2.2 波动方程的有限元解法

有限元法是在差分法和变分法的基础上发展起来的一种数值方法，它吸取了差分法对求解域进行离散处理的启示，又继承了里兹法选择试探函数的合理方法。从实质上看，有限元法与里兹法是等效的，它属于里兹法的范畴，多数问题的有限元方程都是利用变分原理来建立的。但由于有限元法采用了离散处理，能够处理更复杂的问题，并且计算更为简单，因而应用范围更为广泛。对于不同物理性质和数学模型的问题，有限元法求解的基本步骤是相同的，一般典型的分析步骤为：①将求解区域离散化；进行单元划分；②单元分析；③集合成整体；④联立方程组求解。对于不同的物理模型，采用单元类型不尽相同，但分析方法是相同的。采用有限元法分析波动问题的过程^[138-140]如下：

1）根据虚功原理，推导波动方程的变分形式。

2）对求解区域进行离散化、单元分析、整体分析，推导与变分问题等价的线性方程组表达形式。

3）求解线性方程组，将所求得的解作为波动方程的解。

2.2.1 波动方程的变分问题

直接根据边值条件求解波动方程一般很难找到解析解。下面考虑边界条件及初始条件，推导与波动方程等价的变分方程，再导出有限元法的计算公式。弹性固体介质（各向异性或各向同性）中波动方程的初始边界条件为：

边界条件

初始条件

式中，n_i是介质求解区域Ω表面S的法向分量；T_i、f_i表示作用在介质上的面力和体力分量；为质点位移对时间的一阶导数。

用位移分量的变分δu_i乘以式（2-1）的两端并在整个求解区域Ω上积分

由于，同时根据格林公式以及边界条件式（2-9），式（2-11）可写为

式（2-12）就是满足初始边界条件式（2-9）、式（2-10）的波动方程的变分方程，其中Q_k、n_k分别为向量Q、S的法向量n的x_k向分量。只要Ω连续或分片连续，就能满足式（2-12）可积的要求。

2.2.2 波动方程的有限元公式

超声检测过程中，探头向介质发射超声波，从而使机械振动在介质内部质点间由近及远相继连续发生而形成应力波的传播。用有限元法分析超声波传播过程中固体介质内部的应力、变形及运动状态等，需要将固体介质离散成若干个单元，进行单元分析和整体分析，得到有限元计算公式。

因为质点位移是空间坐标的函数，所以位移分量u_i可以表示为

u_i（x₁，x₂，x₃）=N^TU_i=NU^T_i （2-13）

N_x1、N_x2、N_x3为x₁、x₂、x₃方向的节点数，N为形函数向量。

N=N（x₃）ⓍN（x₂）ⓍN（x₁）

N（x_i）=（N₁（x_i），N₂（x_i），…，N_Nxi（x_i））T （2-14）

式中，符号⊗表示矩阵的直积。当采用三线性插值函数时，形函数可表示为

将式（2-13）代入式（2-12）得

以矩阵形式改写式（2-16）得

式中，

式中，、Kgh分别为分块矩阵。

由于式（2-17）中的δU是任意的，因此可以推导出波动方程的有限元方程

式中，、K分别为系统的质量矩阵、刚度矩阵；F为节点载荷向量，超声检测中，F代表探头的激励载荷。

由以上推导过程可以看出式（2-18）与波动方程式（2-1）是等价的，很显然，通过离散得到的有限元方程式（2-18）要比波动方程式（2-1）容易求解，因此可以通过求解有限元方程得到波动方程的解。

2.2.3 波动方程有限元公式的解法及其稳定性

1.波动方程有限元公式求解

目前求解式（2-18）的方法有直接法和振动叠加法两类，直接积分法又包括中心差分法和NewmarK积分法，其中NewmarK积分法是著名的有限元分析软件ANSYS求解动力学问题的方法，利用NewmarK积分法求解波动方程的算法如下。

（1）初始计算

1）通过对求解区域离散化、单元分析、整体分析，把变分问题近似表达成线性方程组，从而形成系统的刚度矩阵K、质量矩阵和阻尼矩阵C（如果忽略阻尼的影响，则C=0）。

2）给定。

3）选择时间步长Δt、参数α和δ，并计算积分常数。

4）形成有效刚度矩阵元：

5）三角分解元：

（2）对于每一时间步长，计算t+Δt时刻作用在系统上的有效载荷

求解t+Δt时刻的位移

求解t+Δt时刻的加速度和速度

根据式（2-22）求得位移向量U（t），进而可求得相应的应变ε（t）和应力σ（t）。

2.波动方程有限元公式求解稳定性探讨

式（2-18）可用与其等价的n个形式相同、不相耦合的微分方程组表示，因此只需讨论该方程组中任意一个方程的稳定性即可，该方程可表示为

由于稳定性讨论的是初始条件或边界条件做微小变动时，方程解的变动情况，因此可令f_i=0，式（2-23）简化为

假设时间步长为Δt，采用NewmarK积分法求解式（2-18），得到如下表达式

将式（2-25）代入式（2-24），可以得到

利用NewmarK积分法，根据式（2-24），将式（2-26）改写为

利用式（2-27）和式（2-24），可将式（2-26）改写为

（2-28）

式中

p_i=Δt²ω²_i （2-29）

假设方程（2-28）的解具有如下形式

（u_i）t+Δt=λ（u_i）t，（u_i）t=λ（u_i）t-Δ_t （2-30）

代入式（2-28），得到关于λ的特征方程

该方程（2-31）的特征值为

其中

方程的解稳定必须满足两个条件：

1）小阻尼情况下具有振荡的特性，即特征根λ_1，2应当为复数，即

4（1+h）²>（2-g）² （2-34）

将式（2-33）代入式（2-34）得

要使式（2-35）在p_i很大时仍然能成立，则

2）绝对值应当是有限大小的，即，要使式（2-35）在p_i很大时仍然能成立，则

δ≥0.5

当式（2-36）成立时，式（2-37）恒成立，因此综合以上分析过程，NewmarK积分法无条件稳定的条件是