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

苗 雨， 孫恬粲， 鄭俊杰

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

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

苗 雨， 孫恬粲， 鄭俊杰

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

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

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

0 引言

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

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

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

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

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

邊界上φ～和q～采用移動(dòng)最小二乘近似

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

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

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

考慮一個(gè)包含n個(gè)子域的問(wèn)題，對(duì)于每一個(gè)子域通過(guò)修正變分原理可得

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

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

對(duì)于粒子k有

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

通過(guò)式（12）－（17）我們最終得到

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

1.1 快速多極雜交邊界點(diǎn)法

快速多極雜交邊界點(diǎn)法主要用來(lái)計(jì)算矩陣與向量相乘，對(duì)于本文的問(wèn)題，矩陣與向量的相乘為如下兩式之一

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

1.1.1 初始快速多極雜交邊界點(diǎn)法

首先將基本解進(jìn)行如下展開(kāi)

上式在｜O1Q｜＞｜O1sJ｜時(shí)成立，其中Rn，m和Sn，m見(jiàn)參考文獻(xiàn)［13］.

將式（23）代入式（21），計(jì)算Cb的貢獻(xiàn)

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

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

于是式（24）可寫(xiě)為

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

1.1.2 新型快速多極雜交邊界點(diǎn)法

以上初始算法中多極矩到局部系數(shù)的傳遞的計(jì)算量為O（189p4），是初始算法的計(jì)算瓶頸，其中p為多極展開(kāi)中的截?cái)囗?xiàng)，新型算法將進(jìn)一步將該部分的計(jì)算量減少至O（p3）［10－11］.

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

其中

ε為殘差值，s（ε）、wk和λk見(jiàn)參考文獻(xiàn)［10－11］.d為立方體格子的棱長(zhǎng).

此外，還存在如下公式

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

將指數(shù)展開(kāi)系數(shù)從O1點(diǎn)傳遞到點(diǎn)有

再將O2點(diǎn)的指數(shù)展開(kāi)系數(shù)傳遞為局部展開(kāi)系數(shù)有

從最終式（38）可以寫(xiě)為

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

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

2 數(shù)值算例

快速多極算法是結(jié)合迭代算法進(jìn)行計(jì)算，文中采用GMRES算法.多極展開(kāi)的截?cái)囗?xiàng)p＝10，指數(shù)傳遞中s＝8，采用C＋＋語(yǔ)言編程，在3.4 GHz的CPU和16.0 GB的RAM個(gè)人電腦上運(yùn)行.

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

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

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

該算例考慮一個(gè)含有1 000個(gè)隨機(jī)分布的球形粒子的立方體模型，見(jiàn)圖4.基體的參數(shù)為：長(zhǎng)L＝200 m，熱傳導(dǎo)率κmatrix＝1 W·m－1·K－1，在左端施加200 K的均勻溫度，右端施加0 K的均勻溫度.考慮了兩種不同的模型，第一種模型中，粒子的半徑r＝5 m，見(jiàn)圖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］進(jìn)行對(duì)比.結(jié)果表明，本文中的算法與Mori-Tanaka方法吻合得比較好，是可行的.

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

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

3 結(jié)論

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

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

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

［3］周宇，錢(qián)煒祺，何開(kāi)鋒，桂業(yè)偉.同時(shí)反演材料熱傳導(dǎo)系數(shù)和比熱的算法［J］.計(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］宋敏，覃先云，張見(jiàn)明，何建平.快速雜交邊界點(diǎn)法在水壩熱傳導(dǎo)分析中的應(yīng)用［J］.水利水電科技進(jìn)展，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

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

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