基于雙向流固耦合的輸流圓管應力應變響應分析

劉 松，吳先梅，彭修乾，陳家熠，王 聰

(1.中國科學院 聲學研究所 聲場聲信息國家重點實驗室，北京 100190;2.中國科學院大學，北京 100049;3.國家超級計算天津中心，天津 300457)

輸流管道廣泛應用于能源環(huán)保、石油化工、儲運交通等領域，在工業(yè)生產(chǎn)及日常生活中發(fā)揮著非常重要的作用。當流體在輸流管道內流動時，其流速與壓強的變化會誘發(fā)管道結構產(chǎn)生形變以及振動，而同時管道結構的形變和振動也會反之影響管道內流體的流動，二者相互作用、相互影響，這種流固交互作用一般稱為流固耦合(fluid-structure interaction,FSI)作用[1]。而這種由于流固間的相互作用所誘發(fā)的輸流管道振動會直接影響到管內流體流動狀況的檢測。因此，對輸流管道的流固耦合作用進行研究，分析其流固耦合特性，這對測量輸流管道內的流體參數(shù)具有重要的指導意義。

早在1933年，Westergaard結合美國地震區(qū)建大壩的需要，發(fā)表了題為“Water Pressure on dams during earthquakes”的文章，首次開展了流固耦合振動的研究，但直到20世紀90年代才引起學者們的廣泛關注。近些年隨著計算機輔助工程(computer aided engineering,CAE)技術的發(fā)展，有很多軟件集成了適用于對結構模型進行流固耦合分析的模塊，這極大地提高了工作效率，促使國內外許多學者在流固耦合問題上進行了積極的探索與研究。柳博瀚等[2]研究了不同管內流速度對彈性管道振動響應的影響，他們通過數(shù)值仿真與模型試驗的相比較，驗證了所采用的雙向流固耦合計算方法的有效性，并得出結論：管內流體是通過增加軸向動力以及改變彈性管道局部曲率的方式影響管道結構的振動頻率；在內流速度相對較小時，管道內流對管道的振動頻率及振幅影響都較小；在內流速度相對較大時，管道內流對管道橫流向振幅及順流向振型影響較為明顯，使管道在低頻段更容易發(fā)生共振。張林等[3]應用ANSYS Workbench軟件構建了管內流體與彎管結構的有限元模型并進行了流體激勵下輸流彎管的流固耦合振動特性研究，在分析管內流體的波動頻率對輸流彎管流固耦合振動幅值的影響后得出結論：橫向振動為輸流彎管流固耦合振動的主要形式；流體的波動頻率與管道的固有頻率越接近，彎管結構的振幅越大。宋學官等[4]書中有一例是對一段兩端為固定約束的柔性管道進行流固耦合研究，其研究結果表明：輸流軟管的最大應力和最大應變發(fā)生在軟管的入口和出口附近，而最大側向位移發(fā)生在軟管的中部附近。袁松松[5]應用ANSYS Workbench軟件對一輸流彎管進行了仿真模擬，研究不同約束方式下的管道振動，其研究結果表明：輸流管道在輸送流體的過程中，流體速度的變化會誘導管道結構產(chǎn)生較大振蕩；不同的約束方式會改變管道結構的最大應力和最大位移，同時管道特征點的位移以及管道結構的最大位移和最大應力都與管內流體加速度的絕對值有關。趙江等[6]利用ANSYS Workbench軟件對流體作用下的T型管進行了雙向流固耦合模擬，通過分析管道在流體作用下的應力和變形情況，對T型管振動劇烈的結合部位進行諧響應分析，其研究表明：基于雙向流固耦合的模態(tài)分析可以較確切地反映內部流體作用下的管道模態(tài)；在各因素中，流體壓強對管道固有頻率的影響最大。Naveen等[7]分析了管道壁厚和彎管長度對流體流速、Mises等效應力以及湍流動能的影響并得出結論：隨著管道壁厚和彎管長度的增加，流體流速、Mises等效應力以及湍流動能都會減小；但隨著湍流強度的增大，流體流速和湍流動能都會隨之增大。Chen[8]對圓形充液彎管進行了流固耦合研究，分析了輸流彎管在不同邊界條件下管內流體流速與彎管穩(wěn)定性的關系，發(fā)現(xiàn)當管內的流體流速超過一定值后，彎管會發(fā)生皺曲失穩(wěn)。Chen等[9]應用ANSYS軟件的多物理場求解器(multi-field solver,MFS)模塊進行了輸流管道流固耦合數(shù)值模擬，仿真模擬PVC管內的流體流動。

上述有關輸流管道流固耦合特性的研究對于測量管內流體參數(shù)有一定的指導意義。而實際上，當管內流體的脈動壓力作用在管道壁面上時，脈動的流體壓力會以脈動的應力應變形式加載到管道內壁面上[10]。因此，為更進一步研究輸流管道的流固耦合作用，可通過對比分析內、外管壁處的應力應變信號與管內流體的脈動信號來探究它們的關系，為測量輸流管道內的流體參數(shù)提供指導。為此，我們分別以內徑d=60 mm、內徑d=40 mm的輸流圓管為研究模型進行雙向流固耦合仿真模擬，主要分析研究了圓管內、外壁處的應力應變信號與管內湍流脈動信號的關聯(lián)。

1 流固耦合基本方程

流固耦合計算可分為單向流固耦合與雙向流固耦合兩種類型。單向流固耦合是將流體運動方程計算得出的流體載荷數(shù)據(jù)通過流固耦合矩陣傳遞給固體運動方程進行計算，但不會將固體結構的計算結果反饋給管內流體運動方程，其計算數(shù)據(jù)的傳遞是單向的。它的應用前提是認為固體結構的形變微小，對流體運動的影響可以忽略不計。而對于雙向流固耦合，其計算數(shù)據(jù)的傳遞是雙向的，即在每個計算時間步內，管內流體計算所得的數(shù)據(jù)會傳遞給固體運動方程進行計算，而固體運動方程的計算結果數(shù)據(jù)又會傳遞回管內流體運動方程進行計算，因為它考慮了管內流體與固體結構兩者之間的相互影響，所以最終計算所得的結果會比單向流固耦合所得的結果更符合實際情況。因而，本次的仿真模擬采用了雙向流固耦合計算。

1.1 流固耦合有限元方程

在研究輸流圓管的流固耦合作用時，首先考慮管內流體的運動，其有限元方程為[11]

(1)

然后為確定圓管結構對管內流體運動的影響，僅考慮在內部流體作用力下的位移和振動，不考慮包括圓管結構外載荷的其它外力，則管道結構有限元方程為

(2)

由上述論述可知，流固耦合作用體現(xiàn)在式(1)中圓管結構對管內流體的載荷項{Ff}與式(2)中管內流體對圓管結構的載荷項{Fs}[12]，其中

(3)

{Fs}=[Rfs]T{p}

(4)

式中：ρf為流體密度；[Rfs]為流固交界面上的耦合矩陣。

在忽略阻尼影響的條件下，聯(lián)立式(1)～式(4)，可得

(5)

式(5)即為圓管結構與管內流體的直接耦合方程。

1.2 流固耦合界面方程

圖1為圓管結構與管內流體在流固耦合交界面處力與位移的傳遞關系示意圖。由于管內流體與管道結構的耦合作用實際上只發(fā)生在流體與管道結構的交界面處，是在交界面處進行物理量的傳遞。因此，在流固耦合交界面處應滿足法向速度(v)連續(xù)、法向應力(τ)連續(xù)以及位移(x)連續(xù)，可分別表示為

圖1 耦合交界面處力與位移傳遞關系示意圖Fig.1 Force and displacement transmit at the interface

vfn=vfnf=vsnf=-vsns=vsn

(6)

τfn=τsn

(7)

xf=xs

(8)

式中：vfn，vsn分別為流固耦合交界面處流體域側的法向速度、固體域側的法向速度；τfn，τsn分別為流固耦合交界面處流體域側的法向應力、固體域側的法向應力；xf，xs分別為流固耦合交界面處流體域側的位移、固體域側的位移。

2 有限元仿真模擬

實際輸流管道內的流體90%以上是以湍流狀態(tài)流動，因此本次模擬的管內流體運動狀態(tài)為湍流。對管內湍流的數(shù)值模擬可以選用直接數(shù)值模擬(direct numerical simulation,DNS)、雷諾平均模擬(reynolds-averaged navier-stokes,RANS)和大渦數(shù)值模擬(large eddy simulation,LES)三種。直接數(shù)值模擬雖然可以模擬所有尺度的湍流運動，具有準確性高的優(yōu)點，但是因其需要非常巨大的計算量導致了難以將之應用于工程實際之中。雷諾平均模擬的計算量較小，因此會常將其應用于工程問題之中[13]，但它的計算準確性較差。而大渦模擬是介于雷諾平均模擬與直接數(shù)值模擬之間的一種數(shù)值模擬方法，它具有兩者的優(yōu)點，非常適合用于本次管內湍流的數(shù)值模擬。

應用ANSYS軟件進行輸流管道的流固耦合模擬計算。取內徑d=60 mm、管長l=300 mm、壁厚δ=4 mm的輸流圓管作為研究模型。首先建立管內流體的有限元模型，因為該流體模型為圓管且形狀規(guī)則，所以采用質量較高的六面體八節(jié)點單元對整個流體模型進行結構化網(wǎng)格劃分，劃分后的流體有限元模型如圖2(a)所示。隨后在Fluent模塊中設置相關參數(shù)：設為瞬態(tài)(Transient)；流體設置為水，其密度為998.2 kg/m3，動力黏度為0.001 Pa·s；考慮管內流體所受重力的影響，在整個流體域Y方向上施加重力載荷，重力加速度為-9.81 m/s2；湍流模型選擇為大渦數(shù)值模擬，并選用WALE亞格子模型[14]，是因為WALE模型的渦黏系數(shù)會在管道壁面附近自然趨近于零，較符合管內流體的實際流動情況，并且由以往的湍流模擬可知它能在很寬的流量范圍內產(chǎn)生令人滿意的流體仿真結果；將流體入口(inlet)設置為速度入口(velocity inlet)，流速為1.6 m/s；將流體出口(outlet)設置為壓力出口(pressure outlet)，靜壓為0 Pa；求解計算的數(shù)值格式設置為有限體積法；壓力-速度耦合方程求解算法設置為Simple算法；對流插值設置為有界中心差分格式。

圖2 有限元模型Fig.2 Finite element model

然后建立圓管結構的有限元模型，同樣選用六面體八節(jié)點單元對圓管結構模型進行結構化網(wǎng)格劃分，劃分后的有限元模型如圖2(b)所示。在Transient Structural模塊中設置該管件材料為聚氯乙烯(PVC)，其楊氏彈性模量為3.4 GPa，泊松比為0.319，密度為1 380 kg/m3；考慮管道結構所受重力的影響，在管道結構Y方向上同樣施加重力載荷，重力加速度為-9.81 m/s2；將管道內壁內表面與管內流體外接觸面均設為流固交界面(fluid solid interface)，并且考慮到管道內壁與管內流體的接觸界面網(wǎng)格，為統(tǒng)一接觸界面的網(wǎng)格尺寸，在進行網(wǎng)格劃分時，已將圓管結構內壁內表面和管內流體外接觸面的網(wǎng)格尺寸均設為2 mm；考慮管道結構的約束狀態(tài)，將圓管結構模型的入口、出口兩端部設為固定約束(fixed support)，而圓管結構的中部無約束。

最后用ANSYS軟件中的系統(tǒng)耦合(system coupling)模塊來實現(xiàn)流體載荷數(shù)據(jù)與固體載荷數(shù)據(jù)的交互。設置仿真總時長為0.200 s，時間步長為0.001 s。經(jīng)迭代運算后，最終便可得到輸流圓管雙向流固耦合的仿真模擬結果。

3 流體脈動與管壁應力應變分析

圖3為輸流圓管雙向流固耦合的仿真模擬圖。由圖3(a)管內流體壓強分布云圖、圖3(b)圓管結構彈性等效應變分布云圖以及圖3(c)圓管結構Mises等效應力分布云圖可以較直觀地看出管內流體壓強能通過流固耦合交界面的傳遞以應力應變形式加載到圓管結構上。但是，若仔細觀察這三個云圖，便發(fā)現(xiàn)管內流體壓強分布與等效應力應變分布存在一定差異。

圖3 雙向流固耦合仿真模擬圖Fig.3 Simulation diagram of the two-way fluid-structure interaction

為了更準確的研究輸流圓管內、外壁處應力應變信號與管內湍流脈動信號的關聯(lián)，我們從計算結果中提取了距輸流圓管模型入口100 mm處橫截面上的數(shù)據(jù)進行對比分析。在圓管結構橫截面處內、外壁面上8等分圓周(360°)設數(shù)據(jù)提取點-載荷點(Point Load，ANSYS中的硬點)，如圖4所示。內壁載荷點在流固耦合交界面(Interface)上，分別標記為I1，I2，…，I8；外壁載荷點在固體(Solid)結構上，分別標記為S1，S2，…，S8。這內、外壁面等間隔45°的8個載荷點在徑向上分別一一對應，即I1對應S1，I2對應S2，…，I8對應S8。

圖4 圓管模型橫截面示意圖Fig.4 The cross section of pipe model

取圓管結構內壁載荷點I1處的等效應變信號、等效應力信號與管內流體的脈動壓強信號作對比，分別如圖5(a)、圖5(b)所示，圖中左側y坐標軸對應藍色的脈動壓強曲線，右側y坐標軸對應綠色的等效應變或等效應力曲線。由圖可以明顯地看出載荷點I1處的等效應力應變信號變化與管內的湍流脈動信號變化不一致。可以用相關系數(shù)來定量描述兩個信號之間的相關程度，相關系數(shù)計算公式為

圖5 載荷點I1處脈動壓強信號與等效應變信號和等效應力信號對比曲線圖Fig.5 The contrast curves of fluctuating pressure and the equivalent strain,the equivalent stress

(9)

式中：R(i,j)為信號向量i與信號向量j的相關系數(shù)；C(i,j)為信號向量i與信號向量j的協(xié)方差矩陣。若信號向量i與信號向量j之間的相似程度越高，則相關系數(shù)值越接近于1。

根據(jù)式(9)求得載荷點I1處等效應變信號、等效應力信號與管內流體脈動壓強信號的相關系數(shù)分別為0.138 7，0.138 6，都與1差距很大，由此可知兩信號之間存在很大差異。為避免偶然性，分別對比分析了圖3中其它內壁載荷點I2，I3，…，I8處的等效應變、等效應力信號與管內流體脈動壓強信號，均發(fā)現(xiàn)存在很大差異，如表1所示。由此可得，難以用等效應力應變信號表征管內的湍流脈動信號。

表1 相關系數(shù)Tab.1 Correlation coefficient

考慮到輸流圓管為圓周結構，根據(jù)式(10)、式(11)將模擬計算所得的應力應變分量(直角坐標系下)轉換為徑向應變εr與周向應變εθ，表達式為

εr=εxcos2θ+εysin2θ+γxysinθcosθ

(10)

εθ=εxsin2θ+εycos2θ-γxysinθcosθ

(11)

式中：εr，εθ分別為徑向主應變、周向主應變；εx，εy分別為直角坐標系下x方向主應變、y方向主應變；γxy為xy切應變；θ為該點所在圓周半徑與x軸正向的夾角。

取載荷點I1處的徑向應變εr、周向應變εθ分別與管內該點處的脈動壓強信號作對比，如圖6(a)、圖6(b)所示。分析后可得圓管結構內壁載荷點I1處的周向應變信號與管內該點處的湍流脈動信號存在很強的正相關關系，其相關系數(shù)為0.998 9；而徑向應變信號與湍流脈動信號存在很強的負相關關系，其相關系數(shù)為-0.997 1。

圖6 載荷點I1處脈動壓強信號與徑向應變信號和周向應變信號對比曲線圖Fig.6 The contrast curves of fluctuating pressure and the radial strain,the circumferential strain

繼續(xù)對比分析圓管結構內壁載荷點I2處的徑向應變信號、周向應變信號與管內該點處的脈動壓強信號，分別如圖7(a)、圖7(b)所示。分析后可得圓管結構內壁載荷點I2處的周向應變信號與管內該點處的湍流脈動信號存在很強的正相關關系，其相關系數(shù)為0.998 0。再取內壁載荷點I3、…、I8處的周向應變信號分別與管內脈動壓強信號作對比，均可得出同樣的結論：管內壁處的周向應變信號與管內壁附近的湍流脈動信號存在很強的正相關關系，如表2所示。這表明可以用圓管內壁的周向應變信號表征管內壁附近的湍流脈動信號。

圖7 載荷點I2處脈動壓強信號與徑向應變信號和周向應變信號對比曲線圖Fig.7 The contrast curves of fluctuating pressure and the radial strain,the circumferential strain

表2 相關系數(shù)Tab.2 Correlation coefficient

接下來，取圓管結構內壁載荷點I1處的徑向應力信號、周向應力信號與管內該點處的脈動壓強信號作對比，分別如圖8(a)、圖8(b)所示。分析后可得圓管結構內壁載荷點I1處的周向應力信號與管內該點處的湍流脈動信號存在很強的正相關關系，其相關系數(shù)為0.999 4；徑向應力信號與湍流脈動信號存在很強的負相關關系，其相關系數(shù)為-0.999 9。對比分析圓管內壁其它載荷點，均可得出同樣的結論：管內壁處的周向應力信號與管內壁附近的湍流脈動信號存在很強的正相關關系，如表2所示。綜上論述可知，能用圓管內壁的周向應力應變信號表征管內壁附近的湍流脈動信號。

圖8 載荷點I1處脈動壓強信號與徑向應力信號和周向應力信號對比曲線圖Fig.8 The contrast curves of fluctuating pressure and the radial stress,the circumferential stress

對圓管結構內壁載荷點與相對應外壁載荷點處的周向應力應變信號進行分析。取內壁載荷點I1與外壁載荷點S1處的周向應變信號、周向應力信號作對比，分別如圖9(a)、圖9(b)所示。分析后可得輸流圓管內壁面與外壁面處的周向應變信號、周向應力信號存在很強的正相關關系，其相關系數(shù)分別為0.999 5、0.999 3，僅幅值大小有差距，但相差很小(約1.4倍)。對比分析其余各對應載荷點，均可得出同樣的結論：管內壁面與外壁面處的周向應變信號、周向應力信號存在很強的正相關關系，如表3所示。綜上可得，可以用輸流圓管外壁處的周向應力應變信號間接表征管內壁附近的湍流脈動信號。

表3 相關系數(shù)Tab.3 Correlation coefficient

圖9 載荷點I1與載荷點S1周向應變應力信號對比曲線圖Fig.9 The circumferential stress-strain signal comparison curve of point I1 and point S1

為避免偶然性，在其它設置都相同的條件下，構建內徑d=40 mm、管長l=300 mm、壁厚δ=4 mm的輸流圓管為研究模型并進行雙向流固耦合仿真模擬，繼續(xù)分析其管壁處應力應變信號與管內湍流脈動信號的關聯(lián)。

取圖10所示圓管結構外壁載荷點P1處的周向應變信號、周向應力信號與管內脈動壓強信號作對比，分別如圖11(a)、圖11(b)所示。分析后可得外壁載荷點P1處周向應變信號、周向應力信號與管內湍流脈動信號存在很強的正相關關系，其相關系數(shù)分別為0.999 6、0.999 5。θ為該點所在圓周半徑與x軸正向的夾角。沿管道外壁抽取不同角度θ處的周向應力應變信號與管內脈動壓強信號進行對比分析，均可得出同樣的結論：外管壁處的周向應力應變信號與管內湍流脈動信號存在很強的正相關關系，如表4所示。

表4 相關系數(shù)Tab.4 Correlation coefficient

圖10 圓管模型(40 mm)橫截面示意圖Fig.10 The cross section of pipe model(40 mm)

圖11 載荷點I1處脈動壓強信號與周向應變信號和周向應力信號對比曲線圖Fig.11 The contrast curves of fluctuating pressure and the circumferential strain,the circumferential stress

綜上所述可知，可以通過檢測輸流圓管外壁處的周向應力應變信號間接測得管內壁附近的湍流脈動信號。

4 結 論

本文研究了輸流圓管管壁處的應力應變信號與管內湍流脈動信號的關系，應用ANSYS軟件進行了輸流圓管雙向流固耦合的仿真模擬。通過分析仿真模擬結果發(fā)現(xiàn)，管內流體壓強能通過流固耦合交界面的傳遞以應力應變形式加載到圓管結構上。管道內壁附近的湍流信號與管道內、外壁面的周向應力應變信號存在很強的正相關關系，而與管道內、外壁面的徑向應力應變信號存在很強的負相關關系。可以通過檢測輸流圓管外壁處的周向應力應變信號間接測量出管內壁附近的湍流脈動信號，這為有效地測量管內流體脈動信號提供了很好的方法。