Orthotropic fitted closures

Model should reproduce Figure 6(a) in Du Hwan Chung and Tai Hun Kwon (2001), 'Improved model of orthotropic closure approximation for flow induced fiber orientation', Polymer Composites, 22(5), 636-649, DOI: 10.1002/pc.10566

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

from fiberoripy.closures import (
    IBOF_closure,
    hybrid_closure,
    orthotropic_fitted_closures,
)
from fiberoripy.orientation import folgar_tucker_ode, integrate_ori_ode

# geometric factor
xi = 1.0

# phenomenological fiber-fiber interaction coefficient
C_I = 0.01

# time steps
t = np.linspace(0, 30, 60)

# 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, 1.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])

# compute solutions by integrating the ODEs
kwargs_ft = {"xi": xi, "Ci": C_I}
sol_ibof = solve_ivp(
    integrate_ori_ode,
    (t.min(), t.max()),
    A0.ravel(),
    t_eval=t,
    args=(L, IBOF_closure, folgar_tucker_ode, kwargs_ft),
)
sol_ft = solve_ivp(
    integrate_ori_ode,
    (t.min(), t.max()),
    A0.ravel(),
    t_eval=t,
    args=(L, hybrid_closure, folgar_tucker_ode, kwargs_ft),
)
sol_orf = solve_ivp(
    integrate_ori_ode,
    (t.min(), t.max()),
    A0.ravel(),
    t_eval=t,
    args=(
        L,
        lambda a: orthotropic_fitted_closures(a, "ORF"),
        folgar_tucker_ode,
        kwargs_ft,
    ),
)
sol_orw = solve_ivp(
    integrate_ori_ode,
    (t.min(), t.max()),
    A0.ravel(),
    t_eval=t,
    args=(
        L,
        lambda a: orthotropic_fitted_closures(a, "ORW"),
        folgar_tucker_ode,
        kwargs_ft,
    ),
)
sol_orw3 = solve_ivp(
    integrate_ori_ode,
    (t.min(), t.max()),
    A0.ravel(),
    t_eval=t,
    args=(
        L,
        lambda a: orthotropic_fitted_closures(a, "ORW3"),
        folgar_tucker_ode,
        kwargs_ft,
    ),
)

# plot components
plt.plot(t, sol_ibof.y[0], linestyle="-", label="$a_{11}$ Hybrid", color="b")
plt.plot(t, sol_orf.y[0], linestyle="-", label="$a_{11}$ ORF", color="r")
plt.plot(t, sol_orw.y[0], linestyle="-", label="$a_{11}$ ORW", color="g")
plt.plot(t, sol_orw3.y[0], linestyle="-.", label="$a_{11}$ ORW3", color="k")

plt.plot(t, sol_ibof.y[1], linestyle="--", label="$a_{12}$ Hybrid", color="b")
plt.plot(t, sol_orf.y[1], linestyle="--", label="$a_{12}$ ORF", color="r")
plt.plot(t, sol_orw.y[1], linestyle="--", label="$a_{12}$ ORW", color="g")
plt.plot(t, sol_orw3.y[1], linestyle="-.", label="$a_{12}$ ORW3", color="k")

# adjust some plot settings.
plt.ylim = [-0.2, 1]
plt.grid(which="major", linestyle="-")
plt.minorticks_on()
plt.grid(which="minor", linestyle="--", alpha=0.2)
plt.xlabel("Time $t$ in $s$")
plt.ylabel("$a_{11}, a_{12}$")
plt.legend(loc="center right")
plt.title(r"Simple shear flow, $\xi = 1$, $C_I = 0.01$")
plt.show()