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

蒲軍平,汪士應(yīng),陳 軼

(1.浙江工業(yè)大學(xué) 建筑工程學(xué)院,浙江 杭州 310014；2.浙江省工程結(jié)構(gòu)與防災(zāi)減災(zāi)技術(shù)研究重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310014；3. 浙江中浩應(yīng)用工程研究院有限公司，浙江 杭州 310011)

?

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

蒲軍平1，2,汪士應(yīng)1,陳 軼3

(1.浙江工業(yè)大學(xué) 建筑工程學(xué)院,浙江 杭州 310014；2.浙江省工程結(jié)構(gòu)與防災(zāi)減災(zāi)技術(shù)研究重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310014；3. 浙江中浩應(yīng)用工程研究院有限公司，浙江 杭州 310011)

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

曲邊四邊形；有限元網(wǎng)格；位移向量；比例劃分

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

1 直邊八節(jié)點(diǎn)等參單元的網(wǎng)格生成

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

(1)

其中：n=1,3時(shí)，ni=n1,m=1,2,…,n1+1;n=2,4時(shí)，ni=n2,m=1,2,…,n2+1.當(dāng)n1=n2=5時(shí)，等份劃分網(wǎng)格示意圖如圖1所示.

圖1 等份劃分直邊四邊形網(wǎng)格Fig.1 Equal parts straight quadrilateral meshes

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

對(duì)于比例系數(shù)變化下的網(wǎng)格劃分同樣采用圖1的四邊形，為了說(shuō)明方便劃分份數(shù)同樣取為n1=5,n2=5.若需要使角點(diǎn)A11附近的網(wǎng)格加密，使1～4條邊的比例系數(shù)分別為R1=0.8，R2=1，R3=1，R4=1.2.比例系數(shù)是按線(xiàn)段前進(jìn)的方向，前一段線(xiàn)段始末點(diǎn)x坐標(biāo)差比上后一段線(xiàn)段始末點(diǎn)x坐標(biāo)的差，或前一段線(xiàn)段始末點(diǎn)y坐標(biāo)的差比上后一段線(xiàn)段始末點(diǎn)y坐標(biāo)的差的比值.前進(jìn)方向?yàn)锳11→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)可得各個(gè)線(xiàn)段上其他點(diǎn)的坐標(biāo)為

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

(4)然后將對(duì)邊各對(duì)應(yīng)點(diǎn)連線(xiàn)，得到有限元網(wǎng)格劃分圖形如圖2所示.從圖2中可以看出：若劃分的份數(shù)加大，則得到在角點(diǎn)A11附近的網(wǎng)格的加密程度會(huì)更高.

圖2 各邊比例系數(shù)不等時(shí)的網(wǎng)格劃分Fig.2 Mesh generation with unequal ratio coefficients

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

由于曲邊單元的線(xiàn)段非線(xiàn)性，為了得到合理的有限元網(wǎng)格，在直邊單元的基礎(chǔ)上，采用位移向量法對(duì)曲邊單元進(jìn)行劃分.

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

n=1,2,3,4

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

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

當(dāng)圓弧邊是等份劃分時(shí)，此時(shí)對(duì)圓心角進(jìn)行等分即可，由此可得各圓弧邊上的劃分點(diǎn)的坐標(biāo)為

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

(6)

其中r為圓弧邊的半徑.

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

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

(7)

對(duì)于中間各點(diǎn)的坐標(biāo)，需要綜合考慮四條曲邊對(duì)應(yīng)點(diǎn)的位移向量的影響，越靠近中間點(diǎn)的圓弧邊對(duì)該點(diǎn)的影響越大.在虛線(xiàn)中所示單元為直邊單元，以任意一條邊為起點(diǎn)，中間點(diǎn)的坐標(biāo)都是固定的，下面以第一條邊為例，說(shuō)明如何求得中間點(diǎn)的坐標(biāo).

當(dāng)圓弧向內(nèi)凹時(shí)，此時(shí)中間點(diǎn)的坐標(biāo)改變值為所在直線(xiàn)首末點(diǎn)對(duì)應(yīng)的位移向量的坐標(biāo)值乘以一個(gè)折減系數(shù)，方向按照位移向量的方向，折減系數(shù)按比例改變，由式(7)并通過(guò)程序語(yǔ)句累加可得

(8)

(9)

同理，通過(guò)程序語(yǔ)句累加可得各中間點(diǎn)的坐標(biāo)為

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

(10)

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

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

對(duì)于圓弧邊，當(dāng)需要按照比例劃分時(shí)，此時(shí)的比例系數(shù)是按照?qǐng)A心夾角進(jìn)行劃分，圓弧邊上劃分點(diǎn)的坐標(biāo)為

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

(11)

各中間點(diǎn)的坐標(biāo)仍由式(10)得到.

3 有限元模型網(wǎng)格生成及計(jì)算

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

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

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

圖6 開(kāi)圓孔的帶裂縫平面單元的網(wǎng)格生成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

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

4 結(jié) 論

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

[1] 關(guān)振群,宋超,顧元憲.有限元網(wǎng)格生成方法研究的新進(jìn)展[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),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] 蒲軍平,姚宇龍.平面有限元網(wǎng)格實(shí)用劃分方法研究[J].浙江工業(yè)大學(xué)學(xué)報(bào),2012,40(1):88-91.

[4] 曾卓,陳家新.掃掠法有限元網(wǎng)格生成方法[J].計(jì)算機(jī)工程與應(yīng)用,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] 梁國(guó)平.變網(wǎng)格的有限元法[J].計(jì)算數(shù)學(xué),1985(4):377-384.

[7] 董秀山,劉潤(rùn)濤.判斷點(diǎn)與簡(jiǎn)單多邊形位置關(guān)系的新算法[J].計(jì)算機(jī)工程與應(yīng)用,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] 趙春風(fēng),李曉杰,閆鴻浩,等.Ⅰ型裂紋尖端圓弧對(duì)應(yīng)力強(qiáng)度因子影響的數(shù)值研究[J].科學(xué)技術(shù)與工程,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] 陳煒,莊順胥,王飛飛.工程結(jié)構(gòu)多裂紋應(yīng)力強(qiáng)度因子分析[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] 劉寶良,鄒廣平,閆相橋.拉伸載荷作用下單邊缺陷-邊裂紋應(yīng)力強(qiáng)度因子研究[J].計(jì)算力學(xué)學(xué)報(bào),2015,32(1):83-88.

(責(zé)任編輯：劉 巖)

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—),男，新疆烏魯木齊人，教授，博士，主要從事固體力學(xué)和結(jié)構(gòu)工程研究，E-mail:pjp@zjut.edu.cn.

TP391.7

A

1006-4303(2016)05-0529-04