卿光輝，王博鰲

（中國民航大學(xué)航空工程學(xué)院，天津 300300）

協(xié)調(diào)位移有限元法采用了一般線性協(xié)調(diào)位移元的插值多項式，但其位移結(jié)果和應(yīng)力結(jié)果的收斂速度相對較差，并存在自鎖問題，而非協(xié)調(diào)位移項的引入[1-2]是改善該問題的有效方法。其次，位移法通常需利用本構(gòu)關(guān)系來計算應(yīng)力，由于應(yīng)變-位移是微商關(guān)系，也是應(yīng)力結(jié)果精度不如位移結(jié)果的原因之一。另一種單變量有限元法是基于最小余能原理的應(yīng)力法，由于應(yīng)力平衡元要求在張量空間中建立所需的平衡關(guān)系和正則條件都比較困難，使其并沒有得到充分的發(fā)展。為了解決該問題，卞學(xué)鐄[3]于1964年從最小余能原理出發(fā)，通過拉格朗日乘子法（也可直接用H-R變分原理）建立了雜交應(yīng)力有限元方法。雜交應(yīng)力元的思想是保持單元內(nèi)部應(yīng)力的平衡，但放松了正則條件，其結(jié)果比位移法更好，且在解決體積剪切自鎖問題和不可壓縮材料問題上有一定的優(yōu)勢。但雜交應(yīng)力元法為了實現(xiàn)位移法求解，要先對應(yīng)力場變分后得到歐拉方程，利用其消去應(yīng)力變量，再對位移變量變分，才能得到位移與外載荷的關(guān)系方程，因而要求對單元級的矩陣求逆。另外，雜交應(yīng)力元的應(yīng)力模式選擇比較煩瑣。

混合法是包含兩類或兩類以上變量的有限元方法。Herrmann[4]和Qukit等[5]最早用該方法分析板殼問題。很多學(xué)者針對混合法中不正定線性方程組的求解問題有一定研究[6]。用混合法分析板殼問題，最明顯的優(yōu)勢是應(yīng)力結(jié)果精度高，且收斂速度快，但一般混合法涉及多類變量時，對計算機的性能要求較高，計算效率不如單變量位移法。另一方面，在網(wǎng)格較粗的情況下，以一般插值多項式表達位移和應(yīng)力變量的混合法數(shù)值結(jié)果存在較明顯的振蕩現(xiàn)象[6]。

文獻[7]在“保辛”理論[8-9]的指導(dǎo)下，從修正的H-R變分原理出發(fā)，建立對偶變量塊體混合元，該方法實質(zhì)上也是位移法，但與雜交應(yīng)力元法不同，沒有應(yīng)力模式的選取要求。

根據(jù)文獻[7]的理論背景和構(gòu)建過程，目前的方法可簡稱為聯(lián)合有限元法。該方法以修正的H-R變分原理和勢能原理為基礎(chǔ)，建立線形求解方程，避免通過本構(gòu)關(guān)系求解應(yīng)力變量，從而提高應(yīng)力計算精度、收斂度及計算效率。同時，不用放松正則條件及應(yīng)力模式的選取，理論上更加嚴謹。通過數(shù)值算例，與ABAQUA的計算結(jié)果進行對比分析，該方法具有一定優(yōu)越性。

1 理論分析

1.1 修正的H-R變分原理與部分混合元列式

其中：σo=[σxzσyzσzz]T為平面外應(yīng)力；Φ11、Φ21和 Φ22為材料參數(shù)矩陣；D1，D2和Dz為微分算子矩陣，即

平面內(nèi)應(yīng)力可由平面外應(yīng)力變量和位移變量表示為

為了實現(xiàn)結(jié)點上位移與平面外應(yīng)力均“保辛”[10-11]的目標，應(yīng)堅持求解各類變量的系數(shù)矩陣對稱原則。設(shè)平面外應(yīng)力σo和位移u在空間域內(nèi)離散后的表達式為

其中：pe=[σexzσeyzσezz]T；qe=[uevewe]T；N=diag（NeNeNe），Ne為單元各結(jié)點組成的向量，Nei=（1+ξiξ）（1+ ηiη）（1+ ζiζ）/8。

將式（4）和式（5）代入式（1）可得

其中

將pe和qe視為相互獨立的變量，則由δΠMR（pe，qe）=0得到歐拉方程如下

文獻[7]用以下公式消除歐拉方程（8）的應(yīng)力與位移的耦合關(guān)系，即

將式（9）代入式（8）中，可得

實際上，式（9）包含了對單元級的K11n求逆。在網(wǎng)格數(shù)增加或單元類型的單元結(jié)點數(shù)較多時，如20結(jié)點實體元，計算效率會明顯下降。同時，也假設(shè)勢能原理的位移條件事先滿足u-=0，即

其中：D為通常位移法中的算子操作矩陣；C為材料本構(gòu)關(guān)系的系數(shù)矩陣。

假設(shè)式（11）和式（1）是對同一工程結(jié)構(gòu)靜力學(xué)問題的表述。將式（5）代入式（11）中，有

則將式（13）替代式（8），得到了關(guān)于位移變量和平面外應(yīng)力變量的單元列式如下

對式（14）或式（15）兩邊分別求和，則問題的控制方程分別為

其中：K11= ∑K11n；K12= ∑K12n；K22= ∑K22n；p= ∑pe；q=∑qe；f= ∑fe。

與位移法同樣，可先通過式（17）獨立求位移的解，然后將結(jié)果代入式（16）求解平面外應(yīng)力。

從以上理論的推導(dǎo)和最終的單元列式（14）和式（15）的形式來看：聯(lián)合有限元法避免了通過較為繁瑣的式（10）求解位移變量，采用較為簡便的位移法求解，從而提高了計算效率；在求解應(yīng)力時，沒有采用計算精度較低的本構(gòu)關(guān)系進行求解，而通過修正的H-R變分原理中應(yīng)力計算公式計算出應(yīng)力結(jié)果；將兩個方程通過線形方程組形式進行同時求解，又進一步提高了計算效率。此方法可稱由H-R變分原理和勢能原理的歐拉方程聯(lián)合得到的單元列式（14）和式（15）為聯(lián)合單元列式（CEF）。

1.2 應(yīng)力求解方法

1.2.1 平面外應(yīng)力的求解方法

對于以聯(lián)合單元列式（14）為基礎(chǔ)的控制方程（16）來說，系數(shù)矩陣K11是對稱“保辛”的。同時，與求解位移相同，便于施加平面外應(yīng)力邊界條件，確保了邊界結(jié)點上的結(jié)果與事先給定的應(yīng)力數(shù)值一致。

1.2.2 平面內(nèi)應(yīng)力的線性方程組方法

為了便于平面內(nèi)應(yīng)力邊界條件的引入，可進一步考慮具體結(jié)構(gòu)中同一材料的線性方程組形式。將式（4）和式（5）代入式（3）中，然后兩邊求和，可分別構(gòu)建關(guān)于同一材料平面內(nèi)應(yīng)力的線性方程組為

其中：K=∑In，In是24×24的單位矩陣（假設(shè)是8結(jié)點塊體元）；pi為待求的平面內(nèi)應(yīng)力向量；Γ=∑（Φ21pe+Φ22（DzN）qe），（DzN）矩陣需根據(jù)具體等參坐標（ξ，η，ζ）的位置代入-1，0，1 等值。構(gòu)造式（18）時無需數(shù)值積分，速度快。

需進一步說明以下兩點：一方面，式（18）中K只有主對角線上存在元素，理所當然對稱“保辛”。另一方面，構(gòu)建式（18）的前提條件是同一材料的單元必須連續(xù)（如呈層狀的復(fù)合材料結(jié)構(gòu)中的同一材料層）。

1.3 非協(xié)調(diào)位移聯(lián)合元

根據(jù)文獻[1-2，15-16]中非協(xié)調(diào)位移元法，可將位移插值表達式分為協(xié)調(diào)部分和非協(xié)調(diào)部分，即

將式（19）代入式（11）后有

對式（20）求和可得

將K′22替代式（17）的K22可得非協(xié)調(diào)位移元的方程。聯(lián)合有限元方法的非協(xié)調(diào)元計算代碼采用了文獻[16]的理論公式。由于式（14）中的系數(shù)矩陣與式（15）的剛度矩陣相比較柔[4-5]，故沒有考慮式（14）的非協(xié)調(diào)形式。

2 數(shù)值算例

已知矩形板的長和寬a=b=1.0，厚度H=0.1。材料參數(shù)[17]為 E11=10E22=10E33，G12=G13=0.6E33，G23=0.5E33，μ12= μ13= μ23=0.25。板上表面（z=0.1）受垂直向下的均布載荷q=1作用，四邊簡支（SSSS）。

計算程序中，每一坐標方向采用兩點高斯積分。網(wǎng)格為12×12×4的四分之一模型，在同一臺集成4核（3.2 GHz）的機器上運算，聯(lián)合有限元法用時3.6 s，而文獻[7]中的方法用時4.2 s。

在以下的數(shù)據(jù)對比曲線圖中，8結(jié)點協(xié)調(diào)聯(lián)合元簡稱CCFEM8，非協(xié)調(diào)聯(lián)合元簡稱NCCFEM8，ABAQUA的8結(jié)點非協(xié)調(diào)位移元，簡稱NC3D8，各計算結(jié)果均進行了無量綱化處理。

圖1～圖6中橫坐標序號與網(wǎng)格劃分密度對應(yīng)關(guān)系如表1所示，其中，網(wǎng)格劃分密度為長×寬方向的單元數(shù)，厚度方向網(wǎng)格劃單元數(shù)均為4。圖7和圖8給出沿著厚度方向應(yīng)力σxy′和σyy′分布及對比結(jié)果，其中橫坐標為沿著厚度方向的坐標。

圖1 u′與精確解比較Fig.1 Comparison between u′and exact solution

圖2 w′與精確解比較Fig.2 Comparison between w′and exact solution

圖3 σxz′與精確解比較Fig.3 Comparison between σxz′and exact solution

圖4 σzz′與精確解比較Fig.4 Comparison between σzz′and exact solution

圖5 σxx′與精確解比較Fig.5 Comparison between σxx′and exact solution

圖 6 σxy′與精確解比較Fig.6 Comparison between σxy′and exact solution

表1 網(wǎng)格劃分的對應(yīng)關(guān)系Tab.1 Corresponding relation between horizontal axis and mesh models

圖7 σxy′沿厚度方向的分布Fig.7 Distribution of stress σxy′along thickness direction

圖8 σyy′沿厚方向的分布Fig.8 Distribution of stress σyy′along thickness direction

以變量靠近最大值與最小值位置為原則，關(guān)于位移和應(yīng)力的比較位置點（參見圖1～圖6）為

圖1～圖6中的誤差值是在網(wǎng)格12×12×4模型下各單元結(jié)果與精確解的比較。從以上數(shù)據(jù)曲線圖可得如下結(jié)果。

1）協(xié)調(diào)聯(lián)合元 CCFEM8 的位移 u′、w′和平面內(nèi)應(yīng)力 σxx′、σxy′的結(jié)果不如 ABAQUS 中非協(xié)調(diào)位移元NC3D8。實質(zhì)上，因為線性的協(xié)調(diào)位移元剛度矩陣的剛度大于真實材料的剛度，而引入了內(nèi)部結(jié)點非協(xié)調(diào)位移項的非協(xié)調(diào)位移元的插值多項式完備，因而其剛度矩陣的剛度更真實。所以，協(xié)調(diào)位移元的應(yīng)力精度是不太可能超越非協(xié)調(diào)位移元的應(yīng)力精度的。但對于協(xié)調(diào)元，聯(lián)合有限元方法的其他應(yīng)力結(jié)果均優(yōu)于NC3D8。CCFEM8元的應(yīng)力結(jié)果收斂穩(wěn)定性相當明顯，且各變量的收斂速度比較均衡。

2）NCCFEM8的位移和應(yīng)力結(jié)果均優(yōu)于NC3D8的結(jié)果，且收斂速度很快，尤其是其平面外應(yīng)力σxz′與NC3D8的結(jié)果比較，收斂優(yōu)勢明顯。從上述圖1～圖6可看出，NCCFEM8的位移和應(yīng)力結(jié)果的收斂速度和精確度基本一致；同時從圖7～圖8可看出，沿厚度方向各結(jié)點，應(yīng)力結(jié)果的精確度也同精確解基本一致。

3）對橫觀材料或正交各異性材料而言，一般有限元法的平面外應(yīng)力σxz′結(jié)果相對較差。實踐表明，若網(wǎng)格進一步細分，NCCFEM8的位移和應(yīng)力結(jié)果均可逼近精確解，而NC3D8的平面外應(yīng)力σxz′要逼近精確解卻可能需要更細的網(wǎng)格劃分（圖3）。

4）σzz′處于中性面的中間位置，3 種單元的 σzz′收斂速度快，其精度沒有明顯的區(qū)別（圖4）。

3 結(jié)語

聯(lián)合有限元方法是通過位移與外載荷關(guān)系的歐拉方程得到較精確的位移結(jié)果后，再用H-R變分原理關(guān)于應(yīng)力變量變分后的歐拉方程替代本構(gòu)關(guān)系求解應(yīng)力。因而，避免直接用本構(gòu)關(guān)系方程求解應(yīng)力，提高了應(yīng)力數(shù)值的精確度。

位移變量求解采用剛度矩陣對稱的線性方程組求解，應(yīng)力變量求解采用系數(shù)矩陣對稱的線性方程組求解。對于求應(yīng)力的線性方程可方便地引入應(yīng)力邊界條件。

從NCCFEM8的數(shù)值結(jié)果來看，其主要意義在于位移和應(yīng)力的快速收斂和結(jié)果的精確性，且位移和應(yīng)力的收斂速度和精確度基本一致。以有關(guān)板殼理論的H-R變分原理和勢能原理為基礎(chǔ)，同樣可簡便地演化出與聯(lián)合有限元方法相對應(yīng)的有限元方法。