A micromechanics-enhanced finite element formulation for modelling heterogeneous materials

# A micromechanics-enhanced finite element formulation for modelling heterogeneous materials

Jan Novák Łukasz Kaczmarczyk Peter Grassl Jan Zeman Chris J. Pearce School of Engineering, University of Glasgow, Glasgow G12 8QQ, UK Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7, 166 29 Praha 6, Czech Republic School of Civil and Environmental Engineering, University of New South Wales, NSW 2052, Sydney, Australia
###### Abstract

In the analysis of composite materials with heterogeneous microstructures, full resolution of the heterogeneities using classical numerical approaches can be computationally prohibitive. This paper presents a micromechanics-enhanced finite element formulation that accurately captures the mechanical behaviour of heterogeneous materials in a computationally efficient manner. The strategy exploits analytical solutions derived by Eshelby for ellipsoidal inclusions in order to determine the mechanical perturbation fields as a result of the underlying heterogeneities. Approximation functions for these perturbation fields are then incorporated into a finite element formulation to augment those of the macroscopic fields. A significant feature of this approach is that the finite element mesh does not explicitly resolve the heterogeneities and that no additional degrees of freedom are introduced. In this paper, hybrid-Trefftz stress finite elements are utilised and performance of the proposed formulation is demonstrated with numerical examples. The method is restricted here to elastic particulate composites with ellipsoidal inclusions but it has been designed to be extensible to a wider class of materials comprising arbitrary shaped inclusions.

###### keywords:
Micromechanics; Equivalent inclusion method; Eshelby’s solution; Heterogeneous materials; Hybrid-stress finite elements; Displacement perturbations
journal: Computer Methods in Applied Mechanics and Engineering

## 1 Introduction

In the analysis of materials with complex microstructures, full resolution of the heterogeneities using classical numerical approaches such as the Finite Element method can be computationally prohibitive. To overcome this, one option is to model the macroscale problem using equivalent properties; however, this can lead to a critical loss of information about the finer scale behaviour and poor understanding of the heterogeneities’ influence on the macroscale response. Numerical approaches such as computational homogenization (often called ) provide an alternative strategy Feyel2000309 ; Geers2010 . These techniques comprise nested Finite Element analyses, where each macroscopic material point response is determined via the numerical solution of an RVE subject to the macroscopic strains. Although such approaches have significant potential for certain classes of problems, they are still computationally demanding and are restricted to situations involving clear separation of scales.

The objective of this work is to develop a Finite Element formulation for modelling the macroscopic mechanical problem that is enhanced to capture the influence of the underlying heterogeneities. In our approach, the Finite Element mesh is not required to explicitly resolve the heterogeneities. Closed-form expressions for the perturbation of the mechanical fields due to the presence of the heterogeneities are determined and these are then utilised to enhance the Finite Element formulation.

The ability to capture the effect of microstructural features independently of the underlying finite element mesh has been an ongoing challenge in computational mechanics research. Partition of Unity methods Sukumar:2001:MHI ; Moes:1999:FEM ; Moes:2003:CAHC ; Wells2002 provide a potential solution to this problem, without mesh refinement, by extending a given solution space with additional functions and has been successfully applied to problems such as cracks and material interfaces. The application of this approach in the context of the current work will be briefly discussed in this paper, whereby the closed-form solutions derived for the mechanical perturbation fields are used to extend the classical finite element method. However, it will be shown that there are some disadvantages to this approach for the particular problem at hand and an alternative approach, centred on the Hybrid-Trefftz stress element formulation Kaczmarczyk20091298 , represents the main focus of this paper. This method does not result in additional degrees of freedom, although it does involve an additional, albeit relatively minor, computational overhead.

The heterogeneities, although currently restricted to simple shapes (ellipsoids), can be randomly sized and randomly distributed without reference to the finite element mesh. Therefore, the proposed approach has the potential to be applicable to a wide range of composite materials, such as fibre reinforced composites Kabele2007194 , porous media nicolaou2004hybrid ; Fritsch2009230 , functionally graded materials sharif2008microstructure , etc. Moreover, it can be extended to general inclusion shapes by evaluating the perturbation functions numerically maz2007approximate ; Novak:2008:CESS .

The paper is structured as follows. The methodology of the proposed strategy is described in Section 2. Construction of the perturbation approximation functions for Finite Element Analysis is derived in Section 3. The implementation into the Hybrid-Trefftz stress element formulation containing an arbitrary number of inclusions is presented in Section 4. Section 5 comprises examples demonstrating the model’s performance. Finally we present the conclusions as well as a discussion on future research directions. An appendix is included that highlights some important, but rather technical, aspects of the proposed technique in order to keep the paper self-contained.

## 2 Micromechanics approach

This section outlines the strategy to calculate the perturbation of mechanical fields due to a heterogeneous microstructure, exploiting the Equivalent Inclusion Method mura1987micromechanics in conjunction with analytical micromechanics. Our objective is to convert the heterogeneous problem into an equivalent homogeneous problem and to derive analytical expressions for the perturbations of the stress, strain and displacement fields that we can then utilise within a finite element formulation.

Consider a body consisting of clearly distinguishable heterogeneities in a matrix (Fig. 1a) subjected to a displacement and traction field.

The stiffness of such a material is decomposed as follows mura1987micromechanics ; Novak:2008:CESS

 \boldmathC=\boldmathC0+V\boldmath% C∗ (1)

where is the stiffness tensor of a homogeneous matrix and is due to the presence of inclusions. is nonzero only within the domain , so that

 V={0 in Ω01 in Ω (2)

As a result of the heterogeneities, the mechanical fields (displacement, strains, stresses) experience a perturbation for which we will derive closed-form expressions based on analytical micromechanics. Symbolically, we can express the decomposition of the mechanical fields as follows:

 u=u0+u∗, \boldmathε=\boldmathε0+\boldmathε∗, \boldmathσ=\boldmathσ0+% \boldmathσ∗ (3)

where, the superscript indicates the macroscopic component of the fields in the absence of heterogeneities and superscript indicates the perturbation (or microscopic) component due to the presence of the heterogeneities. It is worth noting that, traditionally, in analytical micromechanics, the macroscopic fields are assumed to be uniform across the domain, e.g. Vorel:2009:SEM ; pichler2008micron . Here it is assumed that they can be position dependent functions of the Neumann and Dirichlet boundary conditions.

The perturbation fields are determined by employing the equivalent inclusion method for a single heterogeneity embedded in a matrix and then extended here for multiple heterogeneities. In the equivalent inclusion method, the heterogeneous solid is replaced by an equivalent homogeneous solid with uniform material stiffness everywhere (Fig. 1a, b) and suitable stress-free eigenstrains applied in the inclusions so that the homogeneous equivalent solid has the same mechanical fields as the original heterogeneous solid.

### 2.1 Equivalent inclusion method for single heterogeneity problem

Consider first a single heterogeneity embedded in a matrix. Following Eshelby’s fundamental work eshelby1957determination , this problem can be decomposed into two problems of known solution and then assembled back via superposition mura1987micromechanics ; eshelby1957determination , see Fig. 2.

In brief, the solution of the inhomogeneity problem requires the determination of the transformation eigenstrain that induces the identical local mechanical response as the original heterogeneous body. Due to the absence of other inclusions, no strain or stress concentrations are induced and thus remains constant eshelby1957determination .

In the original heterogeneous body, the stress state can be expressed as

 \boldmathσ=\boldmathσ0+\boldmathσ∗=\boldmathC:[\boldmathε0+\boldmathε∗] (4)

In the homogeneous solid, we add a stress-free eigenstrain inside the domain of the inclusion, which has the same material stiffness as the matrix such that

 \boldmathσ=\boldmathC0:[% \boldmathε0+\boldmathε∗−\boldmathετ] (5)

It should be noted that in the matrix.

Given that the macroscopic stress is , it can be see from equating Eqs. (4) and (5), that the stress perturbation in the homogeneous solid can be expressed as

 \boldmathσ∗=\boldmathC0:[% \boldmathε∗−\boldmathετ] (6)

Furthermore, equating Eqs. (4) and (5) also results in the following expression:

 \boldmathC:[\boldmathε0+% \boldmathε∗]=\boldmathC0:[% \boldmathε0+\boldmathε∗−\boldmathετ] (7)

where the transformation eigenstrain is as yet unknown. Eshelby’s solution of the homogeneous inclusion problem eshelby1957determination , relates the eigenstrain to the perturbation strain as follows

 \boldmathε∗=\boldmathS:% \boldmathετ (8)

where denotes the Eshelby tensor and is a function of the heterogeneity’s geometry and the material stiffness of both the matrix and heterogeneity. Substituting this expression into Eq. (7) and rearranging yields

 [\boldmathC−\boldmathC0]:\boldmathε0=[\boldmathC0:\boldmathS−\boldmathC:%\boldmath$S$−\boldmathC0]:% \boldmathετ (9)

This can be recast in a compact form to give an expression for the eigenstrain which depends on the homogeneous strain , the material stiffness of both the matrix and heterogeneity and the Eshelby tensor as

 \boldmathετ=\boldmathB:%\boldmath$ε$0 (10)

where the tensor is provided by:

 \boldmathB=−[\boldmathC∗:\boldmathS+\boldmathC0]−1:\boldmathC∗ (11)

Once the transformation eigenstrain has been determined, the stress perturbation can be computed from Eq. (6) in the form

 \boldmathσ∗=\boldmathC0:[% \boldmathS−\boldmathI]:\boldmathB:\boldmathε0 (12)

It can be seen that the stress perturbation depends on stiffness of the different material phases, the macroscopic strain field and the geometry of the heterogeneity. This closed-form expression for the stress perturbation is at the heart of the proposed finite element enrichment to be discussed later in this paper. It is also useful to derive an expression for the displacement perturbation field as follows

 u∗=\boldmathL:\boldmathετ=\boldmathL:\boldmathB:\boldmathε0 (13)

where the operator is a third order tensor, mapping . For the sake of conciseness and to keep the paper self-contained, the detailed derivation of this operator can be found in Appendix A.

### 2.2 Self-compatibility algorithm for multiple inclusions

In the case of multiple inclusions, the mechanical perturbation fields within individual inclusion domains are no longer uniformly distributed as a result of their mutual interaction. Here, we account for this approximately by assuming that the eigenfields to be piecewice constant within each inclusion. Thus, the perturbations are determined from the Eshelby solution for each individual inclusion, as described above, plus an iterative self-compatibility procedure (Tab. 1) to ensure that the solution correctly reflects the influence of multiple heterogeneities.

This procedure iteratively enforces compatibility (see eg. pichler2010estimation for further reference) between the imposed macroscopic strain and the average eigenstrain inside any given inclusion , so as to account for the influence of the remaining inclusions . An iterative algorithm has been chosen because a closed form solution for multiple inclusions does not exist and a numerical solution would be prohibitively expensive pichler2010estimation .

First, the eigenstrain for each inclusion is calculated (Eq. 10) without reference to the other inclusions (Line 2). Next, the associated perturbation strain for each inclusion is evaluated (Eq. 8) at the centre of all other inclusions (Line 3). The mutual interaction of inclusions is then taken into account via a correction of the eigenstrain (). For each inclusion , this correction is calculated (Line 8) from the inverse of the inclusion’s Eshelby tensor and the perturbation strains of all other inclusions, evaluated at the centroid of inclusion . The perturbation strain resulting from inclusion at the centre of inclusion is denoted as . This is demonstrated in Fig. 3 for a two inclusion problem in 1D. The eigenstrain correction is then used to calculate the correction to the perturbation strains (Line 10). The algorithm continues until a small Euclidean norm between the last two iterations of the total eigenstrains is achieved. At convergence, the corresponding stress and displacement perturbations are recalculated from the corrected transformation eigenstrains.

It is worthwhile noting that the algorithm does not depend on a particular sequence of inclusions, as follows from the elastic reciprocity theorem (pichler2010estimation, , and references therein). Moreover, since the perturbation fields are calculated for the entire macroscopic domain (no RVE is considered), stress admissibility and strain compatibility, in the sense of macro-micro field relations, are fulfilled a priori.

The computational complexity of the Self-compatibility algorithm is . However, this can be improved by taking into account only those inclusions which have a non-negligible influence to the inclusion of interest . Preliminary studies have shown this to give a significant computational speed-up and will be reported in a future paper.

## 3 Construction of perturbation approximation functions for FEA

The above methodology can be utilised to formulate an enhanced Finite Element formulation. The primary task is to determine appropriate approximation functions for the mechanical perturbation fields , and based on the analytical micromechanics developed above and which can then augment the standard macroscopic field approximations. It should be noted that the Voigt-Mandel notation is exclusively used in the forthcoming text.

The perturbation field approximation functions are determined a priori as a linear combination of the perturbation fields evaluated analytically for six load cases, with self-equilibrium enforced by means of the self-compatibility algorithm outlined above (Tab. 1). Each load case corresponds to a unit component of the macroscopic strain vector and the resulting analytically determined stress, strain and displacement perturbation fields are arranged, column-by-column into , and matrices, respectively:

 s∗=[1σ∗…6σ∗],e∗=[1ϵ∗…6ϵ∗] (14)
 u∗=[1u∗…6u∗] (15)

where the left superscript refers to a specific load case, to .

### 3.1 Partition of unity method

Partition of Unity (PU) Methods (for example the eXtended Finite Element Method) extend the underlying basis functions used for interpolating the displacement field by adding an appropriate set of additional functions. Following Melenk:1996:PUM ; Babuska:1997:PUM it has been shown that the displacement field within an element can be interpolated by

 u(x)=n∑i=1(Ni(x)ai+NiNγ(x)bi) (16)

where is the number of nodes per element, is the standard matrix of element shape functions for node , is the identity matrix and the standard displacement degrees of freedom at node . is a matrix containing the additional basis terms and are the associated additional degrees of freedom at node . It is important to recognise that the element shape functions form a partition of unity, i.e.

 n∑i=1Ni(x)=1 (17)

The six analytically derived displacement perturbation functions contained in can be used as the additional functions to augment the standard basis functions. Thus

 Nγ(x)=⎡⎢ ⎢⎣1u∗1(x)…6u∗1(x)0…00…00…01u∗2(x)…6u∗2(x)0…00…00…01u∗3(x)…6u∗3(x)⎤⎥ ⎥⎦ (18)

With this at hand, the PU-based finite element formulation can be derived, see for example Wells2001 . PU methods are particularly useful in problems where the extension of the basis functions is introduced on a node by node basis, so that additional degrees of freedom are only introduced at nodes where the basis is extended. One obvious example of such a local feature that can be modelled in this way is discrete cracks Wells2001 . However, this favourable property is not exploited here because we wish to model a large number of heterogeneities throughout the domain. In 3D problems, there are 3 standard displacement degrees of freedom per node; this would be extended by an additional 18 degrees of freedom per node and per heterogeneity with the proposed approach.

It is also worth noting that for standard finite elements, the volume integration of the discrete system of equations is relatively straightforward. However, extension of the basis functions to include the perturbation functions in Eq. (18) makes this process significantly more arduous. For these reasons, an alternative Finite Element approach using Hybrid Trefftz Stress elements Freitas1998 ; Kaczmarczyk20091298 is considered where the standard basis function is not extended, as with PU methods, but enhanced such that no additional degrees of freedom are introduced. This is described in the next section.

## 4 Hybrid-Trefftz stress element formulation

In this section a finite element formulation based on an enhancement of a hybrid-Trefftz stress (HTS) element formulation Kaczmarczyk20091298 is presented.

The problem requires a solution to the displacement and stress fields as a result of given boundary displacements and tractions on and , respectively. The displacement and stress fields must fulfil the following governing equations:

 LTσ=0in Ωe\ldots Cauchy~{}equilibrium~{}equationLu=ϵin Ωe% \ldots strain-displacement relationshipσ=Cϵin Ωe% \ldots constitutive~{}equationNσ=¯ton Γeσ\ldots static~{}boundary~{}conditionsu=¯¯¯uon Γeu\ldots kinematic% ~{}boundary~{}conditions (19)

where and are the column matrix representation of the second order stress and strain tensor, respectively, represents the displacement vector, is the matrix representation of fourth order stiffness tensor and finally and represent the applied displacements and tractions, respectively. The gradient operator and the matrix of directional cosines of the outward normal to element boundary have the following forms bittnar1996numerical

 L=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∂/∂x000∂/∂y000∂/∂z∂/∂y∂/∂x00∂/∂z∂/∂y∂/∂z0∂/∂x⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,N=⎡⎢⎣nx000nzny0ny0nz0nx00nznynx0⎤⎥⎦ (20)

### 4.1 Stress, strain and displacement approximations

The macroscopic stress field within the HTS element is approximated as

 σ0=S0vv in Ωe (21)

where is the vector of generalised stress degrees of freedom, denotes the matrix of stress approximation functions chosen so as to automatically satisfy the equilibrium conditions Eq. (19). Thus,

 LTS0vv=0 in Ωe (22)

and

 t=NS0vv on Γe (23)

where represents the traction vector induced by the macroscopic stress approximation field.

The macroscopic strain and displacement fields are expressed analogously to Eq. (21) as

 ϵ0=E0vvandu0=U0vv (24)

where and are directly associated with the stress approximation by means of the compatibility equation (19) and constitutive equation (19) as

 S0v=C0E0v=C0LU0v (25)

Since the stress approximation functions are typically polynomial functions, the integration of to get is relatively straightforward.

Rather than extend the solution space to capture the influence of the heterogeneities, as was briefly described in Section 3.1, here we enhance the macroscopic approximations to include the influence of the heterogeneities, thereby not increasing the number of unknowns. The total stress (macroscopic plus perturbation) field within the HTS element is approximated, following Eq. (3), as

 σ=σ0+σ∗=(S0v+S∗v)v in Ωe (26)

where is the perturbation counterpart to . can be constructed from Eq. (14):

 σ∗=s∗ϵ0 (27)

where is the set of analytically defined stress perturbations for six load cases, each one representing a unit component of the macroscopic strain Eq. (14). For the purposes of constructing , we approximate the macroscopic strain field as constant within each finite element and computed as the volume average of the actual macroscopic strain field. From Eq. (24),

 ϵ0,ave=1|Ωe|∫Ωeϵ0dΩe=1|Ωe|∫ΩeEvdΩev=Eavevv (28)

Substituting for into Eq. (27) leads to

 σ∗=s∗Eavevv (29)

Thus, the matrix of stress perturbation approximation functions is:

 S∗v=s∗Eavev (30)

Analogously, the approximation of total strain and displacement fields within the element domain is given by

 ϵ=(E0v+E∗v)vandu=(U0v+U∗v)v (31)

where the perturbation approximation matrices and are, as with their macroscopic counterparts, directly associated with the stress approximation as

 S∗v=C0E∗v=C0LU∗v (32)

It is worthwhile noting that the stress perturbation fields and the corresponding traction perturbation fields, approximated as and , remain in self-equilibrium.

### 4.2 Static boundary conditions

Contrary to general condition in Eq. (19), the equilibrium on the element traction boundary is imposed only in the weighted residual sense as:

 ∫ΓeσWT1[N(σ0+σ∗)−¯t]dΓe=0 (33)

along with representing the matrix of weighting functions. Replacing the total stress field in Eq. (33) by its approximation from Eq. (26), the traction boundary condition becomes

 ∫ΓeσWT1N(S0v+S∗v)vdΓe=∫ΓeσWT1¯tdΓe (34)

### 4.3 Kinematic boundary conditions

Compatibility inside the element domain is also enforced in a weighted residual sense, such that:

 ∫ΩeWT2(ϵ0+ϵ∗−Lu)dΩe=0 (35)

Next, utilising integration by parts and applying Green’s theorem to , Eq. (35) results in

 ∫ΩeWT2(E0v+E∗v)vdΩe + ∫Ωe(LTW2)TudΩe (36) − ∫Γeσ(NW2)TudΓe=∫Γeu(NW2)T¯¯¯udΓe

With the current formulation, it is not necessarily possible to find a solution to both Eqs. (34) & (36) that satisfies both traction and kinematic boundary conditions acting on the element boundary. As a consequence, Eq. (36) is relaxed by introducing an additional and independent approximation of displacements on the element traction boundary:

 uΓ=UΓq in Γeσ (37)

Here, stands for the set of displacement unknowns and is the matrix of boundary displacement approximation functions. Such a formulation of the stress element leads to a hybrid approach  Kaczmarczyk20091298 ; Freitas1998 .

Given the above consideration, introducing Eq. (37) into Eq. (36) results in

 ∫ΩeWT2(E0v+E∗v)vdΩe + ∫Ωe(LTW2)TudΩe (38) − ∫Γeσ(NW2)TuΓdΓe=∫Γeu(NW2)T¯¯¯udΓe

### 4.4 Weighting functions

In order to achieve an energy-consistent formulation, it is required that all weighted terms within the integrals defined above have the dimension of work. The weighting functions then directly follow from the integrands in Eq. (33) and Eq. (35) representing the increment of internal work of strains within and external work of tractions on , respectively. These functions thus admit the following forms:

 W1=UΓ, W2=S0v (39)

First, introducing Eq. (39) into Eq. (38) and taking into account condition (22) yields

 ∫Ωe(S0v)T(E0v+E∗v)vdΩe−∫Γeσ(NS0v)TUΓqdΓe=∫Γeu(NS0v)T¯¯¯udΓe (40)

Second, introducing Eq. (39) into the traction boundary condition (34) yields

 ∫ΓeσUTΓN(S0v+S∗v)vdΓe=∫ΓeσUTΓ¯tdΓe (41)

Combining Eqs. (40) & (41) results in a coupled system of linear equations that can be expressed in compact form as

 [F−AT−(A+A∗)0][cvq]=[pu−pσ] (42)

where the submatrices on the left-hand side follow from Eqs (4041) and are given by the following integrals

 F = ∫Ωe(S0v)T(E0v+E∗v)dΩe=∫ΓeN(S0v)T(U0v+U∗v)dΓe (43) A = ∫ΓeσUΓNS0vdΓeandA∗=∫ΓeσUΓNS∗vdΓe (44)

and for the terms on the right-hand side it holds

 pu=∫Γeu(NS0v)T¯¯¯udΓe,~{}andpσ=∫ΓeσUTΓ¯tdΓe (45)

Note that Eq. (43) illustrates that the  matrix can be evaluated via a boundary rather than volume integral. Thus all terms in Eq. (42) can be evaluated using boundary integrals only.

The size of the system of equations to be solved simultaneously can be reduced via static condensation, representing a significant reduction in computational effort. First, from the first equation in Eq. (42), the generalised stress degrees of freedom are expressed in terms of the displacement degrees of freedom as

 v=F−1(pu+ATq) (46)

This is then substituted into the second equation of Eq. (42) to yield:

 [(A+A∗)F−1AT]q=pσ−(A+A∗)F−1pu (47)

This sparse system of equations is then solved for the displacement degrees of freedom . Subsequently, the stress degrees of freedom can then be calculated on an element-by-element basis.

Our implementation of these HTS elements for composite materials (C-HTS elements) utilises displacement degrees of freedom that are associated with element faces rather than vertices. This has the advantage that the bandwidth of the stiffness matrix is minimised, as is interprocessor communication.

## 5 Numerical Examples

The performance of the key components of the proposed strategy (micromechanical solution, self-compatibility algorithm, finite element analysis convergence, etc.), in terms of efficiency and accuracy, have been explored through two numerical examples.

### 5.1 Three ellipsoidal inclusions in matrix

This example comprises three ellipsoidal inclusions embedded in a cube of matrix. The geometry of the problem is illustrated in Fig. 5 and details of the ellipsoids are given in Table 2, including the semi-axes’ dimensions, centroid coordinates and Euler angles , and , which are successive rotations of the semi-axes , and about global coordinate axes , and , respectively. The cube has side lengths of . The displacement boundary conditions were prescribed on faces , and as , and , respectively. The remaining faces at , and were subject to uniform normal unit tractions. The Young’s modulus for the homogeneous matrix was chosen as and for the heterogeneities as . Poisson’s ratio was chosen as for both materials. All units are consistent.

It is worthwhile noting that the small contrast in stiffness between the two materials was chosen deliberately to maximise the extent of the perturbation fields emanating from the heterogeneities. Large contrasts in stiffness between the matrix and heterogeneities lead to perturbation fields that decay rapidly with distance from the heterogeneities. The close proximity of the three ellipsoidal heterogeneities to each other was also chosen deliberately in order to demonstrate the ability of the formulation to capture the interaction of multiple heterogeneities. Furthermore, the close proximity of one of the ellipsoids to a traction boundary was chosen to demonstrate the ability of the formulation to capture boundary effects.

The problem was analysed using two three-dimensional finite element meshes with different densities, comprising C-HTS elements. The coarse mesh comprised elements and  DOFs (Fig. 6a) and the second, refined, mesh comprised elements and  DOFs (Fig. 6b).

Results from the two enhanced finite element analyses are plotted in the -plane (at ). Fig. 7a and Fig. 8a show the two meshes in this plane. The stress component for both analyses are shown in Fig. 7b and Fig. 8b.

In addition, a reference finite element analysis of the same problem was undertaken for comparison sake. The reference analysis utilised HTS elements but without the proposed enhancement. Unlike the other two analyses, the reference analysis utilised a mesh that explicitly resolved the three ellipsoidal heterogeneities and comprised tetrahedrons with  DOFs (Fig. 9a). The corresponding mesh and stress results of the reference analysis are shown in Fig. 9. Comparison of the stress results from the enhanced formulation and the reference analysis leads to the relative error plots shown in Fig. 7c and Fig. 8c. It can be seen that even the very coarse mesh with the enhanced formulation results in good agreement.

Further comparison of the stress results is shown in Fig. 10, where it can been seen that along the traction boundary, the finer mesh of enhanced elements is able to capture the imposed constant stress field more accurately than the very much finer reference mesh.

The perturbation fields are based on the assumption of a heterogeneity in an infinite medium but the enhanced formulation still exhibits convergence in the regions strongly influenced by the traction boundary.

### 5.2 L-shaped specimen

The proposed modelling strategy is also demonstrated on an example with a large number of inclusion. A 3D L-shaped specimen with fully fixed boundary conditions on the right surface of the right-hand arm and normal traction applied on the top surface of the left-hand arm is analysed, see Fig. 11. The length of the plate is in both and direction, the depth is in direction. The Young moduli were chosen as and for matrix and inclusion respectively. Poisson’s ratio was for both material phases. The microstructure comprised 2,523 spherical inclusions varying in size between and with a uniform spatial distribution (Fig. 11b). All units are consistent.

The solution for three different mesh densities (Fig. 11c, d & e) are compared in Fig. 12. These results are plotted on the - mid-plane. Fig. 12(a–d) shows a plot of the stress component for the three different meshes. These results show that the complex stress distribution resulting from the heterogeneities can be captured and that the solution is converging with mesh refinement. Further local mesh refinement near corners and stress concentrations is possible, although this was not undertaken in this case. As an estimate of the computational overhead of the enhanced formulation for this particular problem on 16 processors, we note that the total solution for C-HTS elements was s, whereas the problem on the identical mesh of HTS elements for the equivalent homogeneous problem consumed s of computer time. This represents a large increase of computational time in comparison to the homogenous problem. However, it should be noted that not all aspects of the solution procedure have been parallelised. It should also be noted that it was not possible to obtain the reference solution for this problem by means of conventional FEA due to the excessive number of degrees of freedom associated with a mesh required to resolve all of the heterogeneities. The mesh generation itself was not possible using our currently available software and hardware facilities. Furthermore, comparison with a PUM based solution, was also not undertaken due to the complexity of the problem.

## 6 Conclusions

A new micromechanics-enhanced finite element formulation has been presented for modelling the influence of a large number of heterogeneities in composite materials in a computationally efficient manner. The strategy exploits closed form solutions derived by Eshelby for ellipsoidal inclusions in order to determine the mechanical perturbation fields as a result of the underlying heterogeneities. Approximation functions for these perturbation fields are then incorporated into a finite element formulation to augment those of the macroscopic fields. A significant feature of this approach is that the finite element mesh does not explicitly resolve the heterogeneities, although the resulting solution still explicitly accounts for their presence. In contrast with traditional homogenization approaches, this method does not rely on separation of scale and does not suffer from loss of information due to averaging or localization.

The proposed technique has been implemented into a hybrid-Trefftz stress (HTS) element formulation and it has been shown that the resulting enhanced elements (C-HTS) require significantly fewer degrees of freedom to capture the detailed mechanical response compared to standard finite elements. The paper also outlines how the proposed micromechanics approach could be used within a Partition-of-Unity (PoU) formulation, although we conclude that this does not fully exploit the advantages of PoU methods and that the proposed hybrid-Trefftz formulation is most appropriate.

A self-compatibility algorithm is used to determine the mutual interactions between inclusions, assuming that the eigenstrain fields are uniform within the domain of each inclusion. It was found that even for topologies exhibiting extremely small distances between the inclusions, this assumption is sufficient. Furthermore, it has been shown that boundary effects, that are not accounted for by the classical micromechanical solution due to the assumption of an infinite medium, can be captured through local mesh refinement.

We have implemented this formulation into our FE code that is optimized for parallel computing. Additional parallelization of the micromechanical aspects of the formulation needs to be investigated for increased efficiency. Further research is required in order to incorporate other improvements such as nonuniform eigenstrains mura1987micromechanics , debonding effects ju2008micromechanical and inclusions of arbitrary shape by evaluating the perturbation functions numerically maz2007approximate ; Novak:2008:CESS .

#### Acknowledgements

Funding by the Glasgow Research Partnership in Engineering (GRPE) under project “Multi-scale modelling of fibre reinforced composites” is gratefully acknowledged. This research was partially supported by the Czech Science Foundation through project GAČR 103/09/P490 and by the Ministry of Education, Youth and Sports of the Czech Republic trough project MSM 684077003. All the analyses were performed by means of YAFFEMS FE code. For more details we refer to the code’s homepage at http://code.google.com/p/yaffems.

## Appendix A Detailed solution of perturbation displacements

The displacement perturbation field in an infinite homogeneous material due to a uniform eigenstrain applied to an ellipsoidal region is provided by the following integral equation (mura1987micromechanics, , Eq. (11.30))

 u∗i=18π(1−ν)[Ψjk,jki−2νΦkk,i−4(1−ν)Φik,k] (48)

where denotes the Poisson’s ratio and the elliptic potentials and are defined as (mura1987micromechanics, , Eq. (11.32))

 Ψij=ετij∫Ω|x−x′|dx′=ετijψ,andΦij=ετij∫Ω1|x−x′|dx′=ετijϕ (49)

The integrals and in Eq. (49) are the harmonic and bi-harmonic potentials respectively. Note that in Eq. (48) and thereafter, standard index notation is employed, together with the generalised summation convection due to Mura mura1987micromechanics . Thus, a repeated index is summed according to the Einstein summation rule (e.g. ), whereas a non-repeated upper-case index equals to the lower-case equivalent (e.g. ). The symbol denotes the partial derivative of with respect to the coordinate .

By combining Eq. (48) and Eq. (49), we obtain

 u∗i=18π(1−ν)[ετjkψ,jki−2νδjkετjkϕ,i−4(1−ν)δijετjkϕ,k] (50)

Similarly to Eshelby’s approach eshelby1957determination , the displacement perturbations are expressed in compact form:

 u∗i=Lεijkετjk, Lεijk=18π(1−ν)[ψ,jki−2νδjkϕ,i−4(1−ν)δijϕ,k] (51)

where the third-order operator maps a transformation eigenstrain to the displacement perturbation field . It is therefore analogous to the well-known Eshelby tensor eshelby1957determination , which relates a transformation eigenstrain to the strain perturbation field.

The operator can be conveniently expressed in terms of the Ferrers-Dyson elliptic integrals, e.g. (mura1987micromechanics, , Eq. (11.36))

 I(λ) = 2πa1a2a3∫∞λdsΔ(s), Ii(λ) = 2πa1a2a3∫∞λds(a2i+s)Δ(s), Iij(λ) = 2πa1a2a3∫∞λds(a2i+s)(a2j+s)Δ(s) (52)

where stands for the -th semi-axis of ellipsoid and is obtained from

 Δ2(s)=3∏i=1(ai+s)2 (53)

The variable is the largest positive root of equation (mura1987micromechanics, , Eq. (11.37))

 xixi(aI+λ)2=1 (54)

Notice that is generally position dependent and non-zero for the points placed outside the inclusion domain , hence called the exterior points. Contrary, for interior points.

All integrals in (52) admit a closed-form expression in terms of the Legendre-Jacobi integrals of the first and second kind, defined as a fuction of an auxiliary angle (mura1987micromechanics, , Eq. (12.17)). It is worth noting that its definition via (mura1987micromechanics, , Eq (11.18))

 θ=sin−1 ⎷1−a23a21 (55)

is valid for interior points only and not everywhere as stated in mura1987micromechanics . Thus, it needs to be replaced with a general formula:

 θ=sin−1 ⎷a21−a23a21+λ (56)

available, e.g. in eshelby1959elastic ; rahman2001newtonian . Moreover, the following identity (mura1987micromechanics, , Eq. (11.40.4))

 [xnxnIi…jN(λ)],p=2xpIi…jP+Ii…j,p(λ) (57)

will be repeatedly proved useful in the sequel.

It follows from Eq. (51) that to express operator , we need to evaluate the first-order derivatives of the potential and the third-order derivatives of . To this end, we start with expressions (mura1987micromechanics, , Eq. (11.38))

 ϕ(λ)=12[I(λ)−xnxnIN(λ)] (58)

and

 ψ,i(λ)=12xi{2ϕ(λ)−a2I[II(λ)−xnxnIIN(λ)]}=12xiQ(λ) (59)

Employing Eq. (57), the first derivative of becomes

 ϕ,i(λ)=12{I,i(λ)−[xnxnIN(λ)],i}=12[I,i(λ)−2xiII(λ)−I,i(λ)]=−xiII(λ) (60)

The third derivative of potential is expressed from Eq. (59) as

 ψ,ijk(λ)=12[δijQ,k(λ)+δikQ,j(λ)+xiQ,jk(λ)] (61)

With the help of Eqs. (60) and (57), the term can be evaluated from

 Q,j(λ)=2ϕ,j(λ)−a2I[II(λ)−xnxnIIN(λ)],j=2xj[a2IIIJ(λ)−IJ(λ)] (62)

This provides the second derivatives of in the form

 Q,jk(λ)=2{δjk[a2IIIJ(λ)−IJ(λ)]+xj[a2IIIJ,k(λ)−IJ,k(λ)]} (63)

After utilising the derivatives of and re-ordering the indices, Eq. (61) becomes

 ψjki(λ) =xiδjk[a2JIJI(λ)−II(λ)]+xkδji[a2JIJK(λ)−IK(λ)] +xjδki[a2JIJK(λ)−IK(λ)]+xjxk[a2JIJK,i(λ)−IK,i(λ)] (64)

with and provided by (mura1987micromechanics, , Eqs. (11.40.1, 11.40))

 Ii…jk,p(λ) = −2πa1a2a3(a2i+λ)…(a2j+λ)(a2k+λ)Δ(λ)λ,p, λ,p = 2xpa2P+λ(a2I+λ)2xixi (65)

Now we are in a position to evaluate in terms of the Ferrers-Dyson integrals. Introducing Eqs. (60) and (A) into Eq. (51) gives

 [8π(1−ν)]Lεijk =xiδjk[a2JIJI(λ)−II(λ)]+(xkδji+xjδki)[a2JIJK(λ)−IK(λ)] +xjxk[a2JIJK,i(λ)−IK,i(λ)]+2νδjkxiII(λ)+4(1−ν)δijxkIK(λ)} (66)

Since for the points inside the inclusion, recall Eq. (54), all derivatives of vanish as well. Therefore, Eq. (A) yields

 [8π(1−ν)]Lε,int =xiδjk[a2JIJI(λ)−II(λ)]+(xkδji+xjδki)[a2JIJK(λ)−IK(λ)] +νδjkxiII(λ)+4(1−ν)δijxkIK(λ) (67)

and Eq. (51) receives its final form

 Lεijk=Lε,intijk+18π(1−ν)xjxk[a2JIJK,i(λ)−IK,i(λ)] (68)

For implementation purposes, it is worth noting that Eq. (51) admits the Voigt representation:

 (69)

## References

• (1) F. Feyel, J.-L. Chaboche, FE multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials, Computer Methods in Applied Mechanics and Engineering 183 (3-4) (2000) 309 – 330.
• (2) M. Geers, V. Kouznetsova, W. Brekelmans, Multi-scale computational homogenization: Trends and challenges, Journal of Computational and Applied Mathematics 234 (2010) 2175 – 2182.
• (3) N. Sukumar, D. Chopp, N. Moës, T. Belytschko, Modeling holes and inclusions by level sets in the extended finite-element method, Computer Methods in Applied Mechanics and Engineering 190 (46–47) (2001) 6183–6200.
• (4) N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, International Journal for Numerical Methods in Engineering 46 (1) (1999) 131–150.
• (5) N. Moës, M. Cloirec, P. Cartraud, J.-F. Remacle, A computational approach to handle complex microstructure geometries, Computer Methods in Applied Mechanics and Engineering 192 (28–30) (2003) 3163–3177.
• (6) G. Wells, L. Sluys, R. D. Borst, Simulating the propagation of displacement discontinuities in a regularised strain-softening medium, International Journal for Numerical Methods in Engineering 53 (5) (2002) 1235–1256.
• (7) L. Kaczmarczyk, C. J. Pearce, A corotational hybrid-Trefftz stress formulation for modelling cohesive cracks, Computer Methods in Applied Mechanics and Engineering 198 (15-16) (2009) 1298 – 1310.
• (8) P. Kabele, Multiscale framework for modeling of fracture in high performance fiber reinforced cementitious composites, Engineering Fracture Mechanics 74 (1-2) (2007) 194 – 209, fracture of Concrete Materials and Structures.
• (9) P. Nicolaou, S. Semiatin, A hybrid micromechanical-macroscopic model for the analysis of the tensile behavior of cavitating materials, Metallurgical and Materials Transactions A 35 (13) (2004) 1141–1149.
• (10) A. Fritsch, C. Hellmich, L. Dormieux, Ductile sliding between mineral crystals followed by rupture of collagen crosslinks: Experimentally supported micromechanical explanation of bone strength, Journal of Theoretical Biology 260 (2) (2009) 230 – 252.
• (11) Z. Sharif-Khodaei, J. Zeman, Microstructure-based modeling of elastic functionally graded materials: One dimensional case, Journal of Mechanics of Materials and Structures 3 (2008) 1773–1796. arXiv:0802.0511.
• (12) V. Maz’ya, G. Schmidt, Approximate approximations, American Mathematical Society, Providence, RI, 2007.
• (13) J. Novák, Calculation of elastic stresses and strains inside a medium with multiple isolated inclusions, in: M. Papadrakakis, B. Topping (Eds.), Proceedings of the Sixth International Conference on Engineering Computational Technology, Stirlingshire, UK, 2008, p. 16 pp, paper 127.
• (14) T. Mura, Micromechanics of Defects in Solids., Martinus Nijhoff Publishers, P. O. Box 163, 3300 AD Dordrecht, The Netherlands, 1987. 587.
• (15) J. Vorel, M. Šejnoha, Evaluation of homogenized thermal conductivities of imperfect carbon-carbon textile composites using the mori-tanaka method, Structural Engineering and Mechanics 33 (4) (2009) 429–446.
• (16) B. Pichler, S. Scheiner, C. Hellmich, From micron-sized needle-shaped hydrates to meter-sized shotcrete tunnel shells: micromechanical upscaling of stiffness and strength of hydrating shotcrete, Acta Geotechnica 3 (4) (2008) 273–294.
• (17) J. D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 241 (1226) (1957) 376–396.
• (18) B. Pichler, C. Hellmich, et al., Estimation of influence tensors for eigenstressed multiphase elastic media with nonaligned inclusion phases of arbitrary ellipsoidal shape, Journal of Engineering Mechanics 136 (2010) 1043–1053.
• (19) J. Melenk, I. Babuška, The partition of unity finite element method: Basic theory and applications, Computer Methods in Applied Mechanics and Engineering 139 (1–4) (1996) 289–314.
• (20) I. Babuška, J. Melenk, The partition of unity method, International Journal for Numerical Methods in Engineering 40 (4) (1997) 727–758.
• (21) G. Wells, Discontinuous modelling of strain localisation and failure, Ph.D. thesis, Delft University of Technology (2001).
• (22) J. Teixeira de Freitas, Formulation of elastostatic hybrid-trefftz stress elements, Computer Methods in Applied Mechanics and Engineering 153 (1998) 127 – 151.
• (23) Z. Bittnar, J. Šejnoha, Numerical methods in structural mechanics, American Society of Civil Engineers, 1996.
• (24) J. Ju, Y. Ko, Micromechanical Elastoplastic Damage Modeling of Progressive Interfacial Arc Debonding for Fiber Reinforced Composites, International Journal of Damage Mechanics 17 (4) (2008) 307.
• (25) J. D. Eshelby, The elastic field outside an ellipsoidal inclusion, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences (1959) 561–569.
• (26) M. Rahman, On the Newtonian potentials of heterogeneous ellipsoids and elliptical discs, Proceedings: Mathematical, Physical and Engineering Sciences 457 (2013) (2001) 2227–2250.