铁路治沙HDPE网的数值模拟计算

程建军，王 瑞，王 连，高 丽，丁录胜，丁泊淞

(1.石河子大学 水利建筑工程学院，新疆 石河子 832003;2.武汉大学 土木建筑工程学院，湖北 武汉 430072；3.新疆铁道勘察设计院有限公司，新疆 乌鲁木齐 830011)

HDPE网(高密度聚乙烯树脂纤维网)作为一种新型的防沙材料，与其他类型的防沙材料相比，具备耐久性持久、防风阻沙效益优良且可大批工业化生产和快速施工等特点，应用前景广阔[1]。随着风沙地区铁路建设的快速推进，大量HDPE网沙障取代了即有的混凝土式挡沙墙和挡沙栅栏。

近年来，针对HDPE网防沙特性的研究逐渐开展起来，屈建军等基于HDPE网的风洞测试发现HDPE网兼有疏透和通风2种形式，其最佳孔隙度为40%～45%，保护区达30h(h为沙障的高度)以上，并通过野外监测验证了HDPE网阻沙、导沙的良好功能[2-3]；张可存等[4]通过对不同孔隙度网的风洞模拟试验，探求输沙量与孔隙度和风速二者之间的关系，最终确定网防沙效益的最佳孔隙度在40%左右；汪万福等[5]通过对敦煌莫高窟顶HDPE网的野外观测发现网的阻沙效益与沙源供给关系密切，沙源充足时以风积为主，相反，则出现较强的风蚀作用。总结以上研究发现，目前针对HDPE网的防风阻沙性能研究主要集中在野外观测和室内风洞测试，都无法得到网式沙障前后的流场演化规律和特征[6-9]。

解决这一问题的另一途径是采用数值模拟计算，而采用数值模拟计算对网的防沙特性进行研究鲜见文献。针对水产养殖方面的渔网，赵云鹏等[10]通过多孔介质模型对水流作用下平面网周围流场的模拟；黄小华、Tang M F等[11-12]基于数学建模的方法对水流作用下网动态变形及受力平衡后空间分布的模拟分析；Bi C W等[13]采用集中质量法和多孔介质模型对渔网的模拟；Chen H等[14]对渔网采用多孔介质模型时阻力系数如何确定进行了研究。以上研究针对水流和渔网的互馈机理，对针对网防风固沙的数值模拟计算研究具有借鉴意义。

针对防风固沙网的数值模拟计算有以下几处难点: ①若建立与野外网式沙障同样尺度的计算模型，则网模型的网格划分及后期计算都较为困难；②采用流体动力学数值模拟计算方法进行网的数值计算又存在诸多设定参数需要先进行风洞测试，而在这方面工作尚未开展；③针对具体的网模拟计算问题，有多种模型可供选择，针对防风阻沙用途的网模型适用性问题尚未解决。

本文对HDPE网的流场特性进行了风洞测试，将其作为基准参考。同时在风洞内进行了网的压降规律测试，获取数值模拟计算所需的压降—速度关系参数。然后进行不同物理模型条件下HDPE网的流场数值模拟计算，研究分析各种模型条件下的参数设置、计算关键步骤，以及各种计算模式的合理选择方法。

1 风洞试验

1.1 风洞测试原理及HDPE网模型布置

风洞试验在中国科学院新疆生态与地理研究所莫索湾野外风洞进行。该风洞为直流吹气式风洞，试验段截面1.3 m×1.0 m，试验段长8 m，试验风速0～20 m·s-1连续可调。试验中采用皮托管和微差压变送器组成的测量系统测定风速，同时通过风速模拟系统软件读取风速数据。测量系统由不同高度的(h0=1，2，3，5，7，9，10，30和50 cm)10个精细皮托管组成(图1)。HDPE网模型高为8 cm，选用网的规格为JZSPE800Lu32，开口率约为55%，丝径为0.32 mm，平均孔目尺寸为0.98×1.14 mm，经纬向强度约为900 N/50 mm。

图1 风洞结构及网的布置示意图

风洞试验结果的可靠性，主要取决于风洞模拟条件与野外实际情况的相似程度，要使风洞流场与野外流场之间具备相似性，即满足几何相似、运动相似和动力相似。在几何相似方面，考虑到风洞湍流附面层与地球贴地边界层的尺度之比，为了避免风洞洞壁对流场干扰，模型布置面积要求小风洞横截面的6%，当采用网模型高度比例尺为1∶20时满足要求；运动相似方面，主要考虑流态和风速廓线相似，流态相似要求考虑雷诺数Re，该风洞在试验风速下Re≈105～106，流态为完全湍流，可达到自模拟要求，经试验验证，风洞内风速廓线满足对数分布；动力相似方面，由于风洞中无温度层结控制装置，仅模拟中性层结条件，试验不涉及自由表面流动，弗劳德数Fr相似可不予考虑。因此，本试验满足几何相似、运动相似、动力相似条件。

在特定入口风速条件下，在水平距离上测试网前10h，5h，1h，0.5h，0h及网后0.5h，2h，3h，20h等共12个位置点的风速剖面，每个风速剖面在垂直高度采用毕托管采集1～50 cm范围内的风速值(见图1)。在数值计算中需要设定各项参数，其中网体前后的压降—速度比是HDPE网数值计算必须参数，为了进行有效数值计算，在风洞里进行了压降测试试验，压降测试采用直径20 cm的圆管，HDPE网置于管中，逐级增大风速，测试HDPE网的迎风侧及背风侧的压力值如图2所示。压降测试是为了求解压降公式中的内部惯性阻力系数和黏性阻力系数，流体动力学数值计算压降公式为

(1)

图2 风洞内网的压降测试

1.2 HDPE网测试结果

为与数值计算对比分析，在风洞测试中仅选择10 m·s-1为入口风速，为突出网体前后局部区域速度变化，在风洞测试数据中仅提取水平距离上网前3h和网后5h范围以内的风速剖面，风速剖面高度范围在1～30 cm内进行数据采集，绘制局部流场如图3所示。

图3 HDPE网的流场分区(局部)

风洞内HDPE网的压降测试选择6，8，10，12，14和16 m·s-1共6组风速，分别测定网前及网后0.1 cm处的压力值，采用多项式拟合方法处理得到不同风速下压降关系，同时对压降数据进行曲线拟合，得到的多项式方程为Δp=0.215 9v2+0.846 9v，相关系数R值0.965 2，表明拟合效果较好(图4)。利用系列压降数据并结合多孔介质模型理论经典式(1)可推算出得内部惯性阻力系数C2为1 305.39；黏性阻力系数1/α为175 298 982.9，这2个参数是本次试验所采用的HDPE网的流体力学性能参数，在后续数值模拟计算中作为输入参数，其数值模拟计算结果才真实可靠。

图4 HDPE网的压降与速度关系

2 HDPE网数值模拟计算方法

2.1 数值模拟计算原理

在运用数值计算方法对一些实际问题进行模拟计算时，设置合理的工作环境、边界条件以及选择恰当的算法是正确计算的前提条件，其中算法的选择对模拟计算效率起到至关重要的作用。目前，对针对复杂工程计算问题，区域离散化算法不失为可靠的选择。区域离散化即采用合理的网格划分方法产生1组有限个离散的点(单元)代替原来连续的空间，常用的离散化方法有有限差分法、有限单元法和有限体积法，以下采用基于完全非结构化网格的有限体积法进行模拟计算。

在数值模拟计算中，必须满足数值计算应遵循的方程组即质量守恒方程与动量守恒方程，同时必须确定κ-ε两方程的计算系数。针对网的数值模拟计算需要进一步对动量守恒方程相关参数进行推算。HDPE网的数值模拟计算中，需要特别阐述动量守恒方程中源项意义，动量守恒方程如下。

(2)

式中：τij为应力张量；p为压力；gi和Fi分别为i方向上的重力体积力和外部体积力。

式(2)中，Fi称为相关计算模型的源项，源项由黏性损失项和惯性损失项2个部分组成，在HDPE网的数值模拟计算中Fi即为式(1)的压降值。

对于κ-ε两方程中常用的5个紊流计算系数，Cμ，C1ε，C2ε，σk和σε，大量气固流体文献显示，其取值较为一致，分别0.09，1.44，1.92，1.0和1.3。κ-ε两方程属于高雷诺数积分到壁面的不可压缩湍流涡黏性模式，要求第1层网格位于湍流充分发展区域，故不必对边界进行网格加密处理，控制方程组采用有限体积法。扩散项采用中心差分格式，对流项采用2阶迎风差分格式，采用分离隐式求解器求解代数方程组，SIMPLEC算法耦合连续性方程和动量方程[15]。

2.2 模拟计算关键设置步骤

在数值计算中，自由边界计算区域的确定对计算结果有很大影响。计算区域选取过小，则无法真实反映实际的开放无界信息，导致计算误差大甚至失真；而选取足够大的计算区域，就会导致计算量庞大，大量数值计算研究证明，过流断面在满足阻塞率小于6%时，模拟计算较为符合真实的流动情况。本次模拟计算为满足过流断面阻塞率要求，流体计算域选用风洞断面尺寸一致，将计算域指定为8 000 mm×1 300 mm×1 000 mm的长方体区域。

选择基于压力的求解器进行求解，而压力—速度的耦合方法应该选择适用于定常不可压缩流体且收敛性较好的SIMPLEC方法，湍流具体表示方法采用流体动力学常用的湍流强度和水力直径表示。在空间离散化方法选择上面，为使结果精度较高，用1阶迎风格式计算收敛后作为初始值，然后采用2阶迎风格式进行计算。

2.3 HDPE网数值模拟计算方法

2.3.1 多孔介质模型方法

在高雷诺数条件下的三维模型进行数值计算是比较困难的，要求网格节点数应在不影响计算速度的情况下尽量多，即网格划分要尽可能的密集，从而提高差值型函数阶数，以达到提高精度的目的。当入口风速为10 m·s-1时，模拟计算显示，在计算域网格节点数不少于104个时，计算精度可达到要求。划分网格前，在DESIGNMODELER(简称DM)中，通过切片操作将计算区域进行划分，定位出多孔介质区域，进行适当的加密处理，以更好地进行网格生成，如图5所示。

图5 DM建模示意图

模拟计算选择标准κ-ε两方程，压力—速度耦合选择定常不可压缩且收敛性较好的SIMPLEC方法。对于多孔介质区域的计算参数设置，通过网的压降测试数据获得，黏性阻力系数1/α为175 298 982.9，内部惯性阻力系数C2为1 305.39。计算模型入口选择速度入口条件(velocity-inlet)，出口选择压力出口条件(pressure-outlet)，由于计算区域满足阻塞率的要求且采用廓线风速入口条件，故其他边界均设为滑移壁面条件。绝对收敛标准小于10-3，并选择入口初始化条件，以满足在开始计算时，每个网格均从以入口条件的状态进行计算。计算收敛结束后，利用TECPLOT进行后处理，速度局部流场图如图6所示。

图6 多孔介质模型条件下的流场

从图6可以看出：在入口廓线风速条件下，气流经过以多孔介质模型代替的实体网后，流场变化较为明显；从流线的变化看，由于多孔介质区域内存在内部阻力和黏性阻力，仅有部分气流可以通过多孔区域，多孔区域的迎风侧流线发生挤压变化，越过多孔区域的气流会有一定速度的增加，并且由于气流向上抬升，会在多孔介质区域后产生回流，从多孔区域透过的气流会抵消这部分回流，导致涡流区消失，即网体将来流分为渗流风和绕流风两部分；多孔介质模型的设置，使气流在网体后变化较大，速度得到较大的削弱，但网体后速度梯度较小。

2.3.2 压力跳跃模型方法

压力跳跃模型建立过程中，通过创建一矩形面域以承载压力跳跃条件。采用四面体网格(Tetrahedrons)对计算区域进行网格划分，网格尺寸设为100 mm,并对压力跳跃面域进行网格加密处理，设置其网格尺寸为10!mm，如图7所示。

图7 采用四面体方法划分的计算域网格

模拟仍计算选择标准κ-ε两方程模型，压力—速度耦合选择定常不可压缩且收敛性较好的SIMPLEC方法。区域条件不需设置，多孔跳跃模型会在边界条件生成baffle条件，此时，编辑baffle，设置压力跳跃系数为1 305.39。计算模型入口选择速度入口条件(velocity-inlet)，出口选择压力出口条件(pressure-outlet)，调用同样以10 m·s-1为梯度风速的廓线风速，其他壁面均设为滑移壁面条件。控制计算时监视残差标准不大于10-3，选择从入口进行初始化计算。模拟计算收敛后，通过TECPLOT进行后处理，其部分流场图如图8所示。

图8 压力跳跃模型条件下的流场

从图8可以看出：压力跳跃模型的设置，同样会使廓线气流发生明显的变化，对实体网的替代具有一定的可行性；压力跳跃面域的存在，阻碍了部分气流的稳定流动，导致压力跳跃面前部分流线发生变化，越过矩形面域顶部后，形成一定的加速区，压力跳跃面域后的部分回流与从压力跳跃面域透过的气流混合后，形成较为均匀的低风速区，越接近底部速度越小。

2.3.3 实体网模型

利用CAD阵列方法建立实体网，如图9所示。

图9 CAD建立的实体模型

实体模型采用便于数值计算的截面为1 mm×1 mm的长方体代替实际的网的经纬丝。网体的实体建模困难在于结构的多重性，其次是网格划分的困难性以及计算的耗时性，文中仅模拟1 mm×400 mm×80 mm的单片网，其计算过程耗时超过了其他模型计算的10倍，完成的网格及相关信息见表1。

表1 网格划分相应参数值

计算选择标准κ-ε两方程，压力—速度耦合选择定常不可压缩且收敛性较好的SIMPLEC方法。入口速度依然采用廓线风速形式。残差监视标准和初始化方法同上。模拟计算收敛后，利用TECPLOT软件进行后处理，得到局部速度云图如图10所示。

图10 实体网模型条件下的流场

由图10可以看出：实体网对廓线风速的影响显著，流场变化比较符合实际情况。实体网的作用，使来流在网前出现挤压，网顶部产生加速区，由于4 mm×4 mm均匀孔的存在，使得网体后出现均匀的气流，速度值呈先减小后增加的微梯度变化趋势。

3种模型对来流的作用效果相似，均具有削弱来流强度，整合网后流场的作用。多孔介质模型和压力跳跃模型都实现了HDPE网的数值计算，体现在采用各种模型对HDPE网进行数值计算后，网背风侧都不存在涡流区，这与实际情况相符(图11)，而HDPE网对来流的消减作用与其开孔率有密切的关系，开孔率越大，消减来流程度越小，当网的开孔率发生变化后，需要针对性地做风洞压降试验获取计算参数后再进行对应数值模拟计算。

以上对比计算发现，3种方法对HDPE网的数值计算结果与风洞试验的规律相同，但从网体前后不同位置上速度分布的具体值分析可知，由于在数值计算中均未考虑网体的弹性模量，导致相应位置速度值有一定的偏差，而压力跳跃模型方法与风洞数据的拟合仍较为理想。由图11可见，实体模型方法在相应位置处的速度值均小于其他模型方法，主要因为实体模型模拟计算无法产生变形，即计算过程中默认网体刚度无限大，实体模型网经纬孔的边界层显著地消减来流风速，同样多孔介质区域和压力跳跃面也没有设置动网格，未能反映网的波动对流场产生的影响效应；多孔介质模型方法能实现对多道网的数值模拟计算，能够考虑区域范围内的流场消减效果，故比较适合三维模型的模拟计算；对于单片HDPE网，压力跳跃模型方法与风洞数据拟合较好，因此，仅考虑单道网对流场的控制作用时，可采用压力跳跃模型方法。

图11 不同模型条件下流场内速度对比

3 结 论

(1)HDPE网在物性参数确定条件下，开孔率为55%，平均孔目尺寸为0.98×1.14 mm的物性条件下进行风洞压降测试，获得内部阻力系数为1 305.39；黏性阻力系数为175 298 982.9，为今后进行类似网的数值计算提供了计算依据。

(2)对HDPE网进行数值模拟计算，采用压力跳跃模型方法和多孔介质模型方法所的计算流场都与实际相符合，即网背风侧不存在涡流区，网将来流分割为两部分，一部分为网上的扰流，一部分为透过网的渗透流，网后最终的流场体现在两部分流体的相互作用结果上。

(3)压力跳跃模型方法用于计算分析已知速度—压降的网，计算收敛性好，实现了网后渗透流的降速，网上扰流的增速，与实际情况相符合，对比3种模型可知，压力跳跃模型的跳跃界面对来流的控制与网的作用规律一致。

(4)多孔介质模型方法与压力跳跃模型方法类似，不同之处在于，多孔介质模型多针对三维计算模型，可以对某区域进行参数设置，而压力跳跃模型方法则是二维计算模型，针对网的数值计算，既可以把网看成二维流场中的1个界面，也可以看作为三维计算模型中的1个降速区域，采用2种模型在参数设定争取的前提下，计算结果都真实可靠，针对网的模拟计算，考虑网厚度可采用多孔介质模型，不考虑网厚度可采用压力跳跃模型。

限于篇幅有限，没有考虑网体随流体改变，将在后续研究中加以说明。