露天堆场防风抑尘网临界孔隙率的数值模拟

潘武轩，何鸿展，宋翀芳，李临平（太原理工大学环境科学与工程学院，山西 太原 030024）

露天堆场防风抑尘网临界孔隙率的数值模拟

潘武轩，何鸿展，宋翀芳*，李临平（太原理工大学环境科学与工程学院，山西 太原 030024）

应用Fluent6.3对防风抑尘网及料堆周围流场进行数值模拟，通过研究不同孔隙率抑尘网与料堆周围湍流特性及料堆表面剪切力分布规律，确定了临界孔隙率.结果显示：高孔隙率（0.3、0.4、0.6）流动状态与无网工况一致，来流风沿迎风面贴附向上，风速逐渐增大，堆顶达到最大；低孔隙率（0、0.2）流态与无网迥异，迎风面处于涡旋中，风速向下，孔隙率为0时涡旋中心高于堆顶3m，孔隙率为0.2时涡旋中心位于堆高2/3处.孔隙率≥0.3时，料堆各表面剪切力变化趋势一致，与无网工况来流风同向，最大剪切力出现在堆顶.孔隙率为0.2时，作为最大起尘量的迎风面，其表面剪切力随高度先增大后减小，最大剪切力处于堆高3/5处.孔隙率为0.25时，湍流结构和剪切力分布发生突变，迎风面为贴附和涡旋复合流动，表面剪切力最小.据此确定来流风速6m/s，运用该几何模型时，临界孔隙率为0.25.

防风抑尘网；数值模拟；临界孔隙率；剪切力

扬尘是大气污染的重要来源，相关研究［1-4］显示，PM10中扬尘排放占39.5％，浓度贡献率为48.7％.防风抑尘网是一种利用空气动力学原理，改变网后微环境，有效抑制开放性露天堆场散尘的方法［5-9］，已在大型堆场及港口等得到成功应用［10-16］.

孔隙率是决定防风抑尘网抑尘效率最重要的参数［17］.胡晋明等［18］，采用数值模拟方法分析了不同孔隙率（0、0.195、0.252、0.33、0.44、0.406、0.49、0.54和0.66）下的风速和湍流度，建议孔隙率取0.3～0.5；姜鑫等［19］应用数值模拟对比了21组不同孔隙率抑尘网后平均风速，得出最佳孔隙率为0.35～0.4；Lee等［20］在风洞试验中测定了网后压力系数，得出孔隙率为0.4～0.5时遮蔽效果最好；Dong等［21］借助粒子图像测速系统PIV（Particle Imaging Velocity）测试了11种不同孔隙率抑尘网后平均速度、湍流度和湍动能，得出抑尘网的遮蔽效果与湍流度直接相关的结论；其他学者［22-24］在研究抑尘网的遮蔽效果时同样以风速、压力和湍流度作为标识指标.

可以看出目前关于抑尘网作用效果的研究多集中以平均风速、平均压力系数和水平湍流度等宏观特性为表征，考虑料堆周围湍流结构和料堆表面微观动力学分布的研究较少.本文通过数值模拟分析了低孔隙率（0、0.2）和高孔隙率（≥0.3）工况绕流与渗流耦合作用时，料堆与抑尘网间湍流结构的微观特性，计算了迎风面和背风面0～17m高度、平顶面0～10m宽度的剪切力分布，根据绕流和渗流主导作用的转换探究了兼具涡旋和贴附流特点的临界孔隙率，为抑尘网孔隙率的优化提供动力学理论依据.

1 数值模拟

1.1 计算区域

数值模拟中料堆计算模型取三维长条棱台.棱台参数设置为：上表面长113m，宽10m；下表面长154m，宽51m；棱台高17m.抑尘网设置在堆前1倍堆高处，网长为1倍堆长.相关研究［25-26］表明，防风网高度大于煤堆高度1.5倍时抑尘效果无明显变化，最佳网高在1～1.5倍的堆高范围内选择比较合适，模型中网高取1.3倍堆高即22m.

计算区域的选择直接影响到结果的可靠性.计算区域选择过大，会增加网格数，增大计算量，耗时多；计算区域选择过小，边界会影响料堆表面流动，影响空气流动的充分发展，影响计算精度.相关研究［27］表明，当计算区域取料堆前6倍的堆宽，料堆后取7倍的堆宽，宽度取2倍的堆长，高度取7倍的堆高时风速和料堆表面剪切力均不再变化.所以计算区域取为714m×119m×308m的长方体.x方向沿堆宽方向，y方向沿堆高方向，z方向沿堆长方向，原点建立在计算区域长方体的内侧左下角，如图1所示.

图1 计算区域（m）Fig.1 Computational domain（m）

1.2 控制方程及边界条件

计算区域选定后建立数学物理模型.料堆周围的空气为不可压缩流体，流动是稳态绝热的.任何流动都须遵守质量守恒和动量守恒基本定律，湍流流动附加湍动能方程和耗散率方程.控制方程组如下：

质量守恒方程：

动量守恒方程：

其中：ρ为空气密度，kg/m3；u、v、w分别为x、y、z方向的速度，m/s；ui、uj分别为i方向和j方向速度分量，m/s；xi、xj分别代表x、y方向的坐标；μ为空气动力黏性系数，kg/（m·s）；δij为克罗内克尔张量.

Si是源项，该值在抑尘网区域外为0，在防风抑尘网区域内可以模拟为多孔介质的作用［28］，源项由黏性损失项和惯性损失项两部分组成：

式中： α为多孔介质的渗透系数，m·s；C2为惯性阻力因子，m-1；Af为孔的总面积，m2；Ap为板的总面积，m2；t为抑尘网厚度，m，取0.002m；C近似等于0.98.

采用标准k-ε方程［29］模拟三维流场，k-ε方程由湍动能k方程和耗散率ε方程组成.

对不可压缩流体，Gb=0，YM=0，G3ε=0.

式中：Gk、Gb为湍动能生成项；Prt为湍动Prandtl数；μt为湍流黏性系数，kg/（m·s）；β为热膨胀系数，K-1；gj为重力加速度在j方向分量，m/s2；T为温度，K；YM为可压缩流中脉动膨胀对总体耗散率的贡献；a为导温系数，m2/s.

紊流系数C1ε、C2ε、C3c、σε和σk的取值分别为1.44、1.92、0.09、1.0和1.3.

边界条件设置：入口边界为速度入口，出口边界为零压力梯度出流，料堆表面及地面为无滑移壁面，抑尘网设置为多孔阶跃边界，模型上表面和前后表面为对称边界.

1.3 网格划分

为探索料堆表面微观特性，有必要将计算空间分区，用一组间距1m的平行平面截料堆表面，使得17m高的料堆迎风面和背风面分别由17个平面组成，10m宽平顶面由10个平面组成.本研究旨在获取料堆和抑尘网周围空气流动特性，远端空气流动对料堆起尘影响甚微，故在网格划分时细划料堆表面，采用三角形网格均匀划分；地面则采用四边形网格，体网格采用六面体网格，如图2料堆表面网格较细密，远离料堆区域网格较疏.抑尘网部分作为多孔阶跃边界，在Fluent6.3中多孔阶跃边界是将抑尘网部分作为动量方程的动量损失源项.网格划分见图2.

图2 网格划分Fig.2 Mesh generation

作为获得数值解的网格应当足够的细密，以致于进一步加密网格已经对数值计算结果基本上没有影响［30］.从疏到密对计算区域进行网格划分，计算出每种网格工况下料堆周围的流场特性，取迎风面剪切力这一宏观量进行网格独立性考核，结果见表1.

表1 网格独立性考核Table 1 Assessment of grid independence

由表1可知当网格数大于1860191时，迎风面表面剪切力变化率仅为0.03％.通过网格独立性考核，取有效计算网格数为1860191.由于抑尘网网格划分影响计算精度，在其他网格划分不变情况下，对抑尘网进行网格独立性考核，发现在网格数大于1860191时，迎风面剪切力不随抑尘网网格数变化，所以不影响计算精度.采用SIMPLE（Semi-Implicit Method for Pressure Linked Equations）算法，标准k-ε方程，收敛误差取10-5，差分格式中压力、动量、湍流能量和湍流耗散率均采用二阶迎风格式.

2 结果及讨论

2.1 速度矢量

图3 不同孔隙率对称面料堆前后速度矢量分布Fig.3 Velocity profiles before and behind the pile on the symmetry plane with deferent porosity

应用FLUENT中的压力基非耦合三维求解器，采用标准k-ε方程对料堆周围空气流场进行数值计算.孔隙率分别为0、0.2、0.25、0.3和0.4，来流风速为6m/s时，对称面料堆周围的速度矢量分布如图3所示.

来流空气流经抑尘网时，一部分从上部绕流，另一部分渗流过网孔.图3显示：低孔隙率（0，0.2）和高孔隙率（0.3，0.4）流动特性截然不同.低孔隙率空气绕流作用大于渗流作用，抑尘网与料堆之间涡旋强度大，紧贴迎风面空气薄层在涡旋作用下速度向下；高孔隙率渗流作用增强，抑尘网与料堆之间无涡旋，迎风面为贴附流，速度向上.

如图3（a）所示，孔隙率为0时，抑尘网和迎风面之间有大尺度涡旋，涡旋中心高度为20m，高于堆高3m；这是由于孔隙率为0时无渗流，空气全部绕流，在网上方形成高速剪切层，网前后压差大，在网后形成大尺度涡旋；整个迎风面处于涡旋影响范围，迎风面最大速度出现在堆顶，沿x、y方向分别为-2.81m/s和-3.21m/s；在涡旋影响下，平顶面依然处于涡旋回流区域，速度为负，与来流方向相反.当孔隙率为0.2时，如图3（b）所示，堆前涡旋强度减弱，涡旋中心靠近迎风面，这是由于孔隙率增大后渗流风增加，网前后压力差较孔隙率为0时减小，涡旋强度减弱，涡旋中心高度为12m，位于迎风面2/3高度处；迎风面最大速度出现在其高度的1/3到2/3高度间，沿x、y方向分别为-0.73m/s和-1.09m/s.

当孔隙率为0.3时，如图3（c）所示，堆前涡旋消失，渗流风进一步增加，网前后压差减小至无法形成涡旋，整个迎风面均为贴附流，迎风面最大速度出现在堆顶，沿x、y方向分别为2.13m/s和2.48m/s；当孔隙率为0.4时，抑尘网与料堆之间无涡旋，同孔隙率为0.3时速度矢量分布基本一致，从速度矢量图3（d）中可看到速度变化较孔隙率为0.3时更平缓，这是因为网前后压差进一步减小，速度从堆底至堆顶逐渐增加.

由此可见：低孔隙率（≤0.2）时，空气绕流占主导，抑尘网与迎风面之间区域处于涡旋流；高孔隙率（≥0.3）时，空气渗流作用强，迎风面为贴附流.低孔隙率与高孔隙率的流动特性明显不同，0.2～0.3之间必然有一个流态突变的临界孔隙率，选择此区间的多个工况进行数值模拟，试图寻找临界孔隙率，其流动兼具涡旋流和贴附流的特性.结果显示，当孔隙率为0.25时，如图3（e）所示，迎风面既非全部处于涡旋区，也非完全贴附流.0～5m空气从地面贴附向上流动，5～15m高度处在绕流和渗流的共同作用下形成紧贴迎风面的强度较小的涡旋，涡旋中心高度为10m，15～17m高度时，又呈现出贴附绕流特性.出现这一复合流动的原因在于：在此孔隙率下，网后渗流与上部绕流共同作用，在垂直压差作用下，中下部形成涡旋，料堆底部涡旋作用衰减，在紧贴地面剪切力作用下空气向上扬起，料堆中上部区域流道面积扩大，压差减小，空气保留来流风特点，贴附料堆向上.为充分利用涡旋和贴附流复合特性从而降低料堆起尘，本研究在来流风速为6m/s下临界孔隙率为0.25.

2.2 剪切力

通过计算料堆周围速度矢量分布，分析了迎风面涡旋和贴附流的流动特性，通过计算料堆表面在不同孔隙率下，迎风面和背风面0～17m高度、平顶面0～10m宽度的表面剪切力，可以定量分析迎风面、背风面和平顶面不同位置剪切力的分布特性.

2.2.1 迎风面剪切力 孔隙率分别为0、0.2、0.25、0.3、0.4、0.6和1.0时，剪切力在迎风坡面不同高度的分布如图4所示，可以看出：未设置抑尘网（孔隙率=1.0）时，迎风面剪切力随高度增加而增大，最大剪切力出现在堆顶17m处，即料堆颗粒在来流风作用下沿迎风面向上扬起，堆顶起尘量最大；孔隙率为0.3、0.4、0.6时，剪切力仍旧随高度增加而增大，只是剪切力大小随孔隙率减小而减小.即在孔隙率≥0.3时，抑尘网的作用在于减小风速和表面剪切力，并没有改变湍流结构.孔隙率≥0.3时，表面剪切力均为正，即沿迎风面向上，与前述贴附流的特性吻合.

孔隙率为0时，湍流结构与高孔隙率时迥然不同.剪切力均为负，即沿迎风面向下，且剪切力绝对值随高度增加而增大；这是由于涡旋中心高于料堆，整个迎风面处于涡旋中，料堆上部靠近涡旋中心，速度梯度和表面剪切力较大，这与前述速度矢量分布也相一致.

孔隙率为0.2时，沿堆高方向剪切力不再呈单调关系，而是随高度增加先增大后减小，1～16m高度间，表面剪切力为负，即涡旋到影响区域为1～16m高度；最大剪切力出现在10m高处，为7.75N；整个迎风面剪切力绝对值之和为74.6N，为无网时的16.7％.

图4 不同孔隙率迎风面剪切力分布Fig.4 Shear force profiles on windward of the pile with deferent porosity

由图4知，孔隙率为0.25时，0～5m高度为向上贴附流，剪切力为正，5～14m高度处于涡旋中，剪切力为负，涡旋强度较孔隙率为0.2时小，14～17m高度处又呈现贴附流特性，剪切力为正.可见该孔隙率下，迎风面为贴附流和涡旋复合流动，表面剪切力随高度有正有负，剪切力随高度变化曲线紧贴横轴，绝对值明显比孔隙率为0.2及0.3时小.综合来看，孔隙率为0.25时迎风面剪切力绝对值之和最小，为26.8N，为无网时的6.0％.2.2.2 背风面剪切力 不同孔隙率剪切力在背风坡面不同高度的分布如图5所示，可以看出：孔隙率为0时和其他孔隙率下剪切力随高度变化趋势明显不同.孔隙率为0时，剪切力随高度变化微小，且均为负，沿x负方向，即沿背风面向上，背风面处于回流涡旋中.孔隙率为0.2、0.25和0.3时，剪切力随高度变化曲线几乎重合，在1～13m高度间，剪切力为正，沿迎风面向下，其他高度处贴附背风面向上且剪切力大小变化不大；当高度大于13m后，剪切力随高度急剧增长，最大剪切力出现在17m处，即堆顶.孔隙率为0时，背风面剪切力绝对值之和最小为12.2N.孔隙率在0.4～1.0时，剪切力变化趋势不变，同一高度，剪切力随孔隙率增大略有增大.背风面剪切力绝对值大小大部分在3N以下相对迎风面和平顶面而言较小，且在孔隙率为0.2～0.3时相差很少，所以料堆表面剪切力主要由迎风面和平顶面剪切力组成，背风面影响很小.

图5 不同孔隙率背风面剪切力分布Fig.5 Shear force profiles on leeward of the pile with deferent porosity

2.2.3 平顶面剪切力 不同孔隙率剪切力在平顶面不同宽度的分布如图6所示，可以看出：孔隙率为0时，由于涡旋中心高度高于料堆，平顶面仍处于涡旋流中，表面剪切力沿平顶面与来流速度相反，均为负值；剪切力随宽度增加而减小，平顶面前半程受涡旋影响显著，剪切力变化快，后半程剪切力变化缓.

孔隙率≥0.2时，剪切力均为正，且随孔隙率减小而减小，这是因为平顶面为贴附流，与来流风速方向相同，且随宽度增加而减小，缘由前半部分靠近迎风面，压力损失影响明显，后半部分压力相对平稳.由图6可见，孔隙率在0.2～0.3时，剪切力随宽度变化微小，同一高度剪切力变化不大，最大相差4.1N.孔隙率为0.2时，涡旋中心高度位于迎风面2/3处，未到达平顶面，所以平顶面为贴附流，剪切力为正，随宽度增加变化甚小，平顶面总剪切力为47.94N，是无网时的16.1％；孔隙率为0.25时，总剪切力为67.12N，是无网时的22.5％；孔隙率为0.3时，总剪切力为85.73N.即孔隙率在0.2～0.3时，平顶面表面剪切力相差不大.由于料堆为长条棱台，两侧面面积相对迎风面、平顶面和背风面面积很小，所以考虑迎风面、背风面和平顶面剪切力之和.通过计算比较，孔隙率为0.25时，总剪切力最小为93.95N.

图6 不同孔隙率平顶面剪切力Fig.6 Shear force profiles on top flat of the pile with deferent porosity

3 结论

3.1 料堆周围速度矢量分布显示：低孔隙率（≤0.2）和高孔隙率（≥0.3）时流动特性迥异.低孔隙率时，迎风面处于涡旋中：孔隙率为0时，涡旋中心高于堆高3m，迎风面和平顶面上均为涡旋流；孔隙率为0.2时，涡旋中心位于迎风面2/3高度（12m）处，平顶面为贴附流.高孔隙率时，堆前无涡旋，迎风面和平顶面均为贴附流.孔隙率为0.25时，迎风面为涡旋和贴附流的复合流动，涡旋中心位于迎风面3/5高度处，5m以下和14m以上为贴附流.

3.2 料堆各表面剪切力分布显示：孔隙率≥0.3时，迎风面同一高度，剪切力随孔隙率增大而增大；同一孔隙率，剪切力随高度增加而增大.比较孔隙率为0，0.2，0.25，0.3，0.4，0.6时的迎风面剪切力，孔隙率为0.25时最小，仅为无网时6.0％.平顶面随宽度增加表面剪切力减小，孔隙率为0.2时，平顶面剪切力之和最小，为无网时的16.1％.

3.3 孔隙率为0.25时，堆前流场特性为涡旋流向贴附流过渡的转折点，且料堆各表面剪切力之和最小.故来流风速为6m/s时，在该几何模型、研究工况下，临界孔隙率为0.25.

［1］郝吉明，程 真，王书肖.我国大气环境污染现状及防治措施研究［J］.环境保护，2012，9：17-20.

［2］李 林，郝吉明，胡京南.北京市能源利用对空气质量的影响分析和预测［J］.中国环境科学，2005，25（6）：746-750.

［3］贺克斌，余学春，陆永祺，等.城市大气污染物来源特征［J］.城市环境与城市生态，2003，16（6）：269-271.

［4］丛晓春，陈志龙，朱晓敏.风蚀性粒子降尘浓度的动力学预报模式［J］.中国环境科学，2012，32（7）：1171-1176.

［5］李建隆，董纪鹏，陈光辉，等.防风抑尘网开孔形式对流场的影响［J］.环境工程学报，2009，3（9）：1725-1728.

［6］陈 凯，朱凤荣，钮珍南.防风网作用效果的风洞实验评估［J］.北京大学学报（自然科学版），2006，42（5）：636-640.

［7］宣 捷，俞学曾.风障减少尘埃飞起的风洞模拟研究［J］.环境科学研究，1997，10（2）：14-18.

［8］Park C W，Lee S J.Experimental study on surface pressure and flow structure around a triangular prism located behind a porous fence［J］.J.Wind Eng.Ind.Aerod.，2003，91（1/2）：165-184.

［9］Wang H，Takle E S.Numerical simulations of shelterbelt effects on wind direction［J］.J.Appl.Meteorol.，1995，34（10）：2206-2219.

［10］丛晓春，詹水芬.防风网多孔动力效应的应用研究［J］.中国矿业大学学报，2009，38（2）：193-196.

［11］赵海珍，梁学功，马爱进，等.防风网防尘技术及其在我国大型煤炭港口的应用与发展对策［J］.环境科学研究，2007，20（2）：68-71.

［12］段振亚，石文梅，郑文娟，等.防尘网抑尘机理研究及工程应用进展［J］.石油化工设备，2010，39（3）：40-44.

［13］金向阳.防（挡）风抑尘网在露天储煤场的应用［J］.煤炭工程，2008（6）：37-38.

［14］李晋旭，李克民，杨 明，等.防风抑尘墙在黑岱沟露天矿中的应用研究［J］.金属矿山，2010（1）：159-162.

［15］唐继臣，孙曙光.防风抑尘网在大型原料场的应用［J］.烧结球团，2010，35（5）：31-34.

［16］张光玉，陈 立，王奇志，等.秦皇岛港煤堆场防风网风洞试验研究［J］.交通环保，2003，24（1）：4-6.

［17］Dong Z B，Luo W Y，Qian G Q，et al.A wind tunnel simulation of the mean velocity fields behind upright porous fences［J］.Agr.Forest Meteorol.，2007，146（1）：82-93.

［18］胡晋明，沈恒根.风障屏蔽作用的数值模拟［J］.污染防治技术，2003，16（4）：1-3.

［19］姜 鑫，宋 力，苏 猛，等.防风抑尘网数值模拟及设置方案分析［J］.内蒙古师范大学学报（自然科学汉文版），2013，42（5）：544-548.

［20］Lee S J，Park C W.Surface-pressure variations on a triangular prism by porous fences in a simulated atmospheric boundary layer［J］.J.Wind Eng.Ind.Aerod.，1998，73（1）：45-58.

［21］Dong Z B，Luo W Y，Qian G Q，et al.A wind tunnel simulation of the turbulence fields behind upright porous wind fences［J］.J.Arid.Environ.，2010，74（2）：193-207.

［22］陈廷国，马思明.不同截面形式防风网流场的3D数值仿真［J］.计算机仿真，2014，31（1）：258-263.

［23］Perera M D A E S.Shelter behind two-dimensional solid and porous fences［J］.J.Wind Eng.Ind.Aerod.，1981，8（1/2）：93-104.

［24］郭 辉.防风网遮蔽效果研究［D］大连：大连理工大学，2008.

［25］周伟朵.用数值模拟方法研究挡风抑尘网高度对抑尘效率影响［J］.电力科技与环保，2012，28（1）：46-47.

［26］Park C W，Lee S J.Verification of the shelter effect of a windbreak on coal piles in the POSCO open storage yards at the Kwang-Yang works［J］.Atmos.Environ.，2002，36（13）：2171-2185.

［27］宋翀芳，彭 林，白慧玲，等.露天堆场防风抑尘网后湍流结构及抑尘效率的数值模拟［J］.中国环境科学，2014，34（7）：1690-1695.

［28］Li G，Ronaldo G M.Numerical simulation of airflow and particle collection by vegetative barriers［J］.Eng.Appl.Comp.Fluid，2012，6（1）：110-122.

［29］Chen G H，Wang W W，Sun C F，et al.3D numerical simulation of wind flow behind a new porous fence［J］.Powder Technol.，2012，230：118-126.

［30］陶文铨.数值传热学［M］.西安：西安交通大学出版社，2001：32.

Numerical simulation for critical porosity of porous fences used to shelter open storage piles.

PAN Wu-xuan，HE Hong-zhan，SONG Chong-fang*，LI Lin-ping（College of Environmental Science and Engineering，Taiyuan University of Technology，Taiyuan 030024，China）.China Environmental Science，2015，35（6）：1638～1644

The flow field around the open storage piles behind a porous fence was numerically simulated by the software Fluent 6.3.The critical porosity was determined based on the datas obtained from the simulation of flow fields for the porous fence with different porosities and the shear stress distribution above the windward side，flat top surface and leeward side of the pile.As a result，the flow structures at high porosities（0.3，0.4，0.6）were similar to that of the unfenced condition.The air velocity was higher from the ground to the top of the windward surface and the velocities were all positive.At low porosities（0，0.2），however，the flow structure differed strongly with that of the high porosities.Large-scale vortices were formed between the fence and the storage pile，thus the velocity of the air above the windward was negative.For the fence with ε=0，the height of vortex center was higher 3m than the height of the pile and for the fence with ε=0.2，the vortex center located at the 2/3 height of the windward surface.In addition，the shear stress distribution was similar to that of the unfenced condition when the fence porosity was exceeded 0.3and the shear stress increased with the increasing porosities.For the ε=0.2 porous fence，however，as the height of the pile increased，the shear stress increased，then decreased and the maximum shear stress was at the 3/5height of the pile.When the fence porosity was 0.25，both vortex and attached flow were formed above the windward surface and the shear stress of the pile was to be the least.Thus the porosity ε=0.25 below or above which turbulence structures and the shear stress distributions differed strongly was determined to be the critical porosity in this study.The method studying aerodynamics microscopic characteristics around the piles provides a new idea for the porous fence investigation.

porous fence；numerical simulation；critical porosity；shear stress

X513

A

1000-6923（2015）06-1638-07

潘武轩（1990-），女，山西运城人，太原理工大学硕士研究生，主要从事空气污染控制和暖通空调系统的数值模拟研究.

2014-10-31

国家自然科学基金资助项目（51108295）

* 责任作者，副教授，scfcindy@163.com