Rheo-SINDy: Finding a constitutive model from rheological data for complex fluids using sparse identification for nonlinear dynamics

Sato, Takeshi; Miyamoto, Souta; Kato, Shota

doi:10.1122/8.0000872

Rheology plays a pivotal role in understanding the flow behavior of fluids by discovering governing equations that relate deformation and stress, known as constitutive equations. Despite the importance of these equations, current methods for deriving them lack a systematic methodology, often relying on sense of physics and incurring substantial costs. To overcome this problem, we propose a novel method named Rheo-SINDy, which employs the sparse identification of nonlinear dynamics (SINDy) algorithm for discovering constitutive models from rheological data. Rheo-SINDy was applied to five distinct scenarios, four with well-established constitutive equations, and one without predefined equations. Our results demonstrate that Rheo-SINDy successfully identified accurate models for the known constitutive equations and derived physically plausible approximate models for the scenario without established equations. Notably, the identified approximate models can accurately reproduce nonlinear shear rheological properties, especially at steady state, including shear thinning. These findings validate the availability of Rheo-SINDy in handling data complexities and underscore its potential for advancing the development of data-driven approaches in rheology. Nevertheless, further refinement of learning strategies is essential for enhancing robustness to fully account for the complexities of real-world rheological data.

I. INTRODUCTION

Mathematical models grounded in physical laws offer profound insights into the behavior of complex systems across science and engineering. These models clarify the underlying mechanisms governing system dynamics and empower predictions and innovations in technology and natural science. Traditionally, model derivation has leaned heavily on theoretical and empirical knowledge, often requiring expert intuition. Data-driven methods have become capable of assisting in developing mathematical models and constructing models that provide advanced predictions [1]. These data-driven methods involve the sparse identification [2–5], symbolic regression [6–11], and physics-informed machine learning methods [12–15]. These methods have emerged as powerful tools for deriving governing equations directly from data.

Rheology is one of the scientific fields that address the properties of flowing matter, which plays a crucial role in many industries, such as designs of chemical processes, by providing insights into the flow behavior of complex fluids. One of the roles of rheology is to discover or derive governing equations that relate deformation and stress, referred to as constitutive equations [16]. Accurate constitutive equations are necessary to predict the flows of complex fluids under complex boundary conditions. Derivations of constitutive equations starting from the principles of continuum mechanics have achieved significant success in the field of rheology, yielding several practical constitutive equations [16]. Additionally, Ilg and Kröger proposed the derivation of constitutive equations following the guidelines of nonequilibrium thermodynamics [17]. Nevertheless, it is generally difficult to theoretically obtain constitutive equations for complex fluids from molecular models. Such cases usually explore mesoscopic coarse-grained models, which are based on molecular theories and suitable for numerical simulations. For example, for polymeric liquids, standard molecular theories have been proposed [18–20] and refined mesoscopic models have been constructed based on them [21–23]. In these models, the motion of individual (coarse-grained) molecules is numerically tracked. Although these models can reproduce rheological data with high accuracy, they require significantly more computational time compared to constitutive equations. Thus, a clear methodology is desired for obtaining constitutive equations from available data with the assistance of the rheological knowledge.

Data-driven methods have addressed the aforementioned challenges and advanced rheological studies such as constitutive modeling, flow predictions of complex fluids, and model selection [24,25]. Some applications have successfully identified constitutive relations of complex fluids or governing equations to predict the dynamics of fluids with knowledge of rheology. These studies have employed neural networks (NNs), including deep NN [26], graph NN [27], recurrent NN [28], physics-informed NN [29–31], multifidelity NN [32], and tensor basis NN [33]. Gaussian process regressions (GPRs) have also been employed, for example, for strain-rate dependent viscosity [34] or for viscoelastic properties [35–38].

Despite the success of NNs and GPRs, with a few exceptions such as the recent work of Lennon and co-workers [33], their black-box nature often obscures the underlying physics, making symbolic regression techniques more appealing for transparency and interoperability [1]. For example, this technique can successfully identify the governing equations for fluid flows [8]. Moreover, the sparse identification of nonlinear dynamics (SINDy) [2], which is one of such methods, has been utilized to track (reduced order) dynamics in the field of fluid mechanics [39]. Inspired by these successes, symbolic regression methods have recently started to be used in the field of rheology as well. For example, Mohammadamin and co-workers relied on SINDy for flexibly identifying the constitutive equations for an elasto-visco-plastic fluid [40]. Generally, the performance of SINDy is significantly influenced by the training data collection method, the candidate terms selected, and the optimization method. However, a comprehensive study to test SINDy for rheological data has not yet been conducted. Before applying SINDy to real-world rheological data, it is highly desirable to investigate fundamental learning strategies such as how to collect training data among various rheological tests and which optimization method to implement.

In this study, we employ SINDy to find constitutive models from rheological data, which we call as Rheo-SINDy. Figure 1 illustrates the schematic diagram of Rheo-SINDy. For the Rheo-SINDy regressions, we prepare a training data set including stress trajectories under simple and oscillatory shear flows and choose the candidate terms based on rheological knowledge of fundamental constitutive equations. Furthermore, multiple optimization methods are compared to find the effective ones for obtaining constitutive equations. This paper demonstrates five case studies. The first four cases verify the performance of Rheo-SINDy to identify the known constitutive equations, while the last one attempts to find the unknown constitutive equation for a coarse-grained mesoscopic model. Through the case studies, we validate the effectiveness of Rheo-SINDy and propose a strategy to find constitutive equations from rheological data. The details are shown below.

FIG. 1.

View large Download slide

Schematic illustration of Rheo-SINDy.

II. METHODS

We use a data-driven method known as a sparse identification of nonlinear dynamics (SINDy), which was originally developed by Brunton and co-workers [2], to obtain constitutive equations of complex fluid dynamics. The SINDy framework considers dynamical systems generally expressed by the following differential equation:

\frac{d x (t)}{d t} = \dot{x} (t) = f [x (t)],

(1)

where the vector

x (t)

represents the state of a system at time

t

and the function

f [x (t)]

determines the dynamics of the state

x (t)

⁠. The basic idea of SINDy is to find dominant terms for describing the dynamics out of numerous candidates using a sparse identification method. One can determine the (sparse) representation of

f

by a dataset including a collection of

x (t)

and

\dot{x} (t)

⁠. The regression to points of

x (t)

and

\dot{x} (t)

is computed with sparsity-promoting techniques, such as

ℓ_{1}

-regularization.

In the rheological community, it is of great importance to determine a relation between stress and strain rate. This relation is a so-called constitutive model or constitutive equation. Most constitutive equations are differential equations that depend on the (extra) stress tensor

τ

and velocity gradient tensor

κ

⁠. We prefer to use the so-called extra stress tensor

τ

as the stress tensor because this tensor satisfies

τ = 0

at equilibrium [16], which is convenient for the SINDy regression. The total stress tensor

σ

can be obtained by the relation

σ = τ + G I

⁠, where

G

is the modulus and

I

is the unit tensor. A general form for such constitutive equations can be written as

\frac{d τ (t)}{d t} = \dot{τ} (t) = f [τ (t), κ (t)] .

(2)

The velocity gradient tensor

κ (t)

is manipulated during rheological measurements. In formal constitutive equations, the time derivative should be frame-invariant, such as the upper-convected derivative [16]. As shown later in Sec. IV, the terms that appeared in the upper-convected derivative are recovered by the SINDy regressions.

We use the SINDy algorithm to find constitutive equations for complex fluids from data, and we refer to this technique as Rheo-SINDy (cf. Fig. 1). Rheo-SINDy requires three types of training data, two of which are transient stress data

T

and those time derivatives

\dot{T}

summarized as the following matrices:

T = [\begin{matrix} t_{x x} & t_{y y} & \dots & t_{z x} \end{matrix}] = [\begin{matrix} τ_{x x} (t_{1}) & τ_{y y} (t_{1}) & \dots & τ_{z x} (t_{1}) \\ τ_{x x} (t_{2}) & τ_{y y} (t_{2}) & \dots & τ_{z x} (t_{2}) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ τ_{x x} (t_{n}) & τ_{y y} (t_{n}) & \dots & τ_{z x} (t_{n}) \end{matrix}]

(3)

and

\dot{T} = [\begin{matrix} {\dot{t}}_{x x} & {\dot{t}}_{y y} & \dots & {\dot{t}}_{z x} \end{matrix}] = [\begin{matrix} {\dot{τ}}_{x x} (t_{1}) & {\dot{τ}}_{y y} (t_{1}) & \dots & {\dot{τ}}_{z x} (t_{1}) \\ {\dot{τ}}_{x x} (t_{2}) & {\dot{τ}}_{y y} (t_{2}) & \dots & {\dot{τ}}_{z x} (t_{2}) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ {\dot{τ}}_{x x} (t_{n}) & {\dot{τ}}_{y y} (t_{n}) & \dots & {\dot{τ}}_{z x} (t_{n}) \end{matrix}],

(4)

where

t_{μ ν}

(⁠

μ ν \in {x x, y y, z z, x y, y z, z x}

⁠) is the stress data for

n

sequential times. In this study, the stress data are collected by applying

κ (t)

to systems of prescribed constitutive equations or mesoscopic models of viscoelastic fluids. The time derivatives of the stress data

\dot{T}

are computed by a numerical differentiation method. When applying Rheo-SINDy with experimental data, we need to address the errors in time derivative data due to limited experimental time resolutions. To address this, it is effective to use a high-accuracy numerical differentiation scheme, such as the differentiation scheme with total variation regularization [41], as reported by Brunton and co-workers [2]. Nevertheless, in this study, we use a simple finite difference method and demonstrate that the correct equations can be regressed with this method. The remaining required data are the velocity gradient data

K

summarized as

K = [\begin{matrix} k_{x x} & k_{y y} & \dots & k_{z x} \end{matrix}] = [\begin{matrix} κ_{x x} (t_{1}) & κ_{y y} (t_{1}) & \dots & κ_{z x} (t_{1}) \\ κ_{x x} (t_{2}) & κ_{y y} (t_{2}) & \dots & κ_{z x} (t_{2}) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ κ_{x x} (t_{n}) & κ_{y y} (t_{n}) & \dots & κ_{z x} (t_{n}) \end{matrix}],

(5)

where

k_{μ ν}

(⁠

μ, ν \in {x, y, z}

⁠) are the velocity gradient data for

n

time steps.

In Rheo-SINDy, we construct a library matrix of functions

Θ

⁠, which includes various nonlinear functions, expressed as

Θ = [\begin{matrix} θ_{1} & θ_{2} & \dots & θ_{N_{Θ}} \end{matrix}] = [\begin{matrix} 1 & T & K & (T \otimes T) & (T \otimes K) & (K \otimes K) & \dots \end{matrix}],

(6)

where

N_{Θ}

is the total number of library functions and

T \otimes K

⁠, for example, denotes all possible combinations of the products of the row components in

T

and

K

for each time

t_{i}

(⁠

1 \leq i \leq n

⁠). The library

Θ

can incorporate not only polynomials but also other functions, such as sinusoidal functions. If the library does not contain necessary functions, the correct expression cannot be obtained by Rheo-SINDy. Thus, functions to be included in the library must be selected carefully by utilizing expertise in rheology and knowledge of the target data. The detailed procedure varies from case to case and is provided in Sec. IV.

Using these expressions, we can rewrite Eq. as

\dot{T} = Θ Ξ,

(7)

where

Ξ

is the sparse coefficient matrix written as

Ξ = [\begin{matrix} ξ_{x x} & ξ_{y y} & \dots & ξ_{z x} \end{matrix}] = [\begin{matrix} ξ_{x x, 1} & ξ_{y y, 1} & \dots & ξ_{z x, 1} \\ ξ_{x x, 2} & ξ_{y y, 2} & \dots & ξ_{z x, 2} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ ξ_{x x, N_{Θ}} & ξ_{y y, N_{Θ}} & \dots & ξ_{z x, N_{Θ}} \end{matrix}] .

(8)

To determine the coefficients

Ξ

⁠, we solve the following optimization problem for each row:

{\hat{ξ}}_{μ ν} = \underset{ξ_{μ ν}}{argmin} ‖ {\dot{t}}_{μ ν} - Θ ξ_{μ ν} ‖_{2}^{2} + R (ξ_{μ ν}),

(9)

where

{\hat{ξ}}_{μ ν}

is the optimized sparse vector,

‖ \dots ‖_{2}

is the

ℓ_{2}

-norm defined as

‖ x ‖_{2} = {(\sum_{i} x_{i}^{2})}^{1 / 2},

(10)

and

R (ξ_{μ ν})

is the regularization term. In SINDy, the appropriate optimization method generally depends on the specific problem, as has already been demonstrated by Fukami and co-workers [39]. To obtain a sparse solution of

Ξ

from rheological data, we apply the following five methods [39]: (i) the sequentially thresholded least square algorithm (STLSQ), (ii) sequentially thresholded Ridge regression (STRidge), (iii) least absolute shrinkage and selection operator (Lasso), (iv) Elastic-Net (E-Net), and (v) adaptive-Lasso (a-Lasso).

These methods employ different regularization terms as shown in Table I. The hyperparameters of

ℓ_{i}

norm (⁠

i = 0, 1, 2

⁠) are denoted as

λ_{i}

(⁠

> 0

⁠). The

ℓ_{0}

and

ℓ_{1}

norms are defined as

‖ ξ_{μ ν} ‖_{0} = \sum_{j} δ (ξ_{μ ν, j})

(11)

and

‖ ξ_{μ ν} ‖_{1} = \sum_{j} | ξ_{μ ν, j} |,

(12)

where

δ (ξ_{μ ν, j})

is the Kronecker delta function, which is equal to

1

if

ξ_{μ ν, j} \neq 0

and

0

otherwise. The vector

ξ_{μ ν}^{'}

in the a-Lasso is defined as

ξ_{μ ν}^{'} = w_{μ ν} \otimes ξ_{μ ν}

⁠, where

\otimes

is the element-wise product and

w_{μ ν}

is the adaptive weight vector and its

j

th element is defined as

w_{μ ν, j} = | ξ_{μ ν, j} |^{- δ}

with

δ

being the positive constant.

TABLE I.

The regularization term $R (ξ_{μ ν})$ for the sparse regression methods.

Method	Regularization term $R (ξ_{μ ν})$
STLSQ	$λ_{0} ‖ ξ_{μ ν} ‖_{0}$
STRidge	$λ_{0} ‖ ξ_{μ ν} ‖_{0} + λ_{2} ‖ ξ_{μ ν} ‖_{2}^{2}$
Lasso	$λ_{1} ‖ ξ_{μ ν} ‖_{1}$
E-Net	$λ_{1} ‖ ξ_{μ ν} ‖_{1} + λ_{2} ‖ ξ_{μ ν} ‖_{2}^{2}$
a-Lasso	$λ_{1} ‖ ξ_{μ ν}^{'} ‖_{1}$

The STLSQ and STRidge were implemented by iteratively conducting the least square regression and the Ridge regression, respectively, while setting the coefficients with smaller absolute values than a certain threshold

α

(⁠

> 0

⁠) to zero based on the original papers [2,42]. Since this coefficient selection by

α

replaces the role of the

ℓ_{0}

norm,

λ_{0}

is set to zero in this implementation. In the STRidge, the hyperparameter

λ_{2}

was set to

0.05

⁠. The Lasso, E-Net, and a-Lasso were implemented using the scikit-learn library [43]. In this library, the loss functions

L (ξ_{μ ν})

for the Lasso and E-Net are, respectively, defined as

L (ξ_{μ ν}) = \frac{1}{2 n} ‖ {\dot{t}}_{μ ν} - Θ ξ_{μ ν} ‖_{2}^{2} + α ‖ ξ_{μ ν} ‖_{1}

(13)

and

\begin{aligned} L (ξ_{μ ν}) = \frac{1}{2 n} ‖ {\dot{t}}_{μ ν} - Θ ξ_{μ ν} ‖_{2}^{2} + α β ‖ ξ_{μ ν} ‖_{1} + \frac{α (1 - β)}{2} ‖ ξ_{μ ν} ‖_{2}^{2}, \end{aligned}

(14)

where

β

is the

ℓ_{1}

ratio and

α

and

β

are the hyperparameters. These two loss functions have the same form when

β = 1

⁠. In this study,

β

for the E-Net was set to

0.5

⁠. According to the original paper of the a-Lasso [44], it can be implemented as the Lasso problem as the following steps:

Define $Θ^{'} = [θ_{1}^{'}, \dots, θ_{N_{Θ}}^{'}]$ ⁠, where $θ_{j}^{'} = θ_{j} / w_{μ ν, j}$ (⁠ $j = 1, \dots, N_{Θ}$ ⁠).
Solve the Lasso problem to obtain ${\hat{ξ}}_{μ ν}^{'}$ using Eq. (13) with $Θ^{'}$ ⁠.
Output ${\hat{ξ}}_{μ ν, j} = {\hat{ξ}}_{μ ν, j}^{'} / w_{μ ν, j}$ (⁠ $j = 1, \dots, N_{Θ}$ ⁠).

The adaptive weight

w_{μ ν, j}

depends on the coefficients, and thereby the output coefficients can be varied in each iteration. To obtain the converged solution, we initialized the weights as unit vectors

w = 1

and repeated the above steps until the coefficients

{\hat{ξ}}_{μ ν, j}

no longer change [39]. Here, the hyperparameter

δ

was set to

3

(see Sec. S1 in the supplementary material for the effect of

δ

⁠). As shown above, each method has a hyperparameter

α

to penalize the solution complexity, which is to be tuned for obtaining good predictive yet parsimonious representations. For this purpose, we test various

α

values and pick the one with the smallest number of terms among the results whose error has the same order as the minimum error (when

α

is sufficiently small).

In this study, we focus on shear rheological measurements that give fundamental rheological properties because they are well-studied and suitable for discussing the applicability of our method to rheology data. Under shear flow, among the components of $κ$ ⁠, only $κ_{x y}$ has nonzero values. Here, $x$ is the velocity direction, and $y$ is the velocity gradient direction. Since the major stress components are $τ_{x x}$ ⁠, $τ_{y y}$ ⁠, $τ_{z z}$ ⁠, and $τ_{x y}$ under shear flow, we only use these components to conduct Rheo-SINDy.

III. CASE STUDIES

For case studies, we first test whether Rheo-SINDy can find appropriate constitutive equations from training data generated by phenomenological constitutive equations, namely, the upper-convected Maxwell (UCM) model and the Giesekus model. Subsequently, we apply Rheo-SINDy to data generated by several dumbbell-based models. This section provides a brief overview of the UCM, Giesekus, and dumbbell models used in this study and the conditions for creating the training datasets. Table II summarizes the conditions for generating training datasets.

TABLE II.

Conditions for generating training datasets.

Model	Parameter	Shear typea	Shear parameters	Δt	Δt_train	Simulation time
UCM	—	S	$\dot{γ} \in {1, 1.7, \dots, 100}$ b	10⁻⁴	10⁻²	10
UCM	—	O	γ₀ = 2, ω = 1	10⁻⁴	10⁻²	100
Giesekus	α_G = 0.5	O	γ₀ = 2, ω ∈ {0.1, 0.2, …, 1}	10⁻⁴	10⁻²	100
Hookean dumbbellc	n_K = 10	O	γ₀ = 2, ω = 0.5	10⁻³	10⁻²	100
FENE-P dumbbelld	n_K = 10	O	γ₀ = 2 or 8, ω ∈ {0.1, 0.2, …, 1}	10⁻⁴	10⁻²	100
FENE dumbbelle	n_K = 10	O	γ₀ = 2 or 8, ω ∈ {0.1, 0.2, …, 1}	10⁻⁴	10⁻²	100

Model	Parameter	Shear typea	Shear parameters	Δt	Δt_train	Simulation time
UCM	—	S	$\dot{γ} \in {1, 1.7, \dots, 100}$ b	10⁻⁴	10⁻²	10
UCM	—	O	γ₀ = 2, ω = 1	10⁻⁴	10⁻²	100
Giesekus	α_G = 0.5	O	γ₀ = 2, ω ∈ {0.1, 0.2, …, 1}	10⁻⁴	10⁻²	100
Hookean dumbbellc	n_K = 10	O	γ₀ = 2, ω = 0.5	10⁻³	10⁻²	100
FENE-P dumbbelld	n_K = 10	O	γ₀ = 2 or 8, ω ∈ {0.1, 0.2, …, 1}	10⁻⁴	10⁻²	100
FENE dumbbelle	n_K = 10	O	γ₀ = 2 or 8, ω ∈ {0.1, 0.2, …, 1}	10⁻⁴	10⁻²	100

^a

“S” means the steady shear flow and “O” means the oscillatory shear flow.

^b

Ten equally spaced values on a logarithmic scale between $\dot{γ} = 1$ and $\dot{γ} = 100$ were used.

^c

Three different numbers of dumbbells (N_p ∈ {10³, 10⁴, 10⁵}) were addressed.

^d

The closed-form expression (N_p → ∞) was examined.

^e

The number of dumbbells was set to N_p = 10⁴.

A. Upper-convected Maxwell (UCM) model

The simplest constitutive equation for viscoelastic fluids is the UCM model [16] shown as

\frac{d τ}{d t} - τ \cdot κ^{T} - κ \cdot τ = - \frac{1}{λ} τ + 2 G D .

(15)

Here, the left-hand side of Eq. is the upper-convected time derivative of

τ

⁠,

κ^{T}

is the transposed

κ

⁠,

λ

is the relaxation time,

G

is the modulus, and

D

is the deformation rate tensor defined as

D = (κ + κ^{T}) / 2

⁠. Using

λ

as the unit time and

G

as the unit stress (i.e.,

λ = G = 1

⁠), we can obtain dimensionless expressions for time

\tilde{t} = t / λ

⁠, velocity gradient tensor

\tilde{κ} = λ κ

⁠, and stress

\tilde{τ} = τ / G

⁠. In what follows, we omit the tilde in dimensionless variables for simplicity. The dimensionless form of the UCM model under shear flow is thus written as

{\dot{τ}}_{x x} = - τ_{x x} + 2 τ_{x y} κ_{x y},

(16)

{\dot{τ}}_{y y / z z} = - τ_{y y / z z} = 0,

(17)

{\dot{τ}}_{x y} = - τ_{x y} + κ_{x y} + τ_{y y} κ_{x y} = - τ_{x y} + κ_{x y} .

(18)

Since the initial conditions for

τ

are set to the values of

τ

at equilibrium, namely,

τ = 0

⁠,

τ_{y y / z z}

of the UCM model is zero under shear flow. Here, “yy/zz” in Eq. (17) means “yy or zz.”

For the UCM model, we generate training data by numerically solving Eqs. (16)–(18) under two shear flow scenarios: steady shear and oscillatory shear tests. For the steady shear test, the shear rate is kept constant (⁠ $κ_{x y} = \dot{γ}$ ⁠) across various values (⁠ $\dot{γ} \in {1, 1.7, 2.8, 4.6, 7.7, 13, 22, 36, 60, 100}$ ⁠) with simulations running from $t = 0$ to $t = 10$ using a time step of $Δ t = 1.0 \times 10^{- 4}$ ⁠. The oscillatory shear test introduces a time-dependent oscillatory shear strain, $γ (t) = γ_{0} \sin (ω t)$ [i.e., $κ_{x y} (t) = \dot{γ} (t) = γ_{0} ω \cos (ω t)$ ], with $γ_{0} = 2$ and $ω = 1$ ⁠, over a period from $t = 0$ to $t = 100$ ⁠, employing the same time step. In both tests, data are collected at intervals of $Δ t_{train} = 1 \times 10^{- 2}$ ⁠, resulting in a total of $10^{4}$ data points for the training data.

B. Giesekus model

The Giesekus model, which is one of the most popular phenomenological constitutive equations [45], shows typical shear rheological properties and is used to fit various complex fluids, including polymer solutions and wormlike micellar solutions. The tensorial form of the Giesekus constitutive equation can be written as

\frac{d τ}{d t} - τ \cdot κ^{T} - κ \cdot τ = - \frac{1}{λ} τ - \frac{α_{G}}{G λ} τ \cdot τ + 2 G D,

(19)

where

α_{G}

is the parameter governing the nonlinear response of the Giesekus model. The Giesekus equation under shear flow is, thus, given by

{\dot{τ}}_{x x} = - τ_{x x} - α_{G} (τ_{x x}^{2} + τ_{x y}^{2}) + 2 τ_{x y} κ_{x y},

(20)

{\dot{τ}}_{y y} = - τ_{y y} - α_{G} (τ_{y y}^{2} + τ_{x y}^{2}),

(21)

{\dot{τ}}_{z z} = 0,

(22)

{\dot{τ}}_{x y} = - τ_{x y} - α_{G} (τ_{x x} + τ_{y y}) τ_{x y} + τ_{y y} κ_{x y} + κ_{x y} .

(23)

Here, all quantities are nondimensionalized by using

λ

as the unit time and

G

as the unit stress. From Eqs. , the total number of collect terms in the Giesekus model is

12

⁠.

We generate the training data by solving Eqs. (20)–(23) numerically with $α_{G} = 0.5$ and $Δ t = 1 \times 10^{- 4}$ ⁠. We note that the Giesekus model with $α_{G} = 0.5$ gives sufficient nonlinear features under shear flow. We applied the oscillatory shear flow with $γ_{0} = 2$ and various $ω$ values (⁠ $ω \in {0.1, 0.2, \dots, 1}$ ⁠) for $0 \leq t \leq 100$ ⁠. From the computed stress data, we collected data at the interval of $Δ t_{train} = 1 \times 10^{- 2}$ ⁠.

C. Dumbbell models

The dumbbell-based models have been widely utilized in numerous previous studies for the computation of viscoelastic fluids and are regarded as the standard mesoscopic model for viscoelastic fluids [18]. As a stepping stone to investigating more complex models, we chose this model, which can be considered as a fundamental model in rheology. As illustrated in Fig. 2, a dumbbell consists of two beads (indexed as

1

or

2

⁠) and a spring that connects them. The Langevin equations for the positions of the two beads

r_{1 / 2} (t)

can be written as

ζ [\frac{d r_{i} (t)}{d t} - κ \cdot r_{i} (t)] = - h (t) {r_{i} (t) - r_{j} (t)} + F_{B i} (t),

(24)

with

(i, j) = (1, 2)

or

(2, 1)

⁠. Here,

ζ

is the friction coefficient,

h (t)

is the spring strength, and

F_{B i} (t)

is the Brownian force acting on the bead

i

⁠. The time evolution equation for the end-to-end vector

R (t)

[

= r_{2} (t) - r_{1} (t)

] of the beads is, thus, obtained as

\begin{aligned} ζ [\frac{d R (t)}{d t} - κ \cdot R (t)] = - 2 h (t) R (t) + {F_{B 2} (t) - F_{B 1} (t)} . \end{aligned}

(25)

The Brownian force is characterized by the first and second-moment averages as

⟨ F_{B i} (t) ⟩ = 0

(26)

and

⟨ F_{B i} (t) F_{B j} (t^{'}) ⟩ = 2 ζ k_{B} T δ_{i j} δ (t - t^{'}) I,

(27)

where

k_{B}

is the Boltzmann constant and

T

is the temperature. From the end-to-end vector

R (t)

⁠, the stress tensor can be expressed as

τ (t) = ρ ⟨ h (t) R (t) R (t) ⟩ - ρ k_{B} T I,

(28)

where

ρ

is the density of dumbbells.

FIG. 2.

View large Download slide

Schematic illustration of the dumbbell model.

There are several expressions for the spring strength

h (t)

⁠. The most basic one is the Hookean spring, defined as

h (t) = h_{eq} = \frac{3 k_{B} T}{n_{K} b_{K}^{2}},

(29)

where

n_{K}

is the number of Kuhn segments per spring and

b_{K}

is the Kuhn length. Reproduction of some nonlinear rheological properties, such as shear thinning under shear flow, necessitates dealing with finite extensible nonlinear elastic (FENE) effects. Although the exact expression for FENE springs is given by the inverse Langevin function, the following empirical expression is widely used [18]:

h (t) = h_{eq} \frac{1 - ⟨ R_{eq}^{2} ⟩ / R_{\max}^{2}}{1 - R^{2} (t) / R_{\max}^{2}},

(30)

where

⟨ R_{eq}^{2} ⟩^{1 / 2} = {(n_{K})}^{1 / 2} b_{K}

is the equilibrium length of the springs and

R_{\max} = n_{K} b_{K}

is the maximum length of the springs. As shown later in Sec. III C 3, a constitutive equation cannot be analytically obtained for the FENE dumbbell model. To address the FENE spring more analytically, the following approximate expression of the FENE spring has been proposed [18]:

h (t) = h_{eq} \frac{1 - ⟨ R_{eq}^{2} ⟩ / R_{\max}^{2}}{1 - ⟨ R^{2} (t) ⟩ / R_{\max}^{2}} = h_{eq} f_{FENE} (t) .

(31)

This spring is referred to as the FENE-P spring. Here, “P” means Peterlin, who proposed the approximate form of the FENE spring law. The average appearing in Eq. makes it possible to obtain the analytical constitutive equation.

We use $λ = ζ / 4 h_{eq}$ as the unit time and $G = ρ k_{B} T$ as the unit stress for the dumbbell models. To simplify the expressions, we omit the tilde representing dimensionless quantities in what follows.

1. Hookean dumbbell model

The most basic dumbbell model is the Hookean dumbbell model, where Hookean springs are employed [cf. Eq. (29)]. From Eqs. (25), (28), and (29), the Hookean dumbbell model reduces to the constitutive equation for the UCM model [cf. Eq. (15)] in the limit of $N_{p} \to \infty$ with $N_{p}$ being the number of dumbbells.

For the Hookean dumbbell model, we generate training data by Brownian dynamics (BD) simulations with the finite numbers of dumbbells (⁠ $N_{p} \in {10^{3}, 10^{4}, 10^{5}}$ ⁠). In the BD simulations of this study, to integrate Eq. (25), we use the explicit Euler method with a small time step. We apply the oscillatory shear flow, $γ (t) = γ_{0} \sin (ω t)$ ⁠, with $γ_{0} = 2$ and $ω = 0.5$ ⁠, over a period from $t = 0$ to $t = 100$ ⁠. The simulations are run with $Δ t = 1 \times 10^{- 3}$ for $0 \leq t \leq 100$ and data are collected at the interval of $Δ t_{train} = 1 \times 10^{- 2}$ ⁠. Each simulation is conducted with five different random seeds, and their average data are used for training. Due to the characteristics of the BD simulation, the training data inherently include noise originating from the finite $N_{p}$ ⁠. We here test whether Rheo-SINDy can find from the noisy data, the constitutive equations for the UCM model shown in Eqs. (16)–(18).

2. FENE-P dumbbell model

We next address the so-called FENE-P dumbbell model, where Eq. (31) is utilized as the spring strength. As shown below, the FENE-P dumbbell model has an analytical solution and is utilized for various flow problems [46,47].

Due to the assumption shown in Eq. , a simple representation of the time evolution for the conformation tensor

C = ⟨ R (t) R (t) ⟩

can be obtained as

\frac{d C}{d t} - C \cdot κ^{T} - κ \cdot C = - f_{FENE} (t) C + \frac{n_{K}}{3} I .

(32)

The stress tensor is, thus, obtained by

τ (t) = ρ h (t) C (t) - ρ k_{B} T I .

(33)

Under shear flow, Eq. reduces to the following expressions:

{\dot{C}}_{x x} = - f_{FENE} C_{x x} + 2 C_{x y} κ_{x y} + \frac{n_{K}}{3},

(34)

{\dot{C}}_{y y / z z} = - f_{FENE} C_{y y / z z} + \frac{n_{K}}{3},

(35)

{\dot{C}}_{x y} = - f_{FENE} C_{x y} + C_{y y} κ_{x y} .

(36)

Using Rheo-SINDy, we test whether or not Eqs. can be discovered from the data.

Although it has not been as widely recognized due to its complexity, the FENE-P dumbbell model can also be expressed in the form of the constitutive equation (i.e., the stress expression) [48]. From the textbook of Bird and co-workers [18], the constitutive equation for the FENE-P model is

\begin{aligned} \frac{d τ}{d t} - τ \cdot κ^{T} - κ \cdot τ = - f_{FENE} (t) τ + 2 D + \frac{D \ln Z}{D t} (τ + I), \end{aligned}

(37)

where

D (\dots) / D t

is the substantial derivative and

Z

is the function expressed as

Z = \frac{1}{1 - ⟨ R^{2} (t) / R_{\max}^{2} ⟩} = 1 + \frac{1}{3 n_{K} Z_{eq}^{- 1}} (tr τ + 3) .

(38)

Here,

Z_{eq}

indicates

Z

at equilibrium. From Eq. , we can see that

tr τ

is tightly related to the (squared) length of dumbbells. Since we do not address the spatial gradient in rheological calculations,

D (\dots) / D t

simply reduces to

d (\dots) / d t

⁠. Using Eqs. , , and , the constitutive equations for the FENE-P dumbbell model under shear flow can be expressed as

\begin{aligned} {\dot{τ}}_{x x} & = - {1 + \frac{1}{3 (n_{K} - 1)}} τ_{x x} - \frac{1}{3 (n_{K} - 1)} (τ_{y y} + τ_{z z}) \\ - \frac{1}{9 n_{K} (n_{K} - 1)} {(tr τ)}^{2} - \frac{1}{3 n_{K}} (2 + \frac{1}{n_{K} - 1}) tr τ τ_{x x} \\ + 2 {1 + \frac{1}{3 (n_{K} - 1)}} τ_{x y} κ_{x y} - \frac{1}{9 n_{K} (n_{K} - 1)} {(tr τ)}^{2} τ_{x x} \\ + \frac{2}{3 (n_{K} - 1)} τ_{x x} τ_{x y} κ_{x y}, \end{aligned}

(39)

\begin{aligned} {\dot{τ}}_{y y / z z} & = - {1 + \frac{1}{3 (n_{K} - 1)}} τ_{y y / z z} - \frac{1}{3 (n_{K} - 1)} (τ_{x x} + τ_{z z / y y}) \\ - \frac{1}{9 n_{K} (n_{K} - 1)} {(tr τ)}^{2} - \frac{1}{3 n_{K}} (2 + \frac{1}{n_{K} - 1}) tr τ τ_{y y / z z} \\ + \frac{2}{3 (n_{K} - 1)} τ_{x y} κ_{x y} - \frac{1}{9 n_{K} (n_{K} - 1)} {(tr τ)}^{2} τ_{y y / z z} \\ + \frac{2}{3 (n_{K} - 1)} τ_{y y / z z} τ_{x y} κ_{x y}, \end{aligned}

(40)

\begin{aligned} {\dot{τ}}_{x y} & = - τ_{x y} + κ_{x y} + τ_{y y} κ_{x y} - \frac{1}{3 n_{K}} (2 + \frac{1}{n_{K} - 1}) tr τ τ_{x y} \\ - \frac{1}{9 n_{K} (n_{K} - 1)} {(tr τ)}^{2} τ_{x y} + \frac{2}{3 (n_{K} - 1)} τ_{x y}^{2} κ_{x y} . \end{aligned}

(41)

For the derivation, please refer to Sec. S2 in the supplementary material. From Eqs. , the constitutive equation for the FENE-P model can be expressed by a polynomial of up to a third degree in

τ

and

κ

⁠. Here, we note that Eqs. become equivalent to the UCM model shown in Eqs. in the limit of

n_{K} \to \infty

⁠.

To generate noise-free training data, we solve Eqs. (33)–(36) with $n_{K} = 10$ and $Δ t = 1 \times 10^{- 4}$ for $0 \leq t \leq 100$ ⁠. We apply the oscillatory shear flow with $γ_{0} = 2$ or $8$ and various $ω$ values (⁠ $ω \in {0.1, 0.2, \dots, 1}$ ⁠). From the computed stress data, we collect data at the interval of $Δ t_{train} = 1 \times 10^{- 2}$ ⁠.

3. FENE dumbbell model

We finally address the FENE dumbbell model, where the spring strength is represented by Eq. (30). Since the FENE dumbbell model does not use any simplification for the spring strength [e.g., Peterlin approximation shown in Eq. (31)], its analytical constitutive equation has not been obtained. We apply Rheo-SINDy to this case to see if an approximate constitutive equation can be obtained. The obtained equations are validated by comparing the data obtained by numerically solving them with the data obtained by BD simulations.

The training data are generated by the BD simulations using Eqs. (25)–(28) and (30) with $n_{K} = 10$ ⁠, $N_{p} = 10^{4}$ ⁠, and $Δ t = 1 \times 10^{- 4}$ for $0 \leq t \leq 100$ ⁠. We apply the oscillatory shear flows with the same parameters as those in the FENE-P dumbbell model. The BD simulation results with five different random seeds are averaged for each condition. Since we do not use any approximation for the spring strength, the values of $h (t)$ differ for each individual dumbbell. From the computed stress data, we collected data at the interval of $Δ t_{train} = 1 \times 10^{- 2}$ ⁠.

IV. RESULTS AND DISCUSSIONS

In this section, we present the results of the case studies for five models explained in Sec. III. Through case studies on phenomenological constitutive equations (i.e., the UCM and Giesekus models), we first investigate two key questions for deriving accurate equations: (i) whether oscillatory shear tests or simple (steady) shear tests are more appropriate and (ii) which of the five optimization methods presented in Sec. II is most suitable. Subsequently, using the experimental and optimization methods identified as appropriate by this investigation, we attempt to obtain the constitutive equations of the dumbbell models. The methods of identifying constitutive models are summarized in Table III. In Secs. IV A–IV F, we examine how the hyperparameter $α$ affects the Rheo-SINDy results. To discuss the effect of the hyperparameter $α$ ⁠, we also examine the results with various $α$ (⁠ $10^{- 9} \leq α \leq 10^{3}$ ⁠). Among the $α$ values examined, we pick $α$ with the smallest number of terms among the results whose error has the same order as the minimum error. For comparison, Rheo-SINDy results with $α$ other than the one chosen by this criterion are also shown.

TABLE III.

Settings of Rheo-SINDy used in case studies.

Model	Exact equations	Optimizationa	Library	$N_{Θ}$
UCM	Eqs. (16)–(18)	All five methodsb	Polynomial (up to third order)	35
Giesekus	Eqs. (20)–(23)	STLSQ, STRidge, a-Lasso	Polynomial (up to second order)	15
Hookean dumbbell	Eqs. (16)–(18)	STRidge, a-Lasso	Polynomial (up to second order)	15
FENE-P dumbbell	Eqs. (34)–(36)	STRidge, a-Lasso	Eq. (42)	26
FENE-P dumbbell	Eqs. (39)–(41)	STRidge, a-Lasso	Eq. (43)	29
FENE dumbbell	N/A	STRidgec, a-Lasso	Eq. (43)	29

Model	Exact equations	Optimizationa	Library	$N_{Θ}$
UCM	Eqs. (16)–(18)	All five methodsb	Polynomial (up to third order)	35
Giesekus	Eqs. (20)–(23)	STLSQ, STRidge, a-Lasso	Polynomial (up to second order)	15
Hookean dumbbell	Eqs. (16)–(18)	STRidge, a-Lasso	Polynomial (up to second order)	15
FENE-P dumbbell	Eqs. (34)–(36)	STRidge, a-Lasso	Eq. (42)	26
FENE-P dumbbell	Eqs. (39)–(41)	STRidge, a-Lasso	Eq. (43)	29
FENE dumbbell	N/A	STRidgec, a-Lasso	Eq. (43)	29

^a

The value of α is picked from sparser solutions whose error is the same order of magnitude as the minimum error.

^b

All five methods are the STLSQ, STRidge, Lasso, E-Net, and a-Lasso explained in Sec. II.

^c

The results are shown in the supplementary material.

A. Upper convected Maxwell model

Through this case study, we first check the appropriate methods to take the shear rheological data for Rheo-SINDy. Figure 3 shows the training data and results for the UCM model. Figures 3(a) and 3(b) are the stress data under the steady shear flows with the various shear rates and those under the oscillatory shear flow.

We conducted the Rheo-SINDy regressions by using the polynomial library that includes up to third order terms of $τ_{x x}$ ⁠, $τ_{y y}$ ⁠, $τ_{x y}$ ⁠, and $κ_{x y}$ ⁠. Thus, there were $35$ candidate terms for each component of the constitutive equation. The terms related to $τ_{z z}$ were excluded because they do not contribute to the UCM dynamics. Figures 3(c) and 3(d) present the numbers of total terms varying with the hyperparameter $α$ obtained by Rheo-SINDy using the training data (a) and (b), respectively. Figure 3(c) indicates that sparse solutions can be obtained by the STLSQ, STRidge, and a-Lasso, but not by the Lasso and E-Net. Moreover, regarding the number of terms, the STLSQ and STRdge exhibit similar behavior. Specifically, the STLSQ and STRidge with $3 \times 10^{- 3} \leq α \leq 3 \times 10^{- 1}$ successfully discovered equations including the correct number of terms. Figure 3(d) indicates that the STLSQ, STRidge, and a-Lasso yielded the correct number of terms, though all five methods gave sparse solutions. In most of the cases where the number of terms obtained was correct, the obtained coefficients were also correct. These results suggest that the oscillatory shear test is more appropriate than the steady shear test to obtain the correct constitutive equations for the UCM model. Figure 3(e) lists the constitutive equations obtained by the STRidge and a-Lasso. We can see that the STRidge and a-Lasso can give the correct constitutive equations, except for the a-Lasso in the steady shear test. Furthermore, we confirmed that the correct equations were obtained even for $α$ values not shown in Fig. 3(e) in the case of the UCM model. These findings show the basic validity of finding the constitutive equations from the rheological data by Rheo-SINDy. Figure 3 indicates that the STLSQ, STRidge, and a-Lasso demonstrate better performance in discovering the correct constitutive equations compared to the Lasso and E-Net; thus, we use the former three methods in the following discussion.

FIG. 3.

Training data obtained by the UCM model [cf. Eqs. (16)–(18)] (a) under steady shear flow ( κ x y = γ ˙) and (b) under oscillatory shear flow [ κ x y = γ 0 ω cos ( ω t )]. The numbers of total terms obtained (c) from the training data (a) and (d) from the training data (b). (e) The constitutive equations obtained by Rheo-SINDy. In (b), x x-, y y-, and x y-components of the stress tensor are plotted with the black solid, red dotted, and blue dashed-dotted lines, respectively. In (c) and (d), the numbers of total terms for five different optimization methods are plotted against the hyperparameter α. The black horizontal lines in (c) and (d) indicate the correct number of the terms in the UCM model. The color version is available online.

View large Download slide

Training data obtained by the UCM model [cf. Eqs. (16)–(18)] (a) under steady shear flow (⁠ $κ_{x y} = \dot{γ}$ ⁠) and (b) under oscillatory shear flow [ $κ_{x y} = γ_{0} ω \cos (ω t)$ ]. The numbers of total terms obtained (c) from the training data (a) and (d) from the training data (b). (e) The constitutive equations obtained by Rheo-SINDy. In (b), $x x$ -, $y y$ -, and $x y$ -components of the stress tensor are plotted with the black solid, red dotted, and blue dashed-dotted lines, respectively. In (c) and (d), the numbers of total terms for five different optimization methods are plotted against the hyperparameter $α$ ⁠. The black horizontal lines in (c) and (d) indicate the correct number of the terms in the UCM model. The color version is available online.

B. Giesekus model

We here explain the results of Rheo-SINDy for the Giesekus model. This case used the polynomial library consisting of up to second-order terms of $τ_{x x}$ ⁠, $τ_{y y}$ ⁠, $τ_{x y}$ ⁠, and $κ_{x y}$ ⁠, which contain sufficient candidate terms to obtain the exact equations. Figure 4 shows (a) the total number of terms and (b) the error rate obtained by Rheo-SINDy for the training data of the Giesekus model. The error rate is defined as the sum of the mean squared errors (MSEs) of ${\dot{t}}_{μ ν} - Θ {\hat{ξ}}_{μ ν}$ ⁠. The MSEs were scaled so that the maximum value of each method was 1. Figures 4(a) and 4(b) indicate that the a-Lasso evidently provides a sparser solution compared to the other two methods when the error rates are comparable. We note that, similar to the number of terms, coefficient values generally depend on $α$ ⁠.

FIG. 4.

View large Download slide

(a) The number of total terms and (b) the error rate of the constitutive equations obtained by Rheo-SINDy for the Giesekus model. The used optimization methods are the STLSQ (green squares), STRidge (red reverse triangles), and a-Lasso (blue triangles). The regressions were conducted with the multiple (⁠ $10$ ⁠) data trajectories of $κ_{x y} = γ_{0} ω \cos (ω t)$ with $γ_{0} = 2$ and $ω \in {0.1, 0.2, \dots, 1}$ for $0 \leq t \leq 100$ ⁠, respectively. The color version is available online.

Figures 5(a) and 5(b) show the constitutive equations found by Rheo-SINDy and the test simulation results, respectively. For test simulations shown in Fig. 5(b), we employed the oscillatory shear flow with $γ_{0} = 4$ and $ω = 0.5$ ⁠, which is outside of the parameters in the training data described in Sec. III B. Figure 5(a) reveals that the STRidge with $α = 3 \times 10^{- 1}$ gave almost exact constitutive equations, including the value of $α_{G}$ [cf. Eqs. (20)–(23)]. As inferred from this, the predictions based on the constitutive equations obtained by the STRidge demonstrate a good agreement with the test data, as shown in Fig. 5(b-ii). In contrast to the success of the STRidge, the STLSQ and a-Lasso seemingly failed to identify the correct solution, as indicated in Fig. 5(a). The constitutive equations obtained by the STLSQ with $α = 3 \times 10^{- 1}$ has a low error rate as shown in Fig. 4(b), but its predicted $τ_{x x}$ significantly deviates from the test data as seen in Fig. 5(b-i). In contrast, although the a-Lasso did not provide the correct solution for $τ_{x x}$ ⁠, the test simulations with the obtained constitutive equations exhibit a good agreement with the test data. These test simulations demonstrate that the STRidge and a-Lasso are promising approaches for Rheo-SINDy.

FIG. 5.

View large Download slide

(a) The constitutive equations obtained by Rheo-SINDy with (i) the STLSQ, (ii) STRidge, and (iii) a-Lasso from the multiple data trajectories, and (b) their test simulation results under the oscillatory shear flow with $γ_{0} = 4$ and $ω = 0.5$ ⁠. The training data are the same as those in Fig. 4. The exact equations for the Giesekus model under shear flow are shown in Eqs. (20)–(23). In (b), the $x x$ -, $y y$ -, and $x y$ -components of the stress tensor are shown with black, blue, and red lines, respectively. The dotted and solid lines in (b) denote the predictions by the equations shown in (a) and those by the exact Giesekus model, respectively. The color version is available online.

C. Hookean dumbbell model

We next explain the results for the Hookean dumbbell model. We here used the polynomial library that includes up to second-order terms of $τ_{x x}$ ⁠, $τ_{y y}$ ⁠, $τ_{x y}$ ⁠, and $κ_{x y}$ ⁠. Thus, the total number of candidate terms was $N_{Θ} = 15$ for each component.

Figure 6 shows the Rheo-SINDy results for the Hookean dumbbell model with the different numbers of dumbbells. We note that the standard deviation of $τ$ in the training data decreases proportionally with $N_{p}^{- 1 / 2}$ ⁠. From Fig. 6(a), as the value of $N_{p}$ increases, sparser solutions are obtained, especially for Rheo-SINDy with the STRidge. Unlike the case of the UCM model (cf. Fig. 3), which can be considered as the “noise-free” case of the Hookean dumbbell model, the STRidge provides the correct number of terms only within a narrow range of $α$ values. Nevertheless, if we choose the appropriate $α$ value, the nearly correct constitutive equations can be found by the STRidge, as shown in the upper part of Fig. 6(b). We note that the terms containing $τ_{y y}$ appear in the time evolution equation for $τ_{x y}$ obtained by the STRidge. Although these terms do not affect the predictions because $τ_{y y} = 0$ ⁠, we speculate that the appearance of these terms is due to the correlation effects of the noise in $R_{x}$ and $R_{y}$ on the stress [cf. Eq. (28)]. When comparing the STRidge and a-Lasso, it is evident that the a-Lasso provides stable and sparse solutions across a broader range of $α$ values, regardless of the $N_{p}$ value. Furthermore, we confirm that the correct equations can be obtained using the a-Lasso, as shown in the lower part of Fig. 6(b). Comparing the results obtained by the a-Lasso with different $N_{p}$ ⁠, we found that the correct number of terms is obtained in almost the same $α$ range for $N_{p} \geq 10^{4}$ ⁠. This indicates that using $N_{p} = 10^{4}$ provides sufficient results in this case. This partially suggests the effectiveness of the a-Lasso in discovering essential terms from noisy data.

FIG. 6.

View large Download slide

(a) The total number of terms and (b) the constitutive equations obtained by Rheo-SINDy with the STRidge (black open symbols) and a-Lasso (red closed symbols) for the Hookean dumbbell model [Eqs. (16)–(18)]. In (a), the horizontal line indicates the correct number of terms, and circle, triangle, and square symbols represent the results for $N_{p} = 10^{3}$ ⁠, $10^{4}$ ⁠, and $10^{5}$ ⁠, respectively. The color version is available online.

D. FENE-P dumbbell model

We next examine whether Rheo-SINDy can find more complex differential equations (i.e., the FENE-P dumbbell model) than the UCM model and the Giesekus model. We utilized Rheo-SINDy with

T

replaced by

C

to discover the differential equations for the conformation tensor

C

of the FENE-P dumbbell model explained in Sec. III C 2. In this case, we prepared the following custom library:

Θ = [\begin{matrix} 1 & Ω (t_{1}) & Ω^{2} (t_{1}) & f_{FENE} (t_{1}) Ω (t_{1}) \\ 1 & Ω (t_{2}) & Ω^{2} (t_{2}) & f_{FENE} (t_{2}) Ω (t_{2}) \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & Ω (t_{n}) & Ω^{2} (t_{n}) & f_{FENE} (t_{n}) Ω (t_{n}) \end{matrix}],

(42)

where

Ω

consists of nonzero components of

C

under shear flow (⁠

C_{x x}

⁠,

C_{y y}

⁠,

C_{z z}

⁠, and

C_{x y}

⁠) and

κ_{x y}

⁠, and

Ω^{2}

is the vector composed of all the multiplied combinations of the

Ω

components. The total number of library functions was, thus,

N_{Θ} = 26

⁠.

Figure 7 displays the results: (a) the total number of predicted terms and (b) the error rate as a function of the hyperparameter $α$ ⁠. The error rate is defined as the sum of the mean squared errors (MSEs) of ${\dot{t}}_{μ ν} - Θ {\hat{ξ}}_{μ ν}$ scaled so that the maximum value of each method is 1. Similar to the results for the phenomenological constitutive equations shown in Secs. IV A and IV B, the a-Lasso provides sparser solutions than the STRidge, and the STRidge gives lower error rates than the a-Lasso. Figure 8 presents the differential equations obtained by the STRidge and a-Lasso for two $α$ values: the one selected by our proposed criterion (i.e., with the smallest number of terms, while the corresponding error is of the same order as the minimum error) and the other has a larger $α$ value (i.e., with a smaller number of terms). From the lower part of Fig. 8, while the a-Lasso provides sparser solutions, they are not correct [cf. Eqs. (34)–(36)]. Specifically, in all cases for $C_{x x}$ ⁠, $C_{y y}$ ⁠, and $C_{z z}$ ⁠, the a-Lasso failed to identify the constant term in Eqs. (34) and (35), which is a possible source of larger errors compared to the STRidge. In the case of the STRidge, we confirmed that by choosing $α$ based on our criterion (⁠ $α = 1 \times 10^{- 1}$ ⁠), nearly correct differential equations can be obtained, as shown in the upper part of Fig. 8. Since the $y y$ -component and $z z$ -component of the stress are equivalent, the exact equations can be recovered by setting $C_{y y} = C_{z z}$ ⁠. Thus, we conclude that the correct differential equations for the FENE-P dumbbell model can be obtained by Rheo-SINDy.

FIG. 7.

View large Download slide

(a) The total number of terms and (b) the error rate of the constitutive equations obtained by Rheo-SINDy with the STRidge (black squares) and a-Lasso (red reverse triangles) for the FENE-P dumbbell model in the conformation expression [Eqs. (34)–(36) with $n_{K} = 10$ ]. The horizontal line in (a) indicates the correct number of terms. The blue circles represent the $α$ values selected for test simulations. Based on our criterion, $α = 1 \times 10^{- 1}$ for the STRidge and $α = 3 \times 10^{- 8}$ for the a-Lasso are picked. The color version is available online.

FIG. 8.

View large Download slide

The constitutive equations obtained by Rheo-SINDy with the STRidge and a-Lasso for the FENE-P dumbbell model in the conformation expression [Eqs. (34)–(36) with $n_{K} = 10$ ]. Based on our criterion, the appropriate $α$ values were determined as $α = 1 \times 10^{- 1}$ for the STRidge and $α = 3 \times 10^{- 8}$ for the a-Lasso.

To validate the obtained equations in Fig. 8, we generated data using the equations under conditions different from those used to generate the training data (namely, the oscillatory shear flow with $γ_{0} = 4$ and $ω = 1$ ⁠) and converted these data to the stress data using the dimensionless form of Eq. (33). Figure 9 compares the generated data with the correct data generated from Eqs. (33)–(36). As shown in this figure, the equations obtained by the STRidge can reproduce the exact solutions even when the equations are not exactly correct (⁠ $α = 1 \times 10^{- 3}$ ⁠). In contrast, the test simulations with the differential equations obtained by the a-Lasso show the deviations from the test data, especially for $τ_{x x}$ ⁠. In particular, although the equations with $α = 3 \times 10^{- 8}$ have the small error rate for the training data [cf. Fig. 7(b)], those predictions deviate from the test data, and this deviation increases over time. These results emphasize the need to choose an appropriate optimization method to obtain reasonable solutions.

FIG. 9.

View large Download slide

Test simulation results under the oscillatory shear flow with $γ_{0} = 4$ and $ω = 1$ using the equations obtained by Rheo-SINDy with (a) the STRidge and (b) a-Lasso for the FENE-P dumbbell model in the conformation expression. The obtained data from the equations were converted to the stress using the dimensionless form of Eq. (33). The black, blue, and red lines show $τ_{x x}$ ⁠, $τ_{y y}$ ⁠, and $τ_{x y}$ ⁠, respectively. The bold, thin dotted, and thin solid lines indicate the exact solutions [Eqs. (33)–(36) with $n_{K} = 10$ ], predictions with smaller $α$ values (⁠ $α = 1 \times 10^{- 3}$ for the STRidge and $α = 3 \times 10^{- 8}$ for the a-Lasso), and predictions with larger $α$ values (⁠ $α = 1 \times 10^{- 1}$ for the STRidge and $α = 1 \times 10^{- 4}$ for the a-Lasso), respectively. The color version is available online.

We then examine whether the stress expression of the constitutive equation for the FENE-P dumbbell model [cf. Eqs. (39)–(41)] can be found by Rheo-SINDy. For such a purpose, we prepared the following custom library:

Θ = [\begin{matrix} 1 & {tr τ (t_{1})}^{p} T_{s} (t_{1}) & {tr τ (t_{1})}^{2} & {T_{s} (t_{1})}^{p} κ_{x y} (t_{1}) \\ 1 & {tr τ (t_{2})}^{p} T_{s} (t_{2}) & {tr τ (t_{2})}^{2} & {T_{s} (t_{2})}^{p} κ_{x y} (t_{2}) \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & {tr τ (t_{n})}^{p} T_{s} (t_{n}) & {tr τ (t_{n})}^{2} & {T_{s} (t_{n})}^{p} κ_{x y} (t_{n}) \end{matrix}],

(43)

where

T_{s}

includes

{τ_{x x}, τ_{y y}, τ_{z z}, τ_{x y}}

and

p

(⁠

= 0, 1, 2

⁠) is the polynomial order. Thus, the total number of library functions was

N_{Θ} = 29

⁠. We designed the library to include all terms present in Eqs. . Moreover, we excluded terms that could potentially become large, such as higher-order terms involving

κ_{x y}

⁠. When such terms are included in the solutions, the differential equations may be unstable, and in worse cases, they may diverge. Furthermore, as proposed by Lennon and co-workers, the principle of frame-invariance can also be used to constrain the candidate terms in

Θ

[33], which will be considered in our future work.

Figure 10 shows (a) the total number of terms and (b) the error rate of the constitutive equation obtained by Rheo-SINDy with the STRidge and a-Lasso for the FENE-P dumbbell model. Similar to the previous case, the a-Lasso yields sparser solutions than the STRidge. Based on the number of terms and the error rates shown in Fig. 10, we chose two $α$ ⁠, one is based on the proposed criterion and the other has a smaller number of terms. Figure 11 presents the equations obtained using the chosen $α$ values. From Fig. 11, the equations predicted by the STRidge with $α = 1$ and the a-Lasso with $α = 1 \times 10^{- 4}$ are almost the same; conversely, the expressions for small $α$ values significantly differ between the two methods. For the STRidge with $α = 1 \times 10^{- 3}$ ⁠, the identified equations are close to the correct equations [cf. Eqs. (39)–(41)]. Furthermore, the coefficient values for the terms obtained correctly are close to the correct values. For the a-Lasso with $α = 1 \times 10^{- 8}$ ⁠, several coefficients for the correctly obtained terms, such as $τ_{x x}$ ⁠, $τ_{x y} κ_{x y}$ ⁠, and $τ_{x x} τ_{x y} κ_{x y}$ in the equation for ${\dot{τ}}_{x x}$ ⁠, are close to the exact values, but for other several terms, such as $tr τ τ_{x x}$ in the equation for ${\dot{τ}}_{x x}$ ⁠, the correct coefficient values are not obtained. Nevertheless, from Fig. 12, which shows the test simulation results, the equations obtained by the STRidge with $α = 1 \times 10^{- 3}$ and the a-Lasso with $α = 1 \times 10^{- 8}$ can well reproduce the exact solutions including the small oscillation of $τ_{y y}$ ⁠. Although the equations obtained by the STRidge and a-Lasso demonstrate similar performance in the test simulations shown in Fig. 12, the difference in predictions is quantified by their MSEs shown in Table IV. When $α$ is small, the STRidge outperforms the a-Lasso. The STRidge, however, generally provides sparse solutions within a narrow range of $α$ values, requiring careful selection of $α$ ⁠. In the case of the FENE-P dumbbell model, as shown in Figs. 7 and 10, the error rate obtained by the STRidge is smaller than that obtained by the a-Lasso over the investigated $α$ region. This indicates that the STRidge is superior to the a-Lasso in representing the training data. However, as shown in Sec. IV E, even with a low error rate, the performance in test simulations for the extrapolation region is not guaranteed. The discussion in Sec. IV D indicates that the optimization method needs to be chosen based on the specific problem.

FIG. 10.

View large Download slide

(a) The total number of terms and (b) the error rate of the constitutive equation obtained by Rheo-SINDy with the STRidge (black squares) and the a-Lasso (red reverse triangles) for the FENE-P dumbbell model in the stress expression [Eqs. (39)–(41) with $n_{K} = 10$ ]. The horizontal short-dashed line in (a) indicates that the number of terms is zero. The blue circles represent the $α$ values selected for test simulations. Based on our criterion, the appropriate $α$ values were determined as $α = 1 \times 10^{- 3}$ for the STRidge and $α = 1 \times 10^{- 8}$ for the a-Lasso. The color version is available online.

FIG. 11.

View large Download slide

The constitutive equations obtained by Rheo-SINDy with the STRidge and a-Lasso for the FENE-P dumbbell model in the stress expression [Eqs. (39)–(41) with $n_{K} = 10$ ]. Based on our criterion, the appropriate $α$ values were determined as $α = 1 \times 10^{- 3}$ for the STRidge and $α = 1 \times 10^{- 8}$ for the a-Lasso.

FIG. 12.

View large Download slide

Test simulation results under the oscillatory shear flow with $γ_{0} = 4$ and $ω = 1$ using the constitutive equations obtained by Rheo-SINDy with (a) the STRidge and (b) a-Lasso for the FENE-P dumbbell model in the stress expression. The black, blue, and red lines represent the $x x$ -, $y y$ -, and $x y$ -components of the stress tensor, respectively. The bold lines show the exact solutions. The thin solid and short-dashed lines indicate the results with smaller $α$ values (⁠ $α = 1 \times 10^{- 3}$ for the STRidge and $α = 1 \times 10^{- 8}$ for the a-Lasso) and with larger $α$ values (⁠ $α = 1$ for the STRidge and $α = 1 \times 10^{- 4}$ for the a-Lasso), respectively. The color version is available online.

TABLE IV.

The mean squared errors (MSEs) between the predicted and correct constitutive equations for the FENE-P dumbbell model in the stress expression [Eqs. (39)–(41) with n_K = 10]. .

Method	α	MSE (τ_xx)	MSE (τ_yy)	MSE (τ_xy)
STRidge	1 × 10⁻³	8.1 × 10⁻³	3.1 × 10⁻⁴	1.1 × 10⁻³
STRidge	1	2.3	8.2 × 10⁻³	8.9 × 10⁻²
a-Lasso	1 × 10⁻⁸	3.2 × 10⁻¹	3.4 × 10⁻³	2.1 × 10⁻²
a-Lasso	1 × 10⁻⁴	2.3	8.2 × 10⁻³	9.0 × 10⁻²

Method	α	MSE (τ_xx)	MSE (τ_yy)	MSE (τ_xy)
STRidge	1 × 10⁻³	8.1 × 10⁻³	3.1 × 10⁻⁴	1.1 × 10⁻³
STRidge	1	2.3	8.2 × 10⁻³	8.9 × 10⁻²
a-Lasso	1 × 10⁻⁸	3.2 × 10⁻¹	3.4 × 10⁻³	2.1 × 10⁻²
a-Lasso	1 × 10⁻⁴	2.3	8.2 × 10⁻³	9.0 × 10⁻²

E. FENE dumbbell model

We next address the FENE dumbbell model. As explained in Sec. III C 3, the FENE dumbbell model does not have an analytical expression of the constitutive equation. Thus, we here aim to identify an approximate constitutive equation using Rheo-SINDy.

To prepare the library $Θ$ for the FENE dumbbell model, we utilize the physical insights obtained from the analytical expression of the FENE-P dumbbell model. We assume the constitutive equation for the FENE-P dumbbell model is similar to that for the FENE dumbbell model. Since the FENE-P dumbbell model is a simplified version of the FENE dumbbell model, we believe that this is a reasonable assumption. Here, we note that the stress expression shown in Eq. (33) is no longer applicable to the FENE dumbbell model since the values of $h (t)$ differ for each individual dumbbell. Thus, it is invalid to obtain stress through the conformation tensor $C$ using Eq. (33). Based on the above considerations, we decided to use the custom library presented in Eq. (43), which was also used to discover the constitutive equation for the FENE-P dumbbell model.

Figure 13 compares (a) the total number of terms and (b) the error rate predicted by the STRidge and a-Lasso. Similar to the previous discussions, we obtain sparse solutions over a wide range of $α$ values with the a-Lasso, whereas the STRidge gives sparse solutions only within a limited range of $α$ ⁠. The left table in Fig. 14 shows the equations obtained by the a-Lasso with two $α$ values chosen in the same way as Fig. 10. The predictions obtained by the STRidge are inferior to those obtained by the a-Lasso shown in Fig. 14; therefore, we focus on the a-Lasso result below (the STRidge results are discussed in Sec. S3 in the supplementary material). From the left table in Fig. 14, if $α$ is appropriately chosen, the a-Lasso can give sparse equations with coefficients of reasonable (not excessively large) magnitudes. Comparing the equations for the FENE-P model obtained by the a-Lasso with $α = 1 \times 10^{- 8}$ (Fig. 11) and those for the FENE model obtained by the a-Lasso with $α = 1 \times 10^{- 6}$ (Fig. 14), the appearing terms are almost identical, which demonstrates the similarity between these models. The difference in the coefficients, thus, represents the difference between these models. The right panels in Fig. 14 show the test simulation results obtained by the equations shown in the left table. With the large $α$ (⁠ $α = 3 \times 10^{- 4}$ ⁠), the identified equation for $τ_{y y}$ becomes ${\dot{τ}}_{y y} = 0$ ⁠, which fails to reproduce the oscillatory behavior of $τ_{y y}$ ⁠. On the contrary, the equations obtained with $α = 1 \times 10^{- 6}$ can reproduce well the BD simulation results outside the range of the training data, including the oscillatory behavior of $τ_{y y}$ ⁠. This success suggests that Rheo-SINDy with the a-Lasso is effective for discovering unknown constitutive equations. Nevertheless, we note that the equations presented in Fig. 14 may fail to predict test data significantly outside the range of the training data. Reproducing such highly nonlinear data would require the nonlinear terms dropped in Fig. 14. In this sense, the constitutive equations for the FENE dumbbell model obtained here are appropriately referred to as the approximate constitutive equations.

FIG. 13.

View large Download slide

(a) The total number of terms and (b) the error rate of the constitutive equation obtained by Rheo-SINDy with the STRidge (black squares) and the a-Lasso (red reverse triangles) for the FENE dumbbell model. The horizontal short-dashed line in (a) indicates that the number of terms is zero. The blue circles represent the $α$ values selected for test simulations. Based on our criterion, the appropriate $α$ values were determined as $α = 1 \times 10^{- 1}$ for the STRidge and $α = 1 \times 10^{- 6}$ for the a-Lasso. The color version is available online.

FIG. 14.

View large Download slide

The constitutive equations obtained by Rheo-SINDy with the a-Lasso for the FENE dumbbell model (left) and the test simulation results under the oscillatory shear flows with $γ_{0} = 3$ and $ω = 1$ (right upper panel) and $γ_{0} = 4$ and $ω = 1$ (right lower panel). The black, blue, and red lines represent the $x x$ -, $y y$ -, and $x y$ -components of the stress tensor, respectively. The bold lines show the exact solutions, and the thin solid and short-dashed lines show the results with the smaller $α$ value (⁠ $α = 1 \times 10^{- 6}$ ⁠), which is chosen based on our criterion, and the larger $α$ value (⁠ $α = 3 \times 10^{- 4}$ ⁠). The color version is available online.

To further investigate rheological quantities under shear flow, we conducted shear simulations under constant shear rates. Figure 15 compares the steady state shear viscosity $η \equiv σ_{x y} / \dot{γ}$ and first normal stress coefficient $Ψ_{1} \equiv (σ_{x x} - σ_{y y}) / {\dot{γ}}^{2}$ obtained from BD simulations with the approximate constitutive equations obtained by Rheo-SINDy. Although the constant shear flows considered in Fig. 15 differ from the oscillatory shear flows used to generate the training data, Rheo-SINDy can reproduce the results of BD simulations, including shear thinning, within the shear rate region of the training data and slightly beyond it, as shown in Fig. 15(a). However, in the high shear rate region, the Rheo-SINDy predictions with $α = 1 \times 10^{- 6}$ exhibit plateaus for $η$ and $Ψ_{1}$ that differ from the results of BD simulations. To improve the Rheo-SINDy predictions in the high shear rate region, we tested the training data with larger $γ_{0}$ (⁠ $γ_{0} = 8$ ⁠) while keeping the same $ω$ values as those used to obtain Figs. 13 and 14. We note that the generated data with $γ_{0} = 8$ show nonlinear stress responses, as shown in Sec. S4 in the supplementary material. Figure 15(b) indicates that the model found by Rheo-SINDy with the training data that include higher shear rates successfully extrapolates continuous shear thinning to the high shear rate regime beyond the training data.

FIG. 15.

View large Download slide

Steady state shear viscosity $η$ (black) and normal stress coefficient $Ψ_{1}$ (red) for the training data (a) with $γ_{0} = 2$ and (b) with $γ_{0} = 8$ ⁠. The symbols are the BD simulation data, and the lines are the Rheo-SINDy predictions. The solid and short-dashed lines show the Rheo-SINDy results with the smaller $α$ value [ $α = 1 \times 10^{- 6}$ for both (a) and (b)] and the larger $α$ value [ $α = 3 \times 10^{- 4}$ for (a) and $α = 1 \times 10^{- 5}$ for (b)]. The shaded regions show the ranges of shear rates in the training data. The color version is available online.

It would be interesting to compare the results obtained by Rheo-SINDy for the FENE-P and FENE dumbbell models. Figure 16 indicates the transient shear viscosity $η^{+}$ and first normal stress coefficient $Ψ_{1}^{+}$ obtained by Rheo-SINDy and those exact solutions for the FENE-P and FENE dumbbell models. We note that, for comparison, the Rheo-SINDy regressions were conducted using the training data with the same parameters (⁠ $γ_{0} = 8$ and $ω \in {0.1, 0.2, \dots, 1}$ ⁠) for both models. For the FENE-P dumbbell model, the regression for the stress expression [cf. Eqs. (39)–(41)] was performed. From the success of the STRidge in this case (cf. Fig. 10), we used the STRidge with $α = 1 \times 10^{- 4}$ ⁠. The Rheo-SINDy results for the FENE dumbbell model were obtained using the a-Lasso with $α = 1 \times 10^{- 5}$ ⁠, which is the same parameter as that determined in Fig. 15(b). Figure 16 shows that Rheo-SINDy for the FENE-P dumbbell model accurately reproduces the transient behavior shown in the exact solutions. However, while Rheo-SINDy for the FENE dumbbell model successfully reproduces the steady state values of $η^{+}$ and $Ψ_{1}^{+}$ ⁠, as expected from Fig. 15(b), it fails to describe the stress overshoot, especially for the large $\dot{γ}$ ⁠. This finding highlights the requirements for a more sophisticated learning strategy for Rheo-SINDy to capture such a highly nonlinear behavior. Although the objective of this study is not to develop a sophisticated constitutive model that approximates the FENE dumbbell model, this discrepancy can be addressed through three possible strategies: (1) adding training data obtained by steady shear flow with high shear rate, (2) using the conformation tensor $C$ as a latent variable, and (3) incorporating higher-order terms into the library for regression. For point (1), it has been reported that including training data obtained from constant shear flow with high shear rate is effective [38]. This strategy enables the development of regression models that can predict transient behavior, even for more complex models. For point (2), the conformation tensor often allows for simpler descriptions, as shown for the FENE-P dumbbell model. It might be effective to use $C$ to obtain regression equations in the FENE dumbbell model. This case will require an additional regression to relate $C$ and $τ$ ⁠, replacing Eq. (33) for the FENE-P dumbbell model. For point (3), Du and co-workers have reported that the approximated model with high-order terms achieves high predictive performance [49]. Following this study, including higher-order terms in the library can enhance the representational ability of Rheo-SINDy. These considerations are interesting subjects for future studies.

FIG. 16.

(a) Transient state shear viscosity η + and (b) first normal stress coefficient Ψ 1 + under the steady shear flow with γ ˙ = 1 (solid lines) and 10 (dotted lines) for the FENE-P (black) and FENE (red) dumbbell models. The thin and bold lines were obtained by Rheo-SINDy and by the equations for the FENE-P and FENE dumbbell models, respectively. The Rheo-SINDy regressions were conducted for the training data using the same parameters ( γ 0 = 8 and ω ∈ { 0.1 , 0.2 , … , 1 }) for both the FENE-P and FENE dumbbell models. While the Rheo-SINDy results for the FENE-P dumbbell model were obtained by the STRidge with α = 1 × 10 − 4, those for the FENE dumbbell model were obtained by the a-Lasso with α = 1 × 10 − 5, which is the same parameter as that in Fig. 15(b). The color version is available online.

View large Download slide

(a) Transient state shear viscosity $η^{+}$ and (b) first normal stress coefficient $Ψ_{1}^{+}$ under the steady shear flow with $\dot{γ} = 1$ (solid lines) and $10$ (dotted lines) for the FENE-P (black) and FENE (red) dumbbell models. The thin and bold lines were obtained by Rheo-SINDy and by the equations for the FENE-P and FENE dumbbell models, respectively. The Rheo-SINDy regressions were conducted for the training data using the same parameters (⁠ $γ_{0} = 8$ and $ω \in {0.1, 0.2, \dots, 1}$ ⁠) for both the FENE-P and FENE dumbbell models. While the Rheo-SINDy results for the FENE-P dumbbell model were obtained by the STRidge with $α = 1 \times 10^{- 4}$ ⁠, those for the FENE dumbbell model were obtained by the a-Lasso with $α = 1 \times 10^{- 5}$ ⁠, which is the same parameter as that in Fig. 15(b). The color version is available online.

Thanks to the equations obtained using Rheo-SINDy, it is possible to provide a physical interpretation with the assistance of rheological knowledge. For example, from the comparison of the equations obtained for the FENE-P dumbbell model (cf. Fig. 11) and those for the FENE dumbbell model (cf. Fig. 14), the equations for larger $α$ value (⁠ $α = 1 \times 10^{- 4}$ for the FENE-P dumbbell model and $α = 3 \times 10^{- 4}$ for the FENE dumbbell model) are similar except for the coefficient values. We note that these $α$ values are not selected based on our criteria. However, to simply analyze the representation obtained by Rheo-SINDy, we consider the sparser solution here. The terms in these equations are the same as those for the UCM model (and, thus, the Hookean dumbbell model). This indicates that all of these models share the same origin based on the dumbbell model. We can also discuss the linear term of stress, which represents the relaxation of stress [see Eq. (15)]. Since the relaxation time at equilibrium (⁠ $λ = ζ / 4 h_{eq}$ ⁠) is taken as the unit time in this study, the coefficient of this term should be $- 1$ at equilibrium [and, thus, for the UCM model, see Eqs. (16)–(18)]. From Figs. 11 and 14, the coefficient of the linear term of stress with the larger $α$ values is smaller than $- 1$ ⁠, which indicates $λ_{sf} < λ_{eq}$ with the subscript “sf” and “eq” standing for “shear flow” and “equiliblium,” respectively. This indicates that under shear flow, the values of spring strength for the FENE-P and FENE dumbbell models become larger than $h_{eq}$ ⁠, which implies the appearance of the FENE effects under flow. From this discussion, it is evident that Rheo-SINDy can provide physically interpretable constitutive equations.

F. Material frame indifference

Finally, we show a method to obtain constitutive models that address the principle of material objectivity. This principle requires that constitutive equations are independent of rotation. However, the results obtained by the Rheo-SINDy regression method presented, thus far do not explicitly satisfy the rotational indifference.

To account for this indifference, Young and co-workers demonstrated that it is effective to separate the rotational components within the velocity gradient tensor from constitutive equations [50]. In this approach, the velocity gradient tensor

κ

in the laboratory frame (LF) is decomposed into the deformation rate tensor

D

and the rotation rate tensor

W

as

κ = D + W

⁠, where

D = (κ + κ^{T}) / 2

and

W = (κ - κ^{T}) / 2

⁠. We then consider a constitutive model in the reference frame (RF) rotated by

W

⁠. Here, we note that the constitutive model in the RF should depend only on

D

⁠. We can compute the rotation matrix

R (θ_{W})

around an axis that relates the LF and RF as

\frac{d R (θ_{W})}{d t} = W \cdot R (θ_{W}),

(44)

where

θ_{W}

is the rotation angle. Using the rotation matrix

R (θ_{W})

⁠, the stress tensor

τ

and the strain-rate tensor

D

are, respectively, transformed into

τ^{'} = R^{T} (θ_{W}) τ R (θ_{W})

and

D^{'} = R^{T} (θ_{W}) D R (θ_{W})

⁠, where prime denotes quantities in the RF. We note that Young and co-workers demonstrated that the constitutive equation of the UCM model in the RF can be written similarly to that in the LF using

τ^{'}

and

D^{'}

[50].

Based on the study of Young and co-workers [50], we here propose a more rigorous approach to account for material objectivity. Note that Lennon and co-workers proposed another method to employ prescribed tensor terms satisfying the frame indifference [33]. Since the current Rheo-SINDy method does not address constraints of inter-coefficients for tensorial terms, the identified equations should be constructed using invariants of the stress and strain-rate tensors. In this approach,

τ^{'}

and

D^{'}

in the RF are diagonalized using rotation matrices around an axis

P (θ_{τ}^{*})

and

P (θ_{D}^{*})

as

τ^{*} = P^{T} (θ_{τ}^{*}) τ^{'} P (θ_{τ}^{*})

and

D^{*} = P^{T} (θ_{D}^{*}) D^{'} P (θ_{D}^{*})

⁠, respectively. Here,

θ_{τ}^{*}

and

θ_{D}^{*}

are the rotation angles for

τ^{*}

and

D^{*}

⁠, respectively. With

τ^{*}

⁠,

D^{*}

⁠, and

ϕ^{*} = θ_{τ}^{*} - θ_{D}^{*}

⁠, the constitutive equations are expressed as follows:

{\dot{τ}}^{*} = f_{τ^{*}} (τ^{*}, D^{*}, ϕ^{*}),

(45)

{\dot{θ}}_{τ}^{*} = f_{θ_{τ}^{*}} (τ^{*}, D^{*}, ϕ^{*}) .

(46)

Here,

ϕ^{*}

indicates the adjustment for using rotation matrices with different rotation angles to obtain

τ^{*}

and

D^{*}

⁠. This type of constitutive equations can cover constitutive relations under deformations in a plane. In the supplementary material (Sec. S5), we demonstrate the derivation of constitutive equations for

τ^{*}

and

θ_{τ}^{*}

in the UCM and Giesekus models. Hereafter, we show that Rheo-SINDy can give Eqs. and for the Giesekus model. From

τ^{*}

and

θ_{τ}^{*}

⁠,

τ

in LF can be determined by the relationship

τ = R (θ_{W}) P (θ_{τ}^{*}) τ^{*} P^{T} (θ_{τ}^{*}) R^{T} (θ_{W})

⁠.

Using the same training data shown in Sec. III B, we here conduct the Rheo-SINDy regressions. In the Giesekus model, the constitutive equations for

τ^{*}

and

θ_{τ}^{*}

are as follows (see the supplementary material for the detailed derivations):

\begin{aligned} {\dot{τ}}_{x x}^{*} = + 2 τ_{x x}^{*} D_{x x}^{*} \cos 2 ϕ^{*} - τ_{x x}^{*} - α_{G} {(τ_{x x}^{*})}^{2} + 2 D_{x x}^{*} \cos 2 ϕ^{*}, \end{aligned}

(47)

\begin{aligned} {\dot{τ}}_{y y}^{*} = - 2 τ_{y y}^{*} D_{x x}^{*} \cos 2 ϕ^{*} - τ_{y y}^{*} - α_{G} {(τ_{y y}^{*})}^{2} - 2 D_{x x}^{*} \cos 2 ϕ^{*}, \end{aligned}

(48)

{\dot{τ}}_{z z}^{*} = - τ_{z z}^{*} - α_{G} {(τ_{z z}^{*})}^{2} = 0,

(49)

{\dot{θ}}_{τ}^{*} = - (\frac{τ_{x x}^{*}}{τ_{x x}^{*} - τ_{y y}^{*}} + \frac{τ_{y y}^{*}}{τ_{x x}^{*} - τ_{y y}^{*}} + \frac{2}{τ_{x x}^{*} - τ_{y y}^{*}}) D_{x x}^{*} \sin 2 ϕ^{*} .

(50)

Similarly to Eqs. , Eqs. are nondimensionalized using

λ

and

G

⁠. We prepared polynomial libraries up to second order consisting of

{τ_{x x}^{*}, τ_{y y}^{*}, τ_{z z}^{*}, D_{x x}^{*} \sin (2 ϕ^{*}), D_{x x}^{*} \cos (2 ϕ^{*})}

for the regression of

τ^{*}

⁠, and

{1 / Δ, τ_{x x}^{*} / Δ, τ_{y y}^{*} / Δ, τ_{z z}^{*} / Δ, D_{x x}^{*} \sin (2 ϕ^{*}), D_{x x}^{*} \cos (2 ϕ^{*})}

with the shorthand

Δ = τ_{x x}^{*} - τ_{y y}^{*}

for the regression of

θ_{τ}^{*}

⁠. Thus, there were

N_{Θ} = 21

and

N_{Θ} = 28

candidate terms for each component of

τ^{*}

and

θ_{τ}^{*}

⁠, respectively. Based on dimensional consideration, it is natural to create the library for the regression of

θ_{τ}^{*}

using terms of

τ^{*} / Δ

⁠.

Figure 17 indicates the Rheo-SINDy results for $τ^{*}$ and $θ_{τ}^{*}$ ⁠. The trends in the number of terms and error rate shown in Figs. 17(a) and 17(b) are the same as those in the regressions with Rheo-SINDy in the LF (cf. Fig. 4). The lower panel in Fig. 17(a) shows that the a-Lasso has a larger error rate than the STRidge, reflecting the inaccuracies in the regression, especially in $τ_{x x}^{*}$ ⁠, as shown in Fig. 17(c). Nevertheless, in the regressions of $θ_{τ}^{*}$ ⁠, the STRidge and a-Lasso have similar error rates, which indicates that these methods have comparable performance in this case. We note that the expression of ${\dot{θ}}_{τ}^{*}$ obtained by the a-Lasso, ${\dot{θ}}_{τ}^{*} ≃ - (1 + 2 τ_{y y}^{*} / Δ + 2 / Δ) D_{x x}^{*} \sin 2 ϕ^{*}$ ⁠, reduces the correct equation shown in Eq. (50). The Rheo-SINDy regressions for $τ^{*}$ using the STRidge have small error rates, but the obtained expressions shown in Fig. 17(c) differ from Eqs. (47) and (48). This is because the Giesekus model with $α_{G} = 0.5$ is equivalent to the Leonov model [16]. Interestingly, by conducting the diagonalization, the equivalent mathematical expression becomes evident through the Rheo-SINDy regression.

FIG. 17.

View large Download slide

The total number of terms and the error rate of the constitutive equations obtained by Rheo-SINDy with the STRidge (red reverse triangles) and the a-Lasso (black triangles) for (a) $τ^{*}$ and for (b) $θ_{τ}^{*}$ ⁠. The blue circles represent the appropriate $α$ values. The horizontal lines in the top figures indicate the correct number of terms. (c) The constitutive equations obtained by Rheo-SINDy with the STRidge and the a-Lasso. The color version is available online.

Through the discussion in this section, we have shown that a constitutive equation that takes into account material objectivity can be obtained from Rheo-SINDy. This observation indicates that the technical preparations are now in place to precisely handle more complex models using Rheo-SINDy.

V. CONCLUDING REMARKS

We tested whether the sparse identification for nonlinear dynamics (SINDy) modified for nonlinear rheological data, which we call Rheo-SINDy, is effective in finding constitutive equations of complex fluid dynamics. We found that Rheo-SINDy can successfully identify correct equations from training data generated from known constitutive equations, as well as provide approximate constitutive equations (or reduced order models) from training data generated by mesoscopic models when constitutive equations are analytically unknown.

Rheo-SINDy was applied for two phenomenological constitutive equations (i.e., the upper-convected Maxwell model and Giesekus model), which revealed the following two things. First, compared to constant shear tests, oscillatory shear tests are appropriate for generating training data. Second, the sequentially thresholded Ridge regression (STRidge) and adaptive-Lasso (a-Lasso) are effective in finding appropriate constitutive equations. We then examined the commonly used mesoscopic model, namely, the dumbbell model with three different representations of spring strength: the Hookean, FENE-P, and FENE springs. Although the Hookean and FENE-P dumbbell models have analytical constitutive equations, for the FENE dumbbell model, there is no analytical expression of the constitutive equation. We confirmed through the Hookean dumbbell model that even in the presence of noise, the a-Lasso provides the correct solution over a wide range of the hyperparameter $α$ ⁠. Rheo-SINDy was also effective in discovering the complex constitutive equations for the FENE-P dumbbell model. This case study revealed that the identification of complex equations requires the preparation of an appropriate custom library based on prior physical knowledge. Using physical insights obtained from the Hookean and FENE-P dumbbell models, we attempted to find approximate constitutive equations for the FENE dumbbell model as combinations of known elemental terms given in a function library. While the approximate equations provided by the a-Lasso successfully predicted steady state shear properties beyond the range of the training data, these equations did not accurately capture transient properties. This finding highlighted the necessity of improving the learning strategy (e.g., how to generate a training dataset) to identify a reasonable equation for a highly nonlinear regime from mesoscopic models without closed-form constitutive equations. Finally, we present a regression method that yields a constitutive model satisfying rotational indifference. These results indicate that Rheo-SINDy is ready to be applied to data obtained from more complex models or experimental data.

From our investigation, Rheo-SINDy with the STRidge or a-Lasso is effective for discovering constitutive equations from nonlinear rheological data. We found that the STRidge is generally superior in terms of retaining correct terms, while the a-Lasso is more robust to the selection of $α$ than the STRidge. To obtain correct constitutive equations, in addition to selecting the appropriate optimization method, we are required to design an appropriate library by using physical insights, namely, domain knowledge. Designing such a proper library necessitates not only including necessary terms but also excluding unnecessary terms. For such a purpose, candidate terms can be chosen to satisfy the principle of frame-invariance [33]. Furthermore, we may conduct regression with physics-informed constraints on the coefficients [51]. However, we cannot entirely exclude the possibility that the proposed Rheo-SINDy alone may fail in regressing models without closed-form constitutive equations or in addressing real-life situations in rheology. Combining Rheo-SINDy with a statistical model (e.g., the Gaussian process regression) may effectively address such complexities.

This research is expected to have an impact on fields such as rheology and fluid mechanics. From a rheological perspective, for several systems such as entangled polymers [52,53] and wormlike micellar solutions [54,55], sophisticated mesoscopic models suitable for numerical simulations under flow have been proposed. These mesoscopic models can generate reasonable training data not only under shear flow but also under extensional flow. Finding new approximate models from the data obtained by these mesoscopic simulations would be an interesting research subject. Furthermore, it would be desirable to conduct Rheo-SINDy for experimental data obtained by large amplitude oscillatory shear (LAOS) experiments [56]. Since the LAOS measurements do not provide all the major stress components under shear flow, exploring methods for discovering the constitutive equations from experimental data would be a future challenge. When approximate constitutive models are identified, those models can be employed for predictions of complex flows, which would deepen our understanding of complex fluids. We will continue our research in these directions.

SUPPLEMENTARY MATERIAL

See the supplementary material for the discussion on the hyperparameter of the a-Lasso (Sec. S1), the derivation of the stress expression for the FENE-P dumbbell model (Sec. S2), the STRidge results for the FENE dumbbell model (Sec. S3), the regression results with larger γ 0 (Sec. S4), and the derivation of the equations to account for material frame indifference (Sec. S5).

ACKNOWLEDGMENTS

We thank Professors Yoshinobu Kawahara and Koji Fukagata for their valuable comments on data-driven methods. We also thank Dr. John J. Molina for carefully reading the manuscript and providing insightful comments. We are grateful to Professor Takashi Taniguchi for encouraging us to start this research. T.S. was partially supported by JST PRESTO (Grant No. JPMJPR22O3). S.M. was financially supported by the Kyoto University Science and Technology Innovation Creation Fellowship (FSMAT) (Grant No. JPMJFS2123).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

T.S. and S.M. contributed equally to this work.

DATA AVAILABILITY

The data that support the findings of this study are available upon reasonable request from the authors.

REFERENCES

1.

Brunton

,

S. L.

, and

J. M.

Kutz

,

Data-Driven Science and Engineering

, 2nd ed. (

Cambridge University

,

Cambridge

,

2022

).

Google Scholar

Crossref

2.

Brunton

,

S. L.

,

J. L.

Proctor

, and

J. N.

Kutz

, “

Discovering governing equations from data by sparse identification of nonlinear dynamical systems

,”

Proc. Natl. Acad. Sci. U. S. A.

113

,

3932

–

3937

(

2016

).

https://doi.org/10.1073/pnas.1517384113

Google Scholar

Crossref

PubMed

3.

de Silva

,

B. M.

,

D. M.

Higdon

,

S. L.

Brunton

, and

J. N.

Kutz

, “

Discovery of physics from data: Universal laws and discrepancies

,”

Front. Artif. Intell.

3

,

25

(

2020

).

https://doi.org/10.3389/frai.2020.00025

Google Scholar

Crossref

PubMed

4.

Kaheman

,

K.

,

S. L.

Brunton

, and

J.

Nathan Kutz

, “

Automatic differentiation to simultaneously identify nonlinear dynamics and extract noise probability distributions from data

,”

Mach. Learn. Sci. Technol.

3

,

015031

(

2022

).

https://doi.org/10.1088/2632-2153/ac567a

Google Scholar

Crossref

5.

Fasel

,

U.

,

J. N.

Kutz

,

B. W.

Brunton

, and

S. L.

Brunton

, “

Ensemble-SINDy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control

,”

Proc. R. Soc. A

478

,

20210904

(

2022

).

https://doi.org/10.1098/rspa.2021.0904

Google Scholar

Crossref

6.

Schmidt

,

M.

, and

H.

Lipson

, “

Distilling free-form natural laws from experimental data

,”

Science

324

,

81

–

85

(

2009

).

https://doi.org/10.1126/science.1165893

Google Scholar

Crossref

PubMed

7.

Bongard

,

J.

, and

H.

Lipson

, “

Automated reverse engineering of nonlinear dynamical systems

,”

Proc. Natl. Acad. Sci. U. S. A.

104

,

9943

–

9948

(

2007

).

https://doi.org/10.1073/pnas.0609476104

Google Scholar

Crossref

PubMed

8.

Reinbold

,

P. A. K.

,

L. M.

Kageorge

,

M. F.

Schatz

, and

R. O.

Grigoriev

, “

Robust learning from noisy, incomplete, high-dimensional experimental data via physically constrained symbolic regression

,”

Nat. Commun.

12

,

3219

(

2021

).

https://doi.org/10.1038/s41467-021-23479-0

Google Scholar

Crossref

PubMed

9.

Udrescu

,

S.-M.

, and

M.

Tegmark

, “

AI Feynman: A physics-inspired method for symbolic regression

,”

Sci. Adv.

6

,

eaay2631

(

2020

).

https://doi.org/10.1126/sciadv.aay2631

Google Scholar

Crossref

PubMed

10.

Cranmer

,

M.

,

A.

Sanchez Gonzalez

,

P.

Battaglia

,

R.

Xu

,

K.

Cranmer

,

D.

Spergel

, and

S.

Ho

, “Discovering symbolic models from deep learning with inductive biases,” in NeurIPS, edited by H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin (Curran Associates Inc., Vancouver, BC, Canada, 2020), Vol. 33, pp. 17429–17442.

11.

Lemos

,

P.

,

N.

Jeffrey

,

M.

Cranmer

,

S.

Ho

, and

P.

Battaglia

, “

Rediscovering orbital mechanics with machine learning

,”

Mach. Learn. Sci. Technol.

4

,

045002

(

2023

).

https://doi.org/10.1088/2632-2153/acfa63

Google Scholar

Crossref

12.

Karniadakis

,

G. E.

,

I. G.

Kevrekidis

,

L.

Lu

,

P.

Perdikaris

,

S.

Wang

, and

L.

Yang

, “

Physics-informed machine learning

,”

Nat. Rev. Phys.

3

,

422

–

440

(

2021

).

https://doi.org/10.1038/s42254-021-00314-5

Google Scholar

Crossref

13.

Raissi

,

M.

,

P.

Perdikaris

, and

G. E.

Karniadakis

, “

Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations

,”

J. Comput. Phys.

378

,

686

–

707

(

2019

).

https://doi.org/10.1016/j.jcp.2018.10.045

Google Scholar

Crossref

14.

Jia

,

X.

,

J.

Willard

,

A.

Karpatne

,

J. S.

Read

,

J. A.

Zwart

,

M.

Steinbach

, and

V.

Kumar

, “

Physics-guided machine learning for scientific discovery: An application in simulating lake temperature profiles

,”

ACM/IMS Trans. Data Sci.

2

,

1

–

26

(

2021

).

https://doi.org/10.1145/3447814

Google Scholar

Crossref

15.

Rosofsky

,

S. G.

,

H.

Al Majed

, and

E. A.

Huerta

, “

Applications of physics informed neural operators

,”

Mach. Learn. Sci. Technol.

4

,

025022

(

2023

).

https://doi.org/10.1088/2632-2153/acd168

Google Scholar

Crossref

16.

Larson

,

R. G.

, Constitutive Equations for Polymer Melts and Solutions, Butterworths Series in Chemical Engineering (Butterworths, Stoneham, MA, 1988).

17.

Ilg

,

P.

, and

M.

Kröger

, “

Molecularly derived constitutive equation for low-molecular polymer melts from thermodynamically guided simulation

,”

J. Rheol.

55

,

69

–

93

(

2011

).

https://doi.org/10.1122/1.3523485

Google Scholar

Crossref

18.

Bird

,

R. B.

,

R. C.

Armstrong

, and

O.

Hassager

,

Dynamics of Polymeric Liquids

, 2nd ed. (

Oxford University

,

New York

,

1987

), Vol. 2.

Google Scholar

19.

Rouse

,

P. E.

, “

A theory of the linear viscoelastic properties of dilute solutions of coiling polymers

,”

J. Chem. Phys.

21

,

1272

–

1280

(

1953

).

https://doi.org/10.1063/1.1699180

Google Scholar

Crossref

20.

Doi

,

M.

, and

S. F.

Edwards

,

The Theory of Polymer Dynamics

(

Oxford University

,

New York

,

1986

).

Google Scholar

21.

Masubuchi

,

Y.

,

J.-I.

Takimoto

,

K.

Koyama

,

G.

Ianniruberto

,

G.

Marrucci

, and

F.

Greco

, “

Brownian simulations of a network of reptating primitive chains

,”

J. Chem. Phys.

115

,

4387

–

4394

(

2001

).

https://doi.org/10.1063/1.1389858

Google Scholar

Crossref

22.

Doi

,

M.

, and

J.

Takimoto

, “

Molecular modelling of entanglement

,”

Philos. Trans. R. Soc. A

361

,

641

–

652

(

2003

).

https://doi.org/10.1098/rsta.2002.1168

Google Scholar

Crossref

23.

Likhtman

,

A. E.

, “

Single-chain slip-link model of entangled polymers: Simultaneous description of neutron spin-echo, rheology, and diffusion

,”

Macromolecules

38

,

6128

–

6139

(

2005

).

https://doi.org/10.1021/ma050399h

Google Scholar

Crossref

24.

Jamali

,

S.

, “

Data-driven rheology: could this be a new paradigm?

,”

Rheol. Bull.

92

,

20

–

24

(

2023

).

Google Scholar

25.

Miyamoto

,

S.

, “

Short review on machine learning-based multi-scale simulation in rheology

,”

Nihon Reoroji Gakkaishi (J. Soc. Rheol. Jpn.)

52

,

15

–

19

(

2024

).

https://doi.org/10.1678/rheology.52.15

Google Scholar

Crossref

26.

Fang

,

L.

,

P.

Ge

,

L.

Zhang

,

W.

E

, and

H.

Lei

, “

DeePN

^{2}

⁠: A deep learning-based non-Newtonian hydrodynamic model

,”

J. Mach. Learn.

1

,

114

–

140

(

2022

).

https://doi.org/10.4208/jml.220115

Google Scholar

Crossref

27.

Mahmoudabadbozchelou

,

M.

,

K. M.

Kamani

,

S. A.

Rogers

, and

S.

Jamali

, “

Digital rheometer twins: Learning the hidden rheology of complex fluids through rheology-informed graph neural networks

,”

Proc. Natl. Acad. Sci. U. S. A.

119

,

e2202234119

(

2022

).

https://doi.org/10.1073/pnas.2202234119

Google Scholar

Crossref

PubMed

28.

Jin

,

H.

,

S.

Yoon

,

F. C.

Park

, and

K. H.

Ahn

, “

Data-driven constitutive model of complex fluids using recurrent neural networks

,”

Rheol. Acta

62

,

569

–

586

(

2023

).

https://doi.org/10.1007/s00397-023-01405-z

Google Scholar

Crossref

29.

Mahmoudabadbozchelou

,

M.

, and

S.

Jamali

, “

Rheology-informed neural networks (RhINNs) for forward and inverse metamodelling of complex fluids

,”

Sci. Rep.

11

,

12015

(

2021

).

https://doi.org/10.1038/s41598-021-91518-3

Google Scholar

Crossref

PubMed

30.

Mahmoudabadbozchelou

,

M.

,

G. E.

Karniadakis

, and

S.

Jamali

, “

nn-PINNs: Non-newtonian physics-informed neural networks for complex fluid modeling

,”

Soft Matter

18

,

172

–

185

(

2022

).

https://doi.org/10.1039/D1SM01298C

Google Scholar

Crossref

31.

Saadat

,

M.

,

M.

Mahmoudabadbozchelou

, and

S.

Jamali

, “

Data-driven selection of constitutive models via rheology-informed neural networks (RhINNs)

,”

Rheol. Acta

61

,

721

–

732

(

2022

).

https://doi.org/10.1007/s00397-022-01357-w

Google Scholar

Crossref

32.

Mahmoudabadbozchelou

,

M.

,

M.

Caggioni

,

S.

Shahsavari

,

W. H.

Hartt

,

G.

Em Karniadakis

, and

S.

Jamali

, “

Data-driven physics-informed constitutive metamodeling of complex fluids: A multifidelity neural network (MFNN) framework

,”

J. Rheol.

65

,

179

–

198

(

2021

).

https://doi.org/10.1122/8.0000138

Google Scholar

Crossref

33.

Lennon

,

K. R.

,

G. H.

McKinley

, and

J. W.

Swan

, “

Scientific machine learning for modeling and simulating complex fluids

,”

Proc. Natl. Acad. Sci. U. S. A.

120

,

e2304669120

(

2023

).

https://doi.org/10.1073/pnas.2304669120

Google Scholar

Crossref

PubMed

34.

Zhao

,

L.

,

Z.

Li

,

B.

Caswell

,

J.

Ouyang

, and

G. E.

Karniadakis

, “

Active learning of constitutive relation from mesoscopic dynamics for macroscopic modeling of non-newtonian flows

,”

J. Comput. Phys.

363

,

116

–

127

(

2018

).

https://doi.org/10.1016/j.jcp.2018.02.039

Google Scholar

Crossref

35.

Zhao

,

L.

,

Z.

Li

,

Z.

Wang

,

B.

Caswell

,

J.

Ouyang

, and

G. E.

Karniadakis

, “

Active- and transfer-learning applied to microscale-macroscale coupling to simulate viscoelastic flows

,”

J. Comput. Phys.

427

,

110069

(

2021

).

https://doi.org/10.1016/j.jcp.2020.110069

Google Scholar

Crossref

36.

Seryo

,

N.

,

T.

Sato

,

J. J.

Molina

, and

T.

Taniguchi

, “

Learning the constitutive relation of polymeric flows with memory

,”

Phys. Rev. Res.

2

,

033107

(

2020

).

https://doi.org/10.1103/PhysRevResearch.2.033107

Google Scholar

Crossref

37.

Seryo

,

N.

,

J. J.

Molina

, and

T.

Taniguchi

, “

Select applications of Bayesian data analysis and machine learning to flow problems

,”

Nihon Reoroji Gakkaishi (J. Soc. Rheol. Jpn.)

49

,

97

–

113

(

2021

).

https://doi.org/10.1678/rheology.49.97

Google Scholar

Crossref

38.

Miyamoto

,

S.

,

J. J.

Molina

, and

T.

Taniguchi

, “

Machine-learned constitutive relations for multi-scale simulations of well-entangled polymer melts

,”

Phys. Fluids

35

,

063113

(

2023

).

https://doi.org/10.1063/5.0156272

Google Scholar

Crossref

39.

Fukami

,

K.

,

T.

Murata

,

K.

Zhang

, and

K.

Fukagata

, “

Sparse identification of nonlinear dynamics with low-dimensionalized flow representations

,”

J. Fluid Mech.

926

,

A10

(

2021

).

https://doi.org/10.1017/jfm.2021.697

Google Scholar

Crossref

40.

Mahmoudabadbozchelou

,

M.

,

K. M.

Kamani

,

S. A.

Rogers

, and

S.

Jamali

, “

Unbiased construction of constitutive relations for soft materials from experiments via rheology-informed neural networks

,”

Proc. Natl. Acad. Sci. U. S. A.

121

,

e2313658121

(

2024

).

https://doi.org/10.1073/pnas.2313658121

Google Scholar

Crossref

PubMed

41.

Chartrand

,

R.

, “

Numerical differentiation of noisy, nonsmooth data

,”

ISRN Appl. Math.

2011

,

164564

(

2011

).

https://doi.org/10.5402/2011/164564

Google Scholar

Crossref

42.

Rudy

,

S. H.

,

S. L.

Brunton

,

J. L.

Proctor

, and

J. N.

Kutz

, “

Data-driven discovery of partial differential equations

,”

Sci. Adv.

3

,

e1602614

(

2017

).

https://doi.org/10.1126/sciadv.1602614

Google Scholar

Crossref

PubMed

43.

Pedregosa

,

F.

,

G.

Varoquaux

,

A.

Gramfort

,

V.

Michel

,

B.

Thirion

,

O.

Grisel

,

M.

Blondel

,

P.

Prettenhofer

,

R.

Weiss

,

V.

Dubourg

,

J.

Vanderplas

,

A.

Passos

,

D.

Cournapeau

,

M.

Brucher

,

M.

Perrot

, and

E.

Duchesnay

, “

Scikit-learn: Machine learning in Python

,”

J. Mach. Learn. Res.

12

,

2825

–

2830

(

2011

).

Google Scholar

44.

Zou

,

H.

, “

The adaptive lasso and its oracle properties

,”

J. Am. Stat. Assoc.

101

,

1418

–

1429

(

2006

).

https://doi.org/10.1198/016214506000000735

Google Scholar

Crossref

45.

Giesekus

,

H.

, “

A simple constitutive equation for polymer fluids based on the concept of deformation-dependent tensorial mobility

,”

J. Non-Newtonian Fluid Mech.

11

,

69

–

109

(

1982

).

https://doi.org/10.1016/0377-0257(82)85016-7

Google Scholar

Crossref

46.

Laso

,

M.

, and

H.

Öttinger

, “

Calculation of viscoelastic flow using molecular models: The connffessit approach

,”

J. Non-Newtonian Fluid Mech.

47

,

1

–

20

(

1993

).

https://doi.org/10.1016/0377-0257(93)80042-A

Google Scholar

Crossref

47.

Graham

,

M. D.

, “

Drag reduction and the dynamics of turbulence in simple and complex fluids

,”

Phys. Fluids

26

,

101301

(

2014

).

https://doi.org/10.1063/1.4895780

Google Scholar

Crossref

48.

Mochimaru

,

Y.

, “

Unsteady-state development of plane Couette flow for viscoelastic fluids

,”

J. Non-Newtonian Fluid Mech.

12

,

135

–

152

(

1983

).

https://doi.org/10.1016/0377-0257(83)80034-2

Google Scholar

Crossref

49.

Du

,

Q.

,

C.

Liu

, and

P.

Yu

, “

FENE dumbbell model and its several linear and nonlinear closure approximations

,”

Multiscale Model. Simul.

4

,

709

–

731

(

2005

).

https://doi.org/10.1137/040612038

Google Scholar

Crossref

50.

Young

,

C. D.

,

P. T.

Corona

,

A.

Datta

,

M. E.

Helgeson

, and

M. D.

Graham

, “

Scattering-informed microstructure prediction during lagrangian evolution (SIMPLE)—A data-driven framework for modeling complex fluids in flow

,”

Rheol. Acta

62

,

587

–

604

(

2023

).

https://doi.org/10.1007/s00397-023-01412-0

Google Scholar

Crossref

51.

Loiseau

,

J.-C.

, and

S. L.

Brunton

, “

Constrained sparse Galerkin regression

,”

J. Fluid Mech.

838

,

42

–

67

(

2018

).

https://doi.org/10.1017/jfm.2017.823

Google Scholar

Crossref

52.

Sato

,

T.

, and

T.

Taniguchi

, “

Rheology and entanglement structure of well-entangled polymer melts: A slip-link simulation study

,”

Macromolecules

52

,

3951

–

3964

(

2019

).

https://doi.org/10.1021/acs.macromol.9b00314

Google Scholar

Crossref

53.

Miyamoto

,

S.

,

T.

Sato

, and

T.

Taniguchi

, “

Stretch-orientation-induced reduction of friction in well-entangled bidisperse blends: A dual slip-link simulation study

,”

Rheol. Acta

62

,

57

–

70

(

2023

).

https://doi.org/10.1007/s00397-022-01378-5

Google Scholar

Crossref

54.

Sato

,

T.

,

S.

Moghadam

,

G.

Tan

, and

R. G.

Larson

, “

A slip-spring simulation model for predicting linear and nonlinear rheology of entangled wormlike micellar solutions

,”

J. Rheol.

64

,

1045

–

1061

(

2020

).

https://doi.org/10.1122/8.0000062

Google Scholar

Crossref

55.

Sato

,

T.

, and

R. G.

Larson

, “

Nonlinear rheology of entangled wormlike micellar solutions predicted by a micelle-slip-spring model

,”

J. Rheol.

66

,

639

–

656

(

2022

).

https://doi.org/10.1122/8.0000426

Google Scholar

Crossref

56.

Hyun

,

K.

,

M.

Wilhelm

,

C. O.

Klein

,

K. S.

Cho

,

J. G.

Nam

,

K. H.

Ahn

,

S. J.

Lee

,

R. H.

Ewoldt

, and

G. H.

McKinley

, “

A review of nonlinear oscillatory shear tests: Analysis and application of large amplitude oscillatory shear (LAOS)

,”

Prog. Polym. Sci.

36

,

1697

–

1753

(

2011

).

https://doi.org/10.1016/j.progpolymsci.2011.02.002

Google Scholar

Crossref

Rheo-SINDy: Finding a constitutive model from rheological data for complex fluids using sparse identification for nonlinear dynamics

I. INTRODUCTION

II. METHODS

III. CASE STUDIES

A. Upper-convected Maxwell (UCM) model

B. Giesekus model

C. Dumbbell models

1. Hookean dumbbell model

2. FENE-P dumbbell model

3. FENE dumbbell model

IV. RESULTS AND DISCUSSIONS

A. Upper convected Maxwell model

B. Giesekus model

C. Hookean dumbbell model

D. FENE-P dumbbell model

E. FENE dumbbell model

F. Material frame indifference

V. CONCLUDING REMARKS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Rheo-SINDy: Finding a constitutive model from rheological data for complex fluids using sparse identification for nonlinear dynamics

I. INTRODUCTION

II. METHODS

III. CASE STUDIES

A. Upper-convected Maxwell (UCM) model

B. Giesekus model

C. Dumbbell models

1. Hookean dumbbell model

2. FENE-P dumbbell model

3. FENE dumbbell model

IV. RESULTS AND DISCUSSIONS

A. Upper convected Maxwell model

B. Giesekus model

C. Hookean dumbbell model

D. FENE-P dumbbell model

E. FENE dumbbell model

F. Material frame indifference

V. CONCLUDING REMARKS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only