平面八節(jié)點有限元網格生成的位移向量法

蒲軍平,汪士應,陳 軼

(1.浙江工業(yè)大學 建筑工程學院,浙江 杭州 310014；2.浙江省工程結構與防災減災技術研究重點實驗室，浙江 杭州 310014；3. 浙江中浩應用工程研究院有限公司，浙江 杭州 310011)

?

平面八節(jié)點有限元網格生成的位移向量法

蒲軍平1，2,汪士應1,陳 軼3

(1.浙江工業(yè)大學 建筑工程學院,浙江 杭州 310014；2.浙江省工程結構與防災減災技術研究重點實驗室，浙江 杭州 310014；3. 浙江中浩應用工程研究院有限公司，浙江 杭州 310011)

在曲邊圖形的有限元計算中，為了保證計算精度及減小計算規(guī)模，提出了一種生成八節(jié)點四邊形單元的位移向量法，以比例漸變的方式綜合考慮了四邊形各曲邊的格柵點對中間各對應點的影響.為了避免出現奇異性單元，在劃分的過程中靈活地應用了比例劃分.利用開發(fā)的VB和FORTRAN程序對一些模型進行了前處理網格劃分和有限元數值計算，結果表明：該方法能簡單、快速地生成有限元網格，并且數值結果與解析解良好吻合.

曲邊四邊形；有限元網格；位移向量；比例劃分

有限元網格生成是數值計算的前提，從簡單的等份劃分到自適應網格的劃分技術已經能解決各種形狀的前處理[1-5].從有限元的前處理到后處理，網格劃分的好壞直接決定了計算的速度和精度.筆者提出了一種生成平面八節(jié)點四邊形單元的位移向量法，將網格生成分為直邊單元和曲邊單元，為了避免出現奇異性單元在劃分的過程中靈活地應用了比例劃分[6]，對關鍵性的局部區(qū)域網格可靈活地進行加密.對于曲邊單元在使用比例劃分的基礎上采用位移向量法進行劃分，該方法使用簡單，通過對一些模型網格劃分和計算，精度滿足實際要求.

1 直邊八節(jié)點等參單元的網格生成

如圖1所示的直邊四邊形，四條邊的編號分別為1，2，3，4，四個角點分別為A11(x11，y11)，A21(x21，y21)，A31(x31，y31)，A41(x41，y41)，其中1，3邊劃分份數為n1，2，4邊的劃分份數為n2.各邊上劃分點的坐標按下式求得

(1)

其中：n=1,3時，ni=n1,m=1,2,…,n1+1;n=2,4時，ni=n2,m=1,2,…,n2+1.當n1=n2=5時，等份劃分網格示意圖如圖1所示.

圖1 等份劃分直邊四邊形網格Fig.1 Equal parts straight quadrilateral meshes

在對一些圖形進行網格劃分時，一些關鍵的局部區(qū)域需要網格加密，而大部分區(qū)域的網格可劃分的較為稀疏，這時需要對每條邊加入比例系數，就可以降低有限元網格的數目，在保證有限元計算精度的前提下提高計算效率.

對于比例系數變化下的網格劃分同樣采用圖1的四邊形，為了說明方便劃分份數同樣取為n1=5,n2=5.若需要使角點A11附近的網格加密，使1～4條邊的比例系數分別為R1=0.8，R2=1，R3=1，R4=1.2.比例系數是按線段前進的方向，前一段線段始末點x坐標差比上后一段線段始末點x坐標的差，或前一段線段始末點y坐標的差比上后一段線段始末點y坐標的差的比值.前進方向為A11→A21→A31→A41，即

Rm=(ymn-ym(n-1))/(ym(n-1)-ym(n-2))=

(xmn-xm(n-1))/(xm(n-1)-xm(n-2))

m=1,2,3,4;n=3,4,…,n1+1

m=1,2,…,n1+1;l=1,2,…,m

(3)

由式(1,3)可得各個線段上其他點的坐標為

n=1,2,3,4;m=2,3,…,n1+1

(4)然后將對邊各對應點連線，得到有限元網格劃分圖形如圖2所示.從圖2中可以看出：若劃分的份數加大，則得到在角點A11附近的網格的加密程度會更高.

圖2 各邊比例系數不等時的網格劃分Fig.2 Mesh generation with unequal ratio coefficients

2 圓弧曲邊八節(jié)點等參單元的網格生成

由于曲邊單元的線段非線性，為了得到合理的有限元網格，在直邊單元的基礎上，采用位移向量法對曲邊單元進行劃分.

對于如圖3所示的以圓弧邊為主的曲邊單元，各邊編號分別為1，2，3，4，四個角點坐標分別為A11(x11，y11)，A21(x21，y21)，A31(x31，y31)，A41(x41，y41)，這里為了說明方便同樣取1，3邊的劃分份數為n1和2，4邊的劃分份數n2均為5，各邊的圓心坐標分別為O1(x1,y1)，O2(x2,y2)，O3(x3,y3)，O4(x4,y4)，各種角度計算式為

n=1,2,3,4

(5)其中：αn為各邊起始點與圓心的連線與x軸正向的夾角；βn為各邊末點與圓心的連線與x軸正向的夾角；θn為各邊起始點與末點直角的圓心角的大小.因此，對于圓弧邊按比例劃分是對所夾圓心角按比例劃分.

圖3 曲邊單元等比例劃分網格圖Fig.3 Dividing grid graph with the same ratio in the curved edge element

當圓弧邊是等份劃分時，此時對圓心角進行等分即可，由此可得各圓弧邊上的劃分點的坐標為

n=1,2,3,4;m=1,2,…,n1-1

(6)

其中r為圓弧邊的半徑.

對于中間點的求取，采用位移向量法.對于每個圓弧邊，可以先分別將各圓弧邊的首尾進行連接，這樣就得到分別以圓弧各角點相連的一個直邊單元，如圖3中虛線所示，對于靠近每條圓弧邊的第一條直邊，可以認為是將此直線進行拉伸將直線拉到圓弧邊所在的位置，這樣以后每條邊都可以認為會受到次圓弧邊的影響，當隨著離圓弧邊越來越遠時，認為受到的影響越來越小，這樣綜合起來看，中間網格上的點會受到四邊形四條邊上對應點的影響，當對四邊形四條邊循環(huán)一次，可以得到個中間點的最終坐標.

n=1,2,3,4;m=2,3,…,n1

(7)

對于中間各點的坐標，需要綜合考慮四條曲邊對應點的位移向量的影響，越靠近中間點的圓弧邊對該點的影響越大.在虛線中所示單元為直邊單元，以任意一條邊為起點，中間點的坐標都是固定的，下面以第一條邊為例，說明如何求得中間點的坐標.

當圓弧向內凹時，此時中間點的坐標改變值為所在直線首末點對應的位移向量的坐標值乘以一個折減系數，方向按照位移向量的方向，折減系數按比例改變，由式(7)并通過程序語句累加可得

(8)

(9)

同理，通過程序語句累加可得各中間點的坐標為

m=2,3，…,n1;n=1,2，…,n2-1

(10)

按照式(10)得出圖3的最終網格劃分圖，如圖4所示.

圖4 圓弧曲邊單元網格劃分Fig.4 The meshing of circular arc curved edge elements

對于圓弧邊，當需要按照比例劃分時，此時的比例系數是按照圓心夾角進行劃分，圓弧邊上劃分點的坐標為

n=1,2,3,4;m=1,2,…,n1-1

(11)

各中間點的坐標仍由式(10)得到.

3 有限元模型網格生成及計算

按以上給出的方法，對一些具體案例進行了網格劃分，利用開發(fā)的VB系統(tǒng)軟件對下述模型進行了前處理網格劃分并應用開發(fā)的FORTRAN程序進行了數值計算.

算例1 帶有裂縫平面板邊長分別為2H和2W，裂縫半寬a=100 mm，H/a=L/a=10，豎向拉應力σ0=100 N/mm2，泊松比v=0.3，彈性模量E=2×105N/mm2.生成平面八節(jié)點等參元有限元計算模型如圖5所示，計算所得裂尖處應力強度因子[8-12]為KI/KI0=0.78.

圖5 帶裂縫平面單元網格劃分圖Fig.5 The plane mesh graph with a crack

圖6 開圓孔的帶裂縫平面單元的網格生成Fig.6 Mesh generation with cracks in a plane element with four circular holes

圖7 裂縫處的局部放大圖Fig.7 Local magnification in the location of a crack

通過對計算結果與解析解進行比較，可看出計算精度滿足實際要求.算例1中單元數取為864，節(jié)點數為2 684，算例2單元數取1 276，節(jié)點數為3 948，算例2較算例1雖然圖形更復雜，但由于采用第2節(jié)中所述的方法，其規(guī)模的增加是可控的，而且這種方法對于更復雜圖形處理的效果將更為明顯.

4 結 論

曲邊平面圖形的有限元計算中，通過對圓弧兩邊的圓心角進行比例劃分來對圓弧邊劃分，可以得到合理的有限元網格.在給出的兩個算例模型中，將板上下邊緣處的平均荷載轉化為集中荷載加在了網格的角點上，由于計算中采用的是八節(jié)點等參元，在每兩個角點之間還有一個中間點，因此，當荷載以每個中間點進行平均時，荷載將更接近均布荷載，計算的精度也將會更高.筆者給出的方法是針對圓弧曲邊進行的網格劃分，從該方法中可以看出：當曲邊不是圓弧邊時，只要該曲邊可以用函數進行模擬，按照此方法在曲邊前進的方向上，根據曲線的性質靈活地劃分格柵點，同樣可以得到合理的有限元網格.

[1] 關振群,宋超,顧元憲.有限元網格生成方法研究的新進展[J].計算機輔助設計與圖形學學報,2003,15(1):1-14.

[2] HU Jiangchun, WANG Hongfang, HE Manchao. Research finite element method mesh generation based on material-discontinuous-nature[J].Advanced materials research,2011,156/157:74-78

[3] 蒲軍平,姚宇龍.平面有限元網格實用劃分方法研究[J].浙江工業(yè)大學學報,2012,40(1):88-91.

[4] 曾卓,陳家新.掃掠法有限元網格生成方法[J].計算機工程與應用,2013,49(2):219-221.[5] KOKO J. A matlab mesh generator for the two-dimensional finite element method[J].Applied mathematics and computation,2014,250:650-664.

[6] 梁國平.變網格的有限元法[J].計算數學,1985(4):377-384.

[7] 董秀山,劉潤濤.判斷點與簡單多邊形位置關系的新算法[J].計算機工程與應用,2009,45(2):185-186.

[8] MAWATARI T,NELSON D V.Method for efficient computation of stress intensity factors from weight functions by singular point elimination[J].Engineering fracture mechanics,2011,78(16):2713-2730.

[9] CAI Gangyi,ZHANG Shixing.Finite element analysis and calculation in cracks of stress intensity factor of thick walled cylinder[J].Advanced mechanical engineering,2010,26/27/28:603-607.

[10] 趙春風,李曉杰,閆鴻浩,等.Ⅰ型裂紋尖端圓弧對應力強度因子影響的數值研究[J].科學技術與工程,2010,10(4):961-965.

[11] MA Youli.Mode I and mode II stress intensity factors based on crack opening and sliding displacements[J].Materials science and engineering,2011,179/180:1417-1422.

[12] 陳煒,莊順胥,王飛飛.工程結構多裂紋應力強度因子分析[J].科技與創(chuàng)新,2014(19):24.

[13] HOSSEINI-TEHRANI P,HOSSEINI-GODARZI A R,TAVANGAR M.Boundary element analysis of stress intensity factor K-I in some two-dimensional dynamic thermoelastic problems[J].Engineering analysis with boundary elements,2005,29(3):232-240.

[14] PU Junping. Effect of defective holes on crack using the boundary element method[J].Engineering analysis with boundary elements,2006,7(30):577-581.

[15] 劉寶良,鄒廣平,閆相橋.拉伸載荷作用下單邊缺陷-邊裂紋應力強度因子研究[J].計算力學學報,2015,32(1):83-88.

(責任編輯：劉 巖)

A displacement vector method for generating planar eight-node finite element mesh

PU Junping1,2, WANG Shiying1, CHEN Yi3

(1.College of Civil Engineering and Architecture, Zhejiang University of Technology, Hangzhou 310014, China; 2.Zhejiang Key Laboratory of Civil Engineering Structures & Disaster Prevention and Mitigation Technology, Hangzhou 310014, China;3.Zhejiang Zhonghao Applied Engineering Technology Institute Co., Ltd., Hangzhou 310011, China)

In the finite element calculation for a geometry with curved edges, a displacement vector method is proposed for generating an eight-node quadrilateral mesh to ensure the computational accuracy and to reduce the amount of computation. In this method, the effect of the grid points on each curved edge on the middle ones on the edge is considered with a gradually changing proportion. To avoid any singular element, the proportional division is flexibly used in the process of mesh generation. With the developed VB and FORTRAN programs, pre-processing mesh generations and finite element computations are conducted for some models. The results show that the finite element mesh can be simply and rapidly generated with this method and that the numerical results are in good agreement with the analytical solution.

curved quadrilateral; finite element mesh; displacement vector; proportional division

2016-01-18

蒲軍平(1962—),男，新疆烏魯木齊人，教授，博士，主要從事固體力學和結構工程研究，E-mail:pjp@zjut.edu.cn.

TP391.7

A

1006-4303(2016)05-0529-04