Comput Mech (2012) 50:413–431 DOI 10.1007/s00466-012-0681-2

ORIGINAL PAPER

Numerical homogenization of heterogeneous and cellular materials utilizing the finite cell method Alexander Düster · Hans-Georg Sehlhorst · Ernst Rank

Received: 7 June 2011 / Accepted: 7 January 2012 / Published online: 26 January 2012 © Springer-Verlag 2012

Abstract We present a new approach for the numerical homogenization of cellular and heterogeneous materials. The procedure is based on the finite cell method, which is applied to efficiently discretize representative volume elements for which effective material properties are computed. The starting point for our homogenization might be either a computer-aided design of a heterogeneous material or a three-dimensional computer tomography (CT-scan) of the specimen of interest. A fully automatic discretization in terms of finite cells, applying a hierarchic extension process to control the discretization error, is utilized to solve the corresponding boundary value problems arising during the homogenization. Special emphasis is placed on the numerical treatment of boundary conditions. To this end we apply the window method, which can be interpreted as a variant of the self-consistency method. Several numerical examples ranging from porous materials to fiber-reinforced composites will be presented, demonstrating the efficiency of our approach. The homogenization procedure will be also applied to a foam, a CT-scan of which is available. Keywords Numerical homogenization · Heterogeneous materials · Porous materials · Composites · Foams · Finite cell method

A. Düster (B) · H.-G. Sehlhorst Numerische Strukturanalyse mit Anwendungen in der Schiffstechnik (M-10), Technische Universität Hamburg-Harburg, Hamburg, Germany e-mail: [email protected] E. Rank Computation in Engineering, Technische Universität München, Munich, Germany

1 Introduction During the past few years, cellular and heterogeneous materials have been the object of growing interest in different fields of application such as civil, mechanical or aerospace engineering as well as in naval architecture. Cellular materials are appealing because they combine properties like low mass density with high stiffness and are therefore utilized to build light-weight structures. Due to the complicated microstructure of cellular materials, which can lead to deformation-induced anisotropic [20,46] and size-dependent effects [2,5,6,23,27,38,50], the modeling and simulation of engineering structures composed of these materials still remains a challenge and accordingly the subject of ongoing research. The situation is very similar with heterogeneous materials. Materials containing pores or particles/fibers with a different mechanical behavior from that of the surrounding matrix are difficult to model and compute. Since a direct numerical approximation of engineering structures composed of cellular or heterogeneous materials is computationally prohibitive, homogenization procedures are of great importance. An overview of analytical homogenization procedures can be found in Aboudi [1], Nemat-Nasser and Hori [36], Zhodi and Wriggers [52] or Barbero [3], for example. Numerical homogenization techniques for micro-macro scale transition have been developed by different research groups: see [13,14,18,24,28,31,37,52], for instance. The FE2 method [12] is a commonly used two-scale numerical homogenization approach, where a microstructure, either a representative volume element (RVE) or a testing volume element (TVE), as introduced by Nemat-Nasser and Hori [36], is embedded into a macroscopic finite element framework via projection and homogenization schemes. The application of the FE2 method to classical continua can

123

414

be found in [12,32,31,33,34,43–45], for gradient-enhanced homogenization in [25,26] or for Cosserat and micromorphic continua in [11,19,23], among others. Although the FE2 method is a very versatile and general approach, which makes it possible to study the mechanical behavior of materials with a complicated microstructure, it has a major drawback which is related to the solution of nested boundary value problems of solid mechanics. Since a suitable microstructure and corresponding boundary value problem has to be solved for each integration point of the finite elements applied on the macro scale, this approach requires a great deal of computation work, particularly when three-dimensional or non-linear problems are involved. The approach adopted for this project was therefore based on the computation of effective material properties, which can be also interpreted as a simplified FE2 method. The effective material properties in this case were computed by means of a numerical differentiation, where the elasticity tensor is obtained via differentiation of the homogenized stresses with respect to the homogenized strains of the microstructure. This approach is similar to a testing method where the microstructure is subjected to several load cases [3,52] resulting in Dirichlet boundary value problems, which need to be solved to determine the effective material properties. Since the quality of the results of the homogenization is largely influenced by the way the boundary conditions are applied to the microstructure [46], special care has be taken in this respect. To this end, we adapted the window method, which Hain and Wriggers [18] successfully applied for the numerical homogenization of hardened cement paste. In order to compute the effective material properties, it is necessary to discretize the RVE representing the microstructure. The most common approach adopted for this purpose is the finite element method. In [24], for example, the numerical homogenization of randomly distributed short cylindrical fiber composites is based on the discretization of the RVE using 10-noded tetrahedral elements with full integration. The meshes resolve the geometry of the particles which are integrated into a surrounding matrix. A different approach was chosen in [29], where hexahedral elements were employed to discretize the RVE as part of the homogenization procedure of heterogeneous materials with inclusions. The meshes are refined locally in order to resolve the particles, leading to elements with hanging nodes, which are constrained to obtain a continuous displacement field. In general, the meshing of complex microstructures is still a challenging problem since it can lead to labor intensive pre-processing, resulting in either very fine meshes or discretizations, which are coarser but composed of distorted elements of reduced approximation quality. These problems are even more pronounced when considering micro-

123

Comput Mech (2012) 50:413–431

structures of increasing complexity, as found in the case of modern composites or foam-like structures. Cartesian grids, which have been used in computational mechanics for a long time, appear to be a promising approach for overcoming—or at least easing—this discretization problem. To homogenize fiber-reinforced materials, Aboudi [1] proposed the method of cells (MOC) where a double periodic array of fibers is discretized with a grid of 2 × 2 sub-cells in two dimensions. This concept was subsequently generalized, so that the geometry of the fibers can be resolved with a Cartesian grid composed of an arbitrary number of sub-cells, where each sub-cell represents a trilinear displacement field and approximates the geometry of the fibers or the embedding matrix with its cubic shape. Matzenmiller and Kurnatowski [30] give an overview of the different versions and enhancements of the MOC and compare them with corresponding finite element computations, leading to the conclusion that the low-order cell approach delivers an accurate approximation of the stiffness of the homogenized composite. The novelty of the proposed approach is the further development of the recently proposed finite cell method (FCM) [9,39] and its application to the efficient solution of homogenization problems. The FCM can be interpreted as a combination of a fictitious domain method [7,35,40] with highorder finite elements [48]. Applying the FCM dramatically simplifies the meshing process, since a simple Cartesian grid is used and the geometry of the microstructure and the different material properties are taken care of during the integration of the stiffness matrices of the cells. In this way, the effort is shifted from mesh generation towards the numerical integration, which can be performed adaptively in a fully automatic way. It will be demonstrated by several numerical examples that this approach can be applied to analyze a computer-aided design of a new heterogeneous material or data obtained by a three-dimensional computer-tomography of the material of interest. The main advantage of the proposed method is that materials with complicated microstructures can be discretized automatically without creating meshes which result in a prohibitive amount of elements. In order to efficiently homogenize the elastic behavior of complex materials, we extend the FCM and combine it with the window method utilized to stipulate boundary conditions for the microstructure. This makes it possible to efficiently compute effective material properties of materials with a complex microstructure. The paper is arranged as follows: following a short introduction to the notation, Sect. 3 describes the FCM and its application to discretize the microstructure of cellular or heterogeneous materials. The procedure applied to perform the homogenization of the RVE is described in Sect. 4. Several numerical examples demonstrating the properties of the proposed homogenization approach are given in Sect. 5, and a

Comput Mech (2012) 50:413–431

415

conclusion plus an outlook on future work round off the paper in Sect. 6. 2 Notation This section serves to summarize the notation of tensors, selected operators and matrices employed in this work. Please note that tensors and tensor products are used according to the Cartesian coordinate system with unit base vectors E , E and E and the Einstein summation convention is ∼1 ∼2 ∼3 applied.

 B(u, v) =

[L v]T C(x) [L u] dΩ Ω

is a bilinear form in which u denotes the displacement, v represents the test function, V is the space of admissible test functions, L is the standard strain-displacement operator and C(x) is the elasticity matrix which is not constant due to the different material properties of the matrix material and the inclusions. The linear form   T (3) F(v) = v f dΩ + v T ¯t dΓ Ω

Tensors Scalar value First-order tensor (vector) Second-order tensor Fourth-order tensor Tensor operations Dyadic product Inner product Tensor product Matrices Local finite element quantities Global quantities

A X = Xi E ∼ ∼i S = Si j E ⊗E ∼ ∼i ∼ j C = C E ⊗E ⊗E ⊗E i jkl ∼ ∼i ∼ j ∼ k ∼l u ⊗ ∼v = u i v j E ⊗E ∼i ∼ j S · F = S F ij ij ∼ ∼ S F = Si j F jk E ⊗E ∼ ∼ ∼i ∼ k

∼

(2)

ΓN

takes the volume loads f and prescribed tractions ¯t into account. A standard finite element analysis would call for a discretization of the domain Ω with finite elements. In the FCM, the main idea is to embed the physical domain Ω into a domain Ωe of simpler shape, see Fig. 1. Therefore the integration of the weak form Be (u, v) = Fe (v)

(4)

has to be extended over Ωe , i.e. the bilinear form reads  Be (u, v) = [L v]T Ce [L u] dΩ (5)

u, J r, R

Ωe

3 The finite cell method

where

The approach utilized for the spatial discretization of the heterogeneous material is an extension of the FCM that was introduced in [9,39] for two and three-dimensional problems of solid mechanics. 3.1 Weak formulation and discretization Let us consider a three-dimensional body consisting of a heterogeneous material. For the sake of clarity, we give a twodimensional sketch of the FCM concept in Fig. 1, although all computations presented in this paper have been conducted for three-dimensional problems. The material to be homogenized may have pores and inclusions with different material properties from that of the matrix. Figure 1 indicates the pores by means of ellipses and inclusions with circles. Note that the proposed method does not impose any restrictions on the topology or geometry of the pores and inclusions, which might have an arbitrary shape. We assume that the material is linear elastic and inclusions to be perfectly bonded to the surrounding material of the matrix. The underlying equation for the FCM is the weak formulation of equilibrium B(u, v) = F(v) ∀ v ∈ V which is defined in the physical domain Ω. In Eq. 1

(1)

 Ce =

C(x) ∀ x ∈ Ω α0 Cm ∀ x ∈ Ωe \ Ω

(6)

is the elasticity matrix of the embedding domain and α0 > 0 is a scalar value which scales the elasticity matrix Cm of the matrix material. Setting α0 = 0 corresponds to a void material, i.e. Ce = 0. In this way points x ∈ Ωe \ Ω which are not located in the physical domain do not contribute to the weak form since the corresponding elasticity matrix vanishes. In order to circumvent conditioning problems which will be addressed in the next section, we set α0 > 0 to a very small value which is close to 0. The linear form of the embedding domain   T (7) Fe (v) = α(x) v f dΩ + v T ¯t dΓ Ωe

ΓN

with  α(x) =

1.0 ∀ x ∈ Ω α0 ∀ x ∈ Ωe \ Ω

(8)

is defined in much the same way. It is easy to see that choosing α0 = 0 and inserting (8) and (6) into Eqs. 5 and 7 yields the weak form defined in the physical domain Ω, i.e.

123

416

Comput Mech (2012) 50:413–431

Fig. 1 A heterogeneous material with pores and inclusions. The physical domain Ω is embedded in Ωe

processing and mesh generation for obtaining meshes that are conform with the boundaries and aligned to the material interfaces becomes very involved. In addition to the laborintensive pre-processing, the computational time needed for the mesh generator and the finite element analysis is likely to be very high, since a complex heterogeneous material would yield a very dense mesh with an enormous amount of finite elements. In the FCM, we apply a Cartesian mesh which consists of hexahedrals based on hierarchic shape functions of arbitrary polynomial order. Performing the discretization  c Ωc with n c cells, the of the embedding domain Ωe = nc=1 bilinear form reads nc   (11) Be (u, v) = [L vc ]T C e [L uc ] dΩ. c=1 Ω c

The linear form (10) can be treated in a similar manner. Applying a standard Bubnov- Galerkin method, the displacement and test function in each cell c are approximated as

Fig. 2 The embedding domain Ωe

 Be (u, v) =

[L v]T Ce [L u] dΩ Ωe



=

[L v]T C(x) [L u] dΩ Ω



+

[L v]T 0 [L u] dΩ = B(u, v)

(9)

Ωe \Ω

and



Fe (v) =

 α(x) v f dΩ + T

Ωe

v T ¯t dΓ = F(v).

(10)

ΓN

The main advantage of embedding the physical domain Ω in Ωe is that the domain Ωe can be easily discretized with a Cartesian mesh, see Fig. 2. Since the elements of this mesh neither correspond to the boundary of Ω nor are they aligned to the material interfaces, we call them finite cells in order to distinguish them from finite elements. The finite cells may be intersected by the boundary of the physical domain or may have different material properties due to the inclusions, which are also disregarded in the Cartesian mesh. Cells which are completely outside of the physical domain Ω are simply discarded, see Fig. 2. An alternative method of applying finite elements, of course, would be to use mesh generators producing unstructured meshes which are usually based on tetrahedrals. Since the geometry of heterogeneous materials might be very complex, it can happen that the relevant pre-

123

uc = N c U c

(12)

vc = N c V c

(13)

where N c denotes the matrix of shape functions, U c is the vector of unknowns and V c represents the coefficients of the test function. The product LN c is denoted as the standard strain-displacement matrix Bc . Our implementation, [8], uses hexahedrals with hierarchical shape functions based on integrated Legendre polynomials, [10,47,48]. This makes it possible to perform either an h-extension with an arbitrary polynomial degree by successively reducing the cell size or a p-extension on a fixed mesh of (coarse) cells. Inserting (12) and (13) into the weak formulation yields the linear equation system KU = F,

(14)

nc K c , the global stiffness matrix, and F = where K = Ac=1 nc Fc the global load vector are obtained by assembling Ac=1 all cell matrices. Due to the similarity between the FCM and the FEM we can take advantage of the existing algorithms and solution methods.

3.2 Numerical integration of cell matrices As discussed in the previous section, the main idea of the FCM is to simplify mesh generation, which in the case of a standard finite element analysis of a heterogeneous material, would produce an enormous amount of elements resulting in a very high number of unknowns. Another, maybe more important fact is that the pre-processing might be quite complicated, resulting in a labor-intensive procedure for generating the finite element model. The FCM significantly reduces

Comput Mech (2012) 50:413–431

the amount of pre-processing required, since mesh generation is significantly simplified. However, the resulting mesh will generally have cells that are intersected by the boundary of the domain Ω or may have different material properties due to the fact that the inclusions in the matrix material are not geometrically resolved by the mesh, see Fig. 2. The effort is accordingly shifted from mesh generation to the task of reliably computing the cell matrices, i.e. the stiffness matrix and load vector of the cells. This calls for a numerical procedure which accurately computes the cell matrices in a fully automatic, error controlled way. One way of improving the accuracy of the integration is to increase the number of Gaussian points employed to compute the stiffness matrix of one element or cell. This approach was coined by Zhodi and Wriggers [52] as the Gauss point oversampling method and utilized in the numerical homogenization of problems with material discontinuities. Due to its simplicity, the Gauss point oversampling method is a very attractive technique. However, in the case of discontinuities introduced by material interfaces, the method does not converge very fast. This is due to the fact that the Gauss quadrature is very efficient for smooth integrands but its convergence rate deteriorates in the case of discontinuous integrands. In order to increase the efficiency of the integration, we apply an adaptive integration scheme based on sub-cells. To be more specific, let us take a look at the computation of the stiffness matrix    BcT (ξ ) C e (x) Bc (ξ ) detJ c dξ dηdζ (15) Kc = ζ

η ξ

of cell c which is either intersected by the boundary of Ω or has different material properties due to inclusions that are geometrically unresolved by the finite cell mesh. In Eq. 15 the integration is carried out in the local coordinate system of the cell, which is the standard approach applied in finite elements. The hierarchic shape functions N c of the hexahedral cell c and accordingly the strain-displacement matrix Bc = LN c are defined in the local coordinate system (ξ, η, ζ ) of the element. Due to the mapping x = Qc (ξ, η, ζ ) of cell c from the local (ξ = [ξ, η, ζ ]T ) to the global (x = [x, y, z]T ) coordinate system, the determinant of the corresponding Jacobian matrix J c = gradT Qc (ξ, η, ζ ) has to be accounted for. As the integrand will be discontinuous, a standard Gaussian integration will not perform well since it presumes that the function to be integrated is a smooth polynomial. We therefore subdivide cells which are broken by the boundary of Ω or which have different material properties, into n sc sub-cells which are introduced just for integration purposes, see Fig. 3. On each of this sub-cells a Gaussian integration is performed with n gp Gaussian points in each direction. The refinement with respect to the number of sub-cells can be interpreted as a kind of h-refinement whereas the increase in the number of

417

Fig. 3 Composed integration of hexahedral elements based on subcells

Gaussian points n gp corresponds to a p-refinement. It should be noted that this h- or p-refinement does not introduce any additional degrees of freedom since it is related to the integration of the stiffness matrix only. Applying the concept of sub-cells to Eq. 15 results in n sc     BcT (ξ (r)) C e (x) Bc (ξ (r)) detJ c Kc = sc=1 t det J˜ sc c

s

r

dr dsdt

(16)

T where the determinant of the Jacobian matrix J˜ sc c = grad sc ˜ c (r, s, t) is introduced due to the additional coordinate sysQ tem r = [r, s, t]T being related to the sub-cell sc. Since the sub-cells are block-shaped, the additional mapping function ˜ sc ξ =Q c (r, s, t) is linear and therefore the Jacobian matrix ˜J sc c is constant. A standard Gaussian quadrature is performed on each of the sub-cells. The refinement of cell c into n sc sub-cells can be carried out in different ways. One option would be to apply spacetrees such as an octree following the geometry of the boundary or inclusion intersecting the cell. However, we opted for a simplified version of the refinement, where each cell c with a discontinuous integrand is refined with an increasing number of sub-cells. To this end, cells that are identified as having a discontinuous integrand are integrated with n sc = n × n × n sub-cells, where n = 1, 2, 3, . . . is the number of sub-cells in each direction being increased until the change of the value of the quadrature falls below a user-defined tolerance. In order to obtain an efficient adaptive integration algorithm, during the refinement n = 1, 2, 3, . . . instead of the stiffness matrix a scalar-valued quantity such as the volume of the element is computed and checked for convergence. The corresponding additional numerical effort to determine the necessary refinement with sub-cells is negligible compared to the integration of the stiffness matrix. However, it makes it possible

123

418

to automatically adjust the effort needed to compute the element matrices in an error-controlled fashion to produce a reliable approximation. Preliminary studies, where the adaptive sub-cell refinement has been compared to other adaptive quadrature schemes, demonstrate that the proposed method is an attractive option from a computationally point of view. A comparison of different adaptive integration schemes applied in the FCM will be presented in a forthcoming paper. In Eq. 6 the elasticity matrix Ce corresponding to the properties of the matrix material has been introduced together with α0 . As mentioned in the previous section, the weak form defined on the embedding domain Ωe equals the one defined on the physical domain Ω provided that α0 = 0. When applying the FCM it is quite likely that the cells are filled with a small fraction of material only. In such a situation, the resulting linear equation system will suffer from a poor condition number. In order to be able to solve the resulting equation system it is advantageous to choose a small value for α0 = 10−q with q in the range q = 8, 9, . . . , 12 in order to prevent ill-conditioning of the overall equation system. Of course, by choosing α0 other than 0, we introduce a model error, since we are solving a different problem. However, in all our computations we set α0 = 10−12 which is very small and so its influence on the accuracy of the solution is negligible. This observation is verified in more detailed investigations of the influence of α0 in our previous work on the FCM, see [9,39]. 4 Numerical homogenization procedure for linear elastic material behavior

Comput Mech (2012) 50:413–431

∼σm  = ∼C eff ∼εm  where

(17)



1 · = Vm

(·) d Vm

(18)

Vm

represents the volume average of the microstructural quantities. In Eq. 17 ∼εm denotes the microstructural small strain tensor and ∼σm is the corresponding microstructural stress tensor. The Dirichlet boundary value problem for the microstructure is defined by applying the projection rule

∼u¯ = ∼ε M X ∼

(19)

where ∼ε M denotes the macroscopic small strain tensor, X ∼ corresponds to the branch vector of the microstructure, see Fig. 4, and ∼u¯ represents the displacements to be applied at the boundary of the microstructure. Once the Dirichlet boundary values for the microstructure have been defined, we can proceed to solving the corresponding boundary value problem resulting in the microstructural strains and stresses. The homogenization procedure has to fulfill the Hill- Mandel condition ∼σm · ∼εm  = ∼σm  · ∼εm  = ∼σ M · ∼ε M

(20)

which enforces the equivalence of the strain energy in the macroscopic and microstructural setting. The Hill- Mandel condition is used to derive the homogenization rule, resulting in the volume average of the microstructural stresses  1 ∼σm  = (∼t ⊗ X ) dΓm (21) ∼ Vm

The numerical homogenization procedure presented in this section is based on the computation of effective material properties of RVEs related to the microstructure of the heterogeneous material under consideration. To this end, the microstructure is subjected to boundary conditions resulting in a Dirichlet boundary value problem which will be discretized by the FCM, as described in Sect. 3. This step is referred to as projection. The computation of the Dirichlet boundary value problem results in microstructural strains and stresses, from which it is possible to calculate volume averaged values by means of homogenization. This enables us to establish a relation between the volume average of microstructural strains and stresses which will be utilized to compute the effective material properties using a numerical differentiation procedure.

where rm represents the nodal forces evaluated at the n nodes on the surface of the corresponding mesh, see Fig. 4. The volume average of the microstructural strains

4.1 Projection and homogenization rule

∼ε m  =

The effective material properties are stated in terms of the fourth-order effective elasticity tensor ∼C eff which maps the volume average of the microstructural strains to the volume average of the microstructural stresses

is of interest for the extraction of the effective material properties. Provided the material under consideration is perfectly bonded, the volume integral in (23) can be transformed into a surface integral

123

Γm

which are computed from the surface integral of the dyadic product of the traction vector ∼t with the branch vector X . In ∼ the discrete form, where the microstructure is discretized by means of the FCM, the corresponding homogenization rule reads σ m  =

n 1  (i) rm ⊗ X(i) Vm

(22)

i=1

(i)

1 (Grad∼u m  + Grad∼u m T ) 2

(23)

Comput Mech (2012) 50:413–431

419

Fig. 6 Two-dimensional illustration of the window method

Fig. 4 Three-dimensional microstructure with boundary nodes

4.2 Computation of effective material parameters

(i)

Fig. 5 Equidistant surface mesh with gray-shaded areas Am

1 ∼ε m  = 2 Vm

 ( ∼u m ⊗ N +N ⊗ ∼u m ) dΓm ∼m ∼m

(24)

Γm

is the outer unit normal vector of the microstrucwhere N ∼m ture, see Fig. 4. In the discrete form the resulting volume average of the microstructural strains reads: εm  =

n   1   (i) (i) (i) (i) um ⊗ N(i) A (25) + N ⊗ u m m m m 2 Vm i=1

(i)

In Eq. 25 Am denotes the area associated with node (i) on the surface of the FCM mesh. Figure 5 depicts an equidistant FCM mesh with cell size h. We distinguish three different cases: interior nodes, edge nodes and corner nodes. In order to overcome problems with a unique definition of the unit normal vector, edge nodes are considered twice and corner nodes three times together with their associated areas and corresponding unit normal vectors.

The microstructure needs to be discretized in order to compute the effective material parameters. To this end, we apply the FCM as described in Sect. 3. The geometry of the microstructure might be defined either by a computer-aided design of a composite material or a three-dimensional computertomography of the specimen of interest. To describe the approach for computing the effective material properties, let us assume that a CT-scan of the material is available which corresponds to the domain ΩCT , as depicted for a two-dimensional situation on the left-hand side of Fig. 6. The RVE is now defined by extracting a subdomain Ωsub of ΩCT . The straight-forward application of the displacements defined by the projection rule (19) onto the boundary of the domain Ωsub would result in a response that overestimates the stiffness of the material, see Zhodi and Wriggers [52]. There are different possibilities of how to apply the boundary conditions on the microstructure and no attempt is being made here to give an overview which can be found in [52], for example. Here we follow Hain and Wriggers’ [18] idea of embedding the domain of the microstructure Ωsub into a so-called window Ωwin with thickness twin , see Fig. 6. This corresponds to the application of a self-consistency method, where the matrix material, into which particles have been incorporated, represents the effective properties. Since the aim of the homogenization procedure is to compute the effective properties, they are, of course, not known to begin with, so the window method defines a nonlinear problem which can be solved iteratively. The bulk value of the underlying microstructural material is a good starting-point for the window properties. The effective parameters of the RVE within the iterative process, i.e. the entries of the material matrix Ceff in Voigt notation are calculated using numerical differentiation. To this end, we choose an arbitrary macroscopic basic state of strain and modify each component at a time using the step size ε = 0.1 Since the symmetric strain tensor has six indi1

Note that index m and M will be skipped in the remainder of this section in order to avoid overloading the notation.

123

420

Comput Mech (2012) 50:413–431

vidual components in the three-dimensional version, we need to solve seven load cases, i.e. one load case corresponding to the basic state and six load cases to modify each strain component by the step size ε. Using the Voigt notation allows us to arrange the strains related to the different load cases  (26) ε˜ = ε˜ 1 . . . ε˜ 7 . We use the projection rule (19) to apply the strains (26) to the window’s boundary Γwin , solve the resulting Dirichlet boundary value problems and, using (22) and (25) respectively, obtain stresses  (27) σ˜  = σ˜1  . . . σ˜7  and strains  ˜ε  = ε˜1  . . . ε˜7 

(28)

for the subdomain Ωsub . The components of C eff can then be computed by ⎧ (σ˜ i −σ˜ i ) ⎪ if j <= 3 ⎨ ε˜ 1−ε˜ j+1

1 j+1 L eff 2 , (29) Ci j ≈ i ) (σ˜ 1i −σ˜ j+1 ⎪ ⎩1 if j > 3 2 ε˜ 1 −ε˜ j+1 

L2 where superscripts denote the individual vector-components of the (˜· )-quantities and the factor 1/2 is due to the Voigt notation. The explicit formulation of the effective elasticity matrix (29) computed by the difference quotient reads ⎤ ⎡ 12 σ 13 σ 14 σ 15 σ 16 σ 17

σ xx xx xx xx xx xx ⎥ ⎢ ⎥ ⎢ ε 12 ε 13 ε 14 2 ε 15 2 ε 16 2 ε 17 ⎥ ⎢ ⎥ ⎢ 13 σ 14 σ 15 σ 16 σ 17 ⎥ ⎢

σ yy yy yy yy yy ⎥ ⎢ ⎥ ⎢

ε 13 ε 14 2 ε 15 2 ε 16 2 ε 17 ⎥ ⎢ ⎥ ⎢ ⎢ 14 σ 15 σ 16 σ 17 ⎥

σzz ⎢ zz zz zz ⎥ ⎥ ⎢ ⎢

ε 14 2 ε 15 2 ε 16 2 ε 17 ⎥ C eff ≈ ⎢ ⎥ ⎢

σx15y σx16y σx17y ⎥ ⎥ ⎢ ⎥ ⎢ ⎢ 2

ε 2

ε 2

ε 16 17 ⎥ 15 ⎥ ⎢ ⎢ 16 σ 17 ⎥

σ yz ⎢ yz ⎥ ⎥ ⎢ ⎢ 2 ε 16 2 ε 17 ⎥ ⎥ ⎢ ⎥ ⎢ ⎣

σx17z ⎦ sym 2 ε 17 with 1j

j

σkl = (σkl1 − σkl ),  

ε 1 j = ˜ε 1  − ˜ε j 

L2

5 Numerical examples The numerical examples we presented in Sect. 5.1 serve to verify the homogenization procedure by considering a porous material, a material with inclusions and a fiber-reinforced material for which we compute effective properties. Having obtained the effective properties, we proceed to apply them via a finite element computation and compare the results to those of FEM and FCM computations, which fully resolve the microstructure of the material of interest. In Sect. 5.2 we apply the homogenization procedure to cellular materials, for which a computer-tomography has been carried out. We consider a metal foam consisting of aluminum and compare the results of the homogenization procedure to reference values obtained from literature. 5.1 Verification of the proposed method The proposed homogenization procedure is verified by considering heterogeneous materials characterized by the volume fraction ϕV =

Vparticles/fibers · 100%, Vtotal

(30)

which can be defined in a similar manner for cellular materials by substituting Vparticles/fibers by Vcellwalls . 5.1.1 Definition of the examples

,

where · L2 denotes the L2 -norm. Once the effective material properties have been obtained, they are applied to the domain of the window Ωwin and the computation of the effective

123

properties is repeated until the Frobenius norm of the difference between two consecutively computed effective material matrices falls within a prescribed tolerance. For the efficiency of the proposed scheme it is important to note that the stiffness matrices of the finite cells applied to discretize the domain Ωsub do not change during the iteration and only have to be calculated once. Since we assume the effective material properties for the domain Ωwin , the corresponding stiffness matrices need to be recalculated every iteration step, so the overall equation system is modified each time. If we are to solve the equation system in an efficient manner, however, we should take advantage of the fact that we need to solve a system with multiple right-hand sides, i.e. seven load cases. We accordingly apply the Parallel Sparse Direct Solver Pardiso [42,41] which factorizes the overall stiffness matrix at each iteration step once only and performs a backward substitution for each of the seven load cases, significantly reducing the overall computational time.

The geometry of the microstructure of the heterogeneous materials is represented by means of an implicit description of the geometry [9]. The dimensions of the domain are l × l × l where l = 10 mm. In the first example, the heteroge-

Comput Mech (2012) 50:413–431

421

neous material to be homogenized consists of a matrix with randomly distributed ellipsoids, see Fig. 7. In the case of the porous material, the ellipsoids represent voids and in the case of inclusions they correspond to a material which is stiffer than the matrix material. In both cases, the volume fraction amounts to ϕV = 3%. The second example is a unidirectional fiber-reinforced composite with a volume fraction of , see ϕV = 18% where the fibers are aligned parallel to E ∼X Fig. 8. In both cases, the linear elastic material parameters of the matrix are Young’s modulus E matrix = 5000 MPa and Poisson’s ratio νmatrix = 0.3. With regard to the ellipsoids, we distinguish between two separate cases: – pores (no material), i.e. E pores = 0.0 MPa and ν pores = 0.0. – stiff material, i.e. E inc = 10 E matrix = 50 GPa and ν inc = 0.3.

Fig. 7 Material with randomly distributed pores

For the fibers we set E fiber = E inc and ν fiber = ν inc . In Fig. 9 the boundary of the physical domain is labelled in order to define the boundary conditions. Two different load cases are defined: uniaxial compression where u¯ 1z = 0, u¯ 4y = 0, u¯ 5x = 0, t¯z6 = −100, and general three-dimensional loading due to body forces f where u¯ 1z = 0, u¯ 4y = 0, u¯ 5x = 0, f x = 300, f y = 200, f z = 100. The numerical examples are solved by applying different approaches:

Fig. 8 Unidirectional fiber-reinforced composite

– with effective materials properties. To this end, Ceff is extracted first and then applied to a macroscopic computation using the classical p-FEM. – directly, i.e. resolving the microstructure – by the FCM and – by classical h-FEM using the commercial FE-Solver Marc.2 5.1.2 Effective material parameters We start with the extraction of effective material parameters. The domain Ωsub , see Fig. 6, is resolved by a grid of c3 cells, whereby c = 8 corresponds to the coarse grid and c = 16 corresponds to the fine one. For the thickness of the window, we choose twin = {0.1, 1, 2, 3}. The polynomial degree of all elements and cells varies from p = 1, . . . , 8 for the coarse grid and from p = 1, . . . , 6 for the fine grid. The trunk space [10,48] is applied in all of the computations. For an accurate 2

Marc is a trademark of MSC Software Corporation, 2 MacArthur Place, Santa Ana, California 92707, US.

Fig. 9 Each surface of the different heterogeneous materials depicted in Figs. 7 and 8 is labelled

integration the number of Gaussian points in one direction is set to 20 for all cases and the finite cells are subdivided automatically into sub-cells, until the error in the volume computed for each cell falls below 0.1%, see Fig. 3. We will now proceed to computing the effective material matrices of the three microstructures as introduced in Sect. 5.1.1, using the proposed homogenization procedure.

123

422

Comput Mech (2012) 50:413–431

pores

versus twin for c = 8 and p = 8

Fig. 12 Cii

pores

versus N for c = 8 and twin = 2 mm

Fig. 13 Ciiinc versus N for c = 8 and twin = 2 mm

Fig. 10 Cii

Fig. 11 Cii

In order to distinguish between the effective material matrices of the different microstructures, they are defined as Cpores , Cinc , and Cfiber . Please note, that the superscript ‘eff’ has been omitted for the sake of simplifying the notation. Before we investigate these matrices, let us begin by taking a look at the influence of the window size. We consider the porous material, which is discretized with c = 8 cells in each direction using a polynomial degree of p = 8 for all cells. In Fig. 10 the coefficients Cii , i = 1, . . . , 3 of the effective material matrix are plotted against the thickness of the window. From this it is evident that the influence of the thickness of the window is moderate when twin ≥ 2 mm. For the remainder of this section we accordingly set twin = 2 mm. Next we study the influence of the discretization where the polynomial degree is successively increased for the two different grids with c = 8 and c = 16. In Figs. 11, 12, 13, 14, 15, 16 the coefficients Cii , i = 1, . . . , 3 are plotted against the number of degrees freedom N for each of the three different heterogeneous material types. A separate figure is provided for each grid, the coarse one (c = 8) and the fine one (c = 16), where a hierarchic extension of the Ansatz is performed by increasing the polynomial degree p.

123

pores

versus N for c = 16 and twin = 2 mm

Fig. 14 Ciiinc versus N for c = 16 and twin = 2 mm

If we look at the convergence of the hierarchic extension, we can conclude that—even with the coarse FCM grid—it is possible to obtain a reasonable degree of accuracy for the effective material properties. In order to study the convergence of the iterative procedure to determine the effective material properties, we

Comput Mech (2012) 50:413–431

Fig. 15 Ciifiber versus N for c = 8 and twin = 2 mm

Fig. 16 Ciifiber versus N for c = 16 and twin = 2 mm

consider C fiber with c = 8 by way of an example.3 The convergence of the iteration is investigated by computing the difference in the effective material matrix between two consecutive iteration steps  new  C − C old Frob   . (31) eFrob = C old  Frob

In the above equation, · Frob denotes the Frobenius norm which is computed for the effective material matrix and plotted in Fig. 17 against the number of iterations. In addition we investigate the convergence of the coefficients Ciifiber , i = 1, . . . , 3 with c = 8 and p = 8 normalized by their corresponding converged values, see Fig. 18. From Figs. 17 and 18 it is evident that the iterative procedure converges quite fast so that accurate results can be obtained for engineering purposes after just 3-4 iterations. Finally, let us consider the effective elasticity matrices obtained by the homogenization process based on the FCM grid with c = 16 cells in each direction, utilizing a polynomial degree of p = 6. The effective material matrices of the 3

Please note that a similar behavior is evident in the other examples, too.

423

Fig. 17 eFrob versus n iter , again c = 8, p = 8 and twin = 2 mm

Fig. 18 Normalized entries of Ciiinc versus n iter , again c = 8, p = 8 and twin = 2 mm

different heterogeneous materials with the pores, inclusions and fiber reinforcement read: C pores ⎡

6.28E3 2.63E3 2.62E3 −1.19E0 6.27E3 2.62E3 −1.44E0 ⎢ ⎢ 6.24E3 −4.55E−1 ⎢ =⎢ 1.81E3

⎣

⎤

−2.76E−1 −1.80E0 −2.24E0 −5.18E−1 1.80E3

−3.06E0 −1.35E0 ⎥ −2.88E0 ⎥ ⎥, 3.32E−1 ⎥ ⎦ −2.13E−1 1.81E3

8.26E−1 −1.08E0 −6.33E−1 −4.46E−2 2.03E3

−7.21E−1 1.78E−1 ⎥ −6.49E−1 ⎥ ⎥, 5.82E−1 ⎥ ⎦ −6.79E−3 2.04E3

sym inc

C ⎡

7.13E3 3.01E3 3.01E3 −4.78E−1 7.12E3 3.01E3 −1.02E0 ⎢ ⎢ 7.09E3 7.30E−2 ⎢ =⎢ 2.04E3

⎣

sym

⎤

and C fiber ⎡

1.39E4 3.79E3 3.79E3 8.88E3 3.64E3 ⎢ ⎢ 8.87E3 ⎢

=⎢ ⎣

sym

⎤

0 −4.53E0 0 0 −9.61E−1 0 ⎥ ⎥ 0 −1.94E1 0 ⎥. 2.72E3 0 −3.88E0 ⎥ ⎦ 2.52E3 0 2.71E3

123

424

Comput Mech (2012) 50:413–431

An anisotropic material behavior can be observed in all of the three cases. However, this anisotropic behavior is rather small for the porous material and the material with inclusions pores when comparing the absolute values like, for example, C13 pores inc with C15 or Cinc 13 with C14 . With regard to the effective elasticity matrix Cfiber we observe that the 11-direction is more pronounced, due to the fiber reinforcement, of course. 5.1.3 Verification of the effective material properties The aim of this section is to verify the effective material properties by comparing them to two different approaches, i.e. the FCM and FEM, where the microstructure of the material is fully resolved. Let us recall the different solution techniques to be applied in more detail: – The effective material properties are utilized within a finite element approximation which is based on a mesh consisting of (2 × 2 × 2) = 8 hexahedral elements applying hierarchic shape functions. This discretization leads to a very low number of degrees of freedom, since the microstructure of the heterogeneous material has already been taken into account by means of the effective material properties, and the material is now treated as homogeneous. – For the FCM approach, accounting for the microstructure during the integration of the stiffness matrix computation, the grid, the adaptive numerical quadrature and the Ansatz space correspond to the setting described in Sect. 5.1.2. This computation does not take advantage of the effective material properties, since the geometry of the microstructure is resolved using the FCM. – The FEM approach likewise resolves the microstructure of the material in full. The discretization is based on 10-noded fully integrated tetrahedral elements and 15-noded fully integrated pentahedrals. The number of elements applied to discretize the porous material varies from 27,582 to 222,719 tetrahedrals. With regard to the material with inclusions, we have between 35,235 and 311,464 tetrahedrals and, in the case of the fiber-reinforced material, 5,824–90,916 pentahedrals are applied. Figures 19, 20, 21, 22 show the coarsest and the finest mesh. Please note that the meshes for the porous domain are not displayed, since they are similar to the ones applied in the case of the material with inclusions. To determine the quality of the different approaches the  1 strain energy U = 2 Ω ∼σ · ∼ε dΩ is plotted against the number of degrees of freedom N in Figs. 23, 24, 25, 26, 27, 28. A very good agreement of the strain energy between the different solution approaches—especially for load case 1— can be observed for the porous material, see Figs. 23 and

123

Fig. 19 h-version: coarsest mesh for the material with inclusions

Fig. 20 h-version: finest mesh for the material with inclusions

Fig. 21 h-version: coarsest mesh for the fiber-reinforced material

Comput Mech (2012) 50:413–431

425 960 950 940 930 920 910 900 890 880 870 860 850

0

500000

1e+06

Fig. 25 Load case 1: strain energy in case of inclusions Fig. 22 h-version: finest mesh for the fiber-reinforced material 520000

1070 510000

1065 1060

500000

1055

490000

1050

480000

1045

470000

1040 460000

1035 450000

1030

0

500000

0

500000

1e+06

1e+06

Fig. 26 Load case 2: strain energy in the case of inclusions Fig. 23 Load case 1: strain energy for the porous material 540000

760 740 720

530000

700 680

520000

660 640 620

510000

600 580

500000

560 540

0

500000

1e+06

Fig. 24 Load case 2: strain energy for the porous material

0

500000

1e+06

Fig. 27 Load case 1: strain energy in case of fiber reinforcement

123

426

Comput Mech (2012) 50:413–431

280000

for efficiently extracting the effective material properties of heterogeneous materials with an accuracy that is definitely acceptable for engineering applications.

270000

5.2 Effective material parameters of cellular materials

260000

In this section we apply the proposed homogenization procedure to determine the effective material properties of cellular materials. The starting point is a CT-scan of the cellular material to be discretized using FCM. At the end of this section we will summarize the results and suggest how to choose the different parameters in order to extract the effective material parameters of cellular materials. Before presenting the numerical results of the homogenization process, let us first briefly recall related approaches familiar from the literature to describe the material behavior of cellular materials. Huber and Gibson described the material behavior of foams by idealized three-dimensional unit cells and validated their results by experiments [21]. Later, Gong et al. published detailed geometric information about PU foams and established models based on Kelvin cells, [15,16]. Their approach was improved by Jang et al. who used 3D CT-scans of polymeric and aluminum foams as a basis for FE-models, [22]. These models are based on different types of cells, including Kelvin cells. A similar investigation was carried out for polymethacrylimide (PMI) foams by Wang et al., who utilized Kelvin cells [49]. The drawback of the aforementioned approaches is that the foam’s real geometry is only approximated. In order to overcome these problems, Wismans et al. converted three-dimensional micro CT-scans of polymeric foams into FE-models in order to study their elastic and hyperelastic properties [51]. However, their approach was restricted to two dimensions, ignoring the importance of the three-dimensional situation. In contrast to [51], we discretize the real 3D geometry which allows for extracting fully three-dimensional linear elastic effective material properties.

290000

250000 240000 230000

0

500000

1e+06

Fig. 28 Load case 2: strain energy in case of fiber reinforcement

24. The strain energy remains almost constant during the p-extension when using the effective material properties. The reason for this is twofold: On the one hand, the effective material matrix C pores is constant within the homogenized material while the anisotropy of the material is not particularly strong. As regards load case 2, only a very small deviation of the solution is observed, based on the effective material parameters as compared with the other solution approaches. This is due to the fact that the overall volume load can only approximately be taken into account during the discretization of the homogenized material. In this case, the volume load is scaled by the volume fraction, i.e. multiplied by a factor of 0.97. In the example with inclusions, the pure FCM and the approach with effective material parameters lead to similar results but seem to deviate slightly from the results obtained with the h-version, see Figs. 25 and 26. This deviation stems from the fact that the material interfaces that lead to kinks in the displacement field cannot be exactly represented by the FCM but are approximated with smooth polynomials, thus introducing a certain discretization error. It is, of course, possible to adjust the cells in such a way that they match the interfaces. However, this approach would somewhat contradict the simplicity of the FCM, which is based on fast, simple grid generation. Taking a look at the absolute values, we see that the deviation amounts to just 1.3% for load case 2. Since the computation of effective material properties is based on volume-averaged quantities, the FCM still yields acceptable results, since the averaging process offsets local effects like those introduced by the material interfaces. For the same reason, there are also some differences between the three approaches in the case of the fiber-reinforced material, see Figs. 27 and 28. Taking the h-version approximation as a reference solution, the maximum error of the approximation based on the effective material properties is about 2.7% when considering load case 2. Summarizing the results of this section, we conclude that the FCM approach is ideal

123

5.2.1 Material under consideration In our investigation, we focus on an equidistant CT-scan of an open-cell aluminum foam with 10 pores per inch (ppi) and a volume fraction, as defined in Eq. 30 of ϕV ≈ 9%. The CT-scan consists of 733 × 729 × 704 voxels where each voxel has the volume of V = 60.331×60.331×60.331 µm3 providing a proper resolution of the individual cell walls, see Fig. 29. In the remainder of this section, all geometrical quantities are expressed in terms of voxels, i.e. we choose a voxel as a unit length in order to derive general guidelines on how to discretize CT-scans by applying the FCM. The mechanical parameters for the cell walls within a linear elastic analysis are assumed to correspond to a Young’s modulus of E = 70 GPa and a Poisson’s ratio of ν = 0.3.

Comput Mech (2012) 50:413–431

427 55 50 45 40 35 30 25

Fig. 29 Resolution of the individual cell walls by the CT-scan 20 1000

10000

100000

1e+06

1e+07

5.2.2 Discretization The goal of this section is to determine a suitable discretization of the foam which can be used for computing its effective material parameters. To this end, we choose a subdomain of size n x × n y × n z = 200 × 200 × 200 = l 3 = 2003 which has its anchor point at position (x, y, z) = (100, 100, 100). Please note that approximately 18 foam cells are included in this subdomain. To determine a reasonable discretization the subdomain is subjected to Dirichlet boundary conditions, i.e. it is clamped at the bottom and loaded by a prescribed displacement of u¯ z = −1.0 at the top. Next, a p -extension with three different finite cell resolutions with 103 , 203 , and 403 cells is performed, resulting in 20, 10, and 5 voxels per finite cell. In the adaptive numerical integration of the finite cells, the number of Gaussian points in one direction is set to 20 and the cell is adaptively subdivided into sub-cells until the change of the computed volume of two consecutive refinements falls below a prescribed tolerance of 0.1%. In Fig. 30 the resulting strain energy U is plotted against the number of degrees of freedom N of a hierarchic p-extension. From the convergence of the strain energy we conclude that a resolution of 10 voxels per finite cell with a polynomial degree of p ≥ 5 yields an acceptable approximation. This kind of discretization with p = 5 results in N = 295,335 degrees of freedom and the corresponding strain energy U ≈ 26.31 deviates by just 1.3% from the finest FCM grid in which we have five voxels per finite cell and p = 6. Based on this preliminary study, we chose 10 voxels per finite cell for the spatial discretization and set p = 5 for the remaining examples. 5.2.3 Window thickness Next, we embed our subdomain into windows with a varying thickness of twin = {1, 5, 10, 20, 40}. The influence of twin onto Cii , i = 1, . . . , 6 is depicted in Fig. 31. As mentioned in the introduction cellular materials exhibit size-dependent boundary layer effects under various load conditions such as shear, tension and bending. It should be noted that shear

Fig. 30 Strain energy U versus number of degrees of freedom N for different FCM resolutions 1800 1600 1400 1200 1000 800 600 400 200 0 0

5

10

15

20

25

30

35

40

Fig. 31 Cii versus twin for the subdomain with 200×200×200 voxels starting at position (100,100,100)

experiments in particular exhibit a strong size dependency. In a two-scale approach, such effects lead to an overestimation of the mechanical properties even when large subdomains (microstructures) are employed. For this reason, we propose using such means as introducing pin supports at the microstructural boundaries or applying periodic boundary conditions in order to relax such local stiffening effects [46]. The window method has a similar effect and there is evidence that, for a constant number of foam cells, the thicker the window the less pronounced the stiffening effect. In the practical application of the resulting effective elasticity matrices to macroscopic problems, however, the distance to the boundary is not known in advance, nor is it constant within the macroscopic domain. It is therefore necessary to find a compromise with respect to the choice of twin . For the foam at hand we choose twin = 10 = l/20. Another important observation is that the z-direction of the foam is more pronounced, indicating anisotropic behavior.

123

428

Comput Mech (2012) 50:413–431

1200

1000

800

600

400

200

8748 12000

26364

40500

49152

Fig. 33 C¯ ii versus Ωsub . The error bars indicate the standard deviation sii as defined in (33)

Fig. 32 Illustration of expansion direction and cell anisotropy

This is not surprising since the manufacturing process of foam materials involves an upward expansion in which the cells are elongated vertically by means of aeration [21,22]. Perpendicular to the vertical direction the cells are nearly equiaxed and we can define A = h 1 / h 2 where h 1 is measured along the vertical axis. For the foam concerned we have A ≈ 1.5 and, as indicated above, the upward expan, see Fig. 32. For obvious sion direction coincides with E ∼Z reasons it is desirable that the subdomain, i.e. the RVE of the foam should have the same number of foam cells in all of the three spatial directions. So A has to be taken into account for generating the subdomains. If we denote the subdomain’s size in the x y-plane by l x y then the expanded thickness emerges as l z = A · l x y and the subdomain’s volume is Ωsub = l x2y · l z = A · l x3y . For the window’s thickness we choose l x y /20 or l z /20 respectively. 5.2.4 Representative subdomain

(32)

k=1

is plotted against Ωsub including error bars indicating the standard deviation

123

   si j = 

n 2 1  ¯ Ci j − Cikj . n−1

(33)

k=1

The next step is to determine a representative subdomain, i.e. RVE. To this end, the subdomain’s size and the starting positions are varied for statistical reasons. In particular, we set l x y = {180, 200, 260, 300, 320} and for each size l x y n = 10 samples with randomly distributed anchor points are determined. In Fig. 33 the arithmetic mean value n 1 k ¯ Ci j Ci j = n

Fig. 34 Two closed faces marked by circles

We notice that the mean values increase slightly in line with the increasing size of the subdomain. This phenomenon can be explained by the higher probability of local stiffening effects such as closed faces [22], see Fig. 34. Almost converged values for C¯ ii are obtained from Ωsub = 40, 500 × 103 voxels on resulting in N ≈ 1.9 × 106 degrees of freedom and corresponding to 4×4×4 foam cells. It should be noted that in [49] a similar number of cells was identified for PMI foams. Figure 35 depicts the corresponding subdomain of the foam together with the FCM mesh consisting of 15153 cells. The computational time needed to calculate the effective material properties on a parallel computer with 8 CPUs with

Comput Mech (2012) 50:413–431

429

E plane =

2 · C + 2 · C2 · C − 2 · C2 · C − C · C2 ) (C11 33 12 11 33 13 13 12 2) (C11 · C33 − C13

.

(35) In our example we have E trans/rise ≈ 725 MPa and E plane ≈ 300 MPa which is in the same order of magnitude as the measurements reported by Jang [22] for an aluminum foam with J ≈ 585 MPa 10 ppi, A ≈ 1.3, and ϕV ≈ 8.7% with E trans/rise J ≈ 363 MPa. and E plane 5.2.5 Summary Fig. 35 Foam related to Ωsub = 40, 500 × 103 and corresponding FCM mesh

6 cores each amounts to approximately 5 h. A total of 9 iteration steps can be conducted during this time-frame to meet the tolerance of eFrob < 10−3 . We wish to point out that this tolerance is quite low and it would be possible to obtain reasonable results with even fewer iteration steps, resulting in a shorter computation time. We suggest contemplating the mean effective elasticity matrix for the purpose of further discussions. To this end, we present the one associated with 40, 500 × 103 voxels ¯ 40,500×103 C ⎡

367.606 120.665 221.621 −0.965 371.579 232.999 −0.588 ⎢ ⎢ 926.978 2.311 ⎢ =⎢ 117.069

⎣

6.136 10.577 24.700 3.565 155.743

⎤

23.670 16.468 ⎥ 42.073 ⎥ ⎥. 2.402 ⎥ ⎦ 3.622 154.299

¯ 40500×103 we can find some Investigating the structure of C similarity to an elasticity matrix of a transverse isotropic material ⎡ ⎤ C11 C12 C13 0 0 0 ⎢ 0 0 0 ⎥ C11 C13 ⎢ ⎥ ⎢ 0 0 0 ⎥ C33 ⎢ ⎥ ⎥, Ctrans−iso = ⎢ (C11 − C12 ) ⎢ 0 0 ⎥ ⎢ ⎥ 2 ⎢ ⎥ ⎣ C44 0 ⎦ sym C44 which is due to the manufacturing process of the foam as mentioned above. Referring to [4] Young’s modulus can be computed for the transverse/expansion direction by E trans/rise =

2 · C + 2 · C2 · C − 2 · C2 · C − C · C2 ) (C11 33 12 11 33 13 13 12 2 − C2 ) (C11 12

(34) and for the plane of isotropy by

We will now proceed to summarize the most important results of the numerical investigations of the foam by way of guidelines for computing effective material properties of foams using the proposed homogenization strategy. – Provided a proper resolution of the individual cell walls by the CT-scan, 10 voxels per finite cell and p = 5 yield accurate results suitable for the computation of the effective material properties. – A representative subdomain should include approximately 4 × 4 × 4 foam cells. – It is known that the window’s size influences the computation of the effective properties [17,18]. For foamed materials this is even more pronounced due to size effects which may lead to stiffening. We therefore recommend investigating the influence of the window size for the foam in question by means of numerical studies. It might be a good idea to compute effective material properties for different values of the window size and to apply the resulting elasticity matrices depending on the distance to the boundary of the macroscopic problem.

6 Conclusions We have presented a numerical homogenization procedure which can be employed to compute the effective material properties of cellular and heterogeneous materials. The key feature of the method is the FCM, which makes it possible to discretize complicated microstructures in a fully automatic way, where the error control is based on a hierarchic extension of the Ansatz space combined with an adaptive quadrature of the cell matrices. Special attention is also given to the way the boundary conditions are applied to the microstructure during the projection step of the homogenization stage. Several numerical examples are presented and demonstrate the performance of the proposed method for materials with pores, inclusions/fibers as well as for cellular materials. Based on the numerical examples, it is evident that

123

430

the proposed method is a very efficient one which produces reliable results, even for complex microstructured materials. Future work will investigate the application of the homogenization procedure to the computation sandwich plates and shells, where two faceplates cover a core of foamed material, for which the effective material properties need to be determined. A further field of application involves composite materials including the modeling and simulation of damage, which is a challenging task. Acknowledgments The support provided by the Deutsche Forschungsgemeinschaft (DFG) is gratefully acknowledged (DU 405/22, RA 624/16-2).

References 1. Aboudi J (1991) Mechanics of composite materials. Elsevier, Amsterdam 2. Anderson WB, Chen CP, Lakes RS (1994) Experimental study of size effects and surface damage in closed-cell polymethacrylimide and open-cell copper foams. Cellular Polym 13(1):1–15 3. Barbero E (2007) Finite element analysis of composite materials. CRC Press, Boca Raton 4. Bower A (2009) Applied mechanics of solids. CRC Press, Boca Raton 5. Brezny R, Green DJ (1990) Characterization of edge effects in cellular materials. J Mater Sci 25:4571–4578 6. Chen C, Fleck NA (2002) Size effects in the constrained deformation of metallic foams. J Mech Phys Solids 50:955–977 7. Del Pino S, Pironneau O (2003) A fictitious domain based general pde solver. In: Kuznetsov PNY, Pironneau O (eds) Numerical methods for scientific computing variational problems and applications. CIMNE, Barcelona 8. Düster A, Bröker H, Heidkamp H, Heißerer U, Kollmannsberger S, Wassouf Z, Krause R, Muthler A, Niggl A, Nübel V, Rücker M, Scholz D (2004) AdhoC 4 —User’s Guide. Lehrstuhl für Bauinformatik, Technische Universität München 9. Düster A, Parvizian J, Yang Z, Rank E (2008) The finite cell method for three-dimensional problems of solid mechanics. Comput Methods Appl Mech Eng 197:3768–3782 10. Düster A, Bröker H, Rank E (2001) The p-version of the finite element method for three-dimensional curved thin walled structures. Int J Numer Methods Eng 52:673–703 11. Feyel F (2003) A multilevel finite element method (FE2 ) to describe the response of highly non-linear structures using generalized continua. Comput Methods Appl Mech Eng 192:3233–3244 12. Feyel F, Chaboche JL (2000) FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fiber SiC/Ti composite materials. Comput Methods Appl Mech Eng 183:309–330 13. Fish J, Belsky V (1995) Multigrid method for periodic heterogeneous media. Part 2: multiscale modeling and quality control in multidimensional case. Comput Methods Appl Mech Eng 126:17– 38 14. Gilat R, Banks-Sills L (eds) (2009) Advances in mathematical modeling and experimental methods for materials and structures— The Jacob Aboudi volume. Springer, New York 15. Gong L, Kyriakides S (2005) Compressive response of open cell foams. Part II: initiation and evolution of crushing. Int J Solids Struct 42:1381–1399 16. Gong L, Kyriakides S, Jang WY (2005) Compressive response of open-cell foams. Part I: morphology and elastic properties. Int J Solids Struct 42:1355–1379

123

Comput Mech (2012) 50:413–431 17. Hain M (2007) Computational Homogenization of micro-structural damage due to frost in hardened cement paste. PhD thesis, Institut für Baumechanik und Numerische Mechanik, Leibniz-Universität Hannover 18. Hain M, Wriggers P (2008) Numerical homogenization of hardened cement paste. Comput Mech 42:197–212 19. Hirschberger CB, Sukumar N, Steinmann P (2008) Computational homogenization of material layers with micromorphic mesostructure. Philos Mag 88:3603–3631 20. Hohe J, Becker W (2003) Effective mechanical behavior of hyperelastic honeycombs and two-dimensional model foams at finite strain. Int J Mech Sci 45:891–913 21. Huber AT, Gibson LJ (1988) Anisotropy of foams. J Mater Sci 23:3031–3040 22. Jang WY, Kraynik AM, Kyriakides S (2008) On the microstructure of open-cell foams and its effect on elastic properties. Int J Solids Struct 45:1845–1875 23. Jänicke R, Diebels S, Sehlhorst HG, Düster A (2009) Two-scale modelling of micromorphic continua. Continuum Mech Thermodyn 12:297–315 24. Kari S, Berger H, Gabbert U (2007) Numerical evaluation of effective material properties of randomly distributed short cylindrical fibre composites. Comput Mater Sci 39:198–204 25. Kouznetsova V (2002) Computational homogenization for the multi-scale analysis of multi-phase materials. PhD thesis, Materials Technology, Eindhoven University of Technology 26. Kouznetsova VG, Geers MGD, Brekelmans WAM (2002) Multiscale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme. Int J Numer Methods Eng 54:1235–1260 27. Lakes RS (1983) Size effects and micromechanics of porous solids. J Mater Sci 18:2572–2580 28. Loehnert S, Wriggers P (2008) Effective behaviour of elastic heterogeneous thin structures at finite deformations. Comput Mech 41:595–606 29. Löhnert S (2005) Computational homogenization of microheterogeneous materials at finite strains including damage. Ph.D. thesis, Institut für Baumechanik und Numerische Mechanik, Universität Hannover 30. Matzenmiller A, Kurnatowski B (2009) A comparison of micromechanical models for the homogenization of microheterogeneous elastic composites. In: Gilat R, Banks-Sills L (eds) Advances in mathematical modeling and experimental methods for materials and structures—The Jacob Aboudi volume. Springer, New York 31. Miehe C, Koch A (2002) Computational micro-to-macro transitions of discretized microstructures undergoing small strains. Arch Appl Mech 72:300–317 32. Miehe C, Schröder J, Schotte J (1999) Computational homogenization analysis in finite plasticity. Simulation of texture development in polycrystalline materials. Comput Methods Appl Mech Eng 171:387–418 33. Miehe C, Schröder J, Bayreuther C (2002) On the homogenization analysis of composite materials based on discretized fluctuations on the micro–structure. Acta Mech 155:1–16 34. Miehe C, Schröder J, Becker M (2002) Computational homogenization analysis in finite elasticity: material and structural instabilities on the micro- and macro-scales of periodic composites and their interaction. Comput Methods Appl Mech Eng 191:4971– 5005 35. Neittaanmäki P, Tiba D (1995) An embedding of domains approach in free boundary problems and optimal design. SIAM J Control Optim 33(5):1587–1602 36. Nemat-Nasser S, Hori M (1999) Micromechanics: overall properties of heterogeneous materials, 2nd edn. North-Holland, Amsterdam

Comput Mech (2012) 50:413–431 37. Oden J, Zhodi T (1997) Analysis and adaptive modeling of highly heterogeneous elastic structures. Comput Methods Appl Mech Eng 148:367–391 38. Onck PR (2002) Cosserat modeling of cellular solids. C R Mecanique 330:717–722 39. Parvizian J, Düster A, Rank E (2007) Finite cell method—h- and p-extension for embedded domain problems in solid mechanics. Comput Mech 41:121–133 40. Saul’ev VK (1963) On solution of some boundary value problems on high performance computers by fictitious domain method. Sib Math J 4:912–925 41. Schenk O, Wächter A, Hagemann M (2007) Matching-based preprocessing algorithms to the solution of saddle-point problems in large-scale nonconvex interior-point optimization. Comput Optim Appl 36(2–3):321–341 42. Schenk O, Bollhöfer M, Römer RA (2008) On large scale diagonalization techniques for the Anderson model of localization. SIAM Rev 50(1):91–112 43. Schröder J (1996) Theoretische und algorithmische Konzepte zur phänomenologischen Beschreibung anisotropen Materialverhaltens. PhD thesis, Institut für Baumechanik und Numerische Mechanik, Universität Hannover 44. Schröder J (2000) Homogenisierungsmethoden der nichtlinearen Kontinuumsmechanik unter Beachtung von Stabilitätsproblemen. Postdoctoral thesis, Institut für Mechanik (Bauwesen) Lehrstuhl I, Universität Stuttgart

431 45. Schröder J (2000) Zu Stabilitätsproblemen bei Mikro–Makro– Übergängen. ZAMM-Zeitschrift für Angewandte Mathematik und Mechanik 80:411–412 46. Sehlhorst HG, Jänicke R, Düster A, Rank E, Steeb H, Diebels S (2009) Numerical investigations of foam-like materials by nested high-order finite element methods. Comput Mech 45:45–59 47. Szabó B, Babuška I (1991) Finite element analysis. Wiley, New York 48. Szabó B, Düster A, Rank E (2004) The p-version of the finite element method. In: Stein E, de Borst R, Hughes TJR (eds) Encyclopedia of computational mechanics, vol 1, chap 5. Wiley, New York pp 119–139 49. Wang J, Wang H, Xiuhua C, Yin Y (2010) Experimental and numerical study of the elastic properties of PMI. J Mater Sci 45:2688–2695 50. Warren WE, Kraynik AM (1991) The nonlinear elastic behavior of open-cell foams. J Appl Mech 58:376–381 51. Wismans JGF, Govaert LE, van Dommelen JAW (2010) X-ray computed tomography-based modeling of polymeric foams: the effect of finite element model size on the large strain response. J Polym Sci B 48:1526–1534 52. Zhodi T, Wriggers P (2005) Introduction to computational micromechanics. Lecture notes in applied and computational mechanics, vol 20. Springer, Berlin

123