等幾何分析方法求解靜電場非齊次邊值問題

張 勇 林 皋 胡志強(qiáng) 徐喜榮

(1.大連理工大學(xué)建設(shè)工程學(xué)部，遼寧 大連 116024； 2.大連理工大學(xué)電子信息與電氣工程學(xué)部，遼寧 大連 116024)

引 言

計(jì)算電磁學(xué)是電磁學(xué)分析的一種重要工具，其中對靜電場求解是最為基礎(chǔ)和重要的環(huán)節(jié)。計(jì)算電磁學(xué)的求解方法主要有解析法、半解析法和數(shù)值解法三大類。解析法和半解析法雖然能夠得到精確或半精確的解析式，但只適用于簡單的求解域、簡單的物理參數(shù)的問題，不具有通用性。數(shù)值解法雖然是在一定條件下的近似解法，但具有較強(qiáng)的通用性，最為常見的有有限差分法、有限元法、邊界元法、矩量法、無網(wǎng)格法[1-8]等。其中，有限元法是應(yīng)用最為廣泛的數(shù)值方法。各類數(shù)值方法都需要對求解域進(jìn)行二次離散或者建模，而離散之后的計(jì)算模型可能與原模型不具有一致性，從而喪失精度。

等幾何分析方法(IGA)是2005年由Hughes[10]率先提出的一種新型數(shù)值方法。該方法避免了傳統(tǒng)數(shù)值方法中求解模型與設(shè)計(jì)模型的非一致性，從而實(shí)現(xiàn)了將問題的分析計(jì)算構(gòu)架于精確模型上，提高了計(jì)算精度與效率。目前這種方法已被用于求解固體力學(xué)和流體力學(xué)問題[11-16]。最近，張勇和林皋等將該方法應(yīng)用于靜電場問題[17]，利用加權(quán)余量法推導(dǎo)了靜電場齊次邊值問題的等幾何分析方程，并實(shí)現(xiàn)了對自然邊界條件和齊次強(qiáng)制邊界條件的施加。另外，也將其應(yīng)用于波導(dǎo)本征問題[18-19]的求解, 但在波導(dǎo)本征問題中，無論是橫向磁場(TM)波還是橫向電場(TE)波，邊界條件都是齊次的。雖然等幾何分析方法具有良好的計(jì)算特性和較高的計(jì)算效率，但是由于非均勻有理B樣條(NURBS)基函數(shù)不具有插值性[20]，使得其求解帶有非齊次強(qiáng)制邊界條件的邊值問題變得復(fù)雜和困難。

自然邊界條件和強(qiáng)制邊界條件的施加方式是不同的，自然邊界條件最終轉(zhuǎn)化為關(guān)于邊界積分，從而對齊次或者非齊次邊界并不敏感。但對于強(qiáng)制邊界條件，需要給定離散自由度值，這對于具有插值特性的數(shù)值方法而言，是極為容易的，比如在有限元方法中，只需通過結(jié)點(diǎn)坐標(biāo)便可計(jì)算出邊界結(jié)點(diǎn)自由度值。但對于不具有插值特性的數(shù)值方法而言，這樣是無法做到的，比如在無網(wǎng)格方法，張淮清等[8]給出了基于拉格朗日乘子方法施加強(qiáng)制邊界條件，文獻(xiàn)[9]則給出了拉格朗日乘子、罰函數(shù)法、混合法等施加邊界條件，其中拉格朗日乘子最為簡潔。本文在文獻(xiàn)[17]基礎(chǔ)上將等幾何分析方法進(jìn)一步拓展到靜電場非齊次邊值問題。建立帶有拉格朗日乘子的泛函，同時(shí)對拉格朗日乘子采用等幾何離散，最終得到關(guān)于控制點(diǎn)電勢及與非齊次邊界相關(guān)的控制點(diǎn)拉格朗日乘子的增廣代數(shù)方程?？紤]到拉格朗日乘子僅賦存于非齊次強(qiáng)制邊界上，其增加的計(jì)算自由度相對于整體自由度而言較小，從而并不會影響等幾何分析方法的計(jì)算效率。

應(yīng)用基于拉格朗日乘子的等幾何分析方法求解靜電場非齊次邊值問題，以頂部帶有正弦分布電勢的方形板以及內(nèi)外壁均為非齊次邊值的圓環(huán)域?yàn)閿?shù)值算例。計(jì)算結(jié)果表明本文方法能很好處理靜電場非齊次邊值問題，通過與傳統(tǒng)方法和解析解的比較，顯示出其低自由度消耗、高效和高精度的特點(diǎn)。

1 基于拉格朗日乘子的靜電場等幾何分析方法推導(dǎo)

1.1 靜電場非齊次邊值問題

靜電場問題是電磁學(xué)分析中的一個(gè)典型邊值問題(BVP)，描述靜電場(域內(nèi)無電荷)特性的控制方程(φ為電勢)為

1.2 曲面的NURBS描述

NURBS是計(jì)算機(jī)輔助設(shè)計(jì)(CAD)中幾何形體的標(biāo)準(zhǔn)數(shù)學(xué)描述，是由NURBS基函數(shù)Ri,p,(i=0,1,…,n)，控制點(diǎn)Pi,(i=0,1,…,n)和節(jié)點(diǎn)向量Ξ={ξ0,ξ1…,ξm}三要素來定義的，其中，p為基函數(shù)次數(shù)，n為控制點(diǎn)數(shù)，m為節(jié)點(diǎn)數(shù)。NURBS基函數(shù)是通過B-樣條有理化得到的，如一維NURBS基函數(shù)可表示為

(3)

式中:ωi是對應(yīng)于Ni,p(ξ)的權(quán)值；Ni,p(ξ)為B-樣條基函數(shù)，可由Cox-de Boor公式[20]遞歸定義為式(4)：

(4)

NURBS基函數(shù)將節(jié)點(diǎn)向量所在的參數(shù)空間中的點(diǎn)映射到控制點(diǎn)所在的物理空間，即NURBS參數(shù)空間的一個(gè)參數(shù)區(qū)間[ξi,ξi+1)與物理空間中的一個(gè)單元Vi相對應(yīng)。

NURBS具有與B樣條一樣的性質(zhì)，如非負(fù)性、單位分解性、緊支性、靈活連續(xù)性、線性無關(guān)性和非插值性等。但與B樣條相比，它可以更準(zhǔn)確地描述圓錐曲線和曲面。當(dāng)權(quán)值相等時(shí)NURBS基函數(shù)則退化為B樣條基函數(shù)。

高維NURBS基函數(shù)可由低維NURBS基函數(shù)張量積得到，如式(5)所示的二維NURBS基函數(shù)形式，p、q分別為兩個(gè)方向基函數(shù)的次數(shù)。

(5)

一般地，NURBS曲面的向量形式可寫為

(6)

式中，Pi,j是NURBS曲面控制點(diǎn)。圖1給出了典型的NURBS曲面，圖1(a)為空間NURBS曲面，圖1(b)為平面內(nèi)NURBS曲面，后者適用于二維靜電場問題。

(a) 二次NURBS曲面

(b) 二次NURBS平面內(nèi)曲面圖1 NURBS曲面(控制點(diǎn)標(biāo)記為·)

1.3 求解域的等幾何離散與細(xì)分

在等幾何分析方法中，求解域的離散模型不是二次建模得到的，而是繼承于求解域的設(shè)計(jì)模型，即 CAD系統(tǒng)中的NURBS描述，如上節(jié)所述，對于二維問題，兩個(gè)張量積方向的節(jié)點(diǎn)向量可表示為Ξ={ξ0,ξ1，…,ξku-1},Η={η0,η1，…,ηkv-1}.參數(shù)域內(nèi)的節(jié)點(diǎn)區(qū)間[ξh,ξh+1)×[ηk,ηk+1)映射成求解域上的一個(gè)單元Vh,k。求解域NURBS單元Vh,k的任一點(diǎn)坐標(biāo)表示為

=NBh,k(ξ,η)xh,k

(7)

從CAD系統(tǒng)導(dǎo)入的初始模型一般并不能滿足解的精度需求，需要對求解域進(jìn)行細(xì)分。而在等幾何分析方法中，求解域的細(xì)分是非通信的、保形的，即細(xì)分的過程無需再與原來的設(shè)計(jì)模型進(jìn)行通信，并且細(xì)分前后的模型是幾何一致的，這是等幾何分析的重要特征，也是其與傳統(tǒng)有限元方法的重要區(qū)別。保形細(xì)分算法有h-細(xì)化方法，p-細(xì)化方法，k-細(xì)化方法等[10]。

1.4 電勢場的等參化

等參方法用于近似求解域的電勢場，求解域中任一點(diǎn)的電勢采用和式(7)相似的表示形式為：

(8)

式中，φi表示求解域控制點(diǎn)的電勢場變量。

1.5 等幾何分析方法離散方程

采用虛功原理[21]將控制微分方程式(1)及非齊次邊界條件式(2)表達(dá)成帶有拉格朗日乘子的泛函

(9)

對泛函Π(φ)取一次變分得

=0

(10)

式(10)可進(jìn)一步離散化為

=0

(11)

對引入的拉格朗日乘子采用NURBS等幾何離散，即

=NB(ξ,η)λ

(12)

將式(8)和(12)代入式(11)，注意到δφ和δλ的任意性，得到等幾何分析方法離散方程的形式為

(13)

(14)

(15)

(16)

2 靜電場非齊次邊值問題的算例分析

2.1 方形槽

通過分離變量法，可得該問題的電勢及電場強(qiáng)度解析解如式(17)所示：

(17)

數(shù)值計(jì)算時(shí)，上述問題簡化為平面問題，如圖2所示，其中紅色點(diǎn)為等幾何分析方法的初始網(wǎng)格控制點(diǎn)。

圖2 方形槽非齊次邊值問題及IGA初始網(wǎng)格示意圖

等幾何分析方法的初始網(wǎng)格只有一個(gè)等幾何單元，不能滿足計(jì)算的精度要求，利用h-細(xì)化方法，可以得到較細(xì)的網(wǎng)格，即網(wǎng)格A(含有36個(gè)控制點(diǎn)自由度和6個(gè)拉格朗日乘子自由度)和網(wǎng)格B(含有100個(gè)控制點(diǎn)自由度和10個(gè)拉格朗日乘子自由度)。表1為本文方法與文獻(xiàn)[8]列出的有限差分法(FDM)、徑向基無網(wǎng)格法(RBF)計(jì)算電勢最大值誤差、相對均方根誤差對比。本文方法中場值計(jì)算點(diǎn)為65×65個(gè)，在方形槽內(nèi)均勻分布。相對均方根誤差ERMS和最大誤差Emax計(jì)算公式為：

(18)

Emax=max(abs(φi_exact-φi_calc))

(19)

式中：φi_exact為解析解；φi_calc為數(shù)值解。

表1數(shù)據(jù)表明：網(wǎng)格A的結(jié)果和有限差分法的精度相當(dāng)，網(wǎng)格 B的結(jié)果和無網(wǎng)格法接近，這表明本文采用拉格朗日乘子處理等幾何分析方法中的非齊次邊值問題是有效的。對于電場強(qiáng)度，由于文獻(xiàn)[8]中沒有給出計(jì)算結(jié)果和誤差，故直接給出本文方法計(jì)算的電場強(qiáng)度相對均方根誤差和最大誤差，見表2.與電勢數(shù)值解的精度相比，電場強(qiáng)度數(shù)值解的精度稍低，這是源于電場強(qiáng)度是電勢的梯度場相關(guān)，也就是它的數(shù)值精度要比電勢低。但網(wǎng)格A到網(wǎng)格B的誤差減小的程度表明，本文方法的數(shù)值解有較快收斂速度。表1和表2中的Emax和ERMS就是最大誤差和相對均方根誤差，計(jì)算式分別為式(18)和(19)。

表1 電勢計(jì)算值的最大值誤差和均方根誤差列表

表2 電場強(qiáng)度計(jì)算值的最大值誤差和均方根誤差列表

上述比較初步確定了本文方法的有效性，為了進(jìn)一步比較其收斂性和收斂速度，采用五種尺寸的單元劃分，對于每一種尺寸的網(wǎng)格，分別對應(yīng)一個(gè)等幾何單元模型和一個(gè)傳統(tǒng)的有限元單元模型。計(jì)算自由度數(shù)量和單元數(shù)量在表3中列出，其中，對于等幾何分析方法，其計(jì)算自由度數(shù)為N+M，N為控制點(diǎn)自由度數(shù)，M為非齊次邊界附件自由度數(shù)，即拉格朗日乘子離散自由度數(shù)。在圖3中用雙對數(shù)坐標(biāo)表示自由度數(shù)和單元數(shù)的關(guān)系，從圖3可以看出等幾何分析方法的單位單元的自由度消耗要比有限元少，并且隨著網(wǎng)格的加密，這種差距在增大。對于非齊次邊值問題，需要附加自由度，但其所占的比重并不大。

表3 各種網(wǎng)格下的單元數(shù)量和自由度數(shù)量

圖3 計(jì)算自由度與單元數(shù)量圖

(a) 電勢相對均方根誤差

(b) x方向電場強(qiáng)度相對均方根誤差

(c) y方向電場強(qiáng)度相對均方根誤差圖4 域內(nèi)電勢和電場強(qiáng)度均方根誤差收斂圖

對每個(gè)網(wǎng)格得出65×65個(gè)計(jì)算點(diǎn)上電勢和電場強(qiáng)度值，并根據(jù)公式(18)和(19)計(jì)算其與解析解(17)的相對均方根誤差。圖4給出了這些誤差隨網(wǎng)格細(xì)化的收斂，圖4(a)、(b)、(c)分別代表電勢、x方向和y方向電場強(qiáng)度的相對均方根誤差收斂圖。從整體上看，電勢場的誤差收斂較快(在網(wǎng)格5中為10-7數(shù)量級)，電場強(qiáng)度收斂較慢(在網(wǎng)格5中為10-4數(shù)量級)；從電勢或者電場強(qiáng)度收斂斜率看，等幾何分析方法的收斂斜率均比有限元大，也就是說等幾何分析方法的收斂速度快。

最后給出整個(gè)域上的等幾何分析方法在網(wǎng)格5中的電勢數(shù)值解分布圖，如圖5(看1063頁)所示，(a)與(b)分別是數(shù)值解及其和解析解的相對誤差分布圖，圖5(b)中的相對誤差在10-6數(shù)量級，表明其數(shù)值解精度高。另外從誤差在域內(nèi)的分布可以看出，非齊次邊界上解的精度要比域內(nèi)解低，這表明對于非齊次邊值，其上拉格朗日乘子的離散數(shù)目不能過少。

從這個(gè)簡單的靜電場非齊次邊值問題可以看出，等幾何分析方法的優(yōu)勢是固有的等幾何性，即計(jì)算模型與設(shè)計(jì)模型一致，以及細(xì)化非通信和保形特點(diǎn)。但是由于其基函數(shù)的非插值性，對于非齊次邊界需要額外的附加自由度，所幸的是邊界自由度數(shù)是O(N)，而域內(nèi)控制點(diǎn)的自由度數(shù)是O(N2)，其中N為一個(gè)方向上自由度的度量，即附加自由度占計(jì)算自由度的比重較小，特別是對于細(xì)網(wǎng)格劃分情況。此外，雖然增加了計(jì)算自由度，但和傳統(tǒng)方法相比，等幾何分析方法在計(jì)算效率上具有單位單元的計(jì)算自由度消耗小、計(jì)算精度高、收斂速度快的特點(diǎn)。

2.2 同心圓柱環(huán)面

電位和電場強(qiáng)度的解析解可由極坐標(biāo)下的分離變量法求得

k(Crk+Dr-k)sin(kθ)eθ

(20)

式中系數(shù)A、B、C、D可根據(jù)邊界條件計(jì)算得出。在計(jì)算中，取Ra=0.5 m，Rb=1.0 m，φa1=1.75 V，φa2=0.5 V，φb1=-0.5 V，φb2=-0.5 V，k=2.

圖6 圓柱環(huán)面非齊次邊值問題及初始網(wǎng)格圖

等幾何分析的初始模型如圖6所示，紅色點(diǎn)表示等幾何模型的控制點(diǎn)，其不在域內(nèi)體現(xiàn)了控制點(diǎn)的非插值性。對于圓柱環(huán)面非齊次邊值問題，其內(nèi)壁和外壁均為非齊次邊界條件，比算例1的齊次和非齊次混合邊值條件要簡單，但本算例模型的邊界是更為復(fù)雜的曲線形式。采用六種尺寸的單元劃分，對于每一種尺寸的網(wǎng)格，分別對應(yīng)一個(gè)等幾何單元模型和一個(gè)傳統(tǒng)的有限元單元模型，各種網(wǎng)格下的單元數(shù)量和自由度數(shù)量如表4所列，對于網(wǎng)格6，有限元模型的自由度數(shù)已經(jīng)超過10 000，而等幾何分析方法只有不到5 000.

表4 各種網(wǎng)格下的單元數(shù)量和自由度數(shù)量

對每個(gè)網(wǎng)格計(jì)算出260(環(huán)向)×65(徑向)個(gè)計(jì)算點(diǎn)上電勢和電場強(qiáng)度值，并根據(jù)式(18)和(19)計(jì)算其與解析解(20)的相對均方根誤差。在圖7(a)、(b)、(c)中分別顯示了電勢和電場強(qiáng)度的相對均方根誤差隨計(jì)算自由度的收斂性，從中可以看出在網(wǎng)格較粗(網(wǎng)格1)的情況下，由于非齊次邊界上附加自由度占的比重較大，使得等幾何分析方法比有限元要消耗更多的自由度數(shù)，但從圖中誤差-自由度的斜率可以看出，隨著網(wǎng)格加密，等幾何分析方法的收斂較快。需指出是圖7為雙對數(shù)坐標(biāo)軸，雖然兩條線的距離不是很大，但體現(xiàn)的差別為倍數(shù)關(guān)系。

圖8(看1064頁)顯示了等幾何分析方法在網(wǎng)格6中的電勢數(shù)值解分布，(a)與(b)分別是數(shù)值解及其和解析解的相對誤差分布圖，圖8(b)中的相對誤差在10-5數(shù)量級，與算例1相比，其精度有所下降，這說明求解域的形狀及非齊次邊界的數(shù)量對結(jié)果的精度有一定的影響。

(a) 電勢相對均方根誤差

(b) x方向電場強(qiáng)度相對均方根誤差

(c) y方向電場強(qiáng)度相對均方根誤差圖7 域內(nèi)電勢和電場強(qiáng)度均方根誤差收斂圖

3 結(jié) 論

本文將等幾何分析方法擴(kuò)展到靜電場非齊次邊值問題，在虛功方程中引入拉格朗日乘子，推導(dǎo)出靜電場非齊次邊值問題的等幾何分析方法求解式。修正后的等幾何分析方法雖然增加了計(jì)算自由度(非齊次邊界上的離散拉格朗日乘子)，但與有限元方法相比，該方法仍然具有計(jì)算自由度少、收斂速度快、計(jì)算精度和計(jì)算效率高的特點(diǎn)，并且非通訊保形細(xì)化方法使得模型的細(xì)分可以自動進(jìn)行，更為簡潔與方便。應(yīng)用此方法求解了方形槽與同心圓環(huán)的靜電場非齊次邊值問題，結(jié)果表明該方法能夠有效的求解非齊次邊值問題，且具有上述的優(yōu)點(diǎn)，可進(jìn)一步在計(jì)算電磁學(xué)中推廣應(yīng)用。

[1] 周 勤,茅乃豐.旋轉(zhuǎn)對稱靜電場的一種新的數(shù)值解法[J].電波科學(xué)學(xué)報(bào),1990,5(4):26-34.

ZHOU Qin,MAO Naifeng.The electrostatic field problems with representations of spline functions and their quasi-charge density analysis method[J].Chinese Journal of Radio Science,1990,5(4):26-34.(in Chinese)

[2] SONG B,FU J.Modified indirect boundary element technique and its application to electromagnetic potential problems[J].IEE proceeding-H of Microwaves Antennas and Propagation,1992,139(3):292-296.

[3] 金建銘.電磁場有限元方法[M].西安電子科技大學(xué)出版社,1998:51-70.

[4] BAMJI S S,BULINSKI A T.Electric field calculations with the boundary element method [J].IEEE Transactions on Electrical Insulation,1993,28(3):420-424.

[5] 甄蜀春,曹 蕾,張繼龍.波動方程差分解法對波導(dǎo)螺釘調(diào)配器的分析[J].電波科學(xué)學(xué)報(bào),2005,20(1):125-127.

ZHEN Shuchun,CAO Lei,ZHANG Jilong.Analysis of screw tuner based on wave equation finite difference method[J].Chinese Journal of Radio Science,2005,20(1):125-127.( in Chinese)

[6] 汪朝暉,廖振方,陳德淑.有限元法分析尖板電極結(jié)構(gòu)的空間靜電場分布[J].重慶大學(xué)學(xué)報(bào),2010,33(5):41-47.

WANG Zhaohui,LIAO Zhenfang,CHEN Deshu.Analysis of spatial electric field with point-plate electrodes configuration using finite element method[J].Journal of Chongqing University,2010,33(5):41-47.(in Chinese)

[7] 梁志偉,趙國偉,徐 杰,等.柱形等離子體天線輻射特性的矩量法分析[J].電波科學(xué)學(xué)報(bào),2008,23(4):749-753.

LIANG Zhiwei,ZHAO Guowei,XU Jie,et al.Analysis of plasma-column antenna using moment method[J].Chinese Journal of Radio Science,2008,23(4):749-753.(in Chinese)

[8] 張淮清.電磁場計(jì)算中的徑向基函數(shù)無網(wǎng)格法研究[D].重慶：重慶大學(xué),2008.

ZHANG Huaiqing.Research on Radial Basis Function Meshless Method in Numercial Computation of Electromagentic Field[D].Chongqing： Chongqing University,2008.(in Chinese)

[9] FERNANDEZ-MENDEZ S,HUERTA A.Imposing essential boundary conditions in mesh-free methods[J].Comput Methods Appl Mech Engrg,2004,193(12/14):1257-1275.

[10] HUGHES T J R,COTTRELL J A and BAZILEVS Y.Isogeometric analysis CAD,finite elements,NURBS,exact geometry and mesh refinement[J].Comput Methods Appl Mech Engrg,2005,194(39/41):4135/4195.

[11] REALI A.An isogeometric analysis approach for the study of structural vibrations[J].Journal of Earthquake Engineering,2006,10(Special Issue 1):1-30.

[12] BAZILEVS Y,CALO V M,ZHANG Y,et al.Isogeometric fluid-structure interaction analysis with applications to arterial blood flow[J].Computational Mechanics,2006,38(4/5):310-322.

[13] COTTRELL J A,REALI A,BAZILEVS Y,et al.Isogeometric analysis of structural vibrations[J].Computer Methods in Applied Mechanics and Engineering,2006,195(41/43):5257-5296.

[14] COTTRELL J A,HUGHES T J R,REALI A.Studies of refinement and continuity in isogeometric structural analysis[J].Computer Methods in Applied Mechanics and Engineering,2007,196(41/44):4160-4183.

[15] ZHANG Y,BAZILEVS Y,GOSWAMI S,et al.Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow[J].Computer Methods in Applied Mechanics and Engineering,2007,196(31/32):2943-2959.

[16] ZHANG Y,LIN G,HU Z Q.Isogeometric analysis based on scaled boundary finite element method[J].IOP Conference Series:Materials Science and Engineering,2010,10(1):012237.

[17] 張 勇,林 皋,劉 俊,等.等幾何分析方法及其應(yīng)用于偏心柱面靜電場問題[J].電波科學(xué)學(xué)報(bào),2012,27(1):177-183.

ZHANG Yong,LIN Gao,LIU Jun,et al.The isogeometric analysis and its application to eccentrically parallel cylinders electrostatic problem[J].Chinese Journal of Radio Science,2012,27(1):177-183.(in Chinese)

[18] 張 勇,林 皋,劉 俊,等.波導(dǎo)本征問題的等幾何分析方法[J].應(yīng)用力學(xué)學(xué)報(bào),2012,29(2):113-119.

ZHANG Yong,LIN Gao,LIU Jun,et al.The isogeometric analysis for eigen problem of waveguide[J].Chinese Journal of Applied Mechanics,2012,29(2):113-119.(in Chinese)

[19] 張 勇,林 皋,胡志強(qiáng),等.基于等幾何分析方法求解任意截面波導(dǎo)本征問題[J].計(jì)算物理,2012,29(4):9-17.

ZHANG Yong,LIN Gao,HU Zhiqiang,et al.Eigenvalue Analysis of Waveguide with Arbitrary Cross-section by Isogeometric Analysis[J].Chinese Journal of Computational Physics,2012,29(4):9-17.(in Chinese)

[20] PIEGL L,TILLER W.The NURBS Book [M].Berlin:Springer Verlag,1997.

[21] HUGHES T J R.The Finite Element Method,Linear Static and Dynamic Finite Element Analysis[M].New York:Dover,2000.