肖映雄，謝凌潔

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

在二維問題有限元分析中，采用四邊形網(wǎng)格可以更好地反映變形體中的位移狀態(tài)和應(yīng)力狀態(tài).為了提高有限元計算精度，往往需要采用相應(yīng)的高階單元.這類單元因具有某些特殊的優(yōu)點，例如，能解決彈性問題的閉鎖現(xiàn)象(Poisson’s ratio locking)，使其在實際計算中被廣泛使用.但與線性元相比，高階單元節(jié)點數(shù)多，單元分析復(fù)雜，相應(yīng)離散化矩陣又具有病態(tài)性，利用通常的方法[1]求解有限元方程時其效率將大大降低.實際計算時，為了減少計算量，提高編程和計算效率，需要考慮采用階譜單元技術(shù)[2-4]，它在生成高階單元特性矩陣時能夠承襲低階單元的特性矩陣.最近，文獻(xiàn)[5]通過適當(dāng)增加棱內(nèi)“虛節(jié)點”和面內(nèi)“虛節(jié)點”的方式設(shè)計了一類四邊形升階譜單元，并結(jié)合矩陣的一維稀疏存儲技術(shù)(CSR 格式)實現(xiàn)了相應(yīng)的分層四邊形高階亞參元方法.

為了提高這種分層有限元分析的整體效率，需要為其提供高效求解器(Solvers).代數(shù)多重網(wǎng)格(Algebraic MultiGrid，AMG)法是求解大規(guī)模有限元離散化線性系統(tǒng)最為有效的代數(shù)方法之一[6-9].在解決實際問題時，如果還能利用易于獲取的某些幾何或分析信息，如節(jié)點坐標(biāo)、單元階次及微分方程的類型等，則可設(shè)計出具有更好計算效率的AMG(簡記GAMG)法[10-12].文獻(xiàn)[13]針對橢圓型定解問題，建立和分析了二次四邊形有限元方程的GAMG 法，但該方法需借助于節(jié)點基函數(shù)的緊支集性和能量極小來建立一系列代數(shù)判據(jù)，導(dǎo)致求解效率降低.對彈性力學(xué)問題，相應(yīng)的GAMG 法研究也取得了很多成果[14-15].最近，文獻(xiàn)[16]針對混凝土隨機(jī)骨料模型的三角形分層二次元離散化系統(tǒng)，設(shè)計了2 種簡單有效的預(yù)條件子，相應(yīng)預(yù)條件共軛梯度(Preconditioned Conjugate Gradient,PCG)法具有很好的計算效率.

本文針對分層四邊形高階元離散化系統(tǒng)，通過利用其系數(shù)矩陣對角分塊矩陣的代數(shù)性質(zhì)和分層結(jié)構(gòu)特性，設(shè)計了一種簡單、有效的預(yù)條件子，獲得了一種內(nèi)迭代效率得到顯著改善的PCG 法，從而提高了相應(yīng)有限元分析的整體計算效率.所采用方法的基本思想是將分層高階元離散化系統(tǒng)本質(zhì)性地化歸為Q4元離散化系統(tǒng)的快速求解，而對Q4元系統(tǒng)，目前已有很多高效的求解器，如GAMG 法.最后，通過對幾類實際問題分層Q8元和分層Q12元離散化系統(tǒng)進(jìn)行求解，驗證了方法的有效性.

1 分層高階四邊形亞參元法

在實際四邊形有限元分析中，為了提高有限元解精度，需要增加單元插值函數(shù)的次數(shù).本節(jié)簡單介紹文獻(xiàn)[5]中通過適當(dāng)增加棱內(nèi)“虛節(jié)點”和面內(nèi)“虛節(jié)點”方式設(shè)計的四邊形升階譜單元及其階譜函數(shù).

考慮如圖1 左圖所示的四邊形單元.為敘述方便，不妨設(shè)Q4單元的四個頂點分別為當(dāng)由Q4元升階為Q8元時，需要分別在4 條棱內(nèi)增加一個虛節(jié)點，記為依次下去，當(dāng)Q8元升階為Q12元時，需要分別在4 條棱內(nèi)再增加一個虛節(jié)點，當(dāng)由Q12元升階為Q17元時，需要分別在4 條棱內(nèi)再增加一個虛節(jié)點，在面內(nèi)增加一個虛節(jié)點當(dāng)由Q17元升階為Q22元時，需要分別在4 條棱內(nèi)再增加一個虛節(jié)點，需要在面內(nèi)再增加一個虛節(jié)點這樣，該單元當(dāng)由Q4元升階為Q22元時，在四條棱內(nèi)增加的虛節(jié)點分別為面內(nèi)增加的虛節(jié)點分別為如圖1 中的左圖所示，其對應(yīng)的基本單元如右圖所示.更高一階單元新增棱內(nèi)“虛節(jié)點”和面內(nèi)“虛節(jié)點”均可類似得到.

圖1 任意四邊形單元及對應(yīng)的基本單元Fig.1 An arbitrary quadrilateral element and the corresponding basic element

設(shè)最終升階后該單元的節(jié)點總數(shù)為Nn個，采用階譜單元(Hierachical element)進(jìn)行分析.二維四邊形階譜單元的形函數(shù)包括頂點形函數(shù)、虛節(jié)點(棱內(nèi)點和面內(nèi)點)形函數(shù)，利用其性質(zhì)可在基本單元上構(gòu)造階譜函數(shù)，列出結(jié)果如下：

(1)頂點(p=1)，這里，p為相應(yīng)插值函數(shù)中完全多項式的最高次數(shù)，對應(yīng)的階譜函數(shù)即為形函數(shù)：

其中，(ξi,ηi)為i點在基本單元中的坐標(biāo).

(2)棱內(nèi)虛節(jié)點(p≥2)，對應(yīng)的階譜函數(shù)分別為：

(3)面內(nèi)虛節(jié)點(p≥4)，對應(yīng)的階譜函數(shù)分別為：

上述各式中，?t為一維階譜單元中可行的階譜基函數(shù)，其形式不唯一.一種常用取法為

為了程序設(shè)計及敘述方便，我們將按以下順序?qū)Ω鞴?jié)點進(jìn)行編號：

相應(yīng)的階譜函數(shù)簡記為Hi，i=1,2,···,Nn.采用亞參元法進(jìn)行分析，即單元位移函數(shù)采用的變換為

而坐標(biāo)的插值采用變換：

設(shè)Q4元的單元剛度矩陣為

其中，Kij為2×2矩陣，i,j=1,2,3,4.

當(dāng)Q4元升階為Q8元時，相應(yīng)的單元剛度矩陣為

在生成高階單元剛度矩陣時始終能夠承襲低階單元的剛度矩陣，這樣可減少計算量，大大提高編程和計算效率.按照通常等參元法將各單元剛度矩陣進(jìn)行組集，可得到如下具有分層結(jié)構(gòu)的整體剛度矩陣：

完全類似地，可得到具有分層結(jié)構(gòu)的整體載荷向量FH,(p).最后，引入位移邊界條件，修改整體剛度矩陣和整體載荷向量，可得相應(yīng)的整體剛度方程分層有限元離散系統(tǒng)

其中，dH,(p)為分層p次四邊形單元位移向量.上述方程組(7)又可寫為

采用高階“階譜單元”亞參元分析的優(yōu)點之一，是它可以為相應(yīng)的分層有限元離散系統(tǒng)（7）或（8）設(shè)計快速求解方法帶來更大方便，以提高相應(yīng)有限元分析的整體計算效率.

2 分層有限元離散系統(tǒng)的PCG 法

2.1 預(yù)條件子的構(gòu)造分層有限元方程(7)的系數(shù)矩陣KH,(p)是對稱正定的，在求解該線性代數(shù)系統(tǒng)時常常采用預(yù)條件共軛梯度(PCG)法.設(shè)預(yù)條件子為BH,(p)，在PCG 算法中，最為關(guān)鍵的環(huán)節(jié)是給出z=BH,(p)r的快速計算方法，其中，r為給定的已知向量.一般而言，可從兩個方面判別一個預(yù)條件子是不是一個好的預(yù)條件子：首先，BH,(p)r的計算復(fù)雜性是否遠(yuǎn)低于 (KH,(p))?1r；其次，Cond(BH,(p)KH,(p))<

下面，基于文獻(xiàn)[16]中算法構(gòu)造思想，將基于標(biāo)量橢圓型問題AMG 的塊對角預(yù)條件子應(yīng)用于分層有限元離散系統(tǒng)(7)的求解.利用系數(shù)矩陣KH,(p)的分層結(jié)構(gòu)特性，不妨取

即可獲得相應(yīng)整體插值算子，記為Pvv，再利用Galerkin方法構(gòu)造下一個粗網(wǎng)格系數(shù)矩陣，即

然后，考慮(11)式中的第二個方程：Kmmzm=rm.通過計算條件數(shù)(參看后面的計算結(jié)果)可知矩陣Kmm是一個良態(tài)矩陣，因此，可直接調(diào)用共軛梯度(CG)法近似求解Kmmzm=rm，并設(shè)其解為，對應(yīng)的預(yù)條件子簡記為

相應(yīng)的z=BH,(p)r的計算方法總結(jié)如下：

算法(計算z=BH,(p)r)

步驟1由Q4元矩陣Kvv得到對角子塊，分別基于，并結(jié)合經(jīng)典RS 粗化算法，獲取GAMG 法所需的各網(wǎng)格層要素：插值算子與限制算子，粗網(wǎng)格矩陣，再通過選取合適的光滑算子，如Gauss-Seidel 迭代，建立GAMG 迭代方法；

步驟2給定初始迭代向量，調(diào)用一次GAMG 法求解Kvvzv=rv，近似解設(shè)為；

步驟3給定初始迭代值向量，調(diào)用s次CG 法求解Kmmzm=rm，近似解設(shè)為；

步驟4取，可得z=BH,(p)r≈.

注1本文所構(gòu)造的預(yù)條件子，充分利用了分層高階四邊形亞參元分析中總剛度矩陣的分塊結(jié)構(gòu)特征，將預(yù)條件子的計算非常自然地轉(zhuǎn)化為2 個解耦后的子系統(tǒng)的近似求解，這樣可分別利用這些子系統(tǒng)矩陣的代數(shù)性質(zhì)(例如，由所有“虛節(jié)點”構(gòu)成的子矩陣的條件數(shù)不是很大，為標(biāo)準(zhǔn)的良態(tài)矩陣)來選取高效求解方法，以提高PCG 內(nèi)迭代的計算效率.

2.2 算例及結(jié)果分析將基于上述預(yù)條件子的PCG(簡記為B-PCG)法應(yīng)用于幾類典型問題分層高階元離散系統(tǒng)(7)的求解，以驗證算法的有效性和魯棒性.在PCG 法中，控制精度取為ε=10?6.本文所有數(shù)值實驗均在Inter(R) Xeon(R)CPU E3-1230V6 處理器的Windows 10環(huán)境下進(jìn)行.

算例1(懸臂梁問題) 考慮如圖2 所示的端部受剪力作用懸臂梁的彎曲問題.設(shè)其長為L=10 m，高為D=1 m，取單位厚度，彈性模量E=2.1×106Pa，泊松比ν=0.167，左端(x=0)剪 力P=1 000 N/m，不考慮自重.對區(qū)域分別進(jìn)行結(jié)構(gòu)和非結(jié)構(gòu)四邊形網(wǎng)格剖分，圖3 分別給出了其中2 種不同類型的網(wǎng)格剖分樣圖.

圖2 端部受剪力作用的懸臂梁Fig.2 Cantilever beam with shear force at the end

圖3 2 種不同形式及規(guī)模的網(wǎng)格剖分樣圖Fig.3 Two samples of mesh with different types and sizes

首先，我們在上述2 類網(wǎng)格下計算了相應(yīng)的分層剛度矩陣KH,(p)及其2 對角塊矩陣的條件數(shù)，其中，p=2,3.總結(jié)如表1 所示，這里，Cond(A)表示給定矩陣A的 條件數(shù)，它可表述為，這里，λmax(A) 和λmin(A)分別為矩陣A的最大和最小特征值，A=KH,(p),p=2,3.

表1 分層剛度矩陣KH,(p)及其分塊矩陣的條件數(shù)(p=2,3)Tab.1 Condition number of the matrices KH,(p)，where p=2,3

表1 分層剛度矩陣KH,(p)及其分塊矩陣的條件數(shù)(p=2,3)Tab.1 Condition number of the matrices KH,(p)，where p=2,3

由此可見，對于分層高階元離散化矩陣KH,(p)，其條件數(shù)基本不隨網(wǎng)格規(guī)模及單元階次的增加而變化，但仍為病態(tài)矩陣.對于KH,(p)的兩個塊對角矩陣，一方面，由各單元頂點組成的Q4元矩陣雖然仍是病態(tài)矩陣，但對應(yīng)條件數(shù)比KH,(p)的條件數(shù)要小得多，可利用已有的高效方法如GAMG 法進(jìn)行求解；另一方面，由棱內(nèi)“虛節(jié)點”和面內(nèi)“虛節(jié)點”組成的矩陣則為一良態(tài)矩陣，其條件數(shù)基本不依賴于網(wǎng)格規(guī)模.

然后，將上述PCG 法應(yīng)用于不同網(wǎng)格規(guī)模下懸臂梁問題分層Q8元和分層Q12元離散系統(tǒng)(7)的求解.為了說明本文方法的優(yōu)越性，將利用工程中應(yīng)用非常廣泛的ILU(0)-PCG[19]和AGMG-FCG[20]進(jìn)行對比實驗，相應(yīng)的數(shù)值試驗結(jié)果總結(jié)如表2 和表3 所示，其中，“#element”表示單元數(shù)，“#dof”表示未知數(shù)總數(shù)，“Iter”表示迭代次數(shù)，“CPU/s”表示計算CPU 時間(單位為s).

表2 幾種PCG 法求解結(jié)構(gòu)網(wǎng)格下懸臂梁問題的迭代次數(shù)和計算CPU 時間Tab.2 Iteration counts and CPU time of several PCG methods for cantilever beam on structured meshes

表3 幾種PCG 法求解非結(jié)構(gòu)網(wǎng)格下懸臂梁問題的迭代次數(shù)和計算CPU 時間Tab.3 Iteration counts and CPU time of several PCG methods for cantilever beam on unstructured meshes

注2對懸臂梁問題，程序中采用的AMG 粗化算法在所有自由邊界上不是標(biāo)準(zhǔn)粗化，在這些節(jié)點上所構(gòu)造的插值算子不能重構(gòu)剛體模式(rigid body modes)，將會導(dǎo)致GAMG 法的收斂速度受到影響.因此，在實施算法中的步驟2 時，改為調(diào)用5 次GAMG法近似求解Q4元 離散系統(tǒng)Kvvzv=rv，以提高PCG法內(nèi)迭代計算效率.在程序?qū)崿F(xiàn)時，我們先將GAMG 迭代所需的各個要素一次性生成并保存好，所以整個迭代時間仍是很快的.

算例2(重力壩問題) 考慮某重力壩問題，其計算斷面尺寸如圖4 所示，作為平面應(yīng)變問題進(jìn)行有限元分析.假設(shè)壩體和壩基均為線彈性材料，其中，壩體的材料常數(shù)為E=20 GPa和 ν=0.167，密度為 γ=2 400 kg/m3，壩基的材料常數(shù)為E=20 GPa和 ν=0.25.邊界條件如圖4 所示，假設(shè)水位與壩頂齊平，計算時只考慮壩體自重和上游水壓力.

圖4 某重力壩截面尺寸及邊界約束(單位：m)Fig.4 Section sizes of a gravity dam and the boundary conditions (unit: m)

對區(qū)域進(jìn)行非結(jié)構(gòu)四邊形網(wǎng)格剖分，圖5 分別給出了2 類不同類型及規(guī)模的網(wǎng)格剖分圖.將上述PCG 法應(yīng)用于不同網(wǎng)格規(guī)模下重力壩問題分層Q8元 和分層Q12元離散系統(tǒng)(7)的求解，相應(yīng)的數(shù)值試驗結(jié)果總結(jié)如表4 所示.

表4 幾種PCG 法求解重力壩問題的迭代次數(shù)和計算CPU 時間Tab.4 Iteration counts and CPU time of several PCG methods for the gravity dam

圖5 2 類不同類型的重力壩網(wǎng)格剖分圖Fig.5 Two types of unstructured meshes with different elements for the gravity dam

算例3(腹拱壩問題) 考慮如圖6 所示的腹拱壩，其中壩體和壩基均采用線彈性模型，壩體的彈性常數(shù)為E=23 GPa 和v=0.17，密度為2 400 kg/m3，壩基的彈性常數(shù)為E=15 GPa 和v=0.17，壩基底部為縱向約束，兩側(cè)均為橫向約束.設(shè)上游水平面與壩頂齊平，計算時只考慮壩體自重和上游對壩體的水壓.

圖6 腹拱壩截面尺寸(單位：m)Fig.6 Section sizes of the arch-abdomen dam (unit: m)

對區(qū)域分別進(jìn)行非結(jié)構(gòu)擬一致四邊形網(wǎng)格剖分和局部加密四邊形網(wǎng)格剖分，圖7 分別給出了不同規(guī)模的兩類網(wǎng)格剖分圖.兩種PCG 法求解不同網(wǎng)格規(guī)模下腹拱壩問題分層Q8元 和分層Q12元離散系統(tǒng)(7)的數(shù)值結(jié)果總結(jié)如表5 所示.

表5 幾種PCG 法求解腹拱壩問題的迭代次數(shù)和計算CPU 時間Tab.5 Iteration counts and CPU time of several PCG methods for the arch-abdomen dam

圖7 2 類不同類型的腹拱壩網(wǎng)格剖分圖Fig.7 Two types of unstructured meshes with different elements for the arch-abdomen dam

由上述數(shù)值結(jié)果可見，我們設(shè)計的B-PCG 法對求解分層Q8元 和Q12元離散化系統(tǒng)均具有很好的計算效率和魯棒性，且隨著單元階次的增加和問題規(guī)模的不斷增加，B-PCG 法的優(yōu)勢更明顯.例如，對腹拱壩問題，當(dāng)單元數(shù)為1 277 時，3 種方法求解分層Q8元離散系統(tǒng)的迭代次數(shù)和計算CPU 時間分別為：ILU(0)-PCG 法為120 次和1.20 s，AGMG-FCG 法為95 次和5.49 s，而B-PCG 法僅為30 次和0.50 s；求解分層Q12元離散系統(tǒng)的迭代次數(shù)和計算CPU 時間分別為：ILU(0)-PCG 法為121 次和3.60 s，AGMG-FCG 法為91 次和52.74 s，而B-PCG 法僅為28 次和1.17 s.在設(shè)計兩水平PCG 法時，由于采用了分層有限元，不需構(gòu)造網(wǎng)格轉(zhuǎn)換算子和線性元矩陣.由于Kvv在PCG算法調(diào)用中保持不變，將基于Kvv的GAMG 法所需的各網(wǎng)格層信息，包括網(wǎng)格轉(zhuǎn)移算子、粗網(wǎng)格系數(shù)矩陣等一次性生成好，再充分利用GAMG 法在求解非結(jié)構(gòu)網(wǎng)格下Q4元離散系統(tǒng)所具的優(yōu)勢，能夠顯著提高內(nèi)迭代的計算效率以及相應(yīng)PCG 法的迭代效率.

3 結(jié)語

四邊形高階單元是能有效提升二維問題有限元分析計算精度的一類常用單元.實際應(yīng)用時，可通過適當(dāng)增加虛節(jié)點的方式增加單元形函數(shù)的次數(shù)，從而提高有限元近似解的計算精度.這些虛節(jié)點僅要求設(shè)置在棱內(nèi)部或面內(nèi)部即可，不需要具體的幾何坐標(biāo)，處理非常靈活，相應(yīng)的階譜函數(shù)也易于構(gòu)造，且其對應(yīng)的有限元矩陣具有分層結(jié)構(gòu)特性，易于為相應(yīng)離散代數(shù)化系統(tǒng)的求解設(shè)計高效求解器.

本文通過利用分層高階四邊形元離散矩陣的分層結(jié)構(gòu)特性及分塊對角矩陣的性質(zhì)，為相應(yīng)離散系統(tǒng)的求解提供了一種簡單有效的預(yù)條件子，所得到的兩水平分層PCG 法具有很好的計算效率.這種預(yù)條件子的構(gòu)造是通過將高次元離散系統(tǒng)方程的求解本質(zhì)性地化歸為Q4元離散系統(tǒng)的求解，再利用現(xiàn)有的高效GAMG 法求解Q4元離散系統(tǒng)，即可獲得內(nèi)迭代計算效率得到顯著提高的PCG 法.本文方法也可以推廣應(yīng)用于三維問題及非線性有限元數(shù)值計算中，或者其它實際工程計算問題的數(shù)值模擬中，以提高相應(yīng)有限元分析的整體效率.