<figure><img src="image_1.png" /></figure>

HAL
open science

# On the reduction of stochastic kinetic theory models of complex fluids

Francisco Chinesta, Amine Ammar, Antonio Falco, Manuel Laso

► To cite this version:

Francisco Chinesta, Amine Ammar, Antonio Falco, Manuel Laso. On the reduction of stochastic kinetic theory models of complex fluids. Modelling and Simulation in Materials Science and Engineering, 2007, 15, pp.639 - 652. 10.1088/0965-0393/15/6/004. hal-01476037

HAL Id: hal-01476037

https://hal.science/hal-01476037

Submitted on 24 Feb 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# On the reduction of stochastic kinetic theory models of complex fluids

F Chinesta¹, A Ammar², A Falco³ and M Laso⁴

¹ Laboratoire de Mécanique des Systèmes et des Procédés, UMR 8106 CNRS-ENSAM-ESEM, 151 Boulevard de l’Hôpital, F-75013 Paris, France

² Laboratoire de Rhéologie UMR 5520 CNRS-INPG-UJF, 1301 rue de la piscine, BP 53 Domaine Universitaire, F-38041 Grenoble Cedex 9, France

³ Universidad CEU Cardenal Herrera, San Bartolomé 55, E-46115 Alfara del Patriarca, Spain

⁴ Laboratory of Non-Metallic Materials, ETSII-UPM, Jose Gutierrez Abascal 2, E-28006 Madrid, Spain

E-mail: francisco.chinesta@paris.ensam.fr, Amine.Ammar@ujf-grenoble.fr, afalco@uch.ceu.es and laso@diquima.upm.es

## Abstract

Kinetic theory models involving the Fokker-Planck equation are usually solved in the framework of stochastic approaches, which allows us to circumvent the difficulties related to the multidimensional character of that equation. In fact, the Fokker-Planck equation governs the evolution of the distribution function that defines the molecular configuration at each point of the physical space and at each time. As the molecular conformation is usually defined by several coordinates, the resulting distribution function will depend on the physical and configuration coordinates and the time. Although different numerical strategies have recently been proposed for solving that equation with efficiency and accuracy (Ammar et al 2006 *J. Non-Newtonian Fluid Mech.* **134** 136–47, Ammar et al 2006 *J. Non-Newtonian Fluid Mech.* **139** 153–76) the stochastic approach is today the most common for solving general kinetic theory models. This paper presents some preliminary results that provide evidence for the potential applicability of model reduction techniques based on the Karhunen-Loève decomposition or on separated representations for reducing the computational efforts related to the solution of such models in the Brownian configuration fields framework.

(Some figures in this article are in colour only in the electronic version)

## 1. Introduction

Many natural and synthetic fluids are viscoelastic materials, in the sense that the stress endured by a macroscopic fluid element depends upon the history of the deformation experienced by that element.
Atomistic modeling is the most detailed level of description that can be applied today in rheological studies, but its simulation requires enormous computer resources, and so they are currently limited to flow geometries of molecular dimensions.

Kinetic theory models provide a coarse-grained description of molecular configurations. They are meant to display in a more or less accurate fashion the important features that govern the flow-induced evolution of configurations.

Micro–macro methods couple the coarse-grained molecular scale of kinetic theory to the macroscopic scale of continuum mechanics (the reader can refer to the review paper [9] and the references therein). This approach is much more demanding in computer resources than more conventional continuum simulations that integrate a constitutive equation to evaluate the viscoelastic contribution of the stress tensor.

Since the early 1990s the field has developed considerably following the introduction of the CONNFFESSIT method by Ottinger and Laso [10].

Kinetic theory provides two basic building blocks: the diffusion or Fokker–Planck equation that governs the evolution of the distribution function (giving the probability distribution of configurations) and an expression relating the viscoelastic stress to the distribution function. The Fokker–Planck equation has the general form

$$ \frac{d\psi}{dt} + \frac{\partial}{\partial \underline{X}} (A\underline{\psi}) = \frac{1}{2} \frac{\partial}{\partial \underline{X}} \frac{\partial}{\partial \underline{X}} : (\underline{D}\underline{\psi}), \quad (1) $$

where vectors are affected by an underline and matrices by a double underline, $d\psi/dt$ is the material derivative, vector $\underline{X}$ defines the coarse-grained configuration and has dimensions $N$. Factor $A$ is an $N$-dimensional vector that defines the drift or deterministic component of the molecular model. Finally $\underline{D}$ is a symmetric, positive definite $N \times N$ matrix that embodies the diffusive or stochastic component of the molecular model. In general both $A$ and $\underline{D}$ (and in consequence the distribution function $\underline{\psi}$) depend on the physical coordinates $\underline{x}$, on the configuration coordinates $\underline{X}$ and on the time $t$.

The second building block of a kinetic theory model is an expression relating the distribution function and the stress. It takes the form

$$ \underline{\tau}_p = \int_C \underline{g}(\underline{X}) \underline{\psi} \, d\underline{X}, \quad (2) $$

where $C$ represents the configuration space and $\underline{g}$ is a model-dependent tensorial function of configuration. In a complex flow, the velocity field is a priori unknown and the stress field are coupled through the conservation laws. In the isothermal and incompressible case the conservation of mass and momentum balance are then expressed (neglecting the body forces) by

$$ \begin{gathered} \mathrm{Div} \underline{v} = 0, \\ \rho \frac{d\underline{v}}{dt} = \mathrm{Div}(-p\underline{I} + \underline{\tau}_p + 2\eta_s \underline{d}), \end{gathered} \quad (3) $$

where $\rho$ is the fluid density, $p$ the pressure and $2\eta_s \underline{d}$ a purely viscous component ($\underline{d}$ being the strain rate tensor). The set of coupled equations (1)–(3), supplemented with suitable initial and boundary conditions in both physical and configuration spaces, is the generic multiscale formulation.

Three basic approaches have been adopted for exploiting the generic multiscale model: the continuum approach; the Fokker–Planck approach and the stochastic approach. This paper focuses on the last approach.

The stochastic approach is based on the mathematical equivalence between the Fokker–Planck equation (1) and the following Ito stochastic differential equation:

$$ d\underline{X} = \underline{A} dt + \underline{B} d\underline{W}, \quad (4) $$
where $\underline{D} = \underline{\underline{B}} \underline{\underline{B}}^T$ and $\underline{W}$ is a Wiener stochastic process of dimension $N$. In a complex flow, the stochastic differential equation (4) applies along individual flow trajectories, the time derivative is thus a material derivation. Instead of solving the deterministic Fokker-Planck equation (1), one solves the associated stochastic differential equation (4) for a large ensemble of realizations of the stochastic process $\underline{X}$ by means of a suitable numerical technique.

The control of the statistical noise is a major issue in stochastic micro-macro simulations based on the stochastic approach. Moreover, to reconstruct the distribution one needs to operate with an extremely large number of particles; however, in general, only the moments of such a distribution are required, which can be computed as an ensemble average using a much more reduced population of particles.

The Brownian configuration fields (BCFs) is a technique proposed in [8] allowing us to reduce the variance as well as to accelerate the numerical simulations. In this paper we present some preliminary results that seem to reveal that BCF models can be reduced by using a Karhunen-Loève decomposition which does not run in the context of the original stochastic approach.

The proposed technique will be illustrated on the kinetic theory model related to a short fiber suspension in a Newtonian fluid. For this reason, we start by introducing such a model in the next section.

### 1.1. Governing equations for a short fiber suspension

In the case of a dilute short fiber suspension, the configuration distribution function (also known as orientation distribution function) gives the probability of finding the fiber in a given direction. Obviously, this function depends on the physical coordinates (space and time) as well as on the configuration coordinates, which taking into account the rigid character of the fibers, are defined on the surface of the unit sphere. Thus, we can write $\psi(\underline{x}, t, \underline{p})$, where $\underline{x}$ defines the position of the fiber center of mass, $t$ the time and $\underline{p}$ the unit vector defining the fiber orientation. The evolution of the distribution function is given by the Fokker-Planck equation:

$$ \frac{d\psi}{dt} = -\frac{\partial}{\partial \underline{p}}(\psi \dot{\underline{p}}) + \frac{\partial}{\partial \underline{p}}\left(D_r \frac{\partial \psi}{\partial \underline{p}}\right), \qquad (5) $$

where $d/dt$ represents the material derivative, $D_r$ is a diffusion coefficient and $\dot{\underline{p}}$ is the fiber rotation velocity. The orientation distribution function must verify the normality condition:

$$ \oint \psi(\underline{p}) d\underline{p} = 1. \qquad (6) $$

When the fibers are assumed ellipsoidal and when the suspension is dilute enough, the rotation velocity can be obtained from the Jeffery's equation

$$ \dot{\underline{p}} = \underline{\underline{\Omega}} \underline{\underline{p}} + k \underline{\underline{D}} \underline{\underline{p}} - k(\underline{\underline{p}}^T \underline{\underline{D}} \underline{\underline{p}}) \underline{\underline{p}}, \qquad (7) $$

where $\underline{\underline{\Omega}}$ and $\underline{\underline{D}}$ are the vorticity and the strain rate tensors, respectively, associated with the fluid flow undisturbed by the presence of the fibers, and $k$ is a scalar which depends on the fiber aspect ratio $\lambda$ (ratio between the fiber length and the fiber diameter)

$$ k = \frac{\lambda^2 - 1}{\lambda^2 + 1}. \qquad (8) $$

In a former work [4] the discretization of the steady-state Fokker-Planck equation was carried out in steady recirculating flows using a particle technique, where the diffusion term was modeled from random motions. It was pointed out that the number of fibers required in this stochastic simulation to describe the fiber distribution increases significantly with the diffusion coefficient $D_r$.
## 2. Reduced order modeling

The model reduction technique that we propose in this work is based on the use of the Karhunen–Loève decomposition, which we summarize in the next section.

### 2.1. The Karhunen–Loève decomposition

We assume that the evolution of a certain field $u(\underline{x}, t)$ is known (being its evolution governed by a PDE). In practical applications, this field is expressed in a discrete form, that is, it is known at the $N_n$ nodes of a spatial mesh and at some times $u(\underline{x}_i, t^n) \equiv u_i^n, \forall n \in [1, \dots, P]$ and $\forall i \in [1, \dots, N_n]$. The main idea of the Karhunen–Loève (KL) decomposition is how to obtain the most typical or characteristic structure $\phi(\underline{x})$ among these $u_i^n(\underline{x})$, $\forall n$. This is equivalent to obtaining a function $\phi(\underline{x})$ that maximizes $\lambda$ defined by

$$ \lambda = \frac{\sum_{n=1}^{n=P} \left[ \sum_{i=1}^{i=N_n} \phi(\underline{x}_i) u_i^n(\underline{x}_i) \right]^2}{\sum_{i=1}^{i=N} (\phi(\underline{x}_i))^2}. \quad (9) $$

By applying standard variational calculus, the maximization of equation (9) leads to

$$ \sum_{n=1}^{n=P} \left[ \left( \sum_{i=1}^{i=N_n} \tilde{\phi}(\underline{x}_i) u_i^n(\underline{x}_i) \right) \left( \sum_{j=1}^{j=N_n} \phi(\underline{x}_j) u_j^n(\underline{x}_j) \right) \right] = \lambda \sum_{i=1}^{i=N_n} \tilde{\phi}(\underline{x}_i) \phi(\underline{x}_i); \quad \forall \tilde{\phi} \quad (10) $$

which can be rewritten in the form

$$ \sum_{i=1}^{i=N_n} \left\{ \sum_{j=1}^{j=N_n} \left[ \sum_{n=1}^{n=P} u^n(\underline{x}_i) u^n(\underline{x}_j) \phi(\underline{x}_j) \right] \tilde{\phi}(\underline{x}_i) \right\} = \lambda \sum_{i=1}^{i=N_n} \tilde{\phi}(\underline{x}_i) \phi(\underline{x}_i); \quad \forall \tilde{\phi} \quad (11) $$

Defining the vector $\underline{\phi}$ such that its $i$-component is $\phi(\underline{x}_i)$, equation (11) takes the following matrix form

$$ \underline{\tilde{\phi}}^{\mathrm{T}} \underline{C} \underline{\phi} = \lambda \underline{\tilde{\phi}}^{\mathrm{T}} \underline{\phi}; \quad \forall \underline{\tilde{\phi}} \Rightarrow \underline{C} \underline{\phi} = \lambda \underline{\phi}, \quad (12) $$

where the two points correlation matrix is given by

$$ C_{ij} = \sum_{n=1}^{n=P} u^n(\underline{x}_i) u^n(\underline{x}_j) \Leftrightarrow \underline{C} = \sum_{n=1}^{n=P} \underline{u}^n(\underline{u}^n)^{\mathrm{T}}, \quad (13) $$

which is symmetric and positive definite. If we define the matrix $\underline{Q}$ containing the discrete field history:

$$ \underline{Q} = \begin{pmatrix} u_1^1 & u_1^2 & \cdots & u_1^P \\ u_2^1 & u_2^2 & \cdots & u_2^P \\ \vdots & \vdots & \ddots & \vdots \\ u_{N_n}^1 & u_{N_n}^2 & \cdots & u_{N_n}^P \end{pmatrix} \quad (14) $$

then, matrix $\underline{C}$ in equation (12) gives

$$ \underline{C} = \underline{Q} \underline{Q}^{\mathrm{T}}. \quad (15) $$

Thus, the functions defining the most characteristic structure of $u_i^n(\underline{x})$ are the eigenfunctions $\phi_k(\underline{x})$, whose discrete expression is $\underline{\phi}_k$, associated with the highest eigenvalues.
## 2.2. A posteriori reduced modeling

If some direct simulations have been carried out, we can determine $u(\underline{x}_i, t^n) \equiv u_i^n$, $\forall i \in [1, \dots, N_n]$, $\forall n \in [1, \dots, P]$, and from that information the $r$ eigenvectors $\phi_k^T = [\phi_k(\underline{x}_1), \dots, \phi_k(\underline{x}_{N_n})]$, $\forall k \in [1, \dots, r]$ (with $r \ll N_n$) related to the $r$-highest eigenvalues $\lambda_1, \dots, \lambda_r$ (that are assumed ordered, being $\lambda_1$ the highest eigenvalue). These eigenvalues verify $\lambda_k > \epsilon\lambda_1$, $\forall k \in [1, \dots, r]$ (with $\epsilon$ a small enough value that in our simulations is set to $10^{-8}$) and $\lambda_k < \epsilon\lambda_1$, $\forall k \in [r+1, \dots, N_n]$.

Now, we can try to use these $r$ eigenfunctions for approximating the solution of a problem slightly different from the one that has served to define $u(\underline{x}_i, t^n)$. For this purpose we need to define the matrix

$$ \underline{B} = \begin{pmatrix} \phi_1(\underline{x}_1) & \phi_2(\underline{x}_1) & \cdots & \phi_r(\underline{x}_1) \\ \phi_1(\underline{x}_2) & \phi_2(\underline{x}_2) & \cdots & \phi_r(\underline{x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_1(\underline{x}_{N_n}) & \phi_2(\underline{x}_{N_n}) & \cdots & \phi_r(\underline{x}_{N_n}) \end{pmatrix}. \quad (16) $$

Now, we consider the linear system of equations resulting from the discretization of a partial differential equation (PDE) in the form

$$ \underline{K} \underline{U}^n = \underline{F}^{n-1}. \quad (17) $$

In the case of transient problems $\underline{F}^{n-1}$ contains the contribution of the solution at the previous time step.

Then, the unknown vector containing the nodal degrees of freedom can be expressed as

$$ \underline{U}^n = \sum_{i=1}^{i=r} \underline{\phi}_i \underline{\xi}_i^n = \underline{B} \underline{\xi}^n, \quad (18) $$

which implies

$$ \underline{K} \underline{U}^n = \underline{F}^{n-1} \Rightarrow \underline{K} \underline{B} \underline{\xi}^n = \underline{F}^{n-1} \quad (19) $$

and multiplying both terms by $\underline{B}^T$ it gives

$$ \underline{B}^T \underline{K} \underline{B} \underline{\xi}^n = \underline{B}^T \underline{F}^{n-1}, \quad (20) $$

which proves that the resulting linear system has a small size, i.e. the dimensions of $\underline{B}^T \underline{K} \underline{B}$ are $r \times r$, with $r \ll N_n$, and the dimensions of both $\underline{x}_i$ and $\underline{B}^T \underline{F}$ are $r \times 1$.

**Remark 2.1.** Equation (20) can also be derived introducing the approximation (18) into the PDE Galerkin form.

See [1] or [11] for more details on the application on model reduction for simulating kinetic theory models based on the use of the Fokker-Planck formalism.

## 3. Reduced Brownian configuration fields (R-BCF)

The BCF technique is based on substituting the solution of equation (4) along individual trajectories by the solution of the evolution of several fields (the so-called Brownian configuration fields). From now on we consider the kinetic theory model related to the flow of a short fiber suspension, with both the flow and the fiber orientation assumed in a 2D physical space (the extension to 3D is straightforward). As we are interested in the orientation solution,
we consider that the flow kinematics and the fiber orientation are uncoupled, and therefore we
are computing the orientation solution for a given velocity field.

We assume a simple shear flow characterized by the following kinematics

$$
\underline{v} = \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} \dot{\gamma} y \\ 0 \end{pmatrix}. \qquad (21)
$$

The fiber orientation can be described from the angle $\varphi$:

$$
\underline{p} = \begin{pmatrix} p_x \\ p_y \end{pmatrix} = \begin{pmatrix} \cos \varphi \\ \sin \varphi \end{pmatrix}. \qquad (22)
$$

Thus, the Fokker-Planck equation (equation (5)) reduces to

$$
\frac{\partial \psi}{\partial t} + u \frac{\partial \psi}{\partial x} = - \frac{\partial}{\partial \varphi} (\psi \dot{\varphi}) + D_r \frac{\partial^2 \psi}{\partial \varphi^2}, \quad (23)
$$

where the diffusion coefficient $D_r$ is assumed scalar and constant.

The fibers are assumed with an infinite aspect ratio ($k = 1$ in equation (7)), leading to

$$
\dot{\varphi} = \frac{d\varphi}{dt} = -\sin^2 \varphi, \qquad (24)
$$

which implies that the flow tends to align the fiber in the flow direction.

The domain in which the Fokker-Planck (FP) equation is defined is defined by Ω = ]0, L] × [−H, H], and the time interval ]0, T]. Due to the advective character of the FP equation in the physical space only a boundary condition is required on the inflow boundary x = 0, which we represent by ψ<sub>inf</sub>(y, t). The initial distribution is given by ψ(x, t = 0) = ψ<sup>0</sup>.

Again, due to the advective character of the FP equation in the physical space, the solution on a flow trajectory y = cte does not depend on the solution on the neighbor trajectories.
Thus, we will restrict our analysis to the one-dimensional trajectory defined by y = 1.

The simplest (explicit and first order) stochastic simulation consists of the following.

1. Define a number of fibers $\varphi_i$ initially distributed in the physical domain $]0, L]$ representing the initial fiber distribution $\psi^0$.

2. Update the position of each fiber according to the flow kinematics (in our case defined on $y = 1$ with $\dot{\gamma} = 1$):

$$
\begin{cases}
x_i^{n+1} = x_i^n + u(x_i^n, y_i^n, t^n) \Delta t = x_i^n + \Delta t, \\
y_i^{n+1} = 1
\end{cases}
\quad \forall i, \forall n \ge 1 \quad (25)
$$

3. Update the orientation of each fiber according to

$$
\varphi_i^{n+1} = \varphi_i^n + \dot{\varphi}(\varphi_i^n, t^n) \Delta t + \Delta W_i^{n,n+1}, \quad \forall i, \quad (26)
$$

where $\dot{\varphi}(\varphi_i^n, t^n) = -\sin^2(\varphi_i^n)$ and $\Delta W$ is a random number with zero mean and a variance of $2D_r \Delta t$.

4. To avoid the loss of fibers we introduce new fibers in the domain through the inflow boundary $x=0$ whose orientation represents the boundary condition $\psi_{\inf}$, in order to keep constant the number of particles into the domain.

Despite the consideration of a first order explicit scheme for the stochastic equation
integration, higher order integration schemes are available.

We prove later that the extraction of representative modes from this kind of analysis is not
evident, a fact that motivates the consideration of the BCF framework.

Because the random term in the stochastic equation is uncorrelated in space, it is possible
to define stochastic fields whose evolutions can be solved by using standard techniques for
PDEs, making use of fixed or moving meshes. In what follows we consider a fixed mesh on
the physical domain on which the evolution of different fields $\varphi_i(x, t)$ is computed. For this purpose we simply substitute the resolution of equations (25) and (26) by the solution of

$$ \varphi_i^{n+1}(x) = \varphi_i^n(x) - u(x) \frac{\partial \varphi_i^{n+\theta}}{\partial x} \Delta t + \dot{\varphi}(\varphi_i^{n+\theta}(x)) \Delta t + \Delta W_i^{n,n+1}, \quad \forall i, \quad (27) $$

where $\partial u/\partial x = 0$.

Different techniques have been used in the literature for solving this equation, all of them stabilizing its advection character: SUPG, discontinuous-Galerkin, etc. In the case of considering an explicit strategy ($\theta = 0$ in equation (27)) no linear system must be solved but extremely small time steps are needed. In the other case, if one considers implicit ($\theta = 1$) or semi-implicit ($0 < \theta < 1$) schemes, the solution of a linear system (whose size corresponds to the number of nodes used in the space discretization, $N_n$) is required at each time step.

In the semi-implicit case, where the flow-induced orientation term is considered at the previous time step, we obtain

$$ \underline{K} \varphi_i^{n+1} = \underline{F}^n (\Delta W_i^{n,n+1}). \quad (28) $$

Obviously, the fully implicit case results in a non-linear system that must be solved using an appropriate linearization schema.

One possibility for reducing the computing time consists of solving the previous model for one field, and then extracting from its time evolution the $r$ characteristic modes via the application of the Karhunen-Loève decomposition, defining the reduced approximation basis. Now, the remaining fields could be computed after projection on the reduced approximation basis. Thus, even if one is using an implicit (or semi-implicit) strategy the size of the linear systems involved is $r \times r$ instead of $N_n \times N_n$.

## 3.1. Numerical results

### 3.1.1. Evaluating the applicability of reduction techniques in the standard stochastic approach.
First, we consider a population consisting of $N$ fibers located on the streamline $y = 0$, isotropically distributed. These fibers are subjected to a unit shear flow ($\dot{y} = 1$) but they do not move along the streamline because the velocity field vanishes on this streamline.

The initial fibers orientation is defined by

$$ \varphi_i^0 = 0. \quad (29) $$

The orientation updating is defined at each time step according to

$$ \varphi_i^{n+1} = \varphi_i^n + \dot{\varphi}(\varphi_i^n) \Delta t + \Delta W_i^{n,n+1}, \quad \forall i. \quad (30) $$

Now, we define the matrix $\underline{Q}$

$$ \underline{Q} = \begin{pmatrix} \varphi_1^1 & \varphi_1^2 & \cdots & \varphi_1^P \\ \varphi_2^1 & \varphi_2^2 & \cdots & \varphi_2^P \\ \vdots & \vdots & \ddots & \vdots \\ \varphi_N^1 & \varphi_N^2 & \cdots & \varphi_N^P \end{pmatrix}, \quad (31) $$

which defines the eigenvalue problem

$$ \underline{Q} \underline{Q}^T \underline{\phi} = \lambda \underline{\phi}. \quad (32) $$

The main difficulty related to the use of these computed modes is that they are associated with a specific particles position, the one that served to define $\underline{Q}$. As soon as the particle trajectories differ from the ones that served to compute the characteristic functions,
the reduced approximation is no more valid. To illustrate this fact, we consider the reduced approximation basis composed of the $r$ eigenvectors related to the highest eigenvalues (the $r$ verifying $\lambda_k > 10^{-8}\lambda_1$, being $\lambda_1$ the highest one). These eigenvectors define the matrix $\underline{B}$ according to equation (16).

Now, equation (30) is written in the matrix form

$$ \underline{\varphi}^{n+1} = \underline{\varphi}^n + \dot{\underline{\varphi}}(\underline{\varphi}^n) \Delta t + \underline{\Delta W}^{n,n+1} \quad (33) $$

that after projection onto the reduced approximation basis results

$$ \underline{B}^T \underline{B} \underline{\xi}^{n+1} = \underline{B}^T \underline{B} \underline{\xi}^n + \underline{B}^T \dot{\underline{\varphi}}(\underline{B} \underline{\xi}^n) \Delta t + \underline{B}^T \underline{\Delta W}^{n,n+1}. \quad (34) $$

In general the initial condition is added to the reduced approximation basis:

$$ \underline{B} \leftarrow [\underline{\varphi}^0 \underline{B}], \quad (35) $$

which allows us to properly define the initial condition associated with the reduced unknown vector $\underline{\xi}^0$:

$$ (\underline{\xi}^0)^T = (1, 0, \dots, 0). \quad (36) $$

The integration of equation (34) from the initial condition (36) allows us to compute the particles orientation at any time:

$$ \underline{\varphi}^n = \underline{B} \underline{\xi}^n \quad (37) $$

and the associated moments of the orientation distribution. In particular the second order moment results

$$ \underline{a}^n = \int_0^{2\pi} \underline{p} \otimes \underline{p} \, \underline{\psi}^n(\underline{\varphi}) \, d\underline{\varphi}, \quad (38) $$

where $\underline{p}$ is defined in equation (22). When one uses the stochastic approach the distribution function is defined from the Dirac masses according to

$$ \underline{\psi}^n(\underline{\varphi}) = \sum_{i=1}^{i=N} \frac{1}{N} \delta(\underline{\varphi} - \underline{\varphi}_i^n), \quad (39) $$

which implies

$$ \underline{a}^n = \frac{1}{N} \sum_{i=1}^{i=N} \begin{pmatrix} \cos^2(\varphi_i^n) & \sin(\varphi_i^n) \cdot \cos(\varphi_i^n) \\ \sin(\varphi_i^n) \cdot \cos(\varphi_i^n) & \sin^2(\varphi_i^n) \end{pmatrix}. \quad (40) $$

Now, we can compare the evolution of the first component of the second order orientation tensor $\underline{a}_{11}$, computed from the particles position determined by using both a standard stochastic procedure (equation (30)) and the reduced modeling defined in equation (34). Figure 1 depicts both solutions from which we can conclude the lack of accuracy related to the reduced modeling. In this simulation we consider $D_r = 0.05$, $N = 300$, $k = 1$, $\dot{\gamma} = 1$, $\Delta t = 0.1$ and $T = 30$. All the fibers were assumed initially aligned in the $x$-direction according to equation (29). The stochastic solution, even if there is a significant noise (which can be reduced by increasing the number of particles), approaches the steady-state exact solution, whereas the computed in the reduced approximation basis (in this case composed of 14 eigenfunctions) evolves towards an isotropic distribution.

In conclusion, the direct reduction of stochastic simulations by applying the Karhunen-Loève decomposition does not run despite the simplicity of the flows considered in the previous analysis. In our opinion the reason for this behavior is that the reduced basis is associated with the particles trajectories involved in the Karhunen-Loève decomposition. Thus, when this reduced basis is applied for describing other particles trajectories significant deviations are expected and, in fact, found.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Comparison between reduced order solution (blue stars) and the standard stochastic one (red stars).</figcaption></figure>

### 3.1.2. Numerical results on the reduced BCF.
Now, we solve the same problem on the streamline related to $y = 1$ using the reduced BCF strategy previously described. On this streamline the kinematics is defined by $u = 1$ and $v = 0$. At present we assume that all the fibers located on the inflow boundary $x = 0$ are aligned in the y-direction, i.e. $\varphi(x = 0, t) = \pi/2$. The same fiber orientation was assumed in the flow domain at the initial time $\varphi(x, t = 0) = \pi/2$.

We compute the evolution of a BCF $\varphi_1(x, t)$ from its initial value $\varphi_1(x, t = 0) = \pi/2$ until the maximum simulation time $T = \pi/2$, being the physical domain length $L = 3$. The diffusion coefficient was set to $D_r = 0.05$ and an infinite fiber aspect ratio.

The number of Brownian configurations fields was set to 200 and 50 nodes were uniformly distributed on the physical domain ($N_n = 50$). A semi-implicit backward finite difference scheme was employed to discretize the transport equation related to the first BCF $\varphi_1(x, t)$, being the time step $\Delta t = T/300$. From this evolution, matrix $\underline{Q}$ was computed, allowing us to extract the characteristic solution information, which in this case consists of five eigenfunctions. Now, the evolution of the remaining BCF $(\varphi_2(x, t) \cdots \varphi_{200}(x, t))$ was computed using the reduced approximation basis consisting of the five eigenfunctions extracted from the complete evolution analysis of the first BCF $\varphi_1(x, t)$. Thus, the size of the linear systems solved at each time step, for all the BCF (except the first one) was $5 \times 5$, allowing a significant CPU time reduction.

Figure 2 depicts the different components of the second order orientation tensor ($a_{22}$ in red, $a_{11}$ in blue and $a_{12}$ in green). The continuous line represents the solution when all the BCFs are computed using the global approximation basis (non-reduced model):

$$ \underline{K} \varphi_i^{n+1} = \underline{F}^n (\Delta W_i^{n,n+1}), \quad \forall i. \quad (41) $$

The evolution of the first BCF served to define the reduced approximation basis ($\underline{B}$) via the Karhunen-Loève decomposition. When this basis is used to span the remaining BCF evolutions

$$ \underline{B}^T \underline{K} \underline{B} \underline{\xi}_i^{n+1} = \underline{B}^T \underline{F}^n (\Delta W_i^{n,n+1}), \quad \forall i > 1, \quad (42) $$
<figure><img src="image_3.png" /><figcaption>Figure 2. Comparison between reduced and non-reduced BCF solutions.</figcaption></figure>

the resulting orientation tensor components are illustrated using the star symbols in figure 2. We notice a slight deviation between both solutions, the non-reduced one and the one just computed using the reduced approximation basis.

The aforementioned deviation was expected because (i) the reduced simulation defines the solution on the approximation basis spanned by the five eigenfunctions extracted from the solution of the first BCF and (ii) the Brownian terms $\Delta W_i^{n,n+1}$ in equations (41) and (42) are different. One could expect that the random effects being the same, i.e. $\Delta W_i^{n,n+1}$ in equations (41) and (42) being the same, the accuracy should increase, because in this case the accuracy only depends on the number of significant eigenfunctions retained from the Karhunen-Loève decomposition of $\varphi_1(x, t)$. In order to prove this, we compute again the evolution of the different orientation tensor components using the non-reduced approach (equation (41)) but considering the random term $W_i^{n,n+1}$ used to integrate equation (42). It can be noticed in figure 2 (discontinuous curves) that in this case both results are very close, proving that a very reduced number of approximation functions are enough to accurately represent the field evolution.

The accuracy can be improved by increasing the number of nodes used in the space discretization, the number of BCF or by reducing the time step. In any case the results are very stable, and similar results were obtained by running the simulation code several times, with deviations that rarely exceed some per cent. As the computed solutions are very accurate there is no necessity of enrichment of the reduced approximation basis. In any case such an adaptation could be carried out using some Krylov subspaces as described in [11] or [1].

## 4. Separated representation of the reduced Brownian configurations fields

In this section we explore the application of a separated representation and the associated tensor product approximation basis for solving the transport equation governing the evolution of the BCFs. Thus, coming back to that equation:

$$ \varphi_i^{n+1}(x) = \varphi_i^n(x) - u(x) \frac{\partial \varphi_i^{n+\theta}}{\partial x} \Delta t + \dot{\varphi}(\varphi_i^{n+\theta}(x)) \Delta t + \Delta W_i^{n,n+1}, \quad \forall i, \quad (43) $$
we can write

$$ \frac{\partial \varphi_i(x, t)}{\partial t} = -u(x) \frac{\partial \varphi_i(x, t)}{\partial x} + (\dot{\varphi}(\varphi_i(x)) + H_i(t)), \quad \forall i, \quad (44) $$

where $H_i(t)$ is a piecewise constant function defined as

$$ H_i(t) = \begin{cases} \frac{\Delta W_i^{0,1}}{\Delta t} & 0 < t < \Delta t \\ \vdots & \vdots \\ \frac{\Delta W_i^{n,n+1}}{\Delta t} & n \Delta t < t < (n+1) \Delta t \\ \vdots & \vdots \\ \frac{\Delta W_i^{P-1,P}}{\Delta t} & T - \Delta t < t < T, \end{cases} \quad (45) $$

where $W_i^{n,n+1}$ is a random number with zero mean and variance $2D_r\Delta t$, $\Delta t$ being the time step that will be considered for the time discretization of equation (44). In any case, with $H_i(t)$ computed according to the previous expression, it becomes a deterministic time function which affects the evolution of the configuration field under consideration.

Now, the separated representation and the associated tensor product approximation basis could be applied to perform a space-time simultaneous resolution, according to the procedures proposed in [2,3] which are based on the following functional approximation:

$$ \varphi_i(x, t) = \sum_j \alpha_j^i \times F_j^i(x) \times G_j^i(t). \quad (46) $$

The construction of this solution requires an iteration scheme involving a projection and an enrichment step at each iteration. The simplest reduction strategy lies in the construction of the tensor product approximation basis, defined by functions $F_j^1(x)$ and $G_j^1(t)$ by solving the evolution problem related to the first configuration field, and then looking for the solution of the remaining configuration fields by a simple projection onto the basis $F_j^1(x)$ and $G_j^1(t)$ which allows us to compute the coefficients $\alpha_j^i$, $\forall i > 1$. Note that due to the non-linear character of equation (44) this projection stage requires an appropriate iteration scheme.

Figure 3 illustrates the results computed by using the reduced approximation basis obtained from the separated representation of the first BCF. In this figure we depict the functions $F_j^1$, $G_j^1$, the components of the second order orientation tensor at the final time $t = T = \pi/2$ and finally the orientation ellipsoids (the axes length represents the intensity of the fiber orientation along the axes direction) in the space-time domain.

Figure 4 compares the just computed components of the second order orientation tensor and the ones computed by solving the separated representation of the 200 BCF.

The computed solution when one solves the reduced representation of all the BCF is very close to the one obtained from the reduced modeling based on the use of the Karhunen-Loève decomposition. For comparison purposes we represent in figure 5 the just computed solution and the one computed using the Karhunen-Loève reduction technique (also depicted in figure 2).

From these preliminary results we notice that the separated representation simulations using the reduced approximation basis computed from a single BCF exhibit lower accuracy than the ones based on the use of the Karhunen-Loève decomposition previously analyzed. A more in-depth analysis on the construction of optimal reduced basis constitutes work in progress.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Reduced BCF using a separated representation: (top-left) functions $F_j^1(x)$; (top-right) functions $G_j^1(t)$, (bottom-left) components of the second order orientation tensor at $t = T = \pi/2$ and (bottom-right) space-time representation of the orientation ellipsoids.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Separated representation of BCF: (left) components of the second order orientation tensor at $t = T = \pi/2$ computed using the reduced approximation basis related to the first BCF and (right) results obtained from a separated representation of the 200 BCF.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Reduced BCF using: (left) a Karhunen–Loève reduction technique and (right) a separated representation of all the 200 BCF.</figcaption></figure>

**5. Conclusions**

In this work we have presented some preliminary results revealing the viability of a model reduction in the context of stochastic simulations. The proposed model reduction is based on a direct numerical simulation of a BCF (or a reduced number of them), from which the more relevant information on the solution can be extracted by applying the Karhunen–Loève decomposition. After that, the evolution of the other BCFs could be computed using the reduced approximation basis just computed, with the associated benefits in the CPU time reduction.

This work opens numerous perspectives. One of them consists of the application of a separated representation and the associated tensor product approximation basis for solving the transport equation governing the evolution of the BCFs, as illustrated in the last section.

At present, due to the low dimensional configuration space as well as to the very simple physical domain considered, no conclusions can be extracted concerning the computing cost and the convergence rate. This work constitutes a first attempt on the reduction of stochastic approaches, and it only proves that a certain reduction can be carried out (in the framework of both the Karhunen–Loève decomposition and the tensor product approximation basis).

**References**

[1] Ammar A, Ryckelynck D, Chinesta F and Keunings R 2006 On the reduction of kinetic theory models related to finitely extensible dumbbells. *J. Non-Newtonian Fluid Mech.* **134** 136–47

[2] Ammar A, Mokdad B, Chinesta F and Keunings R 2006 A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids *J. Non-Newtonian Fluid Mech.* **139** 153–76

[3] Ammar A, Mokdad B and Chinesta F 2007 A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids: II. Transient simulation using space-time separated representations. *J. Non-Newtonian Fluid Mech.* **144** 98–121

[4] Chinesta F, Chaidron G and Poitou A 2003 On the solution of the Fokker–Planck equations in steady recirculating flows involving short fiber suspensions *J. Non-Newtonian Fluid Mech.* **113** 97–125

[5] De Gennes P G 1971 Reptation of a polymer chain in the presence of fixed obstacles. *J. Chem. Phys.* **55** 572–9

[6] Doi M and Edwards S F 1978 *The Theory of Polymer Dynamics* (Oxford: Clarendon)

[7] Fang J, Kroger M and Ottinger H C 2000 A thermodynamically admissible reptation model for fast flows of entangled polymers. II. Model predictions for shear and extensional flows *J. Rheol.* **40** 1293–318
[8] Hulsen M A, van Heel A P G and van der Brule B H A A 1997 Simulation of viscoelastic flows using brownian configuration fields *J. Non-Newtonian Fluid Mech.* **70** 79–101

[9] Keunings R 2004 Micro–macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory *Rheology Reviews* ed D M Binding and K Walters (British Society of Rheology) pp 67–98

[10] Öttinger H C and Laso M 1992 Smart polymers in finite element calculation *Int Congr. on Rheology (Brussels, Belguim)*

[11] Ryckelynck D, Chinesta F, Cueto E and Ammar A 2006 On the *a priori* model reduction: overview and recent developments *Arch. Comput. Methods Eng.* **13/1** 91–128