孫金雄 周瑞

摘 要：管道流致振動現(xiàn)象廣泛存在于核電廠內(nèi)，使管道產(chǎn)生疲勞和噪聲，對管系安全造成威脅。文章對管道流致振動的原理進行探討，并通過對兩端簡支輸流管振動總響應(yīng)計算方法的簡化，得到輸流管線撓度與流速的關(guān)系式，并與流固耦合數(shù)值仿真計算進行對比，兩者結(jié)果非常吻合，這表明用線性流固耦合仿真技術(shù)為核電廠管道設(shè)計與驗證提供參考是可行的。

關(guān)鍵詞：流固耦合；流致振動；管道內(nèi)流

1 概述

核電廠內(nèi)分布大量管道系統(tǒng)，運行時發(fā)生流致振動現(xiàn)象十分普遍，振動會使管道產(chǎn)生疲勞和噪聲，對管系安全造成威脅。常用工業(yè)規(guī)范（ASME、RCCM）對于管道流致振動的控制與驗證未給出明確方法，因此需從原理上探索流致振動的原因，并找到可應(yīng)用于工程的實用方法。雖然工業(yè)規(guī)范不完善，但流致振動問題卻是國內(nèi)外學(xué)者的研究熱點，產(chǎn)生了很多重要研究成果，比如Blevins早在1990年就發(fā)表了經(jīng)典專著[3]，對流致振動原因進行總結(jié)和分析；Paidoussis于2014年發(fā)表的管道內(nèi)流流致振動專著[1]，對管道線性與非線性流固耦合方程進行全面介紹。以上研究都基于理論公式推導(dǎo)，并未與當(dāng)前商用數(shù)值仿真軟件融合。文章通過對輸流管道線性運動方程的推導(dǎo)和簡化，得到管內(nèi)流速與管道最大撓度的關(guān)系，并與用仿真軟件計算的結(jié)果對比，兩者結(jié)果非常吻合，表明數(shù)值仿真方法在輸流管道線性流固耦合分析中精度較高，可用于工程設(shè)計。

2 簡支輸流管線性流固耦合運動方程

2.1 簡支輸流管線計算模型

如圖1所示，兩端簡支的輸流管線偏離平衡位置，形成橫向撓度Y，跨距為L的跨距。我們從管道上切下管道和流體兩個微元，為每個微元進行受力分析。流體密度為？籽；壓力為p；流速為v。管子橫截面積A；長度為L；彈性模量為E；截面慣性矩為I；管道單位長度質(zhì)量為M。

2.2 管道與流體微元平衡方程

由于管道振動，流體產(chǎn)生加速度、同時豎直方向有流體壓力分量、管壁作用在流體上的壓力F，上述力在豎直方向平衡可得：

流體沿管子長度方向壓力梯度由管壁摩擦的剪切應(yīng)力平衡，其方程為：

其中S為管內(nèi)壁周界，q為管內(nèi)表面剪切應(yīng)力。

定義M為無水管道單位長度質(zhì)量，管子橫向剪切應(yīng)力為Q，管子縱向拉應(yīng)力為T，管線豎直方向加速度為：

定義T為管子縱向拉應(yīng)力，Q為管子橫向剪切應(yīng)力，沿管軸線方向力平衡方程為：

聯(lián)立（1）、（2）、（3）消去變量F和Q，剪切應(yīng)力q從（2）（4）中消去，最終得到輸流管道的自由振動的運動方程：

3 方程的解及管道變形的理論近似解

3.1 方程的邊界條件和方程的解

對于兩端簡支輸流管的邊界條件為如下：

當(dāng)x=0和x=L時，管道撓度與加速度為0。

由參考文獻[4]知該方程的解的表達式如下：

其中j=1，2，3，…，wj是j振型的固有頻率，表示為具有正弦和余弦的對稱和反對稱的各個空間振型的總和。

因為系統(tǒng)具有無限個自然振型，而實用解只需考慮頭幾階振型。根據(jù)參考文獻[3]，可以得到管內(nèi)流體流動時管線一階固有頻率w1和管內(nèi)流體靜止時管道固有頻率wN比值近似為：

其中vc為管道靜態(tài)失穩(wěn)時的臨界流速：

3.2 管道最大變形與流速的關(guān)系推導(dǎo)

對于前面計算的兩端簡支管道，單位長度載荷為q。則該管道中心處撓度？啄n計算公式為：

管道固有頻率wn計算公式為：

通常只考慮前兩階振型的響應(yīng)，由參考文獻[3]知第1階振型a1 與第2階振型a2比值為：

因此當(dāng)管道以基波固有頻率振動時，正弦波彎曲的管道第1階振型在管道的響應(yīng)中起主導(dǎo)作用。如果將第1階的響應(yīng)值近似等于管道總響應(yīng)值，則可以得到最大撓度？啄n與第一階固有頻率wn1的關(guān)系如下：

由上式可知，撓度僅與固有頻率的平方成反比關(guān)系。當(dāng)管內(nèi)水流動時，其撓度為？啄i，一階固有頻率為wi1，則：

由公式（10），得出輸流管線撓度？啄i與？啄n的關(guān)系為：

（14）

4 軟件流固耦合計算

4.1 雙向迭代耦合與直接耦合區(qū)別

雙向迭代耦合是指流體和結(jié)構(gòu)用各自求解器在時域積分，通過耦合面在流體和結(jié)構(gòu)之間進行位移與壓力的傳遞。

直接耦合求解是通過數(shù)學(xué)一致化的處理方法，將結(jié)構(gòu)方程組、流體方程組、FSI 界面方程組等不同性質(zhì)的方程組、借助于位移、壓力在不同方程組中的等價性，形成一個方程組耦合系統(tǒng)方程組。

4.2 計算模型與相關(guān)參數(shù)

4.2.1 幾何模型

計算幾何模型與圖1相同，結(jié)構(gòu)為靜態(tài)分析，流體為穩(wěn)態(tài)分析。由于ADINA支持直接耦合求解，所以該算例在ADINA-FSI中進行。

4.2.2 材料參數(shù)

計算采用國際標(biāo)準(zhǔn)單位，材料參數(shù)見表1。

4.2.3 邊界條件

流體為實體單元，管道為殼單元；管道與流體耦合面為有滑移耦合面，選擇滑移面是為匹配理論近似解的假設(shè)，同時滑移邊界對網(wǎng)格粗糙不敏感，無需邊界層。重力加速度為9.81m/s2，在入口處施加線性遞增流速，最大值為臨界失穩(wěn)速度；由公式（8）知vc=149.6931m/s。本次為對比迭代與直接耦合計算差異，進行了兩次計算，耦合方式分別為迭代耦合和直接耦合，其他邊界條件均未發(fā)生變化。

4.3 計算模型與相關(guān)參數(shù)

4.3.1 ADINA計算結(jié)果

（1）迭代耦合與直接耦合。迭代耦合計算出管道最大豎直位移值為-4.668mm；直接耦合計算出最大值為-4.713mm，與迭代耦合幾乎一致。

（2）理論近似解由公式（9）和算例參數(shù)，可以得到該輸流管道靜止時，其中心處的撓度？啄n=2.1112×10-4m。根據(jù)公式（14）可知管內(nèi)流體流動后，管道中心處的撓度表達式為：

其中v為入口處速度。將上式作為ADINA后處理自定義函數(shù)表達式，把管道入口流速作為該表達式自變量，可得管道最大撓度隨入口速度的關(guān)系曲線。

4.3.2 結(jié)果對比

將上述3種分析方法得到的速度-處撓度曲線繪到同一圖內(nèi)進行對比，如圖2所示。

4.3.3 結(jié)果差異分析

（1）計算終止時，入口速度未達到臨界速度149.6931m/s，因為接近臨界速度時，管道剛度為0，產(chǎn)生失穩(wěn)，計算由于不收斂終止。（2）近似精確解與ADINA計算結(jié)果相比，在接近臨界速度時，有一定的差別，但差別較小。造成差別的主要原因是在推導(dǎo)公式時，為簡化計算，只考慮管道的第1階振型對流體激勵的響應(yīng)，后續(xù)的n階響應(yīng)都被忽略，造成估算算出的位移小于實際的響應(yīng)。（3）線性流固耦合問題分析時，迭代耦合與直接耦合結(jié)果差別較小，直接耦合精度更高，但實際計算時，迭代耦合計算效率比直接耦合高很多。

5 結(jié)束語

綜合以上結(jié)論可認為，基于數(shù)值模擬的線性流固耦合仿真技術(shù)，計算精度比簡化公式更高，且當(dāng)結(jié)構(gòu)復(fù)雜，公式無法適用時，該方法優(yōu)勢更加明顯，因此可將其推廣到核電廠的管道設(shè)計應(yīng)用中。

參考文獻

[1]Michael P.Paidoussis.Fluid-structure Interactions Slender structures and Axial Flow：Second Edition VOLUME 1[M].UK：The Academic Press of Elsevier，2014：63-332.

[2]Shigehiko Kaneko，Tomomichi Nakamura. Flow-induced Vibrations Classifications and lessons from Practical Experiences： Second Edition[M]. USA： The Academic Press of Elsevier，2014：157-269.

[3]Robert D.Blevins.Flow-Induced Vibration：Second Edition[M]. Malabar： Krieger Publishing Company，1990.

[4]Housner，G.W.Bending Vibration of a Pipe Line Containing Flowing Fluid：[M].J.Appl.Mech，1954：205-208.

[5]張阿漫.流固耦合動力學(xué)[M].北京：國防工業(yè)出版社，2010.