An efficient homogenization method using the trigonometric interpolation and the fast fourier transform

Vietnam Journal of Mechanics, VAST, Vol. 33, No. 4 (2011), pp. 215 – 223 AN EFFICIENT HOMOGENIZATION METHOD USING THE TRIGONOMETRIC INTERPOLATION AND THE FAST FOURIER TRANSFORM Ngoc-Trung Nguyen1, Christian Licht2 and Jin-Hwe Kweon3 1Kangwon National University, Chunchon, Kangwon-do, South Korea 2LMGC, Université de Montpéllier 2, France 3School of Mechanical and Aerospace Engineering Gyeongsang National University, South Korea Abstract. This study focuses on formulation of the Augmented Lagrangian and ap- plication of the Uzawa’s algorithm to solve the homogenization problem of microscopic periodic media as in composites. Unlike in the finite element model, an equally spaced grid system associated with the microstructure domain is used instead of a finite element mesh topology. Moreover, the trigonometric interpolations for the field variables at every grid point help to handle the periodic conditions. The proposed approach is a compromise between Lagrange multiplier and penalty methods, in that it enables exact representation of constraints while using penalty terms to facilitate the iteration procedure. A typical homogenization problem will be solved using this approach. The results show good con- sistency with those in literatures. Effects of the grid density and the penalty parameter on the convergence have also been investigated. Keywords: Homogenization, Augmented Lagrangian method, trigonometric interpola- tion. 1. INTRODUCTION For composites of complex microstructures, there are two different scales associ- ated with microscopic and macroscopic behaviors to deal with: the slowly varying global variables and the rapidly oscillating local variables. To model a structure of such kind of material using the finite element method (FEM) one should utilize very fine mesh density so that the details at the microscale size can be captured. That leads to a very high com- putational cost and sometimes it is impossible to perform the analysis due to extremely high requirements of computer resources. Instead, a process so-called homogenization is used to characterize the heterogeneous material as a homogenized one and the equivalent material properties are then used in the simulation of the whole structure as in a regular FEM analysis. For the sake of simplicity, the materials can be considered as an assembled body of periodic unit cells as shown in Fig. 1. The assumption of such periodic media has been widely used to develop in both the mathematical analysis and the numerical models [1, 12]. 216 Ngoc-Trung Nguyen, Christian Licht and Jin-Hwe Kweon Among various approaches to predict the effective properties of composites, the mathe- matical homogenization method is preferable due to its systematic background and the ease to implement [1, 10, 11]. One may use a FEM program with a slight modification to solve the homogenization problem to obtain the homogenized coefficients without dif- ficulty [4, 5, 9]. Note that the periodicity conditions must be imposed through constraint equations to reflect the repeatability of the microstructure. Evaluation of several treat- ments of periodicity was discussed in [2]. However, using FEM to solve the homogenization problem may have limitation in meshing the complicated microstructures and handling with a large number of degree of freedom (DOF) in analysis [8]. Fig. 1. The macroscale and microscale of the homogenization problem The objective of this study is to develop an effective approach to solve the homog- enization problem. The variable decomposition technique [3] together with formulation of an Augmented Lagrangian is used leading to a saddle-point problem. The Uzawa’s algo- rithm [6] is then applied to solve the saddle-point problem by iterations. This method is reliable due to the fact that it is proved to be converged [3]. To eliminate the condi- tion on periodicity of the field variables, we employ the trigonometric interpolations. The field variables will then automatically satisfy the periodicity conditions. A grid system associated with the microstructure domain is utilized in this approach instead of a mesh configuration as in the finite element analysis. The meshing and the handling large num- ber of degrees of freedom issues will not be encountered regardless the complexity of the composite microstructure. 2. HOMOGENIZATION OF PERIODIC MEDIA From the asymptotic expansion [1, 11], the homogenized elasticity tensor can be determined explicitly by: ahomijkh = 1 |Y| ∫ Y ( aijkh(y) + aijrs(y)ers(χ kh) ) dy (1) where Y = { y ∈ <n, 0 ≤ yi ≤ Yi, i = 1, n } , |Y| is the volume of the unit cell, An efficient homogenization method using the trigonometric interpolation and the fast fourier transform 217 〈•〉 = 1 |Y| ∫ Y (•) dy denotes the average over the unit cell Y, χkh in (1) is the solution of the cell problem:  − ∂ ∂yj aijrs(y)ers(χ kh ) = ∂ ∂yj aijkh(y) χkh is Y − periodic (2) with the periodicity condition defined by: If y ∈ Y 7→ vi (y) : Y − periodic or vi (y) ∈ Vper (Y) then vi (y) takes equal values on the opposite faces of Y. Generally, we can obtain the solution vE = χ kh by solving the six cell problems and then compute the homogenized elasticity coefficients according to (1). As an alternative, the homogenization problem (2) with periodicity conditions can be formulated in the following forms: For a given macroscopic strain E, (PE) { Find vE ∈ Vper (Y) such that JE (vE) ≤ JE (v) , ∀v∈ Vper (Y)} (3) where JE (v) = 1 2 |Y| ∫ Y a (E+e (v)) (E+e (v))dy Note that the problems (2) and (3) are equivalent. Moreover, the variational formu- lation (3) is equivalent to a problem of minimization with constraints:  Find v ∈ Vper (Y) such that J (v)→ min J (v) = 1 2 a (v, v)− l (v) (4) where a (v, v) = 1 |Y| ∫ Y ae (vE) e (v)dy, l (v) = − 1 |Y| ∫ Y aEe (v)dy Again, by solving 6 problems of formulation (3) or (4) with the imposed macroscopic strains Eij = ( T kh ) ij = 1 2 (δikδjh + δihδjk), where δij is the Kronecker delta symbol; the homogenized coefficients are determined by: ahomijkh = 〈 σkhij 〉 = 1 |Y| ∫ Y σkhij dy = 1 |Y| ∫ Y a ( Tkh + e (vTkh) ) dy (5) In the following sections, we will briefly review the trigonometric interpolation tech- nique and the discrete Fourier transforms. Then, we reformulate the minimization problem (5) by using the Augmented Lagrangian method and the variable decomposition technique. Based on these variational equations, an in-house code will be developed. 3. TRIGONOMETRIC INTERPOLATION AND THE DISCRETE FOURIER TRANSFORMATIONS Consider a unit cell Y = ]0, 1[3, with a given positive number n we define the grid points: (xj, yk, zl) = ( j n , k n , l n ) , where j, k, l ∈ {0, 1, ..., n− 1} as shown in Fig. 2. 218 Ngoc-Trung Nguyen, Christian Licht and Jin-Hwe Kweon Fig. 2. Grid system of the periodic unit cell Due to the periodicity of the field variables, the trigonometric interpolation of a function f (x, y, z) is defined through the values at grid points f (xj , yk, zl) in the system by using the inverse discrete Fourier transform: fˆ (xj, yk, zl) = ∑ r,s,t fˆjkl n3 w−1jkl (xr, ys, zt) = ∑ r,s,t fˆjkl n3 e2ipijxre2ipikyse2ipilzt , r, s, t ∈ {0, 1, ..., n− 1} where, the coefficients fˆjkl of the interpolated function fˆ (x, y, z) are determined by the discrete Fourier transform: fˆjkl = ∑ r,s,t f (xr, ys, zt)wjkl (xr, ys, zt) = ∑ r,s,t f (xr, ys, zt) e −2ipijxre−2ipikyse−2ipilzt The values of the interpolated function fˆ (x, y, z) should be equal to the given values of the function f (x, y, z) at grid points: fˆ (xj , yk, zl) = f (xj, yk, zl) , ∀j, k, l ∈ {0, 1, ..., n− 1} With the above definition of trigonometric interpolation, the field variables in the homogenization problem will then automatically satisfy the periodicity conditions. This benefit helps to eliminate constraints on variables and, as a consequence, to reduce the complexity of the problem. 4. THE AUGMENTED LAGRANGIAN FORMULATION By introducing a supplementary variable q, linked to v through the relation q = ∇v, the problem (4) can be reformulated in the following form: An efficient homogenization method using the trigonometric interpolation and the fast fourier transform 219 For a given macroscopic strain E, find min   1 |Y| ∫ Y 1 2 a (E+ qs) . (E+ qs) dy, ∀v ∈ Vper (Y) and q = ∇v   (6) where qs = 1 2 ( q+ qT ) is the symmetric part of the tensor q = ∇v. The problem (6) is a minimization problem in {q, v} with the supplementary constraint q=∇v. To handle this constraint we use a Largrange multiplier to reduce the problem (6) (and also the problems (2), (3) and (4)) to a saddle-point problem: L (v, q, µ) = 1 |Y| ∫ Y 1 2 a (E + qs) (E + qs) dy+ 1 |Y| ∫ Y µ (∇v− q) dy where µ is a Lagrange multiplier. For r ≥ 0, the Augmented Lagrangian is defined by: Lr (v, q, µ) = 1 |Y| ∫ Y a 2 (E + qs) (E + qs) dy+ 1 |Y| ∫ Y µ (∇v−q)dy+ ∫ Y r 2 |∇v − q|2dy (7) Then, the problem (7) is equivalent to: Find the saddle-point (u,p, λ) such that: Lr (u,p, λ) = min v,q max µ Lr (v, q, µ) (8) Note that if vE is the solution of the problem (3) and (4) then it is equivalent to the fact that {vE, qE = ∇vE, λ} is a saddle-point of Lr. Thus, we have: Lr (vE, qE, µ) ≤ Lr (vE, qE, λ) ≤ Lr (v, q, λ) , ∀ {v, q, µ} Consequently, instead of solving the problem (3) or (4) for vE, we aim to find the saddle-point {vE, qE, λ} of Lr by using the so-called Uzawa’s algorithm [3, 6]. The three- step process at each iteration and the details of the method are given below: Initialization: { q0E, λ 1 } are given arbitrarily. With values of { qn−1E , λ n } calculated at the n-th iteration, vnE, q n E and λ n+1 at the (n+ 1)-th iteration are to be determined successively by: STEP 1: Minimization of Lr in v Taking the variation of Lr in v, the condition of minimization leads to:∫ Y [ (−λn∇ψ) + ( −rqn−1E ∇ψ + r∇v n E∇ψ )] dy = 0, ∀ψ ∈ Vper (Y) (9) By applying the divergence theorem and the periodicity condition of ψ to (9), we yield: ∆vnE = div ( qn−1E + λn r ) STEP 2: Minimization of Lr in q Taking the variation of Lr in q, the condition of minimization leads to: a (qnE)s + rq n E = −aE + r∇v n E − λ n (10) 220 Ngoc-Trung Nguyen, Christian Licht and Jin-Hwe Kweon By solving (10) we can obtain the symmetric and the anti-symmetric parts of qE. The updated values of qnE are: q n E = (q n E)s + (q n E)a STEP 3: Updating λn λn+1 = λn + r (∇vnE − q n E) (11) The solutions of (8) are obtained by iterations. For each iteration, there are three steps, in which (9), (10) and (11) are readily to be processed by using the above-mentioned trigonometric interpolations and the Fourier transforms of displacement field variables. To avoid obtaining a local minimizer solution, the convergence criteria is employed such that the iterating process terminates when the relative error of a variable is smaller than a chosen value ε for a NCONV number of successive iteration. The relative error is taken in L2−norm, i.e. ∥∥∥uk+1i − uki ∥∥∥ L2 /∥∥uki ∥∥L2 ≤ ε. The variable value ui can be the displacement at grid points v, divergence of the displacement q or the Lagrange multiplier λ. Once convergence has been reached, qE coincides with ∇vE. The Uzawa’s algorithm with a suitable choice of the initial and the multiplier values has been proved to achieve better convergence. 5. NUMERICAL EXAMPLE To verify the proposed method, we investigate a typical example of a unidirectional fiber-reinforced composite material using 2D model to compare with results in [7]. Fig. 3. Grid system of the unit cell A grid system of the unit cell for the 2D model is used to predict the properties of the glans fiber/epoxy reinforced composite system as shown in Fig. 3. The volume fraction of fibers is 45%. The plane strain and perfect bonding between constituents assumptions are made. The data on Young’s modulus and Poisson’s ratio of glass fiber and epoxy resin are given in Table 1. Note that for a 2D problem, we use an exchange in notation between the tensor aijkh and its matrix form Cij . Explicitly, in this problem we denote: C11 = a1111, C22 = a2222, An efficient homogenization method using the trigonometric interpolation and the fast fourier transform 221 Table 1. Material properties of glass fiber and epoxy Material E (GPa) ν E-glass fibers 70 0.2 Epoxy resin 3.5 0.35 C12 = a1122 and C66 = a1212. To obtain the homogenized modulus, we will solve 3 cell problems by successively impose the three macroscopic strains: E = ( T11 ) = [ 1 0 0 0 ] ,E = ( T22 ) = [ 0 0 0 1 ] and E = ( 2T12 ) = [ 0 1 1 0 ] Fig. 4. Comparison between the results obtained by the present study and those of three FEM models in [7]: (a) C11 = C22 coefficients, (b) C12 = C21 coefficients, (c) C66 coefficient Calculations of the homogenized values and comparisons with [7] are shown in Fig. 4. A quite good agreement can be observed. Note that [7] used the finite element method to solve the problem with different geometries and sizes of the unit cell. The obtained results of the present study also closely match with those in [9]. With a coarser grid density, the solutions show a small discrepancy with the referred values. However, for finer grid density the solutions are stable and almost identical to those of [7]. The sensitiveness of the solution on the grid density is thus negligible after a certain level. Hence, to employ 222 Ngoc-Trung Nguyen, Christian Licht and Jin-Hwe Kweon the developed program efficiently we can use a rather coarse grid system to obtain results with the same accuracy as in case of fine grid resolution. Fig. 5. Effects of the penalty parameter r on the convergence The influence of the parameter r on the convergence is shown in Fig. 5. We may find that the penalty parameter needs not to tend to infinity to obtain the exact solution as in the ordinary penalization methods, it can be small. That is one of the advantages of the Augmented Lagrangian method, appearance of the last term in (7) improves the convergence properties of the algorithm. 6. CONCLUSION The Augmented Lagrangian approach associated with the trigonometric interpo- lation in this study shows its effectiveness in solving the homogenization problems. The trigonometric interpolation helps to remove the constraint on periodicity. Furthermore, the Augmented Lagrangian method is a compromise between Lagrange multiplier and penalty methods, in that it enables exact representation of constraints while using penalty terms to facilitate the iterative procedure. Combination of these two advantages leads to a quite promising and potential way for handling the homogenization problems, espe- cially for predicting properties of very complicated microstructured materials. Besides, the concept in this approach makes it ready to be used with a microscope image of the real microstructures. The meshing and the handling large number of DOF issues will not be encountered because the grid system is used instead of a topology of the finite element nodes. The applicability of the algorithm discussed in this study is broad, ranging from linear elasticity to nonlinear and/or plasticity applications. An efficient homogenization method using the trigonometric interpolation and the fast fourier transform 223 REFERENCES [1] A. Bensoussan, J. L. Lions, G. Papanicolaou, Asymptotic analysis for periodic structures, North-Holland Publishing Company, (1978). [2] O. Debordes, C. Licht, J. J. Marigo, P. Mialon, J. C. Michel, P. Suquet, Charges limites de milieux fortement heterogenes, Note Technique No. 85-3, Laboratoire de Mecanique Generale des Milieux Continus, Universite des Sciences et Techniques du Languedoc, (1985). [3] M. Fortin, R. Glowinski (Eds.), Augmented Lagrangian Methods: Applications to the Numer- ical Solution of Boundary-Value Problems, North-Holland Publishing Company, (1982). [4] J. M. Guedes, N. Kikuchi, Preprocessing and postprocessing for materials based on the ho- mogenization method with adaptive finite element methods, Computer Methods in Applied Mechanics and Engineering, 83(2) (1990), 143 - 198. [5] B. Hassani, A direct method to derive the boundary conditions of the homogenization equation for symmetric cells, Communications in Numerical Methods in Engineering, 12(3) (1996), 185 - 196. [6] C. Licht, Frottement, Viscoplasticité et Homogénéisation, PhD Thesis, University Montpéllier II, (1987). [7] D. Lukkassen, L. E. Persson, P. Wall, Some engineering and mathematical aspects on the homegenization method, Composites Engineering, 5(5) (1995). [8] H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlin- ear composites with complex microstructure, Computer Methods in Applied Mechanics and Engineering, 157(1-2) (1998), 69 - 94. [9] N. T. Nguyen, K. H. Nguyen, Y. B. Park, J. H. Kweon, Characterization of composite materials with complex microstructure, Proceeding of the KSAS conference, Gyeong-ju, South Korea, (2009). [10] O. A. Oleinik, A. S. Shamaev, G. A. Yosifian, Mathematical problems in elasticity and ho- mogenization, North-Holland, (1992). [11] E. Sanchez-Palencia, A. Zaoui (Eds.), Homogenization Techniques for Composite Media, Springer, (1987). [12] L. Tartar, Quelques remarques sur l’homogénéisation, in: Fujita (Ed.), Functional analysis and numerical analysis, Japan Society for the Promotion of Science, Tokyo, (1978). Received June 20, 2011