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

(武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢 430072)

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 細(xì)分單元Fig.2 Subdivision unit element

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

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

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

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

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

圖4 自適應(yīng)背景網(wǎng)格構(gòu)造過程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ū)分較為困難，涉及點(diǎn)相對于幾何體的位置問題。從以上步驟中，可以知道，在這種情況下該網(wǎng)格的父網(wǎng)格必然為邊界網(wǎng)格，即其包含邊界單元，通過這些邊界單元，便可判斷子網(wǎng)格在域內(nèi)還是域外。在父網(wǎng)格內(nèi)布置m個(gè)臨時(shí)點(diǎn)，這些臨時(shí)點(diǎn)分布于父網(wǎng)格內(nèi)的邊界單元上，那么可采用如下公式判斷點(diǎn)x的位置：

(1)

式中：yi和ni為臨時(shí)點(diǎn)的坐標(biāo)和外法向量。

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

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

(2)

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

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

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

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

(3)

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

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

(4)

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

(5)

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

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

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

二維彈性力學(xué)問題的控制方程為如下幾個(gè)方程：

σij,j+bi=0

(6)

(7)

σij=Cijklεkl

(8)

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

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

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

(9)

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

(10)

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

(11)

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

(12)

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

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

(13)

其中：

(14)

3 邊界元法在重力壩中的應(yīng)用

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

3.1 方法精度驗(yàn)證

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

b1=ex2sinx1+ex2cosx1-

ex1sinx2-ex1cosx2

(15)

b2=ex1sinx2+ex2sinx1-

ex1cosx2-ex2cosx1

(16)

其中位移解析解為：

(17)

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

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

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

圖8 圓孔方板的應(yīng)力計(jì)算結(jié)果Fig.8 Stress calculation results of round hole square plate

3.2 重力壩計(jì)算

該算例考慮如圖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

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

圖10 重力壩位移計(jì)算結(jié)果Fig.10 The displacement calculation results of gravity dam model

圖11 重力壩應(yīng)力計(jì)算結(jié)果Fig.11 Stress calculation results of gravity dam model

4 結(jié) 語

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

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

□

[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.