基于有限元法的箱型梁極限強(qiáng)度影響因素及敏感分析

王保森，馮 亮，耿保陽(yáng)

(1. 山東省船舶控制工程與智能系統(tǒng)工程技術(shù)研究中心，山東 榮成 264300；2. 中國(guó)海洋大學(xué) 工程學(xué)院，山東 青島 266100；3. 威海海洋職業(yè)學(xué)院，山東 榮成 264300)

0 引 言

船體極限強(qiáng)度的研究一般歸結(jié)為對(duì)簡(jiǎn)化箱型梁結(jié)構(gòu)的縱向極限強(qiáng)度計(jì)算，箱型梁的縱向強(qiáng)度可用橫截面承受的最大彎矩來(lái)表示，因此線彈性理論一直被用來(lái)計(jì)算箱型梁的縱向強(qiáng)度。近些年來(lái)，包括屈曲在內(nèi)多種失效形式的出現(xiàn)、結(jié)構(gòu)部件失效產(chǎn)生的相互作用以及橫截面受力重分布等，表明應(yīng)考慮材料非線性及幾何非線性對(duì)箱型梁的極限強(qiáng)度進(jìn)行研究，因此非線性有限元法被廣泛應(yīng)用到船舶極限強(qiáng)度研究中。

國(guó)內(nèi)外許多學(xué)者對(duì)箱形梁結(jié)構(gòu)進(jìn)行了詳細(xì)廣泛的研究，Nshihara，Dowling，Reckling[1-4]對(duì)不同的箱形梁結(jié)構(gòu)進(jìn)行了極限強(qiáng)度研究實(shí)驗(yàn)，通過(guò)對(duì)實(shí)驗(yàn)的結(jié)果進(jìn)行比較分析，改進(jìn)了箱形梁極限強(qiáng)度的理論算法。Hansen[5]采用有限元方法對(duì)Nshihara 的幾個(gè)箱形梁實(shí)驗(yàn)?zāi)Ｐ瓦M(jìn)行了簡(jiǎn)單的研究，其中模型采用的是4 節(jié)點(diǎn)殼單元，通過(guò)不斷增大荷載，使模型發(fā)生整體彎曲，通過(guò)對(duì)對(duì)比5 種箱形梁極限強(qiáng)度實(shí)驗(yàn)結(jié)果和理論算法分別對(duì)應(yīng)的的應(yīng)力-應(yīng)變曲線及彎矩-曲率曲線，發(fā)現(xiàn)結(jié)構(gòu)崩潰段有很大的差異，反映了理論算法的缺陷。Guedes Soares 等[6]詳細(xì)研究了箱形梁在簡(jiǎn)支條件下受純彎矩作用的受力變形過(guò)程，并用Ansys 軟件進(jìn)行了模擬，實(shí)驗(yàn)結(jié)果與有限元結(jié)果符合較好。徐向東等[7]也對(duì)箱形梁進(jìn)行了極限承載能力方面的實(shí)驗(yàn)與理論研究，推出了一個(gè)理論公式，并用實(shí)例進(jìn)行了驗(yàn)證。賀雙元等[8]基于Marc 軟件對(duì)Nshihara 論文中的NST3 箱形梁進(jìn)行極限強(qiáng)度分析。白勇等[9]基于非線性有限元程序Sandy 較全面討論了包括屈服應(yīng)力、楊氏模量、初始缺陷、焊接殘余應(yīng)力、板厚等因素對(duì)船體結(jié)構(gòu)中拱極限強(qiáng)度影響程度，但僅采用一個(gè)算例進(jìn)行分析，普遍性較小。

本文運(yùn)用Abaqus 提供的弧長(zhǎng)法（RIKS）計(jì)算3 個(gè)典型箱型梁極限強(qiáng)度值，并與已有實(shí)驗(yàn)數(shù)據(jù)比對(duì)，驗(yàn)證本文有限元算法的可靠性。通過(guò)計(jì)算考察邊界條件類(lèi)型、網(wǎng)格密度大小和初始缺陷幅值對(duì)數(shù)值計(jì)算結(jié)果的影響，給出具體的誤差值。

1 弧長(zhǎng)法

加弧長(zhǎng)法原理是通過(guò)設(shè)置一個(gè)弧長(zhǎng)參數(shù)來(lái)控制平衡方程的增量迭代和收斂，可以將式（1）寫(xiě)成式（2）增量形式：式中：[KT] 為切線剛度矩陣；{Δu}為位移增量；{ΔP}為載荷增量；{R}為殘差力。

設(shè)第i 步迭代的載荷增量為{ΔP}i，由載荷增量因和參考載荷{Pref}來(lái)控制，即

將式（3）代入式（2）得到第i 步迭代的增量格式：

如圖1 所示，弧長(zhǎng)法在求解中，是把上一步增量計(jì)算的平衡點(diǎn)看做圓心，弧長(zhǎng)增量為半徑，經(jīng)過(guò)牛頓-拉普森迭代找到下一步增量平衡點(diǎn)。每一步的弧長(zhǎng)增量載荷增量因子以及位移增量{Δu}i均通過(guò)下面的約束方程來(lái)控制：

圖 1 弧長(zhǎng)法原理示意圖Fig. 1 Principle diagram of RIKS

經(jīng)過(guò)不斷的迭代，一直到殘差力控制在容差{R}i之內(nèi)。第i 步迭代完成時(shí)有：

2 有限元模型

本文選用Nishihara，Dowling，Reckling 模擬實(shí)船鋼箱梁設(shè)計(jì)的3 種箱型梁結(jié)構(gòu)，3 位學(xué)者均對(duì)模型進(jìn)行了實(shí)驗(yàn)研究，將本文有限元結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，分析這3 種箱型梁極限強(qiáng)度的影響因素及誤差大小。

3 種典型箱型梁即Nshihara 的方形梁MST-3 模型（簡(jiǎn)稱Nshihara3 模型）、Dowling2 箱型梁模型、Reckling23 箱型梁模型的剖面結(jié)構(gòu)形式如圖2 所示。

3 種模型構(gòu)件的詳細(xì)尺寸和材料特性見(jiàn)表1～表3，且3 種模型的跨長(zhǎng)L 分別為540 mm，787 mm，500 mm。

圖 2 三種典型箱型梁剖面尺寸圖Fig. 2 Dimensions of three typical box girders

表 1 Nshihara 3 箱型梁結(jié)構(gòu)尺寸表Tab. 1 Dimensions of Nshihara MST-3 box girder

表 2 Dowling2 箱型梁結(jié)構(gòu)尺寸表Tab. 2 Dimensions of Dowling2 box girder

表 3 Reckling23 箱型梁結(jié)構(gòu)尺寸表Tab. 3 Dimensions of Reckling23 box girder

箱型梁有限元模型采用shell 單元進(jìn)行模擬，單元類(lèi)型S4R5，Nshihara3 箱型梁模型如圖3所示，模型的兩端均使用了相同的延長(zhǎng)段，目的是應(yīng)用St.Venant 原理糾正邊界條件的偏差，在延長(zhǎng)段兩端中心各設(shè)置一個(gè)參考點(diǎn)，兩端節(jié)點(diǎn)采用MPC 多節(jié)點(diǎn)約束，約束條件為：在左端面剛體約束的參考點(diǎn)上約束軸1、軸2、軸3 方向的線位移Ul，U2，U3 以及繞軸1、軸3 的角位移UR1，UR3；右端面剛體約束的參考點(diǎn)上約束軸2、軸3 方向的線位移U2，U3 以及繞軸1、軸3 的角位移URI，UR3，考慮箱型梁表面的初始缺陷，在兩端參考點(diǎn)處施加連續(xù)的垂向彎矩，使梁發(fā)生純彎曲破壞，據(jù)輸出的彎矩-轉(zhuǎn)角曲線可得到箱型梁的極限彎矩值。

圖 3 Nshihara 3 箱型梁模型Fig. 3 Model of Nshihara 3 box girder

為了驗(yàn)證有限元計(jì)算結(jié)果的可信性，以Nshihara3箱型梁模型、Dowling2 箱型梁模型、Reckling23 箱型梁模型為例進(jìn)行計(jì)算，初始缺陷為箱型梁跨長(zhǎng)的0.01，網(wǎng)格密度為5 cm，其中Nshihara3 箱型梁模型結(jié)果如圖4 所示。

圖 4 Nshihara3 箱型梁損傷變形、應(yīng)力云圖及彎矩-轉(zhuǎn)角曲線Fig. 4 Damage deformation and stress cloud diagram bending moment - angle curve of Nshihara3 box girders

3 種箱型梁有限元計(jì)算結(jié)果和實(shí)驗(yàn)值比對(duì)如表4所示，3 種箱型梁有限元解和實(shí)驗(yàn)值的誤差分別為0.92%，0.84%，3.6%，滿足精度要求，驗(yàn)證了有限元算法的可信性。

表 4 三種箱型梁有限元解和實(shí)驗(yàn)值對(duì)比（MPa）Tab. 4 Comparison of finite element solutions and experimental values of three box girders（MPa）

3 有限元計(jì)算

箱型梁極限強(qiáng)度有限元算法的影響因素主要為邊界條件、初始缺陷、網(wǎng)格密度等，因此本節(jié)重點(diǎn)研究這3 種箱型梁模型在有限元方法中，不同邊界條件、初始缺陷、網(wǎng)格密度對(duì)其極限強(qiáng)度的影響程度及誤差大小。

3.1 邊界條件

本文在箱型梁兩端各設(shè)置一段與原模型等長(zhǎng)的延長(zhǎng)段，即1+1+1 模型，延長(zhǎng)段采用較高屈服強(qiáng)度的材料且增加板厚，目的防止延長(zhǎng)段過(guò)早破壞。圖5 為Nshihara3 箱型梁模型在無(wú)延長(zhǎng)段邊界條件、1+1+1 邊界條件下的損傷變形、應(yīng)力云圖，從圖中看出模型延長(zhǎng)段無(wú)明顯變形且應(yīng)力較小，變形與應(yīng)力都集中在箱型梁中段，即設(shè)置延長(zhǎng)段不影響原箱型梁的變形及應(yīng)力分布，說(shuō)明設(shè)置的延長(zhǎng)段比較合理。

圖 5 無(wú)延長(zhǎng)段和1+1+1 邊界條件下Nshihara3 箱型梁損傷變形、應(yīng)力云圖Fig. 5 Damage deformation and stress cloud diagram of no extension and 1+1+1 Boundary conditions

為比較Nshihara3 模型、Dowling2 箱型梁模型、Reckling23 箱型梁模型在上述2 種邊界條件下計(jì)算結(jié)果的差異，分別計(jì)算這2 種邊界條件下3 種模型的極限強(qiáng)度。圖6 為3 種箱型梁模型的在上述2 種邊界條件下計(jì)算得到的彎矩-轉(zhuǎn)角曲線圖。

可以看出，無(wú)延長(zhǎng)段模型比1+1+1 模型計(jì)算結(jié)果偏大，尤其對(duì)Dowling2 箱型梁模型影響更為明顯。表5為3 種箱型梁結(jié)構(gòu)在兩種邊界條件下極限強(qiáng)度有限元值與實(shí)驗(yàn)值對(duì)比?？梢钥闯?，2 種邊界條件下極限彎矩差距為19%左右，且1+1+1 模型計(jì)算結(jié)果更接近實(shí)驗(yàn)值，因此采用具有延長(zhǎng)段的邊界條件可以更好地模擬箱型梁的純彎屈曲變形。

3.2 初始缺陷因子

圖 6 兩種邊界條件下各箱型梁的彎矩-轉(zhuǎn)角曲線Fig. 6 Bending moment - angle curve of each box beam under two boundary conditions

表 5 兩種邊界條件箱型梁極限彎矩值（kN·m）Tab. 5 The ultimate bending moment value of two boundary condition box girders（kN·m）

在研究箱型梁的極限強(qiáng)度時(shí)，需要引入合適的初始缺陷，使數(shù)值模擬過(guò)程與計(jì)算結(jié)果更符合實(shí)際，初始缺陷的處理可采用傅里葉級(jí)數(shù)法，通過(guò)改變單元節(jié)點(diǎn)的位移進(jìn)行施加，但施加過(guò)程比較復(fù)雜，為方便施加初始缺陷，本文引入初始缺陷的方法是采用Abaqus軟件提供的Buckle 分析步對(duì)加筋板模型進(jìn)行特征值屈曲分析，通過(guò)對(duì)特征屈曲模態(tài)按照比例因子進(jìn)行合適的縮放，將縮放后的變形結(jié)果作為初始缺陷引入到后屈曲分析中。

比例因子的取值影響著引入缺陷的大小，因此選取比例因子范圍不應(yīng)過(guò)大，考慮到實(shí)際初始缺陷的量級(jí)，本文按0.01～0.02 L（L 為箱型梁跨長(zhǎng)）大小選取比例因子，即初始缺陷5～10 mm，對(duì)3 種箱型梁進(jìn)行分析。

圖7～圖9 分別為Nshihara3 模型、Dowling2 箱型梁模型、Reckling23 箱型梁模型的屈曲模態(tài)與相應(yīng)的箱型梁極限損傷變形、應(yīng)力云圖。可以看出，屈曲模態(tài)的變形區(qū)域與箱型梁最終變形云圖基本一致，屈曲模態(tài)作為結(jié)構(gòu)初始缺陷，決定模型的屈曲變形位置，并影響模型最終的破壞形態(tài)。

圖10 為Nshihara3 模型的彎矩-轉(zhuǎn)角圖，從圖中可以看出，隨著初始缺陷的增大，極限彎矩呈減小趨勢(shì)，且在初始缺陷為5～10 mm 范圍內(nèi)，極限強(qiáng)度計(jì)算結(jié)果差距不大，其他箱型梁的彎矩-轉(zhuǎn)角曲線差距類(lèi)似。

表6 為3 種箱型梁在5～10 mm 等不同缺陷大小下極限彎矩值，可以看到對(duì)應(yīng)不同初始缺陷，各模型的極限彎矩值，從誤差結(jié)果來(lái)看最大誤差達(dá)到6%，最小誤差為0.71%，說(shuō)明在此范圍內(nèi)選擇初始缺陷的大小是比較合理的，從坡度的變化趨勢(shì)可以看出Nshihara3模型和Reckling23 箱型梁模型均為減小的趨勢(shì)，最適合的缺陷值為8 mm，Dowling2 最適合缺陷為10 mm，均在0.01～0.02 倍箱型梁跨長(zhǎng)范圍內(nèi)。

圖 7 Nshihara3 模型屈曲模態(tài)及損傷變形、應(yīng)力云圖Fig. 7 Buckling mode of Nshihara3 model and Damage deformation、stress cloud diagram

圖 8 Dowling2 箱型梁模型屈曲模態(tài)及損傷變形、應(yīng)力云圖Fig. 8 Buckling mode of Dowling2 model and Damage deformation、stress cloud diagram

3.3 網(wǎng)格密度

一般來(lái)說(shuō)，模型網(wǎng)格越密計(jì)算的結(jié)果更精確。然而隨著網(wǎng)格密度增加，相應(yīng)的計(jì)算時(shí)間會(huì)更長(zhǎng)。因此，合適的網(wǎng)格密度既能得到滿足精度的結(jié)果又可削減計(jì)算時(shí)間，這一規(guī)律在大型結(jié)構(gòu)中尤其明顯。

本文選取網(wǎng)格尺寸為2.5 cm，5 cm，10 cm，15 cm且均設(shè)為正方形網(wǎng)格。圖11 為4 種網(wǎng)格密度下Nshihara3 模型損傷變形、應(yīng)力云圖。

圖 9 Reckling23 箱型梁模型屈曲模態(tài)及損傷變形、應(yīng)力云圖Fig. 9 Buckling mode of Reckling23 model and Damage deformation、stress cloud diagram

圖 10 Nshihara 3 箱型梁不同缺陷下的彎矩-轉(zhuǎn)角曲線Fig. 10 Bending moment - angle curve of Nshihara 3 box girder under different defect

表 6 三種箱型梁在不同缺陷大小下極限彎矩值計(jì)算值（kN·m）Tab. 6 The ultimate bending moment of three kinds of box girders under different defect（kN·m）

圖12 為2.5 cm、5 cm、10 cm、15 cm 4 種網(wǎng)格密度下Nshihara3 模型彎矩-轉(zhuǎn)角曲線，從圖中可以看出網(wǎng)格密度為2.5 cm 與5 cm 極限彎矩大小基本相同，網(wǎng)格密度為10 cm 與15 cm 極限彎矩大小也基本相同，按網(wǎng)格密度從5 cm 變大到10 cm，極限彎矩明顯減小，可見(jiàn)網(wǎng)格粗糙程度對(duì)箱型梁極限彎矩有限元解具有較大影響。

圖 11 不同網(wǎng)格密度下Nshihara3 箱型梁模型損傷變形、應(yīng)力云圖Fig. 11 Damage deformation and stress cloud diagram of Nshihara3 box girder model under different grid densities

圖 12 四種網(wǎng)格密度下Nshihara 箱型梁模型彎矩-轉(zhuǎn)角曲線Fig. 12 Bending moment - Angle curve of Nshihara box girder under four kind of grid densities

表 7 三種箱型梁在不同網(wǎng)格密度下的計(jì)算結(jié)果（kN·m）Tab. 7 The ultimate bending moment of three kinds of box girders under different grid densities（kN·m）

表7 為3 種箱型梁在2.5 cm、5 cm、10 cm、15 cm網(wǎng)格密度下的計(jì)算結(jié)果，從表中可以看出網(wǎng)格密度越稀疏，計(jì)算值越大，在4 種網(wǎng)格尺寸中，Nshihara3模型在網(wǎng)格為5 cm 時(shí)誤差最小，僅與實(shí)驗(yàn)值有0.41%的誤差，但在網(wǎng)格尺寸增加到10 cm 時(shí)，計(jì)算值有較大程度的增長(zhǎng)，彎矩坡度從2.15 增到38.13，說(shuō)明網(wǎng)格大小在5～10 cm 之間時(shí)，計(jì)算結(jié)果受網(wǎng)格密度影響很大，網(wǎng)格尺寸為2.5～5 cm 之間值，彎矩坡度僅為2.15，綜合誤差值與彎矩坡度值可得出，Nshihara3 模型網(wǎng)格大小為2.5～5 cm 之間時(shí)可得到滿足精度且比較穩(wěn)定的結(jié)果，有限元結(jié)果在大于5 cm 時(shí)出現(xiàn)跳躍，網(wǎng)格密度為10 cm 以上時(shí)，Nshihara 箱型梁模型極限強(qiáng)度誤差值超過(guò)30%。

4 結(jié) 語(yǔ)

本文研究箱型梁極限強(qiáng)度有限元法的不穩(wěn)定性，選取Nshihara 的方形梁MST-3 模型、Dowling2 模型、Reckling23 模型等3 個(gè)典型箱型梁模型，通過(guò)對(duì)邊界條件、初始缺陷、網(wǎng)格尺寸等因素進(jìn)行分析，得到以下結(jié)論：

1）在相同荷載條件下，有延長(zhǎng)段模型比無(wú)延長(zhǎng)段模型極限彎矩值更接近實(shí)驗(yàn)值，且無(wú)延長(zhǎng)段的Dowling2箱型梁模型計(jì)算值誤差達(dá)到了20%，說(shuō)明采用延長(zhǎng)段的邊界條件模擬箱型極限破壞得到的計(jì)算值更準(zhǔn)確。

2）比例因子取5～10 mm 時(shí)，3 個(gè)模型在不同初始缺陷下計(jì)算誤差為0.71%～6%，Nshihara 模型和Reckling23 模型最適合缺陷值為8 mm，Dowling2 為10 mm，均在0.01～0.02 倍箱型梁跨長(zhǎng)范圍內(nèi)。

3）在2.5～15 cm 網(wǎng)格密度內(nèi)，有限元結(jié)果在大于5 cm 時(shí)出現(xiàn)跳躍，此時(shí)Nshihara 箱型梁模型極限強(qiáng)度誤差值達(dá)到36%，當(dāng)網(wǎng)格密度在0.1 倍箱型梁跨長(zhǎng)值以內(nèi)時(shí)，計(jì)算結(jié)果更精確。