多重網格求解動力方程方法研究

楊 強，黨亞民，趙文嬌,2，杜 彬,2

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

多重網格求解動力方程方法研究

楊 強1，黨亞民1，趙文嬌1,2，杜 彬1,2

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

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

多重網格；有限元；動力方程；迭代方法

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

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

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

1 多重網格求解過程

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

動力響應分析是一個多自由度系統(tǒng)問題，其動力方程可以表示為

(1)

1.1 初始網格位移解算

利用有限元方法，首先求解初始網格各節(jié)點位移uk-1。初始網格劃分數(shù)量較少，稱為粗網格。

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

通過雙線性插值算子,獲取k層網格下的位移向量uk。定義雙線性插值算子

(2)

1.3 前光滑過程

細網格下的方程組表示為

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

(3)

1.4 粗網格修正過程

}.

(4)

2)限制殘量。利用限制算子以減少誤差影響，通過對細網格的限制

}.

(5)

(6)

3)誤差求解。對于誤差，同樣需要進行解算、插值、迭代的計算過程。通過求解粗網格下的方程組

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

(7)

獲取粗網格下解的誤差vk-1。誤差的求解同樣也采用高斯-賽德爾迭代方法。

(8)

5)細網格校正。最后對k層網格進行精確校正

(9)

1.5 后光滑過程

綜上所述，多重網格算法可描述如下：

1)計算粗網格初始位移;

2)插值計算細網格位移初值;

3)迭代解算細網格位移向量;

4)位移修正;

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

為使結果更加真實可靠，初始網格位移必須相對精確，本文使用長期觀測的GPS觀測站觀測結果檢核初始網格節(jié)點位移。

2 算例分析

本文利用Fortran語言編程，通過算例進行分析。

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

圖1 算例示意圖

圖2 有限元直接解算位移場和最大主應力場

本文利用Fortran語言編程，分別利用普通有限元法和多重網格有限元方法對模型進行計算，得到位移場和應力場，并進行對比。

模型計算得到的位移場及應力場如圖2～4所示。在多重網格計算過程中，首先將結構劃分為較大單元，利用有限元法直接解算位移，將其作為初始網格位移。為更為細致研究結構各部分對荷載的動力響應，利用多重網格方法將各個單元進一步細分。

以圖1中A點為例，以結構力學解析解作為真值[6]，比較有限元直接解算結果與多重網格解算結果見表1。

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

圖3 多重網格初始網格位移場和最大主應力場

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

3 結 論

通過以上研究，本文得出以下結論：

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

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

3)實際解算中，多重網格的運用可以非常靈活，既可以對整個結構進行細分，也可以對某一大單元進行細分，計算簡便且精度較高。

綜上所述，多重網格方法的引入，為地殼運動和形變的研究提供了一個有用的工具。

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

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

[3]張東寧，許忠淮.中國大陸巖石層動力學數(shù)值模型的邊界條件[J].地震學報，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]張淼.結構動力響應分析多重網格方法研究[D].長春:吉林大學，2008.

[6]張森文，曹開彬．計算結構動力響應的狀態(tài)方程直接積分法[J].計算力學學報，2000，17(1)：94-97.

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

[責任編輯：劉文霞]

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

楊 強(1975-)，男，博士.

P227

：A

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