趙大為

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

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

趙大為

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

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

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

0 引言

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

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

1 流固耦合的分析方法

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

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

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

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

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

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

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

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

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

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

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

式中， 、 、 、、為經(jīng)驗(yàn)常數(shù)，對(duì)標(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描述實(shí)際只是與物質(zhì)時(shí)間導(dǎo)數(shù)相關(guān)，反映了物理量時(shí)間上的絕對(duì)變化和相對(duì)變化的關(guān)系，從而容易得到N-S方程的ALE描述形式如下：

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

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

4 迭代法求解雙向耦合

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

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

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

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

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

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

5 算例

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

圖1 幾何形狀和入口條件

各參數(shù)取值為：①流體：密度ρ=1000kg/m3，運(yùn)動(dòng)粘性系數(shù)ν=0.01m2/s，初始?jí)毫镻=500Pa。②渦輪：密度ρ=7800kg/m3，彈性模量E=207GPa，泊松比ν′=0.3。

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

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

圖3 不同時(shí)刻速度矢量圖

6 結(jié)論

上述算例分析了渦輪結(jié)構(gòu)在來(lái)流壓力作用下的旋轉(zhuǎn)，結(jié)構(gòu)在受來(lái)流的沖擊開始旋轉(zhuǎn)，從計(jì)算結(jié)果可以看出，流體的壓力能夠讓流體域中的結(jié)構(gòu)產(chǎn)生運(yùn)動(dòng)和變形，渦輪葉片的受力也極不均勻，轉(zhuǎn)動(dòng)軸心處的應(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)。