基于面積坐標的四邊形雜交有限元法

高何金雨 張世全

摘 要 ???： 本文針對二維線彈性問題提出了一種基于面積坐標的新型雜交應力四邊形有限元法AGQ- LQ ?6. 該方法基于廣義Hellinger-Reissner變分原理，位移逼近采用含內部位移的四節點廣義協調元，應力逼近則采用九參數線性應力模式. 數值算例表明，本文構造的有限元既能保持面積坐標廣義協調元對網格畸變不敏感及粗網格精度較高的優點，又能有效克服泊松閉鎖現象.

關鍵詞 ：線彈性問題； 四邊形面積坐標方法； 雜交應力有限元； 泊松閉鎖現象

中圖分類號 ： O241.82 文獻標識碼 ：A DOI ： ?10.19907/j.0490-6756.2023.041002

A hybrid stress quadrilateral finite element method based on area coordinates

GAO He-Jin-Yu， ?ZHANG Shi-Quan

（School of Mathematics， Sichuan University， Chengdu 610064， China）

In this paper， we propose a new hybrid stress quadrilateral finite element method AGQ- LQ ?6 for the two-dimensional linear elasticity problems by using the area coordinates. In this method， a 4-node generalized conforming element with internal displacements for the displacement approximation and a 9-parameter linear stress mode for the stress approximation are adopted based on the generalized Hellinger-Reissner variational principle. Numerical experiments show that the method is insensitive to mesh distortion. Meanwhile， it is also shown that this method can keep high precision even under coarse mesh and avoid the so-called Poisson-locking phenomenon effectively.

Linear elasticity problem; Quadrilateral area coordinates; Hybrid stress finite element; Poisson-locking phenomenon

（2010 MSC 65M60）

1 引 言

上世紀六十年代以來，如何針對平面彈性問題構造具有較高精度的低階四邊形有限元成為工程計算領域的一個研究熱點. 其中，基于等參坐標的四邊形元法一直有廣泛的應用，其中最為典型的代表是4節點等參雙線性元. 然而，在一些標準的考核算例中，這種最低階的四邊形位移元的數值精度較低，對網格畸變敏感，并且會發生泊松閉鎖現象. 因此，在不增大計算規模的前提下，如何提高低階四邊形元的性能一直廣受關注.

1973年，Wilson等 ?［1］在等參雙線性元基礎上引入非協調內部位移，構造了計算精度更高的Wilson非協調元. 之后，Wilson等 ?［2］又對Wilson元進行了改進，使其在泊松閉鎖測試算例中也能算出一致收斂的結果. 此外，基于廣義變分原理的雜交應力有限元框架也是改進四邊形位移元的性能的有效途徑， 其中的H-R變分原理 ?［3，4］因為將位移和應力分別作為獨立變量的特點而成為研究者們推導混合元/雜交元的一種普遍方法. 之后，在各類不同類型的微分方程有限元求解中，不同的改進方法，如質量集中雜交應力有限元法 ?［5］，組合雜交有限元法 ?［6］，弱Galerkin有限元法 ?［7］，增強混合有限元法 ?［8］等都有合理的運用，并且這些方法都能在不增加計算規模的情況下達到提高數值精度的效果.

1999年，Long等 ?［9-11］提出了采用四邊形面積坐標代替傳統的等參坐標系來構造廣義協調四邊形單元的方法，稱為“四邊形面積坐標（QAC）方法”. 與之前的兩類改進方法不同，QAC方法以坐標變換為切入點，替換掉常用的等參坐標系，轉而使用一種新型面積坐標，其與笛卡爾坐標之間的變換關系是線性的，可相互顯式表出.我們知道，對于等參雙線性元而言，一方面，在等參變換中笛卡爾坐標可以由等參坐標的雙線性形式表示出來，而等參坐標卻無法用笛卡爾坐標的有限項表示出來; 另一方面，其單元剛度矩陣積分式中含有雅可比逆的行列式，只能通過數值積分近似計算，而QAC方法卻可以克服上述缺點.

本文基于雜交應力有限元框架，采用Zhou和Xie在文獻［12，13］中提出的應力模式對四邊形面積坐標下的廣義協調元進行改進，進而提出了一種基于面積坐標的雜交應力元AGQ- LQ ?6. 在該方法中，單元的位移逼近我們采用含內部位移的四節點廣義協調元AGQ6，應力逼近則采用九參數線性應力模式. 最后，我們以標準的數值算例驗證了該方法的數值性能.

后文的安排如下. 在第2節中我們給出四邊形面積坐標的定義及面積元的構造. 在第3節中，我們基于雜交應力有限元框架給出改進的面積元的離散格式. 我們在第4節中給出數值算例. 第5節總結了本文所得結果.

2 四邊形面積坐標

2.1 四邊形的特征參數

四邊形單元面積坐標的構造基于三角形單元面積坐標構造的思想. 我們首先介紹四邊形面積坐標需要用到的形狀特征參數.如圖1所示， A ?是四邊形單元的面積，1，2，3，4 為單元的四個頂點，以圖中的順序標號， ?A′ 和 ?A″ 分別是三角形Δ124和Δ123的面積.參數的定義式為

g ?1= ?A′ A ， g ?2= ?A″ A ， g ?3=1- g ?1， g ?4=1- g ?2

（1）

由式（1） 可以看到，參數 ?g ?3 和 ?g ?4 能分別用 ?g ?1， g ?2 表示，且參數 ?g ?1， g ?2 在一個凸四邊形中的取值范圍是 0< g ?1<1， 0< g ?1<1 .

注1 ??我們稱 ?g ?1 和 ?g ?2 為形狀特征參數. 當 ?g ?1， g ?2 取一些特殊值時，圖1中的單元會轉變為一些有特定形狀的單元：

（1） 當 ?g ?1= g ?2= 1 2 ??時，四邊形單元就會退化為一個平行四邊形；

（2） 當 ?（g ?1 -g ?2） （1-g ?1 -g ?2）=0 時，四邊形單元就會退化為梯形；

（3） 當 ?g ?1 g ?2 （1-g ?1）（ 1-g ?2）=0 時，四邊形單元就會退化為三角形.

2.2 四邊形面積坐標的定義

如圖2所示，令 P ?為四邊形單元內的任意一點，則 P 的面積坐標 （ L ?1， L ?2， L ?3， L ?4） 定義為 ?［10］

L ?i= ?A ?i A ，i=1，2，3，4 ?（2）

其中 A 是該四邊形單元的面積， ?A ?i （ i =1，2，3，4）是由點 P 分別和單元其中兩個相鄰頂點構成的三角形的面積，單元四個頂點的坐標表示為：

節點1： ?（ g ?2，1- g ?2，0，0）;

節點2： ?（ 0，1-g ?1， g ?1，0）;

節點3： ?（0，0，1- g ?2， g ?2） ;

節點4： ?（1- g ?1，0，0， g ?1） .

四邊形面積坐標也可以用笛卡爾坐標表示. 若圖2中單元的四個頂點的笛卡爾坐標依次為 （ x ?1， y ?1） ， （ x ?2， y ?2） ， （ x ?3， y ?3） ， （ x ?4， y ?4） ，那么4個三角形的面積 ?A ?i （ i =1，2，3，4）可以通過行列式計算得到

A ?i= 1 2 ??1 x y 1 ?x ?j ?y ?j 1 ?x ?k ?x ?k ??，

i，j，k=1，2，3 i，j，k=2，3，4 i，j，k=3，4，1 i，j，k=4，1，2 ??.

再將 ?A ?i 代入式（2）可得

L ?i= 1 2A ??a ?i+ b ?ix+ c ?iy ，i=1，2，3，4 ??（3）

其中 A 表示該四角形單元的面積， ??a ?i， b ?i ，c ?i 以及 A 的表達式為

a ?i= x ?j y ?k -x ?k y ?j ，

b ?i= y ?j -y ?k ，

c ?i= x ?k -x ?j ，

A=△124+△324=△132+△134 ，

i，j，k=1，2，3 i，j，k=2，3，4 i，j，k=3，4，1 i，j，k=4，1，2 ?.

進一步，四邊形面積坐標還可以用等參坐標 （ξ，η） 表示. 如圖3所示，笛卡爾坐標（ x ， ?y ）和等參坐標 （ξ，η） 之間存在下述等參映射變換 ?F ?K ： ?K ︿ = ?-1，1 ??2→K ：

x ?y ???=F ??K（ξ，η）= 1 4 ∑ ?4 ?i=1 （1+ ξ ??iξ）（1+ η ??iη） ??x ?i ?y ?i ?，

其中 （ x ?i， y ?i） 表示在四邊形四個頂點的笛卡爾坐標，

ξ ?1 ?ξ ?2 ?η ?1 ?η ?2 ???ξ ?3 ?ξ ?4 ?η ?3 ?η ?4 ??= ??-1 1 -1 -1 ??1 -1 1 1 ????（4）

由式（3）可知，面積坐標 ?L ?i（x，y） 是單元的局部坐標，通過等參變換改寫為以 （ξ，η） 為變量的表達式 ?L ?i（ξ，η） .我們將 x（ξ，η） 和 y（ξ，η） 代入式（3）并整理就可以得到

L ?1= 1 4 （1-ξ） ?g ?2（1-η）+ 1- g ?1 （1+η） ?，

L ?2= 1 4 （1-η）［ （1-g ?2） 1-ξ +

（1- g ?1）（1+ξ）］ ，

L ?3= 1 4 ?1+ξ ??g ?1（1-η）+ 1- g ?2 ?1+η ???，

L ?4= 1 4 （1+η） ?g ?1（1-ξ）+ g ?2（1+ξ）

（5）

文獻［11］中介紹了四邊形單元，以及任意點 P 的面積坐標和笛卡爾坐標之間的一階導數的變換，面積坐標沿邊 ?L ?i=0 （ i =1，2，3，4）的任意冪函數的線積分公式，以及面積坐標的任意冪函數的面積積分公式.

注2 ??與等參坐標相比，面積坐標下的冪函數能給出精確積. 從計算的角度看，相比于等參坐標單元剛度矩陣的數值積分精度更好，結果更準確.

3 基于廣義協調元的雜交應力法

本節基于雜交應力有限元框架對面積坐標四邊形元進行一定的改進.

3.1 記號與基本概念

考慮如下的線彈性問題：

-div ?σ=f ， σ=Dε（u） ?in Ω，

σ·n | ???Γ ??1= T - ?， ?u | ???Γ ??0=0 ， on ??Γ = ?Ω = ??Γ ??1 - ∪ ??Γ ??0 - ???（6）

其中Ω ?∈R ?2 是有界開集; u 代表位移; σ 是應力張量; ε（u）= 1 2 （SymbolQC@

+ SymbolQC@

T）u 表示應變; D 代表彈性模量矩陣; f 是規定的體力，即分布載荷， ?T - ?為給定的表面牽引力; ???Γ ??1 - ?為Ω的邊界 ??Ω中表面牽引力的部分. 在基于雜交變分原理的廣義協調有限元推導中，位移場表示為 v= v ?c+ v ?I ， 其中 ?v ?c 表示單元的協調位移， ?v ?I 表示非協調的內部位移， τ 表示假定的分片獨立應力，位移和應力逼近的有限維空間 ??U ?h= U ?h ?c U ?h ?I 以及 ?V ?h 分別滿足

U ?h ?c∈ U ?c = ?v∈ ?∏ ?K∈ T ?h ?H ?1 K ???2；v | ???Γ ??0=0 ?，

U ?h ?I ??K ?span{bubble functions}，

V ??h∈V=∏ ?K∈ T ??h H（ div ;K），

其中 ?T ?h=∪{K} 表示Ω上四邊形剖分的全體; ??H ?s（K），s≥0 是通常的整數階Sobolev空間; ?H（ div ;K） 的定義為

H（ div ;K）= τ= ??τ ?11， τ ?22， τ ?12 ??T∈ L ?2 ?K ?3 ??.

值得注意的是，文獻［14，15］基于組合雜交變分原理引入了“能量協調性”的概念，對不相容內部位移 ??v ?I 進行了研究，通過能量相容條件

∮ ???Kτ·n· v ??I d s=0， ?τ， ?v ??I ??（7）

并基于修正的Hellinger-Reissner變分原理對面積坐標下的四邊形有限元進行了改進.

3.2 九參數雜交應力格式

當我們采用含內部位移的四節點位移構造雜交應力有限元時，線彈性問題式（6）對應的Hellinger-Reissner變分原理為 ?［12，14］：

Π ??σ ??h， u ??h = ?inf ??v∈ U ??h ??sup ??τ∈ V ??h ?Π （τ，v） ??（8）

其中能量泛函Π （τ，v） 具體表示如下：

Π （τ，v）=∑ ?K ［- 1 2 ∫ ??Kτ· D ?-1τ dΩ +

∫ ??Kτ·ε（v） dΩ -∮ ??Γ ?1∩ ?K T - · v ??c d s-

∮ ???Kτ·n· v ??I d s-∫ ??Kf·v dΩ ］.

注意到

∫ ??Kτ·ε（v） dΩ =∫ ??Kτ·ε（ v ??c） dΩ + ?∫ ??Kτ·ε（ v ??I） dΩ ， v∈ U ??h ?I，

我們借助分部積分公式可得

∫ ??Kτ·ε（ v ??I） dΩ -∮ ???Kτ·n· v ??I d s=

-∫ ??K div τ· v ??I dΩ .

泛函式（8）中的 ?Π ?可以改寫為

Π （τ，v）=∑ ?K ［- 1 2 ∫ ??Kτ· D ?-1τ dΩ +

∫ ??Kτ·ε（ v ??c） dΩ -∫ ??K div τ· v ??I dΩ -

∮ ??Γ ??1∩ K T - · v ??c d s-∫ ??Kf·v dΩ ］.

采用含內部位移的四節點廣義協調元AGQ6 ?［9］的位移場，以及九參數線性應力模式，我們提出一種面積坐標下的雜交應力四邊形有限元AGQ- LQ ?6. ??K∈ T ?h ，用 ?U ?h ?c 表示面積坐標下相應的節點的位移空間，

U ?h ?c={ v ?c； v ?c | ?K∈ span { L ?1， L ?2， L ?3， L ?4，

（ L ?3- L ?1），（ L ?4- L ?2）} ?2} ，

即， ??v ?c∈ U ?h ?c ，具有如下形式

v ?c | ?K= ??N ?1 0 ?0 ?N ?1 ??N ?2 0 ?0 ?N ?2 ??N ?3 0 ?0 ?N ?3 ??N ?4 0 ?0 ?N ?4 ??q ?c= N ?c q ?c ，

其中節點位移

q ?c= ??u ?1 v ?1 u ?2 v ?2 u ?3 v ?3 u ?4 v ?4 ??T ，

形函數 ?N ?i 為

N ?o ?i=- ?g ?k 2 + L ?i+ L ?j+ ξ ?i η ?i g ?kP ?（9）

i，j，k=1，2，3 i，j，k=2，3，4 i，j，k=3，4，1 i，j，k=4，1，2 ?.

其中

P= 3 ?L ?3- L ?1 ??L ?4- L ?2 - ?L ?3- L ?1 ??g ?2- g ?3 ?1+ g ?1 g ?3+ g ?2 g ?4

- ??g ?1- g ?2 ??L ?4- L ?2 - 1 2 ??g ?2 g ?4- g ?1 g ?3 ?1+ g ?1 g ?3+ g ?2 g ?4 .

用 ?U ?h ?I 表示面積坐標下相應的內部非協調位移的空間，

U ?h ?I={ v ?I； v ?I | ?K∈ span ???L ?1 L ?3， L ?4 L ?2 ??2，

K∈ T ?h} ，

即， ??v ?I∈ U ?h ?I 有如下形式

v ?I | ?K= ??N ?λ1 0 ?0 ?Nλ ?1 ??N ?λ2 0 ?0 ?N ?λ2 ??q ?I= N ?I q ?I .

其中內部位移 ?q ?I∈ ??R ?4 ，形函數 ?N ?λi 為

N ?λ1= L ?1 L ?3 ， ??N ?λ2= L ?2 L ?4 ?（10）

那么整個位移空間可以記為 ?U ?h ?W= U ?h ?c ?U ?h ?I . ??v∈ U ?h ?W 具有如下形式

v | ?K= （v ?I+ v ?c） | ?K=［ ?N ?cN ?I］ ?q ?c ?q ?I ?.

注3 ??當單元形狀退化為平行四邊形時，有

g ?1= g ?2= 1 2 ， L ?3- L ?1= ξ 2 ， L ?4- L ?2= η 2 ?.

則形函數式（9）和式（10）會變為

N ?o ?i= 1 4 ?1+ ξ ?iξ ?1+ η ?iη ， i=1，2，3，4，

N ?λ1= 1 16 （1- ξ ?2） ，

N ?λ2= 1 16 （1- η ?2） .

此時節點位移的形函數和內部位移的形函數與Q6元相同，在形式上只相差了常數倍.

當使用面積坐標時，我們用 ?V ?h ?1 表示相應的分片線性應力子空間，

V ?h ?1={τ；τ | ?K∈ span {1，（ L ?3- L ?1），

（ L ?4- L ?2）} ?3} ，

其中， τ∈ V ?h ?1 具有形式 τ | ?K= ?Φ ??β ??β ?（τ） ， ?β ?（τ）∈ ??R ?9 ，

Φ ??β= ?1 ????L ?3- L ?1 1 ???L ?4- L ?2 ?L ?3- L ?1 1 ???L ?4- L ?2 ?L ?3- L ?1 ????L ?4- L ?2 ??.

注4 ??面積坐標中只有兩個是相互獨立的坐標分量，所以在分片應力中， （ L ?3- L ?1）和（ L ?4- L ?2） 可以看作兩個獨立的分量.

AGQ- LQ ?6基于下述修正的Hellinger-Reissner 變分原理：

Π ???A-L ?σ ??h， u ??h = ?inf ??v∈ U ??h ?w ??sup ??τ∈ V ??h ?1 ??Π ???A-L（τ，v） ?（11）

其中

Π ???A-L（τ，v）=∑ ?K ［- 1 2 ∫ ??Kτ· D ?-1τ dΩ +

∫ ??Kτ·ε（v） dΩ -∮ ??Γ ??1∩ K T - · v ??c d s-

∫ ??Kf· v ??c dΩ ］.

因此，與式 （11） 相應的有限元格式為：

求 ??σ ?h， u ?h= u ?c+ u ?I ∈U ?h ?c ?U ?h ?I ?，使得

a ?σ ?h，τ -b τ， u ?h =0， τ∈ V ?h ?1 ，

b ?σ ?h， v ?I+ v ?c =f ?v ?c +g ?v ?c ，

v= v ?I+ v ?c∈ U ?h ?W ?（12）

其中

α（σ，τ）=∑∫ ??Kσ· D ?-1τ dΩ ，

b（τ，v）=∑∫ ??Kτ·ε（v） dΩ ，

f（v）=∑∫ ??Kf·v dΩ

g（v）=∑∮ ??Γ ??1∩ K T - · v ??c d s.

下面我們給出AGQ- LQ ?6單元剛度矩陣的推導，單元應變表示為

ε（v）=ε（ v ?c）+ε（ v ?I）= B ?c q ?c+ B ?c q ?I=

［ ?B ?cB ?I］ ??q ?c ?q ?I ??，

B ?c = ???B ?c1 ???B ?c2 ?????B ?c3 ???B ?c4 ??，

B ?ci = ????N ?0 ?i ?x ?0 0 ???N ?0 ?i ?y ????N ?0 ?i ?y ????N ?0 ?i ?x ??， i=1，2，3，4，

B ?I = ???B ?I1 ???B ?I2 ??，

B ?Ii = ????N ?λi ?x ?0 0 ???N ?λi ?y ????N ?λi ?y ????N ?λi ?x ???i=1，2 .

將式（12）中的積分項寫為

a（σ，τ） ??（K）=∫ ??Kσ· D ?-1τ dΩ =

β ??σ ???T∫ ?1 ?-1 ?∫ ?1 ?-1 ??Φ ??β ??T D ?-1 ?Φ ???β J ?d ξ d η β ??σ

β ??σ ???TH β ??σ，

b（τ，v） ??（K）=∫ ??Kτ·ε（v） dΩ =

β ??σ ???T∫ ?1 ?-1 ?∫ ?1 ?-1 ??Φ ???β ??T［ ?B ??cB ??I］ J ?d ξ d η ??q ??c ?q ??I

（ β ??σ） ??T［ ?G ??cG ??I］ β ??σ，

f（ v ??c） ??（K）+ g（ v ??c） ??（K）=∫ ??Kf· v ??c dΩ +

∮ ??Γ ??1∩ K T - · v ??c d s ［ Q ??c0］ ??q ??c ?q ??I ?.

從式（12）中我們得到單元 K 上應力與位移的關系為：

β ?σ= H ?-1［ ?G ?cG ?I］ ??q ?c ?q ?I ??.

再代入到變分形式（12）的第二個等式中可得

G ?T ?c H ?-1 G ?c ?G ?T ?c H ?-1 G ?I ?G ?T ?I H ?-1 G ?c ?G ?T ?I H ?-1 G ?I ????q ?c ?q ?I ?= ???Q ?c ?T 0 ???（13）

對內部位移參數采取與Wilson 元類似的靜態冷凝方法，在式（13）中消去

q ?I=- ??G ?T ?I H ?-1 G ?I ??-1 G ?T ?I H ?-1 G ?c q ?c ，

則可將式（12） 轉化為僅含有節點參數的常規單元剛度方程 ?K ?8×8 q ?c= ?Q ?c ?T， ?其中8×8的單元剛度矩陣具體表示為

K ?8×8= G ?T ?c H ?-1 G ?c- G ?T ?c H ?-1 G ?I ??G ?T ?I H ?-1 G ?I ??-1

G ?T ?I H ?-1 G ?c .

單元剛度矩陣的顯式表達式可以通過面積積分公式獲得. 因為數值積分法對計算機編碼更方便，所以我們也可以通過等參坐標求出單元剛度矩陣. 與 4 節點等參元的積分形式相比，在面積坐標下形函數中沒有雅可比的逆，通過面積坐標積分公式便可以計算得到 ［K］ 的精確值.

3.3 收斂性分析

我們首先介紹網格的正則細分. 對于每一個四邊形 K∈ T ?h ，若 ?h ?K， h′ ?K 和 ?θ ?K ?i 分別為單元 K 的直徑，最小邊長度和與頂點相關的角，則 ?T ?h 是一個正則細分需滿足以下條件 ?［16］：

存在常數 ?γ′ 和 γ 使得下列不等式對所有 K∈ T ?h 都一致成立：

h ?K≤ ?γ′h′ ?K ，

max ??1≤i≤4 ??cos ?θ ?K ?i ≤γ≤1 .

雜交應力元AGQ- LQ ?6的收斂需要滿足如下更嚴格的細分條件：

在網格參數 h= ?max ??K ?h ?K ?趨于0的情況下，單元 K 的兩條對角線中點連線的距離 ?d ?K 是 ?O（ h ?2） 的.

當AGQ- LQ ?6單元形狀退化為平行四邊形時，有下式成立：

g ?1= g ?2= 1 2 ， L ?3- L ?1= ξ 2 ， L ?4- L ?2= η 2 .

位移形函數式（9）和式（10）以及分片線性應力矩陣 Φ ??β 寫為

N ?o ?i= 1 4 ?1+ ξ ?iξ ?1+ η ?iη ， i=1，2，3，4，

N ?λ1= 1 16 （1- ξ ?2） ，

N ?λ2= 1 16 ?1- η ?2 ，

Φ ??β= ?1 ξ 2 ?η 2 ????1 ξ 2 ?η 2 ????1 ξ 2 ?η 2 ???.

此時，AGQ- LQ ?6元和文獻［17］中的 LQ ?6元的位移逼近空間以及應力逼近空間相同. 換句話說，根據 LQ ?6元的收斂性分析，我們就能推導出AGQ- LQ ?6元對于問題（6）有限元解的存在唯一性證明，并得到與 LQ ?6元一樣的誤差估計，即在 ?L ?2 范數意義下位移逼近能達到2階精度，應變逼近能達到1階精度，應力逼近能達到1階精度.

當AGQ- LQ ?6形狀不是矩形的時候，單元的位移逼近我們采用的是廣義協調元AGQ6的位移場. 針對廣義協調元的收斂，Shi ?［17］通過FEM檢驗給出過證明. 因為非平行四邊形情況下的雜交元AGQ- LQ ?6較為復雜，本文沒有給出嚴格的理論證明，但通過一系列工程算例來驗證了單元的收斂性，并展示的數值算例結果中.結果顯示，AGQ- LQ ?6元的位移逼近可以達到2階精度，應變逼近可以達到1階精度，應力逼近可以達到1階精度，符合預期的理論結果.

注5 ??在后面的數值算例中，我們均采用二分法的細分方法，即連接四邊形對邊中點，這種網格細分方法滿足了雜交元AGQ- LQ ?6的收斂條件.

4 數值算例

我們選取一些代表性數值考核算例，通過比較新方法與幾個相關的四節點四邊形元 ?［9，12，18］的數值結果來驗證了改進效果.

例4.1 ?（網格畸變測試） 在這個測試中，我們將懸臂梁剖分為兩個四邊形單元，如圖4所示. 單元扭曲程度由擾動參數 e 來決定， 尖端點 A 的最大位移數據見表1，其中的楊氏模量 E =1500，泊松比 ν =0.25.

表1中的數據顯示，四邊形面積元 AGQ6不受單元扭曲程度影響，改進的面積元AGQ- LQ ?6也能做到不受參數 e 的影響.

例4.2 ?（Cook平面應力板問題） 板的左端固定，右端施加均勻分布的單位載荷，楊氏模量 E ?= 1，泊松比 ν = 1/3，見圖5. 我們的目的是通過斜網格剖分來檢驗四邊形單元的性能. 鑒于此模型沒有解析解，我們提供的參考值為相對精細網格下的計算結果. 在 C 點的最大位移， A 點的最大主應力和 B 點的最小主應力的數值結果見表2.

表2中的數據顯示，在各種網格情況下，改進的面積元AGQ- LQ ?6的位移和應力結果比AGQ6元和 LQ ?6元更接近參考值.

例4.3 ?（非多項式解析解的平面應變測試） 泊松locking-free測試的算例基于懸臂梁的模型來說明. 我們考慮更為復雜的初值條件和解析解為非多項式的情況. 考慮區域Ω=［0，10］×［-1，1］，楊氏模量 E ?=1500，泊松比 初始值ν =0.3，如圖6. 令體力為

f=E π ?2 ??cos ?πx ?sin ?πy ?- sin （πx） cos （πy） ??.

則問題的精確解為

u= u=（1+ν） cos （πx） sin （πy） ?-2 （1-ν） ?2xy v=- 1+ν ?sin ?πx ?cos ?πy ??+ ?1-ν ??2 x ?2+ν 1+ν ??y ?2-1 ??，

σ=E ?-π sin ?πx ?sin ?πy -2y 0 0 ?sin ?πx ?sin ?πy ??.

我們采用二分法進行網格細分. 表3～表5分別給出了位移的誤差 ???u- u ?h ??0 ??u ??0 ，應變的誤差 ??ε（u）-ε ?u ?h ???0 ??ε（u） ??0 ?和應力的誤差 ???σ- σ ?h ??0 ??σ ??0 ?在不同網格下的數值結果. 值得注意的是， ECQ ?4， LQ ?6和AGQ6的位移和應力誤差階數結果從20×4的網格開始算.因在初始網格下的位移誤差偏大，這對階數的計算有所影響.

從表中數據，我們可以看到：在 ν →0.5時，改進的面積元AGQ- LQ ?6的位移逼近精度能達到2階，應變和應力逼近精度能達到1階，且相比于 LQ ?6和AGQ6的數值結果計算精度有一定的提高.特別地，粗網格下的數值結果比原單元數值性能明顯提升.

5 結 論

本文針對二維線彈性問題提出了一種基于面積坐標的雜交應力四邊形有限元AGQ- LQ ?6. 數值算例顯示，該方法既能保持面積坐標廣義協調元對網格畸變的不敏感性和雜交應力元精度，又能有效克服泊松閉鎖現象.

參考文獻：

［1］ ??Wilson E L， Taylor R L， Doherty W P， ?et al . Incompatible displacement modes ［M］. New York： Academic Press， 1973.

［2］ ?Wilson E L， Taylor R L， Beresford P J. A nonconforming element for stress analysis ［J］. Int J Numer Meth Eng， 1976， 10： 1211.

［3］ ?Hellinger E. Die allgemeinen ansatze der mechanik der kontinua， enz ［J］. Math Wis， 1914， 4： 602.

［4］ ?Reissner E. On a variational theorem in elasticity ［J］. J Math Phys， 1950， 29： 90.

［5］ ?陳豫眉， 周亞婧. 基于 Maxwell 模型的線性粘彈性波傳播問題的質量集中雜交應力有限元法［J］. 四川大學學報： 自然科學版， 2021， 58： 061001.

［6］ ?付衛， 王皓， 張世全. 特征值問題的組合雜交有限元方法［J］. 四川大學學報： 自然科學版， 2017， 54： 708.

［7］ ?劉軒宇， 羅鯤， 王皓. 拋物型積分微分方程的新型全離散弱 Galerkin 有限元法［J］. 四川大學學報： 自然科學版， 2020， 57： 830.

［8］ ?張田田， 張百駒. 非線性彈性問題的增強混合有限元方法［J］. 四川大學學報： 自然科學版， 2021， 58： 011004.

［9］ ?Chen X M， Cen S， Long Y Q， ?et al . Membrane elements insensitive to distortion using the uadrilateral area coordinate method ［J］. Comput Struct， 2004， 85： 35.

［10］ ?Long Y Q， Li J X， Long Z F， ?et al . Area coordinates used in quadrilateral elements ［J］. Comput Mech， 1999， 19： 533.

［11］ Long Z F， Long Y Q， Li J X， ?et al . Some basic formulae for area coordinates used in quadrilateral elements ［J］. Comput Mech， 1999， 19： 841.

［12］ Zhou T X， Xie X P. Optimization of stress modes by energy compatibility for 4-node hybrid quadrilaterals ［J］. Math Numer Sin， 2004， 59： 293.

［13］ Zhou T X. Stabilized hybrid finite element methods based on combination of saddle point principles of elasticity problem ［J］. Math Comp， 2003， 72： 1655.

［14］ Zhou T X. Finite element method based on combination of ‘saddle pointvariational formulations ［J］. Sci China E， 1997， 30： 285.

［15］ Xie X P. An accurate hybrid macro-element with linear displacements ［J］. Comm Numer Meth Eng， 2005， 21： 1.

［16］ Yu G Z， Xie X P， Carstensen C. Uniform convergence and a posteriori error estimation for assumed stress hybrid finite element methods ［J］. Comput Meth Appl Mech Eng， 2011， 200： 2421.

［17］ Shi Z C. The F-E-M-Test for convergence of nonconforming finite elements ［J］. Math Comput， 1987， 49： 391.

［18］ Sumihara K， Pian T H H. Rational approach for assumed stress finite elements ［J］. Numer Meth Eng， 1984， 20： 1685.