含裂紋損傷圓弧曲梁彈性屈曲的有限元網格自適應分析

王永亮

(1. 中國礦業大學(北京)力學與建筑工程學院，北京 100083；2. 中國礦業大學(北京)煤炭資源與安全開采國家重點實驗室，北京 100083)

桿系結構的損傷問題廣泛存在于工程實際中，工程桿系結構大多帶裂紋工作[1-2]，裂紋損傷的存在會改變整個結構的力學性能，影響結構的安全性和適用性[3]。研究含多裂紋損傷梁構件的動力特性、準確預測屈曲承載力可以有效地保障結構在全生命周期內安全使用[4]。曲線形梁構件由于幾何形態復雜，容易誘發彈性屈曲失穩[5]，精確評估各類曲線梁線型、不同曲梁夾角下深梁、淺梁的屈曲荷載成為結構災害分析的重要依據。曲梁中裂紋損傷的存在增加準確預測屈曲失穩承載能力的難度，理論模型、解析方法等往往難以有效分析[6-7]。準確預測不同裂紋損傷位置、大小、數目工況下屈曲荷載承載力以及分析裂紋損傷對屈曲失穩的影響機理[8]，成為理論研究和工程實踐的需求。

有限元法被發展和應用于求解含裂紋損傷曲梁的彈性屈曲荷載和屈曲模態[9-11]，但解答精度依賴于網格劃分質量，解答因網格劃分難免引入誤差[12]。有限元網格自適應分析方法可有效地優化網格分布，在直線梁彈性屈曲[13]、板殼振動[14]、含損傷梁振動[15]、巖體變形和斷裂[16]等問題求解中展示出很好的求解效力。本文將建立圓弧形曲梁裂紋的截面損傷缺陷比擬方案，進行裂紋大小(深度)、位置、數目的模擬，引入變截面Euler-Bernoulli梁的h型有限元網格自適應分析方法[15]，求解含裂紋損傷圓弧曲梁彈性屈曲問題，得到優化的網格和滿足預設誤差限的高精度屈曲荷載和屈曲模態。文中給出求解多種圓弧曲梁彈性屈曲數值算例，對網格自適應劃分以及彈性屈曲解答的收斂性進行了討論，對曲梁夾角、損傷位置、數目、大小等因素影響彈性屈曲荷載和屈曲模態進行了分析，檢驗了網格自適應劃分的有效性。

1 圓弧曲梁彈性屈曲

圖 1 含裂紋損傷曲梁坐標系和符號Fig. 1 Coordinate systems and symbols of cracked curved beam

圖 2 含裂紋損傷圓弧曲梁截面損傷和加載示意圖Fig. 2 Diagram of cross-section damage defect and loading for circularly curved beam with crack damage

表 1 圓弧曲梁根據徑厚比和夾角分類Table 1 Categories of circularly curved beams according to ratio of radius and thickness and subtended angle

本文研究圓弧曲梁彈性屈曲的微分控制方程為[13]：

2 裂紋損傷表征方法

本文采用Euler-Bernoulli梁理論模型研究梁高相對較小薄曲梁，若裂紋損傷發生在梁單一側，則損傷截面中性軸與無損傷截面偏移量較小；若裂紋損傷均勻發生在梁上下兩側，則損傷截面中性軸與無損傷截面完全重合。因此，本文研究的曲梁無裂紋和有裂紋橫截面僅高度不同，二者中性軸考慮為處于相同位置。曲梁中的裂紋損傷，使得梁截面產生弱化、梁的抗彎剛度衰減。本研究采用裂紋的截面損傷缺陷比擬方法[15]，裂紋處的截面剛度為：

3 網格自適應細分加密

有限元計算存在相比當前網格解答具有更高收斂階的超收斂點[19]，利用超收斂點結合單元拼片、高階形函數插值技術，可以提高當前有限元解的精度，得到全域的超收斂解[15,20-21]。本文對于圓弧曲梁的彈性屈曲問題，求得當前網格下屈曲模態(位移)的有限元解后，利用有限元后處理超收斂拼片恢復方法，得到屈曲模態的超收斂解：

利用屈曲模態誤差估計，網格可以進行優化處理來降低和控制屈曲模態的誤差，達到預設的解答精度。本文方法對每個有限元單元e上的振型誤差進行判斷，如果誤差控制式(12)不滿足，則表明該單元上屈曲模態解答的誤差過大，需要通過進行網格優化處理，本文采用單元均勻細分加密的h型網格自適應方式來增加模型自由度、降低單元上解答的誤差[15]。當前單元細分生成的新單元長度與目前誤差和單元階次相關，即利用當前誤差可以估計新單元長度：

4 數值算例

本文方法已經編制相應的Fortran 90語言程序代碼，程序開發實施基于Microsoft Visual Studio和Intel Visual Fortran編程軟件平臺。本節給出求解具有代表性的多種圓弧曲梁彈性屈曲數值算例，對網格自適應劃分以及彈性屈曲解答的收斂性進行了討論，對曲梁夾角、損傷位置、數目、大小等因素影響彈性屈曲荷載和屈曲模態進行了分析，檢驗了網格自適應劃分的有效性。本節所有算例均采用3次元，初始網格采用2個單元，給定的初始誤差限為。

4.1 圓弧曲梁屈曲分析的網格優化

考慮一兩端簡支的圓弧曲梁，梁的幾何和物理參數如下：

本文方法對該曲梁彈性屈曲進行求解，得到屈曲荷載和屈曲模態的解答。文獻[23]采用層狀組合梁模型及有限元方法、文獻[7]采用理論模型解析法分別對該梁彈性屈曲進行分析，得到無裂紋損傷梁屈曲荷載值。為進行對比分析，將上述各方法求得的屈曲荷載值均列于表2。同時，為討論網格數目對解答收斂的影響，本研究采用常規有限元法在4個、5個、6個單元(稀疏均勻網格)以及25個、50個、100個單元(密集均勻網格)上分別進行求解，得到屈曲荷載值；可知隨著單元增多，解答趨于穩定，在接近100個單元時得到收斂解答。使用本文自適應有限元方法進行求解，得到僅為16個單元的優化網格和在此網格下的收斂解答，該屈曲荷載解答(64.406 kN/m)與理論模型解析解(64.966 kN/m)具有很好的吻合度。需要說明的是，組合梁模型使用4個、5個、6個單元得到的解答與解析解相差較大，為梁模型的層狀建模引起。

表 2 單元數目與屈曲荷載結果收斂性Table 2 Convergence for number of elements and buckling loads results

圖3給出本文方法利用自適應網格求解得到的屈曲模態解答，為方便直觀顯示和分析，屈曲模態結果均進行歸一化處理(令最大模態值為1)。由于該梁的物理性質和幾何形態具有左右對稱性，得到圖3(a)所示的屈曲模態亦為左右對稱形式；同時，水平坐標軸上給出了本文的自適應網格(單元端節點)分布，可知單元分布適應模態變化、同樣具有左右對稱性，且在兩端邊界部分使用了相對細密的網格。算法自動優化出非均勻網格，在屈曲模態變化平緩區域使用稀疏網格、在屈曲模態變化劇烈處采用了相對細密的網格，避免了全域使用一致細密網格的冗余性。為了檢驗屈曲模態解答的精確性，本研究使用常規有限元法2500個單元(高密集均勻網格)求解得到高精度屈曲模態作為解析解；該高精度屈曲模態形態與圖3(a)相同，這里不再給出。圖3(b)所示為使用本文自適應方法的屈曲模態與高精度屈曲模態的差值曲線分布，該差值曲線最大值2.17×10-5小于預設誤差限10-4，驗證了本文方法求解解答的精確性。

常規有限元采用4個、5個、6個、25個、50個、100個單元求解得到的屈曲模態解答與高精度屈曲模態解答的差值分布分別如圖4所示，并在圖中水平坐標軸上給出各網格分布。各單元下差值在全域上的最大值分別為1.23×10-1、1.23×10-1、4.12×10-3、7.70×10-3、1.62×10-3、2.63×10-6，可見隨著網格的均勻加密，屈曲模態解答的誤差呈現逐漸降低的趨勢；直到提供足夠多的單元(100個單元)時，才接近獲得滿足預設誤差限10-4的解答。本文自適應方法僅采用16個單元的非均勻分布網格，即可避免高密度均勻分布網格的單元冗余性，提高了計算效率。

圖 3 自適應網格下屈曲模態Fig. 3 Buckling mode on adaptive refinement mesh

4.2 變曲梁夾角下的屈曲荷載

為檢驗本文方法求解不同幾何形式圓弧曲梁彈性屈曲的適用性，本研究對不同圓弧夾角曲梁(淺梁、深度梁)進行分析?？紤]一兩端簡支的圓弧曲梁，梁的幾何和物理參數如下：

4.3 彈性屈曲的裂紋損傷位置影響

圖 4 加密網格下屈曲模態收斂情況Fig. 4 Convergence of buckling modes on refined meshes

表 3 兩端簡支曲梁不同夾角下屈曲荷載值Table 3 Buckling loads of of curved beam with hinged-hinged supports under different subtended angles

表 4 不同裂紋損傷位置下曲梁屈曲荷載值Table 4 Buckling loads of curved beam with crack damage at different locations

圖 5 不同裂紋損傷位置下圓弧曲梁彈性屈曲Fig. 5 Elastic buckling of circularly curved beam with different locations of crack damage

4.4 彈性屈曲的裂紋損傷大小影響

本節分析圓弧曲梁裂紋損傷大小(深度)對彈性屈曲的影響，固定裂紋損傷位于跨中(θc=π/2)，分析裂紋損傷大小hc/h值為0.1、0.2、0.3、0.4、0.5時的彈性屈曲。使用本文方法分別進行求解，得到彈性屈曲荷載解答，結果列于表5?？芍?，隨著裂紋損傷hc增大(hc/h=0.1 →0.5 )，qˉc逐漸降低。

表5 不同裂紋損傷大小下曲梁屈曲荷載值Table5 Buckling loadsof curved beam with crack under different magnitudes

為更直觀顯示和分析變化趨勢，利用上述不同裂紋損傷大小下屈曲荷載結果繪制圖6(a)所示變化曲線?？梢娗奢d隨裂紋損傷加深，出現準線性的降低趨勢；在裂紋擴展到梁高的一半(hc/h=0.5)時，屈曲荷載承載力降低約10%。圖6(b)所示為不同裂紋損傷大小時各屈曲模態與無損傷屈曲模態的差值曲線，可以看出裂紋損傷大小對屈曲模態變化幅度具有重要影響，裂紋損傷程度越大，越容易誘發屈曲模態的大幅度變化。同時，水平坐標軸上給出本文方法求解hc/h=0.5時的自適應最終非均勻網格。

4.5 多裂紋損傷下屈曲模態和網格劃分

為分析多裂紋損傷對曲梁彈性屈曲的影響，采用圖7 所示的兩端簡支圓弧曲梁，該曲梁夾角為θ=π，其余基本幾何和物理參數同式(12)。該曲梁包含3條裂紋損傷，考慮I(裂紋損傷沿曲梁均勻分布)、II(裂紋損傷集中于曲梁端部)、III(裂紋損傷集中于曲梁中部)等3種工況，各工況中裂紋分布角坐標值如表6所示，裂紋分布如圖7所示。

使用本文方法計算該曲梁多裂紋工況的結果如表6所示，工況III的屈曲荷載值最小，即多裂紋損傷越集中于跨中，越容易誘發彈性失穩，這與上文單一裂紋損傷位于跨中易于導致彈性失穩的結論一致。

圖8所示為多裂紋工況各屈曲模態與無損傷屈曲模態的差值曲線，可以看出各裂紋損傷所在局部區域對屈曲模態變化有重要影響，多裂紋損傷缺陷導致屈曲模態劇烈變化。同時，本文算法在各裂紋損傷附近區域使用了相對密集的網格，用于適應裂紋損傷引起屈曲模態的變化，形成優化的非均勻網格、確保解答的可靠性。

表6 多裂紋損傷不同位置下曲梁屈曲荷載值Table 6 Buckling loadsof curved beam with multiple cracks at different locations

圖6 不同裂紋損傷大小下圓弧曲梁彈性屈曲Fig.6 Elastic buckling of circularly curved beam with different magnitudes of crack damage

圖7 各工況多裂紋分布示意圖Fig.7 Diagram of casesfor multiple cracksdistributions

圖8 多裂紋損傷下網格分布和屈曲模態差Fig.8 Mesh distribution and buckling modes of circularly curved beam with multiple cracks

5 結論

本文建立了圓弧形曲梁裂紋的截面損傷缺陷比擬方案，實現裂紋的大小(深度)、位置、數目的模擬；針對含裂紋損傷圓弧曲梁彈性屈曲，引入有限元網格自適應分析方法，得到了優化的網格和滿足預設誤差限的高精度屈曲荷載和屈曲模態解答。經數值算例檢驗，本文方法對淺梁、深度梁等各類變化幾何形式圓弧曲梁的彈性屈曲求解具有良好適用性，解答與解析解具有較高吻合度。研究發現：隨著裂紋損傷增大、接近跨中，均會不同程度降低屈曲荷載，越容易誘發彈性失穩。裂紋損傷將誘發屈曲模態變化，本文自適應算法可劃分出非均勻網格，在裂紋附近區域使用了相對密集的網格來適應裂紋損傷引起屈曲模態的變化。