基于动网格的铁路沿线孔板式沙障流固耦合数值模拟

智凌岩,程建军,2,王 连,辛林桂

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

中国铁路沿线风沙灾害十分严重，风沙灾害严重威胁铁路列车的安全运行[1]。为了减少风沙对铁路的危害，铁路沿线修建了各种挡风沙构筑物。这些挡风沙构筑物均为厚重的混凝土构筑物(图1(a))，诸如斜插板式沙障、挡风墙[2]以及防沙工程系统[3]。此类厚重式挡沙墙自身质量大且体积庞大，不会因风荷载作用而发生倾倒或折断破坏，但由于施工困难且污染环境逐渐被新型的沙障结构所代替。而孔板式沙障正是适应新时期要求出现的新型沙障结构(图1(b))，与传统厚重式风沙构筑物相比，以孔板式沙障为代表的新型沙障结构具有更好的阻风沙功效。但采用孔板式沙障必须考虑其承受风荷载的能力，以及其在风荷载作用下的力学稳定性问题。

图1 挡风沙构筑物

1 流固耦合计算方法

1.1 流固耦合原理

流固耦合问题是流体力学(Computational Fluid Dynamics，CFD)与固体力学(Computational Solid Mechanics，CSM)交叉而生成的一门力学分支，同时也是多学科或多物理场研究的一个重要分支，它是研究可变形固体在流场作用下的各种行为以及固体变形对流场影响这二者相互作用的一门科学。

流固耦合问题可以理解为既涉及固体求解又涉及流体求解，两者又都不能被忽略的模拟问题。它是一个多域问题[4-7]，流固耦合可以有效节约分析时间和成本，同时保证结果更接近于物理现象本身的规律。

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 动网格技术

在孔板式沙障流固耦合计算过程中，沙障受到风的压力荷载会发生变形，沙障的变形也会传递给流场域中的流固耦合面，导致网格变形，流体域网格的变形与更新要应用动网格技术。四面体网格较六面体网格更易于网格运动的更新，因此目前动网格技术大部分应用四面体网格[8]。

动网格计算中网格的运动更新过程可以用多种模型计算，针对孔板式沙障流固耦合计算中流体域网格的运动更新，在保证计算收敛和较少计算时间下，选择弹簧光顺模型[9-11]。

在弹簧光顺模型中，网格的边被理想化为节点间相互连接的弹簧。移动前的网格间距相当于边界移动前由弹簧组成的系统处于平衡状态。以节点位移为自变量，依据胡克定律，经迭代计算得到使各节点上的合力等于零，即新的网格节点。即

(6)

1.3 模型网格及前处理

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的空隙。孔径为9.97 cm的沙障模型如图2所示。

图2 沙障模型示意(单位：cm)

计算域高为20 m，宽度为5 m，长度为100 m。沙障放置在距入口40 m处。孔径为9.97 cm的沙障计算域如图3所示。

图3 计算域示意(单位：cm)

1.3.2 网格划分

流固耦合计算时流体域与固体域要分开计算，所以网格划分时流体域与固体域要分开划分。

划分网格方法采用Tetrahedrons，通过尺寸控制对局部区域进行网格加密，保证较好的网格质量。流体域网格划分时对流固耦合交界面进行加密，沙障固体域网格划分时进行全局加密。孔径为9.97 cm的沙障网格划分结果如图4所示。

图4 网格划分结果

1.3.3 参数设置

流固耦合采用的是System Coupling进行双向流固耦合。

流体域计算边界条件：根据空气动力学原理，当马赫数小于0.3时空气流为不可压缩流，风沙两相流马赫数均小于0.3[12]，故计算模型入口边界条件为Velocity-Inlet(速度入口)；自由出流必须在流态充分发展条件下才能采用，而此模型出口不能确保为自由出流，故模型出口边界条件为Pressure-outlet(压力出口)，其压差为零；因所计算的物理外形以及所期望的流动具有镜像对称的情况，为减小计算量且保证计算结果的准确性模型左右两侧均为Symmetry(对称边界条件)；模型上部边界依据现实情况选取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个沙障模型进行了双向流固耦合数值模拟，以此来分析孔径大小变化对流体域流场的影响以及对沙障受力、变形位移的影响。

2.1 流体域流场结果分析

入口风速为6 m/s，通过三维流固耦合数值模拟结果来观察孔径大小变化对流场的影响。图5为不同孔径大小沙障的流场俯视云图，选取俯视云图高度为1.1 m。

图5 不同孔径沙障流场俯视云图

由不同孔径沙障流场域云图可知，孔板式沙障不同于不透风的厚重式挡沙墙，其障后无涡流区，在障后有大面积的减速区。因沙障孔隙率较大，开孔较多，并且在底部留有一定空隙，所以在障后无法形成涡流区，仅对来流有较大的减速效果。并且由模拟结果发现，障后速度最低区域并不是在障后靠近沙障位置，而是在障后距沙障较远位置。

孔径不同沙障对来流的减速效果也不同。随孔径的减小，沙障对来流的削弱效果越强。障后速度最低区域的速度大小也随孔径减小而降低，但通过结果发现该区域速度最低降至0。由此可知，孔径越小沙障在障后区域对来流的减速效果越强。

沙障孔径变化对流场域的削弱效果仅表现在障后区域，在障前区域随孔径变化对障前流场影响并不明显。沙障流场侧视云图如图6所示。

图6 孔径9.97 cm沙障流场侧视云图

由流场的侧视云图可知，在障前会有大面积的减速区域，但由于沙障孔隙率大，所以对来流的阻碍作用不强仅体现在对来流在障后背风侧有较好的削弱减速效果。模拟相同孔隙率不同孔径的沙障，分析结果可知，不同孔径沙障对来流的阻碍作用即障前的减速区相同，所以可知孔径变化对来流在障前受到的阻碍作用无影响。

2.2 受力结果分析

由双向流固耦合结果分析沙障固体结构所受到的最大剪力可知，沙障立柱受到的剪力较大，主要集中在立柱中下部，且越靠近底部剪力越大，并在底部区域有最大值；沙障面板所受到的剪力主要分布在立柱两侧，且越靠近立柱底部面板所受到的剪力越大、范围越广。在面板中间部位也因变形过大受到较大的剪力。孔径9.97 cm沙障固体模型受力云图如图7所示。

由模拟结果云图可知，沙障模型的立柱为受力最大构件，即立柱最有可能发生受力破坏。所以现以立柱为研究对象，研究孔径变化对立柱受力特点的影响。

标尺单位：10-5N/cm2图7 孔径9.97 cm沙障固体模型受力云图

模拟结果发现不同孔径沙障立柱和面板受力的变化趋势大致相同，仅在受力大小的数值上有所差别。说明孔径的变化不影响沙障的受力分布只影响受力的大小。

图8为孔径9.97 cm沙障立柱受力随高度变化的曲线。分析曲线可知，立柱受到的最大剪力随立柱高度的增加会在底部出现急剧增大，而后随高度增加其最大剪力值会逐步减小，最终达到稳定状态。即立柱受到风力荷载的最大值不是在立柱底部，而是距柱底有一定距离。

图8 孔径9.97 cm沙障立柱受力随高度变化曲线

表1为不同孔径沙障立柱所受到的最大剪力距立柱底端的位置及最大值。由表1数据可知，立柱受到的最大剪力随孔径减小而增大，并且其最大值位置距柱底距离稳定在(4.5±0.025) cm。说明沙障孔径的变化不影响其最大剪力的位置，只影响其最大值的数值，孔径越小沙障立柱所受到的力越大。

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

2.3 位移结果分析

沙障在风力作用下会发生变形，其变形最大区域为孔板的面板中间位置。因为立柱截面尺寸为正方形，相比孔板尺寸比较规则而且厚度也大于孔板厚度，所以立柱产生的位移小于孔板的位移。如图9所示。

由图9可知，每片沙障孔板面板的变形大致呈“U”形分布，其中心线顶部位置变形最大。立柱底部有固定约束所以立柱底部无变形，上部变形次于孔板最大变形，但也有较大变形。孔板中心线与立柱中心线变形特征线如图10所示。

图9 孔径9.97 cm沙障变形位移云图

图10 孔径9.97 cm沙障面板立柱中心位移特征曲线

由图10可知，沙障面板和立柱中心位移随时间的变化特征，明显看出沙障面板位移各时间均明显大于立柱位移。在刚开始时沙障面板与立柱位移会出现瞬时最大值，然后在较短时间内位移会出现反复，称此现象为“冲击效应”，最终随时间增长位移稳定在固定值。

出现该现象的原因是由于刚开始时沙障处于静止稳定状态，风荷载压力值首次作用到沙障结构时，此时沙障受到的风荷载压力值为最大值，沙障面板和立柱会发生瞬时较大位移，此时出现冲击效应最大值。沙障位移的产生反作用于流体域流场，使流场发生变化，即作用在沙障上的风荷载压力值减小。沙障因风荷载压力值的变化沙障本身会发生摆动。经过流体域与沙障固体域的反复耦合作用，其沙障的冲击效应初步消退最终达到稳定状态，沙障面板与立柱位移达到稳定。

沙障孔径的变化并不影响沙障面板与立柱位移变化的特征，仅对其初始冲击效应最大位移与稳定状态下的位移有较大影响。如表2所示。

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

由表2数据可知，随沙障孔径的减小，其面板和立柱冲击效应最大值和稳定值也有所增大。说明相同孔隙率下孔径的减小即孔数量的增加，使沙障与流体的接触面积增加，使沙障对来流有更大的削弱作用，同时沙障受到来流冲击作用和风荷载越大。

3 结论

基于动网格技术，通过对相同孔隙率不同孔径的孔板式沙障进行三维双向流固耦合数值模拟，通过对结果的认真分析对比得出以下结论。

(1)通过三维双向流固耦合对不同孔径的孔板式沙障进行了更为真实的模拟，为阻风沙构筑物的研究提供新的思路和方法，并对孔板式沙障孔径的设计提供了力学基础。

(2)流场结果表明，其障后无涡流区但有大面积的减速区域；沙障孔径的大小对障前会出现的减速区域无影响；障后大面积减速区域随孔径减小其减速区面积越大，对来流的减速效果越强，其减速区速度最低降至0，并未出现涡流区域。

(3)不同孔径沙障在风荷载作用下的受力分布特征相同，仅在最大值数值上有所不同。在风荷载作用下沙障立柱底部区域受到的力较大，并随高度增加受力减小；其面板上所受到的力相较于立柱受的力较小。不同孔径情况下立柱受力最大位置均稳定在距柱底距离(4.5±0.025) cm，但其最大值随孔径减小而增大。

(4)沙障在风荷载作用下的变形主要表现在面板中心线上部和立柱上部，并且面板产生位移较于立柱产生的位移更为显著。沙障在来流首次与其接触时会产生“冲击效应”，在沙障与来流经过反复耦合作用最终“冲击效应”消退，其面板与立柱产生的位移达到稳定状态。

[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] 徐枫.结构流固耦合振动与流动控制的数值模拟[D].哈尔滨：哈尔滨工业大学，2009.

[5] 李树云，谭志洪，姜洋，等.基于流固耦合的滤袋振动的数值模拟[J].环境工程学报，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] 闫清东，刘博深，魏巍.基于动网格的冲焊型液力变矩器流固耦合分析[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] 程建军，庞巧东.戈壁强风区挡风构筑物限制下列车气动力学特性分析[J].铁道标准设计，2013，57(1):1-4.

[13] 辛国伟，程建军，景文宏，等.来流廓线对风沙流场和风沙堆积影响的数值模拟—以挡沙墙为例[J].干旱区研究，2016，33(3):672-679.

[14] 辛国伟，程建军，王连，等.铁路沿线地表条件与风沙流场的互馈规律研究[J].铁道标准设计，2016，60(9):22-27.