基于有限元法的流體與結構的耦合研究

趙大為

（中國艦船研究設計中心， 湖北 武漢430064）

基于有限元法的流體與結構的耦合研究

趙大為

（中國艦船研究設計中心， 湖北 武漢430064）

對有限元法離散流固耦合方程進行了簡要的介紹，給出了粘性流體與彈性結構的全耦合模型。流體采用k-ε湍流模型進行了有限元計算，采用ALE描述實現固體Lagrangian描述到流體Eulerian描述的過渡，使大變形帶來的流體域網格畸變得以避免。對流固耦合方程，采用流固全耦合的迭代算法進行了迭代求解，最后模擬了高雷諾數下湍流流動中沖擊致使渦輪結構旋轉的過程。計算結果表明，流固耦合的算法是可靠的。

有限元法(FEM)；流體沖擊；流固耦合

0 引言

在海洋工程領域中，都會遇到流體和結構的相互作用問題，簡稱為流固耦合問題。流固耦合力學的重要特征是兩相介質間的交互作用，即固體在流體載荷作用下會產生變形或運動，而變形或運動反過來又影響流場的流動，從而改變流體載荷的分布和大小。一般流固耦合的特點是方程的定義域同時有流體域和固體域，未知變量既有流體變量又有固體變量，而且流體域和固體域通常無法單獨求解。

從總體上看，按照耦合機理流固耦合問題可以分為兩大類。第一類是兩相域部分或全部重疊在一起，很難明顯地分開，使描述物理現象的方程，特別是本構方程需要針對具體的物理現象來建立，其耦合效應通過描述問題的微分方程而體現，如滲流問題。第二類問題的特征是耦合作用僅發生在兩相交界面上，在方程上耦合是由兩相耦合面的平衡及協調關系引入的。Zienkiewcz和Bettess[1]曾將上述第二類流固耦合問題分為三種情況，一是流固間有大的相對運動情況，如飛機飛行狀態下的氣動彈性力學問題；二是有限流體運動的短周期情況，如流體受沖擊和水下爆炸問題；三是有限流體位移的長周期情況，如含液容器的流固耦合振動問題。

1 流固耦合的分析方法

流固耦合問題的歷史可以追溯到 1843年，當時Stokes[2]研究了一個無限長圓柱體在無限流體介質中的均勻加速問題。他得出的結論是：流體對圓柱體運動僅有的影響是增加了它的有效質量。Stokes早年的這種引入附加有效質量的基本概念和做法極大地影響了流固耦合分析理論后來的發展。

現代流固耦合分析開始于20世紀50年代。最早進行航天運載器中液體燃料容器的分析研究，而對于核動力系統，最著名的是Fritz和Kiss[3]關于同心圓柱體擺動運動的研究報告，其中流體介質是有限的，固體是剛性體。

自20世紀70年代以來，人們開始研究流體與彈性體間的流固耦合問題。最早的論文是Krajcinovic在1974年發表的。早期的流體-結構耦合分析一般是以核反應堆安全分析為背景的，Belytschko在1977年對這一階段進行了總結[4]，指出當時的耦合分析實際是以解耦的方式進行的，即：首先，忽略結構變形，將結構看作流體的剛性壁，計算流場分布；然后，將流體計算得到的剛性壁壓力加在結構上，計算結構的變形、響應。一般認為，流場作用于剛性壁的壓力比柔性壁要大，所以結果對結構來說是趨于保守的。這種解耦算法一方面導致結構的過于安全的設計；另一方面，當流體—結構系統發生共振時，這種算法是不可信、不安全的。正是因為采用解耦的算法，流體域可以采用有限元，也可以采用有限差分。

在隨后的三十多年中，由于計算固體力學、計算流體力學的不斷發展完善以及計算經驗的不斷積累，最初的分析方法不再被采用。耦合的算法取代了解耦算法，因為有限差分與有限元聯合計算有很大困難以及有限元在流場計算中的不斷完善，有限元成為流體—結構耦合的主導算法。

2 標準k-ε湍流模型的封閉方程組

Reynolds首先認為流體的運動在從層流狀態過渡到湍流狀態后，流體的物理和力學性質沒有變化，流體的連續性沒有受到破外，Navier-Stokes方程仍能描寫湍流的瞬時規律。上述假定雖然至今仍未能得到嚴密論證和證實，但一百多年以來，研究湍流的實踐表明，這些假定并沒有與實際情況發生矛盾，說明Navier-Stokes方程應用于湍流問題是適宜的。如用f表示流動變量的瞬時值，表示f的統計平均值，各量的瞬時值與統計平均值之差即為各量大的脈動值 ′，即：

根據式(1)，將連續方程、動量方程平均化后得(為符號簡單起見，平均值速度u、壓力p符號上的“-”都去掉了)。

方程(3)首先由Reynolds于1895年導出，因而通常稱為 Reynolds方程。其中稱為Reynolds應力。 方程(2)、(3)中只有緩變的平均值，故消除了因脈動而引起的不規則變化，使方程便于應用。由于方程中出現了許多新項，又產生了新的困難。方程中未變量的個數多于方程個數，因而是不封閉的，必須在試驗資料基礎上做出一些假設，建立補充關系后才能求解。對雷諾應力采用與層流類似的模型，可將其表示為：

式中， 為湍流粘性；k為湍流脈動動能的平均值；在式(4)中引進了兩個未知量 和k，為了封閉還要求兩個方程。對于大Reynolds數流動，由量綱分析可得：

式中，ε為脈動動能的耗損速率。關于k和ε的精確方程可有Navier-Stokes方程和Reynolds方程導得，對其中的各項進行分析和某些假設(模式化)后，可以得出能量輸運方程：

方程式中的S為脈動流在單位時間內通過 Reynolds應力從平均流所獲得的能量：

式中， 、 、 、、為經驗常數，對標準模型，Spalding和 Launder[4-5]建議分別取0.09、1.00、1.30、1.44、1.92。

上述方程組成了標準的k-ε湍流模型的封閉方程組。

3 流固耦合系統的方程

ALE描述實際只是與物質時間導數相關，反映了物理量時間上的絕對變化和相對變化的關系，從而容易得到N-S方程的ALE描述形式如下：

從式(10)可以看出，ALE描述并不影響空間導數項。針對邊界問題，相當于調整網格后，使對原問題的描述既不是Lagrangian描述，也不是Eulerian描述，調整網格引入的對流量 也必須考慮在內。

根據公式(10)，ALE描述的k-ε湍流模型封閉方程組可寫為：

4 迭代法求解雙向耦合

在很多耦合問題中，流體的作用力影響結構的變形，同時結構的位移又影響流場的形態。這正是進行流固耦合分析的原因，這種類型的分析叫做“雙向耦合(two-way coupling)”。在某些情況下，結構的變形非常小，它對流體的影響可以忽略。只有流體的應力需要施加到結構上，流體和結構模型之間不需要迭代，這種類型的耦合叫做“單項耦合(one-way coupling)”。

迭代法求解雙向耦合的解法也叫做分離法。流體和結構的求解變量是完全耦合的。流體方程和結構方程是按順序相互迭代求解的，各自在每一步得到的結果提供給另一部份使用，直到耦合系統的解達到收斂，迭代停止。

計算過程可以概括為：為了得到t+Δt時刻的解，在流體模型和結構模型之間開始迭代計算。設初始解為：。對迭代步k=1，2，··，進行下面的求解過程以得到(t+Δt)X。

2)如果只需要滿足應力收斂條件，則只要計算應力殘量并與迭代容差相比較。如果這個標準滿足了，以下3)到5)步驟就可以忽略了。

5)如果只需要滿足位移收斂條件，則要計算位移殘量并與迭代容差相比較。如果應力和位移的標準都要求滿足，則兩個收斂條件都要檢查。如果迭代不收斂，回到第1)步繼續下一個迭代，直到達到FSI迭代的最大數。

6)保存并輸出流體和結構的結果。在這種求解方法中，時間步和求解時間是由流體模型控制的。然而，結構模型中定義的所有時間函數必須覆蓋計算的時間范圍。耦合系統最終控制收斂的參數是由流體模型決定的。這些參數包括應力和位移迭代容差，松弛因子，收斂標準等。流體和結構模型的收斂分別由各自的方程控制。保存和輸出解也同樣由各自的模型控制。

5 算例

此算例模擬在湍流流體作用下致使結構旋轉，流體模型仍然為k-ε湍流模型。如圖1所示，開始分析時，渦輪是靜止的。在渦輪入口處突然施加一個入口邊界的壓力，流體流過渦輪箱，使渦輪旋轉。

圖1 幾何形狀和入口條件

各參數取值為：①流體：密度ρ=1000kg/m3，運動粘性系數ν=0.01m2/s，初始壓力為P=500Pa。②渦輪：密度ρ=7800kg/m3，彈性模量E=207GPa，泊松比ν′=0.3。

采用滑移網格，時間步長 Δt=0.2s。本問題計算所得的流場壓力等值域和速度矢量見圖2、3，圖2中同時給出了渦輪大幅轉動以及應力變化。從圖中可以清楚地看到，在施加突然的壓力之后，隨著渦輪內部流場的變化，渦輪葉片產生了明顯的應力變化和轉動。

圖2 不同時刻流體節點壓力云圖及渦輪有效應力云圖(Pa)

圖3 不同時刻速度矢量圖

6 結論

上述算例分析了渦輪結構在來流壓力作用下的旋轉，結構在受來流的沖擊開始旋轉，從計算結果可以看出，流體的壓力能夠讓流體域中的結構產生運動和變形，渦輪葉片的受力也極不均勻，轉動軸心處的應力變化較大。

[1] Zienkiewicz O C, Bettess P. Fluid-structure dynamic interaction and wave forces: an introduction to numerical treatment[J]. Int. J. Num. Meth. Eng.,1978(13): 1-16.

[2] Karamcheti K. Principle of idea-fluid aerodynamics[M]. New York: John Wiley & Sons,1966.

[3] Fritz R J. The effect of liquids on the dynamic motion of solid [J]. J. Eng. Ind., 1972, 94(1):167-173.

[4] Belytschko T. Method and programs for analysis of fluid-structure systems [J]. Nucl. Eng. Des. ,1977,42:41-52.

[5] Launder B E, Spalding B E. Lectures in mathematical models of turbulence[M]. London: Academic Press,1972: 78-115.

Coupling Research of Fluid and Structure Based on Finite Element Method

ZHAO Da-wei

(China Ship Development and Design Center, Wuhan 430064, China)

The process of the fluid-solid coupling equation discretized by finite element method is simply introduced. It also gives the full coupled model of viscous fluid and elastic structure. The finite element calculation of the fluid is carried out by using k-ε turbulence model. Using ALE description achieves the translation from the lagrangin description of solid to the Eulerian description of the fluid to avoid the mesh distortion of the fluid grid caused by the great deformation. Fluid solid coupling equation iterative solution is carried out by using iterative algorithm fluid-solid coupling. Finally, the rotating process of turbine structure caused by turbulent flow impact in high Reynolds number is simulated. The calculation results show that the fluid solid coupling method is reliable.

finite element method (FEM);fluid impact;fluid-solid coupling (FSI)

U664.5+8

A

趙大為(1982-，男，工程師。研究方向：船舶系統。