基于數(shù)值流形方法的特殊孔緣單元構(gòu)造及開孔板求解

武卓威， 劉 俊*

(1.上海交通大學 海洋工程國家重點實驗室，上海 200240；2.上海交通大學 高新船舶與深海開發(fā)裝備協(xié)同中心，上海 200240)

1 引 言

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

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

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

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

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

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

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

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

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

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

(1)

每個物理片Pi均具有對應(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ù)。在實際應(yīng)用中，權(quán)函數(shù)?；谟邢拊獑卧魏瘮?shù)構(gòu)造。由此，平面問題數(shù)值流形方法的整體逼近表達為

(5)

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

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

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

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

(6)

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

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

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

(7)

(8)

(9)

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

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

由此，本文取[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為自由度向量，每一個特殊局部逼近的自由度數(shù)為10。

(13)

3.2 流形單元類型的劃分

孔緣特殊物理片的基向量項數(shù)增加，增加了求解計算量?？紤]到應(yīng)力集中是一種局部效應(yīng)，僅需對孔緣附近的少數(shù)物理片進行擴展，對于距離孔緣較遠的流形單元，則按照一般平面問題進行處理。

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

圖4 三類流形單元的劃分

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

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

(14)

(15)

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

(16)

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

因此有

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

(17)

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

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

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

(18)

(19)

σ=Dε

(20)

式中D為平面應(yīng)力問題的彈性矩陣。將式(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為總體載荷列陣，分別由流形單元的單元剛度矩陣Ke和單元載荷列陣qe按編號順序組裝而成，其中

(23)

(24)

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

將單元剛度矩陣組裝后求解方程式(22)，即可得到各局部逼近的精確表達式，進而得到問題域上的整體逼近。

4 算例驗證

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

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

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

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

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

圖7給出了n=8時，數(shù)值流形方法與精密網(wǎng)格有限元法計算得到的Von Mises應(yīng)力云圖；圖8給出了點A和點B處切向應(yīng)力隨網(wǎng)格密度變化的對比。可以看出，引入孔緣特殊單元的數(shù)值流形方法在網(wǎng)格密度較低時已經(jīng)能夠獲得較準確的計算結(jié)果，其計算精度顯著高于同等網(wǎng)格尺寸下的有限元法，且計算結(jié)果的收斂速度也有明顯優(yōu)勢。但數(shù)值流形方法的計算結(jié)果并不具備有限元法的下限性，其收斂方向也不是單一的。

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

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

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

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

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

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

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

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

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

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

5 結(jié) 語

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

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

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