充流管道單向流固耦合數(shù)值模擬自動(dòng)化研究

陳堅(jiān)紅，周天情，盛德仁，李 蔚，陳 輝

（1.浙江大學(xué) 熱工與動(dòng)力系統(tǒng)研究所，杭州310027；2.山東電力工程咨詢?cè)河邢薰?，?jì)南250013）

充流管道系統(tǒng)在電力、石油、化工等各類工業(yè)裝置中有著廣泛的應(yīng)用，充流管道的流固耦合（FSI）振動(dòng)是當(dāng)前理論和工程研究的熱點(diǎn)［1-4］.流體在管道中流動(dòng)會(huì)對(duì)管道產(chǎn)生一定的作用力，使管道產(chǎn)生一定的變形，而管道的變形會(huì)促使流體的流場(chǎng)發(fā)生改變，變化后的流場(chǎng)又給管道結(jié)構(gòu)帶來(lái)新的激勵(lì)，這種相互作用的過(guò)程就是流固耦合過(guò)程［5］.近幾十年來(lái)，國(guó)內(nèi)外學(xué)者對(duì)流固耦合的理論和計(jì)算方法開展了廣泛的研究，取得了一些成果.但是，由于流固耦合問(wèn)題存在復(fù)雜性，流固耦合方程的建立還存在許多假設(shè)，遠(yuǎn)遠(yuǎn)不能達(dá)到解決實(shí)際工程問(wèn)題的需要，因此有必要進(jìn)行深入研究.對(duì)于流固耦合過(guò)程的求解，比較簡(jiǎn)單的問(wèn)題可以采用解析法和半解析法，而具有復(fù)雜邊界條件的實(shí)際工程問(wèn)題則很難給出解析解，需要采用數(shù)值計(jì)算方法進(jìn)行反復(fù)迭代計(jì)算才能得到分析結(jié)果［6］.由于采用數(shù)值計(jì)算方法計(jì)算一個(gè)完整的流固耦合過(guò)程時(shí)，計(jì)算量比較大，且實(shí)際工程中大多數(shù)充流管道的計(jì)算只需考慮流體對(duì)管道結(jié)構(gòu)的作用——單向流固耦合過(guò)程，因此，單向流固耦合過(guò)程的求解對(duì)充流管道的計(jì)算有重要的現(xiàn)實(shí)意義.

許多學(xué)者利用數(shù)值模擬方法對(duì)單向流固耦合過(guò)程進(jìn)行了求解，并取得了一定成效.偶國(guó)富［7］運(yùn)用Ansys有限元分析軟件，采用物理環(huán)境順序偶合法進(jìn)行單向流固耦合；楊振國(guó)等［8-9］采用Ansys結(jié)構(gòu)分析軟件與流體分析軟件CFX 進(jìn)行單向流固耦合的模擬；楊濤等［10-13］將Fluent流場(chǎng)分析軟件與Ansys結(jié)構(gòu)分析軟件相結(jié)合，用線性插值、曲線擬合以及取近似值等方法實(shí)現(xiàn)了流場(chǎng)數(shù)據(jù)與結(jié)構(gòu)數(shù)據(jù)的傳遞.然而，以上研究均為對(duì)某個(gè)具體問(wèn)題而進(jìn)行的單向流固耦合數(shù)值模擬，模擬時(shí)要逐步進(jìn)行軟件操作和數(shù)據(jù)處理，過(guò)程極為復(fù)雜，而且在流體計(jì)算與結(jié)構(gòu)計(jì)算之間的數(shù)據(jù)傳遞上也采取了近似算法，計(jì)算的準(zhǔn)確性有待檢驗(yàn).為了解決以上單向流固耦合數(shù)值模擬中操作過(guò)程復(fù)雜、數(shù)據(jù)傳遞失真的問(wèn)題，筆者采用C＋＋語(yǔ)言編寫充流管道流固耦合數(shù)值模擬程序，對(duì)Gambit軟件、Fluent軟件和Ansys軟件進(jìn)行封裝，實(shí)現(xiàn)對(duì)一個(gè)三維充流管道進(jìn)行單向流固耦合數(shù)值模擬的自動(dòng)化計(jì)算分析，程序的計(jì)算流程為：先在Gambit軟件中對(duì)管道和流體進(jìn)行建模，并劃分好網(wǎng)格，然后導(dǎo)入到Fluent軟件中進(jìn)行流體數(shù)值模擬，將經(jīng)過(guò)流體計(jì)算后的管道網(wǎng)格模型連同其所受載荷一起導(dǎo)出，程序?qū)?dǎo)出的數(shù)據(jù)進(jìn)行相應(yīng)整理，并用APDL 參數(shù)化語(yǔ)言編制出Ansys軟件的命令流文件.最后程序調(diào)用Ansys軟件使其自動(dòng)讀取APDL命令流文件，無(wú)失真地加載經(jīng)過(guò)流體計(jì)算后的管道網(wǎng)格模型連同其所受載荷，并進(jìn)行結(jié)構(gòu)分析.

1 數(shù)值計(jì)算方法及程序?qū)崿F(xiàn)

1.1 單向流固耦合的數(shù)值計(jì)算方法

管道中流體的流動(dòng)滿足質(zhì)量守恒定律、動(dòng)量守恒定律和能量守恒定律，這些規(guī)律體現(xiàn)在流體力學(xué)中就是相應(yīng)的連續(xù)性方程和流體運(yùn)動(dòng)方程（即N-S方程）［14］.

按照質(zhì)量守恒定律，對(duì)于某一控制體內(nèi)的流體，其連續(xù)性方程的積分形式為

式中：Vol表示控制體；A表示控制面.根據(jù)奧-高公式，在直角坐標(biāo)系中，對(duì)于不可壓縮的均質(zhì)流體，可將上式化為以下微分形式：

對(duì)于密度和黏性系數(shù)為常數(shù)的流體，其運(yùn)動(dòng)方程（N-S方程）為以下矢量形式：

N-S方程比較準(zhǔn)確地描述了實(shí)際的流動(dòng).黏性流體的流動(dòng)分析可歸結(jié)為對(duì)此方程的研究.

管道結(jié)構(gòu)可以簡(jiǎn)化為線性振動(dòng)系統(tǒng).對(duì)于一個(gè)自由度為n的線性管道振動(dòng)系統(tǒng)，假設(shè)系統(tǒng)中各質(zhì)量、剛度和阻尼已知，根據(jù)牛頓第二定律以及Mykles-tad法，該系統(tǒng)的運(yùn)動(dòng)微分方程為

式中：M、C、K分別為n階對(duì)稱正定質(zhì)量矩陣、n階阻尼對(duì)稱方陣和n階剛度對(duì)稱方陣；F（t）為管道所受外力.

求解彈性機(jī)翼氣動(dòng)力的基本思想有弱耦合法（loose coupling method）和強(qiáng)耦合法（tight coupling method）［15］.同樣，求解充流管道的流固完耦合方法也分為弱耦合法和強(qiáng)耦合法.

弱耦合法即傳統(tǒng)的耦合流體分析模式和結(jié)構(gòu)分析模式的方法，首先完成流體分析，收斂的流體分析結(jié)果轉(zhuǎn)移到結(jié)構(gòu)分析中，用有限元方法計(jì)算結(jié)構(gòu)變形之后，針對(duì)變形的結(jié)構(gòu)再進(jìn)行流體分析，重復(fù)此過(guò)程直至結(jié)構(gòu)變形收斂；強(qiáng)耦合方法中，流體方程和結(jié)構(gòu)方程是同時(shí)求解的，即在流體方程求解迭代期間，間斷地按照還未收斂的流體來(lái)計(jì)算結(jié)構(gòu)變形，再把變形量代入流體計(jì)算的迭代過(guò)程中，直到結(jié)構(gòu)變形和流體力都收斂.

對(duì)于弱耦合法，只須將每次流體計(jì)算收斂后的流體力F（t）代入式（4）中即可進(jìn)行結(jié)構(gòu)分析，即完成了一步單向弱流固耦合計(jì)算.

對(duì)于強(qiáng)流固耦合，流體運(yùn)動(dòng)在結(jié)構(gòu)上產(chǎn)生的作用力須轉(zhuǎn)換為以下形式：

式中：Ma、Ca、Ka分別為流體作用對(duì)結(jié)構(gòu)產(chǎn)生的附加n階對(duì)稱正定質(zhì)量矩陣、附加n階阻尼對(duì)稱方陣和附加n階剛度對(duì)稱方陣，其均為流動(dòng)狀態(tài)的函數(shù)，隨流體流動(dòng)狀態(tài)的變化而變化［5］.

將式（5）代入式（4）中，得

每進(jìn)行一次式（6）的計(jì)算就完成了一步強(qiáng)流固耦合計(jì)算.

由于一方面強(qiáng)耦合算法計(jì)算量太大，且目前還處于基本理論的研究階段；另一方面弱耦合算法能充分利用現(xiàn)有計(jì)算流體力學(xué)和計(jì)算結(jié)構(gòu)力學(xué)的方法和程序，只要增加數(shù)據(jù)交換模塊即可［16］.因此，目前大多數(shù)耦合計(jì)算都建立在弱流固耦合的基礎(chǔ)之上.

本計(jì)算程序利用Fluent軟件和Ansys軟件來(lái)進(jìn)行充流管道單向弱流固耦合的計(jì)算.對(duì)于內(nèi)部充滿流體的管道，管內(nèi)流體的流動(dòng)狀態(tài)可以根據(jù)式（1）和式（3）由Fluent軟件進(jìn)行模擬計(jì)算，然后將流體計(jì)算結(jié)果加載到管道內(nèi)壁上，進(jìn)而根據(jù)式（4）由Ansys軟件進(jìn)行管道結(jié)構(gòu)分析.

1.2 程序?qū)崿F(xiàn)方法

Gambit軟件和Fluent軟件在操作過(guò)程中會(huì)生成相對(duì)應(yīng)的操作命令流文件，這些操作命令流文件能直接被軟件讀取，只要將事先編制好的操作命令流文件導(dǎo)入到軟件中，便可實(shí)現(xiàn)軟件的自動(dòng)化運(yùn)行.而Ansys軟件提供了一種運(yùn)行批處理命令的APDL 格式，只要事先編制好APDL 代碼并導(dǎo)入到Ansys中運(yùn)行，即可實(shí)現(xiàn)Ansys軟件的自動(dòng)化分析.基于上述思路，用C＋＋語(yǔ)言編制出一種用戶界面和計(jì)算過(guò)程控制程序，該程序具有以下功能：（1）管道模型參數(shù)和流體狀態(tài)參數(shù)通過(guò)程序界面上的對(duì)話框輸入；（2）程序根據(jù)輸入的各類參數(shù)自動(dòng)生成相應(yīng)的Gambit操作命令流文件、Fluent操作命令流文件及Ansys的APDL 命令流文件；（3）程序能自動(dòng)將Fluent軟件的計(jì)算結(jié)果導(dǎo)出，并無(wú)失真地轉(zhuǎn)換成Ansys結(jié)構(gòu)分析所需的各類文件，供Ansys軟件調(diào)用；（4）程序能先后自動(dòng)控制調(diào)用Gambit軟件、Fluent軟件和Ansys軟件并執(zhí)行對(duì)應(yīng)的命令流文件；（5）計(jì)算完成后，程序能清晰地顯示出所需的流體分析結(jié)果和結(jié)構(gòu)分析結(jié)果.程序的運(yùn)行流程如圖1所示.

圖1 程序運(yùn)行流程圖Fig.1 Flow chart of the program

2 充流管道數(shù)值模擬分析

2.1 程序計(jì)算結(jié)果驗(yàn)證

為了驗(yàn)證該程序計(jì)算方法的可行性與計(jì)算結(jié)果的準(zhǔn)確性，對(duì)與參考文獻(xiàn)［17］中相類似的充流管道模型進(jìn)行數(shù)值計(jì)算，并將計(jì)算結(jié)果與試驗(yàn)結(jié)果進(jìn)行對(duì)比分析.

文獻(xiàn)［17］中作為對(duì)比的全鋼充流直管段長(zhǎng)為3.81m，內(nèi)徑為0.08m，管壁厚度為0.004 5m，鋼的楊氏彈性模量為21.06×1010Pa，泊松比為0.3，管道內(nèi)液體密度為1 000kg／m3，流速為8.5 m／s.文獻(xiàn)中所進(jìn)行的測(cè)試為模態(tài)測(cè)試.按照?qǐng)D1所示的流程，進(jìn)行相關(guān)參數(shù)的輸入，選擇結(jié)構(gòu)分析類型為“模態(tài)分析”.在Fluent軟件中進(jìn)行流體計(jì)算后可得到管道內(nèi)壁受到的壓力值，將該值導(dǎo)入到Ansys軟件中可得圖2所示的壓力云圖.

通過(guò)Ansys軟件中的模態(tài)分析，可得到該段充流管道的各模態(tài)頻率，對(duì)文獻(xiàn)［17］中圖6的頻響曲線進(jìn)行統(tǒng)計(jì)，也可得出各階頻率值.將本程序的計(jì)算結(jié)果與文獻(xiàn)［17］中的試驗(yàn)結(jié)果進(jìn)行對(duì)比，可得到各階頻率所對(duì)應(yīng)的誤差率，如表1所示.

圖2 全鋼充流直管的壓力云圖（單位：Pa）Fig.2 Pressure distribution of the fluid-filled steel pipeline（unit：Pa）

由表1的數(shù)據(jù)可以看出，本程序計(jì)算出的各階頻率值與文獻(xiàn)［17］中模態(tài)測(cè)試所得的各階頻率值非常接近.因此，可以認(rèn)為采用該計(jì)算程序?qū)Τ淞鞴艿肋M(jìn)行單向流固耦合數(shù)值模擬是可行的，且計(jì)算結(jié)果是準(zhǔn)確的.

表1 充流管道模態(tài)頻率Tab.1 Modal frequency of the fluid-filled pipeline

2.2 充流管道的應(yīng)力分析

在工程管道設(shè)計(jì)中，常用的管道應(yīng)力分析軟件如CAESAR II、GLIF等并不能準(zhǔn)確地將管道內(nèi)部流體的壓力分布施加到管道內(nèi)壁上，從而導(dǎo)致計(jì)算出的管道應(yīng)力結(jié)果不夠準(zhǔn)確.該計(jì)算程序能夠耦合流體計(jì)算和結(jié)構(gòu)計(jì)算，準(zhǔn)確獲知管道內(nèi)壁的壓力分布，使得管道的計(jì)算應(yīng)力值更加準(zhǔn)確.

以某火力發(fā)電廠1 000MW 超超臨界機(jī)組6號(hào)低壓加熱器的疏水管道中兩個(gè)支吊架之間的一段為例進(jìn)行計(jì)算.該管段由半徑相等、兩兩相互垂直的三段直管通過(guò)兩段90°彎頭相連接而成，其內(nèi)徑為0.257m，壁厚為0.008 m，管道內(nèi)壁粗糙度為0.046mm，三段直管的長(zhǎng)度依次為1.8 m、3.7 m和2.0m，兩端彎管的彎曲半徑均為0.381m.

根據(jù)以上結(jié)構(gòu)尺寸，在程序中輸入相應(yīng)的參數(shù)，程序自動(dòng)調(diào)用Gambit軟件進(jìn)行模型建立、網(wǎng)格劃分及邊界的定義，結(jié)構(gòu)及網(wǎng)格劃分如圖3所示.

圖3 疏水管道模型圖Fig.3 Modelling of the drain pipeline

通過(guò)程序自動(dòng)啟動(dòng)Fluent軟件，并在相應(yīng)的數(shù)據(jù)輸入對(duì)話框中輸入疏水流動(dòng)參數(shù)：入口壓力為0.229 MPa、入口溫度為124.6 ℃、入口的流速為1.273m／s.然后Fluent軟件自動(dòng)進(jìn)行流體分析，并將計(jì)算結(jié)果自動(dòng)保存.圖4為管道壁面所受的壓力云圖.

圖4 疏水管道內(nèi)部流體壓力云圖Fig.4 Fluid pressure distribution in the drain pipeline

Fluent計(jì)算完成后，操作程序?qū)luent軟件導(dǎo)出的計(jì)算結(jié)果進(jìn)行自動(dòng)處理，以便無(wú)失真地導(dǎo)入到Ansys軟件中.然后操作程序自動(dòng)啟動(dòng)調(diào)用Ansys軟件、讀入處理過(guò)的結(jié)果文件并自動(dòng)進(jìn)行靜力分析.在靜力分析過(guò)程中，施加在管段上的力包括重力和流體壓力，管段受到的約束包括起始端的剛性固定和末尾段的X方向限位約束.

經(jīng)過(guò)上述流體分析和結(jié)構(gòu)分析可得，該管段所受到的最小應(yīng)力為0.148 MPa，最大應(yīng)力為24.8 MPa，應(yīng)力云圖如圖5所示.

通過(guò)以上實(shí)例可知，該計(jì)算程序完全能夠模擬充流管道內(nèi)部流體壓力分布，并根據(jù)流體壓力進(jìn)行管道應(yīng)力計(jì)算，它可以用于異徑管、周期管等特殊管道的應(yīng)力計(jì)算，從而彌補(bǔ)了通用管道應(yīng)力分析軟件的不足.在閥門、大小頭以及噴嘴等管件的設(shè)計(jì)改造中，該計(jì)算程序也可以用來(lái)進(jìn)行流體流場(chǎng)模擬和結(jié)構(gòu)應(yīng)力計(jì)算，以獲得更好的性能和使用壽命.另外，該計(jì)算程序還可用于對(duì)正在使用中的充液管道進(jìn)行風(fēng)險(xiǎn)預(yù)測(cè)，找出管道中可能出現(xiàn)故障的位置，及時(shí)采取相應(yīng)的安全措施.

圖5 疏水管道應(yīng)力圖（單位：Pa）Fig.5 Stress distribution of the drain pipeline（unit：Pa）

2.3 程序的進(jìn)一步發(fā)展

程序能夠方便地實(shí)現(xiàn)從流體計(jì)算到結(jié)構(gòu)分析的單向流固耦合數(shù)值模擬，在此基礎(chǔ)上，只要稍加改進(jìn)就可以完成完整的流固耦合過(guò)程的數(shù)值分析.具體的改進(jìn)方法如下：獲取結(jié)構(gòu)計(jì)算后各節(jié)點(diǎn)的位移信息，對(duì)變形后的充流管道進(jìn)行重新建模，然后在更新后的模型中再進(jìn)行流體計(jì)算和結(jié)構(gòu)分析，如此往復(fù)循環(huán)，直至最終的管道結(jié)構(gòu)分析計(jì)算結(jié)果收斂.

3 結(jié) 論

采用C＋＋語(yǔ)言編程技術(shù)對(duì)Gambit軟件、Fluent軟件和Ansys軟件進(jìn)行封裝，編制出一個(gè)可對(duì)三維充流管道進(jìn)行單向流固耦合數(shù)值模擬的程序，實(shí)現(xiàn)了從模型建立、流體計(jì)算到結(jié)構(gòu)分析過(guò)程的自動(dòng)化，確保了從流體計(jì)算到結(jié)構(gòu)計(jì)算過(guò)程數(shù)據(jù)傳遞的準(zhǔn)確性.通過(guò)與參考文獻(xiàn)試驗(yàn)結(jié)果的對(duì)比，證明該計(jì)算程序?qū)Τ淞鞴艿肋M(jìn)行單向流固耦合數(shù)值模擬是可行的，且計(jì)算結(jié)果是準(zhǔn)確的.實(shí)例計(jì)算表明，該程序在充流管道單向流固耦合數(shù)值模擬中具有操作簡(jiǎn)便、結(jié)果準(zhǔn)確的特點(diǎn)，可應(yīng)用于多種充流管道的設(shè)計(jì)改造和風(fēng)險(xiǎn)預(yù)測(cè)中.

［1］PITTARD M T，EVANS R P，MAYNES R D，et al.Experimental and numerical investigation of turbulent flow induced pipe vibration in fully developed flow［J］.Review of Scientific Instruments，2004，75（7）：2393-2401.

［2］TIJSSELINGA A S，VARDY A E.Fluid-structure interaction and transient cavitation tests in a T-piece pipe［J］.Journal of Fluids and Structures，2005，20（6）：753-762.

［3］YANG Ke，LI Qiusheng，ZHANG Lixiang.Longitudinal vibration analysis of multi-span liquid-filled pipelines with rigid constraints［J］.Journal of Sound and Vibration，2004，273（1／2）：125-147.

［4］LI Qiusheng，YANG Ke，ZHANG Lixiang，etal.Frequency domain analysis of fluid-structure interaction in liquid-filled pipe systems by transfer matrix method［J］.International Journal of Mechanical Sciences，2002，44：2067-2087.

［5］張立翔，楊柯.流體結(jié)構(gòu)與互動(dòng)理論及其應(yīng)用［M］.北京：科學(xué)出版社，2004：37-42.

［6］朱洪來(lái)，白象忠.流固耦合問(wèn)題的描述方法及分類簡(jiǎn)化準(zhǔn)則［J］.工程力學(xué)，2007，24（10）：92-99.ZHU Honglai，BAI Xiangzhong.Description method and simplified classification rule for fluid-solid interaction problems［J］.Engineering Mechanics，2007，24（10）：92-99.

［7］偶國(guó)富，許根富，朱祖超，等.彎管沖蝕失效流固耦合機(jī)理及數(shù)值模擬［J］.機(jī)械工程學(xué)報(bào)，2009，45（11）：119-124.OU Guofu，XU Genfu，ZHU Zuchao，etal.Fluidstructure interaction mechanism and numerical simulation of elbow erosion failure［J］.Chinese Journal of Mechanical Engineering，2009，45（11）：119-124.

［8］王宇飛，楊振國(guó).氣-液兩相流對(duì)管道性能影響的有限元模擬分析［C］／／第14屆全國(guó)反應(yīng)堆結(jié)構(gòu)力學(xué)會(huì)議論文集.桂林：中國(guó)力學(xué)學(xué)會(huì)，中國(guó)核學(xué)會(huì).2006：99-101.

［9］劉志遠(yuǎn)，鄭源，張文佳，等.Ansys-CFX 單向流固耦合分析的方法［J］.水利水電工程設(shè)計(jì)，2009，28（2）：29-31.LIU Zhiyuan，ZHENG Yuan，ZHANG Wenjia，et al.The method of single-track fluid-structure interaction analysis by Ansys-CFX ［J］.Design of Water Resources ＆Hydroelectric Engineering，2009，28（2）：29-31.

［10］楊濤.液力緩速器流場(chǎng)仿真及有限元分析［D］.武漢：武漢理工大學(xué)汽車工程學(xué)院，2009：32-37.

［11］曾強(qiáng).壓氣機(jī)轉(zhuǎn)子葉片流固耦合計(jì)算及軟件集成研究［D］.南京：南京航空航天大學(xué)能源與動(dòng)力學(xué)院，2006：20-25.

［12］王鄢.大型離心壓縮機(jī)葉輪動(dòng)力分析方法研究［D］.大連：大連理工大學(xué)工程力學(xué)系，2009：30-33.

［13］馬廷衛(wèi).混流式水輪機(jī)轉(zhuǎn)輪結(jié)構(gòu)分析研究［D］.四川：西華大學(xué)能源與環(huán)境學(xué)院，2006：20-22.

［14］韓占忠，王敬，蘭小平.Fluent流體工程仿真計(jì)算［M］.北京：北京理工大學(xué)出版社，2008.

［15］KIM Y，KIM J，JEON Y，etal.Multidisciplinary aerodynamic-structural design optimization of supersonic fighter wing using response surface methodolo-gy［C］／／40th AIAA Aerospace Sciences Meeting.Reno：NV，AIAA，2002.

［16］安效民，徐敏，陳士櫓.二階時(shí)間精度的CFD／CSD 耦合算法研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2009，27（5）：547-552.AN Xiaomin，XU Min，CHEN Shilu.Analysis for second order time accurate CFD／CSD coupled algorithms［J］.Acta Aerodynamica Sinica，2009，27（5）：547-552.

［17］WEN Jihong，SHEN Huijie，YU Dianlong，etal.Theoretical and experimental investigation of flexural wave propagating in a periodic pipe with fluid-filled loading［J］.Chin Phys Lett，2010，27（11）：133-136.