邊界單元法基礎(chǔ)(直接法)

1、邊界單元法基礎(chǔ)（直接法）一、概述近年來在邊界法方面人們發(fā)表了大量的文章和著作。這些方法是以不同的名稱而提出來的，如“邊界積分方程方法”“邊界積分解”，等等。這種方法的數(shù)值解形式是把所考慮的域的邊界劃分為一系列的單元。邊界單元法簡稱BEM是七十年代興起的一種新的計算方法。它將邊界上的廣義位移和廣義力作為獨(dú)立變量且同時用滿足場方程的奇異函數(shù)（源函數(shù)）作為加權(quán)函數(shù)。所以，它是一種特殊格式的加權(quán)余量法。邊界元法只需將求解域的邊界劃分成單元，故使求解問題的維數(shù)降低，如三維問題可轉(zhuǎn)變成二維問題求解。二維問題可化為一維問題。因而，輸入數(shù)據(jù)大為減少，計算時間縮短。由于它只對邊界離散，故離散誤差僅為來源于邊界，

2、而域內(nèi)變量可由解析式的離散形式直接求得。因此，提高了計算精度。求域內(nèi)變量時，只須改變其數(shù)量和坐標(biāo)位置即可。和有限元法一樣，邊界元法可廣泛地用來解決各種工程問題，如彈性力學(xué)、斷裂力學(xué)、塑性力學(xué)、流體力學(xué)、溫度場和電磁場等。邊界元法分為直接法和間接法。直接法是用物理意義明確的變量來建立積分方程，其中未知函數(shù)就是所求的物理量在邊界上的值；間接法是用物理意義不一定很明確的變量來建立積分方程，如位勢問題中用單層位勢和雙層位勢表示物理量。本部分著重敘述直接法。在用加權(quán)余量法建立積分方程時，所使用的權(quán)函數(shù)是數(shù)學(xué)上的“基本解”。基本解在數(shù)學(xué)上是作為微分方程的特殊的非齊次解定義的，它在每個問題上分別具有不同的物

3、理含義。求這個解，特別是便于解析的形式，一般是不容易的，這是數(shù)學(xué)上的難點(diǎn)。然而，除了特殊問題以外，主要微分方程的基本解，數(shù)學(xué)教科書中有所推導(dǎo)，工程技術(shù)人員可直接引用。邊界元法另一個問題是，代數(shù)方程組的系數(shù)矩陣一般是非對稱的，且非零系數(shù)矩陣為滿秩矩陣，這是由于邊界點(diǎn)與全部邊界單元有關(guān)得出的，編程序時需要注意這一點(diǎn)。我們先介紹位勢問題的邊界單元法公式。這些基本概念對任何其它工程問題是類似的。然后再介紹利用邊界元法求解彈性體受力分析問題。這里強(qiáng)調(diào)的是方法在工程中的應(yīng)用。作為有限元法數(shù)值法的一個補(bǔ)充。二、泊松方程的邊界單元法1積分方程的建立和基本解為了說明邊界單元法的積分方程是如何由加權(quán)余量法推導(dǎo)得來

4、的，我們以泊松方程為例來闡述其全部求解過程，這對了解其它問題的邊界元法求解是有益的?？紤]勢函數(shù)。它在域內(nèi)滿足微分方程，即 （在上） （9-1）在邊界上滿足邊界條件，即 （在上） （在上） （9-2）如圖9-1所示，求解域的總邊界。圖9-1 位勢問題的域和邊界可以證明，對于泊松方程或拉普拉斯方程（），一般加權(quán)余量表示式為 （9-3）式中，W是權(quán)函數(shù)，在邊界法中可令W，稱為相應(yīng)于方程（9-1）的基本解，后面描述它。對式（9-3）左邊第二項進(jìn)行分部積分，即 （9-4）式中右邊第一項的被積函數(shù)形式稱愛因斯坦求和約定，即把式（9-4）右端第一項再次分部積分，得所以把上式代入式（9-3）或去（），考慮到，

5、把左右兩邊界積分合并時，得 （9-5）這就是求解泊松方程的積分方程。如果f = 0，即為求解拉普拉斯方程的積分方程。下面來導(dǎo)出相應(yīng)該方程的基本解。假如有一單位勢作用在物體上的P點(diǎn)，那么基本解就是滿足下列方程2*（p,q）+（p-q）=0 （9-6）的二點(diǎn)函數(shù)，其中p，q為無限域中任意兩點(diǎn)，如圖9-2所示。圖9-2 無限域中的p和q點(diǎn)取p點(diǎn)為坐標(biāo)原點(diǎn)，并把式（9-6）寫成極坐標(biāo)表達(dá)式，即（9-7）式中，是狄拉克函數(shù)，它具有如下性質(zhì)： （9-8）r為p點(diǎn)到q點(diǎn)之間的距離。方程（9-6）或（9-7）的解稱為拉普拉斯方程的基本解。在式（9-7）中，當(dāng)時，則為 （9-9）它的解為 （9-10）式中，C和

6、2為積分常數(shù)，可由r = 0時的條件確定。由于 （9-11）應(yīng)用格林公式 （9-12）由圖9-2（b）所示，以p為圓心，以為半徑作一小圓，將式（9-10）代入式（9-12）并使時取極限，于是得由式（9-11）得 （9-13）由于是勢函數(shù)，C1可取作零，所以 （9-14）式（9-14）即為二維拉普拉斯問題的基本解。同理可求得三維拉普拉斯問題的基本解為 （9-15）現(xiàn)在根據(jù)式（9-5）來建立積分方程。由于于是得 （9-16）式中；考慮時，可寫成式（9-16）就是域內(nèi)任一點(diǎn)p的函數(shù)值與邊界的積分關(guān)系，除了微分方程非齊次項f的積分項外，其余都是邊界積分表達(dá)式，故稱積分方程。圖9-3 用半球包圍的邊界點(diǎn)

7、下面我們來建立邊界解的積分表達(dá)式。如果在式（9-16）中把內(nèi)點(diǎn)p取到邊界上P點(diǎn)（三維舉例），如圖9-3所示。這時，為了避免奇性，需要在P點(diǎn)附近稍加改變，即以P點(diǎn)為中心，以為半徑作一小半球包含P點(diǎn)。這樣P點(diǎn)仍為內(nèi)點(diǎn)。此時把邊界分為兩部分：一部分是鼓起的半球部分，用表示；另一部分即。假如點(diǎn)P取在上（上也一樣分析），對于式（9-16）在邊上的積分分為兩部分，則 （9-17）式（9-17）右端第二工面當(dāng)時取極限，并把三維基本解代入時，則對光滑表面有當(dāng)時，而式（9-17）右端第一項仍為上積分，即式（9-16）右端第一項在上積分同樣分為兩部分，即和。而而 如果P點(diǎn)取在上也可推出同樣的結(jié)論。所以，把P點(diǎn)取在

8、邊界上時，則式（9-16）為 （9-18）上式即泊松方程的邊界解公式。它把邊界上的函數(shù)值與邊界上的導(dǎo)數(shù)值通過基本解的邊界積分之間的關(guān)系聯(lián)系起來了。如果，即為拉普拉斯方程的邊界解公式，對于光滑邊界，即為 （9-19）或者寫成一般式為 （9-20）式中，P，Q皆為邊界上的點(diǎn)，C（P）是與P點(diǎn)處的邊界幾何形狀有關(guān)的常數(shù)，如圖9-4所示。圖9-4 C（P）與邊界幾何關(guān)系 （9-21）式中的為邊界點(diǎn)P處的邊界切線之間的夾角。對光滑邊界，所以。2邊界離散為邊界單元并建立代數(shù)方程組由上面推得了兩個關(guān)系式，（i）即式（9-16），得出域內(nèi)點(diǎn)p的函數(shù)值與邊界上P點(diǎn)的函數(shù)值及導(dǎo)數(shù)值之間的關(guān)系。換句話說，邊界上的函

9、數(shù)值和導(dǎo)數(shù)值求得之后，即可求出域內(nèi)任一點(diǎn)的函數(shù)值。（ii）即式（9-20），得到邊界上任一點(diǎn)的函數(shù)值與邊界上其它點(diǎn)函數(shù)值與導(dǎo)數(shù)值之間的關(guān)系。下面對式（9-20）進(jìn)行邊界離散化，即劃分為邊界單元來求邊界上的函數(shù)值及其導(dǎo)數(shù)值。以二維問題為例，離散方案一般有三種情況，如圖9-5所示，即圖9-5（a）常值單元，（b）線性單元，（c）二次單元。所謂常值單元即在單元中函數(shù)值和導(dǎo)數(shù)值均為常數(shù)，用中個單元具有兩個節(jié)點(diǎn)，用節(jié)點(diǎn)值對對函數(shù)作線性插值。二次邊界單元，每個單元有三個節(jié)點(diǎn)，函數(shù)和導(dǎo)數(shù)在單元內(nèi)作二次變化，單元的幾何形狀適應(yīng)于曲線邊界。下面用常值單元將邊界解公式化成節(jié)點(diǎn)值的線性代數(shù)方程組。如圖9-6所示，把

10、邊界劃分為n個邊界單元。其中n1個單元屬于；n2個單元屬于。單元總數(shù)為n = n1 + n2。對于節(jié)點(diǎn)Pi，式（9-20）可寫為 圖9-5 邊界元的不同類型 圖9-6 常值單元（a）常值單元；（b）線性單元；（c）二次單元 （9-22）由于和q在單元內(nèi)為常數(shù)故可提到積分號外面。Pi為單元上的節(jié)點(diǎn)，Q點(diǎn)為單元上的任意點(diǎn)。上的積分與點(diǎn)Pi和單元有關(guān)。 （9-23）由于已知（基本解），也可求得，故式（9-23）亦可求得。于是式（9-22）可寫成 （9-24）進(jìn)一步可寫成 （9-25）式中 （9-26）把式（9-25）寫成矩陣形式，則為式中 （9-27）； （9-28）； （9-29）在式（9-27）

11、中，兩邊都有未知量和已知量，即在，（邊界上函數(shù)值給定）；在上，（邊界上導(dǎo)數(shù)值給定）。給在上有個q值待求，在上有個值待求，把式（9-27）中的已知值和未知值重新整理分別寫出，則為已知已知未知 未知 （9-30）把式（9-30）簡寫為 （9-31）把式（9-31）的未知量移到方程左邊，已知量移到方程的右邊，則得 （9-32）寫成常用的方程組表達(dá)式時，即為 （9-33）式中在式（9-33）中X 中包含有個待求的未知函數(shù)值和個待求的未知導(dǎo)數(shù)值，其未知量仍為n個。F中有個函數(shù)值和個導(dǎo)數(shù)值作為已知的邊界條件給出，n個方程正好求解n個未知數(shù)。矩陣A是一個非稀疏的非對稱的矩陣，與有限元系數(shù)矩陣K是不同的。我們

12、把式（9-16）（這里）改寫為離散形式，即 （9-34）一旦整個邊界上的值和q值都求得了的話，就可利用上式求得域內(nèi)任何點(diǎn)p處的值。3系數(shù)矩陣H，G中元素的計算首先來計算對角線元素和。是的系數(shù)，它是在單元上積分得到的。當(dāng)時，是邊界單元的中點(diǎn)（常值單元），而Q是邊界單元上的任意點(diǎn)，如圖9-7所示。邊界單元的兩端點(diǎn)的坐標(biāo)分別是（x1，y1）和（x2，y2）表示。圖9-7 在上的積分路線由式（9-26）及（9-23）第一式知，對光滑邊界，有而 （9-35）式中，r為Pi到Q的距離，而n是外表面法線，r和n正交，故，所以于是 是的系數(shù)，它也是在上積分的。由式（9-23）知：把拉普拉斯方程基本解式（9-1

13、4）代入，則 （9-36）取元因次局部坐標(biāo)，在單元中點(diǎn)，在兩端點(diǎn)分別±1。令，為單元邊界的長度。代入到式（9-36）中，則 （9-37）再來考慮非對角線中和Gij。如圖9-8所示。當(dāng)時，Pi是單元的中點(diǎn)，Q是單元上的任意點(diǎn)，和Gij的計算，可采用數(shù)值積分（如高斯四點(diǎn)積分）進(jìn)行。我們先導(dǎo)出數(shù)值積分所采用的公式。由式（9-23）第一式，得圖9-8 Pi在上的積分路線由于故 （9-38）圖9-9 Pi到的垂直距離把式（9-38）代入式（9-23），得 （9-39）式中，sj為邊界單元的長度；rij為Pi到中Q點(diǎn)的距離；hij為Pi點(diǎn)到單元的垂直距離。由圖9-9所示，有如果 r和n

14、的夾角為銳角則hij取正號，為鈍角則取負(fù)號，如圖9-10所示。也可由Pi到兩端點(diǎn)的矢量r1和r2的叉積來判別，即逆時針時，hij > 0，順時針時hij < 0，如圖9-11所示。 圖9-10 hij的正負(fù)號 圖9-11 由確定hij的正負(fù)號如果把r1和r2分別寫成 ，則式中，于是hij的正負(fù)號可按c值來決定，即c>0，hij > 0；c<0 hij<0。Gij的計算可由式（9-23）第二式得，即 （9-41）由式（9-39）和（9-41）所表達(dá)Hij和Gij可用高斯四點(diǎn)數(shù)值積分算出，即 （9-42）式中Wk為權(quán)系數(shù)，如圖9-12所表示的積分方案。圖9-12

15、 四點(diǎn)數(shù)值積分的幾何定義在泊松方程中，有一項域積分，如式（9-18）左邊第一項。在離散解時，也要把這項并入到式（9-33）右端項F中去。其作法可把求解域劃分為單元，然后采用數(shù)值積分，即式中，Ae是域單元面積；k從1到7表示采用7點(diǎn)高斯積分方案（見第三章表3-1）；m是域劃分的單元數(shù)，見圖9-13所示。當(dāng)然，這項積分也可化為邊界積分，在彈性全求解時敘述。4域內(nèi)p點(diǎn)函數(shù)值的計算在計算出邊界上的值和q值以后 ，利用式（9-34）式可計算域內(nèi)任一點(diǎn)p的函數(shù)值。如果所求為泊松方程，則需加入域積分項，這時式（9-34）為式中，Hij和Gij雖然也由式（9-39）和（9-41）計算，但其中rij要換成域內(nèi)p

16、點(diǎn)到邊界上的Q點(diǎn)的距離。5線性邊界單元和高階單元的計算首先我們分析線性邊界單元。線性單元具有兩個節(jié)點(diǎn)，函數(shù)和導(dǎo)數(shù)在單元內(nèi)作線性變化。用單元的相交點(diǎn)取為節(jié)點(diǎn)如圖9-14所示。 圖9-13 域積分 圖9-14 線性單元對n個單元，把式（9-20）寫成離散形式為 （9-43）由于和q在單元內(nèi)呈線性變化，不能直接提到積分號外面，故計算要麻煩一些。利用兩點(diǎn)插值公式，對和q進(jìn)行插值，寫成下列形式： （9-44）式中，和分別為單元節(jié)點(diǎn)j和節(jié)點(diǎn)j + 1的函數(shù)值及導(dǎo)數(shù)值。是無因次單元局部坐標(biāo)，即如圖9-15所示。N1和N2當(dāng)?shù)暮瘮?shù)，分別由下式給出，即圖9-15 單元上的坐標(biāo)下面分別來求式（9-43）的積分。把

17、式（9-44）分別代入式（9-43）中去，左邊第二項積分為 （9-45）令； （9-46）于是 （9-47）同理 （9-48）式中 ； （9-49）把式（9-47）和（9-48）代入式（9-43）中去，得 （9-50）上式展開后得 （9-51）對于單連域問題，邊界首尾卸接，這時； （9-52）將式（9-51）整理后得 （9-53）令 ； （9-54）式中 ； 于是，對于節(jié)點(diǎn)Pi，式（9-43）可寫成 （9-55）與常值單元一樣，進(jìn)一步可寫成移項后，得AX = F （9-56）下面導(dǎo)出微元的計算，式（9-46）和（9-49）N1和N2是的函數(shù)而積分微元是對進(jìn)行的，故必須把變換到局部坐標(biāo)系中去，如

18、圖9-16所示。與總體坐標(biāo)（x，y）的關(guān)系為圖9-16 曲線邊界的幾何定義由于所以 （9-57）稱為雅可比行列式。如果對總體坐標(biāo)x，y也進(jìn)行線性插值，即 （9-58）式中，（xj，yj）和（xj+1yj+1）分別為邊界單元兩端點(diǎn)的坐標(biāo)值。由式（9-58）可求得代入式（9-57），得 （9-59）式中sj為單元的長度。于是，對于線性單元，式（9-46）和（9-49）分別由下式表達(dá)， 即 （9-60）式（9-60）中的積分同樣可用四點(diǎn)高斯積分計算。圖9-17 二次單元二次單元是在邊界單元內(nèi)的函數(shù)和其導(dǎo)數(shù)均采用二次插值，如圖9-17所示，其式為 （9-61）式中，為的二次函數(shù)，即 （9-62）與線性

19、單元類似的作法，求下面的積分： （9-63）式中 （9-64）設(shè)，分別為單元1，2，3節(jié)點(diǎn)的總體坐標(biāo)。與線性單元一樣，可寫成 （9-65）于是有 （9-66）將式（9-66）代入式（9-57），得 （9-67）其它推導(dǎo)與線性單元相同，但在每個單元的首尾節(jié)點(diǎn)（如節(jié)點(diǎn)1和2），同樣有兩項系數(shù)疊加以組成總體方程組的系數(shù)。三、彈性體受力分析的邊界元法1基本關(guān)系和預(yù)備知識為了簡化起見，在彈性力學(xué)的表達(dá)式中用張量表示。所以，首先把用到的張量符號簡介如下。指標(biāo)記法：在直角坐標(biāo)系中，通常用x，y，z表示坐標(biāo)。 如果分別用x1，x2，x3來表示上述三個坐標(biāo)，記為xi（i =1，2，3）就稱指標(biāo)記法。這里i稱為下

20、標(biāo)。ai（i =1，2，N）代表a1，a2，aN N個量，則代表N×M個量，則代表N×M×L個量，依此類推。求和約定及啞標(biāo)：某指標(biāo)在某一項中重復(fù)出現(xiàn)，且僅重復(fù)一次，則該項代表一個和式，按重復(fù)指標(biāo)的取值范圍求和，如表示：表示求和重復(fù)的指標(biāo)為啞標(biāo)，用什么字母均可。不求和的指標(biāo)稱為自由標(biāo)。如上例i就是自由標(biāo)?？肆_內(nèi)克符號： 定義為=1（i = j時）；=0（i j時）當(dāng)i, j =1，2，3時，表示可知； ； ； （當(dāng)i, j, k為1，2，3循環(huán)序列）（當(dāng)i, j, k為1，2，3逆循環(huán)序列）（當(dāng)i, j, k有相同指標(biāo)）置換符號eijk定義如下：eijk共代表27個

21、量，有21個量為零。顯然：導(dǎo)數(shù)記號：求導(dǎo)符號用“，”表示。中用指標(biāo)記法表示彈性體的基本關(guān)系時，其平衡方程記為 （在內(nèi)） （9-68）力的邊界條件為 （在上） （9-69）位移邊界條件為 （在上） （9-70）幾何方程為 （9-71）各向幾性彈性體的物理方程為 （9-72）式中 ； （9-73）稱為拉密常數(shù)；G稱剪切模量。或者更簡明地寫為式中 用位移表示的平衡方程為 （9-74）如果考慮熱影響時，則用位移表示的熱彈性方程為 （9-75）考慮熱影響時的應(yīng)力表達(dá)式為 （9-76）式中 式中，a是熱膨脹系數(shù)，T是溫差。以上指標(biāo)記法，對于三維問題，3；對于二維問題，i，j，k=1，2。對于平面應(yīng)力問題需

22、將彈性常數(shù)作相應(yīng)改變，即使E，a相應(yīng)變?yōu)?，?彈性問題邊界積分方程表達(dá)式如前所描述的勢問題一樣，我們把平衡方程和兩類邊界條件寫成加權(quán)余量表達(dá)式為 （9-77）式中，和為權(quán)函數(shù)，分別為加權(quán)場的位移和面力，也即基本解。同樣有 （9-78）現(xiàn)在把式（9-71）和（9-72）應(yīng)用到近似解和中權(quán)場中去。對式（9-77）左邊第一項進(jìn)行分部積分，便得 （9-79）如果考慮溫度影響時，則需把式（9-75）代入到式（9-79）中左邊第一項中去，即 （9-80）有關(guān)溫度影響，可并入體積力的域積分，在后面推導(dǎo)體積力時詳細(xì)論述，在此省略。對式（9-79）再次進(jìn)行分部積分，且用互換定理 （9-81）時，便得到 （9-8

23、2）方程左邊體積力的域積分不含未知函數(shù)，有關(guān)它的計算在以后介紹?，F(xiàn)在的問題是如何把左邊第一項轉(zhuǎn)變?yōu)檫吔绶e分。彈性問題的基本解是滿足平衡方程 （9-83）的解。這里 函數(shù)是表示在i點(diǎn)在l方向施加的單位載荷。如果對上述作如下運(yùn)算，即于是得 把上式代入式（9-82），得 （9-84）式中，表示域內(nèi)任意點(diǎn)i（或p）沿l方向的位移；為基本解，即在無限域中式（9-83）的解。其意義分別為在i點(diǎn)沿l方向有單位力作用時所產(chǎn)生的位移和表面力。如果考慮i點(diǎn)的三個方向時，則式（9-84）為 （9-85）或者把i換p表示時，也可將上式寫成 （9-86）式中，和分別在域內(nèi)p點(diǎn)沿l方向作用單位力時，而在邊界上的Q點(diǎn)在k方

24、向產(chǎn)生的位移分量和面力分量。l，k分別代表p點(diǎn)（或i點(diǎn)）和Q點(diǎn)的有關(guān)坐標(biāo)方向。對三維問題，l，k = 1，2，3。l，k的符號意義如圖9-18所示。對二維問題，l，k = 1，2。式（9-86）是彈性體的邊界積分方程。該式表明區(qū)域內(nèi)任意p點(diǎn)的位移分量ul（p）和邊界上Q點(diǎn)的位移分量u（Q）和表面力分量pk（Q）之間的關(guān)系。如果邊界上所有未知的位移分量和表面力分量已求得，那么域內(nèi)任意點(diǎn)的位移分量可以通過該式求得。圖9-18 在p點(diǎn)在l方向（l = 1）作用單位力而在Q點(diǎn)所產(chǎn)生的3二維彈性體的基本解和邊界元方程的表達(dá)式二維彈性平面問題的基本解，可由無限域的式（9-83）導(dǎo)出。對于平面應(yīng)變問題，位移

25、和面力的基本解分別為 （9-87）對于平面應(yīng)力問題，則分別為 （9-88）式中，r是單位載荷作用點(diǎn)p到所考慮的點(diǎn)Q之間的距離。nj是該點(diǎn)的表面法線的方向余弦。n為外表面法線。當(dāng)l = k時，時，。下面來推導(dǎo)邊界點(diǎn)的積分方程。式（9-86）對區(qū)域內(nèi)任意點(diǎn)都適用。當(dāng)p點(diǎn)取在邊界上P時，基本解會產(chǎn)生奇性。與處理勢問題一樣，以奇點(diǎn)P為圓心以為半徑作小半圓，如圖9-19所示。于是邊界又分為兩部分；一部分是；另一部分是。把基本解代入式（9-86），當(dāng)取極限時，式（9-86）右邊第一項在上積分為 零，而在上積分仍為上的積分。式（9-86）左邊第二項在上積分，當(dāng) 它與邊界形狀a角有關(guān)。積分后代入到式（9-86

26、）中去，得 （9-89）圖9-19 p點(diǎn)取到邊界P時的處理式中，P即所取的邊界點(diǎn)；而為邊界上的位稱移；Clk是與邊界幾何形狀有關(guān)的常數(shù)，表達(dá)式為 （9-90）式中，對于光滑邊界，所以 （9-91）式（9-89）是邊界上的位移分量uk（P），uk（Q）和表面力分量pk（Q）之間的邊界積分關(guān)系式。由該式可以求出邊界上的全部未知位移分量和表面力分量。與勢問題一樣，我們把式（9-89）的邊界積分進(jìn)行離散化，可離散成常值單元、線性單元及二次和高次單元。該式寫成矩陣表達(dá)式為 （9-92）是的2×2矩陣，P*是的2×2矩陣，即； （9-93）； ； （9-94）下面我們把式（9-92）離

27、散成n個單元，并采用常值單元，則該式為 （9-95）式中，Ui為邊界單元的位移矢量，Uj，Pj分別為邊界單元的位移矢量和面力矢量。由于采用的是常值單元，和勢問題一樣節(jié)點(diǎn)取在單元中點(diǎn)，位移變量和力變量在單元內(nèi)是常值，故可提到積分號外面去。n是邊界單元數(shù)。右端第二項為體積力的域積分。令 （9-96）式中，m為域積分單元數(shù)，l為域積分的積分點(diǎn)數(shù)，Wk為權(quán)系數(shù)，Ae為單元e的面積。這是把體積力的域積分用加權(quán)數(shù)值積分得出的。當(dāng)然，也可轉(zhuǎn)變?yōu)檫吔绶e分計算，在后面專門介紹。于是式（9-95）可寫成 （9-97）進(jìn)一步可寫成 （9-98）式中，當(dāng)時，；當(dāng)i = j時，對光滑表面，方程（9-97）或（9-98）

28、對節(jié)點(diǎn)i給出一組方程式（三維是三個，二維是兩個），Uj和Pj都是在節(jié)點(diǎn)j上的未知數(shù)。和Gij是節(jié)點(diǎn)i與物體表面上全部節(jié)點(diǎn)相關(guān)的系數(shù)。這里與勢問題的區(qū)別在于每個節(jié)點(diǎn)位移和面力有兩個分量（三維問題是三個），故和Gij是2×2子矩陣。對所考察的每個節(jié)點(diǎn)都寫出式（9-97），于是有 （9-99）式（9-99）寫成簡化形式為 （9-100）求解方程組（9-100）以前，必須加入邊界條件，此處的邊界條件為在上 （）在上 （）式（9-100）共有2n個位移值和2n個面力值，但有n1個位移值和n2個面力值是作為邊界條件給定的，因而式（9-100）中尚有2n個未知量，即有n1個面力值和n2個位移值待求

29、，把已知位移和已知面力代入式（9-100）中去，使方程組重新組合，即將所有未知項表示矢量X，寫在等式的左邊，其結(jié)果為 （9-101）體力矢量項也包括在X中。4系數(shù)矩陣中元素Hij和Gij的計算首先計算對角線元素Hij和Gii。當(dāng)j = i時，式（9-96）中有關(guān)邊界積分是在的中點(diǎn)，而Q點(diǎn)取在單元中的任意點(diǎn)（在同一個單元中），如圖9-20所示。圖9-20 Pi和Q在中積分由式（9-93）及（9-96）可知 （9-102）把式（9-88）代入上式，考慮到，同時由圖9-20可見，在中，。在中，。代入計算后，便得到。于是得而ci也可由研究剛體運(yùn)動求得。若設(shè)在任一方向上剛體有一個單位的位移，則式（9-1

30、00）變?yōu)槭街?，Il是在l方向上剛體位移的單位矢量，于是H的 對角線項可簡寫為這意味著并不需要用顯式來確定系數(shù)，對光滑表面。Gii的計算是由式（9-93）和（9-96）中第二式進(jìn)行的，即 （9-103）把式（9-88）中的代入上式積分，注意，當(dāng)r = 0時為奇點(diǎn)，故積分限可取為到，使時取極限，計算后得 （9-104）下面介紹非對角線元素Hij和Gij的計算。當(dāng)時，Pi是取在邊界單元上的中點(diǎn)，而是在其它邊界單元上的任意點(diǎn)，如圖9-8所示。計算式仍由式（9-102）和（9-103）計算，只是把下標(biāo)ii換成ij，也是2×2矩陣。把式（9-98）中的和分別代入和Gij中，計算中令，（見圖9-

31、8），于是具體表達(dá)式分別為 （9-105）和 （9-106）假如采用4點(diǎn)高斯積分計算上面積分，只需將微元代入上式，如其它，類似計算。sj為單元邊界長度；Wk為權(quán)系數(shù)；h為Pi點(diǎn)到單元垂直距離；r為Pi到單元積分點(diǎn)的距離；n1，n2分別為的方向余弦。5域內(nèi)點(diǎn)的位移和內(nèi)點(diǎn)、邊界點(diǎn)應(yīng)力的計算一旦邊界上的位移和面力算出后，就可以用式（9-86）求得域內(nèi)點(diǎn)的位移。把該式寫成離散形式時，則為 （9-107）式中，Ui為內(nèi)點(diǎn)Pi的位移矢量。Hij和Gij分別由式（9-105）和（9-106）計算，但式中r為內(nèi)點(diǎn)Pi到 邊界點(diǎn)之間的距離。下面引出內(nèi)點(diǎn)應(yīng)力和邊界點(diǎn)應(yīng)力的計算公式。由式（9-72）可知，應(yīng)力表達(dá)式

32、為 （9-108）將式（9-106）代入幾何方程（9-71），再代入物理方程（9-72），得內(nèi)點(diǎn)應(yīng)力表達(dá)式為 （9-109）上式是全部在內(nèi)點(diǎn)進(jìn)行的。令 （9-110）于是應(yīng)力表達(dá)式為 （9-111）3階張量Dkij和Skij的推導(dǎo)簡要說明如下。首先將在p點(diǎn)對xj微分： （9-112）求得 ； 代入式（9-112），于是得 （9-113）式中，如果將下標(biāo)j換為l，且，則 （9-114）將式（9-113）的下標(biāo)l換為i，便得 （9-115） 如將式（9-113）的下標(biāo)l和j分別換為j和i，就可以得到 （9-116）將式（9-114），（9-115），（9-116）代入式（9-110）整理合并后，便

33、得 （9-117）現(xiàn)在推導(dǎo)Skij表達(dá)式。在p點(diǎn)將對xj微分得 （9-118）將上面關(guān)系以及前面導(dǎo)出，等代入式（9-118），得 （9-119）將上式中下標(biāo)j換為l，便得到；將下標(biāo)l換成i，便得到；將下標(biāo)l和j分別換成j和i，便得到。將其結(jié)果均代入式（9-10）第二式中去，得 （9-120）當(dāng)采用常值單元并把式（9-111）寫成離散后的矩陣形工，得任意內(nèi)點(diǎn)p的應(yīng)力為 （9-121）式中 ； ；中的系數(shù)為，；中的系數(shù)為，；中的系數(shù)為，。它們分別表示如下： （9-122）和 （9-123）式中r為內(nèi)點(diǎn)p到邊界點(diǎn)Q之間的距離。上述公式在編程序計算時，可用4點(diǎn)高斯積分解決，例如：計算和可把式（9-12

34、2）和（9-123）的第一式分別寫成：和其它類推，只是把積分寫成和式，把而已。在數(shù)值積分式中Wk為權(quán)系數(shù)。邊界上的應(yīng)力在工程結(jié)構(gòu)中往往是感興趣的。邊界上的全部位移和面力求出后，就可以計算邊界上的應(yīng)力。由于邊界上的位移uk（Q）沿邊界s上的變化是可以求得的（近似求出可）。再加上物理方程（平面問題為三個） 和力的邊界條件（二個），七個方程正好求出七個未知數(shù)，即三個應(yīng)力和四個位移在邊界上的導(dǎo)數(shù)。具體寫為 （9-124）式中，l = 1，2。sl是邊界的參數(shù)，在上是已知的（差分計算或?qū)k直接微分）。把上式寫成矩陣形式為 （9-125）式中，；是垂直s的單位矢量；s是弧長。四、體積力的邊界積分計算 體

35、積力的計算，在機(jī)械工程中是經(jīng)常遇到的。在熱動力機(jī)械中，體積力包含有重力、離心力和溫度載荷。含有體積力的彈性體的邊界積分方程是式（9-86）或（9-89）。由該二式看出，有一項是域積分。對這項的處理，如果采用在域內(nèi)劃分單元，然后采用數(shù)值積分如式（96-96）那樣，則是不勝其煩的，且所需的數(shù)據(jù)準(zhǔn)備也較多。如果體積力為常量或無熱源的穩(wěn)態(tài)熱載荷，就可以把體積力的域積分轉(zhuǎn)化為邊界積分。這就充分利用了邊界單元法的優(yōu)點(diǎn)?，F(xiàn)在我們先針對二維問題來討論這種算法。域積分的積分項在積分方程中為 （9-126）式中，是整個彈性體的體積域。為了把它轉(zhuǎn)變成邊界積分，我們在這里定義一個張量Gki（P，Q）p是力點(diǎn)，Q是場點(diǎn)

36、。這個張量叫伽遼金張量。二維問題的表達(dá)式為 （9-127）由它組成的基本解為 （9-128）把式（9-127）代入到式（9-128）中去，得 （9-129）可以看出式（9-129）與（9-88）只差一個常數(shù)。常數(shù)只表示剛體位移，對計算結(jié)果是無影響的。把式（9-128）代入到式（9-126）中去，得 （9-130）這個公式似乎復(fù)雜些，但易于轉(zhuǎn)變?yōu)檫吔绶e分。下面我們針對幾種常見的體積力把式（9-130）轉(zhuǎn)化為邊界積分公式。1重力載荷在常重力場中，由于為常數(shù)，為質(zhì)量密度，為加速度。根據(jù)高斯定理，有 （9-131）這就是重力為體積力的體力項轉(zhuǎn)變?yōu)檫吔绶e分的表達(dá)式。把式（9-127）微分后代入式（9-1

37、31）中去，便得 （9-132）于是，內(nèi)點(diǎn)邊界積分方程可以寫為 （9-133）式中 （9-134）在重力場作用下的應(yīng)力表達(dá)式為 （9-135）式（9-135）中的是這樣求得的，首先把式（9-134）代入到應(yīng)力位移關(guān)系（即式（9-72）中去，得到 （9-136）再把式（9-134）求導(dǎo)代入式（9-136）得到 （9-137）2離心力載荷在二維平面情況下，旋轉(zhuǎn)體單位體積所產(chǎn)生的離心力為 （9-138）令 =常數(shù)把式（9-138）代入到式（9-130）中去，得 （9-139）現(xiàn)在對式（9-159）作如下變換：二式相加得到等代入式（9-139），便得到 （9-140）用高斯定理把域積分轉(zhuǎn)變?yōu)檫吔绶e分，

38、即 （9-141）由于是對稱的，方程（9-141）稍加改變，便得 （9-142）如重力載荷情形一樣，民可把式（9-142）簡寫為 （9-143）式中 （9-144）應(yīng)力表達(dá)如同式（9-135），對于離心力所產(chǎn)生的則為 （9-145）3熱載荷對二維穩(wěn)態(tài)熱載荷，體積域積分是由溫差所產(chǎn)生的。根據(jù)熱彈性方程的馬克斯威爾-貝蒂互換定理，得內(nèi)點(diǎn)邊界積分方程為 （9-146）于是 （9-147）式中，由式（9-128）得 （9-148）由改變整數(shù)虛標(biāo)，于是式（9-148）變?yōu)?（9-149）代入式（9-147），得 （9-150）對穩(wěn)態(tài)熱載荷，有 （9-151）于是，我們可以把式（9-150）寫成 （9-1

39、52）由高斯定理，可以把上式轉(zhuǎn)變?yōu)檫吔绶e分，即 （9-153）這就是由于溫度載荷所產(chǎn)生的邊界積分式。該式也可簡寫為 （9-154）式中 （9-155）把式（9-127）代入式（9-155），分別得 （9-156）于是，對于熱彈性方程位移邊界積分表達(dá)式為 （9-157）把式（9-157）微分后代入應(yīng)力位移關(guān)系中去，便得 （9-158）式中 （9-159）五、三維彈性體的邊界元公式1三維彈性體的基本解和三階張量彈性力學(xué)三維問題的基本方程為式（9-68）（9-76）諸式，只是其中下標(biāo)的取值范圍分別為。積分方程表達(dá)式的推導(dǎo)與二維問題是類似的，在此不作重復(fù)。域內(nèi)任意點(diǎn)p的位移ul（p）與邊界積分關(guān)系以及

40、邊界點(diǎn)位移cul（P）與邊界積分關(guān)系仍為式（9-86）和（9-89），只是其中l(wèi)，k =1，2，3。式中的基本解，則分別為 （9-160）（9-161）式中，見圖9-21所示。圖9-21 幾何意義對于光滑表面，系數(shù)c仍為，非光滑表面可由二維推廣。其它如邊界方程離散化，方程組的建立，與二維問題大體是類似的。只是現(xiàn)在自由度增加了，每個邊界節(jié)點(diǎn)有三個位移分量和三個面力分量，n個節(jié)點(diǎn)有3×2n個方程式。積分邊界是三維空間域表面，即面積分?？刹捎糜邢拊ㄓ嘘P(guān)插值公式和數(shù)值積分方案進(jìn)行。在式（9-111）應(yīng)力表達(dá)式中的三階張量Dkij和Skij分別為 （9-162）式中，導(dǎo)數(shù)取在邊界上，。2三維

41、體積力的邊界積分式三維體積力邊界積分公式的導(dǎo)出，與二維是類似的。相應(yīng)于三維基本解的伽遼金張量為 （9-163）將上式微分兩次代入式（9-128），得 （9-164）重力載荷的三維體積力邊界積分式：把式（9-163）微分兩次代入式（9-130）中，得 （9-165）式中 （9-166）用與二維問題類似的方法得應(yīng)力表達(dá)式中的為 （9-167）離心力載荷三維體積力邊界積分式：與二維問題一樣，離心力體積力邊界積分也可寫成 （9-168）式中 （9-169）應(yīng)力表達(dá)式中的為 （9-170）三維穩(wěn)態(tài)熱載荷的邊界積分式：在穩(wěn)態(tài)熱載荷中，體積力的邊界積分項也可寫為式（9-154）的表達(dá)式，即式中 （9-171

42、）應(yīng)力項中的和分別為 （9-172）六、軸對稱體的邊界元公式燃?xì)鉁u輪發(fā)動機(jī)等機(jī)械含有大量的軸對稱體零、部件。這些零、部件用邊界元法分析其結(jié)構(gòu)完整性，會使其輸入數(shù)據(jù)和計算時間大為減少。1軸對稱體基本方程及伽遼金矢量用圓柱坐標(biāo)表示的軸對體的納維葉方程為 （9-173）式中其矢量形式為 （9-174）為了求解方程組（9-173），定義伽遼金矢量為 （9-175）把式（9-175）代入式（9-173）中，得雙調(diào)和方程為 （9-176）式（9-176）寫成標(biāo)量時，則為 （9-177）式中，Gr和Gz是伽遼金矢量的徑向和軸向分量；Fr和Fz是體積力分量。這里的體積力項可通過應(yīng)用Dirac 函數(shù)來表示，即

43、（9-178）式中，R和Z是載荷點(diǎn)的坐標(biāo)；r和z是場點(diǎn)的坐標(biāo)。根據(jù)Dirac 函數(shù)的性質(zhì)，有 （當(dāng)時）為了得到邊界解的公式，需要找到軸對稱體納維葉方程的基本解。由前述，對一般三維問題，其基本解是均勻無限體的點(diǎn)載荷的開爾文問題的解。對于軸對稱問題，基本解的物理意義是軸對稱環(huán)載荷的開爾文解。今考慮一個具有任意橫截面形狀的軸對稱體，在體內(nèi)p點(diǎn)為載荷點(diǎn)，其坐標(biāo)是。邊界點(diǎn)為Q點(diǎn)，其坐標(biāo)為），如圖9-22所示。小寫p表示固定的載荷點(diǎn)，大寫Q表示可改變位置的邊界點(diǎn)。其基本解是描述Q點(diǎn)的狀態(tài)與p點(diǎn)的環(huán)載荷之間的關(guān)系。圖9-22 軸對稱體有兩種方法可以求得軸對稱體的基本解。第一種方法是圍繞對稱軸沿環(huán)載荷路徑積分

44、在維的基本解；第二種方法是直接解出軸對稱形式。這里介紹的是第二種方法。用積分變換方法，可以證明，下列解 （9-179）是滿足無體力情況下方程（9-177）的。式中稱為第二類 階v次的勒讓德函數(shù)。根據(jù)遞推關(guān)系可以用零階的勒讓德函數(shù)來表示-1階的勒讓德函數(shù)，即 （9-180）如果把式（9-180）代入式（9-179）中去，將其結(jié)果代入式（9-175）便得到用零階勒讓德函數(shù)表示的軸對稱體位移核函數(shù)。在具體計算時，可用第一類和第二類橢圓積分來表示勒讓德函數(shù)及其導(dǎo)數(shù)。2軸對稱彈性問題的基本解（位移核和力核）為了方便起見，用完全橢圓積分來表示伽遼金矢量而不用勒讓德函數(shù)時，便得到 （9-181）式中，和分別是第一類和第二類完全橢圓積分并定義為 （9-182）式中把伽遼金矢量分量代入式（9-175）中，得位移矢量為 （9-183）式中，er和ez分別是徑向和軸向的單位矢量。在邊界Q點(diǎn)，徑向和軸向的位移分量為 （9-184）這里pr和pz分別是內(nèi)