武卓威， 劉 俊*

(1.上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240；2.上海交通大學(xué) 高新船舶與深海開(kāi)發(fā)裝備協(xié)同中心，上海 200240)

1 引 言

開(kāi)孔在各類工程結(jié)構(gòu)物中較為常見(jiàn)，如船體結(jié)構(gòu)內(nèi)部通常設(shè)有人孔、流水孔、透氣孔和減輕孔等。這些開(kāi)孔的存在會(huì)破壞板殼結(jié)構(gòu)的完整性，降低臨近區(qū)域結(jié)構(gòu)的強(qiáng)度，同時(shí)可能在孔緣產(chǎn)生局部應(yīng)力集中。

目前，在結(jié)構(gòu)領(lǐng)域應(yīng)用最廣泛的數(shù)值模擬方法為有限元法。在處理開(kāi)孔問(wèn)題時(shí)，由于有限元單元必須與開(kāi)孔區(qū)域較為復(fù)雜的幾何構(gòu)型相適應(yīng)，在網(wǎng)格尺寸較大時(shí)，可能存在形態(tài)較差的單元，使計(jì)算精度受到影響。同時(shí),由于有限元網(wǎng)格收斂速度較慢，為了提高求解精度，將顯著增大前處理的工作量和總計(jì)算量。這一矛盾體現(xiàn)了有限元方法在分析此類結(jié)構(gòu)時(shí)的固有缺陷。

數(shù)值流形方法NMM(Numerical Manifold Me -thod)是一種較為新穎的數(shù)值模擬方法，由石根華[1]首次提出，也稱為有限覆蓋法FCM(Finite Cover Method)[2]，該方法通過(guò)兩套獨(dú)立的覆蓋系統(tǒng)構(gòu)造原問(wèn)題的近似解，能夠統(tǒng)一求解連續(xù)和非連續(xù)問(wèn)題[3]，且具有良好的計(jì)算效率和精度[4]。自提出以來(lái)，數(shù)值流形方法主要應(yīng)用于巖土工程領(lǐng)域裂紋[5-7]、大變形[8]和塊體運(yùn)動(dòng)[9,10]等問(wèn)題的求解，近年來(lái)也逐漸開(kāi)始應(yīng)用于結(jié)構(gòu)分析[11]、流體[12]和熱傳導(dǎo)[13]等其他領(lǐng)域，但總體而言，該方法在各類結(jié)構(gòu)分析中應(yīng)用仍相對(duì)較少。

目前，已有部分有關(guān)數(shù)值流形方法的研究涉及平面圓孔算例[14]，但考慮的圓孔尺寸和載荷條件種類均較少，且未提出針對(duì)開(kāi)孔的特殊流形單元，也未與有限元法計(jì)算結(jié)果進(jìn)行詳細(xì)對(duì)比。本文基于數(shù)值流形方法理論框架和平面圓孔問(wèn)題理論解[15]，構(gòu)造一種適用于平面圓孔問(wèn)題的特殊孔緣單元，從而將數(shù)值流形方法應(yīng)用于該類問(wèn)題的求解。

2 數(shù)值流形方法基本原理

數(shù)值流形方法通過(guò)數(shù)學(xué)覆蓋(Mathematical Cover)和物理覆蓋(Physical Cover)這兩套獨(dú)立的覆蓋系統(tǒng)實(shí)現(xiàn)近似解的構(gòu)造和積分區(qū)域的劃分。其中，自行定義的數(shù)學(xué)覆蓋將求解域劃分為若干數(shù)學(xué)片(Mathematical Patch)，而求解域內(nèi)的邊界、裂紋和交界面等將數(shù)學(xué)片進(jìn)行再次剖分，得到相應(yīng)的物理片(Physical Patch)，物理片的集合即為物理覆蓋。物理片的重疊區(qū)域構(gòu)成數(shù)值流形方法的基本單元，稱為流形單元(Manifold Element)。

圖1 數(shù)值流形方法覆蓋系統(tǒng)

由于數(shù)學(xué)覆蓋的劃分具有任意性，本文采用標(biāo)準(zhǔn)網(wǎng)格數(shù)學(xué)覆蓋，其形式如圖2所示。在標(biāo)準(zhǔn)網(wǎng)格中，每一個(gè)頂點(diǎn)均張成一個(gè)矩形數(shù)學(xué)片，每個(gè)流形單元均是四個(gè)不同物理片的重合區(qū)域，如流形單元E1為P1，P2，P4和P5的重合區(qū)域。

圖2 標(biāo)準(zhǔn)網(wǎng)格下數(shù)值流形方法的覆蓋系統(tǒng)

在NMM覆蓋系統(tǒng)中，每一個(gè)物理片具有各自獨(dú)立的局部逼近(Local Approximation)；各個(gè)物理片的局部逼近通過(guò)權(quán)函數(shù)共同構(gòu)成完整求解域的整體逼近(Global Approximation)。物理片Pi的局部逼近是一個(gè)定義在Pi上的函數(shù)空間，函數(shù)空間的基為一給定的基向量Bi(x,y)。平面問(wèn)題中，物理片Pi的局部逼近記為ui(x,y)，可由基向量中的分量函數(shù)張成，即

(1)

每個(gè)物理片Pi均具有對(duì)應(yīng)的權(quán)函數(shù)wi(x,y)。數(shù)值流形方法屬單位分解 PU(Partition of Unity) 類數(shù)值算法，其權(quán)函數(shù)需滿足以下基本要求，

wi(x,y)≥0 [(x,y)∈Pi]

(2)

wi(x,y)=0 [(x,y)?Pi]

(3)

(4)

式中Ω為求解域，n為物理片總數(shù)。在實(shí)際應(yīng)用中，權(quán)函數(shù)?；谟邢拊獑卧魏瘮?shù)構(gòu)造。由此，平面問(wèn)題數(shù)值流形方法的整體逼近表達(dá)為

(5)

獲得了問(wèn)題域的整體逼近后，即可根據(jù)彈性力學(xué)相關(guān)理論，結(jié)合平衡方程、幾何方程和本構(gòu)方程，利用整體逼近推導(dǎo)應(yīng)力和應(yīng)變，隨后依據(jù)最小勢(shì)能原理建立系統(tǒng)的控制方程，并引入邊界條件，該過(guò)程與有限元方法基本一致。

3 平面圓孔特殊孔緣單元的構(gòu)造

3.1 孔緣局部逼近的構(gòu)造

相較于有限元法，數(shù)值流形方法的優(yōu)越性還體現(xiàn)在可以根據(jù)實(shí)際問(wèn)題，選擇合適形式的基構(gòu)成局部逼近，顯著提高求解的準(zhǔn)確性和收斂速度。局部逼近的基可取為常數(shù)、多項(xiàng)式基或其他級(jí)數(shù)的基，基的維數(shù)決定了該局部逼近的自由度數(shù)。對(duì)于一般的平面問(wèn)題，采用常數(shù)基張成的局部逼近即可取得較好的計(jì)算精度，此時(shí)物理片自由度數(shù)為2，可將式(1)以矩陣形式表達(dá)為

(6)

對(duì)于平面開(kāi)孔問(wèn)題，孔緣區(qū)域在外力作用下可能產(chǎn)生應(yīng)力集中，應(yīng)力變化較為劇烈，采用上述形式的局部逼近計(jì)算效果不好，與有限元法相比也沒(méi)有明顯優(yōu)勢(shì)，因此為提高計(jì)算精度，有必要對(duì)局部逼近的基向量進(jìn)行擴(kuò)展。

對(duì)于如圖3所示平板邊界兩對(duì)邊分別承受均布載荷q1和q2的情況，根據(jù)現(xiàn)有彈性力學(xué)解答[15]，由疊加原理可得極坐標(biāo)系下圓孔附近區(qū)域應(yīng)力場(chǎng)解答如式(7～9)。

圖3 四邊受載的帶圓孔矩形平板

(7)

(8)

(9)

式中σr為徑向應(yīng)力，σθ為切向或環(huán)向應(yīng)力，τr θ為剪切應(yīng)力，R為圓孔半徑。

理論上，根據(jù)坐標(biāo)轉(zhuǎn)換關(guān)系可以進(jìn)行應(yīng)力分量轉(zhuǎn)換，隨后結(jié)合彈性力學(xué)物理方程和幾何方程，可導(dǎo)出上述圓孔問(wèn)題的位移解。然而，該過(guò)程非常繁瑣，且上述解答是特定載荷條件下的解答，與實(shí)際情況常有出入?？紤]到局部逼近更重要的是能反映位移收斂特點(diǎn)，且基向量增加項(xiàng)數(shù)不宜過(guò)多，因此在構(gòu)造基向量時(shí)無(wú)需推導(dǎo)出精確的位移表達(dá)式，只需確定能夠反映位移特點(diǎn)的若干項(xiàng)即可。

由此，本文取[1,x,y,cosθ,sinθ]為基向量，并按式(10)的形式構(gòu)造孔緣區(qū)域局部逼近。

(10)

式中Hi為由基向量Bi(x,y)與二階單位矩陣構(gòu)造的基矩陣，如式(11，12)所示。

(11)

(12)

式(10～12)的(x,y)為以圓孔中心為原點(diǎn)的直角坐標(biāo)值。di為自由度向量，每一個(gè)特殊局部逼近的自由度數(shù)為10。

(13)

3.2 流形單元類型的劃分

孔緣特殊物理片的基向量項(xiàng)數(shù)增加，增加了求解計(jì)算量。考慮到應(yīng)力集中是一種局部效應(yīng)，僅需對(duì)孔緣附近的少數(shù)物理片進(jìn)行擴(kuò)展，對(duì)于距離孔緣較遠(yuǎn)的流形單元，則按照一般平面問(wèn)題進(jìn)行處理。

本文依據(jù)單元包含的物理片類型，將問(wèn)題域中的流形單元?jiǎng)澐譃槿?(1) 孔緣單元,構(gòu)成單元的四個(gè)物理片均為孔緣特殊物理片； (2) 過(guò)渡單元,構(gòu)成單元的物理片中有1～3個(gè)特殊物理片； (3) 常規(guī)單元，構(gòu)成單元的物理片全部為普通物理片。三類單元的分布情況如圖4所示，圖中圓點(diǎn)為張成特殊物理片的頂點(diǎn)。

圖4 三類流形單元的劃分

3.3 權(quán)函數(shù)的構(gòu)造

標(biāo)準(zhǔn)網(wǎng)格下，各數(shù)學(xué)片的權(quán)函數(shù)可按式(14,15)的形式，取為包含于該數(shù)學(xué)片的正方形網(wǎng)格在張成該數(shù)學(xué)片頂點(diǎn)處的有限元形函數(shù)構(gòu)成的分塊函數(shù)。

(14)

(15)

至此，將式(14,10)代入式(5)，即可得到問(wèn)題域上的總體逼近，具體表達(dá)為

(16)

式中Ni由權(quán)函數(shù)及基矩陣得到，其作用類似于有限元法的形函數(shù)矩陣

因此有

u=Nd=[N1N2…Nn][d1d2…dn]T

(17)

式(17)即為總體逼近的具體表達(dá)式。

3.4 數(shù)值流形方法計(jì)算流程

通過(guò)局部逼近和權(quán)函數(shù)，由式(5)可以得到求解域的整體逼近，隨后即可依據(jù)最小勢(shì)能原理構(gòu)建控制方程。由于數(shù)值流形方法覆蓋系統(tǒng)中，數(shù)學(xué)覆蓋與問(wèn)題域邊界并不保持一致，本文采用罰函數(shù)法將位移邊界條件引入勢(shì)能表達(dá)式，故平面問(wèn)題下系統(tǒng)的勢(shì)能如式(18)所示。

(18)

(19)

σ=Dε

(20)

式中D為平面應(yīng)力問(wèn)題的彈性矩陣。將式(17)代入式(19,20)，依次得到

ε=LNd=Bd,σ=DBd=Sd

(21)

式中B為應(yīng)變矩陣，S為應(yīng)力矩陣，其物理意義均與有限元法中的應(yīng)變和應(yīng)力矩陣相似。將式(21)代入式(18)，隨后由?ΠP=0即可得到控制方程，其形式可表示為

Kd=q

(22)

式中K為總體剛度矩陣，q為總體載荷列陣，分別由流形單元的單元?jiǎng)偠染仃嘖e和單元載荷列陣qe按編號(hào)順序組裝而成，其中

(23)

(24)

式中上標(biāo)e為B或N對(duì)應(yīng)當(dāng)前流形單元的部分。式(23，24)的積分運(yùn)算可通過(guò)數(shù)值積分實(shí)現(xiàn)，但在局部逼近存在非多項(xiàng)式的基時(shí)，需要適當(dāng)增加積分點(diǎn)數(shù)以保證計(jì)算精度。

將單元?jiǎng)偠染仃嚱M裝后求解方程式(22)，即可得到各局部逼近的精確表達(dá)式，進(jìn)而得到問(wèn)題域上的整體逼近。

4 算例驗(yàn)證

根據(jù)上述理論編制計(jì)算程序，針對(duì)承受三種不同載荷的帶圓孔正方形板進(jìn)行了分析。圓孔位于板的中心，半徑取為R=2 m。算例1和算例3中板的邊長(zhǎng)為2a，算例2中板的對(duì)角線長(zhǎng)為2a，均取a=10 m。板厚t=0.01 m，彈性模量E=2.1×1011MPa，泊松比μ=0.3。各算例采用數(shù)值流形方法和有限元法求解，得到不同網(wǎng)格精度下的解答，并進(jìn)行對(duì)比。

4.1 單向受拉的帶圓孔矩形板

由于該算例的載荷分布具有對(duì)稱性，可取板的1/4進(jìn)行分析計(jì)算，如圖5所示。圖中AE和BC邊施加對(duì)稱位移約束。取均布載荷q=100 MPa。采用n×n標(biāo)準(zhǔn)網(wǎng)格的數(shù)值流形方法計(jì)算結(jié)果與網(wǎng)格尺寸為a/n、由前處理軟件自動(dòng)生成的有限元計(jì)算結(jié)果進(jìn)行對(duì)比，如圖6所示，并取n=2,4,6,…,20。取n=80的有限元計(jì)算結(jié)果作為近似精確解。本算例中孔半徑與板邊長(zhǎng)之比為1/5，達(dá)到彈性力學(xué)理論中視為無(wú)限大板的判斷標(biāo)準(zhǔn)，但實(shí)際計(jì)算表明，此時(shí)點(diǎn)A和點(diǎn)B處的切向應(yīng)力略大于基于帶圓孔無(wú)限大板理想模型求得的解析解。

圖5 單向受拉的帶圓孔矩形板算例

圖6 有限元網(wǎng)格和數(shù)值流形網(wǎng)格(n =8)

圖7給出了n=8時(shí)，數(shù)值流形方法與精密網(wǎng)格有限元法計(jì)算得到的Von Mises應(yīng)力云圖；圖8給出了點(diǎn)A和點(diǎn)B處切向應(yīng)力隨網(wǎng)格密度變化的對(duì)比?？梢钥闯觯肟拙壧厥鈫卧臄?shù)值流形方法在網(wǎng)格密度較低時(shí)已經(jīng)能夠獲得較準(zhǔn)確的計(jì)算結(jié)果，其計(jì)算精度顯著高于同等網(wǎng)格尺寸下的有限元法，且計(jì)算結(jié)果的收斂速度也有明顯優(yōu)勢(shì)。但數(shù)值流形方法的計(jì)算結(jié)果并不具備有限元法的下限性，其收斂方向也不是單一的。

圖7 Von Mises應(yīng)力云圖對(duì)比(單位:Pa)

圖8 點(diǎn)A和點(diǎn)B切向應(yīng)力隨網(wǎng)格密度變化(算例1)

4.2 四邊受剪的帶圓孔矩形板

本算例中板四邊均受等大、均布的剪切載荷q作用，取其1/4進(jìn)行分析計(jì)算，如圖9所示。圖中邊AD和邊BC施加對(duì)稱位移約束，q=141 MPa，n=2,4,6,…,20。仍以n=80的有限元結(jié)果作為近似精確解。圖10給出了特征點(diǎn)切向應(yīng)力隨網(wǎng)格密度變化，可以看出數(shù)值流形方法仍然具有更高的計(jì)算精度和收斂速度。

圖9 四邊受剪的帶圓孔矩形板算例

圖10 點(diǎn)A和點(diǎn)B切向應(yīng)力隨網(wǎng)格密度變化(算例2)

4.3 一對(duì)邊受線性分布載荷的帶圓孔矩形板

本算例中板的一對(duì)邊受等大和反向的三角形分布載荷，另有一邊固定。如圖11所示，取板的1/2進(jìn)行分析計(jì)算，邊CD為固定邊，邊AF和邊BC施加對(duì)稱位移約束。三角形線性分布載荷在點(diǎn)E處載荷值為0，點(diǎn)D處載荷值為200 MPa。網(wǎng)格數(shù)量按n×2n的方式定義，其中n=2,4,6,…,12，取n=50對(duì)應(yīng)網(wǎng)格尺寸的有限元計(jì)算結(jié)果作為近似精確解。圖12給出了特征點(diǎn)切向應(yīng)力隨網(wǎng)格密度的變化，可以看出數(shù)值流形方法同樣具有更高的計(jì)算精度和收斂速度。

圖11 線性分布載荷帶圓孔矩形板算例

圖12 點(diǎn)A和點(diǎn)B切向應(yīng)力隨網(wǎng)格密度變化(算例3)

5 結(jié) 語(yǔ)

本文針對(duì)有限元法在處理典型開(kāi)孔結(jié)構(gòu)時(shí)的不足之處，基于數(shù)值流形方法的基本框架，提出了一種適用于求解該類問(wèn)題的孔緣單元，這種流形單元具有特殊局部逼近。實(shí)際計(jì)算表明，在網(wǎng)格密度較低的情況下，采用特殊孔緣單元的數(shù)值流形方法的計(jì)算精度相較于有限元法具有明顯優(yōu)勢(shì)，且收斂速度更快。

數(shù)值流形方法是一種較為新穎的數(shù)值模擬方法，目前在結(jié)構(gòu)領(lǐng)域應(yīng)用相對(duì)較少。由于數(shù)值流形方法理論仍在不斷完善，也未有大型通用計(jì)算軟件，因此處理復(fù)雜問(wèn)題仍然有一定困難，但其靈活性強(qiáng)、前處理工作量小以及計(jì)算精度高等突出優(yōu)點(diǎn)均為這種方法帶來(lái)了非常廣闊的工程應(yīng)用前景。

致謝:本文研究開(kāi)展過(guò)程中得到了鄭宏教授的指導(dǎo)和幫助，特此致謝。