船舶開孔梁局部屈曲的無網格MLPG 分析方法

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

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

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

0 引 言

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

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

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

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

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

1 基于MLPG 的梁屈曲問題

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

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

1.1 位移場逼近函數

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

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

式中：w(x?xI)為 權函數，x=[x y]T為空間坐標系內的向量； ?為節點參數值的加權殘值。

泛函 ?中取極小值，可得

其中

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

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

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

1.2 梁平面響應

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

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

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

由于是基于單元的全梁分析，得到的應力僅在局部單元域連續，故其產生的相關誤差對于典型的梁結構來說可以忽略不計。此外，一個連續平面應力場可以在整體梁內通過MLS 來近似實現。單個單元平面應力場的計算公式如下：

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

1.3 平面外響應

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

1.3.1 材料剛度矩陣

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

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

1.3.2 幾何剛度矩陣

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

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

1.3.3 折邊剛度矩陣

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

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

1.4 簡化后的屈曲分析

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

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

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

因此，2 個互補模態uzA和uzB用于計算和解答經降階了的模態特征值，結果經過每一次迭代計算直至收斂。

2 數值模擬算例

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

圖2 4 個孔的單元結構圖Fig. 2 Structural drawing of four hole units

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

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

圖3 無網格MLPG 法計算結果Fig. 3 Calculation results of meshless MLPG method

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

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

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

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

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

圖6 改型方案結構示意圖Fig. 6 Structural diagram of modification scheme

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

圖7 改型算例無網格MLPG 法計算結果Fig. 7 Calculation results of modified example by meshless MLPG method

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

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

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

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

3 結 論

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

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

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

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