基于數值流形方法的特殊孔緣單元構造及開孔板求解

武卓威， 劉 俊*

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

1 引 言

開孔在各類工程結構物中較為常見，如船體結構內部通常設有人孔、流水孔、透氣孔和減輕孔等。這些開孔的存在會破壞板殼結構的完整性，降低臨近區域結構的強度，同時可能在孔緣產生局部應力集中。

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

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

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

2 數值流形方法基本原理

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

圖1 數值流形方法覆蓋系統

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

圖2 標準網格下數值流形方法的覆蓋系統

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

(1)

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

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

(2)

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

(3)

(4)

式中Ω為求解域，n為物理片總數。在實際應用中，權函數常基于有限元單元形函數構造。由此，平面問題數值流形方法的整體逼近表達為

(5)

獲得了問題域的整體逼近后，即可根據彈性力學相關理論，結合平衡方程、幾何方程和本構方程，利用整體逼近推導應力和應變，隨后依據最小勢能原理建立系統的控制方程，并引入邊界條件，該過程與有限元方法基本一致。

3 平面圓孔特殊孔緣單元的構造

3.1 孔緣局部逼近的構造

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

(6)

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

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

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

(7)

(8)

(9)

式中σr為徑向應力，σθ為切向或環向應力，τr θ為剪切應力，R為圓孔半徑。

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

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

(10)

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

(11)

(12)

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

(13)

3.2 流形單元類型的劃分

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

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

圖4 三類流形單元的劃分

3.3 權函數的構造

標準網格下，各數學片的權函數可按式(14,15)的形式，取為包含于該數學片的正方形網格在張成該數學片頂點處的有限元形函數構成的分塊函數。

(14)

(15)

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

(16)

式中Ni由權函數及基矩陣得到，其作用類似于有限元法的形函數矩陣

因此有

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

(17)

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

3.4 數值流形方法計算流程

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

(18)

(19)

σ=Dε

(20)

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

ε=LNd=Bd,σ=DBd=Sd

(21)

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

Kd=q

(22)

式中K為總體剛度矩陣，q為總體載荷列陣，分別由流形單元的單元剛度矩陣Ke和單元載荷列陣qe按編號順序組裝而成，其中

(23)

(24)

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

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

4 算例驗證

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

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

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

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

圖6 有限元網格和數值流形網格(n =8)

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

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

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

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

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

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

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

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

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

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

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

5 結 語

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

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

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