李 響，張浩東，王 橋，周 偉

(武漢大學水資源與水電工程科學國家重點實驗室，武漢 430072)

0 引 言

邊界元法[1]作為有限元法的一個重要補充，對于某些有限元法不能很好解決的問題具有先天的優(yōu)勢。邊界元法只需要對求解問題的邊界進行離散，將求解問題的維度降低一維，并且由于采用了求解問題的基本解，它是一種半解析半數(shù)值的計算方法，因此能得到較高的計算精度。對于裂紋擴展問題，邊界元法只需要在裂紋尖端增加少量網(wǎng)格，極大地減少了工作量。對于一些復雜結構，例如含有細小孔洞的結構，含有細小水管的大壩結構和顆粒增強的復合材料結構等，采用邊界元法只需要對構成其幾何結構的邊界進行離散，減少了網(wǎng)格劃分的困難程度和計算自由度。由于采用了基本解，邊界元法非常適用于求解無限域和半無限域的問題，使得其在巖土工程和地下結構工程中的應用前景非常巨大。

對于齊次方程，邊界元法得到的積分方程只含有邊界積分，然而，對于非齊次方程，即右端項不為零的方程，邊界元法得到的積分方程除了含有邊界積分外，還含有域積分。如果利用劃分體網(wǎng)格的方法來計算域積分，邊界元法便喪失了只對邊界進行離散的巨大優(yōu)勢，因此邊界元法不適合求解非線性問題。而實際工程中的問題最終提煉得到的方程往往是非齊次方程，例如力學問題一般都會考慮重力，這便極大地限制了邊界元法在工程中的應用。

基于體網(wǎng)格的方法主要包括直接域積分法[2]，該方法需要對整個求解域進行離散，即需要劃分體單元，導致邊界類型算法失去了只對邊界進行離散的這一主要優(yōu)勢。對于一些較為復雜的結構，體單元的劃分也比較困難和耗時。

為了消除域積分需要的體網(wǎng)格，國內(nèi)外很多學者提出了各種將域內(nèi)積分轉化為邊界積分的方法[3]，常用的方法主要基于互易定理和散度定理?；诨ヒ锥ɡ淼姆椒ㄖ饕p互易法和多互易法等。雙互易法(dual reciprocity method, DRM)[4]是應用較為廣泛的一種方法，在一些無網(wǎng)格法[5]和邊界面法[6-8]中都有應用。該方法的基本思想是通過再一次的互易定理將域內(nèi)積分轉化為邊界積分。為了達到這一目的，可以將控制方程的非齊次項采用一系列函數(shù)插值，常用的插值函數(shù)有徑向基函(radial basic function, RBF)，需要在域內(nèi)布置節(jié)點，在計算過程中往往需要先求解一個線性方程組，該線性方程組的系數(shù)矩陣一般是一個病態(tài)滿陣，在布置的域內(nèi)節(jié)點不均勻或者節(jié)點數(shù)量很多的情況下，可能導致求解結果的不穩(wěn)定性[9-11]。

基于散度定理的方法是近年來才被提出的處理域積分的新方法，其具有代表性的方法為徑向積分法(radial integration method, RIM)[12,13]，該方法最近得到了較為廣泛的研究。該算法也是一種只需要對邊界進行離散的方法，其基本思想是通過散度定理將域積分轉化為包含一維徑向積分的邊界積分。該方法除了能夠將域積分轉化為邊界積分之外，還能夠非常有效的處理邊界元法中的一些奇異積分。基于散度定理的方法還包括直線積分法(line integration method, LIM)[14]，該方法非常適合與快速多極算法結合求解大規(guī)模問題[15,16]。

本文擬提出一種新的域積分計算方法----自適應背景網(wǎng)格積分法。該方法將域積分近似為背景網(wǎng)格積分，只需要對邊界進行單元離散，非常適合應用于邊界元法中。數(shù)值算例表明，本文提出的方法具有較高的計算精度。本文將該方法應用于二維彈性力學問題，并對大壩結構進行了仿真計算，取得了較好的結果。

1 自適應背景網(wǎng)格積分法

1.1 自適應背景網(wǎng)格

在自適應背景網(wǎng)格積分法中，只需要將求解域的邊界離散為邊界單元，通過邊界單元得到背景網(wǎng)格，在該過程中需要應用到二叉樹結構。在詳細介紹自適應背景網(wǎng)格的構造前，先定義如下參數(shù)：

網(wǎng)格大小：定義網(wǎng)格的大小為過正方形的邊長；單元大小：定義單元大小為單元的長度；單元等級：初始單元的等級設為1，每個單元可以進一步細分為兩個子單元，每個子單元的等級增加1，即等級為l+1的單元由等級為l的單元細分得到；預定義的最大網(wǎng)格大?。汗潭ㄖ?，為事先定義的值來控制葉子網(wǎng)格的大?。活A定義的最大單元等級：固定值，如果邊界單元的等級大于該預定值，停止對該單元的進一步細分。

有了以上定義，構造自適應背景網(wǎng)格的主要步驟如下：

(1)將計算域離散為邊界單元。邊界單元可以直接采用邊界元法中的單元。找到一個盡可能小的正方形覆蓋求解域，該正方形為四叉樹的根節(jié)點，稱為根網(wǎng)格，所有的邊界單元須在根網(wǎng)格之內(nèi)。

(2)細分網(wǎng)格。將每個網(wǎng)格細分為4個相同大小的子正方形，稱為子網(wǎng)格。如果父網(wǎng)格中的邊界單元的部分或全部在子網(wǎng)格內(nèi)，將該單元分配到該子網(wǎng)格內(nèi)。每個單元可以同時在多個網(wǎng)格內(nèi)(見圖1)。

圖1 一個單元在兩個網(wǎng)格中Fig.1 A unit element in two grids

(3)重新細分網(wǎng)格或網(wǎng)格內(nèi)的單元。首先得到每個網(wǎng)格內(nèi)每個單元大小，如果單元的大小不比網(wǎng)格小，細分該單元并將每個新單元的等級增加1(圖2)。重復該步驟，直到該網(wǎng)格內(nèi)每個單元的大小均小于該網(wǎng)格的大小。如果單元的大小比網(wǎng)格大小的一半還小，回到第二步，即細分網(wǎng)格(圖3)。該步驟可確保網(wǎng)格大小和該網(wǎng)格內(nèi)的所有單元的大小較為接近。

圖2 細分單元Fig.2 Subdivision unit element

圖3 細分網(wǎng)格Fig.3 Mesh refinement

(4)判斷網(wǎng)格位置(見圖4)。如果網(wǎng)格在計算域外，刪掉該網(wǎng)格，如果在計算域內(nèi)，進行(5)，否則，進行(6)。

(5)細分域內(nèi)網(wǎng)格。如果域內(nèi)網(wǎng)格大小大于域定義的最大網(wǎng)格大小，將該網(wǎng)格細分為四個相同大小的網(wǎng)格；重復該步驟直到網(wǎng)格大小小于域定義的最大網(wǎng)格大小。

(6)細分邊界網(wǎng)格。如果該邊界網(wǎng)格大小大于域定義的最大網(wǎng)格大小，回到(2)；否則，如果網(wǎng)格內(nèi)單元的最小等級大于預定義的最大單元等級，停止細分，否則，回到(2)。

通過以上步驟即可構造基于邊界單元的自適應背景網(wǎng)格，可以發(fā)現(xiàn)，該方法只需要邊界單元的信息，區(qū)別于其他域積分方法。邊界網(wǎng)格的大小通過邊界單元的大小進行控制，邊界單元大小越小或邊界單元等級越高，得到的邊界網(wǎng)格的大小便越小，而小的邊界網(wǎng)格可以更好地接近原計算域。圖4(b)是一個最終的網(wǎng)格細分例子。

圖4 自適應背景網(wǎng)格構造過程Fig.4 Process of adaptive cell-based domain mesh

在(4)中，需要判斷網(wǎng)格的位置。在該方法中，有3種網(wǎng)格：域外網(wǎng)格，邊界網(wǎng)格和域內(nèi)網(wǎng)格。如果網(wǎng)格中包含邊界單元，稱為邊界網(wǎng)格，該類型的網(wǎng)格容易區(qū)分。域外網(wǎng)格和域內(nèi)網(wǎng)格均不包含邊界單元，區(qū)分較為困難，涉及點相對于幾何體的位置問題。從以上步驟中，可以知道，在這種情況下該網(wǎng)格的父網(wǎng)格必然為邊界網(wǎng)格，即其包含邊界單元，通過這些邊界單元，便可判斷子網(wǎng)格在域內(nèi)還是域外。在父網(wǎng)格內(nèi)布置m個臨時點，這些臨時點分布于父網(wǎng)格內(nèi)的邊界單元上，那么可采用如下公式判斷點x的位置：

(1)

式中：yi和ni為臨時點的坐標和外法向量。

如果K>0，可認為x在域外，否則，可認為x在域內(nèi)。

應用式(1)，在子網(wǎng)格中布置N個點(見圖5)，于是可采用下式來判斷網(wǎng)格的位置。

(2)

式中：xj為網(wǎng)格中的N個點。

圖5 子網(wǎng)格位置Fig.5 Sub-grid position

1.2 背景網(wǎng)格積分法

考慮如下二維問題中的域積分：

(3)

式中：Ω為積分域；f(y)為被積函數(shù)；y{y1,y2}為二維空間中的一點。

采用背景網(wǎng)格積分時，邊界網(wǎng)格往往有一部分在域外，因此可定義如下新的方程：

(4)

式(3)便可采用背景網(wǎng)格積分近似為：

(5)

式中：Ck為第k個葉子網(wǎng)格。

式(5)可采用高斯積分法計算。

2 邊界元法求解包含體力的彈性力學問題

二維彈性力學問題的控制方程為如下幾個方程：

σij,j+bi=0

(6)

(7)

σij=Cijklεkl

(8)

式中：σij為應力張量分量；bi為體力分量；εij為應變張量分量；ui為位移分量；ui,j為位移分量的導數(shù)。

i,j,k,l均為1到2，且：

Cijkl=λδijδkl+μ(δikδjl+δilδjk)

(9)

式中：μ為剪切模量；λ=2vμ/(1-2v)為拉梅常數(shù)，v是泊松比。于是可以得到如下求解彈性力學問題的正則化積分方程：

(10)

式中：x和y分別為場點和源點；u和t分別位移和面力張量；U(x,y)和T(x,y)為基本解，對于平面應變問題，有：

(11)

(1-2v)(rkni-rink)}

(12)

式中：r(x,y)為x和y之間的距離；n為邊界上單位外法向量。

于是式(10)中的域積分可以寫為：

(13)

其中：

(14)

3 邊界元法在重力壩中的應用

本節(jié)給出了兩個數(shù)值算例，驗證了本文提出的方法的可行性和計算精度，并將該方法應用于重力壩的靜力分析中，計算結果與有限元法進行了對比。

3.1 方法精度驗證

該數(shù)值算例是一開孔的方板(圖6)，其尺寸見圖，采用無量綱單位，其體力為：

b1=ex2sinx1+ex2cosx1-

ex1sinx2-ex1cosx2

(15)

b2=ex1sinx2+ex2sinx1-

ex1cosx2-ex2cosx1

(16)

其中位移解析解為：

(17)

圖6 圓孔方板Fig.6 Square plate with round hole

設該問題為平面應變問題，泊松比v=0.25，彈性模量E=1.0。外邊界法向方向位移已知，切向方向和內(nèi)邊界面力已知，在邊界上總共布置400個單元。圖6和圖7為以幾何體中心為原點的極坐標下圓曲線r=1.5上的點的計算結果，橫坐標為計算點與坐標原點組成的向量相對x軸的角度。從圖中可以發(fā)現(xiàn)，計算結果與解析解吻合較好，具有較高的計算精度。

圖7 圓孔方板的位移計算結果Fig.7 The displacement calculation results of Fquare plate with round hole

圖8 圓孔方板的應力計算結果Fig.8 Stress calculation results of round hole square plate

3.2 重力壩計算

該算例考慮如圖9所示的重力壩，其具體尺寸見圖，單位為m，泊松比為v=0.2，彈性模量E=3.0 萬MPa，大壩單位重量為25×103N/m3，左邊水位高45 m，右邊水位高13 m，水的單位重量為9.8×103N/m3，大壩底部完全約束。

圖9 重力壩模型(單位：m)Fig.9 Gravity dam model

該問題可采用平面應變問題求解。由于缺乏解析解，本文與有限元軟件Abaqus的計算結果進行對比。圖10和圖11為直線x1=2.5上的位移u2應力σ2的計算結果，從圖中可以看出，位移計算結果與有限元非常接近，應力計算結果本文方法只需要少量的邊界單元便可得到較好的結果，而有限元法需要大量單元才能得到較好結果。

圖10 重力壩位移計算結果Fig.10 The displacement calculation results of gravity dam model

圖11 重力壩應力計算結果Fig.11 Stress calculation results of gravity dam model

4 結 語

本文提出了一種邊界元法中新的計算域積分的方法，稱為自適應背景網(wǎng)格積分法。該方法只需要對邊界進行單元離散，通過單元構造積分背景網(wǎng)格。在構造過程中需要引入樹結構，對于二維問題為四叉樹結構。本文將該方法應用于二維力學分析中，數(shù)值結果與解析解吻合較好。最后應用該方法對重力壩進行了仿真分析，并與有限元方法進行了對比，結果表明，本文的方法計算精度較高，采用較少的單元便可以獲得較高的精度。

該方法的計算復雜度為O(NM)，其中N為計算節(jié)點數(shù)量，M為背景網(wǎng)格數(shù)量，因此對于大規(guī)模問題計算量較大。為了克服這一問題，可采用快速算法，例如快速多極算法進行加速，該內(nèi)容正在研究過程中，成果將在后續(xù)工作中報道。

□

[1] 姚振漢,王海濤. 邊界元法[M]. 北京:高等教育出版社, 2010.

[2] Dallner R, Kuhn G. Efficient evaluation of volume integrals in the boundary element method[J]. Computer Methods in Applied Mechanics & Engineering, 1993,109(1-2):95-109.

[3] Hsiao S C, Mammoli A A, Ingber M S. The evaluation of domain integrals in complex multiply-connected three-dimensional geometries for boundary element methods[J]. Computational Mechanics, 2003,32(4):226-233.

[4] Nardini D, Brebbia C A. A new approach to free vibration analysis using boundary elements[J]. Applied Mathematical Modelling, 1983,7(3):157-162.

[5] Miao Y, Wang Q, Liao B, et al. A Dual Hybrid Boundary Node Method for 2D Elastodynamics Problems[J]. Computer Modeling in Engineering & Sciences, 2009,53(1):1-22.

[6] Zhou F, Zhang J, Sheng X, et al. Shape variable radial basis function and its application in dual reciprocity boundary face method[J]. Engineering Analysis with Boundary Elements, 2011,35(2):244-252.

[7] Zhang J, Qin X, Han X, et al. A boundary face method for potential problems in three dimensions[J]. International Journal for Numerical Methods in Engineering, 2010,80(3):320-337.

[8] Zhou F, Zhang J, Sheng X, et al. A dual reciprocity boundary face method for 3D non-homogeneous elasticity problems[J]. Engineering Analysis with Boundary Elements, 2012, 6(9):1 301-1 310.

[9] Schaback R. Error estimates and condition numbers for radial basis function interpolation[J]. Advances in Computational Mathematics, 1995,3(3):251-264.

[10] Narcowich F J, Ward J D, Wndland H. Refined error estimates for radial basis function interpolation[J]. Constructive Approximation, 2003,19(4):541-564.

[11] Brownlee R A. Error estimates for interpolation of rough data using the scattered shifts of a radial basis function[J]. Numerical Algorithms, 2005,39(1):57-68.

[12] Gao X W. The radial integration method for evaluation of domain integrals with boundary-only discretization[J]. Engineering Analysis with Boundary Elements, 2002,26(10):905-916.

[13] Gao X W. Evaluation of regular and singular domain integrals with boundary-only discretization-theory and Fortran code[J]. Journal of Computational & Applied Mathematics, 2005,175(2):265-290.

[14] Wang Q, Zhou W, Cheng Y, et al. Line integration method for treatment of domain integrals in 3D boundary element method for potential and elasticity problems[J]. Engineering Analysis with Boundary Elements, 2017,75:1-11.

[15] Wang Q, Zhou W, Cheng Y, et al. A line integration method for the treatment of 3D domain integrals and accelerated by the fast multipole method in the BEM[J]. Computational Mechanics, 2016:1-14.

[16] Wang Q, Zhou W, Cheng Y, et al. The Boundary Element Method with a Fast Multipole Accelerated Integration Technique for 3D Elastostatic Problems with Arbitrary Body Forces[J]. Journal of Scientific Computing, 2016:1-27.