植物被广泛应用在防风固沙工程中。植物能遮蔽地表颗粒防止其发生起动，可以通过吸收风动能从而减小风的剪切应力，能够拦截和捕获沙尘颗粒^［1-2］。

衡量沙生植物防护效益的参数有多种，包括植物种类、植被盖度、植物尺寸、植被群落的配置模式以及植物柔韧性等，如何准确地建立上述参数与植物的防护效能的关系是防沙工程中的关键问题^［3-4］。有关植物防风固沙效益的研究通常把单株植物的孔隙性质或植物群落的配置模式作为研究的重点，而关于其他参数如何影响流场与防护效能目前还没有深入的研究，特别是植物柔韧性如何影响冠层流动与周围微观气候。已有很多野外实验与风洞模拟证明，植物的柔韧性是对于流场变化非常敏感的重要参数，其对于周围流动特征和其整体防护效益有显著的影响^［3］，所以很有必要开展关于植物柔韧性的流体动力学研究。

目前关于水生柔性植物研究较为丰富，主要研究恒定水流与波浪作用下柔性植物的阻力特性^［5-6］，进而研究影响泥沙悬浮与污染物的运输机制，研究方法以水洞模拟和野外实验为主。对柔性植物单元，一般选取真实植物或是将其简化为柔性复合材料制成的植物模型。黄本胜等^［7］关于河滩种树对行洪影响的野外观测表明植物柔韧性对于固滩固岸与防浪护堤的益处。吴福生^［8］和Ghisalberti等^［9］通过水槽实验对淹没的柔性植物进行研究，通过激光雷达测速仪对柔性植物周围流场的信息进行测量，分析不同冠层流动区域的湍流结构以及剖面风廓线信息。

旱生柔性植物相关研究主要关注干旱地区植物沙障柔韧性对地表颗粒物质运输的影响与侵蚀作用^{［2，10-11］}，得益于计算流体力学（Computational Fluid Dynamics，CFD）技术的发展，加上数值模拟通常能够较为全面地反映湍流影响下的柔性植物周围的流场结构，利用CFD研究柔性植物也逐渐成为主流的方法。Liu等^［10］进行了一系列数值模拟实验，确定不同风速下植物沙障的有效高度、叶面积指数与植被盖度。Wu等^［11］通过一系列三维数值模拟，研究气流通过柔性单株防风林的流动情况，比较了不同柔韧性植物的有效防护距离，发现柔性植物比刚性植物具有更好的防护效果。然而，由于描述植物柔性特征的数值模型建模过程较为复杂，需要考虑植物秆茎受到的流体力以及柔性植物反馈给流场的固体力。另外，大气边界层内雷诺数较大，植物周围的流动更为复杂，进行中小尺度的数值模拟也需要耗费巨大的计算资源。He等^［12］使用浸没边界法描述柔性植物在风洞尺度内的流场结构与能量运输特征。但这种对于柔性植物的描述方法由于高昂的计算成本，不能有效应用到真实大气边界层或气候模式的数值计算中。若忽略植物的柔性特征，仅为解决工程实践应用的问题，在植物区域内引入多孔介质模型是目前主流的数值方法，该方法通过定义孔隙介质区域和其中的阻力系数来等效流体经过植物时的扰流作用。刘彦东等^［13］用RANS湍流模型对沉水植物进行三维数值模拟。尽管许多研究表明多孔介质模型能够准确描述植物的孔隙性质^［14-15］，但该方法仍然不能描述植物的柔韧性对于流体的影响。

本文首先以多孔介质模型为基础，将传统模型中的阻力系数修改为时间相关动态阻力系数模型（Varied Porosity Model），通过源项附加法将阻力系数以源项的形式引入到大涡模拟动量方程中进行求解，从而实现柔性植物随着流体不断发生摆动的效果。然后，通过对比流固耦合（Fluid Structure Interaction，FS1）和VPM数值方法描述的柔性植物，证明将柔性植物简化为VPM单元的合理性。最后，将植物柔韧性量化为VPM方程中的摆动频率f与摆动幅度A，分析不同柔韧度的植物在不同冠层流特征区域的流场特性与湍流结构，此外还考虑到工程实践中植物防护效能的评估需求，定量分析植物柔韧性对挡风系数与风速廓线的影响。

本文采用大涡模拟对流场进行求解，其核心思想是将流体运动过程中的大尺度涡旋运动，用空间滤波函数处理后的纳维-斯托克斯方程（N-S）求解，而将尺度小于网格耗散尺度的涡旋运动模化为亚格子模型求解^［16］，对应的不可压缩黏性流体的控制方程为：

式中:t为时间；xi为笛卡尔坐标系下的3个方向；u¯i为过滤后的速度分量；p¯为压力；ρ为流体的密度；ν为流体的动量黏度；τij为亚格子应力张量，表示小尺度涡对大尺度涡的影响；Si为体积力源项。本文采用经典Smagorinsky-Lily亚格子模型对亚格子应力张量进行模化，对应的亚格子混合长度ls为：

式中：k为冯卡门系数；y为网格中心距离壁面的高度；Cs为Smagorinsky常数；Δ为网格的体积。

对于植物区域，本文将单个植物简化为方柱，并以Darcy-Forchhermer（D-F）模型为基础描述流体接触到植物后受到的阻力作用。由于一般真实环境条件下的沙生植物周围的流动雷诺数高，湍流强度大，本文忽略D-F模型中的黏性阻力项，将孔隙植物对流体的衰减作用体现为动量方程（式1）中的附加源项Si：

式中：ρ为流体的密度；νi为流体表观速度在不同方向上的分量（i=x，y，z）；ν为合速度值；F为阻力系数。在传统的多孔介质模型中一般将F定义为孔隙率或是叶面积指数等与植物构型参数相关的常数，一般通过试验获得^［17］。F的大小能够影响空隙介质内部的流速，当F足够大时，内部流速大幅减小，Chang等^［18］在对孔隙圆柱的大涡模拟研究中发现，当体积分数SVF>0.1时，全局时均阻力系数CDG>1.17，相当于实心体圆柱。若将F定义为关于时间的函数，植物内部受F控制的部分会在孔隙介质与固体之间周期性变化，这一特点与柔性植物在风中摇摆的性质相似。因此，本文为了模拟摇摆效应，将F定义为动态的、关于计算时间t周期性波动的正弦函数：

式中：A表示F的振幅，一般与植物的孔隙率相关；f为F波动的频率，一般与植物的柔韧性相关；L为F偏移量，通过该参数可以控制柔性植物在摇摆过程中的倒伏程度。下文将引入式（4）Si的数值模型称为动态阻力系数模型（VPM）。

如图1A所示，单个植株的尺寸是边长为d的正立方体。整体计算域尺度为10d×5d×5d，植物位于中心位置（5d，0，2.5d）。计算域入口处设置为完全发展的湍流入口，湍流度为6.5%，出口处为压力出口，上边界为滑移边界，下边界为非滑移边界，左右边界为循环边界。使用大涡模拟对边界层流动进行求解，不使用壁面函数，该工况的雷诺数Re=8800，第一层网格高度设置小于0.002以保证y⁺<1。设置计算时间步长Δt=0.05以保证每个时间步的库朗数均小于0.5。

统计2 s内的5个计算结果，分别绘制瞬态速度、瞬态压力以及涡结构云图（图2）。植物模型立方体内的速度会随着时间变化产生周期性波动，在100 s时气流能够进入冠层内部，但随着时间发展，内部F逐渐升高，到102 s时可以看到冠层内部已经没有气流进入。冠层前端边界周期性地受到压力，在图中由红色显示，随着时间变化其受力由强变弱，在102 s时又受到较高的压力。该受力特性符合植物随风摇摆机制：植物前端受力向后发生倒伏，然后由于植物自身弹性又逐渐恢复为原来的形态。最后本文使用Q准则对植物周围的涡结构进行捕捉，从100 s时刻可以看到当植物前端受压较大，流体未能进入冠层内部时，植被前端顶部积聚了一个较大尺度的分离涡，当流体逐渐进入冠层内部，这个分离涡破碎并发展为较小尺度的尾流涡，当流体不再能够进入冠层内部时，如102 s时刻所示，这个叫大尺度的分离涡会再次生成（图2C）。

根据柔性植物模型的实验研究，监测植物内部接近剪切层位置的流向速度，流向速度随着时间不断变化并且呈现出周期性波动效应（T≈10 s）。以实验数据为依据，设置VPM模型中f=0.1进行数值模拟，并设置一个不引入VPM的算例进行对照。在VPM模型中，将单棵植物简化为一个矩形的多孔介质区域，本文借鉴先前相关研究将植物秆茎建模为细杆^［19-20］，以直径为0.1d的方柱作为植物秆茎结构，排布3×3的方柱阵列描述单棵植物几何，入口来流条件和计算域尺寸与VPM方案保持一致（图1B）。方柱阵列中每个方柱单元均采用FSI的数值方法模拟柔性植物在流体作用下的摇动效应，设置固体计算域的密度为1 000 kg·m^-3，弹性模量为1×10^-6 N·m^-2，泊松比为0.47。

在相同来流的条件下，对比实验数据与VPM、FS1模型的数值计算结果（图3），VPM模拟了冠层内部流速随时间的周期性波动变化规律，且随着时间的发展，速度有逐渐减小的趋势，在时间370、380、390 s位置处的波峰速度逐渐降低。相较于VPM，FSI能获取植物内部更加细微速度变化，但考虑到FSI模型的稳定性以及对计算资源的需求，在关于中大尺度的植被下垫面数值模拟中，使用VPM对柔性植物进行描述会极大提升计算效率，有利于实际工程应用。对比不引入VPM的刚性植物，尽管流速在湍流作用下也出现了波动，但并未出现明显的周期性波动，与实验结果差距较大。由于阻力系数在冠层内部保持不变，植物被假设为刚体，所以其对于流场的扰动效应会相对较弱，导致其剪切层位置处的平均速度大于柔性植物。

此外，VPM也能反映柔性植物形态发生变化后对流场的影响。在VPM算例中，提取空隙介质内部速度为0.001 m·s^-1的等值面，如图4A所示，灰色区域表示植被内部速度小于0.001 m·s^-1，这时VPM模型中的F处于正弦函数的最高值（F=4），对应图4B中FSI模拟直立状态下的柔性植物，植物的迎风面与来流垂直。随着阻力系数的减小，冠层内部前端的速度逐渐增大，流体能够更多地进入到空隙介质内部，在100.5 s和101 s时刻，流体的迎风面倾斜并后移，这与FSI算例中柔韧方柱产生倒伏的机制类似。

为了避免网格数量影响数值计算结果的准确性，对同一工况下不同数量的网格进行计算。设置3套网格分别为网格1（80×50×60）、网格2（100×70×80）和网格3（140×100×120）。对比不同流向距离位置处的时间平均风速廓线，选择的3个流向位置分别为x=1.25d、x=4d和x=8.75d。结果如图5所示，网格1由于边界层网格较为稀疏，所以对于流场内的湍流计算不准确，导致流向平均流速3个位置处均高估。网格2和网格3的计算结果在位置x=1.25d和x=4d相近，可以认为网格2足以准确解析柔性植物影响下的流场特征。在x=8.75d位置处，最为精细的网格3在z/d的位置出现了轻微的振荡，但在z/d>1时又与网格2的计算结果几乎重合。综上所述，本文可以选择网格2进行下一步计算，在保证计算精度的基础上可以极大提高计算效率。

柔性植物的扰流特性会随着其柔韧性而产生变化，在本文所提到的VPM柔性植物模型中，植物的柔韧特性由孔隙区域的阻力系数F控制。在模型验证部分的VPM相关参数为：振幅A=0.2，频率f=0.1，位移L=0.2，表示孔隙介质内部的F会以周期T=10s在最大F=0.4和F=0之间波动，本文以该实验中的植物为基础，微调VPM模型中的参数并分析柔韧性所带来的变化。假设植物的柔韧性越大，其越容易受到湍流影响而产生摇摆效应，在数学方程上体现为孔隙介质内部的F波动周期越小。为了研究植物摇摆效应对于流场特性的影响，本文设置5种摆动频率，对应不同柔韧性的植物，具体参数设置如表1所列。

对比前3种工况能够分析植物摆动频率的改变对流场的影响；对比工况4与工况1能够分析出同样柔韧性的植物在F增大产生的效应；对比工况5与工况4，可以得出植物区域F在0.75~1波动与F在0~1波动有何异同，前者对应摆动幅度较小的植物，F仅能够在0.75和1之间变化；而后者对应摆动幅度较大，F较小时（F=0），这时植物为完全倒伏状态，对流场没有阻塞效应，随着植物柔性特征作用逐渐恢复为直立的状态，F逐渐恢复为1（图6）。

在植物防风固沙工程中，植物的挡风效应主要由挡风系数D衡量^［21-22］：

式中：Ufront为前端接触区速度；Uback为近尾流区的速度。此外，植物能显著影响近地表的动量输运并有效减少地面沙粒的起动^［4］，本文将对柔性植物影响下的风廓线与地表剪切应力进行定量分析。

随着植物的摆动频率不断增大，瞬态挡风系数D逐渐减小（工况1∶0.81；工况2∶0.78；工况3∶0.72），说明柔韧性较大的植物对于气流的阻塞效应更好。但是工况1在350 s位置处出现了一个最小挡风系数（D=0.08），意味着气流在这一时刻几乎没有障碍地穿过了植物冠层，同时说明近尾流区的速度几乎增大为来流风速，偶尔出现的强阵风会使得植物背侧受保护的区域遭到破坏。

相同柔韧性但F不同的植物挡风系数如图7所示，F的减小会使植物的总体挡风效应减弱，这与实际中孔隙度大、F小的植物挡风能力弱的结论相符。而工况5和工况4对比了F在0.75~1和0~1波动对于其挡风能力的影响，结果显示摆动幅度更大的植物的挡风能力更强。

由于各个算例在100 s时出入口流量已到动态平衡状态，选择用400 s的计算结果进行时间平均。对比瞬态速度与压力可以看出，由于植物摆动频率不同，同一时刻下植物的摆动状态也不同（图8）。例如在300 s时，工况1植物前侧受力，而工况3前侧几乎不受力，证明该时刻植物的F正处于低值。

气流与植物的相互作用在其周围不同位置表现出截然不同的流动特征。Belcher等^［23］在关于植被区的研究中将植物冠层流动大致分为前端接触区、冠层内部与尾流区域。本章选取典型冠层流动区域的速度与压力进行研究，定量分析植物柔韧性如何影响其湍流特性（图9）。

由图10A与图10C可知，植物前端区域的平均和脉动速度与植物摆动频率的关联性较弱，这是由于流体还未接触到植物时，动态阻力系数模型（VPM）对流场的影响较弱。但从速度功率谱密度曲线能够看出，柔韧性越高的植物对应的总体能量越高，且摆动频率最大对应的工况3在2 Hz位置处出现了波峰（图10B），说明在这个工况下，速度在较高的频率下波动更为显著。相较于植物摆动频率对于接触区流速的影响，对比工况1、工况4和工况5（图10D），能够看出流速对于F的变化更敏感。相同柔韧性的植物对应的F越小，接触区流速更大，且F在0.75~1相较于F在0~1对于流速的削弱效应更强。

从图11C可以看出，随着植物摆动频率的增加，其内部平均速度会随之降低，说明植物柔韧性越高，其摇摆效应越强，植物内部的流动会被更多地抑制。脉动效应的强弱可以通过对比其脉动速度均方根来实现，由图11A可得，摆动频率越低，其脉动速度越大，对应的湍动能也越高，说明在冠层内部流动会更加稳定，这有益于植物内部的滞尘效应与防起沙机制。

对比工况1与工况4，在植物柔韧性相同、F减小时，冠层内部平均流速增大，由于F的大小与孔隙率密切相关，孔隙率越大导致流体更多地穿过冠层。对比工况4与工况5，发现冠层内部的均匀性会对其内部流速变化产生很大影响，F在0.75~1对应的植物相较于F在0~1能够消耗更多风场动能，冠层内部的平均风速更低，且脉动速度也更低，可以说明垂向不均匀的植物对风场具有更强的削弱作用。

尾流区的流动特征会显著影响植物的影响范围，对应到工程应用中，尾流区流速的大小会决定植物保护范围与保护效率。对比尾流区与冠层内部的流速，其随时间的变化不再呈现较强的周期性波动规律，在靠近植物的位置，流速在某些时刻出现了负值，这是由于产生回流涡而导致的。从图12C可知，植物的柔韧性对于近尾流区的平均流速影响较小，但仍然可以看出植物摆动频率越大平均流速更小。由图12A可知，摆动频率高的工况3整体能量越高，对应的脉动速度与湍动能会更小。在近尾流区，由于气流刚经过植物冠层，其流动特征与规律与冠层内部相同。

在远尾流区域，流向水平速度已经恢复为来流速度的50%左右，但植物摆动频率对平均流速的影响较小。从图13A可知，摆动频率越大的植物在远尾流区的速度更小。对比功率谱密度曲线，所有工况均没有出现峰值，但整体能量相较于近尾流区增大较为明显，说明植物的摆动使得速度波动效应在时间上更均匀且数值更大。脉动速度相较于近尾流区与冠层内部增大了2倍，这是因为远离植被的保护，远尾流区域的脉动效应更强，湍动能增大了2倍。由于植物的摆动消耗了一部分能量，所以柔韧性高的植物对应的尾流区湍动能更低。

相同柔韧性的植物在F更小时，近尾流区的速度更大，这与冠层内部的规律相同，因为F的大小影响了气流通过植物的速度大小（图14）。但在尾流区较远的位置，F的影响变小。这说明在远尾流区，无论是植物的柔韧性或孔隙率，对于流速的影响都很微弱。同样地，对比工况4与工况5发现，植物的摆动幅度对于远尾流区域的流速影响也很微弱。

监测冠层5个流动区域位置处的风廓线（L1~ L5），流向位置分别为x=1d，x=4d，x=5d，x=6d，x=9d。植物的摆动频率对于风廓线的影响较小，变化比较明显的位置出现在L3和L4，摆动频率越小的植物对应的风廓线在1d高度以内平均风速更小，说明该位置的垂向保护效应更强。而当气流发展到L5，植物的摆动频率几乎对尾流区没有影响（图15）。相较于直接用多孔介质模型描述植物模型，对比L3与L4，若仅使用多孔介质模型描述植物，在设置相同F的情况下，刚性植物内部和近尾流区域的数值结果会高估植物的挡风效应。从L3位置处可知，引入VPM的3种工况冠层内部的相对流速为0.04，而刚性孔隙植物对应的相对速度为0.43，高估的幅度为10.5倍。在近尾流区，植物的阻塞效应导致冠层高度h以下均出现了速度负值，但刚性孔隙植物并不能描述该机制。

在植物防风固沙工程中，地表剪切应力是一个重要的物理参数，其能够在很多方面反映出气流对地表的作用力以及地表抵抗该作用力的能力^［3-4］。首先，地表剪切应力的大小可以直接影响风蚀过程，当风速到达一定阈值时，风开始对地面沙尘颗粒产生剪切作用，使得颗粒发生起动从而引起风沙运动。此外，植物通过其自身结构改变风场的分布，降低近地面的风速，从而减少地表受到的剪切应力，这有助于减少风蚀，保护地表不受侵蚀^［24］。因此本文将对比描述VPM下不同柔韧性植物对地表剪切应力大小与分布的影响。

在102 s时，内部F较大，植物处于直立状态，其左右两侧地表剪切应力增加，前后剪切应力减小（图16），这与亢力强等^［24］的研究结果一致。在103 s时，植物发生倒伏，植物内部前缘的剪切应力增大，两侧剪切应力减小，在植物周围4个角有明显的峰值区域，随着植物由倒伏恢复为直立状态，其两侧的剪切应力又逐渐增大。提取尾流区中轴线上的瞬态地表剪切应力，如图17B所示，地表剪切应力随着F的变化也产生了周期性波动。在102 s时F位于峰值，对应的地表剪切应力最大，相应地，在106 s时，F处于低值，对应的剪切应力最小。

图17A反映出柔韧度不同的植物对地表剪切应力的影响，最为明显的是刚性植物对植物尾流区地表剪切应力的削弱作用较弱，在前端接触区最大相差2.2倍，在尾流区最大相差1.64倍，尽管刚性植物内部的抗风蚀能力较强，但总体来说柔性植物要优于刚性孔隙植物。对比工况1~3可得植物的柔韧性越大，摆动频率越高，尾流区地表剪切应力越小，抗风蚀能力越强，相对剪切应力平均值分别为0.0027（工况1）、0.0025（工况2）、0.0023（工况3）。

本文通过定义柔性植物的摆动频率f与摆幅A来量化植物柔韧度，结合LES计算获得植物周围的流场特性与湍流结构，将气流与植物相互作用的区域划分为前端、内部与尾流区。结果显示，相较于刚性植物，在柔性植物每个区域的平均速度更小，但脉动速度和湍动能更大。这是因为VPM引入了阻力系数F关于时间的脉动特性，从而使得植物周围的湍流效应增强，且f越小，湍动能增幅越大。植物的A是由最大阻力系数控制，因此尾流区平均速度随着A的降低而增大，防风效果更差。从防风固沙工程应用的角度出发，本文分别对挡风系数、风速廓线和地表剪切应力进行了定量分析，结果发现植物的柔韧性会显著影响各类防风效能指标的评估：柔性植物对于气流的削弱作用更强，且对应的地表剪切应力更小。由于目前以植物为研究对象的大部分数值模拟研究仅考虑植物孔隙性质，也就是假设植物为刚性，这导致植物防护效率与其配置参数（高度、排数、间距等）之间的关系预估不准确，工程效益降低。应用流-固耦合（FSI）方法解析柔性植物会使得数值模拟的求解难度与要求的计算资源剧增，特别是需要大尺度求解植被下垫面时^［5］。本文提出的数值方法能够很好地平衡完全解析植物柔性和完全不考虑柔性，从而为提升工程效益和数值模拟效率提供有效参考。在今后的研究中，将植物的柔韧度参数化为材料本身的特征参数将进一步方便工程应用，研究弹性模量、密度与泊松比等于F参数化方案是有价值的课题。此外，本文对于植物的柔性描述为外部模型干预，F没有根据局部流动特征而动态变化，在使用该模型时需要对研究对象的摆动参数与流场特征进行严格标定，因此后续关于植物摆动随着流场动态变化的数值模型研究也很有意义的。

本研究在多孔介质模型基础上，改进了阻力系数模型，这是一种较为高效的描述植物柔韧性的数值方法。本文应用该方法模拟了随气流发生周期性摆动的柔性植物，监测其特征区域的流场信息，分析了植物摆动参数对风场特性的影响。

植物摆动频率对风廓线的影响在冠层内部和近尾流区较为明显。摆动频率越小的植物在冠层内部的垂向保护效应越强，而在近尾流区能越好地抑制速度负值的出现。

植物柔韧性越大，摆动频率越高，挡风系数越低，表明柔韧性大的植物对气流的阻塞效应更好。但高摆动频率可能导致短暂的强阵风现象，使得其防护效果完全失效。

柔性植物对地表剪切应力的削弱作用明显优于刚性植物，特别是在前端接触区和尾流区。植物的柔韧性越大，摆动频率越高，尾流区的地表剪切应力越小，抗风蚀能力越强。