使用有限元方法求解单胞方程

发表于 2024-04-10 更新于 2025-03-06 分类于多尺度阅读次数：

单胞方程的 statement

单胞方程张量场 $E_{kl}^{(mn)}$ 求解

在三维情景, 线弹性条件下, 单胞方程的弱形式为: \[ \int_{\Theta} L_{ijkl} \frac{\partial H_k^{(mn)}}{\partial y_l} \frac{\partial v_i}{\partial y_j} \mathrm{d} y = -\int_{\Theta} L_{ij(mn)} \frac{\partial v_i}{\partial y_j} \mathrm{d} y. \tag{1} \]

场变量 $H_k^{(mn)}$ 满足周期性边界条件.

[!NOTE]

在单胞方程中, $(mn)$ 是对称的自由指标, 使用 $()$ 标注指标 $mn$ 是因为该组指标表示单胞可能的变形模式 (对于三维问题共有 6 种模式, 每一种模式确定一组偏微分方程).

求解出 $H_{k}^{(mn)}$ 后, 可得到均质化的弹性模量为: \[ L_{ij(mn)}^{\mathrm{hom}} = \frac{1}{|\Theta|} \int_{\Theta} L_{ij(mn)} + L_{ijkl}H_{k,l}^{(mn)} \mathrm{d} y. \tag{2} \]

式中，$H_{k,l}^{(mn)} = \frac{\partial H_{k}^{(mn)}}{\partial y_l}$ 。如果记 \[ E_{kl}^{(mn)} \triangleq I_{kl}^{(mn)} + H_{k,l}^{(mn)} \] 式中, $I_{kl}^{(mn)}=\frac{1}{2}(\delta_{km}\delta_{ln}+\delta_{kn}\delta_{lm})$, 那么式 (2) 可以表示为 \[ L_{ij(mn)}^{\mathrm{hom}} = \frac{1}{|\Theta|} \int_{\Theta} L_{ijkl}E_{kl}^{(mn)} \mathrm{d} y. \]

Reduced Order 单胞方程张量场 $P_{kl}^{mn}$ 求解

降阶多尺度的单胞方程弱形式为 \[ \int_{\Theta} L_{ijkl}(y) P_{kl}^{mn(\alpha)}(y) v_{i,j}(y) ~\mathrm{d}y = \int_{\Theta} L_{ijmn}(y) \chi^{(\alpha)}(y) v_{i,j}(y) ~\mathrm{d}y \quad \alpha=1,2, \ldots,N \tag{3} \]

方程 (1) (3) 的区别在于右端项的积分区域，以及正负号。式 (3) 右端对 partition 求和，可以得到式 (1) 的右端项 (相差一个负号) 。方程 (3) 的边界条件与 (1) 相同。

单胞方程的有限元列式

将单胞方程中的张量分量表示为矩阵形式:

\[ L_{ijkl}(y) \rightarrow [\mathbf{D}](y), \quad \frac{\partial v_i}{\partial y_j} \rightarrow [\mathbf{B}][v], \quad \frac{\partial H_k^{(I)}}{\partial y_l} \rightarrow [\mathbf{B}][u]^{(I)}, \quad L_{ij(I)}(y) \rightarrow [\mathbf{D}](y) [\varepsilon]^{(mn)}. \]

式中, $mn = \{11,22,33,23,13,12\}$. 为方便表示向量分量, 我们使用 Voigt 映射 $(\cdot)$ 将二元组 $(mn)$ 映射为一元有序数列 $\{ 1,2,3,4,5,6 \}$. $[\varepsilon]^{(mn)}$ 为 $6\times 1$ 向量, 对应 $(mn)$ 分量为 1, 其它分量为 0. 例如, 当 $mn=11$ 时, 对应的有限元矩阵方程为:

\[ \int_{Y} [B]^{T} [D](y) [B] \mathrm{d} Y [u]^{(11)} = \int_{Y} [B]^{T}[D](y) [\varepsilon]^{(11)} \mathrm{d} Y. \]

式中, $[\varepsilon]^{(11)} = [1,0,0,0,0,0]^{T}$. 记

\[ \int_{Y} [B]^{T} [D](y) [B] \mathrm{d} Y = [K], \quad \int_{Y} [B]^{T}[D](y) [\varepsilon]^{(11)} \mathrm{d} Y = [f]^{(11)}, \]

则有线性方程组

\[ [K][u]^{(11)} = [f]^{(11)}, \]

以及周期性边界条件对边界节点的约束方程

\[ [C_{\mathrm{per}}][u]^{(11)} = [0] \]

和单位平均应变导致的位移边界条件

\[ [C_{\mathrm{disp}}]^{(11)}[u]^{(11)} = [g]. \]

注意, 约束矩阵 $[C_{\mathrm{per}}]$ 与变形模式 $mn$ 无关.

如果考虑材料坐标系, 那么, 在材料坐标系下表示的材料模量 $[\mathrm{D}']$, 需要通过四阶张量的旋转得到宏观坐标系下 $[\mathrm{D}]$: \[ L_{ijkl} = L'_{mnrs} A_{mi} A_{nj} A_{rk} A_{sl} \] 在求解有限元方程之后, 得到的是宏观坐标系下的均质化之后的材料刚度矩阵.

求解 6 组变形模式下的有限元方程, 得到 $[u]^{(mn)}$ 后, 再代入到均质化弹性张量方程中:

\[ [D^{\mathrm{hom}}] [\varepsilon]^{(mn)} = \frac{1}{|Y|} \int_{Y} [D](y) [\varepsilon]^{(mn)} - [D](y)[B][u]^{(mn)} \mathrm{d} Y. \]

若记

\[ [\overline{D}] = \frac{1}{|Y|} \int_{Y} [D](y)\mathrm{d} Y, \quad [\sigma]^{(mn)} = \frac{1}{|Y|} \int_{Y} [D](y)[B] \mathrm{d} Y~[u]^{(mn)}, \]

则有

\[ [D^{\mathrm{hom}}] = [\overline{D}] - [\sigma^{(11)}, \sigma^{(22)}, \sigma^{(33)}, \sigma^{(23)}, \sigma^{(13)}, \sigma^{(12)}]. \]

式中, $\sigma^{(mn)}$ 是由有限元求解的位移场 $u^{(mn)}$ 计算得到应力场, 再求体平均得到. 对于只有纤维和基体两相材料的单胞, 如果它们的体积分数分别为 $c_{f}$ 和 $c_{m}$, 体平均材料矩阵 $[\overline{D}]$ 等于

\[ [\overline{D}] = \frac{1}{|Y|} \int_{Y} [D](y)\mathrm{d} Y = c_f [D_f] + c_m [D_m]. \]

式中, $[D_f]$ 和 $[D_m]$ 分别为纤维和基体材料的弹性模量矩阵.

各个分块上的 E P 张量, 分别对应求解之后的应变场在分块上的平均, 因此 E P 张量的数值形式都是宏观坐标系下的张量. 此外对

单胞方程的求解

使用有限元方法求解单胞方程需要解决如下问题:

生成周期性边界的网格. 比如三维单胞问题, 左右, 前后, 上下面的网格节点要完全一致, 这样才能施加周期性边界条件
施加周期性边界条件. 周期性边界条件包含两部分:
1. 对面主从节点之间的约束关系
2. 用于施加单位平均应变和消除刚体位移的位移边界条件
有限元方程右端项的B矩阵. 单胞变分方程的线性泛函项包含 B 矩阵, 而常用有限元求解器没有与之对应的载荷边界条件

使用 Gmsh 生成周期性网格

Gmsh 带有命令 Periodic, 可以将主面通过平移或旋转的方式与从面相关联, 由此划分网格时可以保证面与面的节点分布完全相同.

施加周期性边界条件

主从节点信息在 Gmsh 网格 .msh 文件中通过关键词 $Periodic进行标识:

从 .msh 文件中获取关联节点以及所属 entity 的信息, 就可以构造周期性边界条件的约束方程了. 在 ABAQUS 中, 一般是通过 equation 关键词确定多个节点自由度之间的约束关系:

\[ A_1 u_i^1 + A_2 u_i^2 + \ldots + A_N u_i^N = 0. \]

式中, $u_i^n$ 上标 $n$ 表示节点编号, 下标 $i$ 表示 $x,y,z$ 三方向自由度编号. 约束方程的变量只能是节点自由度, 并且右端恒等于 0, 所以 equation 只能添加节点自由度的齐次约束. 若希望添加非齐次的约束方程, 那么就需要引入一个额外的参考节点, 并在该节点处施加一个非零的位移条件.

程序采用的方式可以避免设置额外的参考节点, 平均应变条件则通过单胞上特殊节点的位移边界条件来施加, 这可以通过以下推导得到. 以 $x$ 方向为法线方向的两个相对面为例, 我们有

\[ u_i^{1+} - u_i^{1-} = \overline{\varepsilon}_{i1}, \quad \text{对任意面 1 上的主从节点对.} \]

记单胞在 x 轴的顶点坐标 $(1,0,0)$ 处的位移为 $u_i^{x}$, 原点处位移为 $u_i^o$, 于是有

\[ u_i^x - u_i^o = \overline{\varepsilon}_{i1}. \]

所以,

\[ u_i^{1+} - u_i^{1-} = u_i^x - u_i^o, \quad \text{对任意面 1 上的主从节点对}. \]

整理得到面 1 上的周期性边界条件等价为

\[ \begin{aligned} u_i^{1+} - u_i^{1-} - u_i^x + u_i^o &= 0, \quad \text{对任意面 1 上的主从节点对}; \\ u_i^x - u_i^o &= \overline{\varepsilon}_{i1}. \end{aligned} \]

综合其它面的边界条件, 我们得到最终的约束方程

\[ \begin{aligned} u_i^{k+} - u_i^{k-} - u_i^x + u_i^o &= 0, \quad \text{对任意面 k 上的主从节点对}; \\ u_i^x - u_i^o &= \overline{\varepsilon}_{i1}; \\ u_i^y - u_i^o &= \overline{\varepsilon}_{i2}; \\ u_i^z - u_i^o &= \overline{\varepsilon}_{i3}. \end{aligned} \]

第一组方程是关于节点自由度的齐次方程, 且与施加在单胞的单位平均应变无关, 所以在程序实现时可以不用特别考虑边界条件的具体形式; 其它三组方程则通过在节点 $(0,0,0)$, $(1,0,0)$, $(0,1,0)$, $(0,0,1)$ 处施加位移边界条件实现.

在实现第一组约束方程时, 还需要注意剔除冗余的约束方程, 这是因为我们使用 Lagrange 乘子法求解含约束的有限元方程, 如果出现冗余的约束, 则会导致有限元最终得到的刚度矩阵奇异. 我们将用一个二维平面上的四边形单元说明这一点. 如图所示, 单元节点处的位移分别为 $u_1$, $u_2$, $u_3$, $u_4$

沿 $x$ 方向的约束方程有

\[ u_4 - u_2 = u_3 - u_1, \]

沿 $y$ 方向的约束方程为

\[ u_4 - u_3 = u_2 - u_1. \]

但是, 这两个方程是等价的, 因此如果遍历每组对面上每一个顶点,每一条边的节点, 就会出现过约束的问题. 一种避免过约束的方式如图所示, 按顺序读取节点, 并设置约束方程, 跳过图中虚线表示边线上的节点和空心圆圈表示的顶点处的节点, 就可以保证按照周期性边界条件约束所有的节点, 并且不存在冗余的约束方程.

有限元求解结果分析

程序使用三种方式计算均质化材料矩阵:

求解线性位移边界条件下的单胞方程, 对计算结果进行体平均
求解周期性边界条件约束下的单胞方程, 对计算结果进行体平均
求解渐进展开得到的偏微分方程

方式 1,2 又称为 RVE 方法, 即计算单胞在平均应变下的应力场, 对体积积分得到平均应力, 再根据定义得到均质化的材料矩阵.
方式 2 使用周期性边界条件, 是对相对面上的主从节点施加约束关系, 再在角点处施加位移边界条件, 实现所谓的平均应变.
线性位移条件得到的结果应该略大于周期性边界条件的结果.
方式 3 求解新的偏微分方程, 固定单胞的角点等于去除宏观应变 $\varepsilon^c$ 引入的位移量, 求解得到的 $\chi$ 可以看作是对体积积分得到的材料矩阵进行修正:

\[ [D^{\mathrm{hom}}] = [\overline{D}] - [\sigma^{(11)}, \sigma^{(22)}, \sigma^{(33)}, \sigma^{(23)}, \sigma^{(13)}, \sigma^{(12)}] \]