基于復(fù)變函數(shù)解法的桁材開孔腹板彈性屈曲分析

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

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

0 引 言

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

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

1 開孔板的平面應(yīng)力計算

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

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

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

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

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

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

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

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

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

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

2 開孔矩形板的彈性屈曲分析

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

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

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

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

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

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

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

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

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

3 典型算例分析

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

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

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

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

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

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

續(xù)表1

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

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

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

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

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

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

計算結(jié)果如表3 所示，本文方法的屈曲臨界壓力解與有限元模型I 解比較接近，存在-0.56%的誤差。本文方法得到的1 階屈曲變形模態(tài)如圖7 所示，亦與有限元模型I的基本一致。負(fù)偏差則可能來自域離散的數(shù)值積分，這一點通過加密背景網(wǎng)格的驗算可初步驗證。若將M2背景網(wǎng)格的密度分別增加為原來的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

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

4 敏感性對比驗證

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

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

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

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

表5～7分別反映了腰圓孔形狀、矩形板長度a和寬度b的變化對屈曲臨界載荷的影響，以ABAQUS解為基準(zhǔn)，本文方法的結(jié)果誤差率均在1%以內(nèi)?？梢钥闯?，以腰圓孔倍率表示腰圓孔長軸長度與短軸長度之比，腰圓孔倍率越高，彈性屈曲解也越大；彈性屈曲解隨矩形板長度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 矩形板長邊長度a的敏感性效應(yīng)Tab.6 Sensitivity to rectangular plate long-side length

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

5 結(jié) 語

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