細長輸流管內外流耦合振動特性研究

鮑 健，陳正壽

(1.浙江海洋大學 船舶與海運學院，浙江 舟山 316022；2.中國海洋大學 工程學院，山東 青島 266100)

隨著海洋資源的開采逐漸走向更深的海域，大長徑比的深海立管/管線應用越來越廣泛。深海輸流管線主要運輸高溫高壓油氣，沿水深的不均勻流動容易激發細長彈性管的多模態振動響應，另外管內流體高速輸運也會改變彈性管結構的動力特性[1-2]，內外流耦合作用對管線整體運動響應及疲勞損傷影響顯著。內外流與細長柔性管線之間發生流固耦合的機理是不同的，外流以一定的速度流過管線會誘發尾渦脫落，進而引起彈性管體渦激振動(VIV)。當尾渦脫落頻率與管體結構固有頻率接近甚至相等時，“鎖定”現象發生，從而增強管線振動并加速管線的疲勞破壞[3]。當內流沿管線軸向運動時，附加質量力、離心力和科氏力誘發管體流致振動(FIV)。這種雙流固耦合工況(外流—管體，內流—管體)會導致管線動力學特性相較單流固耦合工況更加復雜。

在柔性立管/管線渦激振動方面，國內外學者已經做了相當多探索性的研究工作。Chaplin等[4]通過階梯流下柔性圓柱體振動響應試驗，獲得了沿軸向和橫流向(CF)的一系列激發波數。Trim等[5]對均勻剪切流動下的VIV響應進行了試驗研究，發現順流向的疲勞損傷程度與橫流向相同。Wang和Xiao[6]基于大渦模擬湍流模型對豎直立管在均勻和線性剪切流下的渦激振動過程進行數值研究，得到的結果與試驗數據吻合良好。當軸向內流沿非線性幾何形狀的彎曲管體運動時，額外的離心力和科氏力使管線強非線性內外流耦合響應過程更加復雜，考慮內流的彈性管渦激振動問題在近十幾年來才得到了明顯的關注。Guo和Lou[7]通過彈性管內流的模型試驗，發現隨著內流流速的增加，彈性管的振動幅度將增大，而振動頻率將減小。Duan等[8]通過數值研究發現，行波在輸流管線的橫流向和順流向(IL)VIV響應中都很明顯。Yang等[9]通過參數化的數值研究發現，立管的最大位移和應力可以根據內部流速而增加或減小，并且臨界內流流速會導致不同外流速度下立管的主振模態階數增加。Chen和Kim[10]通過對彈性管在上行、下行內流和剪切外流下振動響應的數值仿真，發現內流的輸送方向在決定大跨度管體的動力響應方面起著不同的作用。柳博瀚等[11]基于計算流體力學(CFD)方法開展管體內流對彈性管振動響應影響的研究，結果表明較高的內流流速能夠顯著影響管體的振動特性。

目前關于內外流耦合下彈性管渦激振動的研究主要關注內流流速的變化對立管振動頻率、振幅，主振模態等特性的影響，如要全面評估彈性管在實際工作環境的動力響應和疲勞損傷，則需對內外流耦合誘發的多模態共振，行、駐波等特性開展更深入研究。文中首先介紹采用的數值方法，然后描述數值計算模型，并驗證數值方法的可靠性，最后基于各工況的數值結果，分析內外流耦合下細長彈性管的動力學特性，并闡述管內渦結構的演變規律。

1 數值方法

1.1 流場控制方程

采用多物理場耦合分析軟件STAR-CCM+[12]，研究內外流耦合下柔性管的強流固耦合特性。大渦模擬(LES)及其動態的Smagorinsky SGS模型[13]被用于求解流體域的水動力。過濾處理后的Navier-Stokes方程為：

(1)

(2)

采用有限體積法對流場進行離散化求解，壓力速度耦合采用SIMPLE算法，動量方程的離散采用有界中心差分格式。

1.2 結構動力學控制方程

將柔性管體簡化為兩端簡支的歐拉—伯努利梁模型[14]，可表示為：

(3)

式中：EI表示管體彎曲剛度；T為立管頂張力，m為管體系統的質量，c為結構阻尼，Fi為立管所受的水動力載荷分量。通過有限元法將式(3)離散化，再采用紐馬克—貝塔法求解。

(4)

通過流場域和固體域之間的數據映射與交換實現了雙向流固耦合，其中徑向基函數方法[15]被用于流場網格邊界和內部節點的移動插值。

2 數值仿真模型

2.1 計算域及網格劃分

選取Chen等[16]開展的剪切流下玻璃纖維增強塑料(FRP)管渦激振動試驗作為計算參考模型，并展開進一步的數值模擬研究。FRP彈性管模型主要參數如表1所示。

表1 FRP管模型的主要參數

圖1為該FRP管數值模型示意，沿管體A端(最小剪切流速，Vmin=0)至B端(最大剪切流速，Vmax)施加線性剪切外流，彈性管模型A端固定，B端鉸接，并且在B端施加一個恒定預張力，其中管內均勻流自A端運輸至B端。如圖2所示，外流場計算域網格劃分采用多面體網格和定向網格技術[12]，管體近壁面流體網格劃分較密，保證了y+<1，計算域遠場的網格較為稀疏，以適應網格的大變形。外流沿x軸正向流動，流動入口距立管橫截面中心15D，流動出口距立管橫截面中心50D，上下邊界分別距立管橫截面中心15D(D為FRP管外徑)。邊界條件設置為：外流入口為速度入口邊界，來流方向與邊界垂直，外流出口設置為參考壓力為0的壓力出口邊界，與立管外表面重合的流體域表面設置為無滑移壁面邊界，前后兩個側面設置為對稱邊界，上下兩個面設置為速度入口邊界，速度大小方向與入口邊界一致。在管體結構的有限元網格劃分方面，立管橫截面為環形結構，沿周長均勻劃分40個網格節點，沿徑向多層分布，并沿軸向以適當的層間距均勻層狀拉伸。

圖1 FRP管數值模型示意

圖2 外流場網格劃分

2.2 數值方法可靠性驗證

通過傅里葉功率譜識別主頻的方法，對采用不同度量單位的模型試驗和數值模擬數據進行對比驗證。圖3分別給出了2種剪切流速下的FRP管跨中處振動位移功率密度譜對比結果。圖3中，線1和線2分別代表模型試驗中測得的橫流向和順流向應變時間歷程功率譜密度，其對應的縱坐標為功率譜密度(模型試驗)，橫坐標為頻率單位Hz。線3和線4分別表示在數值仿真中監測到的橫流向和順流向振動位移時間歷程功率譜密度分布，其對應的縱坐標為功率譜密度(數值仿真)。橫坐標為頻率，單位Hz。

圖3 試驗值和數值模擬結果的頻譜對比

從以上2個頻譜對比分析可以發現，外流流速對確定主振模態的頻率起重要作用。當最大剪切流速度較低時(Vmax=0.3 m/s)，在數值模擬和試驗中觀察到的CF振動過程功率譜密度函數幾乎是單峰分布。當剪切流速增加到Vmax=0.9 m/s時，一些從振模態會被激發。可以得出結論，在數值模擬中觀察到的主頻分布與試驗結果相符。由此證明了文中采用的流固耦合數值方法的可靠性。

其次，使用大渦模擬(LES)及其動態的Smagorinsky SGS模型，進行了雷諾數Re=3 900的單圓柱繞流數值模擬。采用端木玉等[18]的計算模型，圓柱直徑D=0.01 m。計算域長度(x軸的方向)為15D，上游斷面距離圓柱中心為5D，下游斷面距離圓柱中心為15D，方向與來流方向一致；計算域寬度(z軸方向)為πD，方向沿著圓柱方向；計算域高度(y軸方向)為10D，上下面距離圓柱為5D，方向平行于圓柱橫截面。采用多面體非結構化網格對計算域進行劃分，為了滿足大渦模擬對近壁面網格密度的要求，保證了近壁面處y+<1。從表2可知，文中計算的結果與前人計算結果比較接近，驗證了數值模型的有效性。

表2 單圓柱模型驗證

2.3 網格依賴性

網格密度是影響數值結果準確性的關鍵因素，因此在做進一步的數值預報前首先進行了圓柱繞流問題的流體網格依賴性驗證。這里所用網格為多面體非結構化網格，并對柱體近壁區網格局部加密。對于來流速度U=1.5 m/s下剛性圓柱(D=34.8 mm，L=πD)繞流問題，采用了大渦模擬湍流模型，并考慮了沿柱體內壁周長方向節點的數量Nc，無量綱參數Δr/D(Δr是第一層網格高度)和無量綱參數Δz/D(Δz是軸向網格層間距)3種影響網格密度的關鍵參數，計算并比較了不同網格密度下(Mesh-1，Mesh-2，Mesh-3，Mesh-4)的平均阻力系數，均方根阻力系數和均方根升力系數，如表3所示。結果表明：Mesh-3和Mesh-4的均方根升力系數的誤差值減小到5.4%，相應的平均阻力系數和均方根阻力系數隨網格密度的增加幾乎沒有變化。可以得出結論：Mesh-3(其特征為Nc=160，Δr/D=0.000 5，Δz/D=0.005)適用于給定管體模型流體域的數值模擬。因此，在接下來的數值模擬工作中，采用Mesh-3中使用的外流場網格密度。使用STAR-CCM+[12]，在32個處理器(intel(R)Xeon(R)CPU E5-2620 v3 @ 2.40GHz)上，每一種工況耗時超過15×24 h。

表3 網格依賴性測試結果

3 計算結果及分析

3.1 振動強度的時間—空間分布

隨著管體振型及振動軌跡發生改變，沿管體軸線的振動能量也發生變化，通過對比在不同內流流速工況下管體振動位移在時間空間上分布云圖來分析內流對管體展向模態干涉和波形轉變的影響。外流Vmax=1.5 m/s不變，內流分別為0.5 m/s,1.0 m/s,1.5 m/s,2.0 m/s的4個工況下管體順流向與橫流向的位移—時程空間分布如圖4～7所示，橫坐標為標準化時間，縱坐標為標準化的沿管體軸向各監測點的位置。其中，深色區域表示該瞬時振動位置偏離初始位置較大。

圖4 內流Vin=0.5 m/s時，無量綱振動位移時程空間分布

隨著內流流速由0.5 m/s增長到2.0 m/s，管體橫流向振動表現為6階主振模態，順流向振動表現為9階主振模態。這表明了在一定的內流流速范圍內，外流流速對柔性管橫流向和順流向振型起決定性作用。另外，沿管跨方向振動位移的大小不但隨時間變化，而且與內流流速密切相關。當內流流速從0.5 m/s增加到1.0 m/s，振動強度提升，但是當內流流速由1.0 m/s增加到2.0 m/s，振動強度反而有所下降。這些都歸因于流體—固體和模態—模態間的連續能量傳遞。一個有趣的現象是，沿管跨方向局部的位移極值現象很明顯，表現為深色區域。如圖5(a)所示，在t/T=0.35～0.55的時間段中，存在一個位移極大值區域，在t/T=0.55～0.75的時間段中存在一個位移極小值區域。這意味著雖然這個工況下管體橫流向振動表現為6階振型，但是1階模態的振動強度很顯著。

圖5 內流Vin=1.0 m/s 時，無量綱振動位移時程空間分布

另外在順流向上，位移極值區域沿管長交替出現，如圖4(b)至圖7(b)。這表明了2階甚至3階模態開始參與振動能量的轉移，多模態共存及躍遷現象顯著。這種內流—管體和外流—管體的雙重流固耦合系統使管體的振動響應更加復雜。

圖7 內流Vin=2.0 m/s 時，無量綱振動位移時程空間分布

內流不僅可以改變管體振動強度，還能影響管體展向的駐波行波特性[20]。彈性FRP管的渦激振動過程表現為混合駐波行波主導，并且在選取時間范圍內，駐波行波交替出現。駐波模式比行波模式更占優勢，后者在管跨局部位置間歇性出現，尤其是在橫流向的振動中。這一發現表明由于管體長度有限且激發的振動模態不夠高，大部分振動能量仍集中在駐波區域。如圖4所示，在外流Vmax=1.5 m/s，內流0.5 m/s的工況下，管體橫流向和順流向都表現出駐波主導的混合振動模式，表現為在管體的兩端主要由駐波主導，且局部的行波組分在管中段間歇性出現。在t/T=0.65～0.85的時間段內，管體振動位移峰值位置不斷波動，并且振型從6階轉化到了5階。這種相鄰模態間轉變的出現體現出多模態間能量的傳遞。圖5為在外流Vmax=1.5 m/s，內流1.0 m/s的工況下管體順流向與橫流向的位移時程空間分布。相比于內流0.5 m/s的工況，橫流向振動表現出更多的行波組分和更加難以辨識的位移峰谷值位置，這說明此時內流流速的增加使得橫流向振動變得更加強烈。但是在順流向振動上卻發現駐波主導模式更加顯著。圖6為在外流Vmax=1.5 m/s，內流1.5 m/s的工況下管體順流向與橫流向的位移時程空間分布。由圖6可知，此工況下管體的振動情況與內流為0.5 m/s的工況相似，橫流向和順流向振動都表現出駐波主導，局部行波間歇性出現的模式。在t/T=0.50～0.60的時間段內，橫流向振型由6階轉化成了5階，多模態轉變現象間歇性出現。圖7為在外流Vmax=1.5 m/s，內流2.0 m/s的工況下管體順流向與橫流向的位移時程空間分布。由圖7可知，此工況下管體橫流向和順流向的振動情況均為駐波主導，雖然存在較少的行波成分，但不足以導致模態階數的轉變。總的來說，管體的內流流速增加可以明顯改變管體的振型和振幅，并且使其振動響應特性更加復雜化。

圖6 內流Vin=1.5 m/s 時，無量綱振動位移時程空間分布

3.2 順流向的平均偏移及最大變形處的振動軌跡

外流作用下，管體受到拖曳力作用，在順流向上會產生比較大的靜態變形，同時又因為交替變更的升力作用和瞬變的拖曳力分量，產生比較顯著的局部振動變形。柔性管體結構的非線性振動必將影響管內流的運輸路徑。管內流路徑的變化，又必然導致沿管體的科氏力、離心力、附加質量力發生變化，使管體的非線性振動響應更加復雜化。為了更全面立體地揭示大變形輸流管體的振動過程，這里將研究內流流速的變化對于彈性管最大偏移位置及相應運動軌跡的影響。圖8為外流為1.0 m/s，內流分別為0.0 m/s,0.5 m/s，1.0 m/s，1.5 m/s和2.0 m/s的5個工況下管體順流向的平均偏移，由圖可知，管體順流向最大平均偏移量在5D～6D之間，最大偏移位置為y/L=0.375。內流流速的增大對管體順流向最大平均偏移量有一定的影響。圖9選擇性地展示了外流為1.0 m/s，內流分別為0.5 m/s和1.5 m/s時最大偏移位置y/L=0.375處的管體運動軌跡，其中橫坐標為管體順流向上各點的相對位置，縱坐標為管體橫流向上各點的相對位置。由于在開始一段時間內兩端約束的柔性管體受到水流的沖擊會產生一個較大的順流向偏移，圖中數據只截取管體產生穩定靜態偏移后的運動軌跡。由圖9可知，每個周期的“8”字形軌跡發生了變形，且相鄰周期的運動軌跡可能會產生一個較大的位置跳躍，體現了大長徑比管體在較大外流作用下強烈振動的過程。由于各個周期內的“8”字形軌跡變形程度不一，且振動位置跳躍較大，所以多個運動周期疊加之后的軌跡規律較不清晰。這也是細長輸流管在內外流共同作用下多模態交涉的體現。

圖8 外流Vmax=1.0 m/s時，不同內流速度Vin下順流向的平均偏移

圖9 y/L=0.375處管體運動軌跡

圖10為外流為1.5 m/s，內流分別為0.0 m/s,0.5 m/s，1.0 m/s，1.5 m/s和2.0 m/s的5個工況下管體順流向的平均偏移，與圖8對比可知，在較大的外流作用下，管體順流向最大平均偏移量高達13D～15D，且內流流速對最大平均偏移量的影響更加顯著，但是最大偏移位置y/L=0.375幾乎沒有發生變化。圖11選擇性地展示了外流1.5 m/s，內流0.5 m/s和1.5 m/s時最大偏移位置y/L=0.375處的管體運動軌跡。由圖11可知，管體最大偏移位置的運動軌跡在多個運動周期疊加下也變得混亂，反映了管體的振動能量在內外流共同作用下變化顯著。這種顯著的振動能量交換現象也與3.1小結中管體多模態共存和振動模態躍遷息息相關。

圖10 外流Vmax=1.5 m/s時，不同內流速度下順流向的平均偏移

圖11 y/L=0.375處管體運動軌跡

3.3 二次流的渦結構演變

由于扭曲的運輸路徑彎曲和不斷振動的管體，內流在柔性管體中運輸，并隨著管體不斷振動的彎曲和旋轉等因素及其相互交涉使管體內部產生復雜的二次流效應[21]。單獨提取管體內流場數據，對管跨中截面流向渦的時程變化以及沿管長不同截面處的流向渦結構進行比較分析。從圖12可知，在t=0.5～3.0 s的時間范圍內，流向渦一開始附著在管壁上，隨著管體的持續振蕩，流向渦開始脫落，并在內部區域以渦對或單渦結構形式存在并運動。隨著時間的推進，渦的強度、渦對數量以及渦對位置都在發生變換，這對柔性管體的渦激振動起到了顯著的作用。最后可以看到大量離散的渦結構分布在各個位置，這些渦可能是由本身從管壁處脫落，也可能是管體的高頻振動導致渦的分離。對于自由振動的柔性管體，橫流向、順流向以及軸向的彎曲和振動都會對內流場結構產生顯著影響。

圖12 管中跨截面內部渦結構的時間演變歷程

4 結 語

基于CFD方法，對內外流作用下的柔性管體渦激振動流固耦合特性開展數值模擬研究。首先參考不含內流的渦激振動試驗工況開展驗證計算，數值結果與試驗結果吻合良好，驗證了數值方法的可靠性。為研究不同內流流速下大變形柔性管體的渦激振動響應，將不同工況計算得出的標準化振動位移的時—空分布、順流向平均偏移、振動軌跡及內流場橫流向渦的演化進行分析，得出了以下結論：

1)外流速度的增大導致更高的振動模態在管體橫流向和順流向響應中激發出來。隨著內流流速的增加，主振模態階數不變，但振動強度會發生明顯改變，這表明了內流的運動帶來的向心力、科氏力和附加質量力使整個管體系統的振動更加復雜。

2)輸流管體的振動表現出多模態響應特性，除了主振模態，其他潛在的固有模態也被激發出來，并具備一定的振動能量。駐波主導了管體兩個正交方向的振動，行波成分易于在順流向上局部出現，這些都體現出管體渦激振動的模態轉變和躍遷。內流流速對順流向最大偏移位置有一定的影響，在多模態響應的影響下，該位置的運動軌跡有明顯的變形和跳躍。

3)觀察到內流場的二次流效應，其中流向渦沿著管壁不斷生成并脫落，形成對渦或者單渦。隨著持續性的局部振動變形，大量離散的渦被發現，這是渦激振動對管體內流場結構的影響。

文中主要研究了有限長彈性管模型在內外流作用下的振動響應問題，為一定長徑比彈性管體的流固耦合分析提供了研究手段，分析結果可為較大長徑比的彈性管渦激振動預報提供參考，但受限于當前數值計算能力，應用于工程實際中大長徑比彈性管渦激振動計算模擬尚有一定難度。