輸流管道線性流固耦合計算方法在核電工程中的應用

孫金雄 周瑞

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

關鍵詞：流固耦合；流致振動；管道內流

1 概述

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

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

2.1 簡支輸流管線計算模型

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

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

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

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

其中S為管內壁周界，q為管內表面剪切應力。

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

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

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

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

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

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

當x=0和x=L時，管道撓度與加速度為0。

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

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

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

其中vc為管道靜態失穩時的臨界流速：

3.2 管道最大變形與流速的關系推導

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

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

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

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

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

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

（14）

4 軟件流固耦合計算

4.1 雙向迭代耦合與直接耦合區別

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

直接耦合求解是通過數學一致化的處理方法，將結構方程組、流體方程組、FSI 界面方程組等不同性質的方程組、借助于位移、壓力在不同方程組中的等價性，形成一個方程組耦合系統方程組。

4.2 計算模型與相關參數

4.2.1 幾何模型

計算幾何模型與圖1相同，結構為靜態分析，流體為穩態分析。由于ADINA支持直接耦合求解，所以該算例在ADINA-FSI中進行。

4.2.2 材料參數

計算采用國際標準單位，材料參數見表1。

4.2.3 邊界條件

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

4.3 計算模型與相關參數

4.3.1 ADINA計算結果

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

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

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

4.3.2 結果對比

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

4.3.3 結果差異分析

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

5 結束語

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

參考文獻

[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]張阿漫.流固耦合動力學[M].北京：國防工業出版社，2010.