端木玉，陳建平*，易燕，宋博

1 廣州航海學(xué)院 船舶與海洋工程學(xué)院， 廣東 廣州 510725

2 廣州航海學(xué)院 港口與航運管理學(xué)院，廣東 廣州 510725

0 引 言

在船舶結(jié)構(gòu)設(shè)計中，為減輕結(jié)構(gòu)自身重量，便于設(shè)備、管線和人員通行，在結(jié)構(gòu)強度允許的前提下，會使用大量開孔梁和開孔板結(jié)構(gòu)。在進行船舶修理時，如需拆換船舶內(nèi)部設(shè)備，在大型設(shè)備需要出入艙的過程中，有時也需要在結(jié)構(gòu)上開設(shè)大開口。這些孔洞的存在，會削弱結(jié)構(gòu)自身強度，造成局部的應(yīng)力集中甚至是結(jié)構(gòu)失穩(wěn)。因此，有必要對開孔結(jié)構(gòu)重新進行強度校核。

船舶結(jié)構(gòu)上的孔洞除會削弱結(jié)構(gòu)強度外，還容易產(chǎn)生局部屈曲效應(yīng)。如在高剪切應(yīng)力作用下，大部分開有不規(guī)則孔洞的板和梁（如腹板和面板）會因孔洞的原因而受到較高的集中載荷，可能會產(chǎn)生Web-Post 屈曲、Vierendeel Bending 和破壞[1-2]，從而導(dǎo)致明顯不同的結(jié)構(gòu)載荷響應(yīng)和屈曲破壞力[3-5]。對船舶開孔結(jié)構(gòu)的校核，各大船級社都有相關(guān)規(guī)定。針對這些開孔板和開孔梁的設(shè)計及評估方法依舊局限于BS5950：Part 1[6]和Eurocode 3[7]，這些規(guī)則與標(biāo)準考慮到工程安全，大部分是通過逐步校核梁、板的穩(wěn)性及橫截面臨界能力來進行評判的，方法通常偏于保守。針對具有一定幾何形狀（包括布局和尺寸范圍）的有限域的問題，很大程度上是基于簡化的分析方法[2,8-9]，而這些模型計算強度一般較大，在實踐中較難推廣。

在船舶結(jié)構(gòu)分析的傳統(tǒng)方法中，有限元法是最常用也是最為有效的計算方法，其在處理船舶結(jié)構(gòu)變形、應(yīng)力和疲勞分析等方面都有著成熟和廣泛的應(yīng)用[10-12]，但在分析船舶結(jié)構(gòu)場量（如位移場和應(yīng)力場等）劇烈變化的高梯度區(qū)域時，會出現(xiàn)計算精度降低甚至計算中斷等現(xiàn)象。通常，采用加密該區(qū)域網(wǎng)格或者采用高階單元來解決此問題。但是，這給有限元法前、后處理帶來了巨大的工作量，導(dǎo)致計算效率降低，同時，這種方法并不能從根本上消除問題的根源。

與有限元法原理相對應(yīng)的另一種數(shù)值分析方法是無網(wǎng)格法，該方法在近20 年中也得到了很大的發(fā)展[13-16]。無網(wǎng)格法是建立在系列獨立的離散節(jié)點上，通過構(gòu)造這些離散點的近似函數(shù)來對問題域進行求解。無網(wǎng)格法在計算過程中可以根據(jù)需要任意增加節(jié)點，不需要處理節(jié)點之間的拓撲信息，特別適合自適應(yīng)分析計算。同時，由于無網(wǎng)格法的節(jié)點可以由計算機以自主的方式進行創(chuàng)建，故可節(jié)省花費在創(chuàng)建和處理網(wǎng)格上的時間和計算資源。無網(wǎng)格法目前在航空材料、高速碰撞、動態(tài)裂紋擴展、加工成型等諸多領(lǐng)域得到了較為成功的應(yīng)用[17-19]。

鑒于上述背景，基于無網(wǎng)格局部Petrov-Galerkin（MLPG）法，擬提出開孔結(jié)構(gòu)局部屈曲問題的無網(wǎng)格數(shù)值分析方法。首先，基于MLPG 法定義結(jié)構(gòu)的問題域并進行離散，在離散節(jié)點上建立其緊支函數(shù)，采用加權(quán)余量法建立系統(tǒng)位移場（應(yīng)力場）的離散方程；然后，運用移動最小二乘法（MLS）逼近離散節(jié)點的形函數(shù)；最后，通過算例評估該方法對開孔梁板彈性屈曲問題分析的有效性和正確性。

1 基于MLPG 的梁屈曲問題

在目前的研究中，分析板屈曲問題的一般計算公式是在切線剛度矩陣（KT）奇異性的基礎(chǔ)上發(fā)展起來的，一般由問題域材料剛度矩陣（KE）和幾何剛度矩陣（KG）組成。KE只是單純基于Kirchhoff板理論，KG的計算僅僅需要一階運動學(xué)方程的彈簧旋轉(zhuǎn)類比（RSA）原理類推簡化而來[20]。在預(yù)先給定的確定平面應(yīng)力的線性分析中，容易獲得這些參數(shù)。作為大多數(shù)彈性結(jié)構(gòu)問題，KE可以被認為是獨立的內(nèi)部力量。這種方法雖然考慮了材料的型線響應(yīng)，但在發(fā)生變形屈曲之前，通常會忽略小的變形，從而把屈曲問題作為一個線性特征值問題，直接得出其簡化計算公式。

本文采用MLPG 法來離散問題域，然后使用MLS 構(gòu)建問題域位移場的近似函數(shù)。在建立開孔梁的嚴格屈曲模式過程中，文章采用假定模式與互補模式結(jié)合的方法，組成一個降低特征值求解的二階問題，然后通過迭代的方法[21]使多自由度問題能夠得到有效解決。在迭代運算中，局部域的應(yīng)用在求解精度和計算要求方面對兩者有很大的幫助[22-24]。

1.1 位移場逼近函數(shù)

不同于有限元法所采用的分段函數(shù)（位移函數(shù)的形函數(shù)）預(yù)定義在所有單元子域上，MLPG 法位移場函數(shù)是通過MLS 獲得的。這種方法的主要特征是，基于某一問題域 ?內(nèi)一系列給定節(jié)點的曲線擬合來獲得一個連續(xù)逼近函數(shù)，采用多項式組合來完成：

式中：uh(x)為求解區(qū)域場函數(shù)的近似函數(shù)；pT(x)為 多 項 式 基函 數(shù)；a(x)為 待 定系 數(shù)。對于xy平面內(nèi)的問題，pT(x)=[1x y x2xy y2]，可以看做為一個二次基函數(shù)。參考文獻[21]，采用MLS方法規(guī)定的節(jié)點參數(shù)，即節(jié)點位置xI和節(jié)點位移uI(I=1, 2,···,N)來構(gòu)造近似函數(shù)；a(x)可以通過加權(quán)L2 范數(shù)最小化來建立，可以表示為

式中：w(x?xI)為 權(quán)函數(shù)，x=[x y]T為空間坐標(biāo)系內(nèi)的向量； ?為節(jié)點參數(shù)值的加權(quán)殘值。

泛函 ?中取極小值，可得

其中

定義 Φ(x)為由節(jié)點參數(shù)得到的場逼近函數(shù)形函數(shù)矩陣，設(shè) Φ(x)為

由此，式（1）可進一步簡化為

權(quán)函數(shù)在MLS 中具有非常重要的作用，它可以用來控制和調(diào)節(jié)影響域的位置，一般要求其具有緊支性。本文采用三次樣條函數(shù)W(r)來表達權(quán)函數(shù) ：

1.2 梁平面響應(yīng)

為了求得載荷梁的外平面幾何剛度矩陣KG，需要得知與平面響應(yīng)相關(guān)的應(yīng)力分布。在使用MLPG 法離散整個問題域時，有必要簡化離散過程，這樣對于那些具有相同幾何形狀的單元來說，可以避免重復(fù)計算。本文提出將開孔梁分割成梁單元以進行簡化，其中每個單元便構(gòu)成一個被降低了自由度的超級單元。圖1 所示為一個簡化單元，共有4 個質(zhì)心超級節(jié)點，節(jié)點位置如圖所示，每個節(jié)點有3 個平面自由度 (u,v,θ)，每個單元受到集中載荷F和力矩M的作用。在不規(guī)則開口情況下，需要使用更多的單元模型來表示每個不同的開口。

圖1 簡化單元力學(xué)模型Fig. 1 Simplified element mechanical model

通過MLPG 法對單個單元離散建立起其特征單元響應(yīng)，然后再利用一個標(biāo)準的離散裝配過程[23]，把所有單元組裝成整體的梁響應(yīng)。特征單元的響應(yīng)通過考慮單元力交變載荷的情況獲得。在節(jié)點1 處，當(dāng)單元在外載荷作用下達到平衡時，節(jié)點2，3，4 每一個單元對應(yīng)于一個超級單元自由度，如圖1 所示。按照MLPG 方法逼近，這為轉(zhuǎn)化4 個節(jié)點超級單元平面剛度去替代單個單元響應(yīng)提供了一種靈活的方法。

由于是基于單元的全梁分析，得到的應(yīng)力僅在局部單元域連續(xù)，故其產(chǎn)生的相關(guān)誤差對于典型的梁結(jié)構(gòu)來說可以忽略不計。此外，一個連續(xù)平面應(yīng)力場可以在整體梁內(nèi)通過MLS 來近似實現(xiàn)。單個單元平面應(yīng)力場的計算公式如下：

式中：D為 平面應(yīng)力彈性矩陣；ε 為應(yīng)變；Bxy為應(yīng)變位移矩陣；Txy為在MLPG 單元域內(nèi)節(jié)點參數(shù)對應(yīng)超級單元離散載荷和分布載荷的轉(zhuǎn)換矩陣；F為集中載荷 ；ω 為y方向的均布載荷。式（9）封裝了平面離散分析的2 個層次，即MLPG 和超級單元。

1.3 平面外響應(yīng)

由文獻[20]可知，平面外響應(yīng)的切線剛度矩陣KT是 由KE和KG這2 個部分組成的。

1.3.1 材料剛度矩陣

根據(jù)Kirchhoff 薄板理論[22]，運用MLPG 法，KE可 以分解為彎曲剛度矩陣K和施加本質(zhì)邊界條件后得到的罰函數(shù)剛度矩陣Kα：

式中： Γu為 規(guī)定的位移邊界；Rxy為連接節(jié)點位移和節(jié)點旋轉(zhuǎn)自由度的轉(zhuǎn)換矩陣；罰因子 αz和 αθ為常數(shù)，通常其范圍為剛度矩陣最大對角元素值的倍數(shù)，位數(shù)取值范圍為104～1013[23]。從施加約束出發(fā)，罰因子應(yīng)為無窮大，但在實際計算中不可能實現(xiàn)。在實際計算中，如果罰因子取值過小，可能會造成約束不能精確施加；如果取值過大，則計算可能會溢出或者得到病態(tài)解。為簡化計算，本文罰因子的取值為彎曲剛度矩陣K最大對角線元素值的倍數(shù)，且使得 αz=αθ。

1.3.2 幾何剛度矩陣

幾何剛度KG是根據(jù)RSA 方法從平面應(yīng)力場σx， σy和 τxy得 到[20]。使 用 相 同 的 離 散 化 方 法 對KG進行求解，得

式中： Φ(x),x和 Φ(x),y為由MLS 法求得的形函數(shù)的一階偏導(dǎo)數(shù)； σx， σy為x，y方向的應(yīng)力； τxy為xy平面的剪應(yīng)力。

1.3.3 折邊剛度矩陣

折邊對開孔梁的局部屈曲響應(yīng)的影響可以通過一種簡化模型來實現(xiàn)。對于頂部和底部都有折邊的組合工字梁來說，如果這些折邊能夠沿著連接腹板內(nèi)邊緣起到完全約束作用的話，那么通?？梢灾豢紤]它們的材料和幾何剛度；如果有上、下板約束，則可以把板結(jié)構(gòu)進行等效處理從而簡化為梁結(jié)構(gòu)。相對于建立三維問題模型來說，上面提到的簡化方法僅需采用一維模型即可達到簡化折邊響應(yīng)的目的。同時，這種方法還保留了二維平面模型的適用性。

考慮一個一維折邊 Γf，其尺寸為tf×bf，其中tf為折邊厚度，bf為折邊寬度。在忽略次要變形的影響下，折邊材料剛度貢獻來自其扭轉(zhuǎn)剛度GJ，則折邊的材料剛度矩陣為

1.4 簡化后的屈曲分析

目前，針對開孔梁的失穩(wěn)評估模型一般是在僅取梁的某一適當(dāng)部分的基礎(chǔ)上建立起來的。本文采用局部模型，沿梁的長度方向取至少4 個單元，這樣具有更好的分析效果。本文所提方法只需要在無網(wǎng)格MLPG 區(qū)域內(nèi)計算材料和幾何剛度矩陣即可，大幅減少了剛度矩陣的規(guī)模，提高了計算效率。在梁單元局部域內(nèi)，在按照比例施加載荷情況下梁的臨界屈曲載荷因子 λcr可以由相關(guān)的KE和KG的特征值求得，然而相關(guān)矩陣的大小導(dǎo)致求解其特征值的工作量龐大，大幅增加了計算成本。考慮到RSA 方法[20]，當(dāng)梁在屈曲模態(tài)時，λcr也可以作為內(nèi)部能量吸收和轉(zhuǎn)動等效轉(zhuǎn)化的比例系數(shù)。本文所提方法采用實際屈曲載荷因子的近似模態(tài)代替精確模態(tài)的束縛，典型的模態(tài)可以近似地通過對線性結(jié)構(gòu)執(zhí)行屈曲模態(tài)來獲得[20-21]。

考慮一個近似平面外模態(tài)uz由對應(yīng)于外部的平面力fz獲得，忽略屈曲前平面內(nèi)變形的影響，則臨界屈曲載荷因子 λcr可由下式確定。[22]

基于假定的屈曲模態(tài)，對預(yù)定的模式通過進行迭代運算來進一步改進，式（21）可以用來獲得一個初始解的一個近似臨界屈曲載荷因子 λcr。本文方法是通過局部域的手段來降低問題的維數(shù)，從而使得假定的模態(tài)階數(shù)降低，這樣，特征值的求解就變得容易[24-26]。根據(jù)這種方法，相應(yīng)于一個任意的平面外載荷模式fzA，最初的近似模態(tài)uzA在一個新的負載模式fzB作用下可以建立一個互補的模態(tài)uzB。可以定義為：

因此，2 個互補模態(tài)uzA和uzB用于計算和解答經(jīng)降階了的模態(tài)特征值，結(jié)果經(jīng)過每一次迭代計算直至收斂。

2 數(shù)值模擬算例

算例采用一連續(xù)多孔板。在板沿長度方向，等距開有半徑R=250 mm的系列圓形孔，對稱半長L=40 000 mm ，寬度B=2 000 mm， 板厚t=5 mm。這種結(jié)構(gòu)在船舶結(jié)構(gòu)中非常普遍，例如肋板、船底縱桁等，這些結(jié)構(gòu)的一般特征就是開有一定數(shù)量的減輕孔或者人孔。圖2 所示為該連續(xù)多孔板結(jié)構(gòu)靠近兩端的帶有4.5 個圓形孔的部分簡化形式，圖中數(shù)值的單位為mm。算例中用到的計算常量包括：材料的楊氏彈性模量E=3.0×105MPa，泊 松 比 μ=0.3， 密 度 ρ=7.5×103kg/m3。板 的 上、下端為自由邊，左側(cè)邊界剛性固定，右側(cè)平直邊部分受到與板中性面方向一致的均布壓力ω =2×104N。

圖2 4 個孔的單元結(jié)構(gòu)圖Fig. 2 Structural drawing of four hole units

算例采用MLPG 無網(wǎng)格算法，初始節(jié)點42 108個，如圖3（a）所示；采用高斯積分法，積分背景初始網(wǎng)格83 194 個，經(jīng)過五步自適應(yīng)計算，得到如圖3（b）所示的計算結(jié)果。

為了進行比較分析，算例采用有限元法（NASTRAN 2016）進行了驗證性計算。圖4（a）所示為有限元離散網(wǎng)格，共4 900 個節(jié)點，4 583 個網(wǎng)格，計算得到如圖4（b）所示的計算結(jié)果。

圖3 無網(wǎng)格MLPG 法計算結(jié)果Fig. 3 Calculation results of meshless MLPG method

圖4 有限元法計算結(jié)果Fig. 4 Calculation results of finite element method

從2 種方法的計算結(jié)果來看，計算的應(yīng)力分布高度相似，并且應(yīng)力集中區(qū)域幾乎都落在相同的區(qū)域。從無網(wǎng)格方法的計算結(jié)果來看，由于采用的是連續(xù)應(yīng)力場函數(shù)，不僅可以清晰地得到每一點的應(yīng)力分布，甚至在那些邊緣應(yīng)力出現(xiàn)間斷時也都可以獲得。這就是無網(wǎng)格法相對于有限元的重要優(yōu)勢。

為進一步分析文章所提方法對結(jié)構(gòu)局部屈曲載荷因子分析的有效性，對有限元法與本文所提方法的計算結(jié)果進行了比較，結(jié)果如圖5 所示。從圖中可以看出，有限元法基本上對于整個結(jié)構(gòu)的屈曲載荷因子都不敏感，對于結(jié)構(gòu)當(dāng)中由孔洞所帶來的的局部結(jié)構(gòu)屈曲并沒有明顯的響應(yīng)；而本文所提的運用MLPG 方法計算得出在結(jié)構(gòu)2 個端部的屈曲載荷因子較大，結(jié)構(gòu)中部屈曲載荷因子降低，該結(jié)果與原結(jié)構(gòu)的值吻合較好。

圖5 運用MLPG 法和有限元法計算的屈曲載荷因子變化圖Fig. 5 Change diagram of buckling load factor calculated by MLPG method and finite element method

為了更好地確認本文所提方法能夠更好地分析臨界屈曲載荷因子大小隨結(jié)構(gòu)變化的敏感性，即分析結(jié)構(gòu)屈曲載荷因子隨結(jié)構(gòu)變化而變化的情況，將原結(jié)構(gòu)兩側(cè)第2 和第3 個圓形孔按照以直徑作為邊長修改為了正方形（稱為改型方案），如圖6 所示。圖中，數(shù)值單位為mm。

圖6 改型方案結(jié)構(gòu)示意圖Fig. 6 Structural diagram of modification scheme

運用MLPG 無網(wǎng)格算法進行計算，初始節(jié)點44 400 個，如圖7（a）所示；采用高斯積分法，積分背景初始網(wǎng)格為87 923 個，經(jīng)過五步自適應(yīng)計算，得到如圖7（b）所示的計算結(jié)果。

圖7 改型算例無網(wǎng)格MLPG 法計算結(jié)果Fig. 7 Calculation results of modified example by meshless MLPG method

同樣，采用有限元法（使用NASTRAN 2016）進行比較分析。圖8（a）所示為有限元離散網(wǎng)格，共3 585 個節(jié)點，3 127 個網(wǎng)格，得到如圖8（b）所示的計算結(jié)果。

圖8 改型算例有限元法計算結(jié)果Fig. 8 Calculation results of modified example by finite element method

將經(jīng)過修改的結(jié)構(gòu)再次運用2 種方法予以計算。從計算應(yīng)力云圖結(jié)果來看，應(yīng)力分布依然高度相似，并且應(yīng)力集中趨勢也基本一致。重新計算出的屈曲載荷因子如圖9 所示。有限元法計算結(jié)果顯示結(jié)構(gòu)修改后的載荷因子與原結(jié)構(gòu)載荷因子幾乎相同，說明其對于孔洞的細微修改并不敏感，而本文所提出的方法則非常敏感，兩側(cè)第2 和第3 個孔洞的擴大（由半徑R=250 mm 擴大為邊長a=500 mm 的正方形）明顯降低了其所在位置的屈曲載荷因子，同時也對其兩側(cè)結(jié)構(gòu)產(chǎn)生了一定的影響，而對于其外的結(jié)構(gòu)基本不產(chǎn)生影響，這個計算結(jié)果也與原結(jié)構(gòu)有著更高的契合度。

圖9 運用MLPG 法和有限元法計算的改型算例屈曲載荷因子變化圖Fig. 9 Change diagram of modified example buckling load factor calculated by MLPG method and finite element method

3 結(jié) 論

本文采用MLPG 法對開孔梁和開孔板的彈性屈曲問題予以了計算和評估。在建立開孔梁的嚴格屈曲模式過程中，采用假定模式與互補模式相結(jié)合的方法，降階求解了模態(tài)特征值。在迭代運算中，局部域的應(yīng)用在求解精度和計算要求方面均非常有效。

運用有限元法和本文所提方法對算例進行了計算，通過對應(yīng)力計算結(jié)果的分析和比較，得到2 種方法計算的應(yīng)力分布高度相似，應(yīng)力集中區(qū)域基本相同。本文方法由于采用的是連續(xù)應(yīng)力場函數(shù)，所以可以清晰地得到場點的應(yīng)力分布。

通過對算例中結(jié)構(gòu)局部屈曲載荷因子的計算和比較，得出本文所提方法能夠較好地仿真結(jié)構(gòu)2 個端部的屈曲載荷因子，其結(jié)果與原結(jié)構(gòu)吻合較好。

算例還測試了本文方法和有限元法臨界屈曲載荷因子的大小隨結(jié)構(gòu)變化的敏感性，結(jié)果顯示傳統(tǒng)方法（有限元法）對板梁孔洞的細微修改不敏感，而本文所提方法則非常敏感，該計算結(jié)果與原結(jié)構(gòu)的契合度更好。