Finite element solution strategy to analyze heterogeneous structures
Strategija analize heterogenih struktur z metodo končnih elementov
Martin Lamut
1CO Vesolje-SI, Aškerčeva 12, 1000 Ljubljana, Slovenia Corresponding author. E-mail: martin.lamut@space.si
*
Received: September 19, 2011
Accepted: November 2, 2011
Abstract: In this contribution a general strategy for solving a coupled micro-macro problems is presented which enables analyses of modern heterogeneous materials. It provides an efficient problem solving tool to structures with complex microstructures, used in a demanding structural components. The method uses a nested finite element solution strategy called multilevel finite element approach-ML-FEM. Within the ML-FEM framework one conducts an embedded micro-scale computation in order to obtain quantities required at the macroscopic level. The application of ML-FEM circumvents the need to construct an explicit macroscale constitution formulation, considering increased computational costs. Increased computation is linked to detailed microscopic analysis for which the statistical representative volume element-RVE is needed. RVE will be derived based on the convergence criterion. In this work a general method for calculation of the consistent macroscopic stiffness matrix via sensitivity analysis of a micro level is shown. As an example the proposed method is applied on a simple test specimen under compression consisting microstructures with porosities and stiff inclusions.
Povzetek: V tem delu je bila razvita splošna strategija za reševanje vezanih mikro-makro sodobnih heterogenih materialov. Strategija je učinkovito orodje pri reševanju problemov s kompleksno mikrostrukturo, uporabljeno v zahtevnih inženirskih komponentah. Strategija uporablja večnivojski način reševanja problemov, kjer na mikroskopski in makroskopski ravni poteka analiza z metodo končnih elementov
(ML-FEM). Pri tej metodi reševanja makroskopska konstitutivna zveza ni več potrebna, saj je le-ta na račun povečanega računskega časa pridobljena z natančno mikroskopsko analizo. Ta je izvedena na statističnem reprezentativnem volumnu (RVE), katerega velikost določimo s konvergenčnim merilom. Metoda je splošen način reševanja makroskopske togostne matrike preko občutljivostne analize mikroskopskega nivoja. Lastnosti metode so bile preizkušene na enostavnem tlačnem preizkusu za porozno mikrostrukturo in mikrostrukturo s togimi vključki.
Keywords: Heterogeneous materials, multiscale analysis, macroscopic tangent computation, sensitivity analysis
Ključne besede: heterogeni materiali, mikro-makro analize, makroskopska togost, občutljivostna analiza
Introduction
Heterogeneous materials used in engineering sciences have physical properties that vary throughout their microstructures. Heterogeneities, such as inclusions, pores, fibers and grain boundaries, have a significant impact on the observed macroscopic behavior of multi-phase materials. In engineering some typical examples are metal alloy systems, various composites, porous and cracked structures, polymeric blends and polycrystalline materials.
To describe the macroscopic overall characteristics of heterogeneous structures is a vital problem in many engineering applications. The ability to convey information across length scales is essential for a better understanding of the sources of physical behavior observed on higher scales. Using mi-
cromechanical models of the microstructural elements, homogenization techniques allow an efficient and correct transfer of microscale information to the macroscale analysis. The fundamental methodology of homogenization is the characterization of the macroscopic behavior of the heterogeneous material by appropriately identifying and testing a statistically representative micromechanical sample. Once an appropriate sample is identified it can be used in the multiscale analysis methodology. The most straightforward way is to use the multilevel finite element method ML-FEM[1-5]. When analyses at both levels are made in the context of FEM, it can be referred to as the FE2 method[6, 7]. The application of ML-FEM circumvents the need to construct an explicit macroscale constitution formulation, though at an increased computational cost. The con-
stitutive equations are written only on microscopic scale and homogenisation and localization equations are used to compute the macroscopic strains and stresses knowing the mechanical state at microscopic level.
By analyzing the engineering structure, the point of interest is usually localized in the so called critical region, where detailed analyses are needed. So to further increase the efficiency of the computation the structure can be divided into subdomains, critical region and the rest of the structure. In the critical region an embedded ML-FEM computation is conducted, while elsewhere a classical homogenization technique is used. In either case a statistical mi-cromechanical model or representative volume element (RVE) will be needed.
The purpose of this contribution is mainly two fold. First, the statistical RVE size will be derived based on convergence criterion of the several parameters being monitored. The second purpose of this work tackles the efficiency of multilevel computation. Since a conventional way of macroscopic tangent computation in a condensation procedure, necessitate the computation of a Shur complement. It inflicts for increasingly complex microstructure higher memory allocation demands that may not be met by today's computers. Therefore, as an alternative, a tangent computation tech-
nique based on a sensitivity analysis of a microscopic level will be presented.
Methods
Numerical RVE size
In order to estimate the effective properties of heterogeneous material, most of the micro-macro methods assume the existence of a micromechanical sample that is statistically representative of the microstructural features. The usual approach[8] is to determine a relation between averages, E\ defined through (o)RVE = E*(e)RVE . Here o and e are the stress and strain fields within a statistically representative volume element. The RVE is considered both smaller enough than the macro scale media and bigger enough than the heterogeneities on the micro scale, without introducing non-existing properties (e.g. anisotropy).
In this contribution, macroscopically isotropic materials are considered, therefore the two linear elastic constants (bulk and shear moduli) describing the form of E'can be computed using:
(1.1)
¡(^RVE-j^RVE
V (s')xMv
where a' and e' denotes the deviatoric part. Macroscopically isotropic heterogeneous structure is achieved by random particle distribution at the mi-croscale. Therefore, for a given sample size, multiple distributions of particles are possible. In order to capture a statistical measure of the range of responses from different distributions, a simple averaging of three samples per RVE size was used.
To model random porous microstructures a matrix containing randomly distributed pores throughout a square L x L was considered. The size of the particles were determined relatively to unit length of the RVE such that 0.1 < 2r < .15. Mechanical properties of the matrix material was K = 167 GPa and G = 77 GPa. In order to determine a suitable RVE size, one must monitor the range of estimates to E" for successively larger samples, shown on Figure 1. The following sequences of particles per sample are used (N): 2, 4, 15 and 32. Relying on the expectation as RVE size increases indefinitely the effective properties of material constants (K, G) will converge towards E\
For numerical simulation of the response a 2D quadrilateral plane strain 2 x2 Gauss rule elements were used. To determine the effective bulk and shear moduli, since the effective response is assumed isotropic, only one
test loading is necessary en = 0.01005. In Table 1, the perturbation magnitudes are shown for various quantities as a function of pore number in the sample.
Besides convergent material properties, the RVE must be tested upon the influence of the microstructural geometry
Figure 1. A series of test samples with increasing size, the volume fraction of particles is fixed at 0.6 %.
Table 1. Perturbation magnitudes for shear and bulk moduli as a function of particles number (N).
		K	G
N	RVE size	K32 ~ K V	Gn ~GN
		Kn	Gn
2	0.73	0.0 IX	0.018
4	1	0.021	0.023
15	2	0.013	0.013
32	3	0.002	0.005
properties. This can be done by tracking various quantities such as: strain energy function, maximal stresses, averaged stresses in the particles or matrix etc. In this work the maximal effective stress was considered (von Misses). To guarantee the mixed stress fields besides the previously used normal test loading en = 0.01005 the shear loading condition was used e12 = 0.01005 all the rest stays
the same as described previously. Figure 2 is showing the convergence of the max. effective stress in the RVE by increasing its size.
Based on the tests the statistical RVE size 2 (approximately 15 particles) is chosen. This size is used in all subsequent analysis. The outline of the determination of the RVE size and ef-
size
Figure 2. Max. effective stress for two loading cases depending upon RVE size
Table 2. Perturbation magnitudes of max. oeff for normal and shear load condition as a function of particles number (N).
n	RVE size	fb	fs <V.32 -<Vw .. mn0''
		32	°V.32
2	0.73	14,4	16.2
4	1	3.6	3.8
15	2	2.9	1.5
32	3	0.15	1.5
fective material constants are made for particles representing voids. The same procedures are used to determine the size and effective material constants also for microstructure with stiff inclusions. All the values are given in the Results & discussion section.
Multilevel finite element - principle
The material under consideration is assumed to be macroscopically homogeneous, so that continuum mechanics can be used to describe macroscopic behavior. However, at the micro level is the material configuration heterogeneous, consisting of many distinguishable components e.g. grains, cavities, hard inclusions. The microscopic length scale is much smaller than the characteristic length of the macroscopic structure, therefore in this case a hypothesis on periodicity of the microstructure is acceptable. The multilevel finite element strategy may be described by the following subsequent steps: [9, 10]
Figure 3. Shematic diagram of the ML-FEM model.
1.	Determination of a statistical representative volume element (RVE), used in homogenization and in embedded multilevel analysis.
2.	The macroscopic structure to be analyzed is discretized by finite elements. The external load is applied by an incremental procedure.
3.	Macroscopic deformation gradient tensor (FM) is calculated for every integration point of the macrostructure in a multilevel subdomain.
4.	From the macroscopic deformation tensor appropriate boundary conditions are derived to be imposed on the RVE that is assigned to this point.
5.	Upon the solution of the boundary value problem for the RVE, the macroscopic stress tensor (PM) is obtained by averaging the resulting RVE stress field over the volume of the RVE.
6.	Additionally, the local macroscopic consistent tangent is derived from the sensitivity analysis of the RVE.
This framework is schematically illustrated in Figure 3. In the subsequent sections these issues are discussed in more detail.
Linking macroscopic and microscopic levels
The actual coupling between the macroscopic and microscopic scales is based on averaging theorems. The energy averaging theorem, known in
the literature as the Hill condition or macrohomogeneity condition111, 12], requires that the macroscopic volume average of the variation of work performed on the RVE is equal to the local variation of the work on the macroscale. Formulated in terms of a deformation gradient tensor and the first Piola-Kirchhoff stress tensor, the work criterion in differential form is written:
(2.1)
In words, this equality states that in the transition from the microscopic scale to the macroscopic scale, energy is conserved.
It is well known that this criterion is not satisfied for arbitrary boundary conditions (BC) applied to the RVE. Classically three types of RVE boundary conditions are used, i.e. prescribed displacements, prescribed tractions and prescribed periodicity. Periodicity here is referring on an assumption on global periodicity of the microstructure, suggesting that the whole macroscopic
specimen consists of spatially repeated unit cells. Among them the periodic BCs show a more reasonable estimation of the effective properties. This was supported and justified by numbers of authors[13-16]. The periodicity conditions for the microstructural RVE are written in a general format as:
(2.2)
where x and X represents the actual
and initial position vector and p the boundary traction of the RVE. In the equation (2.1) the macroscopic first Piola-Kirchhoff stress tensor (PM) and the macroscopic deformation gradient tensor (Fm) are the fundamental kinetical and kinematical measures which are defined in terms of the volume average of their microscopic counterparts. Every time that the work criterion is satisfied, the volume average of the macroscopic above mentioned measures can be obtained through the knowledge of boundary information only.
(2.3)
(2.4)
Here, is undeformed RVE, P
? RVE	? m
and Fm are microscopic stress tensor and deformation gradient tensor, respectively, r represents the boundary of the RVE, while N represents the normal vector to the surface of RVE.
Macroscopic tangent computation
In the realization of the multilevel FEM approach, the macroscopic constitutive formulation is not explicitly obtained from the experimental data. Instead, the needed stiffness matrix at every macroscopic integration point has to be determined directly from the numerical relation of the macroscopic stress (PM) and macroscopic deformation gradient (FM) at that point[15, 17, 18]. The weak form of the macroscale problem in the absence of body forces and acceleration can be written in variational form as:
m = \SF„'RJV
M 1 Ml
(2.5)
To solve the macroscopic primal problem within ML-FEM setting, at i-th iteration step of a standard New-ton-Raphson solution scheme, the following linearization needs to be computed.
The macrolevel element tangent stiffness matrix and the residual force vector can be obtained with the knowledge of the stress (PM) and macroscopic tangent (SPM/dFM) obtained from the RVE analysis, since FM is explicit function of node displacements. The (PM) can be obtained directly from RVE analysis by using averaging theorem equation (2.4), while for the determination of the macroscopic tangent a RVE sensitivity analysis is performed. For the sensitivity prob-lem[19] the residuals and the vector of unknowns are defined as a function of sensitivity parameters, which are in this case the elements of tensor (FM). The sensitivity problem can then be obtained from the primal problem by differentiating the response functional and the residuals with respect to macroscopic deformation gradient (FM), and the following system on the microlevel has to be solved:
(2.7)
where, Wm represents response functional on the microlevel, a is a set of unknowns (displacements), while ^ represents arbitrary sensitivity param-
eter in our case F
M'
(2.6)
With the assembling of the macroscopic stiffness matrix is the problem on the macro level fully described and can be solved to produce an update of the macroscopic displacement field.
Remark 1: for consistency the particular type of BC employed for the comp. of K must match the type of BC employed in the computation of P.
Results & discussion
In order to evaluate the presented ML-FEM strategy a simple compres-
sion test, Figure 4a, of a homogeneous matrix material with 6 % volume fraction of randomly distributed voids or stiff inclusions has been examined. The material parameters used in the analysis are: shear modulus of matrix material and stiff inclusions are G =
m
77 GPa and Gi = 307 GPa respectively, bulk modulus of matrix material and stiff inclusions are K = 167 GPa
m
and Ki = 667 GPa respectively. For the homogenized part the effective material constants are calculated from RVE tests presented in the Numerical RVE size section: voided microstructure G = 72 GPa and K=156 GPa, microstruc-
Figure 4. a) Axisymmetric model of compresion test. b) Axisymmetric numerical model: two subdomains with critical points indicated.
ture with stiff inclusions G = 90 GPa and K = 194 GPa. On the macroscopic level a 2D quadrilateral plane strain 2 x 2 Gauss rule elements were used. Load was applied incrementally: load displacements A = 1 unit relative to L. The multilevel algorithm has been implemented into computer program AceFEM[20], where a special macroscopic element can be readily defined in open source code.
Figure 4b shows a discretized numerical model where the homogenized and multilevel subdomains are clearly indicated. In the later, two critical points are marked where a detailed RVE analysis has been done. From Figure 4b a
straightforward estimate, regarding the considered test, can be done about the amount of computation needed for each macroscopical load increment: number of elements in the multilevel subdomain times the Gauss points per elements. For the present test it takes, 18 Elements x 4 Gauss points, RVE analyses for each load increment. In order to further speed up the analysis and to make it more useful for complicated engineering applications, multilevel algorithm was set up for parallel computations.
In Figure 5 the contour plots of the equivalent Misses stress considering the two microstructures are compared.
a)	b)
Figure 5. Effective stress in the macrostructure 2D a) stiff inclusions, b) voided microstructure.
As expected the microstructure with stiff inclusions produces higher stresses but the interesting point is that the contour plots are not the same which would be the case for completely homogenized structure. The influence of more realistic model on the microscopic level is clearly visible even in a simple case considered.
The detailed analyses of the RVE at the critical points for both microstructures
Point o
are depicted in Figures 6 and 7. It can be seen that the voids act as a stress concentrator and that some stress concentration regions can be seen between neighboring voids, Figure 6. The RVE taken from macroscopic point (o) is subjected to mainly hydrostatic compression stresses, while RVE from point (-) is exposed to some amount of deviatoric stresses as well. This phenomenon has greater influence on the voided microstructure, where higher
Point-
Figure 6. Effective stress in the RVE from the critical point of the macrostructure voided microstructure.
Point o	Point -
Figure 7. Effective stress in the RVE containing stiff inclusions from the critical point of the macrostructure.
effective stresses are observed in the mainly hydrostatic region. Generally by comparing the microscopic and macroscopic stress fields a substantially higher stresses are observed on the microscopic level. So by simultaneously examining RVE at critical macro points, while deforming the macro structure, a deeper understanding of deformation mechanisms can be obtained, which can be very helpful in studying the damage mechanisms and densification in various engineering materials.
Conclusions
Firstly it was shown how to determine a statistical representative volume element (RVE), which was later used in the homogenization process and as well in embedded multilevel analysis. A multilevel finite element analysis strategy for the simulation of the mechanical behavior of heterogeneous materials has been outlined. The performance of the method was illustrated by the modeling of simple compression test. Two different microstructures were tested: first including voids and second with hard inclusions embedded in the matrix. The presented multilevel finite element analysis strategy provides an efficient approach to determine the macroscopic response of heterogeneous materials with accurate account for microstructural phenomena. In the
ML-FEM strategy the computational efficiency hinges on the correct and effective macroscopic tangent computation, in this work this is done by sensitivity analysis of the microscopic level. It enables a problem solving tool for a variety of different micro-macro problems which includes complex microstructures.
References
[1]	Kouznetsova, V., Brekelmans, W. A.
M., Baaijens, F. P. T. (2001): An approach to micro-macro modeling of heterogeneous materials, Comp. Mech., Vol. 27, pp. 37-48.
[2]	Temizer, I., Wriggers, P. (2008): On
the computation of the macroscopic tangent for multiscale volumetric homogenization problems, Comput. Methods Appl. Mech. En-grg, Vol. 198, pp. 495-510.
[3]	Smit, R. J. M., BREKELMANS, W.
A. M., Meijer, H. E. H. (1998): Prediction of the mechanical behavior of nonlinear heterogeneous systems by multi-level finite element modeling, Comput. Methods Appl. Mech. Engrg., Vol. 155, pp. 181-192.
[4]	Miehe, C., Koch, A. (2002): Compu-
tational micro-to-macro transitions of discretized microstructures undergoing small strains, Archive of Applied Mechanics, Vol. 72 pp. 300-317.
[5]	Miehe, C. (2003): Computational mi-
cro-to-macro transitions for dis-
[9]
[10]
cretized micro-structures of heterogeneous materials at finite strains based on the minimization of averaged incremental energy, Comput. [14] Meth. Appl. Mech. Eng., Vol. 192, pp.559-591.
[6]	Feyel, F., Chaboche, J. L. (2000): FEE2
multiscale approach for modelling the elastoviscoplastic behaviour of [15] long fibre SiC/Ti composite materials, Comput. Methods Appl. Mech. Engrg., Vol. 183, pp. 309330.
[7]	Feyel, F. (1999): Multiscale FE2 elas-
toviscoplastic analysis of composite structures, Comput. Mat. Sci., Vol. 16, pp. 344-354.	[16]
[8]	Zohdi, T. I., Wriggers, P. (2001): Com-
putational Micro-macro Material Testing, Arch. Comput. Methods Eng., Vol. 8, pp. 131-228. Lamut, M. (2010): Mikro in makrome-hansko modeliranje kompozitnih [17] materialov. Ph. D. Thesis. Ljubljana: University of Ljubljana 2010; 87 p. Lamut, M., Korelc, J., Rodič, T. (2011): Multiscale modelling of heterogeneous materials, Math. [18] Tech., Vol. 45, pp.421-426.
[11]	Hill, R. (1963): Elastic properties of
reinforced solids : some theoretical principles, J. Mech. Phys. Solids., Vol. II, pp. 357-372.
[12]	Hill, R. (1972): On constitutive
macro-variables for heterogeneous [19] solids at finite strain, Proc. R. Soc. Lond. A., Vol. 326, pp. 131-147.
[13]	Terada, K., Hori, M., Kyoya, T.,
Kikuchi, N. (2000): Simulation of the multi-scale convergence in [20] computational homogenization
approaches, International Journal of Solids and Structures, Vol. 37, pp. 2285-2311.
Khisaeva, Z., F., Ostoja-Starzewski, M. (2006): On the Size of RVE in Finite Elasticity of Random Composites, J. Elasticity, Vol. 85, pp.153-173.
Sluis, O., Schreurs, P. J. G., Brekelmans, W. A. M., Meijer, H.E.H. (2000): Overall behaviour of heterogeneous elastoviscoplastic materials: effect of microstructural modelling, Mechanics of Materials, Vol. 32, pp.449-462.
Temizer, I., Wriggers, P. (2008): On the computation of the macroscopic tangent for multiscale volumetric homogenization	problems,
Comput. Methods Appl. Mech. Engrg., Vol. 198, pp. 495-510.
Temizer, I., Zohdi, T. I. (2007): A numerical method for homogenization in non-linear elasticity, Comput. Mech., Vol. 40, pp.281-298.
Miehe, C. (1996): Numerical computation of algorithmic (consistent) tangent moduli in large-strain computational inelasticity, Comput. Methods Appl. Mech. Engrg., Vol. 134, pp. 223-240.
Korelc, J. (2009): Automation of primal and sensitivity analysis of transient coupled problems, Comput. Mech., Vol. 44, pp. 631649.
Korelc, J. (2007): AceGen & AceFEM user manual, FGG-UL.