基于動網(wǎng)格的鐵路沿線孔板式沙障流固耦合數(shù)值模擬

智凌巖,程建軍,2,王 連,辛林桂

(1.石河子大學水利建筑工程學院,新疆石河子 832003;2.中鐵西北科學研究院有限公司,蘭州 730000)

中國鐵路沿線風沙災害十分嚴重，風沙災害嚴重威脅鐵路列車的安全運行[1]。為了減少風沙對鐵路的危害，鐵路沿線修建了各種擋風沙構筑物。這些擋風沙構筑物均為厚重的混凝土構筑物(圖1(a))，諸如斜插板式沙障、擋風墻[2]以及防沙工程系統(tǒng)[3]。此類厚重式擋沙墻自身質量大且體積龐大，不會因風荷載作用而發(fā)生傾倒或折斷破壞，但由于施工困難且污染環(huán)境逐漸被新型的沙障結構所代替。而孔板式沙障正是適應新時期要求出現(xiàn)的新型沙障結構(圖1(b))，與傳統(tǒng)厚重式風沙構筑物相比，以孔板式沙障為代表的新型沙障結構具有更好的阻風沙功效。但采用孔板式沙障必須考慮其承受風荷載的能力，以及其在風荷載作用下的力學穩(wěn)定性問題。

圖1 擋風沙構筑物

1 流固耦合計算方法

1.1 流固耦合原理

流固耦合問題是流體力學(Computational Fluid Dynamics，CFD)與固體力學(Computational Solid Mechanics，CSM)交叉而生成的一門力學分支，同時也是多學科或多物理場研究的一個重要分支，它是研究可變形固體在流場作用下的各種行為以及固體變形對流場影響這二者相互作用的一門科學。

流固耦合問題可以理解為既涉及固體求解又涉及流體求解，兩者又都不能被忽略的模擬問題。它是一個多域問題[4-7]，流固耦合可以有效節(jié)約分析時間和成本，同時保證結果更接近于物理現(xiàn)象本身的規(guī)律。

1.1.1 流體控制方程

流體流動要遵守基本的守恒定律，包括質量守恒定律、動量守恒定律、能量守恒定律。對于一般的牛頓流不考慮傳熱問題，其守恒定律通過如下控制方程描述。

質量守恒方程

(1)

動量守恒方程

(2)

式中，t為時間；ff為體積力矢量；ρf為流體密度；v為流體速度矢量；τf為剪切力張量，可表示為

τf=(-p+μ·v)I+2μe

(3)

1.1.2 固體控制方程

固體部分的守恒方程可以由牛頓第二定律導出

(4)

1.1.3 流固耦合方程

流固耦合也遵循最基本的守恒原則，在流固耦合交界面處，應滿足流體與固體應力(τ)、位移(d)等變量的相等或守恒，即滿足如下方程

(5)

注：下標f表示流體，s表示固體。

1.2 動網(wǎng)格技術

在孔板式沙障流固耦合計算過程中，沙障受到風的壓力荷載會發(fā)生變形，沙障的變形也會傳遞給流場域中的流固耦合面，導致網(wǎng)格變形，流體域網(wǎng)格的變形與更新要應用動網(wǎng)格技術。四面體網(wǎng)格較六面體網(wǎng)格更易于網(wǎng)格運動的更新，因此目前動網(wǎng)格技術大部分應用四面體網(wǎng)格[8]。

動網(wǎng)格計算中網(wǎng)格的運動更新過程可以用多種模型計算，針對孔板式沙障流固耦合計算中流體域網(wǎng)格的運動更新，在保證計算收斂和較少計算時間下，選擇彈簧光順模型[9-11]。

在彈簧光順模型中，網(wǎng)格的邊被理想化為節(jié)點間相互連接的彈簧。移動前的網(wǎng)格間距相當于邊界移動前由彈簧組成的系統(tǒng)處于平衡狀態(tài)。以節(jié)點位移為自變量，依據(jù)胡克定律，經(jīng)迭代計算得到使各節(jié)點上的合力等于零，即新的網(wǎng)格節(jié)點。即

(6)

1.3 模型網(wǎng)格及前處理

1.3.1 模型介紹

孔板式沙障孔隙率均為50%，按孔徑大小建立4個三維模型，孔徑分別為9.97、7.98、6.65、5.70 cm。

沙障面板厚為2 cm，高度為100 cm，長度為200 cm；立柱高為105 cm，長寬均為4 cm。沙障底部有5 cm的空隙?？讖綖?.97 cm的沙障模型如圖2所示。

圖2 沙障模型示意(單位：cm)

計算域高為20 m，寬度為5 m，長度為100 m。沙障放置在距入口40 m處?？讖綖?.97 cm的沙障計算域如圖3所示。

圖3 計算域示意(單位：cm)

1.3.2 網(wǎng)格劃分

流固耦合計算時流體域與固體域要分開計算，所以網(wǎng)格劃分時流體域與固體域要分開劃分。

劃分網(wǎng)格方法采用Tetrahedrons，通過尺寸控制對局部區(qū)域進行網(wǎng)格加密，保證較好的網(wǎng)格質量。流體域網(wǎng)格劃分時對流固耦合交界面進行加密，沙障固體域網(wǎng)格劃分時進行全局加密?？讖綖?.97 cm的沙障網(wǎng)格劃分結果如圖4所示。

圖4 網(wǎng)格劃分結果

1.3.3 參數(shù)設置

流固耦合采用的是System Coupling進行雙向流固耦合。

流體域計算邊界條件：根據(jù)空氣動力學原理，當馬赫數(shù)小于0.3時空氣流為不可壓縮流，風沙兩相流馬赫數(shù)均小于0.3[12]，故計算模型入口邊界條件為Velocity-Inlet(速度入口)；自由出流必須在流態(tài)充分發(fā)展條件下才能采用，而此模型出口不能確保為自由出流，故模型出口邊界條件為Pressure-outlet(壓力出口)，其壓差為零；因所計算的物理外形以及所期望的流動具有鏡像對稱的情況，為減小計算量且保證計算結果的準確性模型左右兩側均為Symmetry(對稱邊界條件)；模型上部邊界依據(jù)現(xiàn)實情況選取Pressure-outlet(壓力出口)；下底面和流固耦合交界面采用wall。

流體域計算模型選用標準湍流模型，湍流強度I=0.05，湍流半徑R=1 m。并選取Syamlal-O’Brien曳力模型。方程組求解計算方法采用SIMPLEC算法。模擬風速選取6 m/s[13-14]。

固體結構模型材料為Structural Steel，設定固體結構模型與流體域的接觸面為流固耦合面。在立柱底部添加Fixed Support。流固耦合計算過程中設定Co-Sim，Sequencez中的Sequence:Fluid Flow為1；Transient Structural為2，即設定雙向流固耦合開始的計算順序為先計算流體再計算固體。

2 結果分析

本文對相同孔隙率不同孔徑的4個沙障模型進行了雙向流固耦合數(shù)值模擬，以此來分析孔徑大小變化對流體域流場的影響以及對沙障受力、變形位移的影響。

2.1 流體域流場結果分析

入口風速為6 m/s，通過三維流固耦合數(shù)值模擬結果來觀察孔徑大小變化對流場的影響。圖5為不同孔徑大小沙障的流場俯視云圖，選取俯視云圖高度為1.1 m。

圖5 不同孔徑沙障流場俯視云圖

由不同孔徑沙障流場域云圖可知，孔板式沙障不同于不透風的厚重式擋沙墻，其障后無渦流區(qū)，在障后有大面積的減速區(qū)。因沙障孔隙率較大，開孔較多，并且在底部留有一定空隙，所以在障后無法形成渦流區(qū)，僅對來流有較大的減速效果。并且由模擬結果發(fā)現(xiàn)，障后速度最低區(qū)域并不是在障后靠近沙障位置，而是在障后距沙障較遠位置。

孔徑不同沙障對來流的減速效果也不同。隨孔徑的減小，沙障對來流的削弱效果越強。障后速度最低區(qū)域的速度大小也隨孔徑減小而降低，但通過結果發(fā)現(xiàn)該區(qū)域速度最低降至0。由此可知，孔徑越小沙障在障后區(qū)域對來流的減速效果越強。

沙障孔徑變化對流場域的削弱效果僅表現(xiàn)在障后區(qū)域，在障前區(qū)域隨孔徑變化對障前流場影響并不明顯。沙障流場側視云圖如圖6所示。

圖6 孔徑9.97 cm沙障流場側視云圖

由流場的側視云圖可知，在障前會有大面積的減速區(qū)域，但由于沙障孔隙率大，所以對來流的阻礙作用不強僅體現(xiàn)在對來流在障后背風側有較好的削弱減速效果。模擬相同孔隙率不同孔徑的沙障，分析結果可知，不同孔徑沙障對來流的阻礙作用即障前的減速區(qū)相同，所以可知孔徑變化對來流在障前受到的阻礙作用無影響。

2.2 受力結果分析

由雙向流固耦合結果分析沙障固體結構所受到的最大剪力可知，沙障立柱受到的剪力較大，主要集中在立柱中下部，且越靠近底部剪力越大，并在底部區(qū)域有最大值；沙障面板所受到的剪力主要分布在立柱兩側，且越靠近立柱底部面板所受到的剪力越大、范圍越廣。在面板中間部位也因變形過大受到較大的剪力?？讖?.97 cm沙障固體模型受力云圖如圖7所示。

由模擬結果云圖可知，沙障模型的立柱為受力最大構件，即立柱最有可能發(fā)生受力破壞。所以現(xiàn)以立柱為研究對象，研究孔徑變化對立柱受力特點的影響。

標尺單位：10-5N/cm2圖7 孔徑9.97 cm沙障固體模型受力云圖

模擬結果發(fā)現(xiàn)不同孔徑沙障立柱和面板受力的變化趨勢大致相同，僅在受力大小的數(shù)值上有所差別。說明孔徑的變化不影響沙障的受力分布只影響受力的大小。

圖8為孔徑9.97 cm沙障立柱受力隨高度變化的曲線。分析曲線可知，立柱受到的最大剪力隨立柱高度的增加會在底部出現(xiàn)急劇增大，而后隨高度增加其最大剪力值會逐步減小，最終達到穩(wěn)定狀態(tài)。即立柱受到風力荷載的最大值不是在立柱底部，而是距柱底有一定距離。

圖8 孔徑9.97 cm沙障立柱受力隨高度變化曲線

表1為不同孔徑沙障立柱所受到的最大剪力距立柱底端的位置及最大值。由表1數(shù)據(jù)可知，立柱受到的最大剪力隨孔徑減小而增大，并且其最大值位置距柱底距離穩(wěn)定在(4.5±0.025) cm。說明沙障孔徑的變化不影響其最大剪力的位置，只影響其最大值的數(shù)值，孔徑越小沙障立柱所受到的力越大。

表1 不同孔徑沙障立柱受力最大剪力及位置

2.3 位移結果分析

沙障在風力作用下會發(fā)生變形，其變形最大區(qū)域為孔板的面板中間位置。因為立柱截面尺寸為正方形，相比孔板尺寸比較規(guī)則而且厚度也大于孔板厚度，所以立柱產(chǎn)生的位移小于孔板的位移。如圖9所示。

由圖9可知，每片沙障孔板面板的變形大致呈“U”形分布，其中心線頂部位置變形最大。立柱底部有固定約束所以立柱底部無變形，上部變形次于孔板最大變形，但也有較大變形?？装逯行木€與立柱中心線變形特征線如圖10所示。

圖9 孔徑9.97 cm沙障變形位移云圖

圖10 孔徑9.97 cm沙障面板立柱中心位移特征曲線

由圖10可知，沙障面板和立柱中心位移隨時間的變化特征，明顯看出沙障面板位移各時間均明顯大于立柱位移。在剛開始時沙障面板與立柱位移會出現(xiàn)瞬時最大值，然后在較短時間內位移會出現(xiàn)反復，稱此現(xiàn)象為“沖擊效應”，最終隨時間增長位移穩(wěn)定在固定值。

出現(xiàn)該現(xiàn)象的原因是由于剛開始時沙障處于靜止穩(wěn)定狀態(tài)，風荷載壓力值首次作用到沙障結構時，此時沙障受到的風荷載壓力值為最大值，沙障面板和立柱會發(fā)生瞬時較大位移，此時出現(xiàn)沖擊效應最大值。沙障位移的產(chǎn)生反作用于流體域流場，使流場發(fā)生變化，即作用在沙障上的風荷載壓力值減小。沙障因風荷載壓力值的變化沙障本身會發(fā)生擺動。經(jīng)過流體域與沙障固體域的反復耦合作用，其沙障的沖擊效應初步消退最終達到穩(wěn)定狀態(tài)，沙障面板與立柱位移達到穩(wěn)定。

沙障孔徑的變化并不影響沙障面板與立柱位移變化的特征，僅對其初始沖擊效應最大位移與穩(wěn)定狀態(tài)下的位移有較大影響。如表2所示。

表2 不同孔徑沙障面板與立柱受力特征值

由表2數(shù)據(jù)可知，隨沙障孔徑的減小，其面板和立柱沖擊效應最大值和穩(wěn)定值也有所增大。說明相同孔隙率下孔徑的減小即孔數(shù)量的增加，使沙障與流體的接觸面積增加，使沙障對來流有更大的削弱作用，同時沙障受到來流沖擊作用和風荷載越大。

3 結論

基于動網(wǎng)格技術，通過對相同孔隙率不同孔徑的孔板式沙障進行三維雙向流固耦合數(shù)值模擬，通過對結果的認真分析對比得出以下結論。

(1)通過三維雙向流固耦合對不同孔徑的孔板式沙障進行了更為真實的模擬，為阻風沙構筑物的研究提供新的思路和方法，并對孔板式沙障孔徑的設計提供了力學基礎。

(2)流場結果表明，其障后無渦流區(qū)但有大面積的減速區(qū)域；沙障孔徑的大小對障前會出現(xiàn)的減速區(qū)域無影響；障后大面積減速區(qū)域隨孔徑減小其減速區(qū)面積越大，對來流的減速效果越強，其減速區(qū)速度最低降至0，并未出現(xiàn)渦流區(qū)域。

(3)不同孔徑沙障在風荷載作用下的受力分布特征相同，僅在最大值數(shù)值上有所不同。在風荷載作用下沙障立柱底部區(qū)域受到的力較大，并隨高度增加受力減??；其面板上所受到的力相較于立柱受的力較小。不同孔徑情況下立柱受力最大位置均穩(wěn)定在距柱底距離(4.5±0.025) cm，但其最大值隨孔徑減小而增大。

(4)沙障在風荷載作用下的變形主要表現(xiàn)在面板中心線上部和立柱上部，并且面板產(chǎn)生位移較于立柱產(chǎn)生的位移更為顯著。沙障在來流首次與其接觸時會產(chǎn)生“沖擊效應”，在沙障與來流經(jīng)過反復耦合作用最終“沖擊效應”消退，其面板與立柱產(chǎn)生的位移達到穩(wěn)定狀態(tài)。

[1] Zhang J P， Wang Y S， Jiang F Q. Numerical analysis on the features of sand flow movement around the embankment of Lan-Xin railway in Gobi region[J]. China Railway Science， 2011，32(4):14-18.

[2] Cheng J， Lei J， Li S， et al. Disturbance of the inclined inserting-type sand fence to wind-sand flow fields and its sand control characteristics[J]. Aeolian Research， 2016，21：139-150.

[3] Cheng， J. J.， Xue， C. X. The sand-damage-prevention engineering system for the railway in the desert region of the Qinghai-Tibet plateau[J]. Journal of Wind Engineering & Industrial Aerodynamics， 2014，125:30-37.

[4] 徐楓.結構流固耦合振動與流動控制的數(shù)值模擬[D].哈爾濱：哈爾濱工業(yè)大學，2009.

[5] 李樹云，譚志洪，姜洋，等.基于流固耦合的濾袋振動的數(shù)值模擬[J].環(huán)境工程學報，2016，10(5):2535-2540.

[6] Paik K J， Carrica P M. Fluid-structure interaction for an elastic structure interacting with free surface in a rolling tank[J]. Ocean Engineering， 2014，84(7):201-212.

[7] Jaiman R， Geubelle P， Loth E， et al. Combined interface boundary condition method for unsteady fluid-structure interaction[J]. Computer Methods in Applied Mechanics and Engineering， 2011，200(1-4):27-39.

[8] 閆清東，劉博深，魏巍.基于動網(wǎng)格的沖焊型液力變矩器流固耦合分析[J].華中科技大學學報(自然科學版)，2015，43(12):37-41.

[9] Stein K，Tezduyar T E， Benncy R. Mesh moving techniques for fluid-structure interactions with large displacements[J]. Auris Nasus Larynx，2003，70(1):58-63.

[10] Stein K，Tezduyar T E， Benncy R.Automatic mesh update with the solid-extension mesh moving technique[J]. Computer Methods in Applied Mechanics and Engineering， 2004，193(21):2019-2032.

[11] Pei J，Yuan S， Yuan J，et al. The influence of the flow rate on periodic flow unsteadiness behaviors in a sewage centrifugal pump[J]. Journal of Hydrodynamics， 2013，25(5):702-709.

[12] 程建軍，龐巧東.戈壁強風區(qū)擋風構筑物限制下列車氣動力學特性分析[J].鐵道標準設計，2013，57(1):1-4.

[13] 辛國偉，程建軍，景文宏，等.來流廓線對風沙流場和風沙堆積影響的數(shù)值模擬—以擋沙墻為例[J].干旱區(qū)研究，2016，33(3):672-679.

[14] 辛國偉，程建軍，王連，等.鐵路沿線地表條件與風沙流場的互饋規(guī)律研究[J].鐵道標準設計，2016，60(9):22-27.