高何金雨 張世全

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

關(guān)鍵詞 ：線彈性問題； 四邊形面積坐標(biāo)方法； 雜交應(yīng)力有限元； 泊松閉鎖現(xiàn)象

中圖分類號 ： O241.82 文獻(xiàn)標(biāo)識碼 ：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 引 言

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

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

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

本文基于雜交應(yīng)力有限元框架，采用Zhou和Xie在文獻(xiàn)［12，13］中提出的應(yīng)力模式對四邊形面積坐標(biāo)下的廣義協(xié)調(diào)元進(jìn)行改進(jìn)，進(jìn)而提出了一種基于面積坐標(biāo)的雜交應(yīng)力元AGQ- LQ ?6. 在該方法中，單元的位移逼近我們采用含內(nèi)部位移的四節(jié)點(diǎn)廣義協(xié)調(diào)元AGQ6，應(yīng)力逼近則采用九參數(shù)線性應(yīng)力模式. 最后，我們以標(biāo)準(zhǔn)的數(shù)值算例驗(yàn)證了該方法的數(shù)值性能.

后文的安排如下. 在第2節(jié)中我們給出四邊形面積坐標(biāo)的定義及面積元的構(gòu)造. 在第3節(jié)中，我們基于雜交應(yīng)力有限元框架給出改進(jìn)的面積元的離散格式. 我們在第4節(jié)中給出數(shù)值算例. 第5節(jié)總結(jié)了本文所得結(jié)果.

2 四邊形面積坐標(biāo)

2.1 四邊形的特征參數(shù)

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

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

（1）

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

注1 ??我們稱 ?g ?1 和 ?g ?2 為形狀特征參數(shù). 當(dāng) ?g ?1， g ?2 取一些特殊值時，圖1中的單元會轉(zhuǎn)變?yōu)橐恍┯刑囟ㄐ螤畹膯卧?/p>

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

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

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

2.2 四邊形面積坐標(biāo)的定義

如圖2所示，令 P ?為四邊形單元內(nèi)的任意一點(diǎn)，則 P 的面積坐標(biāo) （ 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）是由點(diǎn) P 分別和單元其中兩個相鄰頂點(diǎn)構(gòu)成的三角形的面積，單元四個頂點(diǎn)的坐標(biāo)表示為：

節(jié)點(diǎn)1： ?（ g ?2，1- g ?2，0，0）;

節(jié)點(diǎn)2： ?（ 0，1-g ?1， g ?1，0）;

節(jié)點(diǎn)3： ?（0，0，1- g ?2， g ?2） ;

節(jié)點(diǎn)4： ?（1- g ?1，0，0， g ?1） .

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

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 的表達(dá)式為

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 ?.

進(jìn)一步，四邊形面積坐標(biāo)還可以用等參坐標(biāo) （ξ，η） 表示. 如圖3所示，笛卡爾坐標(biāo)（ x ， ?y ）和等參坐標(biāo) （ξ，η） 之間存在下述等參映射變換 ?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） 表示在四邊形四個頂點(diǎn)的笛卡爾坐標(biāo)，

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

由式（3）可知，面積坐標(biāo) ?L ?i（x，y） 是單元的局部坐標(biāo)，通過等參變換改寫為以 （ξ，η） 為變量的表達(dá)式 ?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）

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

注2 ??與等參坐標(biāo)相比，面積坐標(biāo)下的冪函數(shù)能給出精確積. 從計(jì)算的角度看，相比于等參坐標(biāo)單元剛度矩陣的數(shù)值積分精度更好，結(jié)果更準(zhǔn)確.

3 基于廣義協(xié)調(diào)元的雜交應(yīng)力法

本節(jié)基于雜交應(yīng)力有限元框架對面積坐標(biāo)四邊形元進(jìn)行一定的改進(jìn).

3.1 記號與基本概念

考慮如下的線彈性問題：

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

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

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

+ SymbolQC@

T）u 表示應(yīng)變; D 代表彈性模量矩陣; f 是規(guī)定的體力，即分布載荷， ?T - ?為給定的表面牽引力; ???Γ ??1 - ?為Ω的邊界 ??Ω中表面牽引力的部分. 在基于雜交變分原理的廣義協(xié)調(diào)有限元推導(dǎo)中，位移場表示為 v= v ?c+ v ?I ， 其中 ?v ?c 表示單元的協(xié)調(diào)位移， ?v ?I 表示非協(xié)調(diào)的內(nèi)部位移， τ 表示假定的分片獨(dú)立應(yīng)力，位移和應(yīng)力逼近的有限維空間 ??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 是通常的整數(shù)階Sobolev空間; ?H（ div ;K） 的定義為

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

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

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

并基于修正的Hellinger-Reissner變分原理對面積坐標(biāo)下的四邊形有限元進(jìn)行了改進(jìn).

3.2 九參數(shù)雜交應(yīng)力格式

當(dāng)我們采用含內(nèi)部位移的四節(jié)點(diǎn)位移構(gòu)造雜交應(yīng)力有限元時，線彈性問題式（6）對應(yīng)的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Ω ］.

采用含內(nèi)部位移的四節(jié)點(diǎn)廣義協(xié)調(diào)元AGQ6 ?［9］的位移場，以及九參數(shù)線性應(yīng)力模式，我們提出一種面積坐標(biāo)下的雜交應(yīng)力四邊形有限元AGQ- LQ ?6. ??K∈ T ?h ，用 ?U ?h ?c 表示面積坐標(biāo)下相應(yīng)的節(jié)點(diǎn)的位移空間，

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 ，

其中節(jié)點(diǎn)位移

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

形函數(shù) ?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 表示面積坐標(biāo)下相應(yīng)的內(nèi)部非協(xié)調(diào)位移的空間，

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 .

其中內(nèi)部位移 ?q ?I∈ ??R ?4 ，形函數(shù) ?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 ??當(dāng)單元形狀退化為平行四邊形時，有

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

則形函數(shù)式（9）和式（10）會變?yōu)?/p>

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） .

此時節(jié)點(diǎn)位移的形函數(shù)和內(nèi)部位移的形函數(shù)與Q6元相同，在形式上只相差了常數(shù)倍.

當(dāng)使用面積坐標(biāo)時，我們用 ?V ?h ?1 表示相應(yīng)的分片線性應(yīng)力子空間，

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 ??面積坐標(biāo)中只有兩個是相互獨(dú)立的坐標(biāo)分量，所以在分片應(yīng)力中， （ L ?3- L ?1）和（ L ?4- L ?2） 可以看作兩個獨(dú)立的分量.

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） 相應(yīng)的有限元格式為：

求 ??σ ?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單元剛度矩陣的推導(dǎo)，單元應(yīng)變表示為

ε（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）中的積分項(xiàng)寫為

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 上應(yīng)力與位移的關(guān)系為：

β ?σ= 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）

對內(nèi)部位移參數(shù)采取與Wilson 元類似的靜態(tài)冷凝方法，在式（13）中消去

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

則可將式（12） 轉(zhuǎn)化為僅含有節(jié)點(diǎn)參數(shù)的常規(guī)單元剛度方程 ?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 .

單元剛度矩陣的顯式表達(dá)式可以通過面積積分公式獲得. 因?yàn)閿?shù)值積分法對計(jì)算機(jī)編碼更方便，所以我們也可以通過等參坐標(biāo)求出單元剛度矩陣. 與 4 節(jié)點(diǎn)等參元的積分形式相比，在面積坐標(biāo)下形函數(shù)中沒有雅可比的逆，通過面積坐標(biāo)積分公式便可以計(jì)算得到 ［K］ 的精確值.

3.3 收斂性分析

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

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

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

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

雜交應(yīng)力元AGQ- LQ ?6的收斂需要滿足如下更嚴(yán)格的細(xì)分條件：

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

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

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

位移形函數(shù)式（9）和式（10）以及分片線性應(yīng)力矩陣 Φ ??β 寫為

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元和文獻(xiàn)［17］中的 LQ ?6元的位移逼近空間以及應(yīng)力逼近空間相同. 換句話說，根據(jù) LQ ?6元的收斂性分析，我們就能推導(dǎo)出AGQ- LQ ?6元對于問題（6）有限元解的存在唯一性證明，并得到與 LQ ?6元一樣的誤差估計(jì)，即在 ?L ?2 范數(shù)意義下位移逼近能達(dá)到2階精度，應(yīng)變逼近能達(dá)到1階精度，應(yīng)力逼近能達(dá)到1階精度.

當(dāng)AGQ- LQ ?6形狀不是矩形的時候，單元的位移逼近我們采用的是廣義協(xié)調(diào)元AGQ6的位移場. 針對廣義協(xié)調(diào)元的收斂，Shi ?［17］通過FEM檢驗(yàn)給出過證明. 因?yàn)榉瞧叫兴倪呅吻闆r下的雜交元AGQ- LQ ?6較為復(fù)雜，本文沒有給出嚴(yán)格的理論證明，但通過一系列工程算例來驗(yàn)證了單元的收斂性，并展示的數(shù)值算例結(jié)果中.結(jié)果顯示，AGQ- LQ ?6元的位移逼近可以達(dá)到2階精度，應(yīng)變逼近可以達(dá)到1階精度，應(yīng)力逼近可以達(dá)到1階精度，符合預(yù)期的理論結(jié)果.

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

4 數(shù)值算例

我們選取一些代表性數(shù)值考核算例，通過比較新方法與幾個相關(guān)的四節(jié)點(diǎn)四邊形元 ?［9，12，18］的數(shù)值結(jié)果來驗(yàn)證了改進(jìn)效果.

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

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

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

表2中的數(shù)據(jù)顯示，在各種網(wǎng)格情況下，改進(jìn)的面積元AGQ- LQ ?6的位移和應(yīng)力結(jié)果比AGQ6元和 LQ ?6元更接近參考值.

例4.3 ?（非多項(xiàng)式解析解的平面應(yīng)變測試） 泊松locking-free測試的算例基于懸臂梁的模型來說明. 我們考慮更為復(fù)雜的初值條件和解析解為非多項(xiàng)式的情況. 考慮區(qū)域Ω=［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 ??.

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

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

5 結(jié) 論

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

參考文獻(xiàn)：

［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 模型的線性粘彈性波傳播問題的質(zhì)量集中雜交應(yīng)力有限元法［J］. 四川大學(xué)學(xué)報(bào)： 自然科學(xué)版， 2021， 58： 061001.

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

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

［8］ ?張?zhí)锾铮?張百駒. 非線性彈性問題的增強(qiáng)混合有限元方法［J］. 四川大學(xué)學(xué)報(bào)： 自然科學(xué)版， 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.