Y.Wang,E.Lü,2,J.Zhaoand J.Guo

Meshfree Method for the Topological Design of Microstructural Composites

Y.Wang1,E.Lü1,2,J.Zhao1and J.Guo1

Meshfreemethodshavefoundgoodapplicationsinmanynewresearches,which show very good potential to be powerful numerical tools.As an alternative to the mesh based methods,meshfree methods have the advantage of not using a predefined mesh for the domain discretization.In this study,a mesh free scheme based on the radial point interpolation method was used to solve the topological design of microstructures for composite materials.The explicit form of the radial point interpolation method(RPIM)interpolation augmented with polynomials is presented,which satisfies range-restricted properties and is applicable to integrate a physically meaningful density interpolation.Meanwhile gauss quadrature scheme was applied in order to calculate the physical meaningful properties of microstructure,while the homogenization method is applied to evaluate the effective macroscopic properties of a periodic microstructure.Typical numerical examples are used to demonstrate the effectiveness of the proposed method for designing composite materials with expected effective elasticity tensor.

Topology optimization,meshfree,RPIM,composite materials.

1 Introduction

Topology optimization is a mathematical approach to optimize material layout within a given design space,under a given set of loads and boundary conditions,so that the resulting layout meets a prescribed set of performance targets.Michell(1904)introduce the idea of the topology optimization method firstly by the work of in frame theory.In the past two decades,structural topology optimization as a new approach in structural optimization has developed rapidly with many new contributions to theory,computational methods and applications[Bendsoe and Sigmund(2003)].By now,various schemes have been developed,such as the homogenization method[Bendsoe and Kikuchi(1988)；Guedes and Kikuchi(1990)]the

SIMP method[Zhou and Rozvany(1991)；Mlejnek(1992)；Bendsoe and Sigmund,(1999)],the ESO method[Xie and Steven(1993)],and the level set-based method[Sethian and Wiegman(2000)；Wang,Wang and Guo(2003)；Allaire,Jouve and Toader(2004)；Luo,Wang,Wang and Wei(2008)].

Mostrecently,severalapproachesarestudiedbasedonnodaldesignvariables[RahmatallaandSwan(2004)；MatsuiandTerada(2004)；Guest,PrevostandBelytschko(2004)；Paulino and Le(2009)].According to these approaches,values of material density at nodes are considered as design variables,and traditional element material densities are obtained from nodal values with interpolation schemes.Meanwhile,a number of researchers concentrated on getting rid of the elements and meshes in the process of numerical treatment,and the concepts of meshfree methods have been shaped up[Liu and Gu(2005)].Meshfree methods use a set of nodes scattered on the boundaries of the domain to represent the problem domain and the boundaries.These sets of nodes which are called field nodes do not form a mesh of elements.In this case,contrary to the Finite Element Method(FEM),there is no relation between nodes and background cells.

By now,there are several mesh free methods being developed,such as the S-moothed Particle Hydrodynamics[Liu and Liu(2005)],the Element Free Galerkin Method[Belytschko,Lu and Gu(1994)],Point Interpolation Method[Liu and Gu(2001)],the hp-clouds method[Duarte and Oden(1996)],the Reproducing Kernel Particle Method[Liu,Sun and Zhang(1995)],Generalized Finite Element Method[Babuska,Banerjee and Osborn(2004)]and the Partition of Unity Method[Babuska and Melenk(1997)].Based on the applied function approximation/interpolation scheme,the Moving Least Square(MLS)approximation and Point Interpolation Method(PIM)are the most widely used methods[Liu,Gu and Dai(2004)].The ElementFreeGalerkin(EFG)andMeshlessPetrov-Galerkin(MLPG)[AtluriandZhu(1998)；Atluri and Shen(2002)；Dong,Alotaibi,Mohiuddine and Atluri(2014)]methods have been developed based on MLS approximation.Compared with the MLS approximation,PIM using interpolation to construct shape functions to possess the Kronecker delta function property is useful for implementing the boundary conditions[Cui,Liu and Li(2010)].

Even though meshfree methods have been applied successfully in vast variety of areas and shown very good potential to be powerful numerical tools,the application of methods are still in their developing stage.On the application of topology optimization for material design,the early work was made by Sigmund by applying an inverse homogenization method[Sigmund(1995)；Sigmund and Torquato(1997)].The microstructural composites microstructural composites periodic microstructures,and the layout including the geometric shape and topology of the microstructure is of great importance for achieving the effective properties of the composite.

It is known that inverse homogenization method(IHM)[Sigmund(1994)；Hassani and Hinton(1998)]searches for an appropriate material architecture in the base cell model,which has been employed to create various microstructures for composite with desired properties[Sigmund(1995)；Sigmund and Torquato(1996),Larsen,Sigmund and Bouwstra(1997)；Diaz and Sigmund(2010)；Zhou,Li and Li(2010)；Zhou,Li,Chen,Sun and Li(2011)；Dong,Gamal and Atluri(2013)；Dong and Atluri(2013)]or extreme properties[Sigmund and Torquato(1996)；Sigmund and Torquato(1997)；Sigmund(2000)；Guest and Prévost(2007)；Huang,Zhou,Xie and Li(2013)；Wang,Luo,Zhang and Kang(2014)].

In this study,a mesh free method based on the RPIM is implemented to solve the topology optimization problem in the design of microstructural composite.The numerical homogenization method,to determine the effective properties of the composite,is combined with the RPIM to implement the inverse design of composite materials.The proposed meshfree method will be employed to optimize the shape and topology of the microstructure.Numerical examples will be used to demonstrate the effectiveness of the proposed method for the design of composites.

2 Radial Point Interpolation Method

Since the mesh free methods get rid of using mesh of elements,the field variable u(x)at a point of interest x is interpolated in the problem domain using the function values at neighbor nodes of the point x.The neighbor nodes are called local support domain,and anon-local support domain of a point x determines the number of nodes n to be used to approximate the function value at x which can be expressed as

Where n is the number of neighbor nodes within the non-local support domain of the arbitrary point x,uiis the nodal field variable at field node i and φiis the RPIM shape function corresponding to field node i.The non-local support domain of a point x determines the number of nodes n to be used to approximate the function value at x.

The RPIM interpolation augmented with polynomials can be written as

where Ri(x)is a compactly supported radial basis function(CSRBF),n is the number of CSRBFs,pj(x)is a polynomial in the space coordinate and m is the number of polynomial basis functions,aiand bjare the coefficients to be determined.In radial basis functions,the variable is only the distance between the point of interest x and a node at xi,for 2D problem

There are a variety of CSRBF types and their characteristics have been widely investigated[Wu(1995)；Wendland(1995)].The CSRBFs are strictly positive definite for all r less than or equal to some fix value,and can be constructed to have desired amount of smoothness of 2k.In a CSRBF,there is a shape parameter,δ,that determines the size of non-local support domain of the CSRBF.When r≥ δ ,their value is regarded as zero[Liu and Gu(2005)].In this study,the Wendland-C6 CSRBF[Wendland(1995)]is applied,which is stated as

The polynomial term in Eq.2 is not always necessary.A number of advantages such as improvement in accuracy and stability of interpolations have been found for adding the polynomial terms[Liu and Gu(2005)].Via the studies,it is concluded that the RPIM shape functions with pure CSRBFs usually can not pass the standard patch tests.It means that the pure CSRBFs can not construct a linear polynomial fi eld exactly.Thus,polynomial terms are added to ensure the C1 consistency in order to pass the standard patch test.Meanwhile,it is also found that adding polynomials is bene fi cial to improve the accuracy of the results and the interpolation stability for CSRBFs,and reduce the sensitivity of the shape parameters.Due to these advantages,in this study CSRBFs augmented with linear order polynomial term have been used to construct the shape functions.

In order to obtain the shape functions,coefficients aiand bjin Eq.2 should be determined.For this purpose,a support domain including a number of field nodes is formed for the point of interest x.Coefficients aiand bjin Eq.2 can be evaluated by enforcing Eq.2 to be satisfied at the n nodes surrounding the point of interest x(non-local support domain).This leads to n linear equations,one for each node.The matrix form of Eq.2 can be expressed as

Where the vector of function values Usis

The moment matrix of CSRBFs can be expressed as:

The polynomial moment matrix can be expressed as

There are n+m variables in Eq.9,and the additional m equations can be added using the following m constraint conditions.

Combining Eq.6 and 10 yields the following set of equations in the matrix form:

The matrix R0is symmetric,and the matrixGwill also be symmetric,thus we can obtain

From Eq.2 and 12 we can obtain

Thus the RPIM shape functions can be expressed as

Radial point interpolation shape function ΦT(x)is the main differences of the Mesh free methods and FEM is at the stage of shape function construction.The RPIM shape functions corresponding to the nodal variables vector are finally expressed as

Using CSRBFs as a basis to construct RPIM shape function has several advantages:1)CSRBFs can solve the singularity problem of the RPIM effectively；2)RPIM shape functions are stable,which is flexible for arbitrary and irregular nodal distribution.The RPIM shape functions satisfy six basic properties:1)Delta function property；2)Partitions of unity；3)Compact support；4)Continuity；5)Reproducibility；6)Compatibility.More details and discussion about the properties of the above shape functions which are called RPIM shape functions can be found by[Liu and Gu(2005)].

3 Implementation of RPIM

3.1 Material Interpolation

Considering the RPIM shape function φi,the density ρ(x)at computational point x can be expressed as

Where i=1,...,n represents the total number of field nodes located within the nonlocal inf l uence demain of the computational point,and ρ(x)are the corresponding nodal densities(design variables).

To solve the material design problem using the meshfree method,this paper introduces the concept of material representation model in SIMP[Gibiansky and Sigmund(2000)]into the mesh free method.Thus The tensors of elasticity C at point x can be calculated based on the nodal variable ρ(x)

Where P is the penalty,and the elastic moduli of the solid material phase isCbase.

3.2 Calculation of Effective Material Properties

Homogenization theory is based on the asymptotic expansion of the governing equation,enabling a separation of the macro-and microscopic length scales,so as to extract homogeneous effective material properties from heterogeneous media.As shown in Fig 1,to evaluate the effective properties of the composite using the homogenization method[Guedes and Kikuchi(1990)],we assume that(1)the metamaterials consist of an assembly of microstructures(base cell)；(2)the size of the microstructure is much smaller than the buck composite to allow scaledecomposition；and(3)the effective homogenized property of the composite can be predicted by a single unit cell.

Figure 1:Schematic of microstructural composite.

According to the theory of homogenization,the effective homogenized properties of the composite can be computed as follows:

In this study,the topological shape optimization will be performed within the unit cell Y that is regarded as the design domain.The period Y is assumed to be very small in comparison with the dimension of the overall domain Ω of the medium.

(1)The locally varying strain fieldsare defined by

in which the displacement fields χklcan be obtained by solving the following equation

where ν is the virtual displacement field.

3.3 Formulation of the optimization problem

Theaimofthisworkistooptimizethetopologyandshapeofthemicrostructureunder specified effective elasticity tensors.An optimization problem including these features can be written as:

where wijklis the weighting factor associated with corresponding component of elasticity tensor,ρminand ρmaxare the lower and upper bounds of the design variables to guarantee a stable iteration,Vminand Vmaxare lower and upper bounds to limit the volume fractions of solid phase.The bilinear energy form a(χ ,v,ρ)and the linear load form l(v,ρ)are described by

3.4 Design sensitivity analysis

In this study,the Method of Moving Asymptotes(MMA)[Svanberg(2005)]is applied to solve the optimization problem.Once the nodal design variables are updated by the MMA,the first-order derivatives of the objective function and constraints with respect to the design variables is required.

The sensitivities of the effective elasticity tensor respect to the design variables can be given by

4 Numerical Results and Discussions

In this section,two classic test problems are used to demonstrate the effectiveness of the RPIM meshfree method.The design domain is discretized with a set of meshfree nodes,and the regular background virtual cells are used only in the procedure of numerical integration.In each cell,4×4 Gauss quadrature is used to evaluate the stiffness matrix of the structure.Fig.2 shows the relationship between the background cells,field nodes and Gauss points.The convergence criteria is the change of the design variables of two consecutive iterations is lower than 1e-4,or the maximum number of iterations is 200 based on our numerical experience.

Figure 2:Relationship of background cells,field nodes and Gauss points.

4.1 Material with Desired Elasticity Tensor

In the case, it aims to achieve the effective elasticity tensor0.2,under the constraint of volume fraction V=35%.For the base material,the artificial material properties are given as:Young’s moduli E=10,Poisson’s ratios ν=0.3.

For representing local geometric details of topology plots,the Gauss points are applied to describe structural topology.The contours of the densities on Gauss points are displayed as Fig.3,in which the red colour indicates the solid material(ρ=ρmax=1),while the blue colour shows the weak material/holes(ρ = ρmin=0.01).

Figure 3:Contour plots of point-wise densities on Gauss points.

The results of Case(a)given in Fig.3 from(1)to(6)as the initial design,four intermediate designs,and the optimized design,respectively.It can be seen that the topology optimization is actually an iterative process to re-distribute the material in the design space until the objective is achieved.When selecting the initial design,as Fig.3(1),the initial structure must include holes to provide an inhomogeneous initial design for the calculation using homogenization method,and the location and shapes of those holes may inf l uence the result.Fig.4 shows the discrete plots of the design distribution of optimized microstructure on field nodes.To display the layout of the optimized microstructural composite,the contour for solid material on Gauss points is plotted and the corresponding array are generated as seen in Fig.5.

Figure 4:Topology plots of nodal material densities.

Figure 5:Contour for solid material on Gauss points and the 3×3 array of base cell.

It is also noticed that all the boundaries related to the topological designs are suffi cient smooth for providing curves and distinct material interface,which is benefi cial to manufacturing procedure.From the iteration history in Fig.6,it can be found that the proposed method can effectively find the optimal structures with less than 200 iterations,and the volume constraints are mass conservative.

Figure 6:Iteration histories of the objective function and volume constraint.

4.2 Material with Desired Poisson’s ratio

In this case,a microstructural composite with negative Poisson’s ratios ν=-0.5,under the constraint of volume fractionV=35%.For the base material,the artificial material properties are given as:Young’s moduli E=10,Poisson’s ratios ν=0.3.

In this paper,the materials designed are subject to the plane stress condition,thus,the elasticity tensor can be written in the following matrix form:

Figure 7:Contour plots of point-wise densities on Gauss points.

For convenience,the factor of the matrix is set as an constant d=E/(1-ν2).Positive energy term in the classic theory of elasticity requires that the Poisson’s ratio can be evaluated within the interval[-1,1]for plane problems.The objective effective elasticity tensor is set as

Figure 8:Topology plots of nodal material densities.

Figure 9:Contour for solid material on Gauss points and the 3×3 array of base cell.

The topological changes during the optimization progress are given in Fig.7.The optimized structure turned out with district interface of weak and solid materials,and the boundaries are sufficient smooth with curves.It can be figure out that to achieve the negative value of Poisson’s ratio,the re-entrant structure are generated within the microstructure.Fig.8 shows the discrete plots of the design distribution of optimized microstructure on field nodes.The contour for solid material on Gauss points is plotted and the corresponding array are given in Fig.9.The iteration history in Fig.10,it can be found that the volume constrains are achieved within 60 iterations and topological changes are also finished with 80 iterations.The rest iterations mainly contribute to the shape changes at the detailed structures.The effective elasticity tensor of the optimized design is shown in Fig.10 as well,which illustrate that the Possion’s ratio satisfies the objective ν=-0.5.Similar as the first case,the proposed method can converge with less than 200 iterations,and the volume constraints are mass conservative.

Figure 10:Iteration histories of the objective function and volume constraint.

5 Numerical Results and Discussion

This study presented a topology optimization method for designing composite with desired elasticity properties using meshfree method.The meshfree field nodes are defined to denote the design variables and RPIM is employed to construct a physical meaningful density interpolation between nodal density variables and practical material properties.The method has provided a new possibility for meaningfully combining the CSRBFs based RPIM with the well-established optimization meth-ods.With the interpolation scheme devised in this work,we can create any artificially microstructural composite under periodicity.Our ongoing research is to extend the proposed topological shape optimization method to design problems of metamaterials.

Acknowledgement:The research is partially supported by National Sci-Tech SupportPlan ProjectinChina(2015BAD18B0301),and theHigh School Outstanding Young Teacher Training Program in Guangdong Province(2014Yq2014025).

Allaire,G.；Jouve,F.；Toader,A.-M.(2004):Structural optimization using sensitivity analysis and a level-set method.Journal of Computational Physics,vol.194,pp.363-393.

Atluri,S.N.；Zhu,T.(1998):A new Meshless Local Petrov-Galerkin(MLPG)approach in computational mechanics.Computational Mechanics,vol.22,pp.117-127.

Atluri,S.N.；Shen,S.(2002):ThemeshlesslocalPetrov-Galerkin(MLPG)method:Crest.

Babu?ka,I.；Banerjee,U.D.；Osborn,J.E.(2012):Generalized finite element methods:Main ideas,results,and perspective.International Journal of Computational Methods,vol.1,pp.67-103.

Belytschko,T.；Lu,Y.Y.；Gu,L.(1994):Element-free Galerkin methods.International Journal for Numerical Methods in Engineering,vol.37,pp.229-256.

Bends?e,M.P.；Kikuchi,N.(1988):Generating optimal topologies in structural design using a homogenization method.Computer Methods in Applied Mechanics and Engineering,vol.71,pp.197-224.

Bends?e,M.P.；Sigmund,O.(1999):Material interpolation schemes in topology optimization.Archive of Applied Mechanics,vol.69,pp.635-654.

Bends?e,M.P.；Sigmund,O.(2003):Topology optimization:theory,methods and applications:Springer.

Cui,X.Y.；Feng,S.Z.；Li,G.Y.(2014):A cell-based smoothed radial point interpolation method(CS-RPIM)for heat transfer analysis.Engineering Analysis with Boundary Elements,vol.40,pp.147-153.

Diaz,A.R.；Sigmund,O.(2010):A topology optimization method for design of negative permeability metamaterials.Structural and Multidisciplinary Optimization,vol.41,pp.163-177.

Dong,L.；Gamal,S.H.；Atluri,S.N.(2013):Stochastic macro material prop-erties,through direct stochastic modeling of heterogeneous microstructures with randomness of constituent properties and topologies,by using trefftz computational grains(TCG).CMC:Computers,Materials&Continua,vol.37,no.1,pp.1-21.

Dong,L.；Atluri,S.N.(2013):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,no.2,pp.111-154.

Dong,L.；Alotaibi,A.；Mohiuddine,S.A.；Atluri,S.N.(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,no.1,pp.1-85.

Duarte,C.A.；Oden,J.T.(1996):An h-p adaptive method using clouds.Computer Methods in Applied Mechanics and Engineering,vol.139,pp.237-262.

Gibiansky,L.V.；Sigmund,O.(2000):Multiphase composites with extremal bulk modulus.Journal of the Mechanics and Physics of Solids,vol.48,pp.461-498.

Guedes,J.M.；Kikuchi,N.(1990):Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods.Computer Methods in Applied Mechanics and Engineering,vol.83,pp.143-198.

Guest,J.K.；Prévost,J.H.；Belytschko,T.(2004):Achieving minimum length scale in topology optimization using nodal design variables and projection functions.International Journal for Numerical Methods in Engineering,vol.61,pp.238-254.

Guest,J.K.；Prévost,J.H.(2007):Design of maximum permeability material structures.Computer Methods in Applied Mechanics and Engineering,vol.196,pp.1006-1017.

Hassani,B.；Hinton,E.(1998):A review of homogenization and topology optimization III-topology optimization using optimality criteria.Computers&Structures,vol.69,pp.739-756.

Huang,X.；Zhou,S.；Xie,Y.；Li,Q.(2013):Topologyoptimizationofmicrostructures of cellular materials and composites for macrostructures.Computational Materials Science,vol.67,pp.397-407.

Larsen,U.D.；Sigmund,O.；Bouwstra,S.(1997):Design and fabrication of compliant micromechanisms and structures with negative Poisson’s ratio.Journal of Microelectromechanical Systems,vol.6,pp.99-106.

Liu,G.R.；Gu,Y.T.(2001):A local radial point interpolation method(LRPIM)for free vibration analyses of 2-D solids.Journal of Sound and Vibration,vol.246,pp.29-46.

Liu,G.R.；Gu,Y.T.；Dai,K.Y.(2004):Assessment and applications of point interpolation methods for computational mechanics.International Journal for Numerical Methods in Engineering,vol.59,pp.1373-1397.

Liu,G.R.；Gu,Y.T.(2005):An introduction to meshfree methods and their programming,Springer.

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.

Luo,Z.；Wang,M.Y.；Wang,S.；Wei,P.(2008):A level set-based parameterization method for structural shape and topology optimization.International Journal for Numerical Methods in Engineering,vol.76,pp.1-26.

Michell,A.G.M.(1904):The limits of economy of material in frame-structures.Philosophical Magazine,vol.6 pp.589-597.

Mlejnek,H.P.(1992):Some aspects of the genesis of structures.Structural and Multidisciplinary Optimization,vol.5,pp.64-69.

Paulino,G.H.；Le,C.H.(2009):A modified Q4/Q4 element for topology optimization.Structural and Multidisciplinary Optimization,vol.37,pp.255-264.

Rahmatalla,S.F.；Swan,C.C.(2004):A Q4/Q4 continuum structural topology optimization implementation.Structural and Multidisciplinary Optimization,vol.27,pp.130-135.

Rozvany,G.I.N.；Kirsch,U.(1995):Layout optimization of structures.Applied Mechanics Reviews,vol.48,pp.41.

Sethian,J.A.；Wiegmann,A.(2000):Structural boundary design via level set and immersed interface methods.Journal of Computational Physics,vol.163,pp.489-528.

Sigmund,O.(1994):Materials with prescribed constitutive parameters:an inverse homogenization problem.International Journal of Solids and Structures,vol.31,pp.2313-2329.

Sigmund,O.(1995):Tailoring materials with prescribed elastic properties.Mechanics of Materials,vol.20,pp.351-368.

Sigmund,O.；Torquato,S.(1996):Composites with extremal thermal expansion coefficients.Applied Physics Letters,vol.69,pp.3203-3205.

Sigmund,O.(2000):A new class of extremal composites.Journal of the Mechanics and Physics of Solids,vol.48,pp.397-428.

Sigmund,O.(1997):On the design of compliant mechanisms using topology optimization.Journal of Structural Mechanics,vol.25,pp.493-524.

Svanberg,K.(1987):The method of moving asymptotes-a new method for structural optimization.International Journal for Numerical Methods in Engineering,vol.24,pp.359-373.

Wang,M.Y.；Wang,X.；Guo,D.(2003):A level set method for structural topology optimization.Computer Methods in Applied Mechanics and Engineering,vol.192,pp.227-246.

Wang,Y.Q.；Luo,Z.；Zhang,N.；Kang,Z.(2014):Topological shape optimization of microstructural metamaterials using a level set method.Computational Materials Science,vol.87,pp.178-186.

Wendland,H.(1995):Piecewise polynomial,positive definite and compactly supported radial functions of minimal degree.Advances in computational Mathematics,vol.4,pp.389-396.

Wu,Z.(1995):Compactly supported positive definite radial functions.Advances in computational Mathematics,vol.4,pp.283-292.

Xie,Y.M.；Steven,G.P.(1993):A simple evolutionary procedure for structural optimization.Computers&Structures,vol.49,pp.885-896.

Zhou,M.；Rozvany,G.I.N.(1991):The COC algorithm,Part II:topological,geometrical and generalized shape optimization.Computer Methods in Applied Mechanics and Engineering,vol.89,pp.309-336.

Zhou,S.；Li,W.；Li,Q.(2010):Level-set based topology optimization for electromagnetic dipole antenna design.Journal of Computational Physics,vol.229,pp.6915-6930.

Zhou,S.；Li,W.；Chen,Y.；Sun,G.；Li,Q.(2011):Topology optimization for negative permeability metamaterials using level-set algorithm.Acta Materialia,vol.59,pp.2624-2636.

1College of Engineering,South China Agricultural University,Guangzhou,Guangdong,China.

2Corresponding author.E-mail:enlilv@scau.edu.cn