趙大為

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

基于有限元法的流體與結(jié)構(gòu)的耦合研究

趙大為

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

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

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

0 引言

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

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

1 流固耦合的分析方法

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

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

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

在隨后的三十多年中，由于計算固體力學(xué)、計算流體力學(xué)的不斷發(fā)展完善以及計算經(jīng)驗的不斷積累，最初的分析方法不再被采用。耦合的算法取代了解耦算法，因為有限差分與有限元聯(lián)合計算有很大困難以及有限元在流場計算中的不斷完善，有限元成為流體—結(jié)構(gòu)耦合的主導(dǎo)算法。

2 標(biāo)準(zhǔn)k-ε湍流模型的封閉方程組

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

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

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

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

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

方程式中的S為脈動流在單位時間內(nèi)通過 Reynolds應(yīng)力從平均流所獲得的能量：

式中， 、 、 、、為經(jīng)驗常數(shù)，對標(biāo)準(zhǔn)模型，Spalding和 Launder[4-5]建議分別取0.09、1.00、1.30、1.44、1.92。

上述方程組成了標(biāo)準(zhǔn)的k-ε湍流模型的封閉方程組。

3 流固耦合系統(tǒng)的方程

ALE描述實際只是與物質(zhì)時間導(dǎo)數(shù)相關(guān)，反映了物理量時間上的絕對變化和相對變化的關(guān)系，從而容易得到N-S方程的ALE描述形式如下：

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

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

4 迭代法求解雙向耦合

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

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

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

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

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

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

5 算例

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

圖1 幾何形狀和入口條件

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

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

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

圖3 不同時刻速度矢量圖

6 結(jié)論

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

[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-，男，工程師。研究方向：船舶系統(tǒng)。