Fluent流固耦合基礎(chǔ)教程（上）

2016-09-04 by:CAE仿真在線來源:互聯(lián)網(wǎng)

最近用fluent UDF做了個彈性梁的流固耦合問題。這里把經(jīng)驗跟大家分享一下。這個帖子主要是介紹基本操作和程序編寫,算是掃盲貼;對于CFD,有限元、梁理論、流固耦合原理不做過多探討,請大家見諒。另外作者水平有限,時間精力也有限,寫得不好也請見諒了。

作者假設(shè)讀者熟悉Fluent,對CFD的基本概念有所了解,具備一定的C/Fortran語言編程技能,并且了解有限元的基本理論和程序編寫。不過即使背景有所欠缺也無妨,點是了解基本流程和操作方法。

1. 流固耦合問題的種類

流固耦合問題可以分為很多種。按照耦合程度分類,可以分為強耦合問題和弱耦合問題。但是這種分類方法有兩個流行的版本。

第一版本根據(jù)流體和固體之間相互影響的程度來劃分,在工業(yè)應(yīng)用中談得比較多。如果固體運動很小,對流場的影響不大,則認(rèn)為是弱耦合。比如金屬管道內(nèi)的水流引起的管壁運動,機翼的振顫等,可以屬于這類。在計算機技術(shù)不夠發(fā)達的時代,這類問題一般采用簡化的流場模型,如簡單幾何形狀的理想流體,解析出由于固體運動而引起的流場作用力,然后再將這些作用力施加到固體控制方程上,從而得到附加質(zhì)量,附加粘度,附加剛度,以及其他非線性項。隨著計算機技術(shù)的發(fā)展,對這類問題的直接數(shù)值模擬逐漸成為可能。如果固體變形或運動比較大,其對流場的影響不可忽略,則屬于強耦合。強耦合問題必須對流場和固體進行細致的計算,充分考慮固體變形對流場的影響。閥門開閉,血液通過血管瓣膜,旗幟在風(fēng)中飄舞,等都屬于這類問題。解決這類問題需要依靠數(shù)值模擬,傳統(tǒng)的簡化模型很難應(yīng)用在這類問題上。

第二個版本根據(jù)流體和固體求解的模式來劃分,在學(xué)術(shù)領(lǐng)域用的比較多。強耦合問題被認(rèn)為是流場和固體變形必須同時求解的問題。帶有滲流的多孔介質(zhì)變形問題可以算這類問題。弱耦合問題則是流場和固體變形可以分開求解,但是二者的信息交換通過循環(huán)迭代完成。多數(shù)工程問題都可以算做這類問題。

這里我要講的是第一類劃分標(biāo)準(zhǔn)中的弱耦合問題,但所用的數(shù)值方法也可以應(yīng)用到強耦合問題中,只是結(jié)構(gòu)體的有限元方法需要采用大變形理論。按照第二種劃分標(biāo)準(zhǔn)來說,這里解決的是流場和固體變形/運動分別求解,兩場之間的相互作用通過迭代實現(xiàn)。

2. 流固耦合問題的數(shù)值模擬

既然是數(shù)值模擬,則需將系統(tǒng)方程在時域上離散。在每個時間步上分別解算流體域和固體域。固體的變形或位移導(dǎo)致流固邊界的運動。這個運動以動網(wǎng)格的形勢傳遞給流場,然后求解流場。流場的解則包括了這個運動所導(dǎo)致的反力。將這個反力反饋給固體,求解其下一時間步的響應(yīng)。如圖1所示。

圖1

以上的概念恐怕大家都清楚,只是談到實現(xiàn)的時候,需要做很多細致的工作。當(dāng)然,如果大家有ANSYS14 WorkBench,則流固耦合問題很容易解決,因為大部分的工作都自動化處理了。但是如果只用Fluent,則需要做一些工作。即便使用先進的數(shù)值模擬工具,如果能弄清楚整個流程以及其中的一些技巧,也會對提高效率以及保證精度有所幫助。下面就講講具體怎么做。

3. 流固耦合例題

既然是具體講,就不如找個例題來說。圖2是一根柔軟的彈性梁,浸泡在在直徑為Dt的圓截面管道里。梁的截面也是圓形,直徑為D,長度為L。梁的兩端為鉸接。整個結(jié)構(gòu)是三維結(jié)構(gòu)。梁可以在三維空間里彎曲。流體為水,梁中點所在截面上的平均流速為U。管道入口處的流速為Uin。上游長度為Lup,下游長度為Ldw。上下游的長度分別為10D。

圖2

4. 流體模型

流體模型的建立是關(guān)鍵的一環(huán),因為流體模型的好壞直接影響到所得解的合理性。建立好的流體模型需要利用流體力學(xué)知識對所處理的問題進行綜合分析。根據(jù)流體的雷諾數(shù)和其它特征建立合適的網(wǎng)格,選擇合理的湍流模型和算法。我們希望流體模型能夠盡量準(zhǔn)確地捕捉到梁表面的力,因此準(zhǔn)確的near-wall處理方法很重要。采用wall function還是精確地模擬邊界層,需要給出適當(dāng)?shù)睦碛?。由于我們處理的是動力問題,希望模型對梁表面力的變化給出準(zhǔn)確的值。采用動態(tài)k-epsilon 或者k-omega湍流模型計算所得的激勵幅值明顯偏小,這是由于雷諾平均造成的。這里我們采用對流體動力特性捕捉比較好的大渦湍流模型(Large Eddy Simulation, LES)。LES對near-wall處理方法要求比較高,需要邊界上的第一層網(wǎng)格厚度滿足y+ ~ 1。在Gambit里可以方便地建立邊界層網(wǎng)格,達到以上要求。如何建立良好的網(wǎng)格是個復(fù)雜的話題,這里且不涉及。

由于固體-液體邊界要移動,流體網(wǎng)格必須改變,需要運用動網(wǎng)格技術(shù)。此時建模的時候需要將動網(wǎng)格區(qū)域和不動網(wǎng)格區(qū)域分開,如圖3所示。

圖3

固體和液體的交界面附近的區(qū)域是動網(wǎng)格區(qū)域(fluid_dyn),遠離固體的區(qū)域可以定義為不動網(wǎng)格區(qū)域。動網(wǎng)格區(qū)域的大小應(yīng)當(dāng)視具體問題而定,但起碼要能夠囊括固體的最大位置變化范圍。梁表面的網(wǎng)格節(jié)點位置隨著梁的運動而變化,因此需要將梁的表面單獨定義為一個區(qū)域 wall_mov_beam。梁的兩個端點也需要進行控制,因此也可以分別定義獨立區(qū)域:wall_mov_end1 和 wall_mov_end2。在劃分流體網(wǎng)格的時候,接近固體表面的地方需要建立邊界層。Fluent幫助文檔上說邊界層應(yīng)該單獨作為一個動網(wǎng)格區(qū)域,這是為了便于控制邊界層的變形,防止出現(xiàn)邊界層的異化。實際上如果將邊界層的厚度增加到物理邊界層厚度的兩倍以上,并且只采用smoothing網(wǎng)格控制,則可以不必將邊界層單獨定義為一個區(qū)域。但是如果變形較大,流體網(wǎng)格需要重新劃分,則最好將邊界層單獨分離為一個區(qū)域,并且在這個區(qū)域內(nèi)使用smoothing控制。對于管壁附近的邊界層,我們并不關(guān)心管壁上的壓力和剪力,而且考慮到出了邊界層以后這個區(qū)域的局部解對整體解影響不是很大,因此可以適當(dāng)降低網(wǎng)格要求,甚至采用wall function近似逼近,這樣可以節(jié)省不少計算量。下面兩張圖分別是圖4幾何模型(部分)和圖5網(wǎng)格劃分(部分)。

圖4

圖5

流體為不可壓縮流體,求解器采用一次隱式瞬態(tài)算法。LES的動網(wǎng)格計算上Fluent 不支持二次格式,這就需要時間步設(shè)定要小一些。時間步長和求解器參數(shù)的選取也是個復(fù)雜的問題,這里也不展開討論了。在這個算例中,我們?nèi)≥^大的步長,以加快計算速度。具體參數(shù)如下:

Time step size: 0.0005 s

Solver settings: unsteady / pressure-based / 1st order implicit

Viscous model: LES / Smagorinsky-Lilly / no dynamic stress

Solution settings:

velocity-pressure coupling: SIMPLE

Relaxation factor: pressure 0.3 / momentum 0.7 / density 1.0 / body force 1.0

Discretization: pressure standard / momentumn bounded central difference

在不加入流固耦合的情況下,計算結(jié)果收斂很好,穩(wěn)定以后每個時間步上的循環(huán)(iteration)為兩步。每個時間步循環(huán)終止時的殘差為:

continuity x-velocity y-velocity z-velocity

4.2946e-04 8.3873e-06 8.1793e-06 1.3922e-04

流體模型解算成功之后就需要考慮動網(wǎng)格和固體變形的問題了。梁表面流體網(wǎng)格節(jié)點的運動需要利用UDF(User Defined Function)來控制。這將在下面的章節(jié)里說明。這里首先設(shè)置動網(wǎng)格的參數(shù)。

動網(wǎng)格區(qū)域fluid_dyn應(yīng)該被設(shè)置為Deforming,采用的網(wǎng)格控制方法為Smoothing,但是對于復(fù)雜結(jié)構(gòu),可能需要Remeshing。具體選項的含義請參見Fluent 幫助文檔中的用戶手冊第11章(Modeling Flows Using Sliding and Deforming Meshes)。這里只說明大體思路。這個區(qū)域的網(wǎng)格將會隨著這個區(qū)域的邊界(梁表面wall_mov_beam)變化。我們只要控制wall_mov_beam,則Fluent可以自動計算出這個區(qū)域內(nèi)的流體網(wǎng)格的新位置。參數(shù)設(shè)置如圖6所示。

圖6

面區(qū)域wall_mov_end 1和wall_mov_end2也設(shè)置為動網(wǎng)格區(qū)域,跟fluid_dyn類似,具體的網(wǎng)格變形由Fluent自動處理。這是由于我們的梁是兩端鉸接的,因此端面上的網(wǎng)格變化不大。對于其他的情況,比如懸臂梁,這兩個面上的點應(yīng)該跟據(jù)梁的運動用UDF來控制,不能交給Fluent自動處理。具體設(shè)置參數(shù)如圖7所示。

圖7

梁表面網(wǎng)格節(jié)點的位移必須通過UDF來控制。因此在選擇的時候應(yīng)當(dāng)選User-Defined。在Motion Attributes里面選擇相應(yīng)的UDF。這一點在下面的章節(jié)中說明。這里先給出Meshing Options的值,參看圖8。這里提醒讀者,wall_mov_beam是驅(qū)動網(wǎng)格變化的“源頭”。我們通過UDF控制wall_mov_beam的運動,而周邊的網(wǎng)格將隨著wall_mov_beam的運動而改變。