基于FVM的应力求解

FVM：有限体积法，作为一种有限元处理方法，在弹性力学领域得到了广泛应用。该方法主要利用Navier-Stocks方程对多面体（polyhedral）网格进行空间离散。本文旨在针对线弹性材料边界应力问题进行分析。

本文主要解决单一材料位移矢量的求解问题，关于多材料边界问题的划分可参考另一篇博文多材料线弹性固体的应力分析。

文章目录

1. 数学模型
2.数值方法
- 2.1 计算域离散
- 2.2方程离散
- - 2.2.1 二阶时间导数项
  - 2.2.2拉普拉斯项（Div-Grad）
  - 2.2.3梯度项（显式）
- 2.3 边界条件
- 2.4 求解过程

1. 数学模型

数学模型：线弹性固体。

根据柯西第一运动定律：
∂2(ρu)∂t2−∇⋅σ=ρf(1-1)\frac{\partial^2(\rho \mathbf u)}{\partial t^2}-\nabla\cdot\sigma=\rho\mathbf f\tag{1-1}∂t2∂2(ρu)−∇⋅σ=ρf(1-1)

ρ\rhoρ：材料密度
u\bf uu：位移矢量
σ\sigmaσ：柯西应力张量（Cauchy stress tensor）
f\bf ff：体积力（body force）

应变张量 ε\varepsilonε 可表述为：
ε=12[∇u+(∇u)T](1-2)\varepsilon=\frac 12\left[\nabla \mathbf u+(\nabla \mathbf u)^T\right]\tag{1-2}ε=21[∇u+(∇u)T](1-2)

根据胡克定律，应力应变表达式如下：

σ=2με+λtr(∇u)I=2με+λtr(∇ε)I(1-3)\begin{aligned} \sigma &=2\mu\varepsilon+\lambda tr(\nabla \mathbf u)\mathbf I\\[1.2ex] &=2\mu\varepsilon+\lambda tr(\nabla \varepsilon)\mathbf I \end{aligned}\tag{1-3} σ=2με+λtr(∇u)I=2με+λtr(∇ε)I(1-3)

I\bf II: 单位张量

参数 μ\muμ 和 λ\lambdaλ 与材料的杨氏模量 E 和泊松比(Poisson’s ratio) ν\nuν 有关：

μ=E2(1+ν)(1-4)\mu=\frac{E}{2(1+\nu)}\tag{1-4}μ=2(1+ν)E(1-4)

并且，
λ={νE1−ν2planestressνE(1−2ν)(1+ν)planestrainand3D(1-5)\lambda=\begin{cases} \frac{\nu E}{1-\nu^2} & \text plane\;stress \\ \frac{\nu E}{(1-2\nu)(1+\nu)} & plane\;strain\;and\;3D\end{cases}\tag{1-5}λ={1−ν2νE(1−2ν)(1+ν)νEplanestressplanestrainand3D(1-5)

注：以上部分的推导可参考应力应变基础理论分析以及剪切应力张量。

综上，式（1-1）可整理为：

∂2(ρu)∂t2−∇⋅[μ∇u+μ(∇u)T+λtr(∇u)I]=ρf(1-6)\frac{\partial^2(\rho \mathbf u)}{\partial t^2}-\nabla\cdot[\mu\nabla\mathbf u+\mu(\nabla \mathbf u)^T+\lambda tr(\nabla \mathbf u)\mathbf I]=\rho\mathbf f\tag{1-6}∂t2∂2(ρu)−∇⋅[μ∇u+μ(∇u)T+λtr(∇u)I]=ρf(1-6)

为了求解上述方程，需要明确求解区域、时间、初始条件以及边界条件\color{#00F}{求解区域、时间、初始条件以及边界条件}求解区域、时间、初始条件以及边界条件。其中，初始条件取决于位移矢量 u\bf uu 以及速度矢量 ∂u/∂t\partial u/\partial t∂u/∂t 在初始时刻的分布。而边界条件则分为以下几种情况：

固定位置
对称面
固定压力
固定的牵引力
自由表面

2.数值方法

在空间上，采用具有二阶精度的单元中心非结构性有限体积法，对积分形式进行离散化；在时间上，采用具有二阶精度的隐性离散方法。这种离散方式分为两部分：

计算域离散
方程离散

2.1 计算域离散

计算域离散分为时间离散和空间离散。

时间离散：对于瞬态问题，将时间间隔划分为有限个时间步长 Δt\Delta tΔt，并采用时间步进 (time-marching) 的方式进行求解；
空间离散：通过一系列凸的多面控制体（convex polyhedral CV）来定义求解域，这些控制体相互不重叠并且需要完全填充整个空间；

大部分物理量存储于控制体（网格）中心，然而仍有一部分物理量存储于网格单元面。网格单元面分为两种：

内表面：用于连接且仅连接两个网格单元面。每个内表面由两个控制体所共有；
边界面：与计算域边界重合，仅从属于一个网格单元。

如图，P为计算域中心点， N 为相邻计算域中心点， fff 为两者所共有的内表面；b\bf bb 为P点到N点的位移矢量；Sf\mathbf S_fSf 为该内表面的面矢量。控制体的几何形状完全取决于各顶点位置。

2.2方程离散

方程离散是将特定的偏微分方程整理成相应的代数方程，用以表示空间离散后的物理量。

根据高斯定理，式（1-6）可整理为：

∫VP∂2(ρu)∂t2dV−∮SdS⋅[μ∇u+μ(∇u)T+λtr(∇u)I]=∫VPρfdV(2-1)\int_{V_P}\frac{\partial^2(\rho \mathbf u)}{\partial t^2}dV-\oint_Sd\mathbf S\cdot[\mu\nabla\mathbf u+\mu(\nabla \mathbf u)^T+\lambda tr(\nabla \mathbf u)\mathbf I]=\int_{V_P}\rho\mathbf fdV\tag{2-1}∫VP∂t2∂2(ρu)dV−∮SdS⋅[μ∇u+μ(∇u)T+λtr(∇u)I]=∫VPρfdV(2-1)

注：对于位移向量，其内部耦合项通常采用显式（explicit）处理，以便产生对角占优的稀疏矩阵（spares matrix)，便于迭代求解。

下面对式（1-7）逐项进行离散化分析：

2.2.1 二阶时间导数项

∂2u∂t2=un−2uo+uooΔt2(2-2)\frac{\partial^2 u}{\partial t^2}=\frac{u^n-2u^o+u^{oo}}{\Delta t^2}\tag{2-2}∂t2∂2u=Δt2un−2uo+uoo(2-2)

其中，un=u(t+Δt)u^n=u(t+\Delta t)un=u(t+Δt)、uo=u(t)u^o=u(t)uo=u(t)、uoo=u(t−Δt)u^{oo}=u(t-\Delta t)uoo=u(t−Δt)。

假定各点处的位移值呈线性变化，则：u(x)=uP+(x−xp)⋅(∇u)P\mathbf u(x)=\mathbf u_P+(\mathbf x-\mathbf x_p)\cdot(\nabla \mathbf u)_Pu(x)=uP+(x−xp)⋅(∇u)P。

2.2.2拉普拉斯项（Div-Grad）

∫VP∇⋅(μ∇u)dV=∮SdS⋅(μ∇u)=∑μfSf⋅(∇u)f(2-3)\int_{V_P}\nabla\cdot(\mu\nabla\mathbf u)dV=\oint_Sd\mathbf S\cdot(\mu\nabla\mathbf u)=\sum\mu_f\mathbf S_f\cdot(\nabla\mathbf u)_f\tag{2-3}∫VP∇⋅(μ∇u)dV=∮SdS⋅(μ∇u)=∑μfSf⋅(∇u)f(2-3)

隐式离散\color{#0F0}{隐式离散}隐式离散：令d=PN‾\mathbf d =\overline{PN}d=PN，如果 Sf\mathbf S_fSf 平行于 d\mathbf dd ,则：

Sf⋅(∇u)f=∣Sf∣uN−uP∣d∣(2-4)\mathbf S_f\cdot(\nabla\mathbf u)_f=|\mathbf S_f|\frac{\mathbf u_N-\mathbf u_P}{|\mathbf d|}\tag{2-4}Sf⋅(∇u)f=∣Sf∣∣d∣uN−uP(2-4)

注：上式只对正交网格成立。

根据式（2-4)，可以将式（2-3）整理成以下代数形式：
∫VP∇⋅(μ∇u)dV=∮SdS⋅(μ∇u)=aPuP+∑aNuN(2-5)\int_{V_P}\nabla\cdot(\mu\nabla\mathbf u)dV=\oint_Sd\mathbf S\cdot(\mu\nabla\mathbf u)=a_P\mathbf u_P+\sum a_N\mathbf u_N\tag{2-5}∫VP∇⋅(μ∇u)dV=∮SdS⋅(μ∇u)=aPuP+∑aNuN(2-5)

其中，
aN=μf∣Sf∣∣d∣(2-6)a_N=\mu_f\frac{|\mathbf S_f|}{|\mathbf d|}\tag{2-6}aN=μf∣d∣∣Sf∣(2-6)

而

aP=∑(−aN)(2-7)a_P=\sum( -a_N)\tag{2-7}aP=∑(−aN)(2-7)

显性处理：若矢量 Sf\mathbf S_fSf与d\mathbf dd不平行，则需要额外添加一个显性的非正交修正项，

(∇u)f=fx(∇u)P+(1−fx)(∇u)N(2-8)(\nabla\mathbf u)_f=f_x(\nabla \mathbf u)_P+(1-f_x)(\nabla \mathbf u)_N\tag{2-8}(∇u)f=fx(∇u)P+(1−fx)(∇u)N(2-8)

其中，fx=fN‾/PN‾f_x=\overline{fN}/\overline{PN}fx=fN/PN，为插值系数。

不同于之前的隐式计算，这一项需要根据当前的 ∇u\nabla\bf u∇u（每个时间步的初始时刻）进行计算。

另外，方程 {1-6} 中， ∇⋅[μ(∇u)T]、∇⋅[λtr(∇u)I]\nabla\cdot[\mu(\nabla \mathbf u)^T]、\nabla\cdot[\lambda tr(\nabla \mathbf u)\mathbf I]∇⋅[μ(∇u)T]、∇⋅[λtr(∇u)I]均采用显式处理。

2.2.3梯度项（显式）

采用最小二乘法进行计算。

根据 P点的值和梯度推算 N 点的值，并与改点处的实际值相比较，计算两者间的误差：
eN=ϕN−(ϕP+d⋅(∇ϕ)P)(2-9)e_N=\phi_N-(\phi_P+\mathbf d\cdot(\nabla \phi)_P)\tag{2-9}eN=ϕN−(ϕP+d⋅(∇ϕ)P)(2-9)

求单位误差梯度的平方和：

eP2=∑N(ωNeN)2(2-10)e_P^2=\sum_N(\omega_Ne_N)^2\tag{2-10}eP2=N∑(ωNeN)2(2-10)

权函数（weighting function） ωN=1∣d∣\omega_N=\frac 1{|\bf d|}ωN=∣d∣1。

当 ePe_PeP 取最小值时，认为此时的梯度是准确的。

(∇ϕ)P=∑NωN2G−1⋅d(ϕN−ϕP)(2-11)(\nabla \phi)_P=\sum_N\omega_N^2\mathbf G^{-1}\cdot\mathbf d(\phi_N-\phi_P)\tag{2-11}(∇ϕ)P=N∑ωN2G−1⋅d(ϕN−ϕP)(2-11)

G\bf GG 是 P 点的张量：

G=∑NωN2dd(2-12)\mathbf G=\sum_N\omega_N^2\mathbf d\mathbf d\tag{2-12}G=N∑ωN2dd(2-12)

2.3 边界条件

边界条件大致可分为以下几类：

规定边界的位移矢量u\bf uu ,此时需要根据相邻网格的中点值计算面梯度；
给定对称面，对于现有的控制体，在其对称面的另一侧有互为镜像的另一控制体；
牵引边界条件，明确边界面上的力：

F=∣Sb∣t−Sbp(2-13)\mathbf F=|\mathbf S_b|\mathbf t-\mathbf S_bp\tag{2-13}F=∣Sb∣t−Sbp(2-13)

其中，

F\bf FF：直接计入平衡中的力；
Sb\mathbf S_bSb：指向外侧的边界面面积矢量；
t\bf tt：规定的牵引力；
ppp：压力。

2.4 求解过程

综上，式（2-1）可整理为以下形式：

aPuP+∑aNuN=RP(2-14)a_P\mathbf u_P+\sum a_N\mathbf u_N=\mathbf R_P\tag{2-14}aPuP+∑aNuN=RP(2-14)

将每个控制体装配为一个整体方程。

[A][U]=[R](2-15)[A][\mathbf U]=[\mathbf R]\tag{2-15}[A][U]=[R](2-15)

[A]：稀疏矩阵（sparse matrix)，系数aPa_PaP在对角线位置，aNa_NaN在非对角线位置；这一矩阵是对称并且对角占优的（求解器：ICCG，incomplete Cholesky conjugate gradient solver)；
[U][\mathbf U][U]：所有控制体的位移矢量集合；
[R][\bf R][R]：源项。

上述的离散系统中包含一些显性项，它们的值取决于上一次迭代的位移。此时方程（2-15）不必收敛到一个很小的公差，因为下一步计算会应用并更新这些显性项。保证误差小于规定值即可认为收敛。对于非稳态计算，对于每个时间步都要进行上述操作，并用求得的值作为初始条件重新计算。

然而，对于隐式项，如时间导数和 ∇⋅(μ∇u)\nabla\cdot(\mu\nabla\bf u)∇⋅(μ∇u) 等，经实践证明只是稍稍收敛。为此，需要对隐式项进行多重网格加速。

对于位移矢量 x轴分量，考虑到∇⋅[μ(∇u)T]、∇⋅[λtr(∇u)I]\nabla\cdot[\mu(\nabla \mathbf u)^T]、\nabla\cdot[\lambda tr(\nabla \mathbf u)\mathbf I]∇⋅[μ(∇u)T]、∇⋅[λtr(∇u)I]，其右侧相邻网格：

aE=(2μ+λ)∣S∣∣d∣(2-15)a_E=(2\mu+\lambda)\frac{|\mathbf S|}{|\mathbf d|}\tag{2-15}aE=(2μ+λ)∣d∣∣S∣(2-15)

其上侧相邻网格：
aN=μ∣S∣∣d∣(2-16)a_N=\mu\frac{|\mathbf S|}{|\mathbf d|}\tag{2-16}aN=μ∣d∣∣S∣(2-16)

y轴位移分量与x轴相反。在考略到面试量与坐标轴方向夹角的情况下，这一理论可以适用于任意非结构网格。但它收敛缓慢，并且矩阵 [A] 对于每个位移分量都不相同。为此，需要对方程{1-6}进行改写。其离散形式如下：

∂2(ρu)∂t2−∇⋅(μ∇u)⏟implicit−∇⋅[μ(∇u)T+λtr(∇u)I]⏟explicit=ρf(2-17)\frac{\partial^2(\rho \mathbf u)}{\partial t^2} -\underbrace{\nabla\cdot(\mu\nabla\mathbf u)}_{implicit} -\underbrace{\nabla\cdot[\mu(\nabla \mathbf u)^T+\lambda tr(\nabla \mathbf u)\mathbf I]}_{explicit} =\rho\mathbf f\tag{2-17}∂t2∂2(ρu)−implicit∇⋅(μ∇u)−explicit∇⋅[μ(∇u)T+λtr(∇u)I]=ρf(2-17)

根据式（2-15）、(2-16）对其进行改写：
∂2(ρu)∂t2−∇⋅[(2μ+λ)∇u]⏟implicit−∇⋅[μ(∇u)T+λtr(∇u)I−(μ+λ)∇u]⏟explicit=ρf(2-18)\frac{\partial^2(\rho \mathbf u)}{\partial t^2} -\underbrace{\nabla\cdot[(2\mu+\lambda)\nabla\mathbf u]}_{implicit} -\underbrace{\nabla\cdot[\mu(\nabla \mathbf u)^T+\lambda tr(\nabla \mathbf u)\mathbf I-(\mu+\lambda)\nabla \bf u]}_{explicit} =\rho\mathbf f\tag{2-18}∂t2∂2(ρu)−implicit∇⋅[(2μ+λ)∇u]−explicit∇⋅[μ(∇u)T+λtr(∇u)I−(μ+λ)∇u]=ρf(2-18)

对于所有位移分量，上式的aP、aNa_P、a_NaP、aN 都是相同的。

注：上式只适用于正交网格。

参考：

Jasak, H. , & Weller, H. G. . (2000). Application of the finite volume method and unstructured meshes to linear elasticity. International Journal for Numerical Methods in Engineering, 48(2), 267-287.
https://www.jianguoyun.com/p/DeZ17XoQ9s3ZBhij_kk