張 吉 陳岱林 王 徽 梁 博

(中國建筑科學(xué)研究院建筑工程軟件研究所,北京 100013)

一種新型有限元網(wǎng)格形狀評價(jià)指標(biāo)

張 吉 陳岱林 王 徽 梁 博

(中國建筑科學(xué)研究院建筑工程軟件研究所,北京 100013)

本文對通用有限元軟件ANSYS與 SAP2000中的單個(gè)網(wǎng)格評價(jià)指標(biāo)進(jìn)行了分析與總結(jié)。在現(xiàn)有單個(gè)網(wǎng)格評價(jià)指標(biāo)計(jì)算的基礎(chǔ)上,建議了一種新的網(wǎng)格評價(jià)指標(biāo),該評價(jià)指標(biāo)不僅能用于單個(gè)網(wǎng)格的評價(jià),而且能用于復(fù)雜的多網(wǎng)格有限元模型。通過 C++編程實(shí)現(xiàn)了整個(gè)評價(jià)指標(biāo)的計(jì)算,對懸臂梁模型采用不同網(wǎng)格劃分后對形狀指標(biāo)進(jìn)行計(jì)算后,發(fā)現(xiàn)按不同網(wǎng)格的形狀指標(biāo)與計(jì)算精確度存在明顯的正相關(guān)性。將其用于多高層建筑結(jié)構(gòu)網(wǎng)格評價(jià),計(jì)算結(jié)果表明該評價(jià)指標(biāo)是高效而又可靠的。

網(wǎng)格;有限單元法;形狀指數(shù)

1 概述

一般說來,單元和網(wǎng)格性能同時(shí)影響有限元數(shù)值計(jì)算結(jié)果的精度。一方面,性能好的單元,在粗糙的或畸形的網(wǎng)格下都能求解出理想的場變量,場變量包括受軸力下的軸向變形、受剪結(jié)構(gòu)中剪應(yīng)變(應(yīng)力)、受彎構(gòu)件中的拉、壓應(yīng)變 (應(yīng)力)等。例如,在經(jīng)典四結(jié)點(diǎn)雙線性位移場 Q4單元基礎(chǔ)上增加含內(nèi)參的附加位移場構(gòu)造的非協(xié)調(diào)位移元Q6單元,在單元形狀規(guī)則或者畸變的情況下,性能都比只具有雙線性位移形函數(shù)的Q4單元的數(shù)值計(jì)算結(jié)果要好[1]。相反,良好的網(wǎng)格劃分能減低對單元性能的要求。不只是精度與網(wǎng)格有關(guān),常用于考核單元性能的小片檢驗(yàn)也與網(wǎng)格有關(guān),有些非協(xié)調(diào)單元在矩形形狀下就能通過小片試驗(yàn),而在非矩形的形狀下就不能通過[2]。

文獻(xiàn)[6,7]中對映射法、基于柵格法等網(wǎng)格劃分方式進(jìn)行了詳細(xì)的介紹,關(guān)于網(wǎng)格劃分更進(jìn)一步的研究是基于幾何形狀的自適應(yīng)網(wǎng)格劃分技術(shù)[3-5],幾何自適應(yīng)網(wǎng)格能識(shí)別所分析區(qū)域的外形,在場變量梯度可能大的區(qū)域自動(dòng)加密網(wǎng)格。但幾何自適應(yīng)網(wǎng)格劃分過程過于復(fù)雜,且與所加荷載形式相關(guān),目前工程中應(yīng)用不多,而廣泛采用的是指定一個(gè)整體網(wǎng)格控制尺寸,對區(qū)域進(jìn)行均勻劃分。不論何種網(wǎng)格方式,都不能保證在一定的控制尺寸下,網(wǎng)格是最優(yōu)的,尤其是針對高層建筑等復(fù)雜結(jié)構(gòu)的網(wǎng)格劃分形狀。因此,本文所研究的問題是,使用相同的單元,在相當(dāng)?shù)那蠼庀南?對網(wǎng)格形狀優(yōu)劣進(jìn)行評價(jià);并得出一個(gè)定量的指標(biāo),以描述單元形狀對數(shù)值計(jì)算結(jié)果的影響程度。

2 單網(wǎng)格形狀評價(jià)指標(biāo)

當(dāng)前流行的通用有限元軟件 ANSYS與SAP2000中都具有強(qiáng)大的自動(dòng)網(wǎng)格劃分功能,并且依據(jù)不同標(biāo)準(zhǔn)給出每個(gè)單元形狀好壞,在ANSYS中還能以圖形或者文本方式查看形狀差單元的位置與指標(biāo)值。

2.1 ANSYS中的網(wǎng)格評價(jià)指標(biāo)

(1)長寬比(aspect ratio)指標(biāo)

如單元具有矩形形狀,直接采取長邊與短邊比值作為長寬比的值,當(dāng)比值大于 20時(shí),給出警告;當(dāng)比值達(dá)到 106時(shí),給出錯(cuò)誤提示。對于任意四邊形(如圖 1),其長寬比的值,通過以下步驟計(jì)算獲得:

圖 1 任意四邊形單元長寬比計(jì)算簡圖

1 找出四邊形四條邊的中點(diǎn)Mp1、Mp2、Mp3和Mp4,分別連接兩對邊中點(diǎn)形成兩條中線Mp1 Mp3,Mp2 Mp4;

④以其中任意的中線為對稱軸線,繪制一矩形I'J'K'L',并且使矩形的各邊均通過原四邊形的中點(diǎn);

㈣以新繪成兩個(gè)矩形的 I'J'K'L'長邊與短邊之比值的較大值作為原四邊形的長寬比。

此法每步計(jì)算的實(shí)現(xiàn)均需要較大的運(yùn)算量;而且,當(dāng)原四邊形為等腰梯形,且梯形高度與兩腰線中點(diǎn)連線長度相等時(shí),計(jì)算的長寬比值為 1,即該指標(biāo)有可能將其誤判為正方形。

(2)平行偏離(parallel deviation)指標(biāo)

用于評價(jià)兩對邊偏離平行線的程度,即偏離矩形的程度。對圖 2中梯形,通過以下方式計(jì)算:

1先求各邊所在直線的單位向量 v1,v2,v3和v4,再求兩對邊單位向量的點(diǎn)積;

④從向量代數(shù)可知,該點(diǎn)積結(jié)果即為兩對邊夾角余弦值,故可找出兩點(diǎn)積中的較大值,通過反余弦求出其夾角以反映網(wǎng)格偏離矩形的程度。當(dāng)夾角為 0°時(shí),兩對邊平行,當(dāng)夾角變大時(shí),其偏離程度越高,性能也會(huì)越差。

對于不含邊中節(jié)點(diǎn)的四節(jié)點(diǎn)單元,ANSYS默認(rèn)設(shè)置為夾角達(dá) 70°時(shí)給出警告,150°時(shí)給出錯(cuò)誤提示。

圖 2 平行偏離指標(biāo)計(jì)算簡圖

(3)雅克比比例法(Jacobian Ratios)

雅克比,即為將任意形狀單元從物理坐標(biāo)中轉(zhuǎn)往基準(zhǔn)坐標(biāo)系的基準(zhǔn)單元時(shí),采用的變換矩陣的行列式。當(dāng)物理坐標(biāo)系中的原始單元為非矩形形狀時(shí),行列式值為一變量,該變量隨著基準(zhǔn)坐標(biāo)值的變化而不斷改變,因而在同一個(gè)原始單元中的四個(gè)角結(jié)點(diǎn)處,雅克比有不同的值。雅克比比例取定為所有四個(gè)結(jié)點(diǎn)中雅克比的最大值與最小值的比值。不難理解,該比值越大,隱含了該坐標(biāo)變換可能帶來越大計(jì)算隱患,數(shù)值計(jì)算結(jié)果的可靠性也就越低。

但是物理單元的形狀為矩形或者平行四邊形時(shí),雅克比為一常數(shù),且等于A/4,A為矩形母單元面積,此時(shí)該法則將失效。

2.2 SAP2000中的網(wǎng)格評價(jià)指標(biāo)

長寬比指標(biāo)

對于任意的四邊形時(shí)網(wǎng)格,建議了兩種計(jì)算方法,第一種計(jì)算方法如圖 3(a)所示。并指出,在關(guān)注分析問題的應(yīng)力時(shí),長寬比不應(yīng)超過小于 1/3,當(dāng)關(guān)注分析問題的位移場為主時(shí),長寬比不應(yīng)小于 1/5。

圖 3 SAP2000中第一類網(wǎng)格評價(jià)指標(biāo)計(jì)算簡圖

第二種計(jì)算方法如圖 3(b)中所示,其中 A1,A2,A3,A4為四邊形對角線剖分成的四個(gè)小三角形的面積。

除要求單元長寬比之外,還要求四邊形的所有內(nèi)角應(yīng)在 45°到 135°之間;然而,SAP2000提供的第一類計(jì)算方式也存在 ANSYS中長寬比計(jì)算時(shí)相同的不足;第二類計(jì)算方法,當(dāng)四邊形單元為矩形時(shí),不論矩形形狀如何狹長,其評價(jià)指標(biāo)保持為一恒定值1。

3 本文建議的網(wǎng)格形狀評價(jià)指數(shù)

采用上面介紹的通用有限元軟件提供的評價(jià)指標(biāo)在某些情況下具有一定的局限性,通常在對某一種形狀網(wǎng)格判定時(shí),能取得較好的效果,而應(yīng)用于另外形狀時(shí),則不甚理想;因此對單元形狀進(jìn)行評價(jià)時(shí),往往需要同時(shí)考慮多個(gè)指標(biāo),如長寬比,大鈍角,尖銳角等的綜合影響;而且,對于建筑結(jié)構(gòu)而言,關(guān)注的問題不是某一個(gè)單元的性態(tài)或某一點(diǎn)的應(yīng)力,而是結(jié)構(gòu)的位移及內(nèi)力等整體層面的物理量,因而能夠?qū)φ麄€(gè)結(jié)構(gòu)網(wǎng)格好壞評價(jià)的指標(biāo)將更有實(shí)用意義?；谏鲜鲈?本文從在學(xué)習(xí)以上一些做法的優(yōu)點(diǎn)基礎(chǔ)上,提出了如下可行的且易于編程實(shí)現(xiàn)的一種網(wǎng)格評價(jià)指標(biāo)。該評價(jià)指標(biāo)既能處理網(wǎng)格全部為四邊形單元情況,也能處理網(wǎng)格同時(shí)包含四邊形與三角形單元的情況:

3.1 四邊形網(wǎng)格的計(jì)算方法

主要計(jì)算步驟如下:

(1)連接任意四邊形的兩對角點(diǎn)形成兩對角線,兩對角線將四邊形分成 4個(gè)小三角形(如圖 4);

(2)分別求出 △AOB、△ BOD、△DOC、△COA周長,并找出最大值,最小值;

(3)以第(2)步中求出的周長最小的三角形的周長與最大周長三角形的周長的比值;

(4)找出對角線交點(diǎn)與四個(gè)角節(jié)點(diǎn)所連線段長度的最小值與最大值,計(jì)算最小值與最大值的比值乘到第(3)步計(jì)算的結(jié)果上,即為該四邊形最終的形狀參考指數(shù)。

為后文陳述方便,不妨命名該形狀指數(shù)為 irreg_ratio,該指數(shù)的計(jì)算式可寫成1m2,其中:

a)為正方形,邊長為 1m;

b)為平行四變形,底邊長 1m,側(cè)邊長m;

c)為矩形,寬 0.5m,高 2m;

d)為等腰矩形,下底寬 1.5m,上底寬 0.5m。

為了初步考證該指數(shù)的性能,即指數(shù)大小與單元形狀好壞之間的相關(guān)性,對圖 5中一組單元的形狀指數(shù)進(jìn)行計(jì)算,所有單元的面積都是相同的,即

初步看來,該指標(biāo)具有以下一些優(yōu)點(diǎn):

a)計(jì)算簡潔,只需簡單的幾何運(yùn)算即能求解;

b)能同時(shí)反映邊長比偏離 1與小銳角與大鈍角的程度,當(dāng)單元形狀為正方形時(shí),即最佳網(wǎng)格時(shí),irreg_ratio=1;

c)從 1到 4網(wǎng)格,rreg_ratio偏離 1越來越大,其值越來越小,單元性能也將越差,與數(shù)值計(jì)算結(jié)論相符[8]的,也說明了該指標(biāo)的合理性;

整個(gè)結(jié)構(gòu)只劃分成一個(gè)單元情況是非常少見的。而在多個(gè)單元組成的分析區(qū)域中,存在極個(gè)別的形狀差的單元對計(jì)算效果影響是比較微小的,作者已做過類似的數(shù)值試驗(yàn)。因此,為了將該形狀指數(shù)應(yīng)用于多網(wǎng)格分析域的網(wǎng)格評價(jià),在單個(gè)單元形狀指數(shù)的計(jì)算方法上再增加以下步驟:

(5)求出每個(gè)單元的面積與整個(gè)分析域的面積;

(6)遍歷每個(gè)單元,求出每個(gè)單元的形狀指數(shù)與該單元的面積之積的總和;

(7)將第(6)步求出的值與整個(gè)分析域的面積之比,作為整個(gè)分析域網(wǎng)格狀況的評價(jià)指標(biāo)。

3.2 三角形網(wǎng)格的計(jì)算方法

網(wǎng)格劃分時(shí),為了降低劃分難度,或者提高劃分成功率,通過會(huì)允許在整體網(wǎng)格中加非常少量的三角形單元。為了適應(yīng)三角形單元形狀的評價(jià),本文建議做法是直接將三角形單元的形狀指數(shù)置為0。當(dāng)圖 6中四邊形角點(diǎn)A不斷與 C點(diǎn)靠近時(shí),單元形狀將趨于三角形。此時(shí)△COA的周長是 4個(gè)三角形周長中最小值,并且接近于 0。由形狀指數(shù)的定義,并且考慮到三角形單元的計(jì)算性能,不難理解可將形狀指數(shù)值取為 0。

基于上述網(wǎng)格指數(shù)的計(jì)算步驟,作者編制了適合于任意網(wǎng)格形狀評價(jià)指標(biāo)計(jì)算的 C++程序。計(jì)算過程中用到的網(wǎng)格信息均來自于ANSYS軟件,在ANSYS中建成幾何模型后,調(diào)用其網(wǎng)格劃分程序自動(dòng)劃分模型,通過簡單后處理操作就可將數(shù)據(jù)文件中網(wǎng)格信息輸出到指定文本文件,后續(xù)過程所需數(shù)據(jù)均從該文本文件中讀取。

3.3 驗(yàn)證實(shí)例

為了進(jìn)一步驗(yàn)證評價(jià)指標(biāo)用于多單元時(shí)的適用性,以MacNeal問題(如圖 7)作為分析對象,探索指標(biāo)值與各種可能的受力條件下的精度關(guān)系。詳細(xì)參數(shù)及模型描述可以參考相關(guān)文獻(xiàn)[8],本文僅做簡單描述:通過改變單元側(cè)邊與水平面得夾角以調(diào)整單元的畸變程度,不論是梯形還是平行四邊形,均使夾角分 15°,30°,45°三個(gè)階次偏離,在自由端施加軸力以及向下剪力兩種荷載,在約束端采取 6自由度全約束的方式,將計(jì)算的結(jié)果與解析解的比值列于表 1中。

從表 1中可以看出,對于受軸力的情形,由于變形簡單,各種單元形狀都基本相同的計(jì)算效果,此時(shí)計(jì)算指數(shù)值與精度之間關(guān)系并不明顯;對于剪彎受力工況情形,形狀指數(shù)與單元的計(jì)算效果具有一致的變化趨勢,隨著單元形狀的畸變程度增加,指數(shù)變得越小。計(jì)算性能較差的梯形單元,其形狀指數(shù)相對于相同畸變角度的平行四邊形要小,因而驗(yàn)證了本文形狀指數(shù)的適用性。

表1 MacNeal問題采用不同形狀單元計(jì)算的形狀指數(shù)與求解精度關(guān)系

4 結(jié)論

本文對結(jié)構(gòu)有限元分析中一個(gè)重要的影響因素——網(wǎng)格剖分方案進(jìn)行一些探討,通過對比試驗(yàn)及數(shù)值算例分析,得出了如下一些結(jié)論:本文建議的網(wǎng)格形狀評價(jià)的指數(shù),計(jì)算過程比較簡便,既適用于單個(gè)網(wǎng)格,也能適用于多網(wǎng)格區(qū)域的網(wǎng)格性能評價(jià)。將該指標(biāo)推廣到多高層建筑結(jié)構(gòu)網(wǎng)格評價(jià)后,計(jì)算表明該指標(biāo)仍然高效可靠。

[1]王勖成,邵敏.有限單元法基本原理和數(shù)值方法[M].北京:清華大學(xué)出版社,1997.

[2]ROBERT D Cook,DAVID SMalkus,MI CHAEL E Plesha,et a.l Concepts and applications of finite ele ment analysis[M].4th ed.Ne w York:JohnWiley&Sons,2001.

[3]黃曉東,杜群貴,葉邦彥.二維幾何特征自適應(yīng)有限元網(wǎng)格生成(一)——幾何特征識(shí)別[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2004(7):923-928.

[4]杜群貴,劉邵宏,黃曉東.二維幾何特征自適應(yīng)有限元網(wǎng)格生成(二)——算法描述[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2004(7):929-934.

[5]傅沛福,吳淑芳,胡平等.全自動(dòng)自適應(yīng)網(wǎng)格細(xì)化[J].計(jì)算力學(xué)學(xué)報(bào),1997(2):114-118.

[6]K Hole.Finite elementmesh generation methods:a review and classification[J].Computer A ided Design, 1988(20):27-38.

[7]關(guān)振群,宋超等.有限元網(wǎng)格生成方法研究的新進(jìn)展[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)報(bào),2003(1):1-14.

[8]R.H.Mac Nea,l R.L.Harder,“A Proposed Standard Set of Proble ms to Test Finite Ele ment Accuracy”,Proceedings,25th SDMFinite Ele mentValidation Forum,1984.

A New Index for Finite ElementMesh Shape Evaluate

Zhang Ji,Chen Dailin,Wang Hui,Liang Bo

(Institution of Building Engineering Sof tware,ChinaA cademy of Building Research,Beijing100013,China)

Single ele m entmesh shape evaluate index ofANSYS and SAP2000 is analyzed and summarized in the present paper.Based on this indexes,a new index is proposed.Benefit from the new index,not only single elementmesh but alsomesh compose ofmore elements can be evaluated.A cantileverbea m problem ism eshed by di-f ferent shape and procedure programmed by C++ language is used to calculate this index of themesh,after compared the solution w ith exact one,a conclusion can be draw that there is a strong correlation be tween the index value and the precision of solution.Multilayer and high-rise buildingsmesh is also evaluated,the results show the index is still efficient and reliable.

Mesh;F inite E le m entMethod;Shape Index

TU375.1

A

1674-7461(2011)02-0016-05

張吉(1982-),男,博士研究生。主要從事建筑結(jié)構(gòu)設(shè)計(jì)軟件的研發(fā)及有限元數(shù)值分析工作。