ARD-RSC in simple shear

The result should match Figure 5 in Phelps, Jay H, Tucker, Charles L: "An anisotropic rotary diffusion model for fiber orientation in short- and long-fiber thermoplastics", Journal of Non-Newtonian Fluid Mechanics, 156, 3, 2009: DOI: 10.1016/j.jnnfm.2008.08.002.

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

# geometric factor
xi = 1.0

# time steps
t0, tf = 0, 4000
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]])

# compute solution by integrating the ODE
kwargs = {
    "xi": xi,
    "b1": 1.924e-4,
    "kappa": 1.0 / 30.0,
    "b2": 5.839e-3,
    "b3": 0.04,
    "b4": 1.168e-5,
    "b5": 0.0,
}

sol = solve_ivp(
    integrate_ori_ode,
    (t0, tf),
    A0.ravel(),
    t_eval=t,
    args=(L, IBOF_closure, ard_rsc_ode, kwargs),
)

# plot components
plt.plot(sol.t, sol.y[0], "-", label="A11")
plt.plot(sol.t, sol.y[4], "--", label="A22")
plt.plot(sol.t, sol.y[8], "-.", label="A33")
plt.plot(sol.t, sol.y[2], ":", label="A31")

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