空间离散
类似地,得到共计 \(3N\) 个位移自由度
\[
u_{x,1},u_{y,1},u_{z,1},u_{x,2},u_{y,2},u_{z,2},\cdots,u_{x,N},u_{y,N},u_{z,N},
\]
和 \(3N\) 个求解方程
\[\begin{split}
\begin{equation}
\begin{aligned}
&\int_{\Omega}\rho\ddot{\mathbf{u}}\cdot\mathbf{v}_{*,i} \, \mathrm{d}V +
\int_{\Omega}c\dot{\mathbf{u}}\cdot\mathbf{v}_{*,i} \, \mathrm{d}V +
\int_{\Omega} \boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) : \boldsymbol{\sigma} \, \mathrm{d}V + \int_{\Gamma_{R}} \alpha\mathbf{u} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S \\
=& \int_{\Omega}\mathbf{f}\cdot \mathbf{v}_{*,i}\,\mathrm{d}V + \int_{\Gamma_{N}} \tilde{\mathbf{p}} \cdot \mathbf{v}_{*,i}\ \mathrm{d}S + \int_{\Gamma_{R}} \mathbf{f}_{R} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S,\quad \forall \, \mathbf{v}_{*,i} \in \, \mathcal{V},\quad *=x,y,z;\ i=1:N.
\end{aligned}
\end{equation}
\end{split}\]
在静力学求解中,空间离散完得到的是非线性的代数方程,而在动力学求解中,上述方程是一个非线性的微分方程
将上述方程记为
(12)\[
\ddot{\mathbf{F}}^{\text{a}}(t)+\dot{\mathbf{F}}^{\text{v}}(t)+\mathbf{F}^{\text{int}} (t)+ \mathbf{F}^{r}(t) = \mathbf{F}^{\text{ext}}(t),
\]
其中
加速度项(二阶导):\(\ddot{\mathbf{F}}^{\text{a}}(t):=\int_{\Omega}\rho\ddot{\mathbf{u}}\cdot\mathbf{v}_{*,i} \, \mathrm{d}V\)
速度项(一阶导):\(\dot{\mathbf{F}}^{\text{v}}(t):=\int_{\Omega}c\dot{\mathbf{u}}\cdot\mathbf{v}_{*,i} \, \mathrm{d}V\)
内力项向量:\(\mathbf{F}^{\text{int}}(t):=\int_{\Omega} \boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) : \boldsymbol{\sigma} \, \mathrm{d}V\)
Robin 边界条件贡献向量:\(\mathbf{F}^{\text{r}}(t):=\int_{\Gamma_{R}} \alpha\mathbf{u} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S\)
外力向量:\(\mathbf{F}^{\text{ext}}(t):=\int_{\Omega}\mathbf{f}\cdot \mathbf{v}_{*,i}\,\mathrm{d}V + \int_{\Gamma_{N}} \tilde{\mathbf{p}} \cdot \mathbf{v}_{*,i}\ \mathrm{d}S + \int_{\Gamma_{R}} \mathbf{f}_{R} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S\)
为表示方便,这里使用 Vogit 形式,计算推导过程与前文类似,在每个单元上(等式最右端为 \(\mathbf{u}_{E}\),仍记为 \(\mathbf{u}\))
\[\begin{split}
\begin{equation}
\begin{aligned}
&\ddot{\mathbf{F}}^{\text{a}} = \int_{E} \rho\ddot{\mathbf{u}}\cdot\mathbf{v}_{*,i} \, \mathrm{d}E = \int_{E}\rho\mathbf{N}^{T}\ddot{\mathbf{u}}\ \mathbf{d}E= \int_{E}\rho\mathbf{N}^{T}\mathbf{N}\ \mathbf{d}E\cdot\ddot{\mathbf{u}},\\
&\dot{\mathbf{F}}^{\text{v}} = \int_{E} c\dot{\mathbf{u}}\cdot\mathbf{v}_{*,i} \, \mathrm{d}E = \int_{E}\rho\mathbf{N}^{T}\dot{\mathbf{u}}\ \mathbf{d}E= \int_{E}c\mathbf{N}^{T}\mathbf{N}\ \mathbf{d}E\cdot\dot{\mathbf{u}},\\
&\mathbf{F}^{\text{int}} = \int_{E} \boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) : \boldsymbol{\sigma} \, \mathrm{d}E = \int_{E}\mathbf{B}^{T}\boldsymbol{\sigma}\, \mathrm{d}E,\\
&\mathbf{F}^{r}=\int_{\Gamma_{R,E}} \alpha\mathbf{u} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S=\int_{\Gamma_{R,E}} \alpha\mathbf{N}^{T}\mathbf{u}\, \mathrm{d}S=\int_{\Gamma_{R,E}} \alpha\mathbf{N}^{T}\mathbf{N}\, \mathrm{d}S\cdot\mathbf{u},\\
&\mathbf{F}^{\text{ext}}=\int_{E}\mathbf{N}^{T}\mathbf{f}\,\mathrm{d}E + \int_{\Gamma_{N,E}} \mathbf{N}^{T}\tilde{\mathbf{p}}\ \mathrm{d}S + \int_{\Gamma_{R,E}} \mathbf{N}^{T}\mathbf{f}_{R}\, \mathrm{d}S,
\end{aligned}
\end{equation}
\end{split}\]
上式通常写为
\[
\mathbf{M}\ddot{\mathbf{u}}+\mathbf{C}\dot{\mathbf{u}}+\mathbf{F}^{\text{int}}+\mathbf{K}_{r}\mathbf{u}-\mathbf{F}^{\text{ext}} = \mathbf{0},
\]
其中
\[
\begin{equation}
\begin{aligned}
&\mathbf{M} = \int_{E}\rho\mathbf{N}^{T}\mathbf{N}\ \mathbf{d}E,\quad \mathbf{C} = \int_{E}c\mathbf{N}^{T}\mathbf{N}\ \mathbf{d}E,\quad \mathbf{K}_{r}=\int_{\Gamma_{R}^{E}} \alpha\mathbf{N}^{T}\mathbf{N}\, \mathrm{d}S.
\end{aligned}
\end{equation}
\]
对于线弹性本构,有
\[
\mathbf{F}^{\text{int}} =\int_{E}\mathbf{B}^{T}\mathbb{C}\mathbf{B}\ \mathbf{d}E\cdot\mathbf{u}=\mathbf{K}\mathbf{u},
\]
时间离散
通常对 \(\ddot{\mathbf{u}}\) 和 \(\dot{\mathbf{u}}\) 使用中心差分
\[\begin{split}
\begin{equation}
\begin{aligned}
\dot{\mathbf{u}}&=\frac{\mathbf{u}(t+\Delta t) - \mathbf{u}(t-\Delta t)}{2\Delta t}\approx\frac{\mathbf{u}_{n+1} - \mathbf{u}_{n-1}}{2\Delta t},\\
\ddot{\mathbf{u}}&=\frac{\mathbf{u}(t+\Delta t) - 2\mathbf{u}(t) + \mathbf{u}(t-\Delta t)}{\Delta t^2}\approx\frac{\mathbf{u}_{n+1} - 2\mathbf{u}_{n} + \mathbf{u}_{n-1}}{\Delta t^2},
\end{aligned}
\end{equation}
\end{split}\]
于是,方程 (12) 变为
\[
\begin{equation}
\begin{aligned}
\mathbf{M}\frac{\mathbf{u}_{n+1} - 2\mathbf{u}_{n} + \mathbf{u}_{n-1}}{\Delta t^2}+\mathbf{C}\frac{\mathbf{u}_{n+1} - \mathbf{u}_{n-1}}{2\Delta t}
+\mathbf{F}^{\text{int}}+ \mathbf{K}_{r}\mathbf{u}-\mathbf{F}^{\text{ext}} = \mathbf{0}.
\end{aligned}
\end{equation}
\]
显式步求解
(13)\[
\mathbf{M}\frac{\mathbf{u}_{n+1} - 2\mathbf{u}_{n} + \mathbf{u}_{n-1}}{\Delta t^2}+\mathbf{C}\frac{\mathbf{u}_{n+1} - \mathbf{u}_{n-1}}{2\Delta t}+
\mathbf{F}^{\text{int}}_{n}+ \mathbf{K}_{r}\mathbf{u}_{n}-\mathbf{F}^{\text{ext}}_{n} = \mathbf{0}.
\]
此时,方程 (13) 是关于 \(\mathbf{u}_{n+1}\) 的线性方程组
\[
\mathbf{M}_{\text{consitent}}\mathbf{u} = \mathbf{F},
\]
矩阵 \(\mathbf{M}_{\text{consitent}}\) 被称为一致质量矩阵,是一个非对角、稀疏的矩阵,反映了节点之间的动力耦合
在显式动力学中,为了避免求解线性方程组,将一致质量矩阵的质量分配到节点上形成对角矩阵,得到集中质量矩阵
\[\begin{split}
\mathbf{M}_{\text{lumped}}=
\begin{bmatrix}
m_{1}&0&\cdots&0\\
0&m_{2}&\cdots&0\\
\vdots&\vdots&\ddots&\vdots\\
0&0&\cdots&m_{3N}
\end{bmatrix}.
\end{split}\]
由于在隐式动力学中,时间步长通常非常小,因此,使用集中质量矩阵不会带来显著误差,且具备以下重要优势
高效性:大幅提升计算效率(无需求解线性方程组)
稳定性:更大的临界时间步长
低内存:无需存储和求解线性方程组
稳健性:一定程度控制沙漏效应
常用的集中质量矩阵构造方法包括
隐式步求解
\[
\mathbf{M}\frac{\mathbf{u}_{n+1} - 2\mathbf{u}_{n} + \mathbf{u}_{n-1}}{\Delta t^2}+\mathbf{C}\frac{\mathbf{u}_{n+1} - \mathbf{u}_{n-1}}{2\Delta t}+
\mathbf{F}^{\text{int}}_{n+1}+ \mathbf{K}_{r}\mathbf{u}_{n+1}-\mathbf{F}^{\text{ext}}_{n+1} = \mathbf{0}.
\]
通常,在求解过程中,\(\Delta t\) 不是均匀分布的,此时需要使用变步长的中心差分公式