-
在众多自然灾害中，风荷载作用所引起的林木风倒是最常见一种，同时风倒也是形成林隙的主要原因，逐渐成为森林生态的一个扰动因子。风倒现象主要表现为干折、根折、连根拔起等3种方式。全世界每年非灾害性风倒造成的木材减产在15%以上[1]。近年来，人类对自然环境的破坏日益严重，导致全球气候环境的急剧变化，风倒引发的森林破坏频繁发生。因此，展开树木风倒机制的相关研究日益受到科学研究人员的重视，研究结果对于森林生态系统和木材生产的理论研究和生产实践都有极其重要的意义。迄今为止，国内外很多专家学者开展了林木风倒机制方面的研究，并取得了丰富的研究成果。其中，李国旗等[2]研究了位于不同风压和高度处树木的应力分布规律。陈少雄等[3-4]、陈士银[5]基于统计和野外调查的方法分别研究了泡桐Paulownia tomentosa，桉树Eucalyptus robusta，木麻黄Casuarina equisetifolia等树种的抗风性与个体特征之间的关系。代力民等[6]从生态学角度出发研究林隙的形成原因。宋晓鹤[7]，王琳[8]，赖秋明[9]分别从静力学和动力学角度建立了风倒的模型，并展开了相应的分析。吴志华等[10]进行了抗风性状对抗风效果的影响规律的研究。段溪[11]运用层次分析法研究了杭州湾地区亚热带泥质海岸的68种防护树种的抗风性能。英国、加拿大、美国等一些学者对樟子松Pinus sylvestris var. mongolica等树种的风倒机制进行了深入广泛的研究，以动力学理论，结合野外试验中林木和土壤力学参数的测定，地形数字高程模型和实物模型风洞试验，获得了复杂地形下和不同林分密度时的风场分布，建立了一系列力学方程来描述林木的抗风性[12]，对风倒的力学机制也进行了相应的阐述，形成若干关于林木风倒机制的模型，如MC2[13]和HWIND[14]等。尽管这些研究代表了风倒机制研究的国际先进水平，但主要是基于材料力学方法的线性力学模型[15]。虽然描述了风倒机制的许多方面，但受研究条件等的限制在系统设计时对关键环节进行了简化，造成模型与实际情况之间的误差。作者在已有研究成果的基础上，建立了林木发生风倒的非线性动力学模型[16]，并且利用Von-Karman风谱模型从频率域内分析了林木风倒机制[17]，验证了模型的有效性。相比其他模型，这个模型更接近实际物理系统且考虑了实际林木发生风倒的一些不确定性影响因素。但模型也存在一些缺陷，首先，并未考虑树冠的形状、树木个体差异、枝条和树叶的作用等对风倒的影响，其次，也未考虑风荷载作用下的“流固耦合”现象，这也是后续研究的方向和内容。本研究在前期的研究基础上，着重考虑孤立木在稳定风速和阵风作用下所引起的倾覆，通过严格的数学推导给出风倒发生的条件，并以西加云杉Picea sitchenrsis为例，进行了数值计算，给出树木发生风倒的标准风速值。研究结果对森林抗风防护建设和经营提供现实的指导意义。
HTML
-
当风荷载作用于树木的时候，空气动力将会对树杆和树枝产生一个对于基础部分的倾覆力矩，由根部所产生的力矩与其抵抗(抵抗力矩)。我们将倾覆力矩定义为mT(θ)，抵抗力矩定义为mR(θ)，并且假设树为刚性体，我们可以得到如下的模型：
式(1)中：I0表示树木相对根部的转动惯量，θ表示树杆在风力作用下相对于垂直位置的倾角。如图 1所示。以西加云杉为例，通过对其根部的抵抗力矩进行测量，发现抵抗力矩取决于根部铰接部分的抵抗矩、土壤的压力、迎风面根系的拉力、根系土壤的质量等等，如图 2所示。
Figure 1. Basic force diagram of trees
Figure 2. Relationship between whole resisting moment and angle θ
由图 2可知：总体抵抗力矩是关于倾角θ的函数，是很多树木拉力试验中的典型曲线。并且当θ很小的时候，总体抵抗力矩急剧上升，在2°到5°左右，抵抗力矩达到最大值，在10°的时候，抵抗力矩降到其最大抵抗力矩的一半左右。基于本研究的目的，可以将抵抗力矩模拟成一个二次方程，在θ=θm的时候有最大抵抗力矩mR(θ)=Mm，在θ=0时，定义mR(θ)=M0。因此可以得到关于mR(θ)的方程：
在0°到0.5°之间，由于其mR(θ)曲线陡然增加，近似用割线代替，从坐标轴上选取(0, M0)作为起点，得抵抗力矩模型曲线如图 3所示。
Figure 3. Mode diagram of resisting moment mR(θ)
-
如图 1所示，空气动力对于树木一个微小的单元作用力为νn2dA，νn即对单元作用的风速，dA表示树干上具有代表性的单元面积。假设风具有连续的水平速度u0，当树木垂直的时候，风力对树木根部产生的弯矩为Fm(u0)，当树木在风力作用下发生倾斜，倾角为θ的时候，其弯矩变为Fm(u0)cos2(θ)，同时重心也将发生偏移，自身重力将产生一个倾覆力矩，定义为Wmsin(θ)。Wm表示θ=90°时，树木自身重力对根部所产生的力矩。则方程变为：
当根部抵抗力矩较小时，可以近似地认为：sinθ≈θ，cos(θ) ≈1-θ2/2，因此式变为：
且θ较小时，角加速度θ″(t)可用旋转角速度ω=θ′(t)表示为：θ′(t)=ω′(t)=ωdω/dθ，则式(4)变为:
式(5)综合了树木旋转时的动能，同时角速度ω也是θ的函数，则有:
由式(6)可知，当θ→0时，树木在风载作用下开始倾斜，此时ω=0，C=0。则式(6)可表示为:
由以上分析可知，式(7)只有当θm＜5°时才成立，此时θ可以用θm线性表示：
则式(7)可变为：
当φ取正值的时候，式(9)右边如果也保持正值，那么角速度ω也将保持正值，这时随着θ的增大，树木将继续倾斜，根据式(9)可以得到如图 4的函数图像。从图 4中可知，当满足条件(10)时式(9)没有正根：
Figure 4. Function diagram of the formula (a show the condition (10) is satisfied, b show the condition (10) is dissatisfied)
当u0足够大能满足不等式(10)时，树木风倒便会发生。
针对图 1中的树木模型，将树干模拟成一个圆柱，其半径为r0，高度为h，则作用于单元dy上的荷载为：ρacDr0u02dy。其中：ρa表示空气密度，cD表示无量纲的阻力系数(与树干的半径有关)，u0表示作用于树干的平均连续稳定风速。则风荷载相对于根部所产生的弯矩为：
实际上，树木在风载作用下产生倾斜，树木重心位置会发生偏移，重力矩Wm将会对弯矩Fm(u0)产生一个附加的作用，这个附加作用可以通过对式(11)进行修正来体现。
式(12)中：ε=0.4ρgh3/(Er02)，ρ表示树木的密度，E表示弹性模量。式(12)相当于式(11)中的阻力系数cD被提高了1+ε倍。
-
基于2.1的理论分析，以西加云杉为例进行数值计算，选取相应参数为：h=19 m，r0=0.1 m，ρ=900 kg·m-3，ρa=1.2 kg·m-3，cD=0.6。可计算得到抵抗力矩：θm=2.5°=0.044 rad，M0=8 × 103 N·m，Mm=12.8 × 103 N·m。通过对根系强度的可变性进行分析，对上述大小的西加云杉树木，其根倒弯矩值可能在10×103 N·m到50 × 103 N·m，并且能给出图 1模型树木对根部的转动惯量:
对于连续稳定风速引起的风倒，满足条件式的最小风速为：
可得u0=29.9 m·s-1，由于式(13)并未考虑根系和树枝惯性矩的影响，因此计算得到风倒发生时的风速可能是被低估了的风倒标准风速。在临界风速u0=29.9 m·s-1作用下，树木至少需要2 s的时间让倾斜达到θm，而需要4.5 s的时间让树木旋转角速度ω达到最小值。总之，树木如果受到连续稳定风速作用并产生弯矩Fm(u0)，当满足条件式(10)，那么树木将在连续的风速作用下发生倾覆。
除西加云杉以外，还对稻类植物、悬铃木Platanus acerifolia也进行了相应的分析，基于各类参数计算得到稻类植物发生风倒的临界风速u0=11.6 m·s-1，悬铃木发生风倒的临界风速u0=22.8 m·s-1，而根据所查资料[18]，3种林木实际发生风倒的数值分别：26.3, 8.8和19.1 m·s-1。注意到，理论计算值和实际数据比较吻合，但是计算值看起来要比实际数据大一些。究其原因是：实际发生风倒与林木个体参数、土壤条件、根系条件以及枝条、树叶与风荷载所形成的“流固耦合”效应有关。以典型的小麦Triticum aestivum来说，一般情况下都有2到6个分支，假设分支数目为n，如果每一个分支所承受的弯矩都有根系来支持，就降低了理论计算所得到的风倒临界风速。再比如“流固耦合”效应导致流体边界的变化，从而对流动产生反作用。不考虑“流固耦合”作用会使得计算结果是偏于安全的，也即“流固耦合”效应导致实际发生风倒的临界风速降低。根系在风倒过程中有至关重要的地位，长期的连续风荷载作用，会降低根系的支持力和抵抗弯矩的能力，以至于林木可能在较小的风速下发生风倒。因此，理论计算值与实验数据结果比较吻合，同时风的动态作用可以得到合理的解释。
2.1. 风倒原理分析
2.2. 数值计算
-
实际上树木更容易在强大的阵风作用下发生风倒，假设其作用时间为τ，在风速为u0连续风载作用下根部产生倾覆弯矩Fm(u0)，则t=0时刻阵风所引起的冲量矩为最大值，树木的角动量为[I0ω]0-0+=(Mg-M0)τ，其中Mgτ表示阵风所产生的冲击荷载引起的冲量矩，M0τ表示阵风引起根部的抵抗力矩。则树木的角速度为：
式(７)变为：
当t=0+时，积分常数C=I0ω02/2，阵风的作用，导致产生了力矩Mg，将Mg表达成连续稳定风速作用下弯矩Fm(u0)的函数：
式(17)中，${G_{\rm{v}}} = \frac{{\hat u}}{{\bar u}} $为阵风系数，${\hat u} $觠指阵风风速极值，u指阵风风速的平均值。让θ=φθm，由式(16)可得:
显然，式(18)右边必须保持正值，树木才会继续发生倾斜，在其他所有参数已知的情况下，式(18)是关于φ的一个三次函数，假设其梯度为0，则退化为关于φ的二次函数，在$\varphi = {\varphi _{\rm{r}}} = 1 + {\left( {\frac{{{M_{\rm{m}}}-{F_{\rm{m}}}\left( {{u_0}} \right)}}{{{M_{\rm{m}}}-{M_{\rm{0}}}}}} \right)^{\frac{1}{2}}} $处可以得到根的最大值。要保证式(18)右边是正值，则必须保证：
成立。由式可以得到相应的判断条件：
可以看出：式(20)是对式(10)的很好补充。其物理含义是：式(10)指出当风速足够大的时候，在连续稳定风速的作用下，风倒便会发生。式(20)指出当连续风速不够大，以至于不能够引起风倒的情况下，连续的风载将和阵风共同作用下引起风倒[阵风系数GV足够大，并且作用时间τ足够长，能满足式(19)]。式(19)取决于如下的无量纲量：Fm=Fm(u0)/M0；Mm=Mm/M0；I0=I0θm/M0τ2。则条件表示的破坏阵风系数必须满足:
-
由无量纲化的风倒破坏式(21)可知，风倒条件取决于所定义的无量纲量，且都是风速u0和阵风作用时间τ的函数。以西加云杉为例进行数值计算，所选参数与2.2中的相同，则无量纲参数变为：
当阵风风速最大值$ {\hat u}$觠和作用时间τ子满足式(21)，风倒便会发生，可得到风倒的风速临界值：
图 5表示出阵风风速${\hat u} $和连续风速u0(可以由阵风风速的平均值u代替)在不同的时间段τ内的关系。从图中可以看到，当u0接近于临界风速u0=29.9 m·s-1的时候，所有曲线都将终止。如果将1 h风速的作用时间按照τf时间段来划分，可以得到这段时间内发生的最大风速。
Figure 5. Function relationship between gust wind and continuous wind about windthrow (in the case of τ=0.2, 0.5, 1, 2, 5, 10s)
式(23)中，σ/u表示紊流强度，对于独立树木，σ/u一般在0.2左右；而对于森林，σ/u会更高一些。依据式(23)，给出不同的τ值，可以得到如图 6所示的一族直线，这些直线穿过图 5中的曲线，相应产生交点，将交点连在一起，如图 7所示，得到风倒发生的风速标准线。在这个直线下方，不可能出现阵风风速${\hat u} $。由图 7可知：风倒可能在较小风速17.7 m·s-1就可能发生。而对于阵风来说，在风速为27.5 m·s-1，作用时间5 s的情况下风倒就可能发生。尽管5 s时间给树木以冲击荷载看起来比较短，但于树木来说，从受风载作用开始倾斜到其抵抗弯矩最大值位置，大约需要2 s。所以，可以寻找作用时间更短的阵风，对于阵风风速为35.1 m·s-1，作用时间τ=1 s叠加在风速为21.0 m·s-1的风速上就会引起风倒。
Figure 6. Diagram of average wind speed when turbulence intensity is 0.2
Figure 7. Average wind speed diagram of damage when turbulence intensity is 0.2, τ=1, 2, 5s
-
林木的破坏模式有2种，一种是破坏值超过林木材料的允许值而发生的折断现象，另外一种破坏值超过根系的支持力而发生的连根拔起现象。本研究综合考虑在稳定风速和阵风作用下树木的倾覆问题。将树木风倒用一个简单的模型来模拟根系，从动力学角度建立了风倒的动态方程，得到在稳定风速下树木发生风倒的条件。同时将阵风作用模拟为一个冲量荷载作用在树木上部，建立阵风作用下树木风倒的条件。以西加云杉为例，计算了树木在稳定风速和阵风作用下破坏的标准值，将阵风风速的作用时间与平均风速作用时间比较，得到树木发生风倒的极限风速。但是，实际林木风倒临界风速值低于文章所计算出的理论值，这是由于林木树冠形状、个体参数差异、枝条、树叶与风荷载之间的“流固耦合”效应以及复杂的非线性动力学行为都导致实际发生风倒的风速值偏低。这也为今后的研究提供了方向和思路。总之，文章从理论方面的研究结果对森林的抗风防护以及木材的生产实践有很重要的现实指导意义。
DownLoad: