三维常系数线弹性方程的有限元程序实现
一、问题描述
求解平衡态下弹性体的位移场。
- 弹性体受体力和面力，满足力的平衡方程和相应的边界条件；
- 均匀各向同性材料，满足特定本构方程（胡克定律）。
1. 微分形式
有界闭区域$\Omega$满足如下力的平衡方程、边界条件和本构方程：
$$
\begin{cases}
-\partial_j \sigma_{ij}(\vec u) = f_i,~\vec x \in \Omega^{\circ} \\
\vec u = 0,~\vec x \in \Gamma_0 \\
\sigma_{ij}(\vec u) n_j = g_i,~~\vec x \in \Gamma_1,
\end{cases}
$$
其中已知量：
- $\vec f \in L_2(\Omega)$ 弹性体所受体力密度场，
- $\vec g \in L_2(\Gamma_1)$ 弹性体边界所受面力场，
未知量：
- $\vec u$ 平衡时初始区域位移场，
$$
\varepsilon_{ij}(\vec u) = \frac{1}{2}(\partial_i u_j + \partial_j u_i),
$$
$$
\sigma_{ij} = \lambda\varepsilon_{kk}\delta_{ij} + 2\mu\varepsilon_{ij},
$$
其中：
- $\varepsilon$ 无穷小应变张量，
- $\sigma$ 应力张量。
2. 变分形式
考虑变分空间$V={\vec v\in (H_1(\Omega))^3|\vec v|_{\Gamma_0}=0}$。任取 $\vec v$$\in V$对微分形式两边做$L_2$内积，并利用分部积分公式得到如下变分形式：
Find $\vec u \in V$ s.t.
$$
\int_\Omega \sigma(\vec u):\epsilon(\vec v) ~dV= \int_\Omega \vec f \cdot \vec v ~dV + \int_{\Gamma_1} \vec g \cdot \vec v ~dS,~~\forall\vec v\in V
$$
3. 适定性
定义双线性型
$$
a(\vec u,\vec v)=\int_\Omega \sigma(\vec u):\epsilon(\vec v) ~dV,
$$
有界线性泛函
$$
l(\vec v) = \int_\Omega \vec f \cdot \vec v ~dV + \int_{\Gamma_1} \vec g \cdot \vec v ~dS.
$$
分别利用Cauchy-Schwarz不等式和Korn不等式可证明$a(\cdot,\cdot)$的连续性和强制性，再由Lax-Milgram引理可得$\vec u$的存在唯一性。
二、数值实验
1. 微分形式
$$
\begin{cases}
-\partial_j \sigma_{ij}(\vec u) = f_i,~\vec x \in \Omega^{\circ} \\
\vec u = 0,~\vec x \in \Gamma_6 \\
\sigma_{ij}(\vec u) n_j = g_{ki},~~\vec x \in \Gamma_k,~k=1,2,..,5
\end{cases}
$$
- $\Omega = [0,2]^3$
- $\Gamma_k,~k=1,2,..,6$分别为立方体的上下前后右左面
- $\vec f = (0, -1, 0)^T$
- $\vec g_1 = (0,0,\frac{1}{3}y)^T,$ $\vec g_2 = (0,0,-\frac{1}{3}y)^T,$
- $\vec g_3 = (\frac{1}{3}y,0,0)^T,$ $\vec g_4 = (-\frac{1}{3}y,0,0)^T,$
- $\vec g_5 = (0,y,0)^T,$
- $\lambda=\mu=1.$
1.1 精确解
$$
\vec u = [0,~\frac{1}{6} y^2,~0]^T
$$
2. 变分形式
Find $\vec u \in V={\vec v\in (H_1(\Omega))^3|\vec v|_{\Gamma_6}=0}$ s.t.
$$
\int_\Omega \sigma(\vec u):\epsilon(\vec v) ~dV= \int_\Omega \vec f \cdot \vec v ~dV + \sum_{i=1}^5\int_{\Gamma_i} \vec g \cdot \vec v ~dS,~~\forall\vec v\in V
$$
3. 有限元形式
Find $\vec u \in V_h = {\vec v\in (H_1(\Omega))^3|\vec v|_{T} \in P_1(T), \forall T \in \mathbf{T_h}, ~\vec v|_{\Gamma_6}=0}$ s.t.
$$
\int_\Omega \sigma(\vec u):\epsilon(\vec v) ~dV= \int_\Omega \vec f \cdot \vec v ~dV + \sum_{i=1}^{5}\int_{\Gamma_i} \vec g \cdot \vec v dS,~\forall\vec v\in V_h
$$
4. 网格生成
按照.smesh面网格文件格式形成初始立方体面网格，然后利用软件Tetgen自动生成四面体网格文件.mesh。
1 | # Part 1 - node list |
5. 线性系统
5.1 刚度矩阵
将$\vec v = \vec v_k \phi_k$代入变分形式可得方程组
$$
\int_\Omega \sigma_{ij}(\vec u) \partial_j \phi_k ~dV= \int_\Omega f_i \phi_k ~dV + \int_{\Gamma} g_i \phi_k dS,~~i=1,2,3,k=0,1,…,N-1.
$$
继续将$\vec u = \vec u_k \phi_k$代入可得关于未知量$(u_{1,1},u_{1,2},…,u_{N,3})$的3N阶线性方程组。矩阵元素
$$
A_{3m+q,~3k+i}=\int_{\Omega} \lambda \partial_q \phi_m \partial_i \phi_k+\mu\partial_i\phi_m\partial_q\phi_k + \delta_{qi}\mu\partial_j\phi_m\partial_j\phi_k~dV.
$$
实际程序实现采用枚举网格每一单元先构造单元刚度矩阵，然后累加至总体刚度矩阵。后续载荷向量和边界条件项的计算实现也采用同样方式。
5.2 载荷向量
$$
\vec b=(…,\int_{\Omega}f_i \phi_k~dV,…)^T.
$$
5.3 边界条件项
- 对于Dirichlet零边值边界条件：将刚度矩阵对应矩阵元素设为1，载荷向量相应位置设为0；
- Neumann边界条件：载荷向量对应位置加上边界项$\int_\Gamma g_i \phi_k ~dS.$
6. 数值结果
6.1误差阶
终端中程序运行输出结果，并将数据整理形成表格计算各种范数的误差阶情况如下：
| 外接球直径 h~max~ | H~1~范数 | 阶 | L~2~范数 | 阶 | L~1~范数 | 阶 |
|---|---|---|---|---|---|---|
| 2.95804 | 0.171261 | / | 0.05619 | / | 0.185037 | / |
| 1.47902 | 0.077058 | 1.15218000735711 | 0.017648 | 1.67080870426018 | 0.055995 | 1.72444386642493 |
| 0.73951 | 0.035939 | 1.10039447556558 | 0.004975 | 1.82673626527301 | 0.015389 | 1.86339852310954 |
| 0.369755 | 0.017469 | 1.04075324191767 | 0.001302 | 1.93396707715274 | 0.00397 | 1.95468857377173 |
| 0.184877 | 0.008706 | 1.00471509936808 | 0.00033 | 1.98019151891982 | 0.000999 | 1.9905824242379 |
6.2 可视化
将位移场数据存储成.vtk文件，之后利用软件paraview作出三维矢量场图像。
利用有限元软件Freefem编写程序求解相同问题。之后利用位移场移动网格，可绘制出平衡时的网格图像。
Freefem有限元软件只需要程序员编写程序表示出问题的变分形式，之后由软件自动解析生成线性方程组求解变分问题，对于一般的问题求解起来更加方便。PHG和Freefem两者可优势互补，可利用Freefem的运算结果对所编写的PHG程序进行对照校验。
