吳 康,周建中,黃顯學
林木風致隨機振動模擬的研究
吳 康,周建中,黃顯學
(北京林業大學 水土保持學院,教育部水土保持與荒漠化防治重點實驗室,北京100083)
采用線性濾波法通過編程模擬了林木的脈動風速時程曲線,并在Ansys內建立其有限元模型,通過施加模擬得到的風荷載獲得了林木風致隨機振動的動力響應。分析表明這種方法對于模擬林木的風致隨機振動是可行的。
林木;風害;隨機振動;動力響應
隨機振動是自然界普遍存在的現象,每當有風襲來,樹木會在風的作用下隨風搖擺,風作用小時擺幅很小,只有樹枝和樹葉在隨風運動,翩翩起舞,樹干幾乎紋絲不動,堅如磐石,風止即止,這是樹隨風的隨機振動現象。而強風來臨時,樹木會產生很大的擺幅,甚至傾倒、折斷。一直以來強風的造成的森林災害對社會產生了重大的經濟損失,嚴重影響了當地的國民經濟增長,同時對森林生態系統的穩定也造成了很大的影響。近年來風害發生的強度和頻率越來越高,引起越來越多的國內外學者注意,國內學者李秀芬[1]、侯倩[2]分析了風害產生的類型及因素,并為森林管理人員提供了防災減災措施。林木風害是由于其風致隨機振動平衡的破壞引起,為防風害有必要了解林木的風致隨機振動,進行風振時程分析是最有效的方法之一。為此在本研究中討論了產生隨機振動的脈動風荷載數值模擬,并建立了林木的有限元模型,通過導入模擬的風荷載得到了其在脈動風荷下的隨機動力響應。
風是空氣分子的運動,由分子從氣壓大的地方流向氣壓小的地方而產生。科學研究者依據得到大量風的實測資料,通過風的順風向時程曲線,得出自然風是由平均風和脈動風這兩部分所組成,平均風周期較長,通常在10分鐘以上,遠遠大于一般的結構自振周期,相當于靜力作用在結構上。脈動風周期短,一般只有幾秒左右,相當于動力作用,激起結構振動[3]。因此在自然風作用下,林木結構的風荷載可認為是由兩部分組成:一是平均風作用下的靜風荷載,二是由于自然風的紊流成分誘發導致的隨機動荷載。林木在隨機振動平衡未被破壞時,主要產生線彈性變形,依據疊加原理,林木的在風載作用下的變形可分為平均風和脈動風分別作用下的疊加。也即:

式(1)中:s為平均風作用下產生的變形,可用靜力學求解。s(t)為脈動風作用下的變形,需要應用動力學求解。林木的隨機振動是由脈動風作用產生,本研究主要討論的就是模擬脈動風荷導致的林木隨機振動。
由于風在流動中與地表摩擦,使得風速隨距離地表高度的增加而增大。Davenport等根據分析實測得到的結果,提出平均風速沿高度變化的規律如下:

式(2)中:v10表示標準高度為10 m處的平均風速(m/s),α為地面粗糙度系數。
Davenport通過分析總結世界上不同地點,不同高度處測得的90多次強風記錄,于1961年提出風速譜的經驗公式:

式(3)中:K為地面粗糙度系數,v10表示標準高度為10 m處的平均風速(m/s),n表示脈動風頻率,x=1 200,其風速譜峰值不隨高度而變化。除此之外還有其他常用風速譜如Kaimal譜和Simiu譜[4],它們的風速譜都隨高度變化,在本研究中脈動風模擬采用Da-venport譜。
脈動風是一種隨機干擾,在模擬之前必須確定其概率分布及功率譜密度,功率譜密度是脈動風最重要的統計特征,能反映出在某一頻率上脈動風的能量大小,又因功率與荷載成平方關系,因此在模擬中必然會存在如何分解功率譜,從而得到荷載表達式的問題[5]。國內外大量研究表明脈動風可作為均值為零的高斯過程及平穩隨機過程來考慮,具有很明顯的各態歷經性,對于n個具有零均值的平穩高斯過程[6],其譜密度函數矩陣可表示為:

式(4)中 s11(ω),s22(ω)等可由式 (2)求得,且對 于 Davenport譜 有:s11(ω)=s22(ω)=…=snn(ω);對 S(ω)做 Cholesky 分 解 有:S(ω)=H(ω)·H*(ω)T,式中

H*(ω)T為 H(ω)的共軛轉置。
脈動風作用下,在林木結構表面風向及風速在不同位置處并不同步,有的甚至是相互獨立無關的,因此,模擬林木結構的脈動風壓須考慮其空間相關性。根據Shiotani的建議,在本次模擬研究中僅考慮豎向相關性,相干函數表達式為

根據空氣動力學原理,風壓與風速的基本關系為:

式(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。
對于建立林木動力學模型主要有兩種,一種是基于理論分析的SIA 方法[7],這種方法用非常簡單的圖形表示法描述林木的形狀,具有應用廣泛、不需要電腦輔助計算等優點,缺點是不能適用于所有樹種,不能使研究者了解其內部動力機理。另外一種基于數學描述語言精確描述幾何形狀的方法,稱之為L-system 方法 ,它能描述出復雜幾何形狀的樹種,可用有限元方法分析,缺點是求解困難甚是得不到解。鑒于以上兩種方法的優缺點,在本研究中采用的是模態分析方法求解,它是介于以上二者之間的方法。模態分析主要用于確定結構體系的振動特性,同時也是其他動力學分析的基礎,如諧響應分析、瞬態動力學分析、譜分析等。進行林木結構模態分析的有限元軟件采用國內外廣泛應用的Ansys。
本研究模擬實例為北京郊區某公路旁的進行過現場試驗的新疆楊,新疆楊樹干和樹枝結構細長,主干直立向上,呈圓柱形。樹根入土深,與土壤有很強的錨固力,抗風性好。常作為風景樹、行道樹及綠化樹。根據新疆楊以上特性,在Ansys模型中其樹干及第一主分支采用Beam188模擬[9],Beam188基于Timoshenko梁結構理論,考慮了剪切變形的影響,適用于分析從細長到中等粗短的梁結構[10]。同時為研究方便,模型中暫不考慮樹葉和其它次分枝的影響,根部和土壤之間考慮為固結。樹干、樹枝質量密度為均勻分布,全樹彈性模量為定值。樹枝作為懸臂梁固結于樹干,而樹干也同樣作為懸臂梁固結于地面。新疆楊有限元模型及主要參數見表1。

表1 新疆楊有限元主要物理參數Table 1 Main physical parameters of P. opulus bolleana

圖1 新疆楊有限元模型Fig.1 FEM of P.opulus bolleana
采用線性濾波法原理模擬脈動風,線性濾波法又稱白噪聲濾波法,它將隨機過程抽象為滿足一定條件的白噪聲,然后經某一假定系統進行適當變換而擬合出具有隨機性、時間相關性、空間相關性的風速時程模型。線性濾波法中的自回歸(Auto—Regressive,AR)模型因其計算量小、速度快的特點,被廣泛應用于隨機振動和時間系列分析中[13]。依據文獻[5],[13-15]并結合林木結構的特點,編制了基于Matlab的模擬程序。本研究中脈動風AR模擬使用的主要參數如表2所示。因篇幅限制,只選取了新疆楊模型中2、3結點風速時程曲線,從圖2、圖3中可以看出不同高度處的脈動風速大小不同,在互不相同的時刻達到各自峰值,體現出空間相關性的影響,其值的變化范圍在0附近波動,表現出平穩隨機過程的特性。圖4中,將節點3的Davenport脈動風速譜與通過模擬獲得的脈動風速譜用雙對數坐標形式比較,可以看出二者相當吻合,表明了對于給定風速譜密度的隨機過程,采用AR模型能夠很好的實現人工模擬,并且具有較高的準確性和精度。

表2 主要模擬參數Table 2 Main parameters of AR model

圖2 節點2風速時程曲線Fig.2 Wind speed-time curve of Node 2

圖3 節點3風速時程曲線Fig.3 Wind speed-time curve of Node 3

圖4 節點3模擬譜與目標譜對比Fig.4 Comparison between simulation spectrum and target spectrum of Node 3

圖5 節點3風壓時程曲線Fig.5 Wind pressure-time curve of Node 3
在Matlab中輸出各節點風荷載時程數據,施加到Ansys中新疆楊有限元模型相應節點上,并進行瞬態動力學分析,就可以獲得新疆楊的各節點處的位移、速度、加速度時程曲線,為使圖清晰可見,速度及加速度只提取了前15 s的時程曲線。從圖6至圖8可以看出,在脈動風荷載作用下,新疆楊風致振動的位移、速度、加速度也呈現出明顯的隨機性特點。

圖6 節點3位移時程曲線(m)Fig.6 Displacement-time curve of Node 3

圖7 節點3速度時程曲線(m/s)Fig.7 Speed-time curve of Node 3

圖8 節點3加速度時程曲線(m/s-2)Fig.8 Acceleration-time curve of Node 3
本研究利用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.
A simulation study on wind-induced random vibration of trees
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.
forest trees; wind damage; random vibration; dynamic response
S712;S761.2;TU312
A
1673-923X (2012)08-0042-04
2012-04-29
國家自然科學基金項目(30872071)