李真有,肖映雄

(湘潭大學(xué)土木工程與力學(xué)學(xué)院， 湖南湘潭411105)

?

三維Wilson元及在近不可壓彈性問題中的應(yīng)用

李真有,肖映雄

(湘潭大學(xué)土木工程與力學(xué)學(xué)院， 湖南湘潭411105)

為克服三維近不可壓縮問題的體積閉鎖現(xiàn)象，建立了兩種基于六面體單元的Wilson非協(xié)調(diào)元計算格式，并將其應(yīng)用于兩類含混合邊界條件的近不可壓縮彈性問題的求解。數(shù)值結(jié)果表明:Wilson非協(xié)調(diào)元能有效克服三維體積閉鎖現(xiàn)象，與相同規(guī)模下的協(xié)調(diào)元相比較，它具有更高的計算精度。在三維有限元分析中，剖分網(wǎng)格的質(zhì)量將對計算精度影響很大，實際計算時若能采用各向同性網(wǎng)格，則對問題的分析將具有更好的收斂性。

近不可壓縮問題；體積閉鎖；Wilson非協(xié)調(diào)元；網(wǎng)格質(zhì)量

利用通常的低階協(xié)調(diào)有限元方法處理三維近不可壓縮問題常出現(xiàn)體積閉鎖現(xiàn)象?？朔@種閉鎖現(xiàn)象的方法很多，如混合有限元法[1-3]，非協(xié)調(diào)有限元法[4-6]，高階協(xié)調(diào)元法[7-8]及減縮積分法[9-10]等。減縮積分使單元剛度矩陣秩降低，常引起多余的零能模式。與基于H-R變分原理的混合元格式相比，基于能量泛函極小的離散變分問題更容易被解決，并且有限元方程的系數(shù)矩陣是正定的，將給求解帶來很大方便。對三維問題，若利用高階協(xié)調(diào)元方法，要求單元階次p≥8，這在實際計算中是不可取的。此時，常采用低階非協(xié)調(diào)元(如二次元、線性元)來克服三維近不可壓縮問題的閉鎖現(xiàn)象[5,11]，這種方法具有自由度少、精度高的特點，可解決三維問題計算規(guī)模過大的困難。

在三維有限元分析中，常采用六面體網(wǎng)格剖分，這種網(wǎng)格在計算精度、劃分?jǐn)?shù)量、抗畸變程度等方面比四面體網(wǎng)格具有明顯的優(yōu)勢?；诹骟w單元的Wilson非協(xié)調(diào)元，它通過在單元內(nèi)部設(shè)置附加自由度，以提高完全多項式的次數(shù)、改善計算精度。關(guān)于Wilson元的研究和應(yīng)用，已有很多研究成果[12-15]。本文建立了兩種基于六面體單元的Wilson非協(xié)調(diào)元，并將其應(yīng)用于懸臂梁和Cook膜兩類近不可壓縮彈性問題的求解。數(shù)值結(jié)果表明：Wilson非協(xié)調(diào)元能有效克服三維體積閉鎖現(xiàn)象，與相同規(guī)模下的協(xié)調(diào)元相比較，它具有更高的計算精度。還研究了有限元網(wǎng)格質(zhì)量對Wilson元計算精度的影響，計算中若能采用各向同性網(wǎng)格，則對問題的分析將具有更好的收斂性。

1　近不可壓縮問題和體積閉鎖現(xiàn)象

考慮下述三維線彈性力學(xué)模型:

(1)

其中，Ω?R3，?Ω=Γ0∪Γ1且Γ0∩Γ1=Φ, u為位移向量，f為外力，g為邊界Γ1上的面力，n為邊界Γ1的單位外法線向量，拉梅常數(shù)λ和μ可分別用楊氏模量E和泊松比ν表示為:

(2)

這里，泊松比ν∈[0,0.5)，當(dāng)ν→0.5時，稱相應(yīng)的彈性力學(xué)問題為近不可壓縮問題。

記Hm(Ω)={v|?αv∈L2(Ω),|α|≤m}，則問題(1)對應(yīng)的變分問題可描述為：

(3)

在有限元分析中，網(wǎng)格的劃分及單元的選取是影響其計算精度的重要因素。四面體網(wǎng)格和六面體網(wǎng)格是三維有限元分析中常采用的兩種網(wǎng)格類型，而六面體網(wǎng)格在計算精度、網(wǎng)格數(shù)量、抗畸變程度等方面比四面體網(wǎng)格具有明顯的優(yōu)勢，已成為三維有限元分析中的首選網(wǎng)格。設(shè)Th是區(qū)域Ω上的六面體網(wǎng)格剖分，其中h為Th上所有剖分單元的最大直徑。引入如下p次拉格朗日有限元空間：

(4)

設(shè)變分問題(4)的p次有限元離散化線性系統(tǒng)的矩陣形式為:

(5)

對近不可壓縮問題，通常的低階協(xié)調(diào)元(如三線性元、三二次元)解不再收斂到原問題的解或達(dá)不到最優(yōu)收斂階，即出現(xiàn)所謂的體積閉鎖現(xiàn)象。下面，舉兩個例子來進(jìn)行驗證。

圖1　懸臂梁幾何結(jié)構(gòu)示意圖

nx×ny×nzv=0.3v=0.4v=0.49v=0.499uAzuBzuAzuBzuAzuBzuAzuBz16×16×16-0.2186-0.6157-0.2113-0.6023-0.1584-0.4900-0.0739-0.272424×24×24-0.2370-0.6663-0.2308-0.6545-0.1893-0.5680-0.1070-0.371240×4×4-0.2455-0.6895-0.2376-0.6710-0.1906-0.5545-0.1257-0.396580×8×8-0.2527-0.7087-0.2485-0.6994-0.2244-0.6418-0.1681-0.5008160×16×16-0.2547-0.7140-0.2518-0.7079-0.2410-0.6832-0.2061-0.5964

表2　懸臂梁軸線上A點和B點20節(jié)點二次元解Tab.2　The triquadratic element solutions of two typical points A and B on the cantilever axis

算例2(Cook膜問題)考慮如圖2所示的三維Cook膜問題，其左側(cè)面(即x=0)固定，右側(cè)面上面力g=(0,0,t)T，其中，t=10N/mm2，彈性模量E=103N/mm2，不計自重。采用六面體網(wǎng)格進(jìn)行剖分，分別采用8節(jié)點線三性元和20節(jié)點三二次元進(jìn)行計算，相應(yīng)的數(shù)值結(jié)果如表3和表4所示。

圖2　三維Cook膜問題Fig.2　Three dimensional Cook’s membrane problems表3　結(jié)構(gòu)網(wǎng)格下特征點C在z方向上的8節(jié)點線性元解Tab.3　The trilinear element solutions of on the structured meshes

表4　結(jié)構(gòu)網(wǎng)格下特征點C在z方向上的20節(jié)點二次元解Tab.4　The triquadratic element solutions of on the structured meshes

由上述數(shù)值結(jié)果可知：對于可壓問題，隨著網(wǎng)格規(guī)模的不斷增加，特征點處的8節(jié)點線性元解均能收斂于理論解，而對近不可壓縮問題，當(dāng)ν→0.5時，8節(jié)點線性元解已不再收斂，20節(jié)點二次元解盡管收斂，但達(dá)不到最佳收斂結(jié)果，出現(xiàn)所謂的體積閉鎖現(xiàn)象。另外，計算中所使用網(wǎng)格單元的質(zhì)量也會影響計算結(jié)果的精度，如采用各向同性網(wǎng)格，則具有更好的計算精度。

克服這種閉鎖現(xiàn)象的方法很多，如混合有限元法，非協(xié)調(diào)有限元法，高階協(xié)調(diào)元法及減縮積分法等。對三維問題，常采用低階非協(xié)調(diào)元(如二次元、線性元)來克服閉鎖現(xiàn)象。這種方法具有自由度少、精度高的優(yōu)點，可有效解決三維實際問題計算量過大的困難。

2　Wilson非協(xié)調(diào)元

六面體單元對應(yīng)的位移模式不是完全多項式，而決定有限元解精度的是完全多項式的次數(shù)，非完全的高次項對改善精度不起作用，有時還可能會起相反的作用。Wilson非協(xié)調(diào)元通過在單元內(nèi)部設(shè)置附加自由度，達(dá)到提高完全多項式的次數(shù)、改善計算精度的目的。這種單元的位移模式中包含了剛體位移和常量應(yīng)變，但在相鄰單元邊界上不再保持連續(xù)性。筆者對如圖3所示的局部坐標(biāo)系o-ξηζ下的8節(jié)點和20節(jié)點六面體單元，分別建立Wilson非協(xié)調(diào)元計算格式。

首先，設(shè)8節(jié)點六面體單元的非協(xié)調(diào)元位移模式為：

(6)

類似地，可得到20節(jié)點六面體單元形如式(6)所示的位移模式，其中：

利用最小位能原理，可建立如下形式的單元剛度方程：

(7)

消去內(nèi)部附加自由度αe，可得凝聚后的單元剛度方程為：

Keae=Pe,

(8)

將上述8節(jié)點和20節(jié)點Wilson非協(xié)調(diào)元應(yīng)用于兩類近不可壓縮問題：懸臂梁和Cook膜，以驗證其有效性。

首先，對懸臂梁問題，當(dāng)ν→0.5時，在不同網(wǎng)格規(guī)模和不同網(wǎng)格特性(各向同性網(wǎng)格和各向異性網(wǎng)格)下列出了懸臂梁軸線上A點和B點的兩種Wilson非協(xié)調(diào)元解，總結(jié)如表5和表6所示。

表5懸臂梁軸線上A點和B點8節(jié)點Witson元解

Tab.5The 8-node Wilson element solutions of two typical pointsAandB

ν0.40.450.490.499nx×ny×nzuAzuBzuAzuBzuAzuBzuAzuBz16×16×16-0.2440-0.6935-0.2387-0.6828-0.2333-0.6718-0.2308-0.666724×24×24-0.2478-0.7008-0.2441-0.6932-0.2402-0.6853-0.2389-0.682640×4×4-0.2505-0.7061-0.2479-0.7008-0.2448-0.6946-0.2409-0.686580×8×8-0.2522-0.7093-0.2503-0.7054-0.2483-0.7014-0.2474-0.6994160×16×16-0.2529-0.7107-0.2512-0.7073-0.2496-0.7039-0.2491-0.7030

表6　懸臂梁軸線上A點和B點20節(jié)點Witson元解Tab.6　The 20-node Wilson element solutions of two typical points A and B

表7結(jié)構(gòu)化網(wǎng)格下特征C在z方向上的8節(jié)點WILSON解

ν0.400.490.4950.4990.499516×16×163.954773.936903.932953.911783.9039124×24×243.977763.699743.967953.962543.957828×2×83.886813.829643.811363.725023.6493416×2×163.957153.941213.937593.922513.9083632×4×323.989303.985683.984773.982863.98116

表8　結(jié)構(gòu)化網(wǎng)格下特征點C在z方向上的20節(jié)點Wilson解Tab.8　The 20-node Wilson element solution of on the structured meshes

(a) 336個單元

(b) 11 856個單元

圖4三維Cook膜問題兩種非結(jié)構(gòu)化網(wǎng)格剖分

Fig.4Two unstructured meshes used for Cook’s membrane

表9非結(jié)構(gòu)網(wǎng)格下特征C在z方向上的8節(jié)點Wilson解

表10　非結(jié)構(gòu)網(wǎng)格下特征C在z方向上的20節(jié)點Wilson解Tab.10　The 20-node Wilson element solution of on the structured meshes

3　結(jié)　語

對三維近不可壓問題，若利用通常的低階協(xié)調(diào)元方法求解會出現(xiàn)所謂的體積閉鎖現(xiàn)象。Wilson非協(xié)調(diào)元法是克服三維體積閉鎖現(xiàn)象的一種有效方法，它具有自由度少、精度高的優(yōu)點，可解決三維問題計算規(guī)模大的難題。Wilson非協(xié)調(diào)元對網(wǎng)格質(zhì)量依賴性強(qiáng)，如何設(shè)計適用于任意六面體網(wǎng)格單元的低階非協(xié)調(diào)元方法，將是我們進(jìn)一步研究的問題之一。另外，針對這種非協(xié)調(diào)有限元分析中形成的大型的、稀疏的和高度病態(tài)的正定方程組，如何設(shè)計高效求解算法是一個非常困難的問題，將決定其有限元分析的整體效率，這也是我們今后將進(jìn)一步研究的問題。

[1]CHEUNG Y K, CHEN W J.Hybrid element method for incompressible and nearly incompressible materials[J]. International Journal of Solids and Structures, 1989, 25(5): 483-495.

[2]MORLEY M.A mixed family of elements for linear elasticity[J]. Numerische Mathematik, 1989(55):633-666.

[3]CHAMA A, REDDY B D.New stable mixed finite element approximations for problems in linear elasticity[J]. Computer Methods in Applied Mechanics and Engineering, 2013, 256: 211-223.

[4]WANG L H, QI H.A locking-free scheme of nonconforming rectangular finite element for the planar elasticity[J]. Journal of Computational Mathematics, 2004, 22: 641-650.

[5]FALK R S.Nonconforming finite element methods for the equations of linear elasticity[J]. Mathematics of Computation, 1991, 57: 529- 556.

[6]BRENNER S C.A nonconforming mixed multigrid method for the pure traction problem in planar linear elasticity[J]. Mathematics of Computation, 1994, 63: 435-460.

[7]SCOTT L R, VOGELIUS M.Conforming finite element methods for incompressible and nearly incompressible continua[C]//Large Scale Computations in Fluid Mechanics. Providence RI. American Mathematical Socity, USA: 1985, 22: 221-244.

[8]STENBERG R, SURI M.Mixed h-p finite element methods for problems in elasticity and Stokes flow[J]. Numerische Mathematik, 1996, 72: 367-390.

[9]HUGHES T J R.Equivalence of finite elements for nearly incompressible elasticity[J]. Journal of Applied Mechanics, 1977, 44: 181-183.

[10]MALKUS D S, HUGHES T J R.Mixed finite element methods-reduced and selective integration techniques: a unification of concepts[J]. Computer Methods in Applied Mechanics and Engineering, 1978, 15 (1): 63-81.

[11]QI H, WANG L H, ZHENG W Y.On locking-free finite element schemes for three dimensional elasticity[J]. Journal of Computational Mathematics, 2005, 23:101-112.

[12]鹿曉陽.Wilson 非協(xié)調(diào)元構(gòu)造機(jī)理研究[J]. 工程力學(xué)，2000, 17(5): 58-62.

[13]張春生，龍馭球，須寅.三維內(nèi)參型附加非協(xié)調(diào)位移基本項[J]. 工程力學(xué)，2000, 18(5): 58-62.

[14]陳震，林府標(biāo).三維Wilson元與特征值下界[J]. 南昌大學(xué)學(xué)報(理科版)，2010, 34(1): 19-23.

[15]王芬玲，石東偉，石東洋.擬線性Sobolev方程Wilson元解的超收斂分析及外推[J]. 工程數(shù)學(xué)學(xué)報，2012, 29(5): 720- 724.

(責(zé)任編輯唐漢民梁碧芬)

Wilson element and its application in nearly incompressible elasticity problems in three dimensions

LI Zhen-you, XIAO Ying-xiong

(Civil Engineering and Mechanics College, Xiangtan University, Xiangtan 411105, China)

In order to overcome the volume locking phenomenon of nearly incompressible problems in three dimensions, two types of Wilson nonconforming finite elements have been presented based on the hexahedral elements, and the resulting methods are then applied to the solution of two nearly incompressible problems with mixed boundary conditions.The numerical results have been shown that Wilson elements can effectively overcome the locking phenomenon, and it has higher accuracy compared with the conforming elements under the same mesh size.In three-dimensional finite element analysis, the quality of the mesh will have a great effect on the accuracy, and if the isotropic grids can be used in the practical calculations, the method will have better convergence.

nearly incompressible problem；volume locking；Wilson nonconforming element；mesh quality

2016-04-22；

2016-05-01

國家自然科學(xué)基金資助項目(10972191);湖南省自然科學(xué)基金資助項目(14JJ2063);湖南省教育廳資助科研項目(15A183)

肖映雄(1970—)，男，湖南城步人，湘潭大學(xué)教授；E-mail: xyx610xyx@xtu.edu.cn。

10.13624/j.cnki.issn.1001-7445.2016.1271

O343.2

A

1001-7445(2016)04-1271-08

引文格式：李真有,肖映雄.三維Wilson元及在近不可壓彈性問題中的應(yīng)用[J].廣西大學(xué)學(xué)報(自然科學(xué)版)，2016，41(4)：1271-1278.