屈曲問題有限元模擬的隨機缺陷法

閆 棣， 蘇 祺， 李四平

(上海交通大學 船舶海洋與建筑工程學院， 上海 200240)

固體力學研究領域中，工程結構的非線性屈曲問題一直是非常活躍的話題[1].在使用通用有限元軟件進行屈曲模擬的過程中，一些文獻[2-3]直接使用軟件中提供的線性屈曲分析模塊進行屈曲臨界載荷的預測與屈曲模態的判斷.線性屈曲分析方法相當于彈性屈曲分析，然而材料的非線性、結構形狀的初始缺陷以及幾何大變形使得實際結構并不在理論彈性范圍內屈曲.因此，特征值屈曲分析所得的解為非保守解，一般不能用于實際工程的分析[4-5].除此之外，線性屈曲分析方法對材料非線性等問題無法解決，有一定的局限性.一些文獻[6-7]使用線性屈曲分析中所得的一階模態或前幾階模態進行缺陷的引入，而在引入缺陷時所用的比例因子參數沒有定論.這種方法通過人為調控比例因子的大小可能得到不同的結果，計算結果仍需要通過進一步的試驗比對來調整參數，而且在線性屈曲方法無法計算的材料非線性等問題上[8]，繼續使用此方法進行屈曲的分析，所得結論是不可信的.

本文提出一種更為普適的模擬方法，即通過隨機的方法對理想結構模型引入初始缺陷進行屈曲的分析.這種方法不僅能夠得到上述文獻中所用方法計算得到的多階屈曲模態和其對應臨界力，而且將特征值問題轉化為幾何非線性靜力變形問題進行計算，克服了線性屈曲方法對材料非線性等問題的局限性，獲得了更多中間加載過程的信息.此外，以彈性地基上梁的屈曲問題和橫向纖維作用下穿透型矩形脫層屈曲問題為例，計算討論隨機缺陷法的普適性和以上文獻中所使用方法的局限性.同時，通過對比不同大小缺陷對計算結果的影響，提出合理的隨機缺陷引入方法.

1 隨機缺陷法

1.1 理論分析

初始缺陷是指與理論計算模型相比，實際構件在尚未受到載荷作用前就已經存在的各種缺陷，其中包含力學初始缺陷以及幾何初始缺陷.力學初始缺陷主要是由各種原因引起的初始殘余應力；而幾何初始缺陷則包含整體初彎曲，初偏心，局部初始缺陷，等等.

考慮如圖1所示的兩端鉸支的細長桿，當它處于微彎曲平衡的臨界狀態時，若桿內應力不超過材料的比例極限，則壓桿任一截面m-n處的彎曲變形與截面彎矩的關系為

(1)

式中：E為桿件彈性模量；I為桿件截面慣性矩；M(x)為桿件在距離下支座x處的由外力F引起的彎矩.

圖1 兩端鉸支細長桿Fig.1 Schematic of hinged-hinged slender beam

考慮模型有初始缺陷，由于初始缺陷是一個隨機的無窮小量，忽略高階小量，做疊加處理得

(2)

式中：M0(x)為初始缺陷引起的初彎矩，是一個隨機的非連續小量.由于M(x)=-Fy，所以式(2)可寫為

(3)

(4)

解二階非線性常系數微分方程，其通解及特解分別為

Y=asinkx+bcoskx

(5)

(6)

式中：D為微分算子.代入桿端邊界條件(x=0時，y=0，M0(0)=0；x=l時，y=0，M0(l)=0)，可以得到式(4)的解為

(7)

(8)

即當M0(x)趨近于0時，式(7)同式(1)的解等價，表示當加入的隨機缺陷足夠小時并不會影響模型整體的屈曲結果.

然而理想不含初始缺陷的構件在外力作用下并不會發生分支點屈曲，分支點屈曲發生的原因也正是由于實際構件中包含一定的初始缺陷.在通用有限元軟件線性屈曲分析中，將屈曲問題作為特征值問題進行求解.而隨機缺陷法，通過隨機的方法在理想構件中引入缺陷將屈曲作為幾何非線性靜力計算問題進行求解，這更符合實際情況.

在對模型施加缺陷時，隨機缺陷法強調隨機性，充分模擬實際構件中缺陷發生的隨機性.缺陷可以是構件的幾何缺陷，也可以是力學初始缺陷，引入這些小缺陷的目的就是觸發構件屈曲.缺陷大小和隨機性是計算結果精度的主要影響因素.

1.2 有限元實現方法的探討

本文算例中采用隨機選取模型網格節點進行初始偏移的方法隨機地添加缺陷，根據計算結果可以得知偏移的大小取梁截面徑向尺寸的千分之一或更小為宜.這是一種簡單方便的幾何缺陷引入方法.除此之外，本文也嘗試通過程序隨機選取結構中某一微小單元并賦予較小的模量來模擬材料的空心幾何缺陷.力學缺陷上，可以使用虛擬的觸發載荷來啟動不穩定性，也可以通過不均勻的溫度場致使屈曲的發生.

在保證隨機性的前提下，引入初始缺陷的方法有很多，其根本的目的是為了觸發屈曲.初始缺陷越小，有限元模型越接近理想模型，計算結果與理論結果也越接近，但此時為觸發屈曲，模型的計算量也會相應增加.所以有限元計算中觸發屈曲時，初始缺陷的大小、模型的計算量以及計算結果的精度需要有一個平衡的過程.

2 兩端鉸支細長桿多階屈曲模擬

設有長l=100 cm，截面面積A=4 cm2，慣性矩Ix=Iy=1.33 cm4的兩端鉸支細長桿件，長細比λ=173.2，如圖2所示.材料的彈性模量E=200 GPa.

圖2 兩端鉸支壓桿Fig.2 Schematic of hinged-hinged compression bar

使用通用有限元軟件進行建模，采用梁單元進行模擬，并將梁沿軸向劃分為100個單元.通過線性屈曲方法得到前5階屈曲模態及臨界載荷，結果以第1階和第2階模態為例，如圖3所示，其中u2表示桿件徑向位移.參照Euler公式

(9)

式中：Fcr為臨界力；n為屈曲模態數.線性屈曲方法模擬結果與Euler公式理論結果(Fcr1=26.32 kN，Fcr2=105.28 kN)較為吻合.

圖3 線性屈曲方法所得壓桿屈曲前兩階模態Fig.3 The first two order modes of compression bar by linear buckling method

圖4 隨機缺陷法計算壓桿屈曲第1階模態Fig.4 The first mode of compression bar by stochastic defect method

使用隨機缺陷法進行模擬.通過編寫程序操作計算文件，隨機選取桿件上節點施加10-5R0(R0為橫截面徑向尺寸)的橫向偏移，加載位移載荷，使其從0開始緩慢增長，其第1階屈曲模態如圖4所示.

為保證隨機性不影響最終結果，進行了多次計算，每次選取的偏移節點均由程序隨機生成，所得結果相同.在使用節點偏移法的同時，嘗試使用空心幾何缺陷，通過實體單元進行建模，隨機選擇模型中某一位置劃分出微小的單元塊，并賦予較低的材料模量來模擬空心，載荷施加不變，所得到的結果與前述方法相同.

在第1階模態基礎上，第2次計算時不改變隨機缺陷位置，使位移載荷在加載開始時跳過第1階的臨界位移載荷，而后緩慢增長，尋找第2階模態的臨界載荷，計算結果如圖5所示.

圖5 隨機缺陷法計算壓桿屈曲第2階模態Fig.5 The second mode of compression bar by stochastic defect method

通過以上算例可以看出，隨機缺陷法不僅可以計算模型屈曲的最低階模態，而且可以通過加載參數的改變來計算高階模態，計算結果與理論結果較為吻合.

3 彈性地基上梁的屈曲

3.1 嵌入于地基中的梁

將上文中的桿件嵌入于彈性地基中.地基使用Winkler地基模型進行模擬，彈簧勁度系數β=20 MN/m，模型如圖6所示.

圖6 Winkler彈性地基梁Fig.6 Schematic of Winkler foundation model

使用通用有限元軟件進行建模，采用梁單元進行模擬，并將梁沿軸向劃分為100個單元，在除兩端外的每個網格節點處添加線性彈簧單元[9-10].

使用線性屈曲方法進行計算，在梁的一端施加單位位移載荷，所得結果如圖7所示.計算特征值為7.655 86×10-2，所以臨界應變εcr1=7.655 86×10-2.

3.2 擱置于地基表面的梁

在以上模型的基礎之上，加入非線性彈簧，即模擬梁未嵌入地基，而是放置在地基表面上的情況.此時，梁在向上彎曲時，與地基表面脫離；而向下彎曲時，則受到來自地基的壓力.模型相當于兩端鉸支的壓桿穩定問題.在建模過程中彈簧壓縮時的勁度系數β=20 MN/m，而拉伸時勁度系數趨近于0.

使用線性屈曲方法進行計算，發現計算結果與使用線性彈簧時的結果相同.線性屈曲方法會選取彈簧初始剛度或者上一分析步結束時的剛度作為線性剛度進行計算，而忽略后續彈簧參數的變化[11]，即不能計算非線性彈簧或者非線性材料等問題.而在此基礎之上通過影響因子的方式將線性屈曲方法的計算結果作為缺陷引入到模型中再進行靜力計算，顯然也有問題，其結果的有效性有待驗證.同時可以認為，對于材料非線性等問題，使用此方法仍然會失效.

使用隨機缺陷法進行計算，同樣地隨機選取梁上節點引入10-5R0的橫向偏移，其計算結果如圖9與10所示.隨機缺陷法計算所得的臨界載荷Fcr=26.31 kN，與Euler公式的計算結果吻合， 說明隨機缺陷法對非線性彈簧問題同樣適用，相比于線性屈曲方法有較強的普適性.

考慮不同大小初始缺陷對最終臨界載荷Fcr計算結果的影響.對模型分別引入不同大小初始缺陷，所得計算結果如表1所示.其中:δ=R/R0，R為初始缺陷.

圖9 擱置地基表面梁受壓端軸向位移與整體支反力關系Fig.9 Relation of axial displacement of foundation surface beam compression end and integral bearing reaction

圖10 擱置地基表面梁中點處橫向位移與整體支反力關系Fig.10 Relation of transverse displacement at midpoint of foundation surface beam and integral bearing reaction

表1 桿初始缺陷與臨界載荷關系Tab.1 Relation table of stochastic defect and critical load

通過對比分析可以發現，當δ=10-3或更小時，所得計算結果的精度較高.

4 橫向纖維搭橋下穿透型矩形脫層屈曲

設有長2a=200 mm，寬2b=100 mm的穿透型矩形脫層如圖11所示，受到兩側均布載荷Px的作用.脫層的厚度h=1 mm，在橫向纖維的影響下，受到均布載荷q(x)的作用抑制屈曲.材料的彈性模量E=2.3 GPa，泊松比μ=0.3.文獻[12]中也進行了類似的計算，其所使用的方法是首先使用線性屈曲模塊，然后將1階屈曲模態位移乘以比例因子來進行初始缺陷的引入，最終計算結果與理論結果沒有很好的對應關系.

圖11 橫向纖維搭橋下矩形脫層Fig.11 Schematic of rectangular delamination bridged by fibers

根據文獻[13]所述理論進行分析，

q(x)=βω(x),β>0

(10)

(11)

(12)

式中：q(x)為纖維約束力；β為纖維假設為彈簧時的勁度系數；ω(x)為屈曲撓度；D為抗彎剛度；a0為屈曲時特征長度.

根據表2中計算結果與理論值的對比可以發現，使用隨機缺陷法所得到的計算結果有較高的可信度.這種方法通過隨機地改變初始缺陷，可充分模擬實際構件的缺陷形式，降低人為因素對計算結果的干擾，對工程應用中分支點屈曲問題的模擬有較高的參考價值.

圖12 屈曲模態圖Fig.12 Schematic of buckling mode

5 結論

通過使用線性屈曲方法和隨機缺陷法分別對壓桿穩定、彈性地基上梁的屈曲以及橫向纖維搭橋下穿透型矩形脫層屈曲問題進行了有限元模擬.得到以下幾點結論：

(1) 隨機缺陷法將屈曲作為幾何非線性靜力計算問題進行有限元模擬，改進了線性屈曲計算中將屈曲作為特征值問題進行求解的不足，與實際情況更為相符.

(2) 討論了隨機缺陷法的理論依據，并給出了通用有限元軟件使用初始隨機缺陷法的幾種思路.

(3) 隨機缺陷法具有更強的普適性，解決了線性屈曲分析無法計算非線性彈簧等問題.同時，隨機缺陷法同樣可以計算模型高階的屈曲模態.

(4) 通過對比不同大小初始隨機缺陷的計算結果，發現初始隨機缺陷取構件截面徑向尺寸千分之一或更小的值時，有限元計算結果與實際情況較為吻合.