Mean-Field Homogenization

Mean-field homogenization:

is a semi-analytical homogenization approach that can be used to compute thermal and mechanical properties of a composite material;
can also predict volume-averaged microlevel solutions in each constituent of the composite; and
uses formulations based on Eshelby's solution.

This page discusses:

Specifying Constituent Properties for a Composite Material
Mean-Field Homogenization for Linear Elastic Composites
Incremental Mean-Field Homogenization for Nonlinear Composites
Multistep Homogenization
Specifying the Microstructure of the Composite
Composites with Thermal Expansion
Composites with Viscoelasticity
Composites with Damage
Example
Elements
Output
References

Specifying Constituent Properties for a Composite Material

In multiscale material modeling, each composite material consists of one matrix material and one or more inclusions including voids. In Abaqus/Explicit the composite material must consist of one matrix material and only one inclusion; the inclusion material must be linear elastic with constant properties, and the matrix material can be only linear elastic or elastic-plastic with isotropic hardening.

Mean-Field Homogenization for Linear Elastic Composites

The mean-field strain and stress in each phase I is defined as

{\bar{σ}}_{I} = {〈 σ 〉}_{I}

{\bar{ϵ}}_{I} = {〈 ϵ 〉}_{I},

where

〈 f 〉 = \frac{1}{V} \int_{V}^{} f d V .

The averaged macro fields can be written as

〈 • 〉 = v_{f}^{M} {〈 • 〉}_{M} + \sum_{I = 1}^{N} v_{f}^{I} {〈 • 〉}_{I},

where $v_{f}^{M}$ is the volume fraction of the matrix phase and $v_{f}^{I}$ is the volume fraction of the inclusion phase. For the single-inclusion case, if all the constituents are linear elastic, the strain in the inclusion is related to the strain in the matrix through a concentration tensor $A$ .

{〈 ϵ 〉}_{I} = A {〈 ϵ 〉}_{M} .

The form of the concentration tensor defines the homogenization methods below.

Voigt and Reuss Models

The simplest mean-field homogenization formulations are the Voigt and Reuss models. These models do not take into account the shape or the orientation of the inclusions; however, they provide upper and lower bounds of the macro stiffness modulus and, therefore, can be used for validation. The Voigt model is also used in a two-step approach to model composites with multiple inclusions.

The Voigt model assumes uniform strain in the composite ${〈 ϵ 〉}_{I} = {〈 ϵ 〉}_{M}$ , which gives

A = I .

The Reuss model assumes uniform stress in the composite ${〈 σ 〉}_{I} = {〈 σ 〉}_{M}$ , which gives

A = C_{I}^{- 1} C_{M},

where $C_{M}$ is the stiffness of the matrix and $C_{I}$ is the stiffness of the inclusion.

Eshelby-Based Models

The more sophisticated mean-field homogenizations are based on Eshelby's solution. In his 1957 paper (Eshelby, 1957), Eshelby solved a single-inclusion problem described in Figure 1.

Eshelby found that the strain inside the constrained ellipsoid is uniform and given by

ϵ (x) = E : ϵ^{*},

where $x$ represents a point inside the inclusion and $ϵ^{*}$ is the eigenstrain. The Eshelby tensor (E) is a function of the stiffness of the matrix and also the shape and orientation of the inclusion. To calculate Eshelby's tensor usually involves numerical integration over the surface of the ellipsoidal inclusion; however, analytical solutions are available (Mura, 1987) for simple shapes such as a sphere, cylinder, penny, oblate, and prolate. Based on Eshelby's solution, the following models are proposed based on different assumptions about the properties of the outer material surrounding the inclusions.

Mori-Tanaka Model

The Mori-Tanaka formulation assumes that each inclusion behaves like an isolated inclusion and the strain in the matrix is considered as the far-field strain. Therefore, the concentration tensor can be written as

A = {[E : (C_{M}^{- 1} C_{I} - I) + I]}^{- 1} .

Inverse Mori-Tanaka Model

The inverse Mori-Tanaka formulation assumes a high volume fraction of the inclusion and permutes the properties of the matrix and inclusion in a single inclusion problem, which gives

A = {[E : (C_{I}^{- 1} C_{M} - I) + I]}^{- 1} .

Balanced Model

The Mori-Tanaka and inverse Mori-Tanaka formulations give the upper and lower bounds of the composite stiffness for low and high volume fractions. The balanced formulation is the interpolative formulation, which can be written as

H_{1} = {[E : (C_{M}^{- 1} C_{I} - I) + I]}^{- 1},

H_{2} = {[E : (C_{I}^{- 1} C_{M} - I) + I]}^{- 1},

A = {[(1 - φ) H_{1}^{- 1} + φ H_{2}^{- 1}]}^{- 1},

where $φ$ is the smooth interpolation function proposed by Lielens (1999)

φ = \frac{1}{2} v_{f}^{I} (1 + v_{f}^{I}) .

User-Defined Concentration Tensor

You can define the concentration tensor $A$ directly by specifying all of its components. The user-defined concentration tensor must be defined in the local coordinate system of the inclusion. You can use homogenization models other than those listed above, calibration from experimental results, or FE-RVE models to obtain the value of the concentration tensor.

Incremental Mean-Field Homogenization for Nonlinear Composites

An incremental linearization approach is used to extend the homogenization models to composites with nonlinear behavior. The tangent stiffness matrix is used in place of the linear elastic modulus to compute the concentration tensor $A$

Δ ϵ_{I} = A Δ ϵ_{M} .

Since the tangent matrix is not constant in the nonlinear case and is usually a function of the strain, iterations are performed to guarantee a converged solution of the strain increment in each constituent.

Isotropization Methods

Analytical expression of Eshelby's tensor is available only if the inclusions are ellipsoidal and the matrix material is isotropic or transversely isotropic. For matrix materials that are anisotropic, using the stiffness directly to compute the Eshelby's tensor can result in stiffer prediction for the composite behavior (Doghri, 2003). The isotropic part of the tangent stiffness matrix can generally be written as

C^{i s o} = 3 K I^{v o l} + 2 G I^{d e v},

where $K$ is the bulk modulus and $G$ is the shear modulus. $I^{v o l}$ and $I^{d e v}$ are the volumetric and deviatoric part of the fourth-order identity tensor. The method used to extract the isotropic part of the stiffness matrix is not unique, and the accuracy of the prediction varies depending on the properties of the matrix material.

General Method

The general method is the most common isotropization method, and it extracts the isotropic part using

K = \frac{1}{3} I^{v o l} : : C^{a n i},

G = \frac{1}{10} I^{d e v} : : C^{a n i} .

The general method can be used for any anisotropic tangent matrix. This method is the default isotropization method.

Spectral Method

The spectral isotropization method can be used if isotropic hardening is defined in the matrix material (Doghri, 2003). The isotropic part of the tangent matrix is given by

K = K_{e},

G = G_{e} (1 - \frac{3 G_{e}}{3 G_{e} + \frac{d σ^{0}}{d {\bar{ϵ}}^{p l}} ({\bar{ϵ}}^{p l})}),

where $K_{e}$ is the elastic bulk modulus, $G_{e}$ is the elastic shear modulus, $σ^{0} ({\bar{ϵ}}^{p l})$ is the yield stress, and ${\bar{ϵ}}^{p l}$ is the equivalent plastic strain.

Modified Spectral Method

For matrix materials with significant softening after yielding, both the general method and the spectral method might give predictions that are too stiff. The modified spectral method proposed by Selmi et. al. (Selmi, 2011) improves the accuracy of the spectral method by evaluating the derivative of the yield stress at a larger equivalent plastic strain

K = K_{e},

G = k_{G} G_{e} (1 - \frac{3 G_{e}}{3 G_{e} + \frac{d σ^{0}}{d {\bar{ϵ}}^{p l}} (k_{p} {\bar{ϵ}}^{p l} + k_{s})}),

where $k_{G}$ , $k_{p}$ , and $k_{s}$ are the isotropization parameters. You can calibrate these parameters by comparing the mean-field solution to experimental data or to FE-RVE solutions. The isotropization parameters are irrelevant if the constituent does not have plastic behavior.

Isotropic Projection

Options are available to use the isotropic projection of the tangent matrix to compute different parts of the concentration tensor. For better predictions, Doghri et al. (Doghri, 2010) recommended the following:

If the general method is used, only Eshelby's tensor is computed with the isotropic projection.
If the spectral/modified spectral method is used, Hill's tensor given by $P = E : C_{M}^{- 1}$ is computed with the isotropic projection.

It is possible that some choices can lead to unphysical prediction of microlevel strain resulting in unstable material behavior in the composite. In addition, in Abaqus/Standard analyses, the tangent stiffness matrix for the composite might lose its symmetry after homogenization; therefore, you must use the unsymmetric solver to achieve convergence.

Multistep Homogenization

For composites with multiple inclusions, a multistep homogenization approach is used, as shown in Figure 2. The composite is decomposed into "grains," with each grain containing one inclusion family and the matrix. The inclusions in each family have the same material properties, aspect ratio, and orientation. In the first step homogenization is performed in each grain using the user-specified formulation; in the second step the Voigt formulation is used to compute the properties of the overall composite.

An alternative approach is to use the Mori-Tanaka scheme in both the first and the second step, assuming the average strain in the matrix is uniform across all grains. This approach is equivalent to the direct Mori-Tanaka approach proposed by Benveniste (Benveniste et al., 1991). In an Abaqus/Standard analysis, the second approach has the drawback that it might result in an unsymmetric effective modulus when the inclusions are misaligned and nonisotropic. In such cases, use of the unsymmetric solver might be necessary.

Abaqus/Explicit analyses only support one matrix and one inclusion material. For composites with dispersed fiber directions (see Specifying a Second-Order Orientation Tensor in Abaqus/Explicit), the direct Mori-Tanaka approach is used.

Specifying the Microstructure of the Composite

You can specify the microstructure of the composite by giving the shape, volume fraction, aspect ratio, and fiber orientation of each constituent. You can use a distribution to define the volume fraction, aspect ratio, and the orientation tensor (see Distribution Definition).

The orientation tensor can be used only with three-dimensional solid elements.

Specifying the Shape, Volume Fraction, and Aspect Ratio

The volume fraction and aspect ratio can be specified as spatial distributions. The aspect ratio is given by $L / D$ as shown in Figure 3.

The aspect ratio of a spherical inclusion is 1.0, and the aspect ratio of a cylindrical inclusion is infinity; therefore, you do not need to specify the aspect ratio in these cases.

Specifying a Fixed Orientation for the Inclusion and Void Phases

By default, the orientation of the inclusion or void phase is fixed and is aligned with the x-axis of the local coordinate system specified in the section definition. You can specify a different orientation by giving the three components of the orientation vector in the local coordinate system, as shown in Figure 4. The components of the orientation unit vector, $p$ , are

\begin{matrix} p_{1} & = & \sin θ \cos φ, \\ p_{2} & = & \sin θ \sin φ, \\ p_{3} & = & \cos θ . \end{matrix}

The local 3-direction of the inclusion (void included) is the projection of the local z-axis defined in the section onto the surface normal to the inclusion orientation (1-direction). If the local z-axis is within 0.1 degree of being parallel to the 1-direction, the local 3-direction is the projection of the local y-axis onto the surface. The local 2-direction is then at right angles to the local 3-direction, so that the local 2-direction, local 3-direction, and the inclusion orientation (1-direction) form a right-handed set. The local 2 and 3 directions are relevant only if the shape of the inclusion (void included) is an elliptical cylinder.

Specifying a Second-Order Orientation Tensor in Abaqus/Standard

For composites containing short fibers, a second-order orientation tensor, $a$ , is usually used to describe the dispersion of the fiber orientation

a_{i j} = \oint p_{i} p_{j} Ψ (p) d p,

where $\oint d p = \int_{φ = 0}^{2 π} \int_{θ = 0}^{π} \sin θ d θ d φ$ is the integral over the surface of the unit sphere (all possible directions of $p$ ) and $Ψ (p)$ is the orientation distribution function (ODF). The spatial distribution of the orientation tensor is usually given by injection molding software and can be used directly to compute the effective moduli of the composite material that is linear elastic. If the material is nonlinear, the average fields given by $〈 f 〉 = \oint f (p) Ψ (p) d p$ must be computed during the analysis and the ODF is required to calculate the integral. When only the second-order orientation is available, Abaqus recovers the ODF using the approach proposed by Onat and Leckie (1988)

Ψ (p) ≃ \frac{1}{4 π} + \frac{15}{8 π} b_{i j} f_{i j} (p) + \frac{315}{32 π} B_{i j k l} F_{i j k l} (p),

for three-dimensional orientations, where $f_{i j}$ and $F_{i j k l}$ are the tensor basis functions

f_{i j} (p) = p_{i} p_{j} - \frac{1}{3} δ_{i j},

F_{i j k l}^{} = p_{i} p_{j} p_{k} p_{l} - \frac{1}{7} (p_{i} {p_{j} δ}_{k l} + {p_{i} p_{k}}_{} δ_{j l} + p_{i} p_{l} δ_{j k} + p_{k} p_{l} δ_{i j} + p_{j} p_{l} δ_{i k} + p_{j} p_{k} δ_{i l}) + \frac{1}{35} (δ_{i j} δ_{k l} + δ_{i k} δ_{j l} + δ_{i l} δ_{j k}),

and $b_{i j}$ and $B_{i j k l}$ are the deviatoric versions of the orientation tensors

b_{i j} = a_{i j} - \frac{1}{3} δ_{i j},

B_{i j k l}^{} = A_{i j k l} - \frac{1}{7} (a_{i j} δ_{k l} + a_{i k} δ_{j l} + a_{i l} δ_{j k} + a_{k l} δ_{i j} + a_{j l} δ_{i k} + a_{j k} δ_{i l}) + \frac{1}{35} (δ_{i j} δ_{k l} + δ_{i k} δ_{j l} + δ_{i l} δ_{j k}) .

The fourth-order orientation tensor is required during the ODF reconstruction. It is computed using a closure approximation, which is a formula that approximates the fourth-order tensor in terms of the known second-order tensor. A popular closure approximation is the hybrid closure approximation proposed by Advani and Tucker (1987), which is constructed using the interpolation of linear and quadratic closure approximation

A_{i j k l}^{h} = (1 - α) A_{i j k l}^{l} + α A_{i j k l}^{q},

where the linear closure is given by

A_{i j k l}^{l} = - \frac{1}{35} (δ_{i j} δ_{k l} + δ_{i k} δ_{j l} + δ_{i l} δ_{j k}) + \frac{1}{7} (a_{i j} δ_{k l} + a_{i k} δ_{j l} + a_{i l} δ_{j k} + a_{k l} δ_{i j} + a_{j l} δ_{i k} + a_{j k} δ_{i l}),

and the quadratic closure is given by

A_{i j k l}^{q} = a_{i j}^{} a_{k l} .

The weight function is defined as $α = 1 - 27 \det (a)$ for three-dimensional orientations.

Other available closure approximations are the orthotropic smooth approximation (Advani and Tucker, 1987) and the fitted approximation (Cintra and Tucker, 1995).

Only one constituent in the composite is allowed to have orientations specified with the second-order orientation tensor.

Taking into account the symmetry of the ODF $Ψ (p) = Ψ (- p)$ , we compute the integral $\oint f (p) Ψ (p) d p$ numerically as follows:

\oint f (p) Ψ (p) d p ≃ {2 \sum}_{I = 1}^{N} \sum_{I = 1}^{N} f (θ_{I}, φ_{J}) Ψ (θ_{I}, φ_{J}) \sin (θ_{I}) Δ θ Δ φ,

where $Δ θ = Δ φ = \frac{π}{N}$ and $θ_{I} = I Δ θ, φ_{J} = J Δ φ$ . The numerical integral is computed by subdividing the macro point into $N \times N$ subdomains, with each subdomain having the same fiber orientation $p = p (θ_{I}, φ_{J})$ . Then a multistep homogenization (as described in Figure 2) is used to compute the macrolevel solutions. The ODF is more accurately recovered if a large value of N is used; however, the computational cost and memory usage also increases substantially.

When the fiber direction is random in the three-dimensional space, the second-order orientation tensor is $a_{i j} = \frac{1}{3} δ_{i j}$ , and the ODF is constant across all subdomains.

When the direct Mori-Tanaka homogenization approach is used, the fiber response can be calculated either individually in each subdomain or just once using the averaged state over all subdomains. In both cases only the averaged fiber state is available for output. The maximum principal direction of the second-order orientation tensor is used as the 1-direction of the local coordinate system for the inclusion.

Specifying a Second-Order Orientation Tensor in Abaqus/Explicit

In Abaqus/Explicit, when a second-order orientation tensor, $a$ , is used to specify the dispersion of the fiber material, the direct Mori-Tanaka approach is used. Because the fiber material can be only linear elastic with constant moduli, homogenization and strain partitioning can be simplified and no numerical integration is necessary. The following equation is used to compute the orientation averaging terms, $C^{a v g} = \oint C (p) Ψ (p) d p$ :

C_{i j k l}^{a v g} = C_{1} a_{i j k l} + C_{2} (a_{i j} δ_{k l} + a_{k l} δ_{i j}) + C_{3} (a_{i k} δ_{j l} + a_{i l} δ_{j k} + a_{j l} δ_{i k} + a_{j k} δ_{i l}) + C_{4} δ_{i j} δ_{k l} + C_{5} (δ_{i k} δ_{j l} + δ_{i l} δ_{j k}),

where $C_{1}$ , $C_{2}$ , $C_{3}$ , $C_{4}$ , and $C_{5}$ are the five constants of the tensor that is transversely isotropic. ODF reconstruction is not necessary and is not performed.

The second-order orientation tensor is specified in the orientation definition of the section where the multiscale material is assigned.

Composites with Thermal Expansion

Composites with thermal expansion are supported only in Abaqus/Standard.

A three-step approach is used to compute the macro response of a thermoelastic composite (Pierard, 2004). The total macro strain is applied to the composite in the first step with no temperature change. The strain in each constituent is computed with the concentration tensor, and the stress is updated accordingly. In Step 2 an equal temperature change is applied to all constituents; by assuming uniform stress increment in all constituents, the strain increment in the second step is computed in each constituent. In Step 3 an equal and opposite strain increment computed in Step 2 is applied to all constituents to recover the total macro strain; there is no temperature change in this step. The total strain and stress in each constituent is computed by superposing the solutions from all three steps. This approach is also applied to the incremental formulation of nonlinear composites with thermal expansion.

Thermal expansion is supported only for three-dimensional solid elements.

Composites with Viscoelasticity

Composites with viscoelasticity are supported only in Abaqus/Standard.

When one or more constituents have viscoelasticity defined in the time domain, the homogenization approach described above is no longer appropriate because it does not take into account the relaxation of stress with time. However, similar homogenization can be performed in the frequency domain, and homogenized properties of the composite in the time domain can be computed using homogenized properties at different frequencies. This homogenization approach is applicable only when the viscoelasticity in the constituents is defined in the time domain using a Prony series expansion and includes the following steps:

The shear and bulk moduli of all the constituents are transformed to the frequency domain using Fourier transforms (see details in Prony Series Parameters).
The mean-field homogenization is performed in the frequency domain the same way as discussed above, with the exception that the moduli are complex numbers instead of real numbers. The real part of the composite modulus is the storage modulus, $C^{s}$ , and the imaginary part is the loss modulus, $C^{l}$ .
Assuming the time response of the composite can be written in a similar form as that in the constituents where viscoelasticity is defined using a Prony series expansion
$g_{i j} (t) = 1 - \sum_{k = 1, N} g_{k} (1 - e^{- t / τ_{k}}),$
where $g_{i j} (t) = C_{i j} (t) / C_{i j}^{0}$ and $C_{i j}^{0}$ is the instantaneous elastic modulus of the composite, the homogenized storage modulus of the composite, $C^{s}$ , can be written in the frequency domain as follows:
$C_{i j}^{s} (ω) = C_{i j}^{0} [1 - \sum_{k = 1, N} \frac{g_{k}}{1 + τ_{k}^{2} ω^{2}}] .$
The parameters $τ_{k}$ are selected to be equally spaced numbers on a logarithmic scale in the range $[τ_{\min}, τ_{\max}]$ . $τ_{\min}$ and $τ_{\max}$ are the minimum and maximum frequencies used to define the viscoelasticity in all constituents. The parameters $g_{k}$ are computed through a least-squares curve fit at frequencies $ω_{l} = 1 / τ_{l}$ $(l = 1, N)$ .

The thermorheologically simple temperature effects are not supported.

Other limitations include:

The inverse Mori-Tanaka and balanced homogenization formulations are not supported when you use the orientation tensor to define the direction of the inclusion.
When you use the orientation tensor to define the direction of the inclusion, you can define only transversely isotropic elasticity in the inclusion material, and no other inclusions can be defined in the composite.
You cannot use the orientation tensor to define the direction of the inclusion in which viscoelasticity is defined.
Distributions other than distributions of the orientation tensor are not allowed in the material definition.
Temperature and field variable dependencies are allowed to define the elasticity of the matrix material only, and only one inclusion is allowed in this case.

Because the multiscale material is replaced by a composite with homogenized viscoelastic properties, the microlevel output is not available.

Composites with Damage

When damage and failure are defined at the constituent level, the damage in each constituent contributes to the overall damage in an indirect way through stress averaging and strain partitioning. As the stiffness of the constituent decreases, the strain increment in this constituent is likely to increase based on the formulation of the homogenization. In reality, the damage behavior of each constituent is likely to interact with each other and affect the overall damage behavior of the composite; therefore, you might need to specify an additional damage variable, $D$ , for the composite. At any time during the analysis, the stress tensor of the composite is given by

σ = (1 - D) \bar{σ},

where $\bar{σ} = \sum_{i} v_{i}^{f} σ_{i}$ , and $σ_{i} = (1 - D_{i}) {\bar{σ}}_{i}$ is the effective stress of the constituent.

You can choose to combine the damage variables, $D_{i}$ , in all constituents to get the overall damage behavior of the composite. The overall damage variable, $D$ , can be computed as the maximum of $D_{i}$ :

D = \max {D_{i}},

or you can choose to combine the damage variables in a multiplicative sense:

D = 1 - \prod_{i} (1 - D_{i}) .

You can also specify $D$ using user subroutine UDAMAGEMF. In addition, you can specify a separate damage variable, $d_{i}^{int}$ , for each constituent. This variable is used to account for the additional contribution to the overall damage behavior due to interaction (for example, a fiber material might buckle more easily if cracks start to develop in the matrix material) as follows:

σ = (1 - D) \sum_{i} v_{i}^{f} (1 - d_{i}^{int}) σ_{i} .

Example

This simple example illustrates how to define a multiscale material.

The name of the multiscale material is COMPOSITE. The material consists of one inclusion material MAT2 embedded in the matrix material MAT1. The shape of the inclusion is prolate. The Mori-Tanaka homogenization method is used. The inclusion has a volume fraction of $v_{f}^{}$ and an aspect ratio of $a r$ . The direction of the inclusion is fixed and defined with a vector $(p_{1}, p_{2}, p_{3})$ . The name of the matrix is MATRIX_MAT1, and the name of the inclusion is FIBER_MAT2.

Input file template

MATERIAL, NAME=MAT1
…
MATERIAL, NAME=MAT2
…
MATERIAL, NAME=COMPOSITE
MEAN FIELD HOMOGENIZATION, FORMULATION=MT
CONSTITUENT, TYPE=MATRIX, MATERIAL=MAT1, NAME=MATRIX_MAT1
CONSTITUENT, TYPE=INCLUSION, MATERIAL=MAT2, NAME=FIBER_MAT2, SHAPE=PROLATE, DIRECTION=FIXED
 $v_{f}^{}$ ,  $a r$ ,  $p_{1}$ ,  $p_{2}$ ,  $p_{3}$

Elements

In Abaqus/Standard the multiscale material model can be used with

three-dimensional solid elements, including stress/displacement elements, diffusive heat transfer elements, and coupled temperature-displacement elements;
two-dimensional solid elements (stress/displacement only);
membrane elements;
three-dimensional conventional shell elements (stress/displacement only); and
continuum shell elements (stress/displacement only).

In Abaqus/Standard additional plane stress iterations are included for two-dimensional plane stress elements and shell elements.

In Abaqus/Explicit the multiscale material model can be used only with three-dimensional solid stress/displacement elements.

Hybrid elements and incompatible mode elements are not supported.

Output

By default, the output variables are macrolevel results. In addition, you can also request microlevel output. The name of the constituent is appended to the end of the output variable.

For composites reinforced with dispersed short fibers, the fiber output variables are computed in the local system with the local fiber direction as the maximum principal direction of the orientation tensor. The matrix output is available only if the direct Mori-Tanaka homogenization approach is used. Otherwise, only averaged stress and strain are available. If plasticity is specified in the matrix material, averaged plastic strains are also available. The same matrix output limitation also applies to composites reinforced with multiple fiber families.

Microlevel strain output is not available in Abaqus/Explicit.

References

Advani, S., and C. Tucker, “The Use of Tensors to Describe and Predict Fiber Orientation in Short Fiber Composites,” Journal of Rheology, vol. 31 (8), pp. 751–784, 1987.
Benveniste, Y., G. J. Dvorak, and T. Chen, “On Diagonal and Elastic Symmetry of the Approximate Effective Stiffness Tensor of Heterogeneous Media,” Journal of the Mechanics and Physics of Solids, vol. 39(7), pp. 927–946, 1991.
Cintra, J., and C. Tucker, “Orthotropic Closure Approximations for Flow-Induced Fiber Orientation,” Journal of Rheology, vol. 39, pp. 1095–1122, 1995.
Doghri, I., I. Adam, and N. Bilger, “Mean-Field Homogenization of Elasto-Viscoplastic Composites Based on a General Incrementally Affine Linearization Method,” International Journal of Plasticity, pp. 219–238, 2010.
Doghri, I., and A. Ouaar, “Homogenization of Two-Phase Elasto-Plastic Composite Materials and Structures. Study of Tangent Operators, Cyclic Plasticity and Numerical Algorithms,” International Journal of Solids and Structures, vol. 40(7), pp. 1681–1712, 2003.
Eshelby, J. D., “The Determination of the Elastic Field of an Ellipsoidal Inclusion and Related Problems,” Proceedings of the Royal Society of London, pp. 376–396, 1957.
Kammoun, S., L. Brassart, G. Robert, I. Doghri, and L. Delannay, “Micromechanical Modeling of Short Glass-Fiber Reinforced Thermoplastics-Isotropic Damage of Pseudograins,” American Institute of Physics Conference Proceedings 1353, pp. 974–977, 2011.
Lielens, G., “Micro-Macro Modeling of Structured Materials,” PhD Thesis Universite Catholique de Louvain Belgium, 1999.
Mura, T., Micromechanics of Defects in Solids, MARTINUS NIJHOFF Publishers, 1987.
Onat, E.T., and F.A. Leckie, “Representation of Mechanical Behavior in the Presence of Changing Internal Structure,” Journal of Applied Mechanics, vol. 55(1), pp. 1–10, 1988.
Pierard, O., and I. Doghri, “An Enhanced Affine Formulation and the Corresponding Numerical Algorithms for the Mean-Field Homogenization of Elasto-Viscoplastic Composites,” International Journal of Plasticity, pp. 131–157, 2006.
Pierard, O., C. Friebel, and I. Doghri, “Mean-Field Homogenization of Multi-Phase Thermo-Elastic Composites: a General Framework and its Validation,” Composites Science and Technology, vol. 64, pp. 1587–1603, 2004.
Selmi, A., I. Doghri, and L. Adam, “Micromechanical Simulations of Biaxial Yield, Hardening and Plastic Flow in Short Glass Fiber Reinforced Polyamide,” International Journal of Mechanical Sciences, vol. 53(9), pp. 696–706, 2011.