风是森林生态系统中的主要干扰因子，我国每年因风灾所造成森林破坏的损失大约占森林损失总量的53%^[1]，尤其是沿海地区林业受损严重^[2]，森林管理治理成本也不断增加^[3]。在我国华北平原黄河入海口地区，土壤盐碱化严重，刺槐（Robinia pseudoacacia）由于其耐盐碱特性成为该地区防护林常用树种，但由于其根系较浅，在风的侵袭下极易风倒。为削弱风速，降低风灾对刺槐林木的破坏程度，提升刺槐林的防风效应，通过优化林带结构，对易发生风倒危险的林木进行及时的加固是最优选择^[4]。

目前的研究表明，防护林的防风效应与林木的高度、行间距、行数、树冠的形状、相互之间的排列方式及林间孔隙度等因素有关^[5-6]。对于防护林防风效应的研究手段主要包含现场风速实测试验、风洞试验和数值模拟试验^[7]。通过在防护林研究场地内搭建风速传感仪可很好地监测到林内风速变化规律的实测数据^[8]，但由于实际场地林木排列固定，不利于对影响防风效应的因素进行分析，单一的现场试验难以得到规律性的研究成果。利用风洞试验将简化的林木模型通过调整树高、树冠的结构可较好地分析林带疏透度对防风效应的影响^[9-10]，但由于风洞试验中通常用到模拟树木的模型与实际树木在尺寸、形态及力学性能上有较大差异，对研究结果的准确性具有一定的影响^[11]。通过有限元法构建林木模型，在形状及材性上更贴近于林木真实情况，可以更好地研究林木防风效应^[12-14]。以线性滤波法中的AR法模拟实际状态下的脉动风速时程作为风荷载^[15]，对林木进行动力学分析更大限度地还原了林木在实际大风作用下的受力情况^[16]。随着计算流体动力学（CFD）技术的愈发成熟，CFD技术已被证实可以很好地进行林木与风场间流固耦合的计算^[17-19]。有学者对黄河入海口流域人工种植的刺槐进行力学性能测定，为数值模拟还原刺槐真实材性提供支持^[20]。对单棵刺槐进行建模并通过与风场的耦合计算分析刺槐的破坏机理，为林场防护提供科学依据^[21]。将刺槐林带以网状结构进行简化并与流场进行耦合计算，探究不同林带宽度的防风效应^[22]。

目前利用数值模拟技术针对刺槐林防风效应的研究已有了一定的进展，但现有的模拟研究对刺槐林带的还原度较低，关于不同林带结构对刺槐林防风效应的影响尚不明确，亦未给出合理的刺槐林带结构配置，亟待进一步深入研究。本研究以我国山东省东营市沿海地区种植的人工刺槐林为研究对象，将现场实测数据与数值模拟技术相结合，通过现场实测风速数据结合线性滤波法生成研究地总风速时程，以实际刺槐尺寸为研究对象进行仿真建模，利用有限元软件Ansys Workbench实现多棵刺槐林木模型与风场的双向流固耦合计算，阐明不同林带结构中风速的变化规律以及刺槐受风场影响枝干位移规律。研究能够为合理配置刺槐林林带结构，提升防护林防风效应和区域防护林合理配置提供理论和科学依据。

1.   现场实测

1.1   研究区概况

研究区位于山东省东营市军马五分场中的人工刺槐防护林内（37°47′N，118°41′E），属温带大陆性季风气候，东南风为主导，累年平均风速为4.2 m/s。土壤为滨海岩土，沙地棕壤，适宜刺槐的生长^[23-24]。受海风的影响，场地部分刺槐已发生风倒或枝干折断破坏，现场刺槐风害情况如图1所示。

图  1  现场刺槐风害情况

Figure  1.  Wind damage of Robinia pseudoacacia on site

下载: 全尺寸图片 幻灯片

1.2   现场试验布置

为采集林前10 m高度处风速情况，计算刺槐林对风速的削弱情况，实测试验在林前入口及风场30 m处设置2组14 m高测风支架，可对2、6、10和14 m高度处的风速时程进行监测，风场10、20、40、50和60 m5个位置处各设置一组2 m高测风支架，对经过林场的风速变化情况进行监测。利用WS-5921数据采集仪与0 ~ 5 V电压输出型风速传感器（量程为0 ~ 30 m/s）可同时监测全部位置处的实时风速。为采集刺槐枝干的风振位移情况，在林中20 m处的刺槐枝干处架设1个位移传感器，测量位移时可同步监测风速数据，现场实测试验简易图如图2所示。

图  2  现场实测示意图

Figure  2.  Schematic diagram of field measurement

下载: 全尺寸图片 幻灯片

1.3   风　速

为充分证明当地实测数据的有效性，选取2017年5月11日24小时整点时刻林前15 m处风速采集数据，并与东营市河口区气象站在当日记录的整点时刻的平均风速进行比较，结果如图3所示。由图可知：在0—16时，实测风速时程与气象数据接近一致，17、18、22和24时的风速偏差稍大，其中24时风速测量差值最大，相差1.51 m/s。偏差的原因是本研究试验场地与气象站所处的地理位置的差异以及周围环境的不同皆会导致风速记录值上的偏差，但整体来看，试验场地实测风速时程情况与当地气象站记录的数据无论在变化趋势还是相应时刻的风速值皆具有较高的一致性，可以反应当地典型的风速时程情况。

图  3  实测数据与气象站数据对比图

Figure  3.  Comparison of measured data and meteorological station data

下载: 全尺寸图片 幻灯片

2.   数值模拟方法

2.1   流固耦合的守恒方程

双向流固耦合的求解模式通常是利用流体求解器和固体求解器同时计算流体与固体的相互作用，目前利用软件进行求解的主要计算模式为分离求解，即两种求解器求解各自的物理量，流体求解器将所求的应力数据作为固体求解器的载荷条件，而固体求解器产生的位移重现传递回流体求解器中作为流体域新的边界条件，二者通过数据交换实现双向流固耦合的过程以此解决瞬态分析问题，流固耦合交界面必须满足守恒方程：

式中：${\tau }_{\mathrm{f}}$为流体的应力；${\mathit{n}}_{\mathrm{f}}$为流体的单位方向向量；${d}_{\mathrm{f}}$为流体的位移；${q}_{\mathrm{f}}$为流体热流量；${T}_{\mathrm{f}}$为流体温度，同理等式另一侧为固体对应的各物理量。

2.2   刺槐模型的参数

根据场地对刺槐尺寸的实测确定刺槐模型的尺寸，模型树高10 m，胸径为0.18 m，冠幅为2 m。刺槐是典型的正交各向异性材料，在计算时需要设置9个弹性常数，其中包含纵向、径向、弦向3个方向的弹性模量${E}_{\mathrm{L}}$、${E}_{\mathrm{T}}$和${E}_{\mathrm{R}}$，剪切弹性模量${G}_{\mathrm{L}\mathrm{R}}$、${G}_{\mathrm{T}\mathrm{L}}$和${G}_{\mathrm{R}\mathrm{T}}$和泊松比${\mu }_{\mathrm{L}\mathrm{R}}$、${\mu }_{\mathrm{R}\mathrm{T}}$和${\mu }_{\mathrm{T}\mathrm{L}}$。相关学者对刺槐径向、弦向抗弯弹性模量（MOE）、及剪切弹性膜量${G}_{\mathrm{T}\mathrm{L}}$和${G}_{\mathrm{L}\mathrm{R}}$进行测定^[25]，现有研究对刺槐一类大型乔木各方向弹性参数及剪切模量的比例研究表明^[26]，${E}_{\mathrm{L}}$∶${E}_{\mathrm{R}}$∶${E}_{\mathrm{T}}$ = 18.5∶1.9∶1，${G}_{\mathrm{L}\mathrm{R}}$∶${G}_{\mathrm{T}\mathrm{L}}$∶${G}_{\mathrm{R}\mathrm{T}}$ = 4.3∶3.2∶1，结合马克斯韦尔定理即${\mu }_{\mathrm{i}\mathrm{j}}$/${E}_{\mathrm{i}}$ = ${\mu }_{\mathrm{j}\mathrm{i}}$/${E}_{\mathrm{j}}$，即可求出刺槐的9个弹性常数，本研究所用刺槐弹性常数见表1所示。

表  1  刺槐的弹性常数

Table  1.  Elastic constant of Robinia pseudoacacia

项目 Project	数值 Numerical value		项目 Project	数值 Numerical value		项目 Project	数值 Numerical value
${E}_{\mathrm{L}}$	1.48 × 1010 Pa		${G}_{\mathrm{L}\mathrm{R}}$	4.62 × 108 Pa		${\mu }_{\mathrm{L}\mathrm{R}}$	0.043
${E}_{\mathrm{R}}$	1.57 × 109 Pa		${G}_{\mathrm{T}\mathrm{L}}$	4.47 × 108 Pa		${\mu }_{\mathrm{R}\mathrm{T}}$	0.710
${E}_{\mathrm{T}}$	8.04 × 108 Pa		${G}_{\mathrm{R}\mathrm{T}}$	1.24 × 108 Pa		${\mu }_{\mathrm{T}\mathrm{L}}$	0.029
注：L、R和T分别代表纵向、径向和弦向，E、G和$\mu$分别代表弹性模量、剪切模量和泊松比。Notes: L, R and T represent longitudina, radial and chord direction, respectively, E, G and $\mu$ represent elastic modulus, shear modulus and Poisson’s ratio, respectively.

下载: 导出CSV 

| 显示表格

2.3   脉动风速时程模拟

运用线性滤波法中的AR法对实测风速进行风速时程模拟，为验证脉动风速模拟的有效性，将模拟得到的风速谱（模拟谱）与选取的Davenport谱（目标谱）进行比较如图4所示，模拟谱与目标谱整体基本吻合，具有很小的差异性，可对目标谱进行良好的脉动风速时程模拟，由此可知本研究应用AR法进行脉动风速时程的模拟是可行的。

图  4  模拟谱与风速谱的比较

Figure  4.  Comparison between simulated spectrum and wind speed spectrum

下载: 全尺寸图片 幻灯片

2.4   计算模型

流场为用Ansys Workbench软件建立的20 m × 20 m × 60 m三维流场模型，刺槐模型的建立采用分形L系统理论，该理论为林木结构特征的构建提供了有效的方法^[27]，基于此理论通过编制的命令流文件生成刺槐计算模型，对刺槐保留一级分枝的枝干部分，以现场监测风速点位及刺槐分布作为参考依据，建立多棵刺槐与脉动风场的双向流固耦合模型示意图如图5所示，本研究中双向流固耦合详细设置参数如表2所示。

图  5  刺槐与风场双向流固耦合模型示意图

Figure  5.  Schematic diagram of two-way fluid-structure coupling between Robinia pseudoacacia and wind field

下载: 全尺寸图片 幻灯片

表  2  双向流固耦合详细设置参数

Table  2.  Detailed setting parameters for two-way fluid-structure coupling

项目 Project	参数 Parameter		项目 Project	参数 Parameter
湍流模型 Turbulence model	标准k-ε湍流模型 Standard k-ε turbulence model		计算时间 Computing time	10秒 10 s
流体材料 Fluid material	25 ℃下模拟生成的风速${v}_{0}$ Wind speed ${v}_{0}$ generated by simulation at 25 ℃		时间步数 Time step number	100步 100 steps
分析类型 Analysis type	瞬态分析 Transient analysis		时间步长 Time step length	0.01秒 0.01 s
刺槐模型根部 Robinia pseudoacacia model root	固结约束 Fixed support		模型网格划分方法 Mesh generation method of model	自由网格划分 Free grid division
风场四周边界条件 Circumferential boundary conditions of wind field	无滑移墙面 No slip wall		风场出、入口边界条件 Boundary conditions at the entrance and exit of wind field	静态网格移动模式 Static grid movement mode

下载: 导出CSV 

| 显示表格

3.   数值模拟的有效性验证

3.1   风速归一化对比验证

图6为本研究通过数值模拟技术生成的林内风速变化数据与现场实测风速数据在2 m高度处的水平方向风速归一化对比图。由图可知，实测试验林场中的风速不断降低，且降低幅度越来越小，林中监测的5处风速归一化值分别为0.69${v}_{0}$、0.57${v}_{0}$、0.53${v}_{0}$、0.52${v}_{0}$和0.53${v}_{0}$（${v}_{0}$为林前入口初始风速值）。模拟风速经过每棵刺槐模型时皆有下降趋势，在两棵刺槐模型间风速得到一定恢复，在经过下一棵刺槐模型时风速又被削弱，整体风速呈先下降后稳定趋势，同现场实测风速变化情况一致，数值模拟林后风速减小到0.66${v}_{0}$，与实测数据相差13%。

图  6  水平方向风速归一化对比图

v/v₀为归一化风速，即被削弱的风速与初始风速的比值。下同。 v/v₀ is the normalized wind speed, i.e, the ratio of the weakened wind speed to the initial wind speed. The same below.

Figure  6.  Wind speed comparison in horizontal direction

下载: 全尺寸图片 幻灯片

图7为林中15 m处对2、6、10、14 m 4个高度下实测和模拟风速数据垂直方向风速归一化对比图。由图可知，模拟试验与实测试验在林中垂直方向风速的变化情况接近，实测试验在4个高度下风速归一化的值分别为0.57${v}_{0}$、0.59${v}_{0}$、0.85${v}_{0}$和0.92${v}_{0}$，而模拟试验在4个高度下的归一化风速值分别为0.63${v}_{0}$、0.62${v}_{0}$、0.82${v}_{0}$和0.96${v}_{0}$，在2 m高度处模拟试验误差值最大为6%，模拟试验风速的变化特征无论在水平方向还是垂直方向与现场实测数据皆具有良好的一致性。

图  7  垂直方向风速归一化对比图

Figure  7.  Wind speed comparison diagram in vertical direction

下载: 全尺寸图片 幻灯片

3.2   刺槐枝干风振位移对比验证

图8为林场20 m处实测刺槐枝干风振位移与模拟值的对比情况，刺槐枝干在初始受风荷载阶段具有较大程度的位移，随之摆动情况趋于稳定，实测值在枝干趋于稳定时的位移约为49.5 mm，模拟值在趋于稳定时枝干的位移约为56.5 mm，两者相差仅为7.0 mm。

图  8  刺槐分枝风振位移对比

Figure  8.  Comparison of wind-induced displacement of Robinia pseudoacacia branches

下载: 全尺寸图片 幻灯片

综上比较，本研究所用双向流固耦合模拟涉及的风速和枝干位移两个物理量与现场实测变化规律一致，误差较小，通过此技术模拟刺槐林与风场的相互作用具有一定的可行性，可为不同林带结构下刺槐林合理配置的研究提供科学依据。

4.   结果与分析

4.1   刺槐行数对防风效果的影响

为探究刺槐林行数的防风效果，分别在林中处设置1行、2行、3行、4行及5行刺槐模型作为5组对比试验，每行刺槐模型固定行间距为5 m。目前通常将风速减弱大于20%的长度作为防护林的有效防护距离^[28]，为了更好地评估防护林的防风效果，本研究将5组试验各自的有效防护距离通过后处理提取后绘制如图9所示。在5 m行间距下，随着刺槐行数的增加，刺槐树后的有效防护距离随之增大，5行刺槐的有效防护距离为5.31 m，比单行刺槐有效防护距离1.88 m提升了182.4%，行数对刺槐有效防护距离的提升较为明显。

图  9  不同行数刺槐的有效防护距离

Figure  9.  Effective protection distance of Robinia pseudoacacia with different rows

下载: 全尺寸图片 幻灯片

为探究刺槐行数对水平风速的影响，对模型记录的水平风速进行归一化处理，将5组对比试验林中风速变化情况绘制如图10所示。由图可知，每组风速在经过刺槐模型后都会有所降低，1行至5行5组试验各自使林后风速降低至0.80${v}_{0}$、0.68${v}_{0}$、0.61${v}_{0}$、0.58${v}_{0}$和0.57${v}_{0}$，其中5组试验风速在经过第1行刺槐模型后降低明显，降低了0.2${v}_{0}$的风速，每增加1行刺槐，林后风速依次降低0.12${v}_{0}$、0.07${v}_{0}$、0.03${v}_{0}$和0.01${v}_{0}$。在第5行刺槐时，风速基本不再发生变化达到最小值，单行防护林的防风效率最高，随着刺槐林木行数的增多，风速会被削弱，但削弱幅度越来越小，第4和第5行刺槐风速减弱不再明显。

图  10  不同行数刺槐的归一化风速对比

Figure  10.  Comparison of normalized wind speed for different numbers of Robinia pseudoacacia

下载: 全尺寸图片 幻灯片

4.2   刺槐行间距对防风效果的影响

为研究相邻刺槐行间距对防风效应的影响，本研究以两行刺槐为一整体对行间距分别为0.4H、0.5H、0.7H、1.0H和1.5H（H为刺槐树高）5种不同行间距的刺槐进行模拟，各自行间距下的有效防护距离如图11所示。由图可知，随着刺槐行间距的增加，有效防护距离会缩短，0.4H行间距时有效防护距离最大，且比1.5H行间距下刺槐的有效距离提升30.9%，通过缩短刺槐的行间距可有效提升整体防风效应。

图  11  不同行间距刺槐的有效防护距离

H为树高，下同。H is tree height, the same below.

Figure  11.  Effective protection distance of Robinia pseudoacacia with different row spacing

下载: 全尺寸图片 幻灯片

为探究刺槐行间距对水平风速的影响，绘制5组不同行间距下刺槐的归一化风速对比图（图12）。风速经过5种行间距下第1排刺槐会降低至0.80${v}_{0}$，其中行间距为0.4H的刺槐防风效果最明显，两行刺槐使风速减小至0.65${v}_{0}$，随着行间距的增大，两行刺槐整体的防风效果减弱，行间距为0.5H、0.7H、1.0H和1.5H下的第2行刺槐将初始${v}_{0}$分别降至0.68${v}_{0}$、0.73${v}_{0}$、0.76${v}_{0}$和0.79${v}_{0}$，5组试验第2行刺槐后风速比第1行刺槐风速分别降低了0.15${v}_{0}$、0.12${v}_{0}$、0.07${v}_{0}$、0.04${v}_{0}$和0.01${v}_{0}$。在1.5H的行间距下，第2行刺槐几乎承受了和第1行刺槐相同风速下的荷载作用。

图  12  不同行间距刺槐的归一化风速对比

Figure  12.  Comparison of normalized wind speed of Robinia pseudoacacia with different row spacing

下载: 全尺寸图片 幻灯片

4.3   刺槐排列方式对防风效果的影响

本研究对多棵刺槐模型进行了两种不同排列方式的模拟，分别是方形和交错排列，刺槐的行、列间距均为4 m，两种排列方式皆采用3行刺槐进行模拟，方形排列后两行刺槐皆位于第1行刺槐的正后方，而交错排列的第2行刺槐位于首行相邻两刺槐正中间位置的后方4 m行间距处，第3行刺槐位于首行刺槐正后方8 m行间距处。两种排列方式得到的刺槐在2 m高度处X-Y面风速变化云图（图13）中，Y向为风向方向。方形排列的刺槐林第2、3行刺槐分别位于前一行刺槐遮蔽区内，没有充分发挥自身的防风作用，导致每行刺槐之间仍有浅绿色的高风速带，而交错排列的第1行刺槐充分发挥了自身的防风效果，弥补了方形排列的不足，减少了刺槐列间高风速带，使林中及出口处形成了一个较完整的低风速遮蔽区，提升了刺槐林整体的防风效应。

图  13  刺槐X-Y面风速变化云图

Figure  13.  Cloud image of wind speed variation on X-Y surface of Robinia pseudoacacia

下载: 全尺寸图片 幻灯片

两种排列方式下末排刺槐枝干的风振位移情况如图14所示。不同排列方式对枝干风振位移影响较大，在风载初期，方形排列的刺槐具有较大的风振位移，而交错排列的刺槐分枝风振位移幅值相对较小，原因就是方形排列的刺槐由于内部含有高风速带，而交错排布的刺槐内部形成了一个良好的低风速遮蔽区。随后两种刺槐分枝位移趋于稳定，方形排列与交错排列的刺槐分枝位移情况相近但位移值稍低，原因是沿风向分布方形排列的刺槐，两者的间距更小。综合考虑，为避免刺槐枝干发生折断，应选用交错排列的方式以降低枝干在大风作用下产生的位移幅值，交错排列是多行刺槐林最佳的排列方式。

图  14  两种排列方式下刺槐枝干位移情况

Figure  14.  Displacement of Robinia pseudoacacia branches under two arrangements

下载: 全尺寸图片 幻灯片

综上分析，为切实可靠降低研究地防护林的风害，提升防护林的整体防风效应，根据本研究结果，对后续刺槐防护林建设有以下两点提议：（1）多行刺槐林进行整体防风时，选用5行刺槐为一个整体，行间距选取0.4倍树高并以交错排列方式进行种植；（2）单行防护林的防风效率最高，在改善例如农田等局部小气候，降低风害作用的情况下，可以考虑单行刺槐进行防护且防护距离不宜大于1.5倍刺槐树高，如防护所需距离过大则需增加刺槐行数且两行刺槐间距不宜大于1.5倍刺槐树高。

5.   结　　论

本研究主要通过数值模拟技术，结合现场实测风速数据对不同林带结构下刺槐林防风效应进行研究，主要结论如下：

（1）模拟试验经实测数据对比，水平归一化风速误差为13%，垂直归一化风速误差为6%，刺槐分枝风振位移相对稳定时误差为7 mm，双向流固耦合试验具有较好的仿真精度。

（2）首行刺槐防风效率最高，可使风速降低至0.80${v}_{0}$，2至5行刺槐的阻风效果将风速依次减弱至0.68${v}_{0}$、0.61${v}_{0}$、0.58${v}_{0}$、0.57${v}_{0}$，在5 m行间距下随着刺槐行数的增加风速逐渐降低，降低程度逐渐缓慢，后4行照前1行总风速依次降低0.12${v}_{0}$、0.07${v}_{0}$、0.03${v}_{0}$和0.01${v}_{0}$，刺槐行数增至5行时风速基本不再变化。

（3）相邻两行刺槐随着行间距的增加，整体的防风效果降低，风在经过0.4H、0.5H、0.7H、1.0H和1.5H行间距下的两行刺槐时，第2行风速分别比第1行减小0.15${v}_{0}$、0.12${v}_{0}$、0.07${v}_{0}$、0.04${v}_{0}$和0.01${v}_{0}$，其中0.4H行间距下的两行刺槐防风效果最佳，可使风速降低至0.65${v}_{0}$，当两行刺槐行间距大于1.5H时，经过第2行刺槐的风速几乎不再受第1行刺槐的影响。

（4）方形排列的刺槐沿风场方向林间有多条高风速带，交错排列的刺槐林使林中、后部形成一个良好的低风速遮蔽区。刺槐枝干的最大位移情况受排列方式影响明显，交错排列的刺槐林可明显降低刺槐分枝风振位移变形幅值，从而减少刺槐枝干风损折断现象，提升刺槐防护林的防风效应。