.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/heterogeneous/composite_parameter_identification.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_heterogeneous_composite_parameter_identification.py: 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 fractions - **Self-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) **Key system**: simcoon ``Parameter`` for generic file templating .. GENERATED FROM PYTHON SOURCE LINES 19-27 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt from scipy.optimize import differential_evolution import simcoon as sim from simcoon import parameter as par import os .. GENERATED FROM PYTHON SOURCE LINES 28-43 Problem Setup ------------- We have a glass-particle / epoxy composite. The matrix properties are known: - :math:`E_m = 2250` MPa, :math:`\\nu_m = 0.19` The reinforcement properties :math:`(E_f, \\nu_f)` are unknown and must be identified from experimental effective Young's modulus measurements at several volume fractions. **Key System**: simcoon's key system allows parameters to be injected into any input file via alphanumeric placeholders (e.g., ``@2p``). This makes the identification workflow generic and applicable to any simulation tool (simcoon solver, Mori-Tanaka, fedoo FE, etc.). .. GENERATED FROM PYTHON SOURCE LINES 43-51 .. code-block:: Python def load_experimental_data(filepath): """Load experimental E_eff vs volume fraction data.""" data = np.loadtxt(filepath) return data[:, 0], data[:, 1] .. GENERATED FROM PYTHON SOURCE LINES 52-58 Forward Model ------------- For each set of candidate properties :math:`(E_f, \\nu_f)`, we compute the effective Young's modulus at all experimental volume fractions using simcoon's ``L_eff`` function. .. GENERATED FROM PYTHON SOURCE LINES 58-173 .. code-block:: Python def compute_E_eff( E_f, nu_f, concentrations, umat_name, param_list, path_keys, path_data, 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") param_list : list of Parameter Parameter objects with keys for file templating path_keys : str Path to template files path_data : str Path to working data directory (must be "data" relative to cwd) props_composite : numpy.ndarray Composite definition array [nphases, num_file, int1, int2, n_matrix] Returns ------- numpy.ndarray Effective Young's modulus at each volume fraction [MPa] """ E_eff = np.zeros(len(concentrations)) for i, c in enumerate(concentrations): param_list[0].value = 1.0 - c param_list[1].value = c param_list[2].value = E_f param_list[3].value = nu_f par.copy_parameters(param_list, path_keys, path_data) par.apply_parameters(param_list, path_data) L = sim.L_eff(umat_name, props_composite, 0, 0.0, 0.0, 0.0) 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, param_list, path_keys, path_data, 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, param_list, path_keys, path_data, props_composite, ) return np.mean((E_pred - E_exp) ** 2) except Exception: return 1e12 def identify_reinforcement( c_exp, E_exp, umat_name, param_list, path_keys, path_data, 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, param_list, path_keys, path_data, 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, } .. GENERATED FROM PYTHON SOURCE LINES 174-185 Main Execution -------------- We compare two identification strategies: 1. **Mori-Tanaka** (``MIMTN``): dilute approximation, best for c < 30% 2. **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 :math:`E_f`. Self-Consistent gives more physically meaningful identified properties. .. GENERATED FROM PYTHON SOURCE LINES 185-338 .. code-block:: Python if __name__ == "__main__": # ----------------------------------------------------------------- # Change to example directory (L_eff reads from data/ in cwd) # ----------------------------------------------------------------- try: script_dir = os.path.dirname(os.path.abspath(__file__)) except NameError: script_dir = os.getcwd() os.chdir(script_dir) path_keys = "keys_ident" path_data = "data" # ----------------------------------------------------------------- # Load experimental data (Wang 2003) # ----------------------------------------------------------------- c_exp, E_exp = load_experimental_data("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") # ----------------------------------------------------------------- # Load parameters (key system) # ----------------------------------------------------------------- # Define parameters inline (no external file needed): # @0p -> matrix concentration # @1p -> reinforcement concentration # @2p -> reinforcement E [MPa] # @3p -> reinforcement nu param_list = [ par.Parameter(0, bounds=(0.0, 1.0), key="@0p", sim_input_files=["Nellipsoids0.dat"]), par.Parameter(1, bounds=(0.0, 1.0), key="@1p", sim_input_files=["Nellipsoids0.dat"]), par.Parameter(2, bounds=(10000, 200000), key="@2p", sim_input_files=["Nellipsoids0.dat"]), par.Parameter(3, bounds=(0.01, 0.45), key="@3p", sim_input_files=["Nellipsoids0.dat"]), ] # Composite definition props_composite = np.array([2, 0, 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", param_list, path_keys, path_data, props_composite, bounds, ) # ----------------------------------------------------------------- # Identification with Self-Consistent # ----------------------------------------------------------------- result_SC = identify_reinforcement( c_exp, E_exp, "MISCN", param_list, path_keys, path_data, props_composite, bounds, ) # ----------------------------------------------------------------- # Summary # ----------------------------------------------------------------- print(f"\n{'=' * 60}") print(f" SUMMARY") print(f"{'=' * 60}") print(f" Reference glass: E_f ~ 73000 MPa, nu_f ~ 0.22") print(f"") 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", param_list, path_keys, path_data, props_composite, ) E_SC = compute_E_eff( result_SC["E_f"], result_SC["nu_f"], c_model, "MISCN", param_list, path_keys, path_data, props_composite, ) # Reference with handbook glass properties E_ref_MT = compute_E_eff( 73000, 0.22, c_model, "MIMTN", param_list, path_keys, path_data, props_composite, ) E_ref_SC = compute_E_eff( 73000, 0.22, c_model, "MISCN", param_list, path_keys, path_data, 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() .. image-sg:: /examples/heterogeneous/images/sphx_glr_composite_parameter_identification_001.png :alt: Composite Reinforcement Identification — Glass/Epoxy Mori-Tanaka vs Self-Consistent + Differential Evolution, Identified Reinforcement Properties, Handbook vs Identified Properties :srcset: /examples/heterogeneous/images/sphx_glr_composite_parameter_identification_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none ============================================================ 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.1728336983 differential_evolution step 10: f(x)= 2841691.1728336983 differential_evolution step 11: f(x)= 2841647.8003525925 differential_evolution step 12: f(x)= 2841602.746182008 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.1478061355 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.1084521133 differential_evolution step 25: f(x)= 2841396.1084521133 differential_evolution step 26: f(x)= 2841396.1084521133 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.396678291 differential_evolution step 31: f(x)= 2841391.396678291 differential_evolution step 32: f(x)= 2841391.396678291 differential_evolution step 33: f(x)= 2841391.396678291 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.0629338934 differential_evolution step 54: f(x)= 2841390.0629338934 differential_evolution step 55: f(x)= 2841390.0629338934 differential_evolution step 56: f(x)= 2841390.060185104 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.33263124695 differential_evolution step 2: f(x)= 52713.18502556835 differential_evolution step 3: f(x)= 52250.08051017501 differential_evolution step 4: f(x)= 52250.08051017501 differential_evolution step 5: f(x)= 51026.38861423012 differential_evolution step 6: f(x)= 50844.44748894559 differential_evolution step 7: f(x)= 50844.44748894559 differential_evolution step 8: f(x)= 50844.44748894559 differential_evolution step 9: f(x)= 50784.99401241427 differential_evolution step 10: f(x)= 50638.777388705355 differential_evolution step 11: f(x)= 50638.777388705355 differential_evolution step 12: f(x)= 50638.777388705355 differential_evolution step 13: f(x)= 50638.777388705355 differential_evolution step 14: f(x)= 50623.994479305074 differential_evolution step 15: f(x)= 50623.783592824555 differential_evolution step 16: f(x)= 50590.65814852624 differential_evolution step 17: f(x)= 50590.65814852624 differential_evolution step 18: f(x)= 50588.314183930845 differential_evolution step 19: f(x)= 50584.761535897684 differential_evolution step 20: f(x)= 50584.761535897684 differential_evolution step 21: f(x)= 50583.00679534065 differential_evolution step 22: f(x)= 50583.00679534065 differential_evolution step 23: f(x)= 50582.404035831416 differential_evolution step 24: f(x)= 50582.404035831416 differential_evolution step 25: f(x)= 50582.404035831416 differential_evolution step 26: f(x)= 50582.16590752159 differential_evolution step 27: f(x)= 50582.11894011469 differential_evolution step 28: f(x)= 50581.99120108114 differential_evolution step 29: f(x)= 50581.99120108114 differential_evolution step 30: f(x)= 50581.88523979694 differential_evolution step 31: f(x)= 50581.88523979694 differential_evolution step 32: f(x)= 50581.88523979694 differential_evolution step 33: f(x)= 50581.88523979694 differential_evolution step 34: f(x)= 50581.88523979694 differential_evolution step 35: f(x)= 50581.88523979694 differential_evolution step 36: f(x)= 50581.88523979694 differential_evolution step 37: f(x)= 50581.88523979694 differential_evolution step 38: f(x)= 50581.88147021726 differential_evolution step 39: f(x)= 50581.88147021726 differential_evolution step 40: f(x)= 50581.88147021726 differential_evolution step 41: f(x)= 50581.88121183913 differential_evolution step 42: f(x)= 50581.88104237749 differential_evolution step 43: f(x)= 50581.88085197378 differential_evolution step 44: f(x)= 50581.880180304695 differential_evolution step 45: f(x)= 50581.880180304695 differential_evolution step 46: f(x)= 50581.88017582471 differential_evolution step 47: f(x)= 50581.88017582471 differential_evolution step 48: f(x)= 50581.88017482083 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%). .. rst-class:: sphx-glr-timing **Total running time of the script:** (2 minutes 30.858 seconds) .. _sphx_glr_download_examples_heterogeneous_composite_parameter_identification.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: composite_parameter_identification.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: composite_parameter_identification.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: composite_parameter_identification.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_