16.4.1 方程

16.4.1.1 线性等温和等熵弹性

考虑一个被边界 $\Omega$ 包围的固体物体 $\sum$。通常情况下，体积力作用在 $\sum$ 上，而表面力和/或位移则施加在 $\Omega$ 上。设 $\underline{u}$ 表示在 $\Sigma$ 上定义的位移矢量，该矢量是由施加在 $\Omega$ 上的力和位移引起的。从 $\underline{u}$ 出发，可以计算出变形张量 $\underline{\epsilon }$，它是 $\underline{u}$ 梯度的对称部分：

$$\underline{\epsilon }=\frac{1}{2}\left(\underline{\nabla u}+{\left(\underline{\nabla u}\right)}^t\right)\text{\hspace{1em}(16.1)}$$

在当前线性等温和各向同性弹性背景下，应力$\underline{\underline{\sigma }}$的本构方程由以下给出：

$$\underline{\sigma }=\frac{E}{1+v}\left(\frac{v}{1-2v}\mathrm{tr}\left(\underline{\epsilon }\right)\underline{I}+\underline{\epsilon }\right)\text{\hspace{1em}(16.2)}$$

其中，$\underline{I}$ 是单位张量。该方程涉及两个物理参数：杨氏模量 $E$ 和泊松比 $v$。

将应力 $\underline{\sigma }$ 的这一表达式引入动量方程，并在用其定义替换 $\underline{\epsilon }$ 之后，最终得到关于位移 $\underline{u}$ 的二阶椭圆方程。替换后，动量方程可写为：

$$0=-\rho \underline{a}+\nabla \cdot \underline{\sigma }=-\rho \underline{a}+\frac{E}{1+\nu }\nabla \cdot \left(\frac{\nu }{1-2\nu }\mathrm{tr}\left(\underline{\nabla u}\right)\underline{I}+\frac{1}{2}\left(\underline{\nabla u}+{\left(\underline{\nabla u}\right)}^t\right)\right)\text{\hspace{1em}(16.3)}$$

其中，$\rho$ 表示固体的密度，$\underline{a}$ 表示加速度。该方程要求在边界条件中给出位移条件或力条件。因此，在固体被固定的边界上，将指定零位移，而在自由悬挂的墙壁上可以施加节点力或连续力。在流固耦合（FSI）模拟中，在固体与周围流体接触的边界上，将指定由流体产生的连续力作为边界条件。

16.4.1.2 冯·米塞斯应力的评估

冯·米塞斯应力在用于预测延性材料（如金属）屈服现象的冯·米塞斯准则中被引用。其计算公式为：

$${\sigma }_e={\left(\frac{1}{2}\left[{\left({\sigma }_{xx}-{\sigma }_{yy}\right)}^2+{\left({\sigma }_{yy}-{\sigma }_{zz}\right)}^2+{\left({\sigma }_{zz}-{\sigma }_{xx}\right)}^2+6\left({\sigma }_{xy}^2+{\sigma }_{yz}^2+{\sigma }_{xz}^2\right)\right]\right)}^{\frac{1}{2}}\text{\hspace{1em}(16.4)}$$

该方程是代数形式的，因此可以直接从先前计算的应力分量中评估出冯·米塞斯应力。

对于二维平面应力状态，关系式变为：

$${\sigma }_e={\left(\frac{1}{2}\left[{\left({\sigma }_{xx}-{\sigma }_{yy}\right)}^2+{\sigma }_{xx}^2+{\sigma }_{yy}^2+6{\sigma }_{xy}^2)\right]\right)}^{\frac{1}{2}}={\left({\sigma }_{xx}^2+{\sigma }_{yy}^2-{\sigma }_{xx}{\sigma }_{yy}+3{\sigma }_{xy}^2\right)}^{\frac{1}{2}}\text{\hspace{1em}(16.5)}$$

冯·米塞斯应力可用于线性和非线性弹性情况。它最初是为延性材料（如金属）开发的。将其应用于其他固体材料则需要验证和常识判断。

16.4.2 有限元表示

除了少数简单情况外，方程16.3（第827页）无法解析求解，因此必须采用数值方法。如前所述，线性弹性问题属于椭圆型问题，即其行为主要由边界值（而非初始值）决定。因此，使用有限元方法来求解固体中的应力和变形是相当自然的。在有限元框架下，被分析的结构被离散化（近似）为由$N$个离散区域（元素）组成，这些区域在有限数量的点（节点）$\left(M\right)$处连接。对于2D情况，这些元素是三角形或四边形，而对于3D情况，则使用四面体、六面体、楔形或金字塔形。

在节点上定义位移未知数，并通过每个元素层面的适当插值函数进行连接。这些插值函数在每个元素层面以及其一阶导数上是连续的；它们被选择得在相邻元素之间的界面处也能保证连续性。在许多应用中，考虑连续的分段线性插值函数。这样的插值函数${\varphi }_i$与未知节点位移向量$u_i$一样多，特别是每个插值函数在节点$i$处${\varphi }_i=1$，而在其他节点处$=0$。

16.4.2.1 系统矩阵的构建

首先，我们暂时忽略方程16.3（第827页）中的加速度项。设$\tilde{\underline{u}}$为位移场的近似解，通过节点值$u_i$和插值函数${\varphi }_i$来描述。插值函数和近似解$\tilde{\underline{u}}$用于将偏微分方程16.3（第827页）转化为包含$M$个向量方程和$M$个向量未知数的线性代数系统。为此，我们采用以下步骤：

• 通过将方程16.3（第827页）与$M$个插值函数相乘，生成$M$个方程。

• 应用部分积分规则和格林定理，降低导数阶数，并便于将力边界条件融入系统中。

• 最终，将基于节点值$u_i$和插值函数描述的位移场近似解引入系统。

经过离散化处理，方程16.3（第827页）的形式变为

$$0=-{\int }_{\Omega }{\varphi }_i{\underline{f}}_{\Omega }d\Omega +{\int }_{\sum }\frac{E}{1+v}\nabla {\varphi }_i\cdot \left(\frac{v}{1-2v}\mathrm{tr}\left(\underline{\nabla \underline{u}}\right)\underline{I}+\frac{1}{2}\left(\underline{\nabla \underline{u}}+{\left(\underline{\nabla \underline{u}}\right)}^t\right)\right)d\sum \text{\hspace{1em}(16.6)}$$

其中，$f_{\Omega }$ 表示施加在边界 $\Omega$ 上的外力向量。为了便于理解，我们可以将公式 16.6（第 828 页）重新排列为：

$$\left[K\right]\{u\}=\{f\}\text{\hspace{1em}(16.7)}$$

在公式16.7（第828页）中，$\left[K\right]$ 是总刚度矩阵，$\{u\}$ 是节点位移未知量 $u_i$ 的向量，而 $\{f\}$ 则来自施加在固体边界上的外部（连续或节点）力。特别是在FSI（流固耦合）模拟中，它们包含了由流体运动产生的力。加法性质允许在公式16.6（第828页）中的积分在每个单独的单元层面进行，以获得局部刚度矩阵，随后将其组装成线性系统的总刚度矩阵 $\left[K\right]$。在完成总刚度矩阵的组装后，通过从系统中移除相应的未知量和方程来施加位移边界条件。这些方法有详细的文档记录（例如，参见Zienkiewicz [732]（第1099页））。在当前背景下解决系统问题时，使用了一种迭代求解器。

16.4.2.2 动态结构系统

在线性结构动力学系统中，内力与节点位移成线性比例关系，且结构刚度矩阵保持不变。然而，动量方程中应考虑惯性项，以及可能的结构阻尼项。离散化的公式16.7（第828页）现在变为

$$\left[M\right]\{\ddot{u}\left(t\right)\}+\left[C\right]\{\dot{u}\left(t\right)\}+\left[K\right]\{u\left(t\right)\}=\{f\left(t\right)\}\text{\hspace{1em}(16.8)}$$

其中，$\left[M\right]$、$\left[\mathcal{C}\right]$ 和 $\left[K\right]$ 分别代表结构质量矩阵、结构阻尼矩阵和刚度矩阵，而 $\{\ddot{u}\left(t\right)\}$、$\{\dot{u}\left(t\right)\}$ 和 $\{u\left(t\right)\}$ 则分别表示节点加速度、速度和位移向量，$\{f\left(t\right)\}$ 是施加的载荷向量。在线性结构动力学系统中，矩阵是常数，而向量是时间的函数。目前，结构阻尼矩阵 $\left[\mathcal{C}\right]$ 并未被考虑，但为了保持方程的普遍性，它仍被保留在方程中。

在数值求解有限元半离散运动方程（如公式16.8（第829页）所示）的直接时间积分方法中，文献中提供了几种方法；其中两种方法，即纽马克方法（Newmark [475]（第1084页））和后向欧拉方法，在Fluent中可用，并在后续章节中进行了描述。

16.4.2.2.1 纽马克方法

纽马克系列时间积分算法（Newmark [475]（第1084页））是最流行的时间积分方法之一，作为一种单步算法。正如Hughes [261]（第1072页）所指出的，半离散公式16.8（第829页）可以重写为：

$$\left[M\right]\left\{{\ddot{u}}_{n+1}\right\}+\left[C\right]\left\{{\dot{u}}_{n+1}\right\}+\left[K\right]\left\{u_{n+1}\right\}=\left\{f_{n+1}\right\}\text{\hspace{1em}(16.9)}$$

其中，$\left\{{\ddot{u}}_{n+1}\right\}$、$\left\{{\dot{u}}_{n+1}\right\}$ 和 $\left\{u_{n+1}\right\}$ 分别表示在时间 $t_{n+1}$ 时的节点加速度、速度和位移向量，而 $\left\{f_{n+1}\right\}$ 则表示在时间 $t_{n+1}$ 时施加的载荷向量。除了方程 16.9（第 829 页）之外，Newmark 系列时间积分算法还要求更新位移和速度，如下所示：

$$\left\{{\dot{u}}_{n+1}\right\}=\left\{{\dot{u}}_n\right\}+\left[\left(1-\delta \right)\left\{{\ddot{u}}_n\right\}+\delta \left\{{\ddot{u}}_{n+1}\right\}\right]\Delta t\text{\hspace{1em}(16.10)}$$

$$\left\{u_{n+1}\right\}=\left\{u_n\right\}+\left\{{\dot{u}}_n\right\}\Delta t+\left[\left(\frac{1}{2}-\alpha \right)\left\{{\ddot{u}}_n\right\}+\alpha \left\{{\ddot{u}}_{n+1}\right\}\right]\Delta t^2\text{\hspace{1em}(16.11)}$$

其中，$\alpha$ 和 $\delta$ 是 Newmark 积分参数，$\left\{{\ddot{u}}_n\right\}$、$\left\{{\dot{u}}_n\right\}$ 和 $\left\{u_n\right\}$ 分别表示在时间 $t_n$ 时刻的节点加速度、速度和位移向量，而 $\Delta t=t_{n+1}-t_n$ 则是时间步长。利用这三个代数方程（方程 16.9（第 829 页）、方程 16.10（第 829 页）和方程 16.11（第 829 页）），可以针对未知量 $\left\{u_{n+1}\right\}$ 以及三个已知量，构建一个单步时间积分器。

$$\begin{array}{rl}& \left({\alpha }_0[M]+{\alpha }_1[C]+[K]\right)\left\{u_{n+1}\right\}\\ & =\left\{f_{n+1}\right\}+\left[M\right]\left\{{\alpha }_0\left\{u_n\right\}+{\alpha }_2\left\{{\dot{u}}_n\right\}+{\alpha }_3\left\{{\ddot{u}}_n\right\}\right\}+\left[C\right]\left({\alpha }_1\left\{u_n\right\}+{\alpha }_4\left\{{\dot{u}}_n\right\}+{\alpha }_5\left\{{\ddot{u}}_n\right\}\right)\end{array}\text{\hspace{1em}(16.12)}$$

各系数由以下公式给出：${\alpha }_0=1/\alpha \Delta t^2.{\alpha }_1=\delta /\alpha \Delta t.{\alpha }_2=1/\alpha \Delta t.{\alpha }_3=1/\Delta t-1$ . ${\alpha }_4=\delta /\alpha -1$ . 以及 ${\alpha }_5=\Delta t/2\left(\delta /\alpha -2\right)$ .

计算过程如下。首先，利用公式16.12（第830页）计算节点位移向量 $\left\{u_{n+1}\right\}$。接着，程序通过以下关系式计算节点速度向量 $\left\{{\dot{u}}_{n+1}\right\}$ 和节点加速度向量 $\left\{{\ddot{u}}_{n+1}\right\}$：

$$\left\{{\dot{u}}_{n+1}\right\}={\alpha }_1\left(\left\{u_{n+1}\right\}-\left\{u_n\right\}\right)-{\alpha }_4\left\{{\dot{u}}_n\right\}-{\alpha }_5\left\{{\ddot{u}}_n\right\}\text{\hspace{1em}(16.13)}$$

$$\left\{{\ddot{u}}_{n+1}\right\}={\alpha }_0\left(\left\{u_{n+1}\right\}-\left\{u_n\right\}\right)-{\alpha }_2\left\{{\dot{u}}_n\right\}-{\alpha }_3\left\{{\ddot{u}}_n\right\}\text{\hspace{1em}(16.14)}$$

在选择适用于有限元半离散运动方程（方程16.7，第828页）的适当时间积分方案时，最重要的因素是精度、稳定性和耗散性。在条件稳定的时间积分算法中，稳定性受所选时间步长的影响；而在无条件稳定的时间积分算法中，时间步长可以独立于稳定性考虑来选择。

在新马克方法中，数值算法耗散的量可以通过新马克参数之一（$\delta$）来控制，具体如下：

$$\begin{array}{l}\delta \ge \frac{1}{2}\\ \\ \alpha \ge \frac{1}{4}(\frac{1}{2}+\delta {)}^2\end{array}\text{\hspace{1em}(16.15)}$$

(16.15)

当Newmark参数满足上述条件时，Newmark系列方法是无条件稳定的（Hughes [261]（第1072页））。通过引入幅值衰减因子$\gamma$，上述条件可以表述为：

$$\begin{array}{rl}& \delta =\frac{1}{2}+\gamma \\ & \alpha =\frac{1}{4}{\left(1+\gamma \right)}^2\\ & \gamma \ge 0\end{array}\text{\hspace{1em}(16.16)}$$

因此，Fluent 为您提供了 Newmark 积分过程，该过程通过输入幅度衰减因子 $\gamma$ 实现无条件稳定。

16.4.2.2.2 后向欧拉法

后向欧拉法要求速度和加速度按以下方式更新：

$$\left\{{\dot{u}}_{n+1}\right\}=\frac{\left\{u_{n+1}\right\}-\left\{u_n\right\}}{\Delta t}\text{\hspace{1em}(16.17)}$$

$$\left\{{\ddot{u}}_{n+1}\right\}=\frac{\left\{{\dot{u}}_{n+1}\right\}-\left\{{\dot{u}}_n\right\}}{\Delta t}\text{\hspace{1em}(16.18)}$$

结合方程16.17（第830页）和方程16.18（第831页）与方程16.9（第829页）给出的半离散运动方程，可以得到一个关于未知量 $\left\{u_{n+1}\right\}$ 的单步时间积分器，如下所示：

$$\left(\frac{1}{\Delta t^2}\left[M\right]+\frac{1}{\Delta t}\left[C\right]+\left[K\right]\right)\left\{u_{n+1}\right\}=\left\{F_{n+1}^a\right\}+\left[M\right]\left(\frac{\left\{u_n\right\}}{\Delta t^2}+\frac{\left\{{\dot{u}}_n\right\}}{\Delta t}\right)+\left[C\right]\frac{\left\{u_n\right\}}{\Delta t}\text{\hspace{1em}(16.19)}$$

16.4.2.2.3 瑞利阻尼

在结构模型中，瑞利阻尼适用于瞬态本征流固耦合（FSI）案例。如《Newmark方法》（第829页）所述，瞬态运动方程的最终形式可以写为：

$$\begin{array}{rl}& \left({\alpha }_0[M]+{\alpha }_1[C]+[K]\right)\left\{u_{n+1}\right\}\\ & =\left\{f_{n+1}\right\}+\left[M\right]\left\{{\alpha }_0\left\{u_n\right\}+{\alpha }_2\left\{{\dot{u}}_n\right\}+{\alpha }_3\left\{{\ddot{u}}_n\right\}\right\}+\left[C\right]\left({\alpha }_1\left\{u_n\right\}+{\alpha }_4\left\{{\dot{u}}_n\right\}+{\alpha }_5\left\{{\ddot{u}}_n\right\}\right)\end{array}\text{\hspace{1em}(16.20)}$$

其中，各系数由以下公式给出：${\alpha }_0=1/\alpha \Delta t^2$，${\alpha }_1=\delta /\alpha \Delta t$，${\alpha }_2=1/\alpha \Delta t$，${\alpha }_3=1/\Delta t-1$，${\alpha }_4=\delta /\alpha -1$，以及${\alpha }_5=\Delta t/2\left(\delta /\alpha -2\right)$。

此外，结构阻尼矩阵的Rayleigh模型如下：

$$\left[C\right]={\alpha }_d\left[M\right]+{\beta }_d\left[K\right]\text{\hspace{1em}(16.21)}$$

其中，${\alpha }_d$ 是质量比例阻尼系数，而 ${\beta }_d$ 是刚度比例阻尼系数。