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積分法無條件穩定的條件是