多重網(wǎng)格求解動(dòng)力方程方法研究

楊 強(qiáng)，黨亞民，趙文嬌,2，杜 彬,2

(1.中國(guó)測(cè)繪科學(xué)研究院，北京 100830；2.山東科技大學(xué)，山東 青島 266590)

多重網(wǎng)格求解動(dòng)力方程方法研究

楊 強(qiáng)1，黨亞民1，趙文嬌1,2，杜 彬1,2

(1.中國(guó)測(cè)繪科學(xué)研究院，北京 100830；2.山東科技大學(xué)，山東 青島 266590)

在求解動(dòng)力方程過(guò)程中，引入多重網(wǎng)格方法，通過(guò)多次迭代、逐漸細(xì)分網(wǎng)格的方法簡(jiǎn)化有限元模型的計(jì)算，并進(jìn)行實(shí)例分析。研究結(jié)果表明，多重網(wǎng)格方法運(yùn)用靈活、減少運(yùn)算量、計(jì)算簡(jiǎn)便且具有降噪的優(yōu)勢(shì)，這種方法的引入，為地殼運(yùn)動(dòng)和形變的研究提供了一個(gè)有用的工具。

多重網(wǎng)格；有限元；動(dòng)力方程；迭代方法

利用大地測(cè)量技術(shù)研究地表形變特征，并進(jìn)一步確定地殼的應(yīng)變特征是大地測(cè)量和地球動(dòng)力學(xué)主要研究?jī)?nèi)容之一。而地殼運(yùn)動(dòng)作為一個(gè)動(dòng)態(tài)演化過(guò)程，其位移、應(yīng)力隨時(shí)間不斷演化發(fā)展。地殼應(yīng)力、應(yīng)變場(chǎng)的動(dòng)態(tài)演化過(guò)程，體現(xiàn)了地殼內(nèi)部各物理因素的變化與影響程度[1-2]。

一個(gè)動(dòng)態(tài)系統(tǒng)的動(dòng)力特性主要包括：幾何結(jié)構(gòu)、邊界條件和材料性質(zhì)，其解析表達(dá)式則是將其質(zhì)量、剛度和阻尼分布分別利用質(zhì)量矩陣、剛度矩陣和粘滯矩陣表示出來(lái)，以此來(lái)確定系統(tǒng)運(yùn)動(dòng)的動(dòng)力特性。解算方法很多，比較典型的方法是有限元方法[3]。

傳統(tǒng)有限元方法的離散化過(guò)程和數(shù)值求解過(guò)程是相互獨(dú)立、互無(wú)作用的，其離散化過(guò)程不能預(yù)測(cè)合適的解，往往造成計(jì)算上的浪費(fèi)。而將網(wǎng)格劃分很細(xì)，又使得代數(shù)方程組階數(shù)過(guò)大，而精度沒(méi)有提高。為了克服這種缺陷，本文將多重網(wǎng)格理論與有限元方法相結(jié)合，提出一種有效求解地殼動(dòng)力方程的方法。

1 多重網(wǎng)格求解過(guò)程

多重網(wǎng)格方法的主要特點(diǎn)就是在粗網(wǎng)格層上迭代求解其誤差方程直到誤差收斂的基礎(chǔ)上，通過(guò)插值的方法將誤差校正結(jié)果返回到細(xì)網(wǎng)格層上進(jìn)行校正。通過(guò)這種迭代、校正的過(guò)程，減少了計(jì)算量、提高了解算精度，在一定程度上彌補(bǔ)了傳統(tǒng)有限元方法的不足[4-5]。

動(dòng)力響應(yīng)分析是一個(gè)多自由度系統(tǒng)問(wèn)題，其動(dòng)力方程可以表示為

(1)

1.1 初始網(wǎng)格位移解算

利用有限元方法，首先求解初始網(wǎng)格各節(jié)點(diǎn)位移uk-1。初始網(wǎng)格劃分?jǐn)?shù)量較少，稱(chēng)為粗網(wǎng)格。

1.2 定義向下延拓雙線性插值算子

通過(guò)雙線性插值算子,獲取k層網(wǎng)格下的位移向量uk。定義雙線性插值算子

(2)

1.3 前光滑過(guò)程

細(xì)網(wǎng)格下的方程組表示為

{Kk}u(t+Δt)={Qt+Δt}.

(3)

1.4 粗網(wǎng)格修正過(guò)程

}.

(4)

2)限制殘量。利用限制算子以減少誤差影響，通過(guò)對(duì)細(xì)網(wǎng)格的限制

}.

(5)

(6)

3)誤差求解。對(duì)于誤差，同樣需要進(jìn)行解算、插值、迭代的計(jì)算過(guò)程。通過(guò)求解粗網(wǎng)格下的方程組

{Kk-1}{vk-1}={dk-1}

(7)

獲取粗網(wǎng)格下解的誤差vk-1。誤差的求解同樣也采用高斯-賽德?tīng)柕椒ā?/p>

(8)

5)細(xì)網(wǎng)格校正。最后對(duì)k層網(wǎng)格進(jìn)行精確校正

(9)

1.5 后光滑過(guò)程

綜上所述，多重網(wǎng)格算法可描述如下：

1)計(jì)算粗網(wǎng)格初始位移;

2)插值計(jì)算細(xì)網(wǎng)格位移初值;

3)迭代解算細(xì)網(wǎng)格位移向量;

4)位移修正;

5)后光滑消除高頻誤差。

為使結(jié)果更加真實(shí)可靠，初始網(wǎng)格位移必須相對(duì)精確，本文使用長(zhǎng)期觀測(cè)的GPS觀測(cè)站觀測(cè)結(jié)果檢核初始網(wǎng)格節(jié)點(diǎn)位移。

2 算例分析

本文利用Fortran語(yǔ)言編程，通過(guò)算例進(jìn)行分析。

算例：如圖1所示，結(jié)構(gòu)上部y方向固定和右部上半部x方向固定，每單元橫向跨度為5 m，縱向?yàn)? m，厚度為0.5 m，介質(zhì)彈性模量為90 GPa，泊松比為0.25，密度為2 750 kg/m3，粘滯系數(shù)為5.0×1021Pa·s，左側(cè)加載向右荷載P=10 MPa。加載時(shí)間長(zhǎng)度為10 s，步長(zhǎng)Δt=0.004 s。本文采用彈性模型進(jìn)行計(jì)算，A點(diǎn)假設(shè)當(dāng)時(shí)間為t=2 s時(shí)發(fā)生斷裂，斷裂帶發(fā)生右旋走滑，滑動(dòng)量為2.5e-4 m。

圖1 算例示意圖

圖2 有限元直接解算位移場(chǎng)和最大主應(yīng)力場(chǎng)

本文利用Fortran語(yǔ)言編程，分別利用普通有限元法和多重網(wǎng)格有限元方法對(duì)模型進(jìn)行計(jì)算，得到位移場(chǎng)和應(yīng)力場(chǎng)，并進(jìn)行對(duì)比。

模型計(jì)算得到的位移場(chǎng)及應(yīng)力場(chǎng)如圖2～4所示。在多重網(wǎng)格計(jì)算過(guò)程中，首先將結(jié)構(gòu)劃分為較大單元，利用有限元法直接解算位移，將其作為初始網(wǎng)格位移。為更為細(xì)致研究結(jié)構(gòu)各部分對(duì)荷載的動(dòng)力響應(yīng)，利用多重網(wǎng)格方法將各個(gè)單元進(jìn)一步細(xì)分。

以圖1中A點(diǎn)為例，以結(jié)構(gòu)力學(xué)解析解作為真值[6]，比較有限元直接解算結(jié)果與多重網(wǎng)格解算結(jié)果見(jiàn)表1。

表1 結(jié)果比較(Δt=0.004 s)

圖3 多重網(wǎng)格初始網(wǎng)格位移場(chǎng)和最大主應(yīng)力場(chǎng)

從表1中可以看出，有限元直接解誤差并不穩(wěn)定，而多重網(wǎng)格解誤差隨著時(shí)間的增加逐漸降低。因此，使用多重網(wǎng)格求解動(dòng)力模型有其精度上的優(yōu)勢(shì)。

3 結(jié) 論

通過(guò)以上研究，本文得出以下結(jié)論：

1)針對(duì)有限元方法運(yùn)算量較大，大量網(wǎng)格往往需要專(zhuān)業(yè)軟件進(jìn)行解算的問(wèn)題，本文提出了多重網(wǎng)格解算方法，以數(shù)量較少的初始網(wǎng)格直接采用有限元解算，然后通過(guò)細(xì)分網(wǎng)格，采用迭代方法就可以實(shí)現(xiàn)大量網(wǎng)格的解算，大大簡(jiǎn)化了運(yùn)算量，可以自主編程實(shí)現(xiàn)；

2)隨著時(shí)間的增加，多重網(wǎng)格精度逐漸提高，而有限元解則沒(méi)有明顯提高，顯示出多重網(wǎng)格具有降噪、提高精度的優(yōu)勢(shì)；

3)實(shí)際解算中，多重網(wǎng)格的運(yùn)用可以非常靈活，既可以對(duì)整個(gè)結(jié)構(gòu)進(jìn)行細(xì)分，也可以對(duì)某一大單元進(jìn)行細(xì)分，計(jì)算簡(jiǎn)便且精度較高。

綜上所述，多重網(wǎng)格方法的引入，為地殼運(yùn)動(dòng)和形變的研究提供了一個(gè)有用的工具。

[1]黨亞民，晁定波，許才軍，等.利用GPS形變資料確定地殼形變的應(yīng)變特征[J].測(cè)繪科學(xué)，2001，26(3)：10-13.

[2]黨亞民，陳俊勇，張燕平，等.利用GPS資料分析南天山地區(qū)的地殼形變特征[J].測(cè)繪科學(xué)，2002，27(4)：13-15.

[3]張東寧，許忠淮.中國(guó)大陸巖石層動(dòng)力學(xué)數(shù)值模型的邊界條件[J].地震學(xué)報(bào)，1999，21(2)：133-139.

[4]ZHANG J.Accelerated multigrid high accuracy solution of the convection diffusion equation with high Reynolds number[J].Numerical Methods for Partial Differential Equations，1997，13:77-92.

[5]張淼.結(jié)構(gòu)動(dòng)力響應(yīng)分析多重網(wǎng)格方法研究[D].長(zhǎng)春:吉林大學(xué)，2008.

[6]張森文，曹開(kāi)彬．計(jì)算結(jié)構(gòu)動(dòng)力響應(yīng)的狀態(tài)方程直接積分法[J].計(jì)算力學(xué)學(xué)報(bào)，2000，17(1)：94-97.

[7]楊強(qiáng)，黨亞民.利用GPS速度場(chǎng)估算青藏高原地殼韌性層等效粘滯系數(shù)分布的研究[J].測(cè)繪學(xué)報(bào)，2010,39(5)：497-502.

[責(zé)任編輯：劉文霞]

Researchonthesolutiontothedynamicequationwithmultigridmethod

YANG Qiang1,DANG Ya-min1,ZHAO Wen-jiao1,2,DU Bin1,2

(1.Chinese Academy of Surveying and Mapping，Beijing 100830,China; 2.Shandong University of Science and Technology，Qingdao 266590，China)

The method is introduced. Multigri mothod is used during solving the dynamic equations，of which the several times iterative and gradual grid subdivision are used to simplify the calculation of FEM.Its case analysis shows the multigrid method has the advantages of using flexibly，reducing operand，simplifying calculation and restraining noise.The introduction of multi-grid method provides a useful tool for the study of crustal movement and deformation.

multigrid；finite element；dynamic equation；iterative method

2012-09-11

楊 強(qiáng)(1975-)，男，博士.

P227

：A

：1006-7949(2014)01-0017-04