張 吉 陳岱林 王 徽 梁 博

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

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

張 吉 陳岱林 王 徽 梁 博

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

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

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

1 概述

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

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

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

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

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

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

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

圖 1 任意四邊形單元長寬比計算簡圖

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

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

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

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

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

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

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

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

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

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

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

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

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

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

長寬比指標(biāo)

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

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

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

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

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

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

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

主要計算步驟如下:

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

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

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

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

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

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

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

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

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

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

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

a)計算簡潔,只需簡單的幾何運算即能求解;

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

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

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

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

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

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

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

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

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

3.3 驗證實例

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

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

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

4 結(jié)論

本文對結(jié)構(gòu)有限元分析中一個重要的影響因素——網(wǎng)格剖分方案進行一些探討,通過對比試驗及數(shù)值算例分析,得出了如下一些結(jié)論:本文建議的網(wǎng)格形狀評價的指數(shù),計算過程比較簡便,既適用于單個網(wǎng)格,也能適用于多網(wǎng)格區(qū)域的網(wǎng)格性能評價。將該指標(biāo)推廣到多高層建筑結(jié)構(gòu)網(wǎng)格評價后,計算表明該指標(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)格生成(一)——幾何特征識別[J].計算機輔助設(shè)計與圖形學(xué)學(xué)報,2004(7):923-928.

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

[5]傅沛福,吳淑芳,胡平等.全自動自適應(yīng)網(wǎng)格細化[J].計算力學(xué)學(xué)報,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].計算機輔助設(shè)計與圖形學(xué)報,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è)計軟件的研發(fā)及有限元數(shù)值分析工作。