新型快速多极杂交边界点法在复合材料热传导中的应用

苗 雨， 孙恬粲， 郑俊杰

（华中科技大学土木工程与力学学院，湖北武汉 430074）

新型快速多极杂交边界点法在复合材料热传导中的应用

苗 雨， 孙恬粲， 郑俊杰

（华中科技大学土木工程与力学学院，湖北武汉 430074）

将杂交边界点法应用于复合材料的热传导模拟，推导一种求解复合材料的方程，该方程减少计算自由度，效率更高.将新型快速多极算法与杂交边界点法结合进行大规模计算，数值算例中对包含大量粒子的复合材料进行模拟，结果表明快速多极杂交边界点法可行，具有一定的应用前景.

杂交边界点法；三维复合材料；方程；快速多极算法

0 引言

在工程计算中［1－3］，主要存在两大难题，①是网格生成，②是大规模计算.传统的有限元法等基于网格的算法需要在划分网格上耗费大量时间，对于复杂的问题，如复合材料等问题尤为严重.对于网格划分这一难题，本文采用近年来发展起来的一种无网格法——杂交边界点法［4－7］来解决，该方法基于移动最小二乘与修正变分原理，不仅可以将求解问题的维数降低一维，而且无需划分网格.该算法已经应用于弹性力学和弹性动力学等多个领域，具有很大的发展潜力.

然而，传统的杂交边界点法的系数矩阵为密集非对称矩阵，无法在普通微机上求解大规模问题.为了解决大规模计算这一难题，本文将新型快速多极算法（new fast multipole method）［8－11］与杂交边界点法结合，推导了相应的公式.并成功应用于复合材料热传导问题的大规模模拟中，数值算例表明，该新型快速多极杂交边界点法是有效的.

首先推导杂交边界点法求解复合材料热传导问题的公式，然后推导新型快速多极杂交边界点法求解三维热传导问题的公式，最后给出数值算例.与张见明等人［12］的方法相比，本文求解复合材料的热传导问题的公式可以减少计算自由度，且在初始快速多极算法的基础上，新型快速多极算法应用于杂交边界点法，可进一步减少计算时间.

1 复合材料中的杂交边界点法

在热传导问题中，杂交边界点法的三个独立变量分别为：域内的温度φ，边界上的温度φ～和边界上的法向温度流q～．

边界上φ～和q～采用移动最小二乘近似

其中，ΦJ（s）为节点sJ的移动最小二乘形函数.

域内的φ和q采用基本解近似

其中κ为热传导率，r（Q，sJ）为两点间的距离.

考虑一个包含n个子域的问题，对于每一个子域通过修正变分原理可得

考虑如图1中所示的复合材料模型，对于基体材料有

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

对于粒子k有

在基体与粒子的交界面上存在连续性条件，于是有

通过式（12）－（17）我们最终得到

其中，A和d通过组装边界条件得到.如果包含的粒子形状和材料相同，则Ck只需要计算一次.

1.1 快速多极杂交边界点法

快速多极杂交边界点法主要用来计算矩阵与向量相乘，对于本文的问题，矩阵与向量的相乘为如下两式之一

快速多极算法通过将远端节点的相互作用化为格子与格子之间的相互作用来减少计算量.考虑两个格子Ca和Cb分别包含Na和Nb个粒子，见图2.方便起见，我们首先只计算式（21）.

1.1.1 初始快速多极杂交边界点法

首先将基本解进行如下展开

上式在｜O1Q｜＞｜O1sJ｜时成立，其中Rn，m和Sn，m见参考文献［13］.

将式（23）代入式（21），计算Cb的贡献

其中Mn，m（O1）为多极矩，且

假设Ca和Cb属于两个大一点的格子和并且他们依旧相距很远，见图2，则有多极矩到多极矩间的传递

于是式（24）可写为

以上便是初始快速多极杂交边界点法的基本传递公式，式（22）的计算方法类似.

1.1.2 新型快速多极杂交边界点法

以上初始算法中多极矩到局部系数的传递的计算量为O（189p4），是初始算法的计算瓶颈，其中p为多极展开中的截断项，新型算法将进一步将该部分的计算量减少至O（p3）［10－11］.

将格子B的相互作用列表分为6部分，分别为上列表、下列表、北列表、南列表、西列表和东列表.假设格子C在B的下列表中，见图3，点O（x0，y0，z0）和Q（x，y，z）分别在格子C和B中，由于z＞z0成立，于是

其中

ε为残差值，s（ε）、wk和λk见参考文献［10－11］.d为立方体格子的棱长.

此外，还存在如下公式

假设O1为C的中心，O2为B的中心，于是式（24）可以写为

将指数展开系数从O1点传递到点有

再将O2点的指数展开系数传递为局部展开系数有

从最终式（38）可以写为

以上便是新型算法的主要传递公式，在实际计算过程中，由于格子C可能在B的其他列表中，因此需要对坐标系进行变换，保证以上公式的成立.

可以发现，新型算法将初始算法中的多极矩到局部展开系数的传递过程分解为了三个新的传递过程，该三个新的传递过程可以进一步减少计算量［10－11］.

2 数值算例

快速多极算法是结合迭代算法进行计算，文中采用GMRES算法.多极展开的截断项p＝10，指数传递中s＝8，采用C＋＋语言编程，在3.4 GHz的CPU和16.0 GB的RAM个人电脑上运行.

文中的模型以面为单元，每个面利用其曲面参数方程建立，在每个面的参数空间布置节点.如，对于一个三维曲面，其参数方程为

在建立该曲面时，以参数u和v建模，离散时对参数u和v进行离散，映射至三维实际空间.

2.1 含大量球形粒子的复合材料模拟

该算例考虑一个含有1 000个随机分布的球形粒子的立方体模型，见图4.基体的参数为：长L＝200 m，热传导率κmatrix＝1 W·m－1·K－1，在左端施加200 K的均匀温度，右端施加0 K的均匀温度.考虑了两种不同的模型，第一种模型中，粒子的半径r＝5 m，见图4（a），粒子的热传导率从1 W·m－1·K－1变化到10 W·m－1·K－1.在第二种模型中，粒子的热传导率为κinclusion＝6 W·m－1·K－1，粒子的半径从2 m变化到8 m.

图5和图6分别显示了两种模型得到的等效热传导率，为了显示本文算法的正确性，将所得结果与Mori-Tanaka方法［14］进行对比.结果表明，本文中的算法与Mori-Tanaka方法吻合得比较好，是可行的.

2.2 含大量纤维复合材料的模拟

在本算例中模拟量如图7所示的短纤维复合材料.在所有模型中，基体参数为L＝800 m，H＝200 m，W＝200 m，热传导率κmatrix＝1 W·m－1·K－1，在左端施加100 K均匀温度，右端施加0 K的均匀温度.纤维粒子的模型见图8.在本算例中考虑了三种情况，且纤维粒子都是均匀分布，纤维粒子的参数可见表1，第四列表示纤维粒子的热传导率，第五列表示纤维粒子的个数.纤维粒子的热传导率、半径以及长度对复合材料的等效热传导率的影响可见计算结果图9、图10和图11.从这三个图中可以看出，随着纤维粒子的热传导率、半径以及长度的增加，复合材料的等效热传导率也随之增加.

3 结论

推导了复合材料热传导问题的杂交边界点法求解公式，并将新型快速多极算法与杂交边界点法结合，对包含大量粒子的复合材料进行了计算模拟，该方法结合了无网格不需要划分网格的优点和快速算法的优点，非常适用于计算大规模问题.同时，该算法还是一种边界类型的数值算法，类似于边界元一样将问题的维数降低了一维，本文的方法还可以应用于复合材料的力学性质的模拟，这也是本文以后进一步的工作之一.

［1］袁光伟，杭旭登，盛志强，岳晶岩.辐射扩散计算方法若干研究进展［J］.计算物理，2009，26（4）：475－500.

［2］曾敏，王刚，谢公南，陈秋炀，王秋旺.复杂多孔介质腔体内自然对流换热的数值研究［J］.计算物理，2008，25（4）：445－449.

［3］周宇，钱炜祺，何开锋，桂业伟.同时反演材料热传导系数和比热的算法［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］宋敏，覃先云，张见明，何建平.快速杂交边界点法在水坝热传导分析中的应用［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

国家自然科学基金（51378234）及中央高校基本科研业务费专项资金（2013QN027）资助项目

苗雨（1979－），男，山东菏泽，博士后，从事计算力学及其工程应用研究，E-mail：my_miaoyu＠163.com