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

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

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

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

圖1 擋風(fēng)沙構(gòu)筑物

1 流固耦合計(jì)算方法

1.1 流固耦合原理

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

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

1.1.1 流體控制方程

流體流動(dòng)要遵守基本的守恒定律，包括質(zhì)量守恒定律、動(dòng)量守恒定律、能量守恒定律。對(duì)于一般的牛頓流不考慮傳熱問(wèn)題，其守恒定律通過(guò)如下控制方程描述。

質(zhì)量守恒方程

·(ρfv)=0

(1)

動(dòng)量守恒方程

·(ρfv-τf)=ff

(2)

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

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

(3)

1.1.2 固體控制方程

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

·σs+fs

(4)

1.1.3 流固耦合方程

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

(5)

注：下標(biāo)f表示流體，s表示固體。

1.2 動(dòng)網(wǎng)格技術(shù)

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

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

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

(6)

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

1.3.1 模型介紹

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

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

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

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

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

1.3.2 網(wǎng)格劃分

流固耦合計(jì)算時(shí)流體域與固體域要分開(kāi)計(jì)算，所以網(wǎng)格劃分時(shí)流體域與固體域要分開(kāi)劃分。

劃分網(wǎng)格方法采用Tetrahedrons，通過(guò)尺寸控制對(duì)局部區(qū)域進(jìn)行網(wǎng)格加密，保證較好的網(wǎng)格質(zhì)量。流體域網(wǎng)格劃分時(shí)對(duì)流固耦合交界面進(jìn)行加密，沙障固體域網(wǎng)格劃分時(shí)進(jìn)行全局加密。孔徑為9.97 cm的沙障網(wǎng)格劃分結(jié)果如圖4所示。

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

1.3.3 參數(shù)設(shè)置

流固耦合采用的是System Coupling進(jìn)行雙向流固耦合。

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

流體域計(jì)算模型選用標(biāo)準(zhǔn)湍流模型，湍流強(qiáng)度I=0.05，湍流半徑R=1 m。并選取Syamlal-O’Brien曳力模型。方程組求解計(jì)算方法采用SIMPLEC算法。模擬風(fēng)速選取6 m/s[13-14]。

固體結(jié)構(gòu)模型材料為Structural Steel，設(shè)定固體結(jié)構(gòu)模型與流體域的接觸面為流固耦合面。在立柱底部添加Fixed Support。流固耦合計(jì)算過(guò)程中設(shè)定Co-Sim，Sequencez中的Sequence:Fluid Flow為1；Transient Structural為2，即設(shè)定雙向流固耦合開(kāi)始的計(jì)算順序?yàn)橄扔?jì)算流體再計(jì)算固體。

2 結(jié)果分析

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

2.1 流體域流場(chǎng)結(jié)果分析

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

圖5 不同孔徑沙障流場(chǎng)俯視云圖

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

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

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

圖6 孔徑9.97 cm沙障流場(chǎng)側(cè)視云圖

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

2.2 受力結(jié)果分析

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

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

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

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

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

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

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

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

2.3 位移結(jié)果分析

沙障在風(fēng)力作用下會(huì)發(fā)生變形，其變形最大區(qū)域?yàn)榭装宓拿姘逯虚g位置。因?yàn)榱⒅孛娉叽鐬檎叫?，相比孔板尺寸比較規(guī)則而且厚度也大于孔板厚度，所以立柱產(chǎn)生的位移小于孔板的位移。如圖9所示。

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

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

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

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

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

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

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

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

3 結(jié)論

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

(1)通過(guò)三維雙向流固耦合對(duì)不同孔徑的孔板式沙障進(jìn)行了更為真實(shí)的模擬，為阻風(fēng)沙構(gòu)筑物的研究提供新的思路和方法，并對(duì)孔板式沙障孔徑的設(shè)計(jì)提供了力學(xué)基礎(chǔ)。

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

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

(4)沙障在風(fēng)荷載作用下的變形主要表現(xiàn)在面板中心線上部和立柱上部，并且面板產(chǎn)生位移較于立柱產(chǎn)生的位移更為顯著。沙障在來(lái)流首次與其接觸時(shí)會(huì)產(chǎn)生“沖擊效應(yīng)”，在沙障與來(lái)流經(jīng)過(guò)反復(fù)耦合作用最終“沖擊效應(yīng)”消退，其面板與立柱產(chǎn)生的位移達(dá)到穩(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] 徐楓.結(jié)構(gòu)流固耦合振動(dòng)與流動(dòng)控制的數(shù)值模擬[D].哈爾濱：哈爾濱工業(yè)大學(xué)，2009.

[5] 李樹(shù)云，譚志洪，姜洋，等.基于流固耦合的濾袋振動(dòng)的數(shù)值模擬[J].環(huán)境工程學(xué)報(bào)，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] 閆清東，劉博深，魏巍.基于動(dòng)網(wǎng)格的沖焊型液力變矩器流固耦合分析[J].華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版)，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] 程建軍，龐巧東.戈壁強(qiáng)風(fēng)區(qū)擋風(fēng)構(gòu)筑物限制下列車(chē)氣動(dòng)力學(xué)特性分析[J].鐵道標(biāo)準(zhǔn)設(shè)計(jì)，2013，57(1):1-4.

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

[14] 辛國(guó)偉，程建軍，王連，等.鐵路沿線地表?xiàng)l件與風(fēng)沙流場(chǎng)的互饋規(guī)律研究[J].鐵道標(biāo)準(zhǔn)設(shè)計(jì)，2016，60(9):22-27.