Mori Tanaka: Why symmetrize?

The following notebook is intended to explain the various Mori-Tanaka implementations and solution strategies used in HomoPy. The goal is not to derive the equations, but to provide a learning by doing example. A prerequisite is a basic understanding of homogenization methods in the field of elasticity.

The starting point is the formulation after Benveniste (1987), which can also be found in Brylka (2017)

\[\bar{\mathbb{C}}^{\text{MT}} = \mathbb{C}_m + c_f \left<(\mathbb{C}_{f,\alpha} - \mathbb{C}_m)\mathbb{A}^{\text{SIP}}_{f,\alpha}\right>_f (c_m\mathbb{I}^\text{S}+c_f\left<\mathbb{A}^\text{SIP}_{f,\alpha}\right>_f)^{-1},\]

where \(\bar{\mathbb{C}}^{\text{MT}}\) is the effective stiffness, \(\mathbb{C}_m\) is the matrix stiffness, \(\mathbb{C}_{f,\alpha}\) is the stiffness of the inclusion of type \(\alpha\), \(c_m\) and \(c_f\) are the volume fractions of matrix and fiber, respectively, \(\mathbb{A}^{\text{SIP}}_{f,\alpha}\) is the strain localization tensor based on Eshelby’s solution for the inclusion shape of type \(\alpha\) and \(\mathbb{I}^\text{S}\) is the symmetric identity of order four. The brackets \(\left< \cdot \right>_f\) indicate the fiber volume average.

The problem with this formulation is that depending on the different shapes of inclusions and elastic symmetries (isotropy, transverse isotropy, …) of the constituents, the effective stiffness is not guaranteed to be major symmetric, which violates thermodynamic principles.

To demonstrate this, we start with a simple example:

First, we assume a single inclusion type, namely carbon fibers aligned with the x-axis in a polyamide-6 matrix. For simplicity, the carbon fiber is assumed to be isotropic, while in reality it should be transverse isotropic.

[1]:

import numpy as np

np.set_printoptions(linewidth=100, precision=4)

from homopy.methods import *
from homopy.elasticity import *
from homopy.stiffness_plot import *

# define fiber and matrix properties
carbon_fiber = Isotropy(242e9, 0.1)
v_frac_carbon = 0.25
a_carbon = 347
polyamid6 = Isotropy(1.18e9, 0.35)

# define the MT homogenization and print it's effective stiffness
mt_carbon = mt_carbon = MoriTanaka(
    polyamid6, carbon_fiber, v_frac_carbon, a_carbon, shape="ellipsoid"
)
print(mt_carbon.effective_stiffness66)

[[6.0989e+10 1.1464e+09 1.1464e+09 0.0000e+00 0.0000e+00 0.0000e+00]
 [1.1464e+09 2.7483e+09 1.4048e+09 0.0000e+00 0.0000e+00 0.0000e+00]
 [1.1464e+09 1.4048e+09 2.7483e+09 0.0000e+00 0.0000e+00 0.0000e+00]
 [0.0000e+00 0.0000e+00 0.0000e+00 1.3435e+09 0.0000e+00 0.0000e+00]
 [0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 1.4507e+09 0.0000e+00]
 [0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 1.4507e+09]]

The result is a symmetric, i.e. thermodynamically consistent, transversely isotropic effective stiffness. The symmetry can also be validated by using

[2]:

mt_carbon.is_symmetric()

Left minor symmetry: passed
Right minor symmetry: passed
Major symmetry: passed

So far, there is nothing to complain about. This is due to the fact, that the Mori Tanaka formulation always results in a symmetric stiffness, when the matrix is isotropic, and the single inclusion is isotropic.

When we add glass fibers (or any other type of inclusion in that sense) to the mix, we suddenly get a different behaviour, as the next example will demonstrate.

[3]:

# define fiberproperties
glass_fiber = Isotropy(80e9, 0.22)
v_frac_glass = 0.25
a_glass = 225

# define the MT homogenization and print it's effective stiffness
mt_hybrid = MoriTanaka(
    polyamid6,
    [carbon_fiber, glass_fiber],
    [v_frac_carbon / 2, v_frac_glass / 2],
    [a_carbon, a_glass],
    2 * ["ellipsoid"],
)
print(mt_hybrid.effective_stiffness66)
print("\n")
mt_hybrid.is_symmetric()

[[4.1198e+10 1.2154e+09 1.2154e+09 0.0000e+00 0.0000e+00 0.0000e+00]
 [1.2155e+09 2.7383e+09 1.3995e+09 0.0000e+00 0.0000e+00 0.0000e+00]
 [1.2155e+09 1.3995e+09 2.7383e+09 0.0000e+00 0.0000e+00 0.0000e+00]
 [0.0000e+00 0.0000e+00 0.0000e+00 1.3388e+09 0.0000e+00 0.0000e+00]
 [0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 1.4436e+09 0.0000e+00]
 [0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 1.4436e+09]]


Left minor symmetry: passed
Right minor symmetry: passed
Major symmetry: failed
The rel. residuum for major sym. is: res = 6.406621895094511e-06

We see that the effective stiffness does not contain the major symmetry and thus is not thermodynamically consistent. The relative residuum indicates that the major symmetry is only violated ‘a bit’, but it is violated. Therefore, a motivation is given to use an alternative algorithm which ensures that the effective stiffness will be symmetric.

By abbreviating the classic formulation above with

\[\bar{\mathbb{C}}^{\text{MT}} = \mathbb{C}_m + c_f \underbrace{\left<(\mathbb{C}_{f,\alpha} - \mathbb{C}_m)\mathbb{A}^{\text{SIP}}_{f,\alpha}\right>_f}_{\mathbb{X}} \underbrace{\left(c_m\mathbb{I}^\text{S}+c_f\left<\mathbb{A}^\text{SIP}_{f,\alpha}\right>_f\right)^{-1}}_{\mathbb{Y}^{-1}},\]

Segura et al. (2023) uses the relation

\[\mathbb{X}\mathbb{Y}^{-1} = \left(\mathbb{Y}\mathbb{X}^{-1}\right)^{-1}\]

to derive the alternative form

\[\bar{\mathbb{C}}^{\text{MT}} = \mathbb{C}_m + c_f \left( c_m\left<(\mathbb{C}_{f,\alpha} - \mathbb{C}_m)\mathbb{A}^{\text{SIP}}_{f,\alpha}\right>_f^{-1} +c_f\underbrace{\left<\mathbb{A}^\text{SIP}_{f,\alpha}\right>_f \left<(\mathbb{C}_{f,\alpha} - \mathbb{C}_m)\mathbb{A}^{\text{SIP}}_{f,\alpha}\right>_f^{-1}}_\mathbb{Z}\right)^{-1},\]

where the left term in the paranthesis is always symmetric. Therefore, only the right term, i.e. \(\mathbb{Z}\), needs to be symmetrized, to recover the major symmetry, i.e.

\[\bar{\mathbb{C}}^{\text{MT}} = \mathbb{C}_m + c_f \left( c_m\left<(\mathbb{C}_{f,\alpha} - \mathbb{C}_m)\mathbb{A}^{\text{SIP}}_{f,\alpha}\right>_f^{-1} +\frac{c_f}{2}(\mathbb{Z}+\mathbb{Z}^\mathsf{T})\right)^{-1}.\]

HomoPy uses this formulation, when the ‘symmetrize’ flag is set to True

[4]:

mt_hybrid_sym = MoriTanaka(
    polyamid6,
    [carbon_fiber, glass_fiber],
    [v_frac_carbon / 2, v_frac_glass / 2],
    [a_carbon, a_glass],
    2 * ["ellipsoid"],
    symmetrize=True,
)
print(mt_hybrid_sym.effective_stiffness66)
print("\n")
mt_hybrid_sym.is_symmetric()

[[4.1198e+10 1.2155e+09 1.2155e+09 0.0000e+00 0.0000e+00 0.0000e+00]
 [1.2155e+09 2.7383e+09 1.3995e+09 0.0000e+00 0.0000e+00 0.0000e+00]
 [1.2155e+09 1.3995e+09 2.7383e+09 0.0000e+00 0.0000e+00 0.0000e+00]
 [0.0000e+00 0.0000e+00 0.0000e+00 1.3388e+09 0.0000e+00 0.0000e+00]
 [0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 1.4436e+09 0.0000e+00]
 [0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 1.4436e+09]]


Left minor symmetry: passed
Right minor symmetry: passed
Major symmetry: passed