多材料线弹性固体的应力分析

本文在基于FVM的应力求解一文基础上，进一步分析多材料线弹性固体的应力问题。

1. 数学模型

前提条件：多材料线弹性固体、绝热并忽略体积力

根据高斯定理以及柯西第一运动定律，在忽略体力的情况下，有：

$${\int }_{V_P}\rho \frac{{\partial }^2\mathbf{u}}{\partial t^2}dV={\oint }_S\mathbf{n}\cdot \sigma dS\text{\hspace{1em}(1-1)}$$

牵引矢量 $\mathbf{t}$ ：

$$\mathbf{t}=\mathbf{n}\cdot \sigma =\mu \mathbf{n}\cdot \nabla \mathbf{u}+\mu \nabla \mathbf{u}\cdot \mathbf{n}+\lambda tr(\nabla \mathbf{u})\mathbf{n}\text{\hspace{1em}(1-2)}$$

将牵引矢量划分为法向(normal)与切向(tangential)两部分。

法向：

$${\mathbf{t}}_n=(2\mu +\lambda )\mathbf{n}\cdot \nabla {\mathbf{u}}_n+\lambda \mathbf{n}tr({\nabla }_t{\mathbf{u}}_t)\text{\hspace{1em}(1-3)}$$

切向：

$${\mathbf{t}}_t=\mu \mathbf{n}\cdot \nabla {\mathbf{u}}_t+\mu {\nabla }_t{u}_n\text{\hspace{1em}(1-4)}$$

其中，${\nabla }_t=(\mathbf{I}-\mathbf{n}\mathbf{n})\cdot \nabla$ —切向梯度算子。

结合式（1-3) 、(1-4）可得：
$$\mathbf{t}=(2\mu +\lambda )\mathbf{n}\cdot \nabla \mathbf{u}-(\mu +\lambda )\mathbf{n}\cdot \nabla {\mathbf{u}}_t+\mu {\nabla }_t{u}_n+\lambda \mathbf{n}tr({\nabla }_t{\mathbf{u}}_t)\text{\hspace{1em}(1-5)}$$

其在数学形式上等价于
$$\mathbf{t}=(2\mu +\lambda )\mathbf{n}\cdot \nabla \mathbf{u}-(\mu +\lambda )\mathbf{n}\cdot \nabla \mathbf{u}+\mu (\nabla \mathbf{u}{)}^T+\lambda tr(\nabla \mathbf{u})\mathbf{n}\text{\hspace{1em}(1-6)}$$

注：可参考基于FVM的应力求解一文中的式（2-18)。

2. 离散方式

离散方式主要分为两种：

• 计算域离散
• 方程离散

2.1 计算域离散

计算域离散分为时间离散和空间离散。 相应内容可参考基于FVM的应力求解 2.1节

假定不同材料之间的界面与控制体内表面重合。则每个控制体可以由普通内表面（f）、交界面（i）、边界面（b）来表示：

$$S_P=\partial V_P=\sum _f{S}_f+\sum _i{S}_i+\sum _b{S}_b\text{\hspace{1em}(2-1)}$$

$S_P$ 为控制体 $V_P$ 的边界面。

2.2 数学模型离散

对于控制体体积为 $V_P$ 的动量方程，其完全离散对应项为：

$${\rho }_p\frac{{\mathbf{u}}_p^n-2{\mathbf{u}}_p^o+{\mathbf{u}}_p^{oo}}{(\nabla t{)}^2}=\sum _f{\mathbf{t}}_f^n{S}_f+\sum _i{\mathbf{t}}_i^n{S}_i+\sum _b{\mathbf{t}}_b^n{S}_b\text{\hspace{1em}(2-2)}$$

• 下角标（subscript）P：网格中心值
• 下角标 $f、i、n$：分别代表内表面(internal)、交界面(interface)、边界面(boundary)的面中心值；
• 上角标（superscript） $n、o、oo$ ：分别代表新值（当前时间步时间）、旧值、旧-旧值。

下面将分别介绍上述三种不同表面的牵引力计算。

2.2.1 内表面牵引力

材料内表面与多材料边界面不重合，其牵引力可通过对式（1-5）离散求得：

$${\mathbf{t}}_f^n=(2{\mu }_f+{\lambda }_f){\mathbf{n}}_f\cdot (\nabla \mathbf{u}{)}_f^n-({\mu }_f+{\lambda }_f){\mathbf{n}}_f\cdot (\nabla {\mathbf{u}}_t{)}_f^n+{\mu }_f({\nabla }_t{u}_n{)}_f^n+{\lambda }_f{\mathbf{n}}_{f}tr({\nabla }_t{\mathbf{u}}_t{)}_f^n\text{\hspace{1em}(2-3)}$$

控制体内表面法向位移导数 ${\mathbf{n}}_f\cdot (\nabla \mathbf{u}{)}_f^n$ 离散形式如下：

$${\mathbf{n}}_f\cdot (\nabla \mathbf{u}{)}_f^n=\underset{orthogonal}{\underbrace{|{\mathbf{\Delta }}_f|\frac{{\mathbf{u}}_N^n-{\mathbf{u}}_P^n}{|\mathbf{d}|}}}+\underset{non-orthogonal}{\underbrace{({\mathbf{n}}_f-{\mathbf{\Delta }}_f)\cdot (\nabla \mathbf{u}{)}_f^n}}\text{\hspace{1em}(2-4)}$$

其中，${\mathbf{\Delta }}_f=\mathbf{d}/(\mathbf{d}\cdot {\mathbf{n}}_f)$。

方程（2-4）中正交部分进行隐式处理，而非正交修正部分则采用显性方法。另外，该方程式也可以应用于控制体内表面切向位移 ${\mathbf{n}}_f\cdot (\nabla {\mathbf{u}}_t{)}_f^n$。

式（2-3）右侧最后两项与前两项的非正交修正部分相同，需要计算控制体位移场的表面梯度，故需要采用显性离散方式。

位移的面中心梯度可通过与相邻网格体中心值线性插值（linear interpolation）求得：

$$(\nabla \mathbf{u}{)}_f=f_x(\nabla \mathbf{u}{)}_P+(1-f_x)(\nabla \mathbf{u}{)}_N\text{\hspace{1em}(2-5)}$$

其中，$f_x$ 为插值系数，

$$f_x=\overline{fN}/\overline{PN}\text{\hspace{1em}(2-6)}$$

位移的网格中心梯度采用离散的高斯定理进行计算：

$$(\nabla \mathbf{u}{)}_P=\frac{1}{V_P}\sum _f{\mathbf{n}}_f{\mathbf{u}}_f{S}_f\text{\hspace{1em}(2-7)}$$

这里的面中心位移矢量 ${\mathbf{u}}_f$ 可与相邻网格进行插值求得。为保证位移矢量具有二阶精度，线 $\overline{PN}$ 必须与面 $f$ 相交于面心位置，否则，必须引入一个偏度（skewness）修正的线性插值。

$${\mathbf{u}}_f=f_x{\mathbf{u}}_P+(1-f_x){\mathbf{u}}_N+{\mathbf{m}}_f\cdot (\nabla \mathbf{u}{)}_f\text{\hspace{1em}(2-8)}$$

如图所示，${\mathbf{m}}_f$ 为偏度修正矢量，是线 $\overline{PN}$ 与面 $f$ 的交点到面心的位移矢量。

注：在离散后的动量方程中，用于偏度、非正交修正以及其他显性项的网格中心偏移梯度均由之前的外迭代求得。

2.2.2 交界面牵引力

求解这一类问题时，需要保证交界面处位移矢量以及牵引力的连续性。

如图，$V_P$ 和 $V_{Ni}$ 为两种不同弹性体的控制体体积，交界面为 $i$，下角标 ${}_{ia}$ 和 ${}_{ib}$ 分别用于指定 $V_P$ 和 $V_{Ni}$ 的特征参数。${\delta }_{an}$ 和 ${\delta }_{bn}$ 分别为控制体$V_P$ 和 $V_{Ni}$ 到交界面的距离。 这里同样将交界面的中心值分为法向和切向两部分。

根据式 (1-3) 以及式（2-4)，可求出交接面两侧法向牵引力分别为：

$$({\mathbf{t}}_n{)}_{ia}=(2{\mu }_{ia}+{\lambda }_{ia})\frac{({\mathbf{u}}_n{)}_i-({\mathbf{u}}_n{)}_P}{{\delta }_{an}}+{\lambda }_{ia}{\mathbf{n}}_{i}tr({\nabla }_t{\mathbf{u}}_t{)}_{ia}\text{\hspace{1em}(2-9)}$$

$$({\mathbf{t}}_n{)}_{ib}=(2{\mu }_{ib}+{\lambda }_{ib})\frac{({\mathbf{u}}_n{)}_{Ni}-({\mathbf{u}}_n{)}_i}{{\delta }_{bn}}+{\lambda }_{ib}{\mathbf{n}}_{i}tr({\nabla }_t{\mathbf{u}}_t{)}_{ib}\text{\hspace{1em}(2-10)}$$

注：为简化运算，假设网格正交

为保证牵引力在界面的连续性，必须保证 $({\mathbf{t}}_n{)}_{ia}=({\mathbf{t}}_n{)}_{ib}$。结合式（2-9）与式（2-10）可求出交界面位移矢量的切向部分 $(u_n{)}_i$：

$$\begin{array}{rl}（{\mathbf{u}}_n{)}_i& =\frac{(2{\mu }_{ia}+{\lambda }_{ia}){\delta }_{bn}({\mathbf{u}}_n{)}_P+(2{\mu }_{ib}+{\lambda }_{ib}){\delta }_{an}({\mathbf{u}}_n{)}_{Ni}}{(2{\mu }_{ia}+{\lambda }_{ia}){\delta }_{bn}+(2{\mu }_{ib}+{\lambda }_{ib}){\delta }_{an}}\\ & +\frac{{\delta }_{an}{\delta }_{bn}\left[{\lambda }_{ib}{\mathbf{n}}_{i}tr({\nabla }_t{\mathbf{u}}_t{)}_{ib}-{\lambda }_{ia}{\mathbf{n}}_{i}tr({\nabla }_t{\mathbf{u}}_t{)}_{ia}\right]}{(2{\mu }_{ia}+{\lambda }_{ia}){\delta }_{bn}+(2{\mu }_{ib}+{\lambda }_{ib}){\delta }_{an}}\end{array}\text{\hspace{1em}(2-11)}$$

将方程（2-11）代入到（2-9）中，有：

$$\begin{array}{rl}({\mathbf{t}}_n{)}_i& =\overline{(2\mu +\lambda {)}_i}\frac{({\mathbf{u}}_n{)}_{Ni}-({\mathbf{u}}_n{)}_P}{{\delta }_{in}}\\ & +\frac{(2{\mu }_{ia}+{\lambda }_{ia}){\delta }_{bn}{\lambda }_{ib}{\mathbf{n}}_{i}tr({\nabla }_t{\mathbf{u}}_t{)}_{ib}+(2{\mu }_{ib}+{\lambda }_{ib}){\delta }_{an}{\lambda }_{ia}{\mathbf{n}}_{i}tr({\nabla }_t{\mathbf{u}}_t{)}_{ia}}{(2{\mu }_{ia}+{\lambda }_{ia}){\delta }_{bn}+(2{\mu }_{ib}+{\lambda }_{ib}){\delta }_{an}}\end{array}\text{\hspace{1em}(2-12)}$$

其中，
$${\delta }_{in}={\delta }_{an}+{\delta }_{bn}\text{\hspace{1em}(2-13)}$$

由调和插值（harmonic interpolation）得到的表面系数：

$$\overline{(2\mu +\lambda {)}_i}=\frac{(2{\mu }_{ia}+{\lambda }_{ia})(2{\mu }_{ib}+{\lambda }_{ib})}{\frac{{\delta }_{bn}}{{\delta }_{in}}(2{\mu }_{ia}+{\lambda }_{ia})+\frac{{\delta }_{an}}{{\delta }_{in}}(2{\mu }_{ib}+{\lambda }_{ib})}\text{\hspace{1em}(2-14)}$$

对牵引力切向部分采用同样的方法，得到切向位移矢量：

$$(u_t{)}_i=\frac{{\mu }_{ia}{\delta }_{bn}({\mathbf{u}}_t{)}_P}{{\mu }_{ia}{\delta }_{bn}+{\mu }_{ib}{\delta }_{an}}+\frac{{\delta }_{an}{\delta }_{bn}\left[{\mu }_{ib}({\nabla }_t{u}_n{)}_{ib}-{\mu }_{ia}({\nabla }_t{u}_n{)}_{ia}\right]}{{\mu }_{ia}{\delta }_{bn}+{\mu }_{ib}{\delta }_{an}}\text{\hspace{1em}(2-15)}$$

切向牵引力：
$$\begin{array}{rl}({\mathbf{t}}_t{)}_i& =\overline{{\mu }_i}\frac{({\mathbf{u}}_t{)}_{Ni}-({\mathbf{u}}_t{)}_P}{{\delta }_{in}}+\\ & +\frac{{\mu }_{ia}{\mu }_{ib}{\delta }_{bn}({\nabla }_t{u}_n{)}_{ib}+{\mu }_{iab}{\mu }_{ia}{\delta }_{an}({\nabla }_t{u}_n{)}_{ia}}{{\mu }_{ia}{\delta }_{bn}+{\mu }_{ib}{\delta }_{an}}\end{array}\text{\hspace{1em}(2-16)}$$

以及表面系数 $\overline{{\mu }_i}$：
$$\overline{{\mu }_i}=\frac{{\mu }_{ia}{\mu }_{ib}}{\frac{{\delta }_{bn}}{{\delta }_{in}}{\mu }_{ia}+\frac{{\delta }_{an}}{{\delta }_{in}}{\mu }_{ib}}\text{\hspace{1em}(2-16)}$$

综上，界面处牵引力为：
$$\begin{array}{rl}{\mathbf{t}}_i& =\overline{(2\mu +\lambda {)}_i}\frac{(\mathbf{u}{)}_{Ni}-(\mathbf{u}{)}_P}{{\delta }_{in}}-\left[\overline{(2\mu +\lambda {)}_i}-\overline{{\mu }_i}\right]\frac{({\mathbf{u}}_t{)}_{Ni}-({\mathbf{u}}_t{)}_P}{{\delta }_{in}}\\ & +\frac{(2{\mu }_{ia}+{\lambda }_{ia}){\delta }_{bn}{\lambda }_{ib}{\mathbf{n}}_{i}tr({\nabla }_t{\mathbf{u}}_t{)}_{ib}+(2{\mu }_{ib}+{\lambda }_{ib}){\delta }_{an}{\lambda }_{ia}{\mathbf{n}}_{i}tr({\nabla }_t{\mathbf{u}}_t{)}_{ia}}{(2{\mu }_{ia}+{\lambda }_{ia}){\delta }_{bn}+(2{\mu }_{ib}+{\lambda }_{ib}){\delta }_{an}}\\ & +\frac{{\mu }_{ia}{\mu }_{ib}{\delta }_{bn}({\nabla }_t{u}_n{)}_{ib}+{\mu }_{ib}{\mu }_{ia}{\delta }_{an}({\nabla }_t{u}_n{)}_{ia}}{{\mu }_{ia}{\delta }_{bn}+{\mu }_{ib}{\delta }_{an}}\end{array}\text{\hspace{1em}(2-17)}$$

上式RHS第一项采用隐式处理，其余各项为显性处理。

• 边界面切向梯度的计算

2.2.3 边界面牵引力

边界面上的牵引力与规定的边界条件有关。对于固定位移边界条件，牵引力为：

$$\begin{array}{rl}{\mathbf{t}}_b^n& =(2{\mu }_b+{\lambda }_b)\frac{{\mathbf{u}}_b^n-{\mathbf{u}}_P^n}{{\delta }_n}-({\mu }_b+{\lambda }_b)\frac{({\mathbf{u}}_t{)}_b^n-({\mathbf{u}}_t{)}_P^n}{{\delta }_n}\\ & +{\mu }_b({\nabla }_t{u}_n{)}_b^n+{\lambda }_b{\mathbf{n}}_{b}tr({\nabla }_t{\mathbf{u}}_t{)}_b^n\end{array}\text{\hspace{1em}(2-18)}$$

其中，下角标 ${}_b$ 代表边界面上的值，${\delta }_n$ 代表相邻网格中心到边界面的距离。

2.2.4 线性代数方程

通过对控制体网格中心与相邻控制体网格中心处的值进行线性插值，可将式（2-2）进行改写：

$$a_P{\mathbf{u}}_P^n+\sum _N{a}_N{\mathbf{u}}_N^n+\sum _{Ni}{a}_{Ni}{\mathbf{u}}_{Ni}^n={\mathbf{R}}_P^n\text{\hspace{1em}(2-19)}$$

注意区分多材料交界面与普通内表面。这里，

$$a_P=\frac{{\rho }_P{V}_P}{(\Delta t{)}^2}+\sum _f(2{\mu }_f+{\lambda }_f)\frac{{\mathbf{\Delta }}_f}{\mathbf{d}}{S}_f+\sum _i\overline{(2\mu +\lambda {)}_i}\frac{1}{{\delta }_{in}}{S}_i\text{\hspace{1em}(2-20)}$$

$$a_N=-(2{\mu }_f+{\lambda }_f)\frac{{\mathbf{\Delta }}_f}{{\mathbf{d}}_f}{S}_f\text{\hspace{1em}(2-21)}$$

$$a_{Ni}=-\overline{(2\mu +\lambda {)}_i}\frac{1}{{\delta }_{in}}{S}_i\text{\hspace{1em}(2-23)}$$

$$\begin{array}{rl}{\mathbf{R}}_P^n& ={\rho }_P{V}_P\left(\frac{2{\mathbf{u}}_P^o}{(\Delta t{)}^2}-\frac{2{\mathbf{u}}_P^{oo}}{(\Delta t{)}^2}\right)+\sum _f(2{\mu }_f+{\lambda }_f)({\mathbf{n}}_f-{\mathbf{\Delta }}_f)\cdot (\nabla \mathbf{u}{)}_f^n{S}_f\\ & -\sum _f({\mu }_f+{\lambda }_f)\left[|{\mathbf{\Delta }}_f|\frac{({\mathbf{u}}_t{)}_N^n-({\mathbf{u}}_t{)}_P^n}{|\mathbf{d}|}+({\mathbf{n}}_f-{\mathbf{\Delta }}_f)\cdot (\nabla {\mathbf{u}}_t{)}_f^n{S}_f\right]\\ & +\sum _f{\mu }_f{\mathbf{n}}_f\cdot ({\nabla }_t{u}_n{)}_f{S}_f+\sum _f{\lambda }_f{\mathbf{n}}_{f}tr[({\nabla }_t{\mathbf{u}}_t{)}_f]S_f\\ & -\sum _i\left[\overline{(2\mu +\lambda {)}_i}-\overline{{\mu }_i}\right]\frac{({\mathbf{u}}_t{)}_{Ni}-({\mathbf{u}}_t{)}_P}{{\delta }_{in}}\\ & +\sum _i\frac{(2{\mu }_{ia}+{\lambda }_{ia}){\delta }_{bn}{\lambda }_{ib}{\mathbf{n}}_{i}tr({\nabla }_t{\mathbf{u}}_t{)}_{ib}+(2{\mu }_{ib}+{\lambda }_{ib}){\delta }_{an}{\lambda }_{ia}{\mathbf{n}}_{i}tr({\nabla }_t{\mathbf{u}}_t{)}_{ia}}{(2{\mu }_{ia}+{\lambda }_{ia}){\delta }_{bn}+(2{\mu }_{ib}+{\lambda }_{ib}){\delta }_{an}}\\ & +\sum _i\frac{{\mu }_{ia}{\mu }_{ib}{\delta }_{bn}({\nabla }_t{u}_n{)}_{ib}+{\mu }_{ib}{\mu }_{ia}{\delta }_{an}({\nabla }_t{u}_n{)}_{ia}}{{\mu }_{ia}{\delta }_{bn}+{\mu }_{ib}{\delta }_{an}}\end{array}\text{\hspace{1em}(2-24)}$$

3. 求解过程

根据式（2-19)，将控制体的每个代数方程组合成在一起，形成代数方程组。其矩阵形式如下：

$$[A][\mathbf{U}]=[\mathbf{R}]\text{\hspace{1em}(2-25)}$$

• [A]：稀疏矩阵（sparse matrix)，系数$a_P$在对角线位置，$a_N、a_{Ni}$ 在非对角线位置；这一矩阵是对称并且对角占优的（求解器：ICCG)；
• $[\mathbf{U}]$：所有控制体的位移矢量集合；
• $[\mathbf{R}]$：源项。

求解一般步骤：

1. 更新时间步，根据上一时间步长的值，初始化因变量；
2. 利用线性插值公式（2-8）计算内表面面心处的位移矢量；用式（2-11）、（2-15）计算交界面处的位移矢量；
3. 根据上一时间步求得的面心位移矢量，利用式（2-7）计算控制体体心位移矢量；
4. 采用线性插值方法计算面心位移梯度；
5. 通过式（2-25）求出控制体位移 $\mathbf{u}$;
6. 判断位移场是否收敛，若收敛，返回步骤1.；否则返回步骤2.。

参考：

1. Tukovi, Ivankovi, A. , & Kara, A. . (2012). Finite-volume stress analysis in multi-material linear elastic body. , aop(aop).
2. https://blog.csdn.net/hanbingchegu/article/details/108378493