iARD in simple shear

The result should match Figure 1 in Huan-Chang Tseng, Rong-Yeu Chang, Chia-Hsiang Hsu: "An objective tensor to predict anisotropic fiber orientation in concentrated suspensions", Journal of Rheology 60, 215 (2016): DOI: 10.1122/1.4939098

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, iardrpr_ode, integrate_ori_ode

# geometric factor
xi = 1.0

# time steps
t0, tf = 0, 2000
t = np.linspace(t0, tf, 100)

# 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]])

# compute solutions by integrating the ODEs
kwargs_ref = {
    "xi": xi,
    "b1": 3.842e-4,
    "kappa": 1.0 / 30.0,
    "b2": -1.786e-3,
    "b3": 5.25e-2,
    "b4": 1.168e-5,
    "b5": -5.0e-4,
}
sol_ref = solve_ivp(
    integrate_ori_ode,
    (t0, tf),
    A0.ravel(),
    t_eval=t,
    args=(L, IBOF_closure, ard_rsc_ode, kwargs_ref),
)

kwargs = {"xi": xi, "Ci": 0.025, "Cm": 1.0, "alpha": 0.967}
sol = solve_ivp(
    integrate_ori_ode,
    (t0, tf),
    A0.ravel(),
    t_eval=t,
    args=(L, IBOF_closure, iardrpr_ode, kwargs),
)

# plot components
plt.plot(sol_ref.t, sol_ref.y[0], "k-", label="A11 (ARD-RSC)")
plt.plot(sol_ref.t, sol_ref.y[4], "r-", label="A22 (ARD-RSC)")
plt.plot(sol_ref.t, sol_ref.y[8], "b-", label="A33 (ARD-RSC)")
plt.plot(sol_ref.t, sol_ref.y[2], "g-", label="A31 (ARD-RSC)")
plt.plot(sol.t, sol.y[0], "k.", label="A11 (iARD-RPR)")
plt.plot(sol.t, sol.y[4], "r.", label="A22 (iARD-RPR)")
plt.plot(sol.t, sol.y[8], "b.", label="A33 (iARD-RPR)")
plt.plot(sol.t, sol.y[2], "g.", label="A31 (iARD-RPR)")

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