物质点法基础（一）：基本原理与控制方程

1. 物质点法

物质点法(MPM)将连续体离散成一组质点. 每个质点代表一块材料区域并携带该区域材料的所有物质信息, 如质量, 速度, 应力和应变等. 计算网格仅用于动量方程的求解和空间导数的计算, 它不携带任何物质信息. 在每一个时间步中, 质点与计算网格完全固连(拉格朗日计算步), 可以用标准的有限元法在计算网格上求解物体的运动方程. 计算网格结点的运动方程可以通过将质点的运动量映射到计算网格得到, 求解后将计算网格结点的运动量映射回各质点, 从而得到这些质点在下一时刻的运动量(对流计算步).

物质点法是一种完全的拉格朗日质点类方法, 在每步中质点和计算网格没有相对运动, 避免了欧拉法因非线性对流项产生的数值困难. 质点已经携带了连续体的所有物质信息, 因此物质点法在每个时间步结束时抛弃变形后的计算网格, 在新的时间步中仍可以采用未变形的计算网格.

在并行算法方面, 主要基于消息传递模型(MPI)和信息共享模型(OpenMP)实现.

2. 控制方程

2.1 物体运动和变形的描述

物体在 t = 0 时刻所占据的空间区域称为初始构形. 物体在 t 时刻所占据的空间区域成为现时构形. 为了度量物体的运动, 需要选取一个特定的构形作为基准, 称为参考构形.

在参考构形中质点的矢径 $X$可以表示为

$$X = X_i e_i, i = 1,2,3 \qquad (2.1)$$

其中 $e_i$为直角坐标系的基矢量, $X_i$为参考构形中质点矢径 $X$在 $e_i$上的投影. 质点在参考构形中的矢径 $X$不随时间 $t$变化, $X_i$称为物质坐标或拉格朗日坐标, 它可以作为该质点的标记.

在现时构形中质点 $X$的矢径 $x$可以表示为

$$x = x_i e_i , i = 1, 2, 3 \qquad (2.2)$$

其中 $x_i$为矢径 $x$在 $e_i$上的投影, 坐标 $x_i$给出了质点 $X$在空间中的位置,称为空间坐标或欧拉坐标. 质点 $X$的运动方程可以表示为

$$x_i = x_i(X,t) \qquad (2.3)$$

在拉格朗日描述中, 质点 $X$的位移为

$$u_i = x_i (X, t) - X_i \qquad (2.4)$$

在欧拉描述中, 质点 $X$的位移为

$$u_i = x_i - X_i(x,t) \qquad (2.5)$$

质点的速度等于其矢径 $x$的变化率, 即令 $X$保持不变时矢径 $x$对时间的偏导数. 令 $X$保持不变时物理量对时间的导数称为物质导数, 也称为全导数或拉格朗日导数. 由式 $(2.4)$可以得到质点的速度为

$$v_i = \frac{\partial x_i(X, t)}{\partial t} = \frac{\partial u_i(X, t)}{\partial t} \equiv \frac{\partial u_i}{\partial t} \qquad (2.6)$$

质点的加速度为其速度的物质导数, 即

$$\alpha_i = \frac{\partial v_i(X,t)}{\partial t} = \frac{\partial^2 u_i (X, t)}{\partial t^2} \equiv \ddot{u}_i \qquad (2.7)$$

如果物理量 $F$是用空间坐标 $x$和时间 $t$描述的, 即 $F = F(x,t)$, 可以先利用式 $(2.3)$把它变换为复合函数 $F = F(x(X, t), t)$, 它的物质导数为

$$\begin{aligned} \frac{DF(x, t)}{Dt} &= \frac{\partial F(x, t)}{\partial t} + \frac{\partial F(x, t)}{\partial x_i}\frac{\partial x_i(X, t)}{\partial t} \\ &= \frac{\partial F(x, t)}{\partial t} + \frac{\partial F(x, t)}{\partial x_i}v_i \end{aligned} \qquad (2.8)$$

$F(x, t)$描述的是时刻 $t$空间点 $x$处的物理量, 因此 $\partial F(x, t)/\partial t$表示的是物理量在空间固定点 $x$处的变化率, 称为空间时间导数, 也称为局部导数或欧拉导数. 它反映了物理量的非定常性.

2.2 变形梯度

质点的现时坐标 $x_i$相对于物质坐标 $X_j$的偏导数 $F_{ij} = \partial x_i / \partial X_j$称为变形梯度, 它是一个非对称的二阶张量. 初始构形中由相邻质点 $X$和 $X+dX$构成的无限小线元 $dX$在现时构形中变为

$$dx_i = x_i(X + dX, t) - x_i(X, t) \qquad (2.9)$$

对 $x_i(X + dX, t)$在 $X$处作泰勒展开, 并略去高阶小项得

$$dx_i = \frac{\partial x_i}{\partial X_j} dX_j \qquad (2.10)$$

变形梯度的行列式

$$J = \frac{\partial x_i}{\partial X_j} = \begin{vmatrix} \frac{\partial x_1}{\partial X_1} &\frac{\partial x_1}{\partial X_2} &\frac{\partial x_1}{\partial X_3}\\ \frac{\partial x_2}{\partial X_1} &\frac{\partial x_2}{\partial X_2} &\frac{\partial x_2}{\partial X_3}\\ \frac{\partial x_3}{\partial X_1} &\frac{\partial x_3}{\partial X_2} &\frac{\partial x_3}{\partial X_3}\\ \end{vmatrix} \qquad (2.11)$$

称为雅可比行列式.

引入排列张量 $e_{ijk}$, 可将式 $(2.11)$展开为

$$J = e_{ijk}\frac{\partial x_i}{\partial X_1} \frac{\partial x_j}{\partial X_2} \frac{\partial x_k}{\partial X_3} \qquad (2.12)$$

由行列式的性质和 $J$的定义式 $(2.11)$可得

$$e_{plm}J = e_{ijk}\frac{\partial x_i}{\partial X_p} \frac{\partial x_j}{\partial X_l} \frac{\partial x_k}{\partial X_m} \qquad (2.13)$$

设初始构形中的平行六面体体元的三个线元分别是 $dX$, $\delta X$和 $\Delta X$, 它们在现时构形中分别为 $dx$, $\delta x$和 $\Delta x$, 且

$$dx_i = \frac{\partial x_i}{\partial X_p}dX_p,\ \delta x_j = \frac{\partial x_j}{\partial X_l}\delta X_l,\ \Delta x_k = \frac{\partial x_k}{\partial X_m} \Delta X_m$$

现时构形中由 $dx$, $\delta x$和 $\Delta x$构成的六面体体元的体积为

$$dV = \begin{vmatrix} dx_1 &dx_2 &dx_3\\ \delta x_1 &\delta x_2 &\delta x_3\\ \Delta x_1 &\Delta x_2 &\Delta x_3\\ \end{vmatrix} = JdV_0 \qquad (2.14)$$

因此变形梯度矩阵的行列式 $J$可用来表示变形过程中体元的体积变化, 即

$$J = \frac{dV}{dV_0} \qquad (2.15)$$

可见 $J$表示变形前后体元体积之比.

2.3 变形率

考察一个典型的质点 $P$, 它在时刻 $t$的坐标为 $x_j$, 瞬时速度为 $v_i(x_j, t)$. 与 $P$点相邻的质点 $P'$的坐标为 $x_j+dx_j$, 它相对于 $P$点的相对速度为

$$dv_i = v_i(x_j + dx_j, t) - v_i(x_j, t) = \frac{\partial v_i}{\partial x_j} dx_j \qquad (2.16)$$

式中 $\partial v_i/\partial x_j$为速度梯度张量. 它可分解为对称部分和反对称部分之和, 即

$$\begin{aligned} \frac{\partial v_i}{\partial x_j} &= \frac{1}{2}\left( \frac{\partial v_i}{\partial x_j} - \frac{\partial v_j}{\partial x_i} \right) + \frac{1}{2}\left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right) \\ &= \Omega_{ij} + D_{ij} \end{aligned} \qquad (2.17)$$

式中

$$\Omega_{ij} = \frac{1}{2}\left( \frac{\partial v_i}{\partial x_j} - \frac{\partial v_j}{\partial x_i} \right) \qquad (2.18)$$

$$D_{ij} = \frac{1}{2}\left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right) \qquad (2.19)$$

分别为旋律张量和变形率张量.

与速度梯度分解相似， 相对速度式 $(2.16)$也可分解为

$$dv_i = dv_i^* + dv_i^{**} \qquad (2.20)$$

其中

$$dv_i^* = \Omega_{ij}dx_j, \quad dv_i^{**} = D_{ij}dx_j \qquad (2.21)$$

旋率张量 $\Omega_{ij}$是反对称张量, 它只有 3 个独立的分量, 总存在一个对偶矢量(旋度矢量) $\omega$, 使得

$$\Omega_{ij} = -e_{ijk}\omega_k \qquad (2.22)$$

将上式代入式 $(2.21)$的第一式得

$$dv_i^* = -e_{ijk} \omega_k dx_j = (\omega \times dx)_i \qquad (2.23)$$

上式表明相对速度 $dv_i^*$等于 $\omega$和 $dx$的矢量积. 在质点 $P$的邻域内, 相对速度 $dv_i^*$对应于这个邻域绕过 $P$点的某轴的一个刚体运动, 矢量 $\omega$表示转动的角速度.

可以证明, 变形张量 $D_{ij}$是相对于现时构形定义的柯西应变的速率, 也就是真实应变 $\epsilon_{ij}$的速率，即

$$D_{ij} = \frac{\partial \epsilon_{ij}}{\partial t} \qquad (2.24)$$

2.4 柯西应力

考虑物体在时刻 $t$的现时构形内的一个有向面元 $n\Delta A$, 其法向单位向量 $n$在坐标轴 $e_k$上的投影记为 $n_k$. 面元 $n\Delta A$两侧的介质通过面元相互作用以力元 $\Delta T$, 这个力元除以面元的面积就定义了该面元上的应力矢量 $t^{(n)}$:

$$t_i^{(n)} = \lim_{\Delta A \rightarrow 0} \frac{\Delta T_i}{\Delta A} = \frac{dT_i}{dA} \qquad (2.25)$$

面元 $n\Delta A$与另外三个垂直于坐标轴的面元 $n^{(1)}\Delta A_1$, $n^{(2)}\Delta A_2$和 $n^{(3)}\Delta A_3$构成一个四面体, 且有 $\Delta A_k = \Delta An_k$. 面元 $n^{(k)}\Delta A_k$上的应力矢量记为 $t^{(k)}$, 由平衡条件可得

$$\Delta T_i = t_i^{(n)} \Delta A = t_i^{(k)} \Delta A_k = t_i^{(k)} \Delta An_k = \sigma_{ki} n_k \Delta A \qquad (2.26)$$

其中 $\sigma_{ki} = t_i^{(k)}$是垂直于坐标轴 $x_k$的面元 $n^{(k)} \Delta A_k$上的应力矢量 $t^{(k)}$在坐标轴 $x_i$上的分量. 由式$(2.25)$和式$(2.26)$得

$$t_i^{(n)} = \sigma_{ki}n_k \qquad (2.27)$$

(箭头朝外)

由现时构形中垂直于坐标轴的三个面元上的应力矢量的九个分量 $\sigma_{ij}$定义了一个张量, 称为柯西应力张量. 由微元体关于力矩的平衡条件可以证明柯西应力张量是对称的, 即

$$\sigma_{ij} = \sigma_{ji} \qquad (2.28)$$

柯西应力是定义在现时构形每单位面积上的接触力, 它是与变形相关的真实应力.

2.5 焦曼应力率

考察一承受单向应力 $\sigma$作用并绕 $x_3$轴作刚体转动的杆, 如图2.4所示. 在杆平行于 $x_1$的瞬时，$\sigma_{11} = \sigma$, $\sigma_{22} = 0$, 而当杆转动到平行于 $x_2$轴的位置时, $\sigma_{11} = 0$, $\sigma_{22} = \sigma$. 相对于在空间固定的坐标系 $Px_1x_2$而言, 虽然杆内应力状态并未变化, 但刚体转动改变了柯西应力张量的分量, 导致柯西应力 $\sigma_{ij}$的变化率,无论是空间导数 $\partial \sigma_{ij}/\partial t$还是物质导数的时间偏导数都不等于零. 因此在本构方程中使用与变形率 $D_{ij}$相关联的应力率时, 无论 $\partial \sigma_{ij}/\partial t$还是应力的物质导数都不是一个适当的度量. 用于本构方程中的应力率必须不受刚体转动的影响, 即必须是客观张量.

考虑包含 $P$点的一个作瞬时转动的物质元, 动坐标系 $P\bar{x}_1\bar{x}_2\bar{x}_3$与物质元固连, 跟随物质元一起作瞬时刚体转动. 虽然柯西应力张量在空间固定坐标系 $Px_1x_2x_3$中的变化率不为零, 但它在动坐标系 $P\bar{x}_1\bar{x}_2\bar{x}_3$中的变化率为零, 是一个客观张量.

在所考察的时刻 $t$, 动坐标系与固定坐标系 $x_i$重合, 如图 2.5 所示. 在 $P$点附近的任一相邻质点 $Q$, 其在动坐标系中的坐标 $d\bar{x}_i$不受刚体转动影响, 而它在固定坐标系中的坐标以速率

$$dv_i = \Omega_{ij}d\bar{x}_j \qquad (2.29)$$

变化. 在时刻 $t+dt$, 有

$$dx_i = d\bar{x}_i + dv_i dt = (\delta_{ij} + \Omega_{ij}dt) d\bar{x}_j \qquad (2.30)$$

由式($2.30$)可知, 动坐标系和固定坐标系之间的坐标变换的变换系数为

$$c_{ij} = \delta_{ij} + \Omega_{ij} dt \qquad (2.31)$$

在动坐标系内定义应力率

$$\sigma_{ij}^{\nabla} (t) = \lim_{dt \rightarrow 0} \frac{1}{dt} \left[\bar{\sigma}_{ij}(t+dt) - \bar{\sigma}_{ij}(t)\right] \qquad (2.32)$$

显然物质元作刚体运动时, $\sigma_{ij}^{\nabla}(t) = 0$. 质点 $P$在时刻 $t+dt$相对于固定坐标系的应力张量是

$$\sigma_{ij}(t+dt) = \sigma_{ij}(t) + \frac{\partial \sigma_{ij}}{\partial t}(t) dt \qquad (2.33)$$

利用坐标变换将这个应力张量的分量变换到随体坐标系 $P\bar{x}_1\bar{x}_2\bar{x}_3$中, 有

$$\begin{aligned} \bar{\sigma}_{ij}(t + dt) &= c_{pi} c_{qj} \sigma_{pq}(t + dt) \\ &= (\delta_{pi} + \Omega_{pi} dt)(\delta_{qj} + \Omega_{qj} dt)\left[\sigma_{pq}(t) + \frac{\partial \sigma_{pq}}{\partial t}(t) dt\right] \\ &= \sigma_{ij}(t) + \left(\frac{\partial \sigma_{ij}}{\partial t} + \sigma_{ip}\Omega_{pj} + \sigma_{pj}\Omega_{ip}\right) dt + O(dt^2) \end{aligned} \qquad (2.34)$$

将上式代入式$(2.32)$, 并利用 $\bar{\sigma}_{ij}(t) = \sigma_{ij}(t)$和 $\Omega_{ij}$的反对称性, 最后得

$$\sigma_{ij}^{\nabla} = \frac{\partial \sigma_{ij}}{\partial t} - \sigma_{ip}\Omega_{jp} - \sigma_{pj}\Omega_{ip} \qquad (2.35)$$

按上式定义的应力率叫焦曼应力率, 它是不受刚体转动影响的客观张量, 可以在本构方程中使用.