Note
Go to the end to download the full example code.
Composite Parameter Identification using Mori-Tanaka and Self-Consistent
This example demonstrates inverse identification of reinforcement elastic properties in a two-phase composite material. Given experimental effective Young’s modulus data at various volume fractions (Wang 2003), we identify the Young’s modulus and Poisson’s ratio of the reinforcement phase.
Two micromechanical schemes are compared:
Mori-Tanaka (
MIMTN): accurate at low-to-moderate volume fractionsSelf-Consistent (
MISCN): better at high volume fractions (c > 0.3)
Forward model: simcoon mean-field homogenization (L_eff)
Optimization: scipy.optimize.differential_evolution (global optimizer)
Phases: described in memory, so an evaluation touches no file at all
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import differential_evolution
import simcoon as sim
from simcoon.solver.micromechanics import Ellipsoid
import os
Problem Setup
We have a glass-particle / epoxy composite. The matrix properties are known:
\(E_m = 2250\) MPa, \(\\nu_m = 0.19\)
The reinforcement properties \((E_f, \\nu_f)\) are unknown and must be identified from experimental effective Young’s modulus measurements at several volume fractions.
The phases are built as objects and passed to L_eff. A candidate
\((E_f, \\nu_f)\) is written into the reinforcement,
Note that simcoon’s key system (simcoon.parameter) remains the tool of choice when
the forward model is an external code driven by input files, such as an
Abaqus job; it is not needed when simcoon itself is the forward model.
def load_experimental_data(filepath):
"""Load experimental E_eff vs volume fraction data."""
data = np.loadtxt(filepath)
return data[:, 0], data[:, 1]
Forward Model
For each set of candidate properties \((E_f, \\nu_f)\), we compute the
effective Young’s modulus at all experimental volume fractions using simcoon’s
L_eff function.
def compute_E_eff(E_f, nu_f, concentrations, umat_name, props_composite):
"""
Compute effective Young's modulus at given volume fractions.
Parameters
----------
E_f : float
Reinforcement Young's modulus [MPa]
nu_f : float
Reinforcement Poisson's ratio [-]
concentrations : array-like
Volume fractions of reinforcement
umat_name : str
Micromechanical scheme ("MIMTN" or "MISCN")
props_composite : numpy.ndarray
Composite definition array [nphases, unused, int1, int2, n_matrix]
Returns
-------
numpy.ndarray
Effective Young's modulus at each volume fraction [MPa]
"""
matrix = Ellipsoid(
number=0,
umat_name="ELISO",
save=1,
nstatev=1,
props=np.array([2250.0, 0.19, 8.8e-5]),
)
reinforcement = Ellipsoid(
number=1,
umat_name="ELISO",
save=1,
nstatev=1,
props=np.array([E_f, nu_f, 0.5e-6]),
)
E_eff = np.zeros(len(concentrations))
for i, c in enumerate(concentrations):
matrix.concentration = 1.0 - c
reinforcement.concentration = c
L = sim.L_eff(
umat_name,
props_composite,
0,
phases=[matrix, reinforcement],
)
iso_props = sim.L_iso_props(L).flatten()
E_eff[i] = iso_props[0]
return E_eff
def cost_function(params_opt, c_exp, E_exp, umat_name, props_composite):
"""MSE cost function for reinforcement property identification."""
E_f, nu_f = params_opt
try:
E_pred = compute_E_eff(E_f, nu_f, c_exp, umat_name, props_composite)
return np.mean((E_pred - E_exp) ** 2)
except Exception:
return 1e12
def identify_reinforcement(
c_exp,
E_exp,
umat_name,
props_composite,
bounds,
verbose=True,
):
"""
Identify reinforcement properties using differential evolution.
Returns
-------
dict
Identified parameters and optimization result
"""
if verbose:
print(f"\n{'=' * 60}")
print(f" Identification with {umat_name}")
print(f"{'=' * 60}")
result = differential_evolution(
cost_function,
bounds=bounds,
args=(c_exp, E_exp, umat_name, props_composite),
strategy="best1bin",
maxiter=200,
popsize=15,
tol=1e-8,
mutation=(0.5, 1.0),
recombination=0.7,
seed=42,
polish=True,
disp=verbose,
)
E_f_opt, nu_f_opt = result.x
if verbose:
print(f"\n E_f = {E_f_opt:.0f} MPa")
print(f" nu_f = {nu_f_opt:.3f}")
print(f" MSE = {result.fun:.2e} MPa^2")
return {
"E_f": E_f_opt,
"nu_f": nu_f_opt,
"mse": result.fun,
"result": result,
"umat_name": umat_name,
}
Main Execution
We compare two identification strategies:
Mori-Tanaka (
MIMTN): dilute approximation, best for c < 30%Self-Consistent (
MISCN): accounts for percolation, better at high c
At high volume fractions (c = 0.5), Mori-Tanaka underestimates stiffness and compensates by overestimating \(E_f\). Self-Consistent gives more physically meaningful identified properties.
if __name__ == "__main__":
try:
script_dir = os.path.dirname(os.path.abspath(__file__))
except NameError:
script_dir = os.getcwd()
# -----------------------------------------------------------------
# Load experimental data (Wang 2003)
# -----------------------------------------------------------------
c_exp, E_exp = load_experimental_data(os.path.join(script_dir, "data", "E_exp.txt"))
print("=" * 60)
print(" COMPOSITE REINFORCEMENT IDENTIFICATION")
print(" Glass/Epoxy — Mean-field + Differential Evolution")
print("=" * 60)
print(f"\nExperimental data (Wang 2003): {len(c_exp)} points")
for c, E in zip(c_exp, E_exp):
print(f" c = {c:.1f} -> E_eff = {E:.0f} MPa")
# Composite definition: [int1, int2, n_matrix] (Eshelby integration points, matrix index)
props_composite = np.array([50, 50, 0], dtype="float")
# Bounds for identification (E_f, nu_f)
bounds = [(10000, 200000), (0.01, 0.45)]
# -----------------------------------------------------------------
# Identification with Mori-Tanaka
# -----------------------------------------------------------------
result_MT = identify_reinforcement(
c_exp,
E_exp,
"MIMTN",
props_composite,
bounds,
)
# -----------------------------------------------------------------
# Identification with Self-Consistent
# -----------------------------------------------------------------
result_SC = identify_reinforcement(
c_exp,
E_exp,
"MISCN",
props_composite,
bounds,
)
# -----------------------------------------------------------------
# Summary
# -----------------------------------------------------------------
print(f"\n{'=' * 60}")
print(" SUMMARY")
print(f"{'=' * 60}")
print(" Reference glass: E_f ~ 73000 MPa, nu_f ~ 0.22")
print("")
print(f" {'Scheme':<20} {'E_f [MPa]':>12} {'nu_f':>8} {'MSE':>14}")
print(f" {'-' * 56}")
print(
f" {'Mori-Tanaka':<20} {result_MT['E_f']:>12.0f} "
f"{result_MT['nu_f']:>8.3f} {result_MT['mse']:>14.2e}"
)
print(
f" {'Self-Consistent':<20} {result_SC['E_f']:>12.0f} "
f"{result_SC['nu_f']:>8.3f} {result_SC['mse']:>14.2e}"
)
print()
# Note on model limitations
print(" Note: Mori-Tanaka underestimates stiffness at high volume")
print(" fractions and compensates by overestimating E_f. Self-Consistent")
print(" accounts for phase connectivity and gives more physically")
print(" meaningful results for this dataset (c up to 50%).")
# -----------------------------------------------------------------
# Visualization
# -----------------------------------------------------------------
c_model = np.arange(0.0, 0.51, 0.01)
E_MT = compute_E_eff(
result_MT["E_f"],
result_MT["nu_f"],
c_model,
"MIMTN",
props_composite,
)
E_SC = compute_E_eff(
result_SC["E_f"],
result_SC["nu_f"],
c_model,
"MISCN",
props_composite,
)
# Reference with handbook glass properties
E_ref_MT = compute_E_eff(73000, 0.22, c_model, "MIMTN", props_composite)
E_ref_SC = compute_E_eff(73000, 0.22, c_model, "MISCN", props_composite)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
# Left: Identified fits
ax1.plot(
c_exp, E_exp, "kx", markersize=10, markeredgewidth=2, label="Exp. (Wang 2003)"
)
ax1.plot(
c_model,
E_MT,
"b-",
linewidth=2,
label=f"MT identified (E_f={result_MT['E_f']:.0f})",
)
ax1.plot(
c_model,
E_SC,
"r-",
linewidth=2,
label=f"SC identified (E_f={result_SC['E_f']:.0f})",
)
ax1.set_xlabel("Reinforcement volume fraction $c$", fontsize=12)
ax1.set_ylabel("Effective Young's modulus $E_{eff}$ [MPa]", fontsize=12)
ax1.set_title("Identified Reinforcement Properties", fontsize=13)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
# Right: Reference vs Identified
ax2.plot(
c_exp, E_exp, "kx", markersize=10, markeredgewidth=2, label="Exp. (Wang 2003)"
)
ax2.plot(
c_model,
E_ref_MT,
"b--",
linewidth=1.5,
alpha=0.6,
label="MT with E_f=73 GPa (handbook)",
)
ax2.plot(
c_model,
E_ref_SC,
"r--",
linewidth=1.5,
alpha=0.6,
label="SC with E_f=73 GPa (handbook)",
)
ax2.plot(c_model, E_MT, "b-", linewidth=2, alpha=0.8, label="MT identified")
ax2.plot(c_model, E_SC, "r-", linewidth=2, alpha=0.8, label="SC identified")
ax2.set_xlabel("Reinforcement volume fraction $c$", fontsize=12)
ax2.set_ylabel("Effective Young's modulus $E_{eff}$ [MPa]", fontsize=12)
ax2.set_title("Handbook vs Identified Properties", fontsize=13)
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)
fig.suptitle(
"Composite Reinforcement Identification — Glass/Epoxy\n"
"Mori-Tanaka vs Self-Consistent + Differential Evolution",
fontsize=14,
fontweight="bold",
)
plt.tight_layout()
plt.show()
============================================================
COMPOSITE REINFORCEMENT IDENTIFICATION
Glass/Epoxy — Mean-field + Differential Evolution
============================================================
Experimental data (Wang 2003): 6 points
c = 0.0 -> E_eff = 2250 MPa
c = 0.1 -> E_eff = 2580 MPa
c = 0.2 -> E_eff = 3390 MPa
c = 0.3 -> E_eff = 4480 MPa
c = 0.4 -> E_eff = 6310 MPa
c = 0.5 -> E_eff = 10500 MPa
============================================================
Identification with MIMTN
============================================================
differential_evolution step 1: f(x)= 2845441.8584308363
differential_evolution step 2: f(x)= 2845441.8584308363
differential_evolution step 3: f(x)= 2844384.527162666
differential_evolution step 4: f(x)= 2843372.44679017
differential_evolution step 5: f(x)= 2843130.1492564585
differential_evolution step 6: f(x)= 2842473.3377612284
differential_evolution step 7: f(x)= 2842473.3377612284
differential_evolution step 8: f(x)= 2841803.683566658
differential_evolution step 9: f(x)= 2841691.1728336997
differential_evolution step 10: f(x)= 2841691.1728336997
differential_evolution step 11: f(x)= 2841647.8003525925
differential_evolution step 12: f(x)= 2841602.7461820096
differential_evolution step 13: f(x)= 2841453.417909387
differential_evolution step 14: f(x)= 2841453.417909387
differential_evolution step 15: f(x)= 2841453.417909387
differential_evolution step 16: f(x)= 2841453.417909387
differential_evolution step 17: f(x)= 2841419.1478061364
differential_evolution step 18: f(x)= 2841401.2186515722
differential_evolution step 19: f(x)= 2841401.2186515722
differential_evolution step 20: f(x)= 2841401.2186515722
differential_evolution step 21: f(x)= 2841401.2186515722
differential_evolution step 22: f(x)= 2841401.2186515722
differential_evolution step 23: f(x)= 2841401.2186515722
differential_evolution step 24: f(x)= 2841396.108452113
differential_evolution step 25: f(x)= 2841396.108452113
differential_evolution step 26: f(x)= 2841396.108452113
differential_evolution step 27: f(x)= 2841393.3240744397
differential_evolution step 28: f(x)= 2841393.3240744397
differential_evolution step 29: f(x)= 2841393.3240744397
differential_evolution step 30: f(x)= 2841391.3966782917
differential_evolution step 31: f(x)= 2841391.3966782917
differential_evolution step 32: f(x)= 2841391.3966782917
differential_evolution step 33: f(x)= 2841391.3966782917
differential_evolution step 34: f(x)= 2841391.1774388435
differential_evolution step 35: f(x)= 2841390.819151416
differential_evolution step 36: f(x)= 2841390.2110353014
differential_evolution step 37: f(x)= 2841390.2110353014
differential_evolution step 38: f(x)= 2841390.2110353014
differential_evolution step 39: f(x)= 2841390.2110353014
differential_evolution step 40: f(x)= 2841390.2110353014
differential_evolution step 41: f(x)= 2841390.2110353014
differential_evolution step 42: f(x)= 2841390.2110353014
differential_evolution step 43: f(x)= 2841390.161741854
differential_evolution step 44: f(x)= 2841390.101437967
differential_evolution step 45: f(x)= 2841390.101437967
differential_evolution step 46: f(x)= 2841390.101437967
differential_evolution step 47: f(x)= 2841390.101437967
differential_evolution step 48: f(x)= 2841390.101437967
differential_evolution step 49: f(x)= 2841390.078150217
differential_evolution step 50: f(x)= 2841390.07126537
differential_evolution step 51: f(x)= 2841390.07126537
differential_evolution step 52: f(x)= 2841390.065127143
differential_evolution step 53: f(x)= 2841390.0629338943
differential_evolution step 54: f(x)= 2841390.0629338943
differential_evolution step 55: f(x)= 2841390.0629338943
differential_evolution step 56: f(x)= 2841390.060185105
Polishing solution with 'L-BFGS-B'
E_f = 200000 MPa
nu_f = 0.450
MSE = 2.84e+06 MPa^2
============================================================
Identification with MISCN
============================================================
differential_evolution step 1: f(x)= 54451.33267154567
differential_evolution step 2: f(x)= 52713.18502556835
differential_evolution step 3: f(x)= 52250.080494583985
differential_evolution step 4: f(x)= 52250.080494583985
differential_evolution step 5: f(x)= 51026.38867824256
differential_evolution step 6: f(x)= 50844.447549089666
differential_evolution step 7: f(x)= 50844.447549089666
differential_evolution step 8: f(x)= 50844.447549089666
differential_evolution step 9: f(x)= 50784.99401181517
differential_evolution step 10: f(x)= 50638.77745195191
differential_evolution step 11: f(x)= 50638.77745195191
differential_evolution step 12: f(x)= 50638.77745195191
differential_evolution step 13: f(x)= 50638.77745195191
differential_evolution step 14: f(x)= 50623.99454003979
differential_evolution step 15: f(x)= 50623.783659892244
differential_evolution step 16: f(x)= 50590.658213934024
differential_evolution step 17: f(x)= 50590.658213934024
differential_evolution step 18: f(x)= 50588.31424864657
differential_evolution step 19: f(x)= 50584.76160028943
differential_evolution step 20: f(x)= 50584.76160028943
differential_evolution step 21: f(x)= 50583.00686011609
differential_evolution step 22: f(x)= 50583.00686011609
differential_evolution step 23: f(x)= 50582.40410017907
differential_evolution step 24: f(x)= 50582.40410017907
differential_evolution step 25: f(x)= 50582.40410017907
differential_evolution step 26: f(x)= 50582.1659720676
differential_evolution step 27: f(x)= 50582.11900479212
differential_evolution step 28: f(x)= 50581.99126579074
differential_evolution step 29: f(x)= 50581.99126579074
differential_evolution step 30: f(x)= 50581.88530440981
differential_evolution step 31: f(x)= 50581.88530440981
differential_evolution step 32: f(x)= 50581.88530440981
differential_evolution step 33: f(x)= 50581.88530440981
differential_evolution step 34: f(x)= 50581.88530440981
differential_evolution step 35: f(x)= 50581.88530440981
differential_evolution step 36: f(x)= 50581.88530440981
differential_evolution step 37: f(x)= 50581.88530440981
differential_evolution step 38: f(x)= 50581.88153485019
differential_evolution step 39: f(x)= 50581.88153485019
differential_evolution step 40: f(x)= 50581.88153485019
differential_evolution step 41: f(x)= 50581.881276461645
differential_evolution step 42: f(x)= 50581.88110700669
differential_evolution step 43: f(x)= 50581.88091661406
differential_evolution step 44: f(x)= 50581.8802449384
differential_evolution step 45: f(x)= 50581.8802449384
differential_evolution step 46: f(x)= 50581.88024045628
differential_evolution step 47: f(x)= 50581.88024045628
differential_evolution step 48: f(x)= 50581.88023945281
Polishing solution with 'L-BFGS-B'
E_f = 46367 MPa
nu_f = 0.010
MSE = 5.06e+04 MPa^2
============================================================
SUMMARY
============================================================
Reference glass: E_f ~ 73000 MPa, nu_f ~ 0.22
Scheme E_f [MPa] nu_f MSE
--------------------------------------------------------
Mori-Tanaka 200000 0.450 2.84e+06
Self-Consistent 46367 0.010 5.06e+04
Note: Mori-Tanaka underestimates stiffness at high volume
fractions and compensates by overestimating E_f. Self-Consistent
accounts for phase connectivity and gives more physically
meaningful results for this dataset (c up to 50%).
Total running time of the script: (1 minutes 57.594 seconds)