楊 源，莫中華，唐文勇，孫啟榮，沈亞明

（1.南通中遠(yuǎn)海運(yùn)川崎船舶工程有限公司，江蘇南通226005；2.上海交通大學(xué)海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海200240）

0 引 言

在船體梁總縱彎矩作用下，甲板縱桁和船底縱桁承受較高的軸向壓力，結(jié)構(gòu)設(shè)計(jì)時(shí)要確保它們具有足夠的屈曲強(qiáng)度。但是，出于人員、管路穿行和結(jié)構(gòu)減重的目的，縱桁的腹板上不可避免地會(huì)布置一些開(kāi)孔。開(kāi)孔形狀以腰圓形和圓形較為常見(jiàn)，孔尺寸相比板格短邊長(zhǎng)一般較大，此時(shí)開(kāi)孔對(duì)板格屈曲強(qiáng)度的影響值得關(guān)注。由于開(kāi)孔的存在，結(jié)構(gòu)分析域的形狀變得不規(guī)則，面內(nèi)應(yīng)力發(fā)生了重分配，依靠純解析方法難以準(zhǔn)確處理該問(wèn)題。

早期，Schlack[1]采用Ritz法研究了中心開(kāi)圓孔的矩形板格受壓屈曲問(wèn)題，在小孔徑、方形板條件下的計(jì)算結(jié)果與試驗(yàn)結(jié)果大體接近，這種模型過(guò)于理想化，實(shí)用價(jià)值有限。近年來(lái)，El-Sawy[2]、褚洪[3]、李政杰[4]等采用有限元法研究了開(kāi)孔板彈性屈曲壓力的各種敏感性因素。潘祖興[5]以復(fù)變應(yīng)力函數(shù)重構(gòu)法配合區(qū)域劃分法研究了帶矩形開(kāi)孔的矩形板屈曲，結(jié)果準(zhǔn)確性高，只是區(qū)域分解法更適于通過(guò)簡(jiǎn)單分解即可獲得較規(guī)則子區(qū)域形狀的開(kāi)孔結(jié)構(gòu)域，不易推廣到其他開(kāi)孔形狀的情形。本文聚焦船體結(jié)構(gòu)中常見(jiàn)的開(kāi)腰圓孔矩形板格的板邊受壓彈性屈曲問(wèn)題，以復(fù)變函數(shù)方法完成開(kāi)孔板平面應(yīng)力計(jì)算，以基于背景網(wǎng)格的Ritz法計(jì)算彈性屈曲載荷，形成了一種有別于有限元的半解析算法，可應(yīng)用在船體結(jié)構(gòu)設(shè)計(jì)校核中。

1 開(kāi)孔板的平面應(yīng)力計(jì)算

待分析問(wèn)題的力學(xué)模型如圖1所示，在矩形板的中部存在一個(gè)腰圓形的開(kāi)孔，開(kāi)孔板的外板邊可能受到單向或兩向的壓應(yīng)力作用，不考慮體力作用。開(kāi)孔板的長(zhǎng)、寬、厚分別標(biāo)記為a、b、t，孔的長(zhǎng)軸、短軸的長(zhǎng)度分別標(biāo)記為lx、ly，孔心位置以坐標(biāo)( xh,yh)來(lái)表示，E、μ 為板材料的彈性模量和泊松比。計(jì)算開(kāi)孔板格屈曲首先要確定板上任意一處位置的平面應(yīng)力狀態(tài)量( σx,σy,τxy)。

根據(jù)Muskhelishvili提出的復(fù)變函數(shù)解法[6]，體力為0 時(shí)，均勻各向同性材質(zhì)下的彈性平面問(wèn)題的應(yīng)力狀態(tài)解可由兩個(gè)復(fù)解析函數(shù)φ( z )和ψ( z )決定。各向應(yīng)力、面內(nèi)位移與φ( z )、ψ( z )之間存在如下的關(guān)系：

圖1 分析模型示意圖Fig.1 Diagram of analytical model

對(duì)于非圓形孔口的計(jì)算一般要進(jìn)行保角變換處理，將彈性體在z平面上所占的區(qū)域變換為平面上以ξ = 0 為原點(diǎn)，1 為半徑的單位圓內(nèi)部區(qū)域，此時(shí)的保角變換關(guān)系可表示為z = f( ξ )。將其代入φ1( z)和ψ1( z )中，z的函數(shù)變?yōu)棣蔚暮瘮?shù)，函數(shù)的單值解析性質(zhì)不變，有

根據(jù)復(fù)變函數(shù)理論，在環(huán)形區(qū)域內(nèi)單值解析的函數(shù)可展開(kāi)成Laurent 級(jí)數(shù)，且這種展開(kāi)應(yīng)是唯一的。故可設(shè)：

（10）-（11）式中，Ak和Ck為復(fù)常數(shù)，k= 0對(duì)應(yīng)的常數(shù)項(xiàng)對(duì)分析無(wú)影響，可不考慮。為編程方便，上式進(jìn)一步轉(zhuǎn)化為函數(shù)g( ξ )的線性組合，未知量變?yōu)閷?shí)常數(shù)αk和βk，共8n個(gè)。

采用最小二乘邊界配點(diǎn)法[5-9]，構(gòu)建邊界條件的限制方程。邊界面力tx，ty與應(yīng)力σx,σy,τxy的關(guān)系如下：

式中，nx和ny為邊界法矢的方向余弦。本模型包括內(nèi)外兩部分邊界，內(nèi)邊界為腰圓形孔邊，邊界面力為0，外邊界為矩形，邊界面力為法向壓力。將內(nèi)、外邊界等距離散為M個(gè)點(diǎn)，在每個(gè)點(diǎn)j上，用resj表示x,y兩個(gè)方向上應(yīng)力殘差的平方和：

最接近精確解答的參數(shù){ αk} ，{βk} 應(yīng)使各點(diǎn)處的resj之和最小，故有

上式形成8n個(gè)線性方程，以求解未知參數(shù){αk} ，{ βk} ，k為1至4n的整數(shù)，從而進(jìn)一步確定結(jié)構(gòu)域內(nèi)每一點(diǎn)處的平面應(yīng)力狀態(tài)量。

2 開(kāi)孔矩形板的彈性屈曲分析

假設(shè)在板邊的面內(nèi)壓力作用下，開(kāi)孔板正處于中性平衡狀態(tài)。擾動(dòng)引起結(jié)構(gòu)發(fā)生偏離初始構(gòu)形的微彎，結(jié)構(gòu)域Ω內(nèi)任一點(diǎn)( x,y )的撓度記為w，假設(shè)：

式中，Aj為未知常數(shù)，qj( x,y )為滿足結(jié)構(gòu)位移邊界條件的基函數(shù)。基函數(shù)的選取與板的面外邊界條件有關(guān)，本文主要關(guān)注四邊簡(jiǎn)支情況，不妨繼續(xù)采用經(jīng)典矩形板的基函數(shù)構(gòu)造形式：

此時(shí)，腰圓孔的出現(xiàn)增加了一道內(nèi)圈的孔邊界，模型的孔邊是完全自由的，撓度未受限制，因此（18）式是滿足位移邊界條件的，但孔邊彎矩和剪力均為0，這兩項(xiàng)力邊界條件則不再能夠通過(guò)假設(shè)的基函數(shù)自動(dòng)滿足。因此為獲得較準(zhǔn)確的結(jié)果，需要增加N的數(shù)目，后面的算例說(shuō)明22項(xiàng)的基函數(shù)組合已可獲得足夠準(zhǔn)確的解。

開(kāi)孔板的彎曲變形能Vε可表達(dá)為

而在板邊單位壓力作用下外力功W為

式中，D = Et3/[ 12( 1 - μ2)]，表示板的撓曲剛度。結(jié)構(gòu)位能U = Vε- λW，λ代表施加的板邊壓力載荷。根據(jù)位能駐值原理，有

中性平衡狀態(tài)下板結(jié)構(gòu)撓度解的不唯一，表現(xiàn)為方程組系數(shù)矩陣的奇異，使屈曲臨界載荷的計(jì)算轉(zhuǎn)化為矩陣特征值的計(jì)算。對(duì)開(kāi)腰圓孔的矩形板，板內(nèi)各處的平面應(yīng)力σx，σy和τxy不再是簡(jiǎn)單的均勻分布或線性分布關(guān)系，積分域Ω 的形狀也是不規(guī)則的。顯然，直接從整體的積分域出發(fā)來(lái)計(jì)算（22）-（23）式的積分比較困難，因此將結(jié)構(gòu)域離散為三角形來(lái)實(shí)施積分計(jì)算。

這時(shí)結(jié)構(gòu)域離散生成的網(wǎng)格屬于背景網(wǎng)格，不影響板結(jié)構(gòu)撓曲面的連續(xù)性。在離散得到的三角形積分域內(nèi)可使用二維Gauss積分方法計(jì)算積分，矩陣MA，MB的元素可表示為

上式采用二維四點(diǎn)Gauss 積分公式[10]，ws為第s 個(gè)積分點(diǎn)的權(quán)重，Ar為離散得到的編號(hào)為r 的三角形面積，(,)為編號(hào)r 的三角形中第s 個(gè)積分點(diǎn)在總體坐標(biāo)系下的坐標(biāo)。每個(gè)積分點(diǎn)位置處的面內(nèi)應(yīng)力σx、σy，τxy由上一節(jié)所述的平面應(yīng)力計(jì)算方法確定。

3 典型算例分析

參照?qǐng)D1 取一典型的船體桁材板格，ls= 0.8 m，ly= 0.6 m，a = 2.55 m，b = 0.85 m，t = 0.01 m，孔的中心與板的中心重合，矩形板的長(zhǎng)板邊不受面力，短邊受1 MPa的均布?jí)毫ΑＦ矫鎽?yīng)力計(jì)算時(shí)，短板邊等分為20段，長(zhǎng)板邊等分為60段，腰圓形孔邊按1/4圓弧長(zhǎng)等分20段的密度離散，然后取每一段的中點(diǎn)位置作為邊界計(jì)算的配點(diǎn)，Laurent展開(kāi)的參數(shù)n取30。計(jì)算在Matlab 7.7下編程實(shí)現(xiàn)。開(kāi)孔腰圓孔的保角變換參照文獻(xiàn)［11］的做法進(jìn)行，獲得的保角變換關(guān)系如（26）式，變換前后的內(nèi)外邊界線見(jiàn)圖2。

圖2 保角變換示意圖Fig.2 Diagram of conformal transformation

Matlab 編程計(jì)算開(kāi)孔板的平面x 向應(yīng)力如圖3 所示?？紤]到應(yīng)力分布的對(duì)稱性，在如圖4 所示的1/4 對(duì)稱腰圓孔邊上，取弧長(zhǎng)等分點(diǎn)處的孔邊切向應(yīng)力σθ結(jié)果，與ABAQUS 6.10 下的有限元平面應(yīng)力解對(duì)比，顯示于表1。作為對(duì)比基準(zhǔn)用的開(kāi)孔板格有限元模型，腰圓孔邊網(wǎng)格尺度取0.012 m，外矩形板邊網(wǎng)格尺度取0.025 m，在兩者間的區(qū)域網(wǎng)格尺度逐漸過(guò)渡，以確保得到的孔邊應(yīng)力足夠準(zhǔn)確。表1中除了序號(hào)為7、8 的兩個(gè)計(jì)算點(diǎn)由于應(yīng)力結(jié)果的絕對(duì)值較小而反映出較大的相對(duì)誤差外，其他11 個(gè)計(jì)算點(diǎn)的相對(duì)誤差率均在4%以內(nèi)?？紤]到應(yīng)力集中效應(yīng)在孔邊最大，開(kāi)孔板上其他位置處應(yīng)力的誤差將更小，可見(jiàn)復(fù)變函數(shù)解法的開(kāi)孔板格應(yīng)力計(jì)算是比較準(zhǔn)確的。

圖3 短邊受壓下的開(kāi)孔板σy應(yīng)力云圖 Fig.3 σy stress contour of perforated plate under longitudinal axial compression

圖4 孔邊應(yīng)力σθ計(jì)算點(diǎn)示意圖Fig.4 Diagram of calculation points for tangential stress σθ around the oval hole

表1 孔邊應(yīng)力σθ計(jì)算對(duì)比Tab.1 Comparison of tangential stress σθ around the oval hole

續(xù)表1

屈曲載荷計(jì)算環(huán)節(jié)的域網(wǎng)格離散采用較粗的三角形網(wǎng)格，基于Delauney 三角化方法實(shí)現(xiàn)，如圖5所示。計(jì)算采用（24）式的構(gòu)造形式，以22 項(xiàng)基函數(shù)的線性組合來(lái)近似板的撓曲面，具體參數(shù)見(jiàn)表2。

有限元特征值屈曲作為對(duì)比方法，在ABAQUS 6.10 下建立2 種網(wǎng)格的有限元模型。第一種如圖6 所示，為global mesh size=0.05 m 下以Free 方式劃分得到的有限元網(wǎng)格模型，以四邊形單元（S4R）為主，輔以三角形單元（S3），記作有限元模型I。這種網(wǎng)格較密且均勻，經(jīng)過(guò)驗(yàn)證更密的網(wǎng)格對(duì)計(jì)算結(jié)果影響已經(jīng)比較微弱，可認(rèn)為這種網(wǎng)格密度下的分析結(jié)果準(zhǔn)確度是足夠的，以此作為對(duì)比的基準(zhǔn)模型。第二種采用如圖4 所示的三角形背景網(wǎng)格，生成三角形單元（S3）離散的有限元模型，記作有限元模型II。兩種有限元模型在加載時(shí)均采用shell edge load 方式施加板邊的面內(nèi)壓力，并在矩形的3 個(gè)頂點(diǎn)處施加適當(dāng)?shù)拿鎯?nèi)約束，限制模型的剛體位移，特征值屈曲計(jì)算選用蘭索斯迭代算法完成。

圖5 算例采用的積分背景網(wǎng)格Fig.5 Integration background mesh of the example

圖6 對(duì)比用的有限元模型Fig.6 Finite element model for comparison

表2 算例采用的基函數(shù)序列Tab.2 Basis functionse quence of the example

表3 算例的對(duì)比結(jié)果Tab.3 Comparison of results for the examples

計(jì)算結(jié)果如表3 所示，本文方法的屈曲臨界壓力解與有限元模型I 解比較接近，存在-0.56%的誤差。本文方法得到的1 階屈曲變形模態(tài)如圖7 所示，亦與有限元模型I的基本一致。負(fù)偏差則可能來(lái)自域離散的數(shù)值積分，這一點(diǎn)通過(guò)加密背景網(wǎng)格的驗(yàn)算可初步驗(yàn)證。若將M2背景網(wǎng)格的密度分別增加為原來(lái)的2 倍和3 倍，相應(yīng)的彈性屈曲解分別為116.92 MPa 和116.89 MPa，與有限元模型I 的解已非常接近。而有限元模型II 的解則出現(xiàn)了7%以上的明顯結(jié)果偏差，這反映出背景網(wǎng)格和有限元網(wǎng)格存在著明顯差別。

圖7 本文方法獲得的一階屈曲模態(tài)變形Fig.7 First-order buckling mode deformation obtained by this method

有限元方法采用分片位移場(chǎng)假設(shè)，網(wǎng)格離散對(duì)位移場(chǎng)連續(xù)性的削弱在網(wǎng)格較稀疏時(shí)會(huì)引起明顯的誤差，需要通過(guò)足夠密的網(wǎng)格來(lái)保證精度。而本文方法在平面應(yīng)力和屈曲載荷計(jì)算中均采用了完全連續(xù)的位移場(chǎng)假設(shè)，網(wǎng)格的離散只是為了不規(guī)則形狀域下的數(shù)值積分計(jì)算，故本文方法對(duì)于網(wǎng)格離散的依賴度較低。

4 敏感性對(duì)比驗(yàn)證

在展示開(kāi)孔板格彈性屈曲載荷敏感性的同時(shí)，本章進(jìn)一步展示本文方法的準(zhǔn)確性。敏感性對(duì)比計(jì)算中，均采用與上一節(jié)相似的網(wǎng)格離散。

首先，通過(guò)設(shè)定不同的x、y 向板邊壓應(yīng)力比，可獲得如圖8 所示的彈性屈曲結(jié)果曲線。以ABAQUS 解為基準(zhǔn)，本文方法的結(jié)果誤差率在2%以內(nèi)。誤差率隨著長(zhǎng)邊受壓比重的增加有微弱增加。板格在短邊受壓為主時(shí)，可獲得較高的屈曲承載能力，故這種情況在船體結(jié)構(gòu)校核上也更常見(jiàn)，下面的計(jì)算主要針對(duì)短邊受壓情形展開(kāi)對(duì)比。

圖8 雙向軸壓作用的彈性屈曲解Fig.8 Elastic buckling solution under bi-axial compression

將腰圓孔沿x軸方向朝左側(cè)逐漸移動(dòng)，分別計(jì)算xh從1.275 m逐步減小至0.575 m的彈性屈曲壓力解，結(jié)果如表4。當(dāng)xh大于0.675 m時(shí)，本文方法與ABAQUS有限元解的吻合度較好，以ABAQUS解為基準(zhǔn)，本文方法的結(jié)果誤差率保持在1%以內(nèi)。當(dāng)開(kāi)孔距離孔邊較近時(shí)，本方法的誤差開(kāi)始增大，這與有限的基函數(shù)選取不能反映過(guò)于局部化的撓曲變形有關(guān)。

表5～7分別反映了腰圓孔形狀、矩形板長(zhǎng)度a和寬度b的變化對(duì)屈曲臨界載荷的影響，以ABAQUS解為基準(zhǔn)，本文方法的結(jié)果誤差率均在1%以內(nèi)?？梢钥闯?，以腰圓孔倍率表示腰圓孔長(zhǎng)軸長(zhǎng)度與短軸長(zhǎng)度之比，腰圓孔倍率越高，彈性屈曲解也越大；彈性屈曲解隨矩形板長(zhǎng)度a的增大而逐漸減小，但a>2.85 m后，彈性屈曲解的減小幅度已十分微弱；彈性屈曲解隨矩形板寬度b的增大而逐漸減小。

表4 腰圓孔位置的敏感性效應(yīng)Tab.4 Sensitivity to oval cut-out location

表5 腰圓孔倍率的敏感性效應(yīng)Tab.5 Sensitivity to oval cut-out shape

表6 矩形板長(zhǎng)邊長(zhǎng)度a的敏感性效應(yīng)Tab.6 Sensitivity to rectangular plate long-side length

表7 矩形板短邊長(zhǎng)度b的敏感性效應(yīng)Tab.7 Sensitivity to rectangular plate short-side length

5 結(jié) 語(yǔ)

對(duì)船體桁材開(kāi)孔腹板的屈曲問(wèn)題，本文提出以改進(jìn)的復(fù)變函數(shù)方法計(jì)算板內(nèi)的平面應(yīng)力，結(jié)合背景網(wǎng)格積分下的Ritz 法，確定開(kāi)孔板格的彈性屈曲載荷。該方法基于完全連續(xù)的面內(nèi)和面外位移場(chǎng)假設(shè)，對(duì)網(wǎng)格的依賴性相比有限元已有明顯的削弱，可以在較粗的域離散網(wǎng)格下獲得十分準(zhǔn)確的彈性屈曲解。通過(guò)與ABAQUS 屈曲有限元結(jié)果的系列對(duì)比，本方法能夠較好地反映雙向受壓、孔位移動(dòng)、腰圓孔倍率、矩形板長(zhǎng)度和寬度變化的敏感性影響，具有良好的結(jié)果精度，可用于設(shè)計(jì)階段船體桁材開(kāi)孔板格的屈曲校核。