Model comparison in simple shear

The result should match Figure 2 in Favaloro, A.J., Tucker III, C.L.: "Analysis of anisotropic diffusion models for fiber orientation", Composites Part A, 126 (2019): DOI: 10.1016/j.compositesa.2019.105605

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

from fiberoripy.closures import IBOF_closure
from fiberoripy.orientation import (
    ard_rsc_ode,
    folgar_tucker_ode,
    iard_ode,
    integrate_ori_ode,
    mrd_ode,
    pard_ode,
)

# geometric factor
xi = 1.0

# time steps
t0, tf = 0, 80
t = np.linspace(t0, tf, 500)

# initial fiber orientation state
A0 = 1.0 / 3.0 * np.eye(3)

# define a function that describes the (time-dependend) velocity gradient
def L(t):
    """Velocity gradient."""
    return np.array([[0.0, 0.0, 1.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])

kwargs_ft = {"xi": xi, "Ci": 0.0311}
sol_ft = solve_ivp(
    integrate_ori_ode,
    (t0, tf),
    A0.ravel(),
    t_eval=t,
    args=(L, IBOF_closure, folgar_tucker_ode, kwargs_ft),
)

kwargs_iard = {"xi": xi, "Ci": 0.0562, "Cm": 0.9977}
sol_iard = solve_ivp(
    integrate_ori_ode,
    (t0, tf),
    A0.ravel(),
    t_eval=t,
    args=(L, IBOF_closure, iard_ode, kwargs_iard),
)

kwargs_pard = {"xi": xi, "Ci": 0.0169, "Omega": 0.9868}
sol_pard = solve_ivp(
    integrate_ori_ode,
    (t0, tf),
    A0.ravel(),
    t_eval=t,
    args=(L, IBOF_closure, pard_ode, kwargs_pard),
)

kwargs_mrd = {"xi": xi, "Ci": 0.0198, "D1": 1.0, "D2": 0.4796, "D3": 0.0120}
sol_mrd = solve_ivp(
    integrate_ori_ode,
    (t0, tf),
    A0.ravel(),
    t_eval=t,
    args=(L, IBOF_closure, mrd_ode, kwargs_mrd),
)

kwargs_wpt = {
    "xi": xi,
    "b1": 0.0504 * (1.0 - 0.995),
    "kappa": 1.0,
    "b2": 0.0,
    "b3": 0.0504 * 0.995,
    "b4": 0.0,
    "b5": 0.0,
}
sol_wpt = solve_ivp(
    integrate_ori_ode,
    (t0, tf),
    A0.ravel(),
    t_eval=t,
    args=(L, IBOF_closure, ard_rsc_ode, kwargs_wpt),
)

kwargs_pt = {
    "xi": xi,
    "b1": 1.924e-4,
    "kappa": 1.0,
    "b2": 5.839e-3,
    "b3": 0.04,
    "b4": 1.168e-5,
    "b5": 0.0,
}
sol_pt = solve_ivp(
    integrate_ori_ode,
    (t0, tf),
    A0.ravel(),
    t_eval=t,
    args=(L, IBOF_closure, ard_rsc_ode, kwargs_pt),
)

# plot components
plt.plot(sol_ft.t, sol_ft.y[0], "k-", label="FT")
plt.plot(sol_iard.t, sol_iard.y[0], "b-.", label="iARD")
plt.plot(sol_pard.t, sol_pard.y[0], "b--", label="pARD")
plt.plot(sol_wpt.t, sol_wpt.y[0], "b:", label="WPT")
plt.plot(sol_mrd.t, sol_mrd.y[0], "r:", label="MRD")
plt.plot(sol_pt.t, sol_pt.y[0], "r--", label="PT")

# adjust some plot settings.
plt.xlabel("Time $t$ in s")
plt.ylim([0.3, 0.9])
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()