王保森，馮 亮，耿保陽

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

0 引 言

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

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

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

1 弧長法

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

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

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

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

圖 1 弧長法原理示意圖Fig. 1 Principle diagram of RIKS

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

2 有限元模型

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

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

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

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

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

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

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

表 4 三種箱型梁有限元解和實驗值對比（MPa）Tab. 4 Comparison of finite element solutions and experimental values of three box girders（MPa）

3 有限元計算

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

3.1 邊界條件

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

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

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

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

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）

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

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

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

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

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

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

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

3.3 網(wǎng)格密度

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

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

圖 9 Reckling23 箱型梁模型屈曲模態(tài)及損傷變形、應力云圖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 三種箱型梁在不同缺陷大小下極限彎矩值計算值（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，極限彎矩明顯減小，可見網(wǎng)格粗糙程度對箱型梁極限彎矩有限元解具有較大影響。

圖 11 不同網(wǎng)格密度下Nshihara3 箱型梁模型損傷變形、應力云圖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)格密度下的計算結(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)格密度下的計算結(jié)果，從表中可以看出網(wǎng)格密度越稀疏，計算值越大，在4 種網(wǎng)格尺寸中，Nshihara3模型在網(wǎng)格為5 cm 時誤差最小，僅與實驗值有0.41%的誤差，但在網(wǎng)格尺寸增加到10 cm 時，計算值有較大程度的增長，彎矩坡度從2.15 增到38.13，說明網(wǎng)格大小在5～10 cm 之間時，計算結(jié)果受網(wǎng)格密度影響很大，網(wǎng)格尺寸為2.5～5 cm 之間值，彎矩坡度僅為2.15，綜合誤差值與彎矩坡度值可得出，Nshihara3 模型網(wǎng)格大小為2.5～5 cm 之間時可得到滿足精度且比較穩(wěn)定的結(jié)果，有限元結(jié)果在大于5 cm 時出現(xiàn)跳躍，網(wǎng)格密度為10 cm 以上時，Nshihara 箱型梁模型極限強度誤差值超過30%。

4 結(jié) 語

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

1）在相同荷載條件下，有延長段模型比無延長段模型極限彎矩值更接近實驗值，且無延長段的Dowling2箱型梁模型計算值誤差達到了20%，說明采用延長段的邊界條件模擬箱型極限破壞得到的計算值更準確。

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

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