一種新型有限元網格形狀評價指標

張 吉 陳岱林 王 徽 梁 博

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

一種新型有限元網格形狀評價指標

張 吉 陳岱林 王 徽 梁 博

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

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

網格;有限單元法;形狀指數

1 概述

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

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

2 單網格形狀評價指標

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

2.1 ANSYS中的網格評價指標

(1)長寬比(aspect ratio)指標

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

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

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

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

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

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

(2)平行偏離(parallel deviation)指標

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

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

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

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

圖 2 平行偏離指標計算簡圖

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

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

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

2.2 SAP2000中的網格評價指標

長寬比指標

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

圖 3 SAP2000中第一類網格評價指標計算簡圖

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

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

3 本文建議的網格形狀評價指數

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

3.1 四邊形網格的計算方法

主要計算步驟如下:

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

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

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

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

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

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

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

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

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

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

初步看來,該指標具有以下一些優點:

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

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

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

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

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

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

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

3.2 三角形網格的計算方法

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

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

3.3 驗證實例

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

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

表1 MacNeal問題采用不同形狀單元計算的形狀指數與求解精度關系

4 結論

本文對結構有限元分析中一個重要的影響因素——網格剖分方案進行一些探討,通過對比試驗及數值算例分析,得出了如下一些結論:本文建議的網格形狀評價的指數,計算過程比較簡便,既適用于單個網格,也能適用于多網格區域的網格性能評價。將該指標推廣到多高層建筑結構網格評價后,計算表明該指標仍然高效可靠。

[1]王勖成,邵敏.有限單元法基本原理和數值方法[M].北京:清華大學出版社,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]黃曉東,杜群貴,葉邦彥.二維幾何特征自適應有限元網格生成(一)——幾何特征識別[J].計算機輔助設計與圖形學學報,2004(7):923-928.

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

[5]傅沛福,吳淑芳,胡平等.全自動自適應網格細化[J].計算力學學報,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]關振群,宋超等.有限元網格生成方法研究的新進展[J].計算機輔助設計與圖形學報,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-),男,博士研究生。主要從事建筑結構設計軟件的研發及有限元數值分析工作。