苗 雨， 孫恬粲， 鄭俊杰

（華中科技大學(xué)土木工程與力學(xué)學(xué)院，湖北武漢 430074）

新型快速多極雜交邊界點法在復(fù)合材料熱傳導(dǎo)中的應(yīng)用

苗 雨， 孫恬粲， 鄭俊杰

（華中科技大學(xué)土木工程與力學(xué)學(xué)院，湖北武漢 430074）

將雜交邊界點法應(yīng)用于復(fù)合材料的熱傳導(dǎo)模擬，推導(dǎo)一種求解復(fù)合材料的方程，該方程減少計算自由度，效率更高.將新型快速多極算法與雜交邊界點法結(jié)合進行大規(guī)模計算，數(shù)值算例中對包含大量粒子的復(fù)合材料進行模擬，結(jié)果表明快速多極雜交邊界點法可行，具有一定的應(yīng)用前景.

雜交邊界點法；三維復(fù)合材料；方程；快速多極算法

0 引言

在工程計算中［1－3］，主要存在兩大難題，①是網(wǎng)格生成，②是大規(guī)模計算.傳統(tǒng)的有限元法等基于網(wǎng)格的算法需要在劃分網(wǎng)格上耗費大量時間，對于復(fù)雜的問題，如復(fù)合材料等問題尤為嚴(yán)重.對于網(wǎng)格劃分這一難題，本文采用近年來發(fā)展起來的一種無網(wǎng)格法——雜交邊界點法［4－7］來解決，該方法基于移動最小二乘與修正變分原理，不僅可以將求解問題的維數(shù)降低一維，而且無需劃分網(wǎng)格.該算法已經(jīng)應(yīng)用于彈性力學(xué)和彈性動力學(xué)等多個領(lǐng)域，具有很大的發(fā)展?jié)摿?

然而，傳統(tǒng)的雜交邊界點法的系數(shù)矩陣為密集非對稱矩陣，無法在普通微機上求解大規(guī)模問題.為了解決大規(guī)模計算這一難題，本文將新型快速多極算法（new fast multipole method）［8－11］與雜交邊界點法結(jié)合，推導(dǎo)了相應(yīng)的公式.并成功應(yīng)用于復(fù)合材料熱傳導(dǎo)問題的大規(guī)模模擬中，數(shù)值算例表明，該新型快速多極雜交邊界點法是有效的.

首先推導(dǎo)雜交邊界點法求解復(fù)合材料熱傳導(dǎo)問題的公式，然后推導(dǎo)新型快速多極雜交邊界點法求解三維熱傳導(dǎo)問題的公式，最后給出數(shù)值算例.與張見明等人［12］的方法相比，本文求解復(fù)合材料的熱傳導(dǎo)問題的公式可以減少計算自由度，且在初始快速多極算法的基礎(chǔ)上，新型快速多極算法應(yīng)用于雜交邊界點法，可進一步減少計算時間.

1 復(fù)合材料中的雜交邊界點法

在熱傳導(dǎo)問題中，雜交邊界點法的三個獨立變量分別為：域內(nèi)的溫度φ，邊界上的溫度φ～和邊界上的法向溫度流q～．

邊界上φ～和q～采用移動最小二乘近似

其中，ΦJ（s）為節(jié)點sJ的移動最小二乘形函數(shù).

域內(nèi)的φ和q采用基本解近似

其中κ為熱傳導(dǎo)率，r（Q，sJ）為兩點間的距離.

考慮一個包含n個子域的問題，對于每一個子域通過修正變分原理可得

考慮如圖1中所示的復(fù)合材料模型，對于基體材料有

其中，0表示基體子域，其他表示粒子子域.

對于粒子k有

在基體與粒子的交界面上存在連續(xù)性條件，于是有

通過式（12）－（17）我們最終得到

其中，A和d通過組裝邊界條件得到.如果包含的粒子形狀和材料相同，則Ck只需要計算一次.

1.1 快速多極雜交邊界點法

快速多極雜交邊界點法主要用來計算矩陣與向量相乘，對于本文的問題，矩陣與向量的相乘為如下兩式之一

快速多極算法通過將遠(yuǎn)端節(jié)點的相互作用化為格子與格子之間的相互作用來減少計算量.考慮兩個格子Ca和Cb分別包含Na和Nb個粒子，見圖2.方便起見，我們首先只計算式（21）.

1.1.1 初始快速多極雜交邊界點法

首先將基本解進行如下展開

上式在｜O1Q｜＞｜O1sJ｜時成立，其中Rn，m和Sn，m見參考文獻［13］.

將式（23）代入式（21），計算Cb的貢獻

其中Mn，m（O1）為多極矩，且

假設(shè)Ca和Cb屬于兩個大一點的格子和并且他們依舊相距很遠(yuǎn)，見圖2，則有多極矩到多極矩間的傳遞

于是式（24）可寫為

以上便是初始快速多極雜交邊界點法的基本傳遞公式，式（22）的計算方法類似.

1.1.2 新型快速多極雜交邊界點法

以上初始算法中多極矩到局部系數(shù)的傳遞的計算量為O（189p4），是初始算法的計算瓶頸，其中p為多極展開中的截斷項，新型算法將進一步將該部分的計算量減少至O（p3）［10－11］.

將格子B的相互作用列表分為6部分，分別為上列表、下列表、北列表、南列表、西列表和東列表.假設(shè)格子C在B的下列表中，見圖3，點O（x0，y0，z0）和Q（x，y，z）分別在格子C和B中，由于z＞z0成立，于是

其中

ε為殘差值，s（ε）、wk和λk見參考文獻［10－11］.d為立方體格子的棱長.

此外，還存在如下公式

假設(shè)O1為C的中心，O2為B的中心，于是式（24）可以寫為

將指數(shù)展開系數(shù)從O1點傳遞到點有

再將O2點的指數(shù)展開系數(shù)傳遞為局部展開系數(shù)有

從最終式（38）可以寫為

以上便是新型算法的主要傳遞公式，在實際計算過程中，由于格子C可能在B的其他列表中，因此需要對坐標(biāo)系進行變換，保證以上公式的成立.

可以發(fā)現(xiàn)，新型算法將初始算法中的多極矩到局部展開系數(shù)的傳遞過程分解為了三個新的傳遞過程，該三個新的傳遞過程可以進一步減少計算量［10－11］.

2 數(shù)值算例

快速多極算法是結(jié)合迭代算法進行計算，文中采用GMRES算法.多極展開的截斷項p＝10，指數(shù)傳遞中s＝8，采用C＋＋語言編程，在3.4 GHz的CPU和16.0 GB的RAM個人電腦上運行.

文中的模型以面為單元，每個面利用其曲面參數(shù)方程建立，在每個面的參數(shù)空間布置節(jié)點.如，對于一個三維曲面，其參數(shù)方程為

在建立該曲面時，以參數(shù)u和v建模，離散時對參數(shù)u和v進行離散，映射至三維實際空間.

2.1 含大量球形粒子的復(fù)合材料模擬

該算例考慮一個含有1 000個隨機分布的球形粒子的立方體模型，見圖4.基體的參數(shù)為：長L＝200 m，熱傳導(dǎo)率κmatrix＝1 W·m－1·K－1，在左端施加200 K的均勻溫度，右端施加0 K的均勻溫度.考慮了兩種不同的模型，第一種模型中，粒子的半徑r＝5 m，見圖4（a），粒子的熱傳導(dǎo)率從1 W·m－1·K－1變化到10 W·m－1·K－1.在第二種模型中，粒子的熱傳導(dǎo)率為κinclusion＝6 W·m－1·K－1，粒子的半徑從2 m變化到8 m.

圖5和圖6分別顯示了兩種模型得到的等效熱傳導(dǎo)率，為了顯示本文算法的正確性，將所得結(jié)果與Mori-Tanaka方法［14］進行對比.結(jié)果表明，本文中的算法與Mori-Tanaka方法吻合得比較好，是可行的.

2.2 含大量纖維復(fù)合材料的模擬

在本算例中模擬量如圖7所示的短纖維復(fù)合材料.在所有模型中，基體參數(shù)為L＝800 m，H＝200 m，W＝200 m，熱傳導(dǎo)率κmatrix＝1 W·m－1·K－1，在左端施加100 K均勻溫度，右端施加0 K的均勻溫度.纖維粒子的模型見圖8.在本算例中考慮了三種情況，且纖維粒子都是均勻分布，纖維粒子的參數(shù)可見表1，第四列表示纖維粒子的熱傳導(dǎo)率，第五列表示纖維粒子的個數(shù).纖維粒子的熱傳導(dǎo)率、半徑以及長度對復(fù)合材料的等效熱傳導(dǎo)率的影響可見計算結(jié)果圖9、圖10和圖11.從這三個圖中可以看出，隨著纖維粒子的熱傳導(dǎo)率、半徑以及長度的增加，復(fù)合材料的等效熱傳導(dǎo)率也隨之增加.

3 結(jié)論

推導(dǎo)了復(fù)合材料熱傳導(dǎo)問題的雜交邊界點法求解公式，并將新型快速多極算法與雜交邊界點法結(jié)合，對包含大量粒子的復(fù)合材料進行了計算模擬，該方法結(jié)合了無網(wǎng)格不需要劃分網(wǎng)格的優(yōu)點和快速算法的優(yōu)點，非常適用于計算大規(guī)模問題.同時，該算法還是一種邊界類型的數(shù)值算法，類似于邊界元一樣將問題的維數(shù)降低了一維，本文的方法還可以應(yīng)用于復(fù)合材料的力學(xué)性質(zhì)的模擬，這也是本文以后進一步的工作之一.

［1］袁光偉，杭旭登，盛志強，岳晶巖.輻射擴散計算方法若干研究進展［J］.計算物理，2009，26（4）：475－500.

［2］曾敏，王剛，謝公南，陳秋煬，王秋旺.復(fù)雜多孔介質(zhì)腔體內(nèi)自然對流換熱的數(shù)值研究［J］.計算物理，2008，25（4）：445－449.

［3］周宇，錢煒祺，何開鋒，桂業(yè)偉.同時反演材料熱傳導(dǎo)系數(shù)和比熱的算法［J］.計算物理，2011，28（5）：719－724.

［4］Zhang J，Yao Z，Li H.A hybrid boundary node method［J］.International Journal for Numerical Methods in Engineering，2002，53：751－763.

［5］Zhang J，Yao Z.The regular hybrid boundary node method for three-dimensional linear elasticity［J］.Engineering Analysis with Boundary Element，2004，28：525－534.

［6］Miao Y，Wang Y，Yu F.Development of hybrid boundary node method in two-dimensional elasticity［J］.Engineering Analysis with Boundary Elements，2005，29：703－712.

［7］Miao Y，Wang Y H.Meshless analysis for three-dimensional elasticity with singular hybrid boundary node method［J］. Applied Mathematics and Mechanics，2006，27（6）：673－681.

［8］Rokhlin V.Rapid solution of integral equations of classical potential theory［J］.Journal of Computational Physics，1985，60：187－207.

［9］Greengard L，Rokhlin V.A fast algorithm for particles simulations［J］.Journal of Computational Physics，1987，73：325－348.

［10］Greengard L，Rokhlin V.A new version of the fast multipole method for the Laplace equation in three dimensions［J］.Acta Numerica，1997，6：229－269.

［11］Cheng H，Greengard L，Rokhlin V.A fast adaptive multipole algorithm in three dimensions［J］.Journal of Computational Physics，1999，155：468－498.

［12］宋敏，覃先云，張見明，何建平.快速雜交邊界點法在水壩熱傳導(dǎo)分析中的應(yīng)用［J］.水利水電科技進展，2010，30：27－34.

［13］Yoshida K，Nishimura N，Kobayashi S.Application of new fast multipole boundary integral equation method to crack problems in 3D［J］.Engineering Analysis with Boundary Elements，2001，25：239－247.

［14］Karris A.An examination of the Mori-Tanaka effective medium—Approximation for multiphase composites［J］.Journal of Applied Mechanics，1989，56：83－88.

New Fast Multipole Hybrid Boundary Node Method for Thermal Analysis of 3D Composites

MIAO Yu，SUN Tiancan，ZHENG Junjie
（School of Civil Engineering＆Mechanics，Huazhong University of Science and Technology，Wuhan 430074，China）

Hybrid boundary node method（hybrid BNM）is applied for thermal analysis of 3D composites.A new formulation is derived for inclusion-based composites，which reduces degrees of freedom（DOFs）compared with conventional multi-domain solver.A version of fast multipole method（FMM）is coupled with new formulation for large scale analysis.Numerical examples are presented to analyze thermal behavior of composites with many inclusions.

hybrid boundary node method；3D composites；formulation；fast multipole method

date： 2013－06－06；Revised date： 2013－10－11

O24

A

1001-246X（2014）03-0335-08

2013－06－06；

2013－10－11

國家自然科學(xué)基金（51378234）及中央高校基本科研業(yè)務(wù)費專項資金（2013QN027）資助項目

苗雨（1979－），男，山東菏澤，博士后，從事計算力學(xué)及其工程應(yīng)用研究，E-mail：my_miaoyu＠163.com