林木风致随机振动模拟的研究

(北京林业大学水土保持学院，教育部水土保持与荒漠化防治重点实验室，北京100083)

采用线性滤波法通过编程模拟了林木的脉动风速时程曲线，并在Ansys内建立其有限元模型，通过施加模拟得到的风荷载获得了林木风致随机振动的动力响应。分析表明这种方法对于模拟林木的风致随机振动是可行的。

随机振动是自然界普遍存在的现象，每当有风袭来，树木会在风的作用下随风摇摆，风作用小时摆幅很小，只有树枝和树叶在随风运动，翩翩起舞，树干几乎纹丝不动，坚如磐石，风止即止，这是树随风的随机振动现象。而强风来临时，树木会产生很大的摆幅，甚至倾倒、折断。一直以来强风的造成的森林灾害对社会产生了重大的经济损失，严重影响了当地的国民经济增长，同时对森林生态系统的稳定也造成了很大的影响。近年来风害发生的强度和频率越来越高，引起越来越多的国内外学者注意，国内学者李秀芬[1]、侯倩[2]分析了风害产生的类型及因素，并为森林管理人员提供了防灾减灾措施。林木风害是由于其风致随机振动平衡的破坏引起，为防风害有必要了解林木的风致随机振动，进行风振时程分析是最有效的方法之一。为此在本研究中讨论了产生随机振动的脉动风荷载数值模拟，并建立了林木的有限元模型，通过导入模拟的风荷载得到了其在脉动风荷下的随机动力响应。

1 风的基本组成

风是空气分子的运动，由分子从气压大的地方流向气压小的地方而产生。科学研究者依据得到大量风的实测资料，通过风的顺风向时程曲线，得出自然风是由平均风和脉动风这两部分所组成，平均风周期较长，通常在10分钟以上，远远大于一般的结构自振周期，相当于静力作用在结构上。脉动风周期短，一般只有几秒左右，相当于动力作用，激起结构振动[3]。因此在自然风作用下，林木结构的风荷载可认为是由两部分组成：一是平均风作用下的静风荷载，二是由于自然风的紊流成分诱发导致的随机动荷载。林木在随机振动平衡未被破坏时，主要产生线弹性变形，依据叠加原理，林木的在风载作用下的变形可分为平均风和脉动风分别作用下的叠加。也即：

式（1）中：s为平均风作用下产生的变形，可用静力学求解。s(t)为脉动风作用下的变形，需要应用动力学求解。林木的随机振动是由脉动风作用产生，本研究主要讨论的就是模拟脉动风荷导致的林木随机振动。

1.1 平均风速

由于风在流动中与地表摩擦，使得风速随距离地表高度的增加而增大。Davenport等根据分析实测得到的结果，提出平均风速沿高度变化的规律如下：

式（2）中：v10表示标准高度为10 m处的平均风速（m/s），α为地面粗糙度系数。

1.2 脉动风速谱

Davenport通过分析总结世界上不同地点，不同高度处测得的90多次强风记录，于1961年提出风速谱的经验公式：

式（3）中：K为地面粗糙度系数，v10表示标准高度为10 m处的平均风速（m/s），n表示脉动风频率，x=1 200，其风速谱峰值不随高度而变化。除此之外还有其他常用风速谱如Kaimal谱和Simiu谱[4]，它们的风速谱都随高度变化，在本研究中脉动风模拟采用Da-venport谱。

1.3 功率谱函数

脉动风是一种随机干扰，在模拟之前必须确定其概率分布及功率谱密度，功率谱密度是脉动风最重要的统计特征，能反映出在某一频率上脉动风的能量大小，又因功率与荷载成平方关系，因此在模拟中必然会存在如何分解功率谱，从而得到荷载表达式的问题[5]。国内外大量研究表明脉动风可作为均值为零的高斯过程及平稳随机过程来考虑，具有很明显的各态历经性，对于n个具有零均值的平稳高斯过程[6]，其谱密度函数矩阵可表示为：

式（4）中 s11(ω)，s22(ω)等可由式 (2)求得，且对于 Davenport谱有：s11(ω)=s22(ω)=…=snn(ω)；对 S(ω)做 Cholesky 分解有：S(ω)=H(ω)·H*(ω)T，式中

1.4 脉动风压相干函数

脉动风作用下，在林木结构表面风向及风速在不同位置处并不同步，有的甚至是相互独立无关的，因此，模拟林木结构的脉动风压须考虑其空间相关性。根据Shiotani的建议，在本次模拟研究中仅考虑竖向相关性，相干函数表达式为

1.5 脉动风荷载

式（5）中：V = v + v (t )，ρ为空气质量密度。W=ρV2=ρ[v+v(t)]2=ρ[v2+2w(t)+ v2(t)]；平均风压：w=ρv2；脉动风压：w(t)=ρ[2v2v(t)+ v2(t)]；

根据每个节点的受风荷面积Si，乘以节点处的脉动风压时程wi(t)就可得到脉动风荷载的时程曲线也即：Fi=wi(t)·Si。

2 林木有限元模拟方法

对于建立林木动力学模型主要有两种，一种是基于理论分析的SIA 方法[7]，这种方法用非常简单的图形表示法描述林木的形状，具有应用广泛、不需要电脑辅助计算等优点，缺点是不能适用于所有树种，不能使研究者了解其内部动力机理。另外一种基于数学描述语言精确描述几何形状的方法，称之为L-system 方法，它能描述出复杂几何形状的树种，可用有限元方法分析，缺点是求解困难甚是得不到解。鉴于以上两种方法的优缺点，在本研究中采用的是模态分析方法求解，它是介于以上二者之间的方法。模态分析主要用于确定结构体系的振动特性，同时也是其他动力学分析的基础，如谐响应分析、瞬态动力学分析、谱分析等。进行林木结构模态分析的有限元软件采用国内外广泛应用的Ansys。

3 实例模拟

3.1 Ansys有限元模型

本研究模拟实例为北京郊区某公路旁的进行过现场试验的新疆杨，新疆杨树干和树枝结构细长，主干直立向上，呈圆柱形。树根入土深，与土壤有很强的锚固力，抗风性好。常作为风景树、行道树及绿化树。根据新疆杨以上特性，在Ansys模型中其树干及第一主分支采用Beam188模拟[9]，Beam188基于Timoshenko梁结构理论，考虑了剪切变形的影响，适用于分析从细长到中等粗短的梁结构[10]。同时为研究方便，模型中暂不考虑树叶和其它次分枝的影响，根部和土壤之间考虑为固结。树干、树枝质量密度为均匀分布，全树弹性模量为定值。树枝作为悬臂梁固结于树干，而树干也同样作为悬臂梁固结于地面。新疆杨有限元模型及主要参数见表1。

表1 新疆杨有限元主要物理参数Table 1 Main physical parameters of P. opulus bolleana

3.2 Matlab脉动风荷模拟

采用线性滤波法原理模拟脉动风，线性滤波法又称白噪声滤波法，它将随机过程抽象为满足一定条件的白噪声，然后经某一假定系统进行适当变换而拟合出具有随机性、时间相关性、空间相关性的风速时程模型。线性滤波法中的自回归(Auto—Regressive，AR)模型因其计算量小、速度快的特点，被广泛应用于随机振动和时间系列分析中[13]。依据文献[5],[13-15]并结合林木结构的特点，编制了基于Matlab的模拟程序。本研究中脉动风AR模拟使用的主要参数如表2所示。因篇幅限制，只选取了新疆杨模型中2、3结点风速时程曲线，从图2、图3中可以看出不同高度处的脉动风速大小不同，在互不相同的时刻达到各自峰值，体现出空间相关性的影响，其值的变化范围在0附近波动，表现出平稳随机过程的特性。图4中，将节点3的Davenport脉动风速谱与通过模拟获得的脉动风速谱用双对数坐标形式比较，可以看出二者相当吻合，表明了对于给定风速谱密度的随机过程，采用AR模型能够很好的实现人工模拟，并且具有较高的准确性和精度。

图4 节点3模拟谱与目标谱对比Fig.4 Comparison between simulation spectrum and target spectrum of Node 3

3.3 新疆杨风致振动随机响应

在Matlab中输出各节点风荷载时程数据，施加到Ansys中新疆杨有限元模型相应节点上，并进行瞬态动力学分析，就可以获得新疆杨的各节点处的位移、速度、加速度时程曲线，为使图清晰可见，速度及加速度只提取了前15 s的时程曲线。从图6至图8可以看出，在脉动风荷载作用下，新疆杨风致振动的位移、速度、加速度也呈现出明显的随机性特点。

4 结论及展望

本研究利用AR模拟方法实现了林木脉动风荷载的模拟，而修改AR模拟参数还可以对处于不同地区、不同地貌的林木进行风荷模拟，同时对新疆杨的有限元模型修改物理参数及建模的几何形状后又可以推广到其他的树种，将以上二者结合起来就可以实现不同地区、不同树种的风致随机振动模拟。在有限元模型中考虑树叶、其它次要树枝后，能使模拟更加趋于真实。对于实际为林木采用的防灾减灾措施，如除叶剪枝、施加减震器等，都能用是否建立枝，叶模型、改变林木阻尼(质量阻尼、刚度阻尼)等来实现。在模型中设置生死单元后还可以将风致振动过程中断枝的对林木整体运动的影响体现出来，这种情况在实际试验中是比较难实现的。综上所述进行林木风振时程分析可以使研究者能够全面的了解林木的风振响应特征以及直观的反应出风振控制的效果，而本研究采用的人工模拟风荷载时程曲线方法和程序适用于林木结构的脉动风荷模拟，为进行林木风振时程分析提供了新思路。

[1] 李秀芬,朱教君,王庆礼,等.森林的风/雪灾害研究综述[J].生态学报,2005.1,25(1):148-157.

[2] 侯倩,李意德,康文星,等.海南热带滨海城市防台风防护林树种的选择[J].中南林业科技大学学报,2011,31(5):184-191.

[3] 张相庭.结构风压和风振计算[M].上海：同济大学出版社,1985.

[4] Davenport A G. The Relationship of W-ind Structure to Wind Loading[M].WE-BE,1963.

[5] 陆飞,李爱群,程文瀼,等.脉动风荷模拟中的几点问题的探讨[J].特种结构,2002,19(3):18-20.

[6] 王之宏.风荷载的模拟研究[J].建筑结构学报,1994,15(1):44-52.

[7] Wessolly L, Erb M. Handbuch der bau-mstatik und baumkontrolle[M]. Patzer,1998.

[8] Przemyslaw Prusinkiewicz, Aristid Lind-enmayer, Hanan J S, et al.The Algorithmic Beauty of Plants[M]. Springer, 1996.

[9] Jonsson M J, Foetzkia, Kalberer M, et al. Root-soil rotation stiffness of Norway spruce (Picea abies(L.) Karst) growing on subalpine forested slopes[J]. Plant Soil,2006,285:267-277.

[10] Ansys Inc. ANSYS Theory Reference[M]. Canonsburg, 2002.

[11] Richard B, Cai Zhiyong, Charlie G C, et al. Wood Handbook[M].Forest Products Laboratory, 2010.

[12] Moore J.R, Maguire D.A. Natural sway frequencies and damping ratios of trees: concepts, review and synthesis of previous studies[J].Trees,2004, 18:195-203.

[13] 刘锡良,周颖.荷载的几种模拟方法[J].工业建筑,2005,35(5):81-84.

[14] 李元齐,董石麟.大跨度空间结构风荷载模拟技术研究及程序编制[J].空间结构, 2001,7(3):3-11.

[15] 康文星,赵仲辉,邓湘雯.杉木林冠层的动力效应及动能传递规律的研究[J],中南林业科技大学学报,2007,27(2):1-6.

[16] 朱瑞兆.风压计算的研究[M].北京：科学出版社,1976.

WU Kang, ZHOU Jian-zhong, HUANG Xian-xue
(Key Lab. of Soil and Water Conservation & Desertification Combating , Ministry of Education, School of Soil and Water Conservation,Beijing Forestry University , Beijing 100083 , China)

The turbulent wind velocity time-history curve of forest trees has been simulated by using linear filtering method and programming, the finite element model（FEM） of the forest trees was built up within the Ansys, and the dynamic response data of the forest trees’ wind-induced random vibration were obtained by applying the simulated wind load. It shows that this method is feasible to simulate forest trees wind-induced random vibration.

吴康（1986—），男，湖南常德人，硕士，主要从事混凝土结构、振动控制、仿生结构方面的研究

周建中（1970—），男，湖南衡东人，副教授，博士，主要从事混凝土结构、异形结构及振动控制、仿生结构、木结构方面的研究