An approach to approximate and fem-Model the conductivity and elasticity of multicomponent material

available numerical or experimental reference conductivity data for particular composites. We use the eXtended Finite Elements Method (XFEM) to estimate the effective conductivity of 2D macroscopically-isotropic composites containing elliptic inclusions and the equivalent ones with circular inclusions for comparisons with the approximations. 3. Scope The thesis focuses on conductivity and elasticity of multi- component materials, the Finite Element Method (FEM) and approximation schems 4 4. Research methods  Near interaction approximations has been constructed from the minimum energy for the macroscopic conductivity and elasticity of the multi-component matrix composites with spherical (circular) inclusions. Equivalent replacement of complex-geometry inclusions by the equivalent spherical, circular, disk and needle ones with equivalent properties using polarization approximation, dilute solutions, and experimental referemce results.  Numerical method: use Matlab program to homogenize some periodic material models in the framework of FEM method (XFEM). The results of FEM are considered as the accurate reference results for comparisons with the approximation ones. 5. The contributions of the thesis Beside Introduction section, the thesis contains 3 Chapters, a Conclusion section and a list of publications relevant to the thesis. References cited in the thesis are listed at the end of the thesis. CHAPTER 1. OVERVIEW 1.1. Opening Multi-component materials have complex structures, different individual mechanical properties. Many authors offered different evaluation methods, including the effective medium approximations and the variational ones. Geometric parameters have bên added to improve the étimates. In this chapter, the author presents the concept of hômgenization and an overview of the constructions of approximation methods for complex multi-component materials. The stress field ( ) x is related to the strain field ( ) x by Hook’s law: ( ) ( ) : ( ), x C x x (1.1) The average values of the stress and strain on V is defined as: 1 1 , . V V d d V V      x x (1.2) 5 Assume homogeneous boundary conditions for displacements: ( ) . 0u x x (1.3) Or the respective ones for the tractions    0n n (1.4) With the solutions σ, ε on V, the relationship between the averaged stress and strain on V is presented through the effective elastic tensor Ceff: : , ( , ).eff effk   eff effC C T (1.5) effk and eff are effective elastic bulk and shear moduli. Another approach is to determine the effective elasticity coefficients by finding the infimum of the energy function on V (the fields  need to be compatible): 0 0 0: : inf : : ,eff V d          C C x (1.6) Or through the dual principle (the fields  need to be equilibrated): 0 0 1 0 1: ( ) : inf : ( ) : .eff V d           C C x (1.7) Similarly, the equations for the conductivity problem: The flux J must satisfy the equilibrium condition: · ( ) 0 J x With the solutions J,  E T on V, the thermal conductivity coefficient (effective) ceff is determined as: .eff effJ c E c T     (1.8) The minimum energy principles are also the main tools to find the macroscopic conductivity: 0 0 0· inf · ,eff V c d    E EE E cE E x (1.9) and: 6 1( ) · inf · ,eff V c d     0 0 0 1 J J J J c J J x (1.10) 1.2. Overview of approximation methods for multi-component materials 1.2.1. Dilute solutions The effective conductivity ceff of the dilute solution of ellipsoidal inclusions with axes ratio a: b: c, randomly oriented in a continuous matrix is expressed in the form: ( ) ( , ) , 1 ,eff M I I M c I M Ic c v c c D c c v    (1.11) 1 1 1 ( , ) , 3 (1 ) (1 ) (1 ) [ ]Mc I M I M I M I M c D c c c A c A c B c B c C c C          The general fomula of ( , )c I MD c c for spherical (d=3) and circular (d=2) inclusions is: ( , ) . ( 1) M c I M I M dc D c c c d c    1.2.2. Maxwell Approximation Maxwell approximation is built for 2-phase material from the matrix + spherical inclusions with any volume ratios, not limited by dilute distribution case (M - matrix symbol, I – inclusion symbol).       1 1 * * * * 2 1 * * * * 1 , 1 2( 1) ; , 1 2( 1)( 2) ( ) ; . 2 4 eff I M M I M M eff I M M M M I M M M eff I M M M MA M M M I M M M M M v v c d c c d c dc v v d K K K K d K K K d v v d K d d dK d                                                 (1.12) 1.2.3. Differential Approximation - DA we obtain the following differential equations for the effective conductivity ceff = c(1) of the composite 7 1 1 1 ( ) ( , ), 1 (0) , 0 1 , , n I I c I I n M I I dc v c c D c c dt v t c c t v v                   (1.16a) For elastic coefficient 1 1 1 1 ( ) ( , , , ), 1 1 ( ) ( , , , ), 1 (0) , (0) 0 1 , , n I I K I I I n I I K I I I n M M I I dK v K K D K K dt v t d v D K K dt v t K K t v v                                         (1.16b) 1.2.4. Self-consistent approximations - SA The Self-consistent approximation method (SA) for composite materials n components, is cSA=c solution of the following equation: 1 ( ) ( , ) 0 . n I I c Iv c c D c c        (1.17) SA for the moduli of elasticity are the solutions KSA=K and  SA=  of a system of two equations 1 1 ( ) ( , , , ) 0, ( ) ( , , , ) 0 . n I I K I I n I I M I I v K K D K K v D K K                           (1.18) 1.2.5. Mori-Tanaka Approximation (MTA) Mori-Tanaka type approximation (MTA), based on the assumption that the fields in an circular inclusion are determined as if it is embedded in the matrix with remote average strain of the matrix, yields [Le Quang] for two-phase composites 1 1( )·{ [ · ·( )] } .MTA M I I M M M I M Iv v v       c c c c I p c c c I (1.19) 8 While for the multi-phase ones (matrix + n inclusions) 1 1 1 1 1 1 1 { ·[ · ·( )] } ·{ [ · ·( )] } . n MTA M M I I M I M n M I M I M v v v v                          c c c I p c c c I I p c c c (1.20) MTA for the effective conductivity of d-dimentional multi- component isotropic materials with spherical inclusions (circular) has the following form 1 1 ( ) / [ ( 1) ] . / [ ( 1) ] n I I M M I M MTA M n M I M I M v c c dc c d c c c v v dc c d c                    (1.21) 1.3. Three-point correlation estimates of Phạm ĐC Three-point correlation estimates of Phạm ĐC are for the effective thermal conductivity of the multi-component materials. The bounds have been built from the minimum energy principles. The general expression of the upper bound for ceff is eff 0 **(2 ) ,cc P c c  (1.22) where c0 is a positive parameter, 1 ** ** * ( ) ,c v P c c c c             (1.23) 2 ** 0 1 , , 10 3( ) ( ) . 2 n nv c c c A X X c c                  (1.24) Similarly, the respective expression of the dual principle is written in the form: 2 2 1 1 ** 0 0 0 1 , , 10 3 (1 2 ) ( ) . 2 n nv c c c c c A X X c c                     (1.25) 9 Where 1 0 0 1 0 0 1 2 2 1 2 2 n n v X c c c c v X c c c c                       (1.26) We choose the value c0 to eliminate the component ** **,c c to make the inequality stronger to get the respective bounds. 1.3. Finite element method for homogenization solution The finite element method (FEM), is a numerical method for solving problems of engineering and mathematical physics. Typical problem areas of interest include structural analysis, heat transfer, fluid flow, mass transport, and electromagnetic potential. The analytical solution of these problems generally require the solution to boundary value problems for partial differential equations. The finite element method formulation of the problem results in a system of algebraic equations. The method approximates the unknown function over the domain. To solve the problem, it subdivides a large system into smaller, simpler parts that are called finite elements. The simple equations that model these finite elements are then assembled into a larger system of equations that models the entire problem. FEM then uses variational methods from the calculus of variations to approximate a solution by minimizing an associated error function. A typical work out of the method involves dividing the domain of the problem into a collection of subdomains, with each subdomain represented by a set of element equations to the original problem, followed by systematically recombining all sets of element equations into a global system of equations for the final calculation. The global system of equations has known solution techniques, and can be calculated from the initial values of the original problem to obtain a numerical answer. 10 CHƯƠNG 2. FINITE ELEMENT METHOD 2.1. Introduce 2.1.1. FEM for thermal solution The finite element method with fine meshes is now used for reference comparisons. Let us consider a periodic cell Ω with the external boundary ∂Ω. The strong form of the conductivity problem is written as:                 0 , , , periodicin , . antiperiodicin , in c in T in T          q x q x x E x E x x x q x n (2.1) The weak form associated with the above equations is given by: .i i ic N N T d     (2.2) Using linear form functions for triangular elements with 3 joint as:  N x, y ax by c,   (2.3) Equation in matrix form: ( ) .e e eE x B T           (2.4) The effective thermal conductivity in x1 direction can be computed as: 1 1 1 . eff V L c q d V T     The domain V is the periodic cell, so the effective thermal conductivity of RVE is given by E 1 eff eff F Mc c 11 We solve the temperature at each node element position, and find the effective conductivity coefficient according to the heat transfer equation: . 1 avg eff Tq k X    (2.5) 2.1.2. Model of thermal solution We consider the three-dimensional cubic periodic microstructure as given in Fig.2.1. The finite element mesh for a periodic cell is given in Fig. 2.2. Fig 2.1: Mesing up for Body center cubic Some results about thermal solution: Fig 2.2: Thermal distribution in the model 12 2.1.2. Model for elastic solution The displacement field according to the degrees of freedom at the joint element qe      . . e e u N q (2.6) The deformation state of the joint elements will be:             . e e e e u N q B q      (2.7) The equation for stress of elements:      . e e D  (2.8) The full potential of the element:    .e eeu U A  (2.9)             1 . 2 T T e e e ee e e q q K q q P  (2.10) As a result, we get the equation to solve by finite element method:                 1 1 0 . e eN N T T e e ee e e e L K L q L P q              (2.11) Or    .K q P    (2.12) 2.2. Extended – Finite Element Method (XFEM) 2.2.1. LevelSet function The XFEM displacement approximation can be expressed by *( ) ( ) ( ) ( ) . e h i i j j i i u x N x u N x x a      (2.13) One important example of such a function would be the signed distance function: 13 ( ) ( ) ,x s x x x   (2.14) Fig 2.3: LevelSet function An example of level set function for one and four inclusions is provided in Figure 3: Fig 2.4 Level set function defining a circular interface 2.2.2. Using eXtend – Finite element method for thermal solution The field equations of thermo-static problem are given by ( ) ( ) 0 in ( ) ( ) ( ) q x r x q x C x T x        (2.15) Above, q(x) denotes the heat flux, r(x) is a heat source term and c(x) the conductivity tensor. The periodic boundary conditions are: q.n is antiperiodic on  , and T is periodic on  . (C) is the conductivity matrix. More precisely, C(x) = CI for anisotropic inclusions is defined as ,TIC R C R   (2.16) where C∗ is conductivity matrix in local coordinate system, The weak form associated is given by finding such that 14 ( ) d 0,q T d r T          By substituting the temperature field defined in (8) into the weak form (16) we obtain the discrete system of linear ordinary equations  d , T K Q d T a  where d are nodal unknowns and K and Q are the global stiffness matrix and external flux, respectively. More precisely, the matrix K and vector Q are defined by ( ) d , d .T TK B c x B Q N r        where B and N are the matrices of shape function derivatives and shape functions associated with the approximation scheme (8). The periodic boundary condition is finally introduced to (17) by mean of multiplicator Lagrange. 2.2.2. Results of thermal solution, using XFEM The calculation results indicate the temperature change in the calculation model Hình 2.5: Temperature distribution in the model CHAPTER 3. POLARIZATION APPROXIMATION (PA) 3.1. Introdution Consider a representative volume element (RVE) of an isotropic n-icomponent material that occupies spherical region V of Euclidean space. The center of the sphere is also the origin of the Cartesian system of coordinates {x}. The RVE consists of n components occupying regions Va V of volumes va and having conductivities (thermal, electrical, etc.) ca(a = 1,. . .,n; the volume of V is assumed 15 to be the unity). Starting from the minimum energy principles and using Hashin–Shtrikman polarization trial fields, one derives the following three-point correlation bounds on the effective conductivity ceff of the composite (Le & Pham, 1991; Pham, 1993) 1 1 0 ** 0 **(2 ) [ (2 ) ] . eff c cP c c c P c c      (3.1) where 1 * * 1 * ( ) ( ) . n c v P c c c c          (3.2) , , , 3 ij ij V v A d             x (3.3) where conventional summation on repeating Latin indices (but not on the Greek indices) is assumed; Latin indices after comma designate differentiation with respective Cartesian coordinates; the arbitrary positive constant c0 is often referred to as the conductivity of a comparison material; the harmonic potentials ua(x) appear in the expressions of Hashin– Shtrikman polarization trial fields; the three- point correlation parameters Ab ac relate the microgeometries of the three phases Va, Vb, Vc. Simple property (polarization) functions P, being monotonously increasing functions of their arguments and sharing the same structure, shall take a central place in our bounds and estimates for an easy qualitative comparisons between them. If one takes c0 = cmax = max{c1,. . .,cn} (or c0 = cmin = min{c1,. . .,cn}), then ** 0c  (or ** 0c  ) and can be neglected to strengthen the inequalities in (1), and subsequently one obtains Hashin–Shtrikman bounds (2 ) (2 ) .effc max c minP c c P c  (3.4) With (3.5) in hands, one substitutes c0 = cM into (3.7) to obtain ** ** 0c c  . Then, the bounds (1) converge to the unique value of the effective conductivity of the model (2 ) .eff c Mc P c (3.5) 16 Then taking c0 = cM, one finds that ** ** 0c c  and deduces the polarization approximation (PA) for the effective conductivity of our n-component matrix-inclusion composite (generally in d dimensions) expressed through a property (polarization) function P: (( 1) ).eff PA c Mc c P d c   (3.6) 3.2. Result Using Ansys softwate mesh model for FEM, and author build a program to calculate by Matlab, show results as graphs. 3.2.1. Two-dimensional periodic three-component composites formed We examine two-dimensional periodic three-component composites formed from a continuous matrix phase and two inclusion phases, which have isotropic effective isotropic properties. The first one is bodycentered square periodic cell as given in the figure 1 and the second one is body-centered hexagonal microstructure as shown in the figure 2. The diameter of inclusions of each phases is taken such that vI2=vI3. The effective thermal conductivity ceff is computed with theparameters shown in the table 1. The obtained results are reported in the figure 3.a), 3.b), 3.c), 3.d) corresponding the data in the table 1.a) 1.b), 1.c), 1.d) respectively. Figure 3.1: body-centered square and hexagonal periodic three- component microstructure The thermal conductivity of inclusions and matrix: cM c1 c2 1 10 3 17 Figure 3.2: approximations and finite element results for the effective conductivity of the three-phase matrix mixtures. 3.2.2. Three-dimensional periodic three-component composites formed Figure 3.3: 3D cubic periodic three-component BCC Figure 3.4: 3D cubic periodic three-component FCC Example of calculation according to the data in the table (a) cM = 1 c1 = 3 c2 = 10 (b) cM = 3 c1 = 1 c2 = 10 (c) cM = 3 c1 = 10 c2 = 1 18 (d) cM = 10 c1 = 1 c2 = 3 Figure 3.5: Grahp of results for 3D solution 3.2.1. Effective medium approximations for the elastic moduli Building models for elastic problems according to Body-Center Cubic Example (a) KM = 4  M=2 KI2 = 1  I2=0.4 KI3 = 20  I3=12 (b) KM = 4  M=2 KI2 = 20  I2=12 KI3 = 1  I3=0.4 (c) KM = 1  M=0.4 KI2 = 4  I2=2 KI3 = 20  I3=12 (d) KM = 20  M=12 KI2 = 4  I2=2 KI3 = 10  I3=0.4 19 Fig 3.6: Graph of effective elastic modulu results CHAPTER 4. EQUIVALENT APPROXIMATION 4.1 Equivalent inclusion approach 4.1.1 Dilute solution for equivalent circle inclusions Presume one has particles of certain shapes from a particular component material, and the effective conductivity of a dilute suspension of those randomly oriented particles, having conductivity cα and volume proportion vα (α = 2,...,n) in a matrix of conductivity c1 = cM, has the form ( ) ( , ), 1.eff M M Mc c v c c D c c v       (4.1) 20 In the meantime, the dilute suspension of d-dimensional spherical particles having conductivity c¯α and volume proportion vα in the matrix of the same conductivity cM has the particular expression ( ) 1. ( ) , 1 eff M M M M dc c c v c c v c d c           (4.2) Equalizing (1) and (2), one finds 2 ( 1) ( ) ( , ) . ( ) ( , ) M M M M M M M dc d c c c D c c c dc c c D c c            (4.3) In special case the anisotropic inclusions have the elipse shape, has the particular expression 2( )(1 ) ( , ) . 2( )( ) M M M M M c c c r D c c c r c r c c             (4.4) and inclusions have ellipsoid (3D) shape 1 1 1 ( , ) . 3 (1 ) (1 ) (1 ) M M M M M c D c c c A c A c B c B c C c C                        (4.5) 2 2 0 0 2 2 2 2 0 ˆ ˆˆ ˆ ˆ ˆ , , ˆˆ2 ( ) ( ) 2 ( ) ( ) ˆˆ ˆ ˆ, ( ) ( )( )( ). ˆ2 ( ) ( ) a b c a b cdt dt A B a t t b t t a b c dt C t a t b t c t c t t                                       4.1.2 Materials with anisotropic inclusions We consider the two-component 2D square-periodic suspension of anisotropic inclusions having conductivity cI1 and cI2 in a matrix of conductivity cM     1 1 ,..., , ( 1) 2 . 2 ,..., , I Id M I M I Id M D c c c d c c D c c c     (4.6) 4.1.3 Equivalent inclusion approach spherical inclusions (platelet, fibrous) 21 Our spherical equivalent inclusion polarization approximation (SEIPA) for the effective conductivity of the composite would have particular expression: 1 2 . 2 3 ( )eff I MSEIPA M I M M v v c c c c c     (4.7) Our platelet equivalent inclusion polarization approximation (PEIPA) for the effective conductivity of the composite would have particular expression: 1 2 . 3 2 ( )eff I MPEIPA I I M I v v c c c c c     (4.8) Our fibrous equivalent inclusion polarization approximation (FEIPA) for the effective conductivity of the composite would have particular expression: 12 / 3 2 3 . 5 2 ( )eff I M I MFEIPA I M M I v v c c c c c c c       (4.9) 4.2 Result 4.2.1 Examples in 2D As the first numerical example, we consider some two-component 2D square periodic suspensions of elliptic inclusions having conductivity cI in a matrix of conductivity cM (Figure. 4.1) 22 Figure 4.1: Square periodic cell with elliptic inclusions Figure 4.2: Graphics of the effective conductivity Hashin– Shtrikman upper (lower) bound. (a) cM = 1, cI = 10; (b) cM = 10, cI = 1 4.2.2 Random suspension of elliptic inclusions in a continuous matrix The space of material can be entirely filled by spheres and fibres distributed randomly with dimensions varying to infinitely small such that the inclusion volume proportion can approach. Fìgure 4.3: (a) A random elliptic inclusion configuration; equivalent circular inclusion configuration, with the same inclusion volume 23 proportion. 4.2.3 Example in 3D Hình 4.4: Cubic periodic suspensions of prolate and oblate spheroid inclusions Figure 4.5: Graphics of the effective conductivity. 4.2.4 Anisotropic inclusions As a numerical example we consider the two-component 2D square- periodic suspension of anisotropic inclusions having conductivity cI1 and cI2 in a matrix of conductivity cM 24 Hình 4.6: The periodic cells of disorderly anisotropic inclusions having circular shape (a) The square cell; (b) The hexagonal cell; (c) The random cell Figure 4.7: The graphics of the effective conductivity for square cell. Hình 4.8: The graphics of the effective conductivity for hexagonal cell. Hình 4.9: Graphics of the effective conductivity with random cell. 25 4.2.5 Experimental (EXP) data Figure 4.10: Spherical equivalent inclusion polarization approximation (SEIPA) Figure 4.11: Platelet equivalent inclusion polarization approximation (PEIPA) Figure 4.12: Fibrous equivalent inclusion polarization approximation (FEIPA) CONCLUSION Main conclusions of the thesis are 1) The author has studied the extended FEM method, to overcome the difficulties of meshing with complex phase geometry (such as random distribution ellipse) to solve the problem of shifting 26 grid system, possibly. There is no need to pay attention to the intersection between phases with different conductivity. 2) Using Ansys software and building FEM programs with Matlab. The author has successfully developed the computing programs according to the eXtended Finite Element Method. 3) Construction of the near-interaction approximations for the conductivity and elastic moduli of macroscopically isotropic composites with spherical inclusions based on the polarization bounds of Pham (1995), in which the near-interactions between the particles and the surrounding matrix have been estimated exactly, while the far-interctions between the different particles are approximated. The approximations are simplr, always satisfy HS bounds, close to the numerical and experimental results. 4) The thesis has constructed equivalent inclusion approximatation. Using the available approximations (such as the near interaction or the polarization ones), I determine the inclusion equivalent conductivity based on a comparison between the dilute solutions for the idealistic inclusion and real inclusion composites. That would bring the complex problems to the equivalent simple ones for all the ranges of inclusions’ proportions. I also extend the approximation to materials with inclusions of spherical, platelet, and fiber forms, using experimental reference at certain volume proportions of the inclusions. FURTHER DEVELOPMENTs After this thesis, the author would like to have more time and conditions to continue research on anisotropic materials, as well as to solve more complex inclusion composites, close to the reality. I mant to refine the XFEM method in applications to ensure that the results would be reliable in case no experimental results are available. Meanwhile, the experimental results are still needed for comparisons, whenever they are available. 27 PUBLICATIONS OF THE AUTHOR