1. 动力学方程
基本方程
平衡方程：
\[\sigma_{ij,j} + f_i - \rho u_{i,tt} - \mu u_{i,t} = 0 \quad 在V域内
\]
几何方程：
\[\epsilon_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i}) \quad 在V域内
\]
物理方程：
\[\sigma_{ij} = D_{ijkl}\epsilon_{kl} \quad 在V域内
\]
边界条件：
\[u_i = \bar{u}_i \quad 在 s_u边界上
\]
\[\sigma_{ij}n_j = \bar{T}_i \quad 在s_{\sigma}边界上
\]
初始条件：
\[u_i(x,y,z,0) = \bar{u}_i(x,y,z)
\]
\[u_{i,t}(x,y,z,0) = \bar{u}_{i,t}(x,y,z)
\]
解释
$ \sigma_{ij,j}$ \(\sigma_{ij,j}\)分量
对于每个方向\(i\)（\(i=1,2,3\)），应力张量的散度\(\sigma_{ij,j}\) 可以写为：
\[\sigma_{ij,j} = \frac{\partial \sigma_{i1}}{\partial x_1} + \frac{\partial \sigma_{i2}}{\partial x_2} + \frac{\partial \sigma_{i3}}{\partial x_3}
\]
具体地，我们可以将其展开为三个分量：
对于\(x\) 方向 (\(i=1\))：
\[\sigma_{1j,j} =
\sigma_{11,1} + \sigma_{12,2} + \sigma_{13,3} \\
\qquad \qquad =\frac{\partial \sigma_{11}}{\partial x_1} + \frac{\partial \sigma_{12}}{\partial x_2} + \frac{\partial \sigma_{13}}{\partial x_3}
\]
这表示作用在\(x\) 方向上单位体积内的净内力。
同理，可以展开y，z，然后
\[\begin{equation}
\sigma_{ij,j} =
\left\{
\begin{aligned}
&\sigma_{11,1} + \sigma_{12,2} + \sigma_{13,3} &&= \frac{\partial \sigma_{11}}{\partial x_1} + \frac{\partial \sigma_{12}}{\partial x_2} + \frac{\partial \sigma_{13}}{\partial x_3} && \text{(对于 $x$ 方向)} \\
&\sigma_{21,1} + \sigma_{22,2} + \sigma_{23,3} &&= \frac{\partial \sigma_{21}}{\partial x_1} + \frac{\partial \sigma_{22}}{\partial x_2} + \frac{\partial \sigma_{23}}{\partial x_3} && \text{(对于 $y$ 方向)} \\
&\sigma_{31,1} + \sigma_{32,2} + \sigma_{33,3} &&= \frac{\partial \sigma_{31}}{\partial x_1} + \frac{\partial \sigma_{32}}{\partial x_2} + \frac{\partial \sigma_{33}}{\partial x_3} && \text{(对于 $z$ 方向)}
\end{aligned}
\right.
\end{equation}
\]
平衡方程
给出的方程是连续介质力学中描述材料内部应力、外力和加速度之间关系的偏微分方程，具体来说它是一个线性动弹性理论中的运动方程。
让我们逐步解析这个方程：
\[\sigma_{ij,j} + f_i - \rho u_{i,tt} - \mu u_{i,t} = 0
\]
这里：
- $ \sigma_{ij} $是应力张量，描述了作用在材料内的内力分布。
- 下标 $ ij,j $表示对第二个下标的求和，并且对这个下标进行空间偏导数运算（根据Einstein求和约定）。所以 $ \sigma_{ij,j} $实际上表示的是应力张量的散度，它给出了单位体积上的净内力。
- $ f_i $是体力（如重力或电磁力）的分量，作用于单位体积的材料上。
- $ \rho $是材料的密度。
- $ u_i $是位移向量的分量，表示材料中某点相对于其未变形位置的移动。
- $ u_{i,tt} $表示位移向量的时间二阶导数，即加速度。
- $ u_{i,t} $表示位移向量的时间一阶导数，即速度。
- $ \mu $是一个阻尼系数，用于描述材料内部摩擦导致的能量耗散。
方程右边等于零意味着在考虑的所有效应（内力、外力、惯性和阻尼）平衡时，系统的净加速度为零。
\(\epsilon_{ij}\) 分量描述
应变张量 \(\mathbf{\epsilon}\) 是一个二阶张量，它描述了物体内部各点相对于未变形状态的位移梯度。
在数学上，应变张量可以表示为：
\[\epsilon_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)
\]
其中，\(u_i\) 是物体内部某点在 \(i\) 方向上的位移分量，\(x_j\) 是该点的坐标分量。
应变张量的分量解释
- \(\epsilon_{ij}\) 表示在 \(j\) 方向单位长度上，由于形变而在 \(i\) 方向产生的相对位移量的一半。
- 当 \(i = j\) 时，\(\epsilon_{ii}\)（如 \(\epsilon_{11}\)，\(\epsilon_{22}\)，\(\epsilon_{33}\)）表示物体在对应方向（如 \(x_1\)，\(x_2\)，\(x_3\)）上的正应变或线应变。它描述了物体在该方向上长度的相对变化。
- 当 \(i \neq j\) 时，\(\epsilon_{ij}\)（如 \(\epsilon_{12}\)，\(\epsilon_{23}\)，\(\epsilon_{31}\)）表示物体在 \(i\) 和 \(j\) 方向之间的剪应变。它描述了物体在这两个方向之间相对角度的变化或剪切变形。
应变张量的性质
- 对称性：由于应变张量的定义中包含 \(\frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)\)，因此 \(\epsilon_{ij} = \epsilon_{ji}\)，即应变张量是对称的。
- 无量纲：应变张量的分量是无量纲的，因为它们表示的是相对位移或相对变化。
2. 求解基本步骤
现以三维实体动力分析为例，用有限元法求解的基本步骤如下：
(1）连续区域的离散化
在动力分析中，因为引入了时间坐标，处理的是四维（x，y，z，t）问题。在有限元分析中一般采用部分离散的方法，即只对空间域进行离散，这样一来，此步骤和静力分析时相同。
(2）构造插值函数
由于只对空间域进行离散，所以单元内位移u、v、w的插值分别表示为：
\[u(x,y,z,t)= \sum_{i=1} ^ n N_i(x,y,z) * u_i(t)
\]
\[v(x,y,z,t)= \sum_{i=1} ^ n N_i N_i(x,y,z) * v_i(t)
\]
\[w（(x,y,z,t)= \sum_{i=1} ^ n N_i(x,y,z) * w_i(t)
\]
或者写成矩阵形式：
\[u = N * a
\]
其中，u、v、w分别为x、y、z方向上的位移向量，N为形函数矩阵，a为节点位移向量（包含时间函数）。
\[N = [N_1, N_2, \cdots, N_n], \quad \quad N_i = NI_{3 \times 3}, \quad i = 1,2,\cdots,n
\]
\[a^e = \begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_n \end{bmatrix}, \quad a_i = \begin{bmatrix} u_i(t) \\ v_i(t) \\ w_i(t) \end{bmatrix}, \quad i = 1,2,\cdots,n
\]
(3) 形成系统的求解方程
平衡方程式及力的边界条件的等效积分形式的伽辽金提法可表示如下：
\[\int_V
(\delta u_i (\sigma_{ij,j}+f_i - \rho u_{i,tt}
- \mu u_{i,t})dV -
\int_{S_2}\delta u_i\left(\sigma_{ij}n_j - \overline{T}_i\right)ds = 0
\]
对上式的第1项\(\int_V\delta u_i\sigma_{ij,j}dV\)进行分部积分，并代入物理方程，则从上式可以得到：
\[\int_V\left(\delta \epsilon_{ij}D_{ijkl}\epsilon_{kl}+\delta \sigma_{ij}\mu_{ij,k}+\delta \sigma_{ij}u_{i,j}\right)dV = \int_V\delta u_i f_i dV + \int_{S_2}\delta u_i\overline{T}_i ds
\]
将空间离散后的位移表达式（现在情况下，\(u_1=u\)，\(u_2=v\)，\(u_3=w\)）代人上式，并注,意到结点位移变化$ \delta a $的任意性，最终得到系统的求解方程（在动力学问题中，又称运动方,程）如下：
\[M\ddot{a}(t)+C\dot{a}(t)+Ka(t)=Q(t)
\]
如果忽略阻尼的影响，即令阻尼矩阵\(C\)为零，则原运动方程，可以简化为：
\[M\ddot{a}(t) + Ka(t) = Q(t)
\]
其中，\(\ddot{a}(t)\)是系统的结点加速度向量，\(a(t)\)是系统的结点位移向量，\(M\)是系统的质量矩阵，\(K\)是系统的刚度矩阵，\(Q(t)\)是系统的结点载荷向量。这些矩阵和向量分别由各自的单元矩阵和向量集成，即
\[M = \sum M^e, \quad C = \sum C^e, \quad K = \sum K^e, \quad Q = \sum Q^e
\]
其中，\(M^e\)、\(C^e\)、\(K^e\)和\(Q^e\)分别是单元的质量矩阵、阻尼矩阵、刚度矩阵和载荷向量，它们分别由以下公式计算：
\[M^e=\int_V \rho N^TN dV
\]
\[C^e=\int_V \mu N^TN dV
\]
\[K^e=\int_V B^T DB dV
\]
\[Q^e=\int_V N^T f dV + \int_{S_2} N^T T ds
\]
在这里，\(\rho\)是材料的密度，\(N\)是形函数矩阵，\(B\)是应变-位移矩阵，\(D\)是弹性矩阵，\(f\)是体积力向量，\(T\)是面力向量，\(V\)是单元体积，\(S_2\)是单元受面力作用的面。
(4)求解方程
求解方程在后续展开
(5)计算结构应力和应变
从以上步骤可以看出，和静力分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。其他的计算步骤和静力分析是完全相同的。
关于二阶常微分方程组的解法,原则上可利用求解常微分方程组的常用方法(例如Runge-Kutta方法)求解,但是在有限元动力分析中,因为矩阵阶数很高,用这些常用算法一般是不经济的,所以只对少数有效的方法感兴趣,和前一章关于一阶常微分方程组的求解方法相同,这些方法可分为两类,即直接积分法和振型叠加法。
3. 质量矩阵和阻尼矩阵
3.1 协调质量矩阵和集中质量矩阵
上节所表达的单元质量矩阵
\[ M^{\prime}=\int_{v}\rho N^{T}NdV
\]
称为协调质量矩阵或一致质量矩阵，这是因为导出它时，和导出刚度矩阵所根据的原理（伽辽金方法）及所采用位移插值函数是一致的。此外，在有限元法中还经常采用所谓集中（或团集）质量矩阵。它假定单元的质量集中在结点上，这样得到的质量矩阵是对角线矩阵。
将单元协调质量矩阵\(M^{\prime}\)转换为单元集中质量矩阵\(M_i\)，即对\(M^{\prime}\)进行对角化的方法，有多种方案可供选择。以下分实体单元和结构单元进行讨论。
-
实体单元
现在介绍两种常用的方法，即
（1）第一种方法
\[\left(M_i\right)_{ij}=\begin{cases}
\sum_{k=1}^{n}\left(M^{\prime}\right)_{kk}=\sum_{k=1}^{n}\int_{v}\rho N_i^TN_kdV&(j=i)\\
0&(j\neq i)
\end{cases}
\]
其中，\(n\)是单元的结点数。该式的力学意义是：\(M_i\)中每一行的主元素等于\(M^{\prime}\)中该行所有元素之和，而非主元素为零。
（2）第二种方法
\[\left(M_i\right)_{ij}=\begin{cases}
a\left(M^{\prime}\right)_{ii}=a\int_{v}\rho N_i^TN_idV&(j=i)\\
0&(j\neq i)
\end{cases}
\]
此式的力学意义是：\(M_i\)中每一行的主元素等于\(M^{\prime}\)中该行主元素乘以缩放因子\(a\)，而非主元素为零。因子\(a\)根据质量守恒原则确定，即\(M_i\)中对应于每一方向的所有自由度的元素之和应等于整个单元的质量。对于实体单元，则有
\[\sum_{i=1}^{n} \left(M_i\right)_z = a \sum_{i=1}^{n} \left(M'\right)_{zz} = W L_z = \rho V L_z
\]
例 13.1 计算平面应力（应变）单元的协调质量矩阵 \(M'\) 和集中质量矩阵 \(M_i\)。单元
形式采用第 2 章 2.2 节所述的 3 结点三角形单元。
（1）协调质量矩阵
位移插值函数是
\[ N = [N_1, N_2, N_3] I
\]
其中 \(I\) 是 \(2 \times 2\) 单位矩阵。
\[ N_i = (a + b x_i + c y_i)/2A \quad (i = 1, 2, 3)
\]
系数 \(a\), \(b\), \(c\) 见（2.2.7）式；\(A\) 是三角形单元面积。
按式可以算得单元的协调质量矩阵
\[M' = \frac{W}{3} \begin{bmatrix}
\frac{1}{2} & 0 & \frac{1}{4} & 0 & \frac{1}{4} & 0 \\
0 & \frac{1}{2} & 0 & \frac{1}{4} & 0 & \frac{1}{4} \\
\frac{1}{4} & 0 & \frac{1}{2} & 0 & \frac{1}{4} & 0 \\
0 & \frac{1}{4} & 0 & \frac{1}{2} & 0 & \frac{1}{4} \\
\frac{1}{4} & 0 & \frac{1}{4} & 0 & \frac{1}{2} & 0 \\
0 & \frac{1}{4} & 0 & \frac{1}{4} & 0 & \frac{1}{2}
\end{bmatrix}
\]
其中，$$W=\rho tA$$是单元的质量，t是单元的厚度。
（2）集中质量矩阵
①得到集中质量矩阵为
\[M_{i}^{e}=\frac{W}{3}\begin{bmatrix}
1&0&0&0&0&0\\
0&1&0&0&0&0\\
0&0&1&0&0&0\\
0&0&0&1&0&0\\
0&0&0&0&1&0\\
0&0&0&0&0&1
\end{bmatrix}
\]
此式的力学意义是：在单元的每个结点上集中\(\frac{1}{3}\)的质量。
按形成的集中质量矩阵，可以得到因子\(α = 2\)，这样得到的集中质量矩阵。虽然对于3结点三角形单元，形成的质量矩阵完全相同。但是对于一般形式的单元并不是如此。
例如对于8结点矩形单元，对\(M^{\prime}\)进行对角化，所得到的集中质量矩阵\(M\)就不相同。对于角结点：
第1种
\[(M^{e})_{ii}=-\frac{1}{12}W_{i}^{e}\quad(i = 1,2,3,4)
\]
第2种
\[(M^{e})_{ii}=\frac{1}{36}W_{i}^{e}\quad(i = 1,2,3,4)
\]
对于边中点：
由方法1得到的是
\[(M^{e})_{jj}=\frac{1}{3}W_{j}^{e}\quad(j = 5,6,7,8)
\]
由方法2式得到的是
\[(M^{e})_{jj}=\frac{2}{9}W_{j}^{e}\quad(j = 5,6,7,8)
\]
方法1在角结点上\((M^{e})_{ii}\)是负值，这在力学上是不合理的，同时也将对计算精度产生不良影响。因此在实际分析中，更多的是推荐用方法2来计算集中质量矩阵。
- 结构单元
对\(M^{\prime}\)进行对角化常用的简便方法是忽略\(M^{\prime}\)中对应于转动自由度的元素。对于\(M^{\prime}\)中与位移自由度相关的元素则采用方法1、方法2进行处理。例如对于9.2节中所讨论的2结点梁单元、协调质量矩阵和集中质量矩阵如下所示：
(1) 协调质量矩阵
位移插值函数是
\[N = [N_{1},N_{2},N_{3},N_{4}]$ (13.2.7)
\]
其中
\[N_{1}=1-\frac{3x^{2}}{l^{2}}+\frac{2x^{3}}{l^{3}}
\]
\[N_{2}=x-\frac{2x^{2}}{l}+\frac{x^{3}}{l^{2}}
\]
\[N_{3}=\frac{3x^{2}}{l^{2}}-\frac{2x^{3}}{l^{3}}
\]
\[N_{4}=-\frac{x^{2}}{l}+\frac{x^{3}}{l^{2}}
\]
按3式可以算得单元的协调质量矩阵为
浙公网安备 33010602011771号