結構極限強度分析的三種有限元解法研究

彭大煒, 張世聯

(上海交通大學,上海200030)

結構極限強度分析的三種有限元解法研究

彭大煒, 張世聯

(上海交通大學,上海200030)

介紹了弧長法、阻尼因子法和準靜態法等3種計算結構極限強度的非線性有限元解法及其求解思路。通過Reckling No.23模型數值計算,對上述3種解法各自的計算過程、求解特點和要點進行了比較和歸納。為更準確、高效地運用非線性有限元程序計算分析結構極限強度提出了合理的建議。

弧長法;阻尼因子法;準靜態法;極限承載能力

0 前言

近年來,結構設計的思路已經逐步從傳統的許用應力設計向極限狀態設計發展[1]。研究結構的極限強度水平,確定其與設計載荷水平之間的確切裕度,已成為結構理性設計的核心和基礎。

目前,船體結構極限承載能力研究的方法主要有[2]:實船事故調查和模型試驗法、直接法(解析法)、以Smith法為代表的逐步崩潰法和非線性有限元法。隨著計算機技術的不斷發展,非線性有限元法現已成為計算和評估結構極限承載能力最理想的方法,其分析問題的主要流程如圖1所示。

在非線性有限元法中,有3種不同的分析結構極限強度的解法。它們分別是弧長法、阻尼因子法和準靜態法。本文基于ABAQUS大型通用非線性有限元程序,以Reckling-No.23試驗模型為對象,分別采用上述3種解法進行極限強度的數值分析,并研究不同解法的求解過程,對不同解法的應用特點進行了歸納和總結,且還針對如何更準確、高效地運用非線性有限元法計算分析結構極限強度的問題提出合理的建議。

1 理論介紹

從求解本質而言,弧長法和阻尼因子法是以結構非線性靜態平衡方程式(1)的求解(牛頓拉-普森迭代)為基礎。而準靜態方法是以結構非線性運動方程式(2)的顯式求解(中心差分法)為基礎。

式中:{P}為載荷列陣;{I}為內力列陣;[M]為質量矩陣;{¨u}為加速度列陣。

以下分別簡要介紹這3種解法的基本思路。

1.1 弧長法

弧長法的基本思路是通過設置一個參數(弧長l)來控制平衡方程的增量迭代和收斂。將(1)式寫成增量格式:

式中:[KT]為切線剛度矩陣;{Δu}為位移增量;{ΔP}為載荷增量;{R}為殘差力。

在弧長法中,第i步迭代的載荷增量{ΔP}i,由載荷增量因子Δ λi和參考載荷{Pref}來控制,即

將(4)式代入(3)式即得到弧長法的第i步迭代的增量形格式:

圖2 弧長法示意圖

弧長法求解時,是以前一步增量計算得到的平衡點為圓心,以弧長增量Δli為半徑,通過牛頓拉-普森迭代搜索下一個平衡點,如圖2所示。每一步的弧長增量Δli、載荷增量因子Δ λi和位移增量{Δu}i通過下面的約束方程來控制:

通過迭代,直至殘差力在容差{R}i范圍內。當第i步迭代完成時,有

由于在弧長增量Δli中同時包含了載荷增量{ΔPi}和位移增量{Δui}的信息,運用弧長法能夠全程跟蹤結構在“加載”(Δ λi＞0)和“卸載”(Δ λi＜0)過程中的載荷-位移平衡路徑。但在實際運用時,存在跟蹤失敗[3]和收斂性問題。目前弧長法還在不斷地研究改進。

1.2 阻尼因子法

阻尼因子法仍采用傳統的載荷步長增量來進行非線性迭代求解。其思路是通過在平衡方程式(1)中引入阻尼力項來求解結構不穩定崩潰的問題,求解的平衡方程為

其中{Fv}為阻尼力列陣,取決于求解時結構變形的速度,由廣義節點速度{v}來表征:

式中:c為阻尼系數[4];M＊為人工質量矩陣;廣義節點速度{v}={Δu}/Δt。

加載的初始階段,結構尚處于穩定狀態,此時廣義節點速度{v}很小,故阻尼力項{Fv}對平衡方程式(10)幾乎沒有影響,可以忽略。隨著載荷的不斷增加,結構趨向于不穩定,當外載{P}已經不能完全由結構內力{I}來平衡時,結構達到極限狀態,相當部分的應變能將釋放轉化為動能,廣義節點速率{v}迅速增大。此時,阻尼力項{Fv}通過做功消耗釋放的應變能,在平衡方程式(10)中起到維持求解系統的“平衡”作用。

由于采用傳統的載荷步長來進行非線性迭代,阻尼因子法無法繼續有效跟蹤結構在“卸載”過程中的載荷-位移路徑,取而代之的是一條幾乎水平(斜率為0)的直線。通過考察阻尼力項{Fv}為維持系統“平衡”所做的功所占結構應變能百分比的歷史變化曲線,就能夠確定結構的極限承載能力。相比弧長法、阻尼因子法的數值收斂性要好一些,能夠解決的問題也更廣泛。

1.3 準靜態法

與弧長法和阻尼因子法的靜態求解方法不同,準靜態法從本質上講是一個結構動態求解的過程。在求解時,對式(2)采用中心差分法進行顯示的時間積分,由一個增量步的動力學條件計算下一個增量步的動力學條件,直至求解時間結束。準靜態法的基本思路就是用慢速加載的動態分析來模擬靜態問題,所以,求解關鍵是要設置合適加載速率。加載速率過快會導致求解結果的局部性(劇烈的結構局部變形),使計算結果偏離“準靜態”的要求;而過慢的加載速率意味著較長的加載時間,從而使計算時間大幅度增加。故分析時,通常會取從快到慢多個加載速率進行分析比較,以選定一個合適的加載速率。判斷加載速率是否合適的一個重要標準就是分析過程中結構模型的動能與其應變能之比,一般準靜態的要求是小于5%。

由于中心差分法是條件穩定的算法[5],在分析時,時間步長Δt必須小于穩定性限制Δtstable才能保證求解的穩定性:

式中:ωmax為結構最高階固有頻率;Δte為結構模型中最小尺寸的單元的穩定時間步長。Δte與單元的特征尺度Le、彈性模量E和材料密度ρ有關,是穩定性限制Δtstable的上界[6]:

在計算中,由于結構最高階頻率ωmax不易求得,故時間步長Δt就取為Δte。

相比弧長法和阻尼因子法,準靜態方法最大的優勢在于采用中心差分法進行顯式時間積分不存在收斂性的問題。因此,準靜態法能夠很好地求解更復雜結構崩潰問題,如結構的自接觸和材料的失效問題。在求解極限狀態問題時,時間步長Δt往往較小,準靜態“緩慢”加載可能導致求解的機時很長,這時可采用質量放大等方法進行調整。

2 數值計算分析

2.1 Reckling No.23模型試驗

Reckling(1979)[7]采用箱型剖面模擬實船剖面進行了系列總縱極限承載能力試驗,以研究船體梁在極限狀態下各種崩潰模式和剖面應力分布。本文選取加強筋較多的Reckling No.23號模型作為分析對象,模型跨長l=500 mm。表1列出了模型截面尺寸和材料屬性。

圖3 Reckling有限元模型一階屈曲模態

表1 Reckling No.23模型截面尺寸和材料屬性

在進行極限強度分析時,可將如圖3所示模型的第一階線性屈曲模態的變形形式作為結構的初始缺陷引入模型網格。

2.2 弧長法解

采用弧長法時,在模型兩端設置大小相同、方向相反的參考載荷(中垂彎矩),對應各弧長增量步的載荷與位移的表達式如式(7)～式(9)所示。得到的端面彎矩-轉角曲線如圖4所示:考察曲線峰值點對應的彎矩,即為采用弧長法得到的極限載荷。

2.3 阻尼因子法解

圖4 載荷-位移平衡路徑曲線(弧長法)

采用阻尼因子法時,在模型端面設置的載荷應超過結構的極限載荷,一般可取為端面的塑性彎矩。在結構達到極限狀態前,采用阻尼因子法和采用弧長法得到的端面彎矩-轉角曲線幾乎重合;當結構達到極限狀態時,采用阻尼因子法得到的端面彎矩-轉角曲線變成一條水平線,如圖5(a)所示。圖5(b)、圖5(c)表達了中垂加載過程中,阻尼力項所做的功和結構應變能隨端部施加的中垂彎矩和端部轉角的變化曲線。通過綜合考察圖5所示的端面彎矩-轉角曲線、彎矩-能量曲線和轉角-能量曲線來確定結構的極限承載能力。

圖5 阻尼因子法

2.4 準靜態法解

采用準靜態法時,在模型兩端施加隨時間光滑變化的強迫位移(轉角)。當總的強迫位移值一定時,加載速率取決于加載的總時間。當加載時間為1 s時的端部強迫位移加載曲線如圖6(a)所示。這里取加載的總時間為1/100 s、1/20 s和1/1 s由快到慢三個加載速率進行計算,輸出端面的彎矩反力,得到如圖6(b)所示的彎矩-轉角曲線。同時考察如圖6(c)所示的結構動能和應變能歷史曲線,確定合適的加載速率。

圖6 準靜態法

2.5 討論

表2列出了三種數值解法的計算結果及其相對試驗值的誤差。計算表明,雖然三種有限元解法具體的計算思路和過程有所不同,但最終得出的結果基本一致,且與試驗值相吻合。下面對不同解法的應用特點與分析的過程進行討論:

表2 數值計算結果與試驗值比較

(1)采用弧長法應對結構模型的規模加以控制,以避免出收斂性問題。

(2)采用阻尼因子法應同時考察載荷-位移曲線和阻尼力項所做的功占應變能百分比的變化曲線。當結構到達極限狀態時,阻尼力項所做的功將迅速增加,甚至超過結構的應變能。這表明阻尼力項在平衡方程中占據了主導作用,此時對應的載荷即為結構的極限載荷。

(3)采用準靜態方法應取由快到慢若干個加載速率進行分析。加載曲線要求其一階和二階導數都是光滑的,從而避免由于加載不連續引起的求解波動。通過考察加載過程中結構動能與應變能之比來判斷加載速率是否合適。計算表明,采用合適的加載速率得到的準靜態結果與靜態計算的結果接近。

3 總結

本文對求解結構極限承載能力的非線性有限元法中的3種解法:弧長法、阻尼因子法和準靜態法進行了理論介紹,并通過對Reckling No.23模型試驗的數值計算,對3種解法的求解思路、計算特點和關鍵點進行了歸納。計算分析表明,3種解法都是有效的數值方法,可以相互作為輔助和補充。當求解模型規模較大時,可以考慮采用阻尼因子法和準靜態方法進行求解;當求解涉及復雜的結構接觸、材料失效等不容易收斂因素時,采用準靜態方法可能更為有效。表3總結了3種非線性有限元解法各自的特點。

表3 3種非線性有限元解法的特點

[1] Jeom Kee Paik,Anil Kumar Thayamballi.Ultimate Limit State Design of Steel-Plated Structures[M].Chichester,U K:Wiley,2003.

[2] ISSC.Report of ISSC Special Task Committee VI.2-Ultimate Hull Girder Strength[R].Proceedings of ISSC 2000[C].Nagasaki,Japan,2000.

[3] 李元齊,沈祖炎.弧長控制類方法使用中若干問題的探討與改進[J].計算力學學報,1998,15(4):414-422.

[4] ABAQUS Analysis User’s Manual.7.1.1 SOLVING NONLINEAR PROBL EM[S].Using the damping factor.

[5] 王勖成.有限單元法[M].北京:清華大學出版社,2003.

[6] 莊茁,張帆,岑松.ABAQUS非線性有限元分析與實例(ABAQUS數碼工程師系列叢書)[M].北京:科學出版社,2005.

[7] Reckling K A.Behavior of box girder under bending and shear[R].Proceeding of the International Ship&Offshore Structures Congress(ISSC),Pairs,1979,2:46-49.

Study on Three FEA Methods of Structure Ultimate Strength Analysis

PENG Da-wei, ZHANG Shi-lian

(Shanghai Jiao Tong University,Shanghai 200030,China)

Three algorithms,Arc-Length method,Damping Factor method and Quasi-Static method in nonlinear FEA method of ultimate strength analysis are presented in this paper.The three algorithms have different analysis processes,features and key points,which are discussed and compared with one another through a numeral calculation by Reckling No23 test model.The result of the comparisons gives a reasonable suggestion for better use of nonlinear FEA method of ultimate strength analysis.

Arc-Length method,Damping Factor method,Quasi-Static method,ultimate strength

U661.41

A

圖1 非線性有限元分析流程

2009-11-03;修改稿收到日期:2009-12-25

彭大煒(1984-),男,碩士研究生,主要從事船體結構極限強度分析研究。

1001-4500(2010)02-0001-05