Mechanical Analysis of 3D Composite Materials by Hybrid Boundary Node Method

Yu Miao,Zhe Chen,Qiao Wang,2and Hongping Zhu

1 Introduction

Composite materials become more and more important in industrial projects.There has been much interest in the modeling and simulations of composites in order to characterize their mechanical properties for engineering applications.Much effort,time and expense would be saved if the properties of composites could be predicted accurately.Many numerical models based on finite element method(FEM)and boundary element method(BEM)have been developed to simulate composite materials so far.Hou and Wu(1997)studied a multiscale finite element method for solving a class of elliptic problems arising from composite materials and flows in porous media.A numerical homogenization technique based on the FEM with representative volume element(RVE)was used to evaluate the effective material properties with periodic boundary conditions by Kari,Berger,Rodriguez-Ramos and Gabbert(2007).Segurado and LLorca(2004)developed a new 3D quadratic interface finite element to simulate fracture in composite materials Hybrid/Mixed finite elements and mixed methods[Bishay and Atluri(2012);Bishay,Sladek,Sladek and Atluri(2012);Dong,Alotaibi,Mohiuddine and Atluri(2014);Dong,El-Gizawy,Juhany and Atluri(2014)]were developed recently and applied for mechanical modeling of composites.Eischen and Torquato(1993)applied the BEM to determine the effective elastic moduli of continuum models of composite materials.An advanced 3D boundary element method was proposed by Chen and Liu(2005)for characterizations of composite materials.Liu,Nishimura,Otani,Takahashi,Chen and Munakata(2005)developed a new BEM for 3D analysis of fiber-reinforced composites based on a rigid-inclusion model.A new BEM for the numerical analysis of mechanical properties in 3D particle-reinforced composites was proposed by Wang and Yao(2005).Yao,Xu,Wang and Zheng(2009)applied the BEM on the simulation of carbon nanotube(CNT)composites,including the simulation of elastic,thermal and electric properties.Since the BEM has a dense matrix in the final system equation,some researchers[Nishimura and Liu(2004);Liu,Nishimura and Otani(2005);Liu,Nishimura,Otani,Takahashi,Chen and Munakata(2005);Wang and Yao(2005);Wang and Yao(2008);Yao,Xu,Wang and Zheng(2009)]coupled the BEM with fast multipole method(FMM)[Rokhlin(1985)]in order to simulate composites with large scale computation.

Although the FEM and BEM are widely investigated and have been applied in many areas,they may not be convenient and efficient to simulate composites or fracture problems[Dong and Atluri(2013a);Dong and Atluri(2013b);Dong and Atluri(2013c)]because of the cells caused by the mesh procedure.In order to overcome the difficulty caused by meshing,meshless or meshfree methods are widely investigated in recent years.Many meshless methods have been proposed so far,including the smoothed particle hydrodynamics(SPH)method[Gingold and Monaghan(1977);Lucy(1977)],the reproducing kernel particle methods(RKPM)[Liu,Jun and Zhang(1995)],the hp-clouds method[Liszka,Duarte and Tworzydlo(1996)],the element free Galerkin method(EFG)[Belytschko,Lu and Gu(1994)],the meshless local Petrov-Galerkin(MLPG)approach[Atluri and Zhu(1998)],the boundary node method(BNM)[Mukherjee and Mukherjee(1997)],the Galerkin boundary node method(GBNM)[Li and Zhu(2009)],the boundary face method(BFM)[Zhang,Qin,Han and Li(2009)],the Method of Finite Sphere[De,Hong and Bathe(2003)],the boundary point interpolation method[Gu and Liu(2002)]and the hybrid boundary node method(Hybrid BNM)[Zhang and Yao(2001);Zhang,Yao and Li(2002)].

The Hybrid BNM[Zhang,Tanaka and Matsumoto(2004a);Zhang and Yao(2004)]is based on the moving least squares(MLS)approximation and the hybrid displacement variational formula.It not only reduces the spatial dimensions by one like boundary element method(BEM)or BNM,but also does not require any cells neither for interpolation of the solution variables nor for the boundary integration.In fact,the Hybrid BNM requires only discrete nodes located on the surface of domain and its parametric representation.The Hybrid BNM has been developed by Miao et al[Miao,Wang and Yu(2005);Miao and Wang(2006)]and applied to some other areas,such as elastodynamics problems[Miao,Wang,Liao and Zheng(2009)],Helmholtz problems[Miao,Wang and Wang(2009)]and multi-domain problems[Wang,Zheng,Miao and Lv(2011)].

The meshless methods are very suitable for simulation of composite materials and some literatures have reported relative work.Zhang et al.[Zhang,Tanaka and Matsumoto(2004b);Zhang and Tanaka(2008)]simulated the CNT based composites by Hybrid BNM based on a simplified model.Singh,Tanaka and Endo(2007)applied the EFG in the thermal analysis of CNT based composites.Wang,Miao and Zhu(2013a)used a fast multipole Hybrid BNM(FM-HBNM)to simulate the composites in 3D elasticity.There are mainly two models while simulating composites by Hybrid BNM so far.One is the multi-domain solver[Wang,Zheng,Miao and Lv(2011);Wang,Miao and Zhu(2013a)],which is based on the continuity and equilibrium conditions across the interface.The other one is the simplified approach[Zhang,Tanaka and Matsumoto(2004b)],which treats the CNTs as heat superconductors.The multi-domain solver can be applied to general composites since no assumption is added and the unknown vectors for each material are assembled in the final system equation.The simplified approach is very suitable for the case when the inclusions are heat superconductors or rigid bodies and it can reduce the total degrees of freedom(DOFs)substantially,however,it is not suitable to simulate general composites as multi-domain solver.

In this paper,an improved multi-domain model for 3D composites with inclusions in elasticity problems base on Hybrid BNMis proposed.It has already been applied to thermal analysis of composites[Wang,Miao and Zhu(2013b)]successfully and we extend it to elasticity problems in this paper.In the improved model,no other assumption is added and it can be applied to general case like the multi-domain solver.Much important,it can reduce the total DOFs in the final system equation as the simplified approach,and we will find out only the matrix domain is needed to be modeled in some special case.

This paper is organized as follows.The hybrid BNM is introduced in the second section and the multi-domain solver for 3D elasticity problems is reviewed in the third section.The formulation for the improved multi-domain model is derived in the fourth section.Then numerical examples are given in the fifth section.

2 The single-domain Hybrid BNM for 3D elasticity problems

In this section,the Hybrid BNM for3D elasticity problems[MiaoandWang(2006)]is reviewed.The Hybrid BNM is based on a modified variational principle.In 3D elasticity problems,functions in the modified variational principle that assumed to be independent are:displacementsand tractionson the boundary and displacementsulinside the domain.Consider a domain ? enclosed by Γ=Γu+Γtwithandare the prescribed displacements and tractions,respectively.

whereNis the number of nodes located on the surfaceandare nodal values,and ΦJ(s)is the shape function of the MLS approximation,corresponding to node.

The u and t inside the domain can be approximated by fundamental solutions as

wherer=r(sJ,Q)andQis the field point while sJis the source point.

As the modified variational principle holds both in the whole domain and any subdomain,the local sub-domain around each node can be taken into consideration.

Then the following set of equations,expressed in matrix form,are given as[Zhang and Yao(2004);Miao and Wang(2006)]

wherehI(Q)is a weight function,ΓIis a regularly shaped local region around node sIin the parametric representation space of the boundary surface.Therefore,the integrals in Equations(9),(10)and(11)can be computed without using boundary elements.

For a general problem,eitherorcan be known at each node on the boundary and by rearranging Equations(7)and(8),a final algebraic equation in terms of x only can be obtained as below:

For the node sI,ifis known,select the correspond row in U to A,otherwise,select the correspond row in T to A,and the corresponding term of d comes from the matrix-vector product ofor.Then the unknown vector x is obtained by solving the final algebraic equation.The nodal valuesandon the boundary can be computed by the back-substitution of x into Equations(7)and(8)

3 The multi-domain Hybrid BNM for 3D elasticity problems

In this section,the multi-domain Hybrid BNM for 3D elasticity problems[Wang,Zheng,Miao and Lv(2011)]is introduced.Consider a problem withnsub-domains andNinodes are distributed on the bounding surface of subdomain-i.Then for each sub-domain,we can have the following equations from Equations(7)and(8)

where superscriptiindicates the subdomaini.

The boundary nodes of subdomain-ican be sorted intongroups according to the locations of the nodes.Groupicontains nodes that belong exclusively in subdomain-iand groupjcontains nodes that are on the interface with subdomain-j,wherejis from 1 tonandj/=i.Accordingly,Equations(16)and(17)are partitioned into blocked matrix equations as

where the superscriptidenotes the subdomain-i;the subscriptk=1,2,...ndenotes that the prescribed quantities are associated with the nodes in groupk.The double subscriptsij,iandj=1,2,...nis used to convey the U T in Equations(16)and(17),by which the prescribed coefficient matrix blocks are computed,belong to groupiandj,respectively.

At the interfaces between sub-domainiandj,the following continuity conditions exist.For the displacements

and for the tractions

Consider a problem with two sub-domains,then for subdomain-1,we can have

For subdomain-2,we can have

Using Equations(20)and(21),Equations(22),(23),(24)and(25)can be assembled to an overall matrix equation as

Similarly,for three subdomains,the final equation for the entire domain can be written as

Equation(28)can be rewritten as

where A2stands for the coefficient matrix of the two subdomains in the left hand side of Equation(27)and D2stands for the vector in the right hand side of Equation(27).Then fornsubdomains,the final equation for the entire domain can be written as

where An?1and Dn?1stand for the coefficient matrix and right hand vector of then-1 subdomains.

The unknown vector x can be solved by the final matrix equation,and then by backsubstitution x to matrix equations for each sub-domain,the boundary unknowns can be obtained either on the interfaces or the external boundary surfaces.

4 Improved multi-domain Hybrid BNM for 3D elasticity problems

For composite materials,the multi-domain solver introduced in section 3 is a natural way to be chosen.Now consider a matrix with inclusions as showed in Figure 1(a).S0is the sub-domain of matrix andS1,S2,...,Snare sub-domains of the inclusions.For the matrix,Equations(18)and(19)can be rewritten as

Figure 1:(a)The model of matrix with inclusions;(b)Three kinds of inclusions.

where scripts 0 andk,k=1,...,nindicate the subdomain of matrix and thek-th inclusion,respectively.

Three types of inclusions can be considered as shown in Figure 1(b)and the inclusion type 2 is totally in the matrix and it is hollow.Then for the general case,the system equations for thek-th inclusion domain can be written as

If the inclusions are solid and totally embed in the matrix(the inclusion type 3 in Figure 1(b)),then Equations(33)and(34)can be rewritten as

One can obtain the final system equation for the whole domain as

Equation(37)is the conventional multi-domain solver for the composites of matrix with inclusions.All the unknown parametersxabout the matrix and inclusions are solved in this equation.

Now we derive the improved formulation for composites with inclusions,which is based on the multi-domain solver.

From Equation(31)one can have

From Equation(33)one can have

where

Then from Equation(39)one can obtain

And from Equation(33)one can obtain

At the interface between the matrix and inclusion,both the displacement and traction must be continuous.If we use the same set of nodes distributed on the interface,we can obtain

Substituting Equations(42)and(43)into Equation(41)and using Equation(38),one can have

where I is the identity operator.

Using the boundary conditions,for the matrix domain one can have

For thek-th inclusion domain one can have

Submitting Equation(45)into Equation(47)and using Equation(44)we can obtain

Then by assembling Equations(48)and(49)into Equation(46)one can obtain the final system equation as

For the case when the inclusions are solid and totally embed in the matrix,Equation(50)can be rewritten as

If the inclusions are of the same size and material,the corresponding incidence matrices Ckare also identical.Therefore,Ckis required to be formed only once.If the inclusions are of the same shape,Ckcan also be calculated only once because of the similar relationship between inclusion phases of different size and material.If we assume the Poisson’s ratios of the inclusions are the same,one can have the following relationship

whereEkandrkare Young’s module and radius of thek-th inclusion,respectively.In Equations(50)and(59),the unknowns on the interfaces are assembled only once compared with Equation(37),thus,the total unknowns in the final system equation can be reduced by the improved multi-domain solver.In equation(50),if the inclusions are of the same shape and the same boundary condition,Equations(52),(54),(55),(56),(57)and(58)can also be computed once and the computational resource,both the time and memory required,can be saved further.Compared with the simplified approach proposed by Zhang,Tanaka and Matsumoto(2004b),no assumption is needed and it can be applied to analyze more general composites.The improved multi-domain solver is also very suitable to couple with FMM for large scale computation.

5 Numerical examples

The proposed technique has been implemented in C++program.Three numerical examples are studied in this section.For the purpose of error estimation,a formula is defined as

whereandrefer to the exact and numerical solutions respectively and|u|maxis the maximum value ofuoverNnodes.

In all the examples,Gaussian elimination method is used to solve the final system equations and both the multi-domain solver and improved multi-domain solver are used.

5.1 Validation of the method

In this example,a hollow sphere shown in Figure 2 is considered,the outer radius of the sphere is 10 and the inner radius is 1,the Young’s moduleE=1.0 and the Poisson’s ratiov=0.25 Direchlet boundary conditions on all the surfaces are imposed,according to the following exact solution:

Figure 2:A hollow sphere with three subdomains.

In order to apply the multi-domain and improved multi-domain solver,the hollow sphere is divided into three subdomains with the same materials,see Figure 2.For the subdomain 2,the outer radius is 2 and the inner radius is 1.For subdomain 3,the radius is 2 and its center located at point(6,0,0).

In the model,the spherical face is considered as one face and all the spherical faces have the same number of nodes distributed on them.Several different node arrangements are used to show the efficiency and accuracy of the technique.Figure 3 shows the relative error ofuxcomputed by the nodes on line(1,0,0)to(10,0,0).From Figure 3 one can observe that the improved multi-domain solver has the same precision as the conventional multi-domain solver.The CPU time for solving the final system equation is shown in Figure 4 and we can observe that the improved multi-domain is more efficient than the conventional multi-domain solver since it can reduce the total DOFs in the final system equation.

5.2 Cube with eight spherical inclusions

A cube with eight spherical inclusions shown in Figure 5 is considered in this example.The parameters of the cube are:the lengtha=2,the Young’s moduleE=1.0,Poisson’s ratiov=0.25.The cube is roller supported on the left end and a uniform distributed traction is loaded on outer normal direction of the right end.These boundary conditions allow us to estimate the effective Young’s module in the axial direction as

Figure 3:Relative error.

Figure 4:CPU time.

whereEeffis the estimated effective Young’s module,Lis the length,σ is average stress of the two ends and Δuis the difference of displacement between the two ends.

Two models are considered.In the first model,the parameters of the inclusions are:the radiusr=0.2 and the Young’s moduleEinchanges from 1 to 10.In the second model,the Young’s moduleEin=8 and the radius of the inclusion changes from 0.1 to 0.4.In all the models,the Poisson’s ratiov=0.25,600 nodes are distributed on the outer faces and 172 nodes are distributed on each interface.Table 1 shows the information of total unknowns in the final system equation and the computational time.From Table 1,one can indicate that the number of DOFs in the improved multi-domain solver is reduced to nearly one half of that in the conventional multidomain solver.Thus,the improved multi-domain solver is more efficient and From Table 1,one can observe that the CPU time has been saved much by improve solver.Figure 6 and Figure 7 show the effective Young’s modules while the Young’s module of the inclusion and radius of the inclusion changes,respectively.In order to show the validation of the method,the Mori-Tanaka method[Karris(1989)]is also used.From Figure 6 and Figure 7 one can observe that the results obtained by the improved multi-domain solver and conventional multi-domain solver are the same and both of them have good agreement with the Mori-Tanaka method.It should be noted that Mori-Tanaka method can only be used to estimate the effective module of composite for some special cases,it cannot give the field of stress or displacement under boundary conditions.However,by applying numerical methods like introduced in this paper,these aims can be easily achieved.

If the matrix has more inclusions,the Hybrid BNM can also be applied easily.However,without fast algorithms the technique cannot solve large scale problems[Wang,Miao and Zheng(2010);Wang,Miao,Zhu and Zhang(2012);Wang,Miao and Zhu(2013a);Wang,Miao and Zhu(2013c);Miao,Wang,Zhu and Li(2014)]

Table 1:Time information of cub with 8 spherical inclusions.

Figure 5:Model of cub with 8 uniformly distributed spherical inclusions.

Figure 6:Effective Young’s module with Young’s module of the inclusion.

Figure 7:Effective Young’s module with radius of the inclusion.

5.3 A square RVE

A square RVE model with an inclusion embedded is considered in this example.The geometry is shown in Figure 8.The parameters of the matrix are:the lengthL=100,the widthW=20,the heightH=20,the Young’s moduleEmatrix=1.0 and the Poisson’s ratiov=0.25.The RVE is roller supported on one end and applied with a uniform normal displacementuy=100 on the other end.The parameters of the inclusion are:the lengthLc=50,the outer radiusr=5,the thicknesst=0.4 and the Poisson’s ratiov=0.25.If the Young’s module of the inclusion is supposed to be much higher than the matrix,this model can be considered as a CNT based composites.

The ratio of the Young’s modules between the inclusion and the matrix is defined asc=Einclusion/Ematrix.Numerical results for the displacements at linesx=0,z=4.8(through the inclusion)andx=0,z=5.1(very near the inclusion)with different values of ratiocare presented in Figure 9.In this example,we useMSandIMSdenote themulti-domain solverandimproved multi-domain solver,respectively.The results obtained by the multi-domain solver and improved multi-domain solver are the same.From Figure 8,we can observe that the displacements at the two lines are very close with a given ratioc.The difference of the displacements within the inclusion(r=4.8)and near the inclusion(r=5.1)becomes smaller as the increase of ratioc.The displacements become almost constant along the length of the inclusion whenc=1000 and the results are very similar to a RVE with rigid inclusion.

Figure 8:A square RVE model.

Figure 10 shows the displacements at linex=0,z=4.8 with four different thicknesstof the inclusion whilec=200.Whilet=5.0,the inclusion is solid and only its outer face is needed.In Figure 10,we can observe that thickness of the inclusion has little influence on the displacements at the evaluated line.One of the reasons is that the inclusion is nearly a rigid body in this case and the displacements of the inclusion is the same.Figure 11 shows the effective Young’s module of the RVE with various thicknesstand ratio of the Young’s modules between the inclusion and the matrix.From Figure 11 one can observe that the difference of the effective Young’s modules becomes littler and littler as the increase of the ratiocand the effective Young’s modules become very close for the various thicknesstwhen the ratiocis more than 6000.From this example we can find out that the thickness of the inclusion has little influence on the effective Young’s module of the RVE when the Young’s module of the inclusion is much higher than the matrix,and only the outer face of the inclusion is needed be modeled,which can reduce much resource further.

Figure 9:Displacements along the axial lines(t=0.4).

Figure 10:Displacements along the axial lines(c=2000).

Figure 11:Effective Young’s module of the RVE with various t and c.

6 Conclusions

In this paper,an improved multi-domain solver is proposed for the mechanical analysis of composites with inclusions by the Hybrid BNM.Continuity conditions are used in the improved multi-domain solver as the conventional multi-domain solver.The unknown vectors on the interfaces are needed to be assembled only once in the final system equation of improved solver,thus,the computational time and memory required can be reduced.

The improved multi-domain model proposed can be applied to more general composites than the simplified approach,and especial suitable for the composites when the inclusions are solid and totally embedded in the matrix domain.Composites which similar as rigid inclusion based composites are also discussed in the paper and we can find out that the thickness of the rigid inclusion has little influence on the effective Young’s module.

The improved multi-domain solver is very suitable for large-scale computation and coupling it with fast algorithm,such as FMM is a subject of future research.

The difference between the numerical and experimental results has been observe in the work of Yao,Xu,Wang and Zheng(2009)It means that the numerical method has to be improved by considering more practical microstructural factors,including the imperfect interface conditions and thermal contact between fibers.These researches will be investigated in the future.

Acknowledgement:Financial support for the Project from the Natural Science Foundation of China(No.51378234)the National Basic Research Program of China(973 Program:2011CB013800)and the Fundamental Research Funds for the Central Universities(No.2014YQ008 and No.2013QN027).

Atluri,S.;Zhu,T.(1998):A new meshless local Petrov-Galerkin(MLPG)approach in computational mechanics.Computational mechanics,vol.22,pp.117-127.

Belytschko,T.;Lu,Y.Y.;Gu,L.(1994):Element-free Galerkin methods.International journal for numerical methods in engineering,vol.37,pp.229-256.

Bishay,P.;Atluri,S.(2012):High-Performance3DHybrid/Mixed,andSimple3D Voronoi Cell Finite Elements,for Macro-&Micro-mechanical Modeling of Solids,Without Using Multi- field Variational Principles.Computer Modeling in Engineering&Sciences(CMES),vol.84,pp.41-97.

Bishay,P.L.;Sladek,J.;Sladek,V.;Atluri,S.N.(2012):Analysis of Functionally Graded Magneto-Electro-Elastic Composites Using Hybrid/Mixed Finite Elements and Node-Wise Material Properties.CMC:Computers,Materials&Continua,vol.29,pp.213-262.

Chen,X.;Liu,Y.(2005):An advanced 3D boundary element method for characterizations of composite materials.Engineering analysis with boundary elements,vol.29,pp.513-523.

De,S.;Hong,J.-W.;Bathe,K.(2003):On the method of finite spheres in applications:towards the use with ADINA and in a surgical simulator.Computational Mechanics,vol.31,pp.27-37.

Dong,L.;Alotaibi,A.;Mohiuddine,S.;Atluri,S.(2014):Computational methods in engineering:a variety of primal&mixed methods,with global&local interpolations,for well-posed or ill-Posed BCs.CMES:Computer Modeling in Engineering&Sciences,vol.99,pp.1-85.

Dong,L.;Atluri,S.N.(2013a):Fracture&fatigue analyses:SGBEM-FEM or XFEM?Part 1:2D structures.CMES:Computer Modeling in Engineering&Sciences,vol.90,pp.91-146.

Dong,L.;Atluri,S.N.(2013b):Fracture&fatigue analyses:SGBEM-FEM or XFEM?Part 2:3D solids.CMES:Computer Modeling in Engineering&Sciences,vol.90,pp.379-413.

Dong,L.;Atluri,S.N.(2013c):SGBEM Voronoi Cells(SVCs),with embedded arbitrary-shaped inclusions,voids,and/or cracks,for micromechanical modeling of heterogeneous materials.CMC:Computers,Materials&Continua,vol.33,pp.111-154.

Dong,L.;El-Gizawy,A.S.;Juhany,K.A.;Atluri,S.N.(2014):A simple locking-alleviated 4-node mixed-collocation finite element with over-integration,for homogeneous or functionally-graded or thick-section laminated composite beams.CMC:Computers,Materials&Continua,vol.40,pp.49-77.

Eischen,J.;Torquato,S.(1993):Determining elastic behavior of composites by the boundary element method.Journal of applied physics,vol.74,pp.159-170.

Gingold,R.A.;Monaghan,J.J.(1977):Smoothed particle hydrodynamics theory and application to non-spherical stars.Monthly notices of the royal astronomical society,vol.181,pp.375-389.

Gu,Y.;Liu,G.(2002):A boundary point interpolation method for stress analysis of solids.Computational Mechanics,vol.28,pp.47-54.

Hou,T.Y.;Wu,X.-H.(1997):A multiscale finite element method for elliptic problems in composite materials and porous media.Journal of computational physics,vol.134,pp.169-189.

Kari,S.;Berger,H.;Rodriguez-Ramos,R.;Gabbert,U.(2007):Computational evaluation of effective material properties of composites reinforced by randomly distributed spherical particles.Composite structures,vol.77,pp.223-231.

Karris,A.(1989):An Examination of the Mori-Tanaka Effective Medium–Approximation for Multiphase Composites.Journal of Applied Mechanics,vol.56,pp.83-88.

Li,X.;Zhu,J.(2009):A Galerkin boundary node method and its convergence analysis.Journal of computational and applied mathematics,vol.230,pp.314-328.

Liszka,T.;Duarte,C.;Tworzydlo,W.(1996):hp-Meshless cloud method.Computer Methods in Applied Mechanics and Engineering,vol.139,pp.263-288.

Liu,W.K.;Jun,S.;Zhang,Y.F.(1995):Reproducing kernel particle methods.International Journal for Numerical Methods in Fluids,vol.20,pp.1081-1106.

Liu,Y.;Nishimura,N.;Otani,Y.(2005):Large-scale modeling of carbon-nanotube composites by a fast multipole boundary element method.Computational Materials Science,vol.34,pp.173-187.

Liu,Y.;Nishimura,N.;Otani,Y.;Takahashi,T.;Chen,X.;Munakata,H.(2005):A fast boundary element method for the analysis of fiber-reinforced composites based on a rigid-inclusion model.Journal of applied mechanics,vol.72,pp.115-128.

Lucy,L.B.(1977):A numerical approach to the testing of the fission hypothesis.The astronomical journal,vol.82,pp.1013-1024.

Miao,Y.;Wang,Q.;Liao,B.;Zheng,J.(2009):A dual hybrid boundary node method for 2D elastodynamics problems.Computer Modeling in Engineering and Sciences(CMES),vol.53,pp.1.

Miao,Y.;Wang,Q.;Zhu,H.;Li,Y.(2014):Thermal analysis of 3D composites by a new fast multipole hybrid boundary node method.Computational Mechanics,vol.53,pp.77-90.

Miao,Y.;Wang,Y.-h.(2006):Meshless analysis for three-dimensional elasticity with singular hybrid boundary node method.Applied Mathematics and Mechanics,vol.27,pp.673-681.

Miao,Y.;Wang,Y.;Wang,Y.(2009):A meshless hybrid boundary-node method for Helmholtz problems.Engineering analysis with boundary elements,vol.33,pp.120-127.

Miao,Y.;Wang,Y.;Yu,F.(2005):Development of hybrid boundary node method in two-dimensional elasticity.Engineering analysis with boundary elements,vol.29,pp.703-712.

Mukherjee,Y.X.;Mukherjee,S.(1997):The boundary node method for potential problems.International Journal for Numerical Methods in Engineering,vol.40,pp.797-815.

Nishimura,N.;Liu,Y.(2004):Thermal analysis of carbon-nanotube composites using a rigid-line inclusion model by the boundary integral equation method.Computational Mechanics,vol.35,pp.1-10.

Rokhlin,V.(1985):Rapid solution of integral equations of classical potential theory.Journal of Computational Physics,vol.60,pp.187-207.

Segurado,J.;LLorca,J.(2004):A new three-dimensional interface finite element to simulate fracture in composites.International journal of solids and structures,vol.41,pp.2977-2993.

Singh,I.;Tanaka,M.;Endo,M.(2007):Thermal analysis of CNT-based nanocomposites by element free Galerkin method.Computational Mechanics,vol.39,pp.719-728.

Wang,H.;Yao,Z.(2005):A new fast multipole boundary element method for large scale analysis of mechanical properties in 3D particle-reinforced composites.Computer Modeling in Engineering&Sciences,vol.7,pp.85-95.

Wang,H.;Yao,Z.(2008):Arigid-fiber-based boundary element model for strength simulation of carbon nanotube reinforced composites.CMES:Computer Modeling in Engineering&Sciences,vol.29,pp.1-13.

Wang,Q.;Miao,Y.;Zheng,J.(2010):The hybrid boundary node method accel-erated by fast multipole expansion technique for 3D elasticity.Computer Modeling in Engineering and Sciences,vol.70,pp.123-151.

Wang,Q.;Miao,Y.;Zhu,H.(2013a):A fast multipole hybrid boundary node method for composite materials.Computational Mechanics,vol.,pp.885-897.

Wang,Q.;Miao,Y.;Zhu,H.(2013b):A new formulation for thermal analysis of composites by hybrid boundary node method.International Journal of Heat and Mass Transfer,vol.64,pp.322-330.

Wang,Q.;Miao,Y.;Zhu,H.(2013c):Numerical Simulation of Heat Conduction Problems by a New Fast Multipole Hybrid Boundary-Node Method.Numerical Heat Transfer,Part B:Fundamentals,vol.64,pp.436-459.

Wang,Q.;Miao,Y.;Zhu,H.;Zhang,C.(2012):An O(N)Fast Multipole Hybrid Boundary Node Method for 3D Elasticity.Computers Materials and Continua,vol.28,pp.1-25.

Wang,Q.;Zheng,J.;Miao,Y.;Lv,J.(2011):The multi-domain hybrid boundary node method for 3D elasticity.Engineering Analysis with Boundary Elements,vol.35,pp.803-810.

Yao,Z.-H.;Xu,J.-D.;Wang,H.-T.;Zheng,X.-P.(2009):Simulation of CNT composites using fast multipole BEM.Journal of Marine Science and Technology,vol.17,pp.194-202.

Zhang,J.;Qin,X.;Han,X.;Li,G.(2009):A boundary face method for potential problems in three dimensions.International journal for numerical methods in engineering,vol.80,pp.320-337.

Zhang,J.;Tanaka,M.(2008):Fast HdBNM for large-scale thermal analysis of CNT-reinforced composites.Computational Mechanics,vol.41,pp.777-787.

Zhang,J.;Tanaka,M.;Matsumoto,T.(2004a):Meshless analysis of potential problem sin three dimensions with the hybrid boundary node method.International Journal for Numerical Methods in Engineering,vol.59,pp.1147-1166.

Zhang,J.;Tanaka,M.;Matsumoto,T.(2004b):A simplified approach for heat conduction analysis of CNT-based nano-composites.Computer methods in applied mechanics and engineering,vol.193,pp.5597-5609.

Zhang,J.;Yao,Z.(2001):Meshless regular hybrid boundary node method.CMESComputer Modeling in Engineering and Sciences,vol.2,pp.307-318.

Zhang,J.;Yao,Z.(2004):The regular hybrid boundary node method for three dimensional linear elasticity.Engineering analysis with boundary elements,vol.28,pp.525-534.

Zhang,J.;Yao,Z.;Li,H.(2002):A hybrid boundary node method.International Journal for Numerical Methods in Engineering,vol.53,pp.751-763.