ac

1348 Louvain-la-Neuve, Belgium Received 12 February 2003; accepted 12 May 2003 periodic linear media, this theory can be rigorously justified, through convergence results, see e.g. [3,22]. Thus, as the ratio of the period length to the structure length goes to zero, it can be shown that the actual solution of the heterogeneous structure tends to the solution of a problem written for a structure with an Comput. Methods Appl. Mech. Engrg. 192 (2003) 3163–3177 * Corresponding author. Address: Laboratoire de M�eecanique et Mat�eeriaux, Ecole Centrale de Nantes, 1 Rue de la Noe, 44321 1. Introduction The extended finite element method (X-FEM) simplifies greatly the analysis of structures with complex geometry since the mesh is not required to match this geometry. This property is particularly useful when homogenizing the behavior of microstructures, complex by nature. We focus here our attention to the application of X-FEM to microscale problems involved in the homogenization theory. The homogenization theory is a powerful mean for analyzing and modeling heterogeneous structures. Moreover, in the case of Abstract In multiscale analysis of components, there is usually a need to solve microstructures with complex geometries. In this paper, we use the extended finite element method (X-FEM) to solve scales involving complex geometries. The X- FEM allows one to use meshes not necessarily matching the physical surface of the problem while retaining the ac- curacy of the classical finite element approach. For material interfaces, this is achieved by introducing a new enrichment strategy. Although the mesh does not need to conform to the physical surfaces, it needs to be fine enough to capture the geometry of these surfaces. A simple algorithm is described to adaptively refine the mesh to meet this geometrical requirement. Numerical experiments on the periodic homogenization of two-phase complex cells demonstrate the ac- curacy and simplicity of the X-FEM. � 2003 Elsevier B.V. All rights reserved. Keywords: Homogenization; Surface of discontinuity; X-FEM; Periodicity A computational approach to handle complex microstructure geometries N. Mo€ees a,*, M. Cloirec a, P. Cartraud a, J.-F. Remacle b a Laboratoire de M�eecanique et Mat�eeriaux, Ecole Centrale de Nantes, 1 Rue de la Noe, 44321 Nantes, France b D�eepartement de G�eenie Civil et Environnement, Universit�ee Catholique de Louvain, Ba^atiment Vinci, Place du Levant 1, www.elsevier.com/locate/cma Nantes, France. Tel.: +33-2-40-37-68-22; fax: +33-2-40-37-25-73. E-mail address: nicolas.moes@ec-nantes.fr (N. Mo€ees). 0045-7825/03/$ - see front matter � 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0045-7825(03)00346-3 homogenized behavior. Furthermore, it turns out that the effective properties are obtained from the so- lution of a boundary value problem to be solved on a period of the structure, which will be called the basic cell problem in the sequel. Since the basic cell problem is a boundary value problem, classical numerical methods can be used for computing their solution, see [4] for a state-of-the-art in this matter. The majority of studies have been based on the finite element method (FEM), see [8,12] and references herein, and only few papers present a boundary element method implementation, see e.g. [10,18]. Let us also mention the fast Fourier transforms based numerical method [17]. This method is well suited for periodic homogenization, and is very efficient compared to FEM, especially in the linear range [12], pro- vided that no voids or rigid inclusions are present in the basic cell. The advantages of the FEM in micromechanical analysis are indeed the same as in standard engineering problems: its flexibility, and its applicability to non-linear problems, anisotropic materials, and arbitrary geometries. Following [4], one may distinguish, in FEM based periodic homogenization studies, four main groups. In most of the studies, the mesh models the heterogeneities boundaries in the unit cell, i.e. only one phase is present within each finite element. The main drawback of such an approach is the explicit microstructure modeling, which leads to problems for generating the mesh for complex geometries, and requires the use of sophisticated tools, see e.g. [26]. One way to overcome these difficulties is to use a digital image based FEM technique, as initially proposed in [9]. It consists in using a uniform mesh that has the same resolution as the digital image, and then to identify each pixel or voxel as a finite element. Such an approach, however, leads to models which are computationally expensive. A third approach also uses regular meshes, but the interfaces between constituents may be represented independently of element boundaries. For an element cut by an interface, the heterogeneity is treated at the integration point level [27]. As we shall see in the numerical experiments, this strategy yields a reasonable rate of convergence for the homogenized parameters of the basic cell but suffers a slow rate of convergence for the quality of the overall stress distribution over the basic cell (energy norm in the stresses). Another technique is the Voronoi cell finite element model [7], in which an element incorporates the effect of an embedded inclusion in the matrix. Both stress and displacement fields need to be discretized with this technique. The X-FEM provides yet another way to solve the basic cell problem. The advantage of this approach is to retain all the advantages of the finite element approach (applicability to non-linear and anisotropic constitutive laws, wide range of codes already written, robustness, . . .) while considerably easing the meshing step. Indeed, the physical surfaces of the problem do not need to be meshed. They are taken into account by enriching the finite element approximation space through the partition of unity technique [11]. The adequate choice of the enrichment function for cracks was discussed in [1,14], for intersecting cracks and voids in [6]. Concerning material interfaces, an enrichment was introduced in [24]. This latter paper was also the first one to use a level set representation of surfaces in the X-FEM context. This representation not only reduces the surface storage to a usual finite element field but also provides a natural way to express the enrichment. Note that in order to build the level set representation of the material interfaces, the CAD geometries of these surfaces are not necessarily needed. All that is needed is a function returning for one ðx; y; zÞ point the distance to the closest material interface. In other words, the X-FEM approach coupled to the level set representation of the surfaces not only simplifies the mesh step but also the geometrical data pre-processing. In this paper, we improve the enrichment functions for material interfaces and obtain a convergence rate very close to the one obtained with regular finite elements (i.e. conforming meshes). The remainder of this paper is organized as follows. The basic cell problem is described in Section 2. The X-FEM approach to solve the basic cell problem is detailed in Section 3 and a new enrichment function for 3164 N. Mo€ees et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3163–3177 material interface is introduced and compared to previous approaches. The fourth section uses the X-FEM For a linear elastic medium, the complete determination of a requires the solution of problem (3) for six independent data E. One way to proceed is to consider successively the elementary macroscopic strain N. Mo€ees et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3163–3177 3165 states corresponding to one component of E equal to unity, the others being zero. Thus, introducing the second order tensor I ij defined by ij to homogenize three different basic cells: a cell with orthogonal non-woven fibers, a cell with woven fiber and, finally, a cell with randomly placed spherical inclusions. 2. The basic cell problem The formulation of the basic cell problem for composites with periodic microstructure can be derived in a systematic way using the two-scale asymptotic expansion method [3,22], or following a process which is also valid for random media [12,25]. The latter approach is used in this paper. The microstructure description is achieved through the definition of a representative volume element, containing all the geometrical and material data. In the case of a periodic medium, the representative volume element is simply a basic cell X, which forms the composite by spatial repetition. At the macroscopic scale, stress and strain fields are denoted by R and E, which are the average of the corresponding microscopic field, r and e, on the basic cell (the macroscopic and microscopic scales are denoted by x and y, respectively): RðxÞ ¼ hrðx; yÞi ¼ 1jXj R X rðx; yÞdX; EðxÞ ¼ heðx; yÞi: ( ð1Þ The homogenized behavior relating R to E is built using a localization process which consists in submitting the composite structure to a given macroscopic stress or strain state. In the case of a periodic micro- structure, such a problem can be studied on a single basic cell, thanks to periodic boundary conditions. Moreover, the local periodicity of the strain eðx; yÞ with respect to the microscopic scale y, and the relation EðxÞ ¼ heðuðx; yÞÞi leads to the following decomposition of the displacement field: uðx; yÞ ¼ EðxÞ � yþ umðyÞ; um periodic on oX: ð2Þ Therefore, with the data E, the basic cell (or localization) problem is defined by Find r; e; um such that div rðyÞ ¼ 0 rðyÞ ¼ aðyÞ : eðyÞ ¼ aðyÞ : E þ eðumðyÞÞð Þ um periodic r � n anti-periodic; 8>>>< >>>: ð3Þ where ‘‘periodic’’ (‘‘anti-periodic’’) means that values on opposite sides of the boundary oX are equal (opposite). The problem (3) admits a unique solution, up to a rigid body translation. Once problem (3) is solved, the strain localization tensor D is obtained eðyÞ ¼ DðyÞ : E ð4Þ and the effective stiffness ahom may be computed as ahom ¼ haðyÞ : DðyÞi: ð5Þ hom I kh ¼ 12ðdikdjh þ dihdjkÞ ð6Þ Since in this paper a displacement approach is used, the right-hand side of (8) is known. Therefore, the only specific feature of problem (8) is the periodicity condition [12]. 3. X-FEM discretization In this paper, we consider the homogenization of two-phase basic cells. The material interface separating the two phases is described by a level set function. The level set representation is then used to express a jump in the strain field within the elements crossed by the material interface. 3.1. Level set representation Consider a two-phase basic cell. The two phases are separated by an interface which we may locate either explicitely (set of CAD entities) or implicitly by assigning to each node I of the mesh, the distance /I to the interface (with a positive sign if node I is located in one phase and a negative sign if node I is located in the other phase). Next, we may interpolate these nodal informations using the finite element shape functions, NIðxÞ, yielding the level set expression: /ðxÞ ¼ X I /INIðxÞ; ð11Þ where the sum over I indicates a sum over the nodes of the mesh (more precisely only the nodes belonging to the element containing point x). The iso-zero of the level set function / approximates the true location of the interface. A strategy for constructing the level set function in the case where the interface is known only by a set of points is described in [2]. With the X-FEM, the mesh does not need to conform to the interface but the mesh has to be fine enough to locate precisely the interface and to accurately resolve the displacement field i.e. two kinds of error and denoting umij the solution of problem (3) with the data I ij, one has Dkhij ¼ eðumijðyÞÞð Þkh: ð7Þ The problem (3) corresponds to the strain approach, the data being the macroscopic strain E. One can also use a stress approach, with the data R, thus leading to effective compliance Ahom, which is equal to the inverse of ahom. It is important to mention that, in the framework of periodic homogenization, non-periodic materials (i.e. with some randomness in the geometrical shape of the heterogeneities and their spatial distribution) may be studied. In that case, however, one has to use a sufficiently large basic cell, in order to obtain a statistically representative response and accurate effective properties. In this paper, such examples of basic cells will be treated. The variational formulation of the basic cell problems (3), for a data E ¼ I ij is given by Find umij periodic such that aðumij ; vÞ ¼ �aðI ij; vÞ 8v periodic; � ð8Þ where we use the bilinear form: aðu; vÞ ¼ 1jXj Z X eðuÞ : a : eðvÞdX: ð9Þ The effective properties are determined from ahomijkl ¼ aðI kl þ umkl ; I ij þ umijÞ ¼ aðI kl; I ijÞ � aðumkl ; umijÞ: ð10Þ 3166 N. Mo€ees et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3163–3177 appear: geometrical and numerical errors. Since a level set interpolated on the finite element mesh is used to locate the interface, one can notice that both geometrical and numerical errors are reduced as the mesh gets finer (as with the FEM). As an example, Fig. 1 shows the geometry of the basic cell for a landing helicopter grid located on the deck of an aircraft carrier. Using the coarse mesh, Fig. 2 (left), the iso-zero level set obtained is shown Fig. 2 (right). It can be seen that the geometrical representation is poor. A finer mesh, Fig. 3 (left), gives a much better representation of the surface location, Fig. 3 (right). Fig. 1. Unit cell geometry for an helicopter landing grid. Fig. 2. A coarse mesh and the corresponding iso-zero level set. N. Mo€ees et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3163–3177 3167 Fig. 3. A refined mesh and the corresponding iso-zero level set. The finer mesh was obtained using a simple adaptive mesh refinement strategy. The idea is to reduce the size of the elements located close to the physical surface if the local radius of curvature of the surface is smaller than the characteristic element size. The radius of curvature q on an element is computed by the second gradient of the level set interpolation over the element [15,23] and the characteristic element size lc is computed by lc ¼ ðV � d!Þ1=d where V is the element volume and d the problem dimension. An element is tagged for refinement if lc > aq and j/j < bltot; ð12Þ where / is the level set value at the centroid of the element and ltot the side length of the basic cell. The dimensionless parameters a and b govern the extent and level of refinement. Based on the elements tagged for refinement, the mesh can be refined where needed. In order to meet the condition lc > aq starting from one initial coarse mesh, several intermediary meshes are usually needed. As meshes become finer, the evaluation of the curvature q becomes more accurate. The process ends when no element is tagged. The mesh shown Fig. 3 was obtained with a ¼ 0:1 and b ¼ 0:03. The gmsh meshing tool [19] and the aomd mesh database library [20,21] were used to refine the mesh based on the elements tagged. 3.2. Enrichment strategy Once the level set function is defined, we can express the X-FEM approximation. The classical finite approximation over the basic cell is given by uFEM ¼ X I uINIðxÞ: ð13Þ If the finite element mesh conforms to the material interface, the (linear) approximation above yields an OðhÞ convergence rate in the energy norm (provided the solution is smooth) [8]. On the contrary, as will be shown in the numerical experiments, if the finite element does not conform to the interface, a poor rate of convergence is obtained. For instance, a two-phase 1D problem with a mesh non-conforming to the in- terfaces gives an asymptotic rate of convergence as poor as Oð ffiffiffihp Þ [16]. In order to avoid poor rates of convergence, the X-FEM adds an enrichment to the classical finite element approximation through a partition of unity technique [11]: uX-FEM ¼ X I uINIðxÞ þ uenr; uenr ¼ X J aJNJðxÞFðxÞ: ð14Þ The additional degrees of freedom aJ are added at the nodes for which the support is split by the interface. The function F is called the enrichment function. It is discontinuous in its derivative across the interface. In [24], the enrichment is defined as the absolute value of the level set function which indeed has a discon- tinuous first derivative on the interface: F 1ðxÞ ¼ X I /INIðxÞ ����� �����: ð15Þ A similar choice [2], is to define uenr as uenr ¼ X J aJNJðxÞðF1ðxÞ � F1ðxJÞÞ: ð16Þ A closer look at the two enrichments above indicates that they will yield the same approximation space. However, the conditioning number of the two resulting matrices will be different in general. It was shown in 3168 N. Mo€ees et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3163–3177 [2,24] that a smooting of F 1 away from the element layer containing the interface somewhat improved the convergence. The effect of this smoothing in 1D is illustrated Fig. 4. For 2D problems, the smoothing strategy is detailed in [24]. In this paper, we suggest another choice for the enrichment function: F 2ðxÞ ¼ X I j/IjNIðxÞ � X I /INIðxÞ ����� �����: ð17Þ This enrichment is shown Fig. 4 in the 1D case. For two- and three-dimensional problems, the enrichment function F 2 is a ridge centered on the interface and has zero value on the elements which are not crossed by the interface. Remark: If the mesh conforms to the material interface, no nodes are enriched since no element is cut by the material interface. The X-FEM then does behave as the FEM. On an element cut by the interface, the integration needs to be performed carefully. First, we divide the element into subdomains matching the material interface as described in [14]. Then, the number of integration points over each subdomain is chosen so that the integration is ‘‘exact’’ (note that F ðxÞNJ ðxÞ is polynomial). On the elements, not cut by the interface but neighbouring an element cut (for instance element A, Fig. 4), no subdivision needs to be performed but the order of the integration needs to be chosen with care since the enrichment acts on the element. 3.3. Convergence tests In order to compare the convergence rate of the classical and X-FEM, we consider a 2D circular and a 3D spherical inclusion under tension in an infinite domain. These benchmark problems were already considered in [24] and [2] in 2D and 3D, respectively. In order to compare the numerical and exact solution, we use the energy norm. Let uh be an approximate solution and uex the exact solution, the ‘‘exact’’ error is defined by ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiRq F + smoothing F 1 1 2F Interface element A Nodes Enriched Nodes Fig. 4. Several choices for the enrichment function. N. Mo€ees et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3163–3177 3169 e ¼ X eðuh � uexÞ : a : eðuh � uexÞdXffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiR X eðuexÞ : a : eðuexÞdX q : ð18Þ Concerning the material properties, the Poisson ratio is set to 0.25 for the inclusion and to 0.3 for the matrix. The Young�s modulus ratio between the inclusion and matrix is set to 10. Figs. 5 and 6 show the decrease in the energy error with respect to the exact solution in four different cases: • A finite element mesh conforming to the material interface (denoted FEM). • A finite element computation with a non-conforming mesh (denoted FEM-non-conforming). On the ele- ments containing the material interface, the integration of the bilinear form is performed using the ap- propriate constitutive law at each integration point. 3170 N. Mo€ees et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 3163–3177 • The X-FEM using the enrichment F 1 improved by a smoothing (denoted X-FEM-1+ smoothing). • The X-FEM using the new enrichment, F 2, given Eq. (17) (deno