基于ABAQUS有限元法的空洞樹干強度分析

吳貽軍，姚劍飛，吳 俊，王福利，邵卓平

(1.黃山風(fēng)景區(qū)管委會 園林局，安徽 黃山 245800；2.安徽農(nóng)業(yè)大學(xué) 林學(xué)與園林學(xué)院，安徽 合肥 230036)

古樹名木和高大行道樹在經(jīng)歷了長期的風(fēng)雨洗禮后，均可能存在因腐朽造成的空洞現(xiàn)象，它們在強風(fēng)雪載荷下可能發(fā)生斷裂倒伏而造成行人生命安全隱患和財產(chǎn)損失。自然界中，受風(fēng)雪載荷作用，樹干會發(fā)生彎曲、扭轉(zhuǎn)甚至彎扭組合變形，樹干的彎曲斷裂就是一種常見的破壞類型，尤其是在中齡林和成熟林中表現(xiàn)更為突出[1-3]。另一方面，樹木是一種多胞層狀天然復(fù)合材料，樹干在某種程度上是由同心圓筒狀(年輪層次)的薄層組成，在宏觀力學(xué)行為上可視之為一種典型的圓柱對稱性復(fù)合材料，并在強度和剛度方面展示了強烈的各向異性性質(zhì)，表現(xiàn)為木材沿順紋方向的抗拉強度可達150～300 MPa，但橫向抗拉強度和順紋抗剪強度僅為前者的1/40～1/20和1/20左右[4]，因此，由風(fēng)力雪載引起的橫向張力與剪切力極易引發(fā)樹干沿順紋方向的劈裂，但當(dāng)樹木偏冠嚴重或樹干因腐朽存在大空洞時，樹干常常是先發(fā)生扭轉(zhuǎn)劈裂或橫裂使之喪失了整體剛度后再折斷[5-6]。

目前，采用力學(xué)方法分析風(fēng)雪載荷下的樹木強度越來越受到重視[7-10]，同時認為，當(dāng)受風(fēng)雪載荷下樹干最外層的軸向應(yīng)力超過了材料的臨界應(yīng)力強度時就會發(fā)生破壞[11-13]。有限元方法因其數(shù)值模擬便于多工況加載分析以及精準(zhǔn)的計算精度等優(yōu)勢而逐漸被研究者應(yīng)用于樹木抗風(fēng)性研究和指導(dǎo)林業(yè)生產(chǎn)[14-18]。本研究利用已測定的樹木生材力學(xué)性質(zhì)數(shù)據(jù)為基礎(chǔ)，采用有限元數(shù)值模擬方法，利用Abaqus有限元軟件建立樹干的三維模型，分析樹干因腐朽空洞造成的不同內(nèi)外徑比，以及受不同彎扭矩比下的樹干最外層應(yīng)力變化，判斷不同程度空洞樹木的破壞行為和方式，為樹木安全性評估提供一定的理論支持和依據(jù)。

1 材料與方法

1.1 彎、扭載荷作用的空心樹干受力模型建立

樹干由樹皮、形成層、木質(zhì)部等組成，樹干木質(zhì)部內(nèi)部含有節(jié)子、樹脂道、管孔等不同的結(jié)構(gòu)組分，應(yīng)用有限元軟件分析時，需在軟件中對實物建立數(shù)值模擬的幾何模型。樹皮主要為運送養(yǎng)料，形成層多為單列細胞層起分生作用，對樹干的力學(xué)強度影響微弱。用于測定樹木樹干材料力學(xué)性質(zhì)的試件取之于木質(zhì)部，包含了木質(zhì)部內(nèi)部的節(jié)子、樹脂道等不同結(jié)構(gòu)組分，基本反映了樹干呈現(xiàn)的真實力學(xué)性質(zhì)。因此，本研究以力學(xué)分析為目的時，可將其簡化處理并假設(shè)木材細胞腔等間隙填滿物質(zhì)，木材是連續(xù)均勻的，而且心邊材及樹干不同高度處相應(yīng)的材料力學(xué)性質(zhì)一致。另外，為促使樹干受到不同的彎矩與扭矩作用，在樹干模型頂端上增加一加載桿，建立樹干基本受力模型(圖1)，即Γ形空心圓樹干在自由端受F力作用，高為H、偏心距為e、內(nèi)徑為d、外徑為D，樹干基部為固定端約束。

1.2 有限元靜力學(xué)分析

材料定義時，由于木材是正交各向異性材料，定義材料屬性需建立柱坐標(biāo)系對應(yīng)樹干的力學(xué)性質(zhì)。樹干軸向L用1表示，弦向T用2表示，徑向R用3表示，Ei為彈性模量，uij為泊松比，Gij為剪切彈性模量(i、j為1、2、3，表示不同方向)，則在Wood材料屬性中，E1=EL，E2=ET，E3=ER，u12=uTL，u13=uRL，u23=uRT，G12=GTL，G13=GRL，G23=GRT。以黃山松(Pinustaiwanensis)為研究對象，其生材的彈性常數(shù)見表1，彎曲強度為33.59 MPa，順紋剪切強度7.19 MPa。定義加載桿材料屬性時，因其僅是為便于施加集中載荷，起輔助作用，不需要分析其力學(xué)變化，同時，為消除低彈性模量而產(chǎn)生大的彈性變形對加載效果的影響，定義其彈性模量為具有高彈性模量屬性合金鋼[19]材料的10倍，即令其彈性模量為2 000 GPa，泊松比為0.3。

表1 黃山松生材的彈性常數(shù)

網(wǎng)格劃分時，選擇C3D8R單元作為網(wǎng)格劃分單元，該單元為三維的8結(jié)點六面線性減縮積分單元，相比較完全積分單元，其在每個方向上少用1個積分點，僅在單元中心包含1個積分點，優(yōu)點是在承受彎曲載荷作用下不發(fā)生剪切自鎖現(xiàn)象而引發(fā)計算錯誤和可能造成的計算不收斂，且當(dāng)受扭曲變形時，分析的精度不受影響。結(jié)果查看中，為得到樹干指定高度處截面的最外層應(yīng)力變化情況，可在指定高度截面的應(yīng)力值云圖(圖2)中選擇樹干最外層一圈路徑結(jié)點應(yīng)力值導(dǎo)出，再進行相關(guān)數(shù)據(jù)分析。

為研究不同程度空洞樹干樹木的強度和斷裂破壞方式，本研究建立了高H=9.1 m，外徑D=0.77 m的黃山松樹干基本模型，分析：1)因d/D=α(α=0.1～0.9)變化，即空洞大小變化，扭矩不變，指定樹干高度處截面最外層上單元應(yīng)力的變化情況；2)因e/h=β(β=0.1～0.5)變化，即扭矩發(fā)生變化，樹干空洞大小不變，指定樹干高度處截面最外層上單元應(yīng)力的變化情況；3)比較分析因α與β的變化，指定高度處截面最外徑上，彎曲應(yīng)力與扭轉(zhuǎn)剪應(yīng)力的變化情況。

2 結(jié)果與分析

2.1 樹干破壞系數(shù)與載荷的關(guān)系

破壞系數(shù)指樹干受載荷作用下樹干外層的最大應(yīng)力與其相應(yīng)臨界應(yīng)力強度值的比值。圖3是黃山松內(nèi)外徑比，α=0.5，樹干指定高度1.3 m處截面扭矩與彎矩比，β=0.1時，樹木破壞系數(shù)隨施加載荷力F的變化情況。可知彎曲破壞系數(shù)與扭轉(zhuǎn)破壞系數(shù)均隨著F的增大而增加，即樹干最外層最大彎曲應(yīng)力和扭轉(zhuǎn)剪應(yīng)力均隨著載荷F的增大而增加，且在斷裂破壞前呈線性關(guān)系。

2.2 樹干破壞系數(shù)隨樹干內(nèi)外徑比和扭矩與彎矩比的變化情況

圖4是黃山松自由端受同一載荷作用，樹干指定截面高度1.3 m處樹干破壞系數(shù)隨樹干內(nèi)外徑比(α)和扭矩與彎矩比(β)的變化情況。可以發(fā)現(xiàn)：1)樹干可能發(fā)生的彎曲破壞系數(shù)在內(nèi)外徑比α=0.7時，隨著樹干空洞的繼續(xù)增大其彎曲破壞系數(shù)會急劇增加，在α<0.7時，彎曲破壞系數(shù)有一定的增加，但增加幅度較小。即樹木因空洞造成的樹干內(nèi)外徑比為0.7是樹木可能發(fā)生斷裂的臨界值，這與前人野外統(tǒng)計得到的結(jié)論相一致[5,20]。2)樹干可能發(fā)生的彎曲破壞系數(shù)并不會隨著樹木指定截面高度處承受的扭矩與彎矩比值(β)變化而發(fā)生變化，即不受樹木偏冠的程度大小影響。3)樹干可能發(fā)生的扭轉(zhuǎn)破壞系數(shù)在內(nèi)外徑比為0.7，即α=0.7時，隨著樹干空洞繼續(xù)增大其扭轉(zhuǎn)破壞系數(shù)也急劇增加，在α<0.7時，扭轉(zhuǎn)破壞系數(shù)有一定的增加，但變化幅度較小，與彎曲破壞系數(shù)趨勢一致。4)樹干可能發(fā)生的扭轉(zhuǎn)破壞系數(shù)會隨著樹木指定截面高度處承受的扭矩與彎矩比值(β)的增大而增加。

黃山松樹干指定高度1.3 m處截面受集中載荷破壞系數(shù)隨扭矩與彎矩比(β)和樹干內(nèi)外徑比(α)的變化情況見圖5。

由圖5可見：1)對于樹干內(nèi)外徑比α≤0.5時，當(dāng)樹干指定截面高度處扭矩與彎矩比β≤0.2時，樹干可能發(fā)生的扭轉(zhuǎn)破壞系數(shù)均小于彎曲破壞系數(shù)，即樹干主要發(fā)生彎曲破壞，也就是說，當(dāng)樹木空洞形成的樹干內(nèi)外徑比在0.5以下，且指定高度處樹干扭矩與彎矩比在0.2以下，樹木主要發(fā)生彎曲破壞；而當(dāng)β>0.2時，樹干可能發(fā)生的扭轉(zhuǎn)破壞系數(shù)均大于彎曲破壞系數(shù)，即樹干首先發(fā)生扭轉(zhuǎn)破壞。2)對于樹干內(nèi)外徑比α=0.7時，當(dāng)樹干指定截面高度處扭矩與彎矩比β<0.2，約等于0.16時，樹干可能發(fā)生的扭轉(zhuǎn)破壞系數(shù)小于彎曲破壞系數(shù)，即主要發(fā)生彎曲破壞；而當(dāng)β>0.16時，樹干發(fā)生的扭轉(zhuǎn)破壞系數(shù)大于彎曲破壞系數(shù)，即樹干首先發(fā)生扭轉(zhuǎn)破壞。3)對于樹干內(nèi)外徑比α=0.9時，當(dāng)樹干指定截面高度處扭矩與彎矩比β<0.15，約等于0.13時，樹干可能發(fā)生的扭轉(zhuǎn)破壞系數(shù)小于彎曲破壞系數(shù)，即主要發(fā)生彎曲破壞；而當(dāng)β>0.13時，樹干發(fā)生的扭轉(zhuǎn)破壞系數(shù)大于彎曲破壞系數(shù)，即樹干首先發(fā)生扭轉(zhuǎn)破壞。

由上述分析可知，隨著樹干空洞造成的內(nèi)外徑比增加，樹木首先發(fā)生扭轉(zhuǎn)破壞的偏冠程度會逐漸降低，再次證明了當(dāng)樹木生長到一定的偏冠程度，會因樹干空洞的增加，樹木先發(fā)生扭轉(zhuǎn)破壞后，造成樹干整體剛度的降低再引發(fā)樹干折斷。

3 結(jié)論與討論

3.1 結(jié)論

利用有限元方法分析了樹木指定截面高度處不同內(nèi)外徑比和不同扭矩與彎矩比下樹干最外層的應(yīng)力變化情況，探討樹木破壞的方式。結(jié)果表明，樹木發(fā)生彎曲和扭轉(zhuǎn)破壞系數(shù)隨載荷增加而增大，樹干最外層的最大彎曲和扭轉(zhuǎn)剪應(yīng)力隨載荷增加呈線性增加；同時，當(dāng)樹干的空洞發(fā)展到樹干內(nèi)外徑比α=0.7時，隨著樹干空洞的繼續(xù)增大其彎曲和扭轉(zhuǎn)破壞系數(shù)會急劇增加，在α<0.7時，彎曲和扭轉(zhuǎn)破壞系數(shù)有一定的增加，但增加幅度較小。這為林業(yè)工作者對含有腐朽空洞的古樹名木、高大行道樹及庭蔭樹采取支撐拉纖等保護措施的介入時機提供了參考和依據(jù)。此外，研究還發(fā)現(xiàn)，隨著樹干空洞造成內(nèi)外徑比的增加，樹木首先發(fā)生扭轉(zhuǎn)破壞的偏冠程度會逐漸降低，這與前人的研究結(jié)果相一致[5-6]，再次證明了樹木達到一定的偏冠程度后，會因樹干空洞的增加，樹木先發(fā)生扭轉(zhuǎn)破壞，造成樹干整體剛度的降低再引發(fā)樹干折斷。

3.2 討論

利用有限元對樹木強度分析主要集中于樹木的動力學(xué)研究，D.Sellieretal[21]基于枝條的動力學(xué)效證明了樹木達到一定的偏冠程度后，會因樹干空洞的增加，樹木先發(fā)生扭轉(zhuǎn)破壞，造成樹干整體剛度的降低再引發(fā)樹干折斷。應(yīng)建立樹木模型并預(yù)測樹木的風(fēng)激響應(yīng)，J.R.Mooreetal[22]建立了花旗松有限元模型分析樹冠質(zhì)量對其自振頻率的影響，艾曉秋等[14]建立香樟行道樹的有限元模型分析城市行道樹動力學(xué)特性與風(fēng)致破壞影響，而基于樹木樹干作為各向異性材料，從靜力學(xué)方向利用有限元分析樹干破壞方式則鮮有研究。

本研究基于測定的黃山松力學(xué)性質(zhì)，就不同形狀的樹干以及樹洞在樹干內(nèi)部不同位置利用有限元建立相應(yīng)的數(shù)值模擬模型進行了數(shù)據(jù)分析。圖6為圓形樹干偏心圓形樹洞樹高1.3 m處彎曲應(yīng)力和扭轉(zhuǎn)剪應(yīng)力云圖，顯示最薄處、最外層的彎曲應(yīng)力和扭轉(zhuǎn)剪應(yīng)力最大，當(dāng)外載荷增加，即可提取出最薄處、最外層的節(jié)點應(yīng)力值與材料的彎曲強度和順紋剪切強度進行比較，判定可能發(fā)生的破壞方式。圖7為橢圓樹干圓形樹洞樹高1.3 m處彎曲應(yīng)力和扭轉(zhuǎn)剪應(yīng)力云圖，圖8為橢圓樹干橢圓樹洞樹高1.3 m處彎曲應(yīng)力和扭轉(zhuǎn)剪應(yīng)力云圖，圖9為橢圓樹干偏心橢圓樹洞樹高1.3 m處彎曲應(yīng)力和扭轉(zhuǎn)剪應(yīng)力云圖，從應(yīng)力云圖中均能直觀地發(fā)現(xiàn)樹干承受的最大彎曲應(yīng)力和扭轉(zhuǎn)剪應(yīng)力所處位置和大小，幫助分析樹干可能發(fā)生的破壞方式。

自然環(huán)境中，樹干形狀并非規(guī)則的圓形、橢圓形，而且樹干腐朽空洞也呈現(xiàn)不同形狀，基于純理論計算外載荷下樹干指定位置處的彎曲和扭轉(zhuǎn)剪應(yīng)力將非常復(fù)雜，有限元方法的出現(xiàn)能很好地幫助解決計算難題。同時，隨著科技的發(fā)展，各種樹木應(yīng)力波掃描儀[23-24]可以對樹體內(nèi)部的腐朽空洞大小、形狀進行確認，同時高精度激光三維掃描技術(shù)[25-26]可以建立與真實樹干仿真性極高的樹干三維模型，下一步，結(jié)合上述技術(shù)建立個體真實情況的腐朽空洞樹干有限元模型并進行力學(xué)強度分析，將對古樹名木、大樹、行道樹安全性評估具有更加積極的現(xiàn)實意義。