# On the Reduction of Kinetic Theory Models Related to Finitely Extensible Dumbbells

Amine Ammar, David Ryckelynck, Francisco Chinesta, Roland Keunings

► To cite this version:

Amine Ammar, David Ryckelynck, Francisco Chinesta, Roland Keunings. On the Reduction of Kinetic Theory Models Related to Finitely Extensible Dumbbells. Journal of Non-Newtonian Fluid Mechanics, 2006, 134 (1-3), pp.136-147. 10.1016/j.jnnfm.2006.01.007. hal-01007150

HAL Id: hal-01007150

https://hal.science/hal-01007150

Submitted on 2 Mar 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
On the reduction of kinetic theory models related to finitely extensible dumbbells

A. Ammar<sup>a</sup>, D. Ryckelynck<sup>b</sup>, F. Chinesta<sup>b</sup>, R. Keunings<sup>c</sup>

<sup>a</sup> Laboratoire de Rhéologie, INPG, UJF, CNRS (UMR 5520), 1301 rue de la piscine, BP 53 Domaine Universitaire, F-38041 Grenoble Cedex 9, France

<sup>b</sup> 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

<sup>c</sup> CESAME, Université Catholique de Louvain, Bat. Euler, Av. Georges Lemaitre 4, B-1348 Louvain-la-Neuve, Belgium

Stochastic simulation for finitely extensible non-linear elastic (FENE) dumbbells has been successfully applied (see the review paper of Keunings [R. Keunings, Micro–macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory, in: D.M. Binding, K. Walters (Eds.), Rheology Reviews, British Society of Rheology, 2004, pp. 67–98] and the references therein). The main difficulty in these simulations is related to the high number of realizations required for describing accurately the microstructural state due to Brownian effects.

The discretisation of the Fokker–Planck equation with a mesh support (finite elements, finite differences, finite volumes, spectral techniques, ...) allows to go beyond the difficulty related to Brownian effects. However, kinetic theory models involve physical and conformation spaces. Thus, the molecular distribution depends on time, space as well as on the molecular orientation and extension (conformation coordinates). In this form the resulting Fokker–Planck equation is defined in a space of dimension 7.

In the reduction technique proposed in this paper, a reduced approximation basis is constructed. The new shape functions are defined in the whole domain in an appropriate manner. Thus, the number of degrees of freedom involved in the solution of the Fokker–Planck equation is significantly reduced. The construction of those new approximation functions is done with an ‘a priori’ approach, which combines a basis reduction (using the Karhunen–Loève decomposition) with a basis enrichment based on the use of some Krylov subspaces. This numerical technique is applied for solving the FENE model of viscoelastic flows.

**Keywords:** Kinetic theory; Non-Newtonian fluids; Numerical modelling; Model reduction; Karhunen–Loève decomposition; Krylov’s subspaces; FENE model; Fokker–Planck equation

# 1. Introduction

In the last two decades, great progress has been attained in the numerical simulation of complex fluid flows such as polymer flow for example. The new tendency to imply microscopic description of the molecule behavior allows to make a fine description of the physical state of the polymer. Although the development of computer performances allows to improve the numerical description, finer descriptions require new numerical schemes in order to be efficient. The purpose of this work is to present a robust and simplified model to carry

out micro-simulations of molecular behavior within the kinetic theory framework. This framework is useful when it allows to obtain constitutive equation for polymer. But in most cases the explicit derivation of a constitutive equation implies the introduction of some closure approximations. The impact of these approximations can be significant [7].

The present paper focuses on the kinetic theory model related to a dilute polymer model known as finitely extensible non-linear elastic dumbbell (FENE) for which a constitutive equation can- not be derived explicitly. The polymer solution is described as a suspension of non-interacting dumbbells in a Newtonian solvent. The Fokker–Planck formalism is used to describe the kinetic theory model. It allows to describe the evolution of the configu- ration distribution function, which represents the probability of finding a molecule in a particular configuration. Obviously this function depends on time, space and configuration coordinates. The first works concerning its resolution concern the Brownian
dynamics approach used in the CONNFESSIT approach [12]. This approach was considered in [15] for treating MBS (multi-bead-spring) models. A similar technique was proposed in [16] in the context of MBS kinetic models, which introduces a change of variable and uses a Monte-Carlo technique for treating the diffusion term. In these techniques a high number of particles are introduced in the stochastic simulation to account for Brownian effects. A multi-scale approach using deterministic particles for treating the advection and a different set of particles to account for diffusion effects, which leads in fact to a multi-scale approach, was considered in [6]. The same idea was used in the case of short fibre suspensions flows in [5]. In that work, the discretisation of the advection dominated Fokker-Planck equation governing the fibre orientation was carried out using a particle technique, where the diffusion term was modelled from random walks. It was pointed out that the number of fibres required in this stochastic simulation to describe the fibre distribution increases significantly with the diffusion coefficient. Thus, it was argued that for practical applications the use of the particle method in the framework of a stochastic simulation is restricted to very slight diffusion effects. Another deterministic particle approach, very close to that proposed in [3], was analyzed in [1] using smooth particles, but it was noticed that the impact of smoothing on the solution can be significant. Moreover, the fact that the Fokker-Planck equation is defined in a multi-dimensional space induces the necessity of using an extremely large number of particles, with the associate unfavorable incidence on the methods efficiency. The combination of Brownian dynamics with a macroscopic treatment of motion equations is at the basis of the micro-macro approaches deeply reviewed in [8].

Thus, when diffusion effects vary in a large interval, continuous approximations using a fixed or moving mesh seem to be suitable. In this case accurate stabilizations are required for dealing with small diffusion effects. Some attempts at solving the Fokker-Planck equation using a fixed mesh discretisation exist [9,4]. The main difficulties in this approach are related to the multidimensional character of the problem. Thus, the linear systems obtained after usual implicit or semi-implicit space-time discretisations are extremely large for a practical inversion. On the other hand, explicit discretisations, which do not require matrix inversions, have the constraint of too small time steps.

The purpose of this work is to propose an efficient and accurate discretisation technique able to solve the multidimensional Fokker-Planck equations with a substantial reduction in the number of degrees of freedom involved. This technique operates by extracting automatically, and in a way completely transparent for the user, the most relevant information of the unknown solution to construct the functional approximation from the information just extracted. The new shape functions are defined in the whole domain in an appropriate manner (the optimal functions for describing the model solution). Thus, the number of degrees of freedom involved in the solution of the Fokker-Planck equation is significantly reduced. The construction of those new approximation functions is done with an 'a priori' approach, which combines a basis reduction (using the Karhunen-Loève

decomposition) with a basis enrichment based on the use of some Krylov subspaces [13,14].

The plan of the paper is as follows. The next section introduces the mechanical model related to the finitely extensible dumbbells model. The dimensionless problem will be discretised in Section 3 using the finite element method on a fixed mesh support associated with the configuration space. Section 4 presents a general overview of the "a priori" model reduction technique described in detail in our former works [13,14]. The continuous description is followed by the discrete counterpart that will be used to reduce the FENE numerical model. Finally, Section 5 analyzes the performances of the proposed technique for solving the 1D, 2D and 3D FENE models in start-up of elongation flow.

## 2. Mechanical model using a FENE representation

The dumbbell model consists of two beads connected by a spring connector. The bead serves as an interaction point with the solvent and the spring contains the local stiffness depending on local stretching (see [2] for more detail).

The dynamic of the chain is governed by viscous, Brownian, and connector forces. If we denote by $\dot{r}_1$ and $\dot{r}_2$ the velocities of the two beads, these three contributions can be easily identified in the three terms in

$$ -\zeta(\dot{r}_2 - u_0 - \text{grad}(u)r_2) - k_b T \frac{\partial}{\partial r_2} (\ln \Psi) - F^c = 0 \quad (1) $$

$$ -\zeta(\dot{r}_1 - u_0 - \text{grad}(u)r_1) - k_b T \frac{\partial}{\partial r_1} (\ln \Psi) + F^c = 0 \quad (2) $$

where $\zeta$ is the drag coefficient, $u$ the velocity field, $u_0$ the average velocity, $k_b$ the Boltzman constant, $T$ the absolute temperature and $\Psi$ is the probability distribution function. Using the definition of the connector vector $q = r_2 - r_1$ we can derive the following equation:

$$ \dot{q} = \text{grad}(u)q - \frac{2}{\zeta} \left( k_b T \frac{\partial}{\partial q} (\ln \Psi) + F^c(q) \right) \quad (3) $$

The connector force can take different forms leading to different kinetic models. In this work a non-linear extensible dumbbell (FENE) is considered. The connector force is then given by

$$ F^c(q) = \frac{H}{1 - \tilde{q}^2 / q_0^2} q \quad (4) $$

where $\tilde{q} = \|q\|$, $H$ is a connector constant and $q_0$ is the maximum spring length. A particularity of this model is that it has no equivalent constitutive macroscopic equation [7]. The associated evolution equation for the distribution function can be written as

$$ \frac{\partial \Psi}{\partial t} = -\frac{\partial}{\partial q} \left\{ \left( \text{grad}(u)q - \frac{2}{\zeta} F^c(q) \right) \Psi \right\} + \frac{2k_b T}{\zeta} \frac{\partial^2 \Psi}{\partial q^2} \quad (5) $$

The problem defined by Eq. (5) has a characteristic relaxation time $\theta = \zeta/4H$ and a dimensionless finite extensibility parameter $b = Hq_0^2/k_bT$. Thus vector $q$ can be made dimensionless with $\sqrt{k_bT/H}$, $\text{grad}(u)$ with $1/\theta$ (so it can be viewed as a Weissenberg number $We$), time with $\theta$ and the polymer stress tensor
with $n_c k_b T$ where $n_c$ is the number of chains per unit volume. Consequently, the dimensionless form of (5) is

$$ \frac{\partial \Psi}{\partial t} = -\frac{\partial}{\partial q} \left\{ \left( \text{grad}(u)q - \frac{1}{2}f(q)q \right) \Psi \right\} + \frac{1}{2} \frac{\partial^2 \Psi}{\partial q^2} \quad (6) $$

where $f(q)$ becomes the dimensionless connector force, that in the FENE model results:

$$ f(q) = \frac{1}{1 - \tilde{q}^2/b} \quad (7) $$

Moreover, a normalisation condition is associated with the probability distribution:

$$ \int \Psi(q) dq = 1 \quad (8) $$

Finally, the relation between statistical distribution of dumbbell configurations and the polymer stress $\tau_p$ is provided by the Kramers' expression:

$$ \tau_p = \langle f(q)q \otimes q \rangle - \underline{1} = \int \Psi(q) f(q)q \otimes q dq - \underline{1} \quad (9) $$

with $\underline{1}$ the unit tensor.

We must notice with respect to Eq. (6) that this equation defines the time evolution of the distribution function, whose integration requires to specify the initial distribution that we denote by $\Psi_0$. A reasonable choice lies in taking as initial distribution the equilibrium steady state related to a null velocity gradient. That distribution can be obtained by solving

$$ \frac{\partial}{\partial q} \left\{ \left( \frac{1}{2} f(q)q \right) \Psi_0 \right\} + \frac{1}{2} \frac{\partial^2 \Psi_0}{\partial q^2} = 0 \quad (10) $$

By symmetry considerations $\Psi_0 = \Psi_0(\tilde{q})$, and taking into account that $\Psi(\tilde{q} \to \infty) = 0$, the solution of the previous equation results

$$ \Psi_0(\tilde{q}) = \frac{f(\tilde{q})^{-b/2}}{\int_0^\infty f(\tilde{q})^{-b/2} d\tilde{q}} \quad (11) $$

Although in Eq. (6) $\frac{\partial \Psi}{\partial t}$ represents the material derivative (defined in the physical space) in non-homogeneous flows, in this paper, we consider only simple flows with homogeneous solutions. Thus, the material derivative reduces to the partial derivative.

In the case of complex flows with non-homogeneous solutions, out of the scope of the present work, we should consider the physical (related to the fluid element position) and configuration (related to the molecule orientation and extension) variables. When an updated Lagrangian technique is used to discretize the problem in the physical variables, the advection term of the material derivative can be integrated along the flow trajectories using the method of characteristics. Thus, no stabilization will be required for treating advection terms in physical variables. However, if one uses the finite element method for discretizing the problem in the physical space, the mesh becomes too distorted when the nodes move along the nodal trajectories, which implies the necessity of frequent remeshing. To alleviate this drawback, and avoid also the field projections required after each remeshing, one could proceed in a meshless framework. This question

have been deeply addressed in the context of non-Newtonian flows in [11,10].

### 3. Finite element discretisation

In this section we establish the discrete form related to the partial differential equation governing the evolution of the probability distribution (6) in the finite element framework. Taking into account the solution homogeneity, no discretisation in physical variables is required, and therefore, we proceed by discretizing with respect to the conformation coordinates. Indices $v$ affecting a vector variable refer to their components. The components of a second order tensor will be referred using two indices. In this way, Eq. (6) can be written as

$$ \frac{\partial \Psi}{\partial t} + E_0(q) \Psi + \sum_v \left( E_{1v}(q) \frac{\partial \Psi}{\partial x_v} + E_2 \frac{\partial^2 \Psi}{\partial x_v^2} \right) = 0 \quad (12) $$

with

$$ \begin{aligned} E_0(q) &= \sum_v \left( u_{v,v} - \frac{1}{2} f(q) - \frac{1}{2} x_v f_{,v} \right), \\ E_{1v}(q) &= \sum_{v'} \left( u_{v,v'} x_{v'} \right) - \frac{1}{2} f(q) x_v, \quad E_2 = -1/2 \end{aligned} \quad (13) $$

where $f_{,v} = \frac{\partial f(q)}{\partial x_v}$ and the velocity gradient is given by $u_{v,v'} = \frac{\partial u_v}{\partial x_{v'}}$.

The configuration domain is bounded. This domain is chosen such that the distribution function can be assumed vanishing on its boundary. In the case of the FENE model, this domain can be defined by

$$ \Omega = [-b, b]^N $$

being $\mathcal{N}$ is the dimension of the space where the molecular conformation is defined. Of course, with such definition of the domain $\Omega$, the dimensionless connector force is not defined in some regions, in particular for points $x$ such that $\|x\| > \sqrt{b}$ according to Eq. (7). At those points, we simply consider the connector force given by a constant value large enough to avoid a molecular extension higher than its maximum physical value.

Firstly, the problem is formulated in the finite element framework using a weight function $\Psi^*$.

$$ \int_{\Omega} \Psi^* \frac{\partial \Psi}{\partial t} dq + \int_{\Omega} \Psi^* E_0(q) \Psi dq + \sum_{v} \left( \int_{\Omega} \Psi^* E_{1v}(q) \frac{\partial \Psi}{\partial x_v} dq + \int_{\Omega} \Psi^* E_2 \frac{\partial^2 \Psi}{\partial x_v^2} dq \right) = 0 \quad (14) $$

The computational domain $\Omega$ is partitioned into a collection of non-overlapping finite elements. The number of nodes defining the $C^0(\Omega)$ interpolation in each element are 2 in 1D, 4 in 2D (usual $Q_1$-rectangular element) and 8 in 3D (usual $Q_1$-hexahedron element). Thus, the number of nodes in each element
is equal to $2^N$, $\mathcal{N}$ being the space dimension. Thus, we can write

$$ \Psi^e(q) = \sum_{i=1}^{2^N} N_i(q) \Psi_i^e \qquad (15) $$

and

$$ \Psi^{e*}(q) = \sum_{i=1}^{2^N} N_i(q) \Psi_i^{e*} \qquad (16) $$

where $\Psi_i^e$ and $\Psi_i^{e*}$ are the values at node $i$ of $\Psi$ and $\Psi^*$ respectively, and $N_i(q)$ is the associated shape function which takes a unit value at the node $i$, and vanishes at the other nodal positions.

Assuming that the distribution vanishes on the domain boundary $\partial\Omega$, Eq. (14) can be written after integration by part, as

$$ \int_{\Omega} \Psi^* \frac{\partial \Psi}{\partial t} dq + \int_{\Omega} \Psi^* E_0(q) \Psi dq + \sum_{\nu} \left( \int_{\Omega} \Psi^* E_{1\nu}(q) \frac{\partial \Psi}{\partial x_{\nu}} dq - \int_{\Omega} \frac{\partial \Psi^*}{\partial x_{\nu}} E_2 \frac{\partial \Psi}{\partial x_{\nu}} dq \right) = 0 \qquad (17) $$

Due to the advection–diffusion character of Eq. (17) an appropriate stabilization is needed to avoid numerical instabilities induced by the convection term. A non-consistent upwinding formulation is considered here, which modifies the weighting function related to the advection term $\bar{\Psi}^*$ as described later. Thus, the stabilized variational formulation results:

$$ \int_{\Omega} \Psi^* \frac{\partial \Psi}{\partial t} dq + \int_{\Omega} \Psi^* E_0(q) \Psi dq + \sum_{\nu} \left( \int_{\Omega} \bar{\Psi}_{\nu}^* E_{1\nu}(q) \frac{\partial \Psi}{\partial x_{\nu}} dq - \int_{\Omega} \frac{\partial \Psi^*}{\partial x_{\nu}} E_2 \frac{\partial \Psi}{\partial x_{\nu}} dq \right) = 0 \qquad (18) $$

where in each element

$$ \bar{\Psi}^{e*}(q) = \sum_{i=1}^{2^N} \bar{N}_i(q) \Psi_i^{e*} \qquad (19) $$

with

$$ \bar{N}_i(q) = N_i(q) + \frac{\beta h}{2 \|E_1(q)\|} \frac{\partial N_i(q)}{\partial x_{\nu}} E_{1\nu}(q) \qquad (20) $$

where $h$ is the characteristic element length (in the advection direction) and $\beta$ is the upwinding parameter given by

$$ \beta = \coth(Pe) - \frac{1}{Pe} \qquad (21) $$

where the Peclet number $Pe$ is given by

$$ Pe = \frac{\|E_1(q)\|h}{2|E_2|} \qquad (22) $$

Integration of Eq. (18) leads to the system of ordinary differential equations

$$ \underline{\Psi}^{*T} \underline{M} \underline{\Psi} + \underline{\Psi}^{*T} \underline{G}_0 \underline{\Psi} + \underline{\Psi}^{*T} \underline{G}_1 \underline{\Psi} + \underline{\Psi}^{*T} \underline{G}_2 \underline{\Psi} = 0 \qquad (23) $$

**Remark 3.1.** From now on, for the sake of simplicity, vector will be identified by an underline, which is doubled in the case of a matrix. However, $q$ continues to refer the connector vector, for which the previous rule is not applied.

Eq. (23) can be rewritten by introducing the matrix $\underline{G}$, with $\underline{G} = \underline{G}_0 + \underline{G}_1 + \underline{G}_2$, as

$$ \underline{M} \underline{\Psi} + \underline{G} \underline{\Psi} = \underline{0} \qquad (24) $$

This equation can be solved, under the normalisation condition, by applying an appropriate time discretisation preserving the stability and assuring the convergence.

A fully implicit (backwards Euler) scheme proceeds as follows: *From $\underline{\Psi}_n$ just computed ($\underline{\Psi}_0$ being known), compute $\underline{\Psi}_{n+1}$ by solving:*

$$ \underline{\Psi}_{n+1} = (\underline{M} + \Delta t \underline{G})^{-1} \underline{M} \underline{\Psi}_n $$

The update requires the solution of a linear system whose size is directly related to the number of degrees of freedom used in the space discretisation. For this reason, the higher is the space dimension, the larger is the computing cost. In the next section we propose a strategy which allows for significant computing savings. We will prove that the size of the resulting problem depends only on the solution regularity.

# 4. A priori model reduction

## 4.1. Introducing the main ideas

Let us consider the evolution of a field $\Psi_{\text{ref}}(q, t)$ defined at any time $t \in [0, T]$ and for any configuration $q$ in the domain $\Omega$. The purpose of the Karhunen–Loève (KL) expansion of $\Psi$ is to find a set of basis functions $(\gamma_k(q))_{k=1,...,r}$ defined over $\Omega$ and a set of reduced state variables $(\nu_k(t))_{k=1,...,r}$ defined over $[0, T]$ such that

(1) Preliminary stage. $\gamma_k$ is defined by the stationarity of $\lambda(\gamma)$:

$$ \lambda(\gamma) = \frac{\int_{0}^{T} (\int_{\Omega} \Psi_{\text{ref}}(q, t) \gamma(q) dq)^2 dt}{\int_{\Omega} \gamma^2(q) dq} \qquad (25) $$

$$ \delta\lambda(\gamma_k) = 0 \qquad (26) $$

$$ \lambda(\gamma_k) \neq 0 \qquad (27) $$

$$ \int_{\Omega} \gamma_k^2(q) dq = 1 \qquad (28) $$

(2) Projection stage. $(\nu_k(t))_{k=1,...,r}$ minimize $\eta(\nu_1, ..., \nu_r)$:

$$ \eta(\nu_1, \dots, \nu_r) = \int_0^T \int_{\Omega} \left( \Psi_{\text{ref}}(q, t) - \sum_{k=1}^{k=r} \gamma_k(q) \nu_k(t) \right)^2 \times dq dt \qquad (29) $$

The adaptive method that we propose allows to avoid the preliminary construction of basis functions. Thus, during the evolution of $\Psi(q, t)$, the approximation basis $(\phi_k^{(m)})_{k=1,...,r(m)}$ is improved. The superscript $(m)$ indicates that the approximation
basis has been updated *m* times until now. We can prove that the basis functions $\phi_k^{(m)}$ are approximations of the eigenfunction $\gamma_k$ [14]. The evolution of $\Psi(q, t)$ can be approximated by

$$ \Psi_{\phi}^{(m)}(q, t) = \sum_{k=1}^{k=r(m)} \phi_k^{(m)}(q) a_k^{(m)}(t) \quad (30) $$

The adaptive procedure that we describe later includes the enrichment of the approximation basis and the selection of the most significant basis functions in order to represent the state evolution.

Let us consider the weak form of the equations defining the evolution of $\Psi$:

$$ \Psi(q, 0) = \Psi_0(q) \quad (31) $$

$$ R \left( \Psi^*, \Psi, \frac{\partial \Psi}{\partial t}, t \right) = 0, \quad \forall \Psi^* \quad (32) $$

where $R$ is the residual of the governing equations. We denote by $\Psi_{\text{ref}}$ the exact solution of these governing equations.

The norm of the residual $R(\Psi^*, \Psi_{\phi}^{(m)}, \frac{\partial \Psi_{\phi}^{(m)}}{\partial t}, t)$ provides an error estimator to check the quality of $\Psi_{\phi}^{(m)}$. An approximation is appropriate if the residual is lower than a given value $\varepsilon_R$ (small enough) for all time $t \in [0, T]$.

Different adaptive strategies can be developed, all of them including:

(1) A set of basis functions always defined in the whole time interval.

(2) An extension of the approximation basis from the residual of the governing equations.

(3) A selection of the most significant approximation functions from a Karhunen–Loève expansion of the reduced state variables evolution.

The strategy that we propose is an incremental approach. Thus, if we assume that the evolution approximation was done properly over the time interval $[0, t_\alpha]$ thanks to the reduced state variables $a_k^{(m)}$ associated to the basis functions $\phi_k^{(m)}$, $\forall k \in [1, \dots, r^{(m)}]$, the first step of the algorithm is to find the next time $t_\beta$ ($t_\beta \ge t_\alpha$) such that the residual computed at $t_\beta$ does not satisfy the quality criterion. When $t_\beta$ is known, an adaptation of the basis functions is performed for satisfying the quality criterion at $t_\beta$, and then the incremental computation can continue until the new basis adaptation. Obviously, to start the computation we can choose $\phi_1^{(0)}$ collinear to $\Psi_0$.

After this general overview, we now focus on the adaptation procedure.

The adaptation is carried out in two steps for satisfying the quality criterion:

(1) The first step is the KL expansion of the reduced state variables over the time interval $[0, t_\beta]$. Let $C^{(m)}$ be the averaged autocorrelations matrix

$$ C_{kp}^{(m)} = \int_{0}^{t_{\beta}} a_{k}^{(m)}(t) a_{p}^{(m)}(t) dt, \quad k, p \in \{1, \dots, r^{(m)}\} \quad (33) $$

Let $\mu_j$ be the eigenvalue related to the eigenvector $\underline{V}_j^{(m)}$ such that:

$$ \mu_1 \ge \mu_2 \ge \dots \ge \mu_{r(m)} \quad (34) $$

with

$$ \underline{V}_i^{(m)} \cdot \underline{V}_j^{(m)} = \sum_{k=1}^{k=r(m)} V_{ki}^{(m)} V_{kj}^{(m)} = \delta_{ij} \quad (35) $$

The selection of the most significant basis functions is performed by taking the $\tilde{r}^{(m+1)}$ eigenvectors related to the eigenvalues greater than $\varepsilon_{\text{KL}}\mu_1$, i.e.:

$$ \mu_i \ge \varepsilon_{\text{KL}}\mu_1, \quad i \in [1, \dots, \tilde{r}^{(m+1)}] \quad (36) $$

where $\varepsilon_{\text{KL}}$ is a small enough parameter ($\varepsilon_{\text{KL}} = 10^{-8}$ in our simulations). The selected eigenvectors define the basis reduction matrix $\underline{V}^{(m)}$

$$ \underline{V}^{(m)} = [\underline{V}_1^{(m)}, \dots, \underline{V}_{\tilde{r}^{(m+1)}}^{(m)}] \quad (37) $$

This matrix provides new basis functions $(\phi_k^{(m+1)})_{k=1,\dots,\tilde{r}^{(m+1)}}$ such that

$$ \phi_k^{(m+1)}(q) = \gamma_k \sum_{p=1}^{p=r(m)} \phi_p^{(m)}(q) V_{pk}^{(m)}, \quad k = 1, \dots, \tilde{r}^{(m+1)} \quad (38) $$

where $\gamma_k$ is computed in order to normalize the basis functions

$$ \int_{\Omega} \phi_k^{(m+1)}(q) \phi_k^{(m+1)}(q) dq = 1 \quad (39) $$

(2) The second step consists of the extension of the subspace spanned by the functions $(\phi_k^{(m+1)})_{k=1,\dots,\tilde{r}^{(m+1)}}$ by adding the function $\phi_{\tilde{r}^{(m+1)}}(r^{(m+1)} = \tilde{r}^{(m+1)} + 1)$ collinear to the residual at time $t_\beta$, i.e.

$$ \int_{\Omega} \Psi^*(q) \phi_{\tilde{r}^{(m+1)}}^{(m+1)}(q) dq = \gamma_{\tilde{r}^{(m+1)}} R \left( \Psi^*, \Psi_{\phi}^{(m)}, \frac{\partial \Psi_{\phi}^{(m)}}{\partial t}, t \right), \quad \forall \Psi^*, t = t_{\beta} \quad (40) $$

where $\gamma_{\tilde{r}^{(m+1)}}$ is taken to verify

$$ \int_{\Omega} \phi_{\tilde{r}^{(m+1)}}^{(m+1)}(q) \phi_{\tilde{r}^{(m+1)}}^{(m+1)}(q) dq = 1 \quad (41) $$

Because of the basis adaptation, the reduced state variables must be updated over $[0, t_\beta]$ [14]:

$$ a_k^{(m+1)}(t) = \frac{1}{\gamma_k} \sum_{p=1}^{p=r(m)} a_p^{(m)}(t) V_{pk}^{(m)}, \quad \forall t < t_\beta \forall k < r^{(m+1)} \quad (42) $$

and

$$ a_{\tilde{r}^{(m+1)}}^{(m+1)}(t) = 0, \quad \forall t < t_{\beta} \quad (43) $$
For $t \ge t_β$, the evolution process is restarted using the just updated reduced approximation basis:

$$ \Psi^*(q) = \sum_{k=1}^{k=r^{(m+1)}} \phi_k^{(m+1)}(q) a_k^* \quad (44) $$

$$ \Psi_{\phi}^{(m+1)}(q, t) = \sum_{k=1}^{k=r^{(m+1)}} \phi_k^{(m+1)}(q) a_k^{(m+1)}(t) \quad (45) $$

$$ R \left( \Psi_{\phi}^{*}, \Psi_{\phi}^{(m+1)}, \frac{\partial \Psi_{\phi}^{(m+1)}}{\partial t}, t \right) = 0, \quad \forall \Psi_{\phi}^{*}, \forall t \in [t_{\alpha}, t_{\beta}] \quad (46) $$

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

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

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

## 4.2. Reduction of kinetic theory models

In this section we apply the ideas just described to the discrete problem defined by Eq. (24).

We consider that the probability distribution has been accurately described in [0, $t_α$]. We assume that at time $t_β$ the reduced approximation basis is given by $\underline{B}^{(m)}$. The i, j-component of $\underline{B}$ represents the value at node i of the j-eigenvector ($j \in [1, r^{(m)}]$). Moreover, at certain times $t_p, p \in [1, \dots, P]$, the solution is assumed properly computed and defined by the reduced vectors $\underline{a}_p^{(m)}$. Knowing $\underline{a}_p^{(m)}$, the finite element description of $\Psi$ at time $t_p$ results: $\underline{\Psi}_p = \underline{B}^{(m)} \underline{a}_p^{(m)}$. We can assume that the first approximation basis $\underline{B}^{(0)}$ contains a single vector that corresponds to the initial probability distribution $\underline{\Psi}_0$.

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

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

<figure><img src="image_7.png" /><figcaption>Fig. 1. 1D start-up elongation $We = 5$, $b = 10$: (a) distribution function evolution; (b) difference between the exact and computed steady solutions; (c) functions defining the reduced approximation basis; (d) evolution of the reduced order approximation coefficients; (e) evolution of the number of approximation functions; (f) evolution of the polymer stress.</figcaption></figure>
Now, we can compute the evolution of $\Psi$ in $[t_\alpha, t_\beta]$ by solving

$$ \underline{\Psi}_{n+1} = (\underline{M} + \Delta t \underline{G})^{-1} \underline{M} \underline{\Psi}_n \quad (47) $$

that can be written in the reduced approximation basis as

$$ (\underline{B}^{(m)})^{\mathrm{T}} \underline{B}^{(m)} \underline{a}_{n+1}^{(m)} = (\underline{B}^{(m)})^{\mathrm{T}} (\underline{M} + \Delta t \underline{G})^{-1} \underline{M} \underline{B}^{(m)} \underline{a}_{n}^{(m)} \quad (48) $$

This results in

$$ \underline{a}_{n+1}^{(m)} = ((\underline{B}^{(m)})^{\mathrm{T}} \underline{B}^{(m)})^{-1} (\underline{B}^{(m)})^{\mathrm{T}} (\underline{M} + \Delta t \underline{G})^{-1} \underline{M} \underline{B}^{(m)} \underline{a}_{n}^{(m)} \quad (49) $$

At time $t_\beta$ we compute the residual according to

$$ \underline{R}_{\beta} = \underline{B}^{(m)} \underline{a}_{\beta+\Delta t}^{(m)} - (\underline{M} + \Delta t \underline{G})^{-1} \underline{M} \underline{B}^{(m)} \underline{a}_{\beta}^{(m)} \quad (50) $$

<figure><img src="image_8.png" /><figcaption>Fig. 2. 1D start-up elongation We = 10, b = 50: (a) distribution function evolution; (b) difference between the exact and computed steady solutions; (c) functions defining the reduced approximation basis; (d) evolution of the reduced order approximation coefficients; (e) evolution of the number of approximation functions; (f) evolution of the polymer stress.</figcaption></figure>

If $\|\underline{R}_\beta\| < \epsilon_R$, we can put $t_\alpha = t_\beta$ and continue the resolution of the evolution problem using the reduced approximation basis $\underline{B}^{(m)}$ for $t \ge t_\alpha$. When $\|\underline{R}_\beta\| > \epsilon_R$ we need to look for a smaller time $t_\beta$ ($t_\alpha + \frac{t_\beta - t_\alpha}{2} \to t_\beta$) such that $\|\underline{R}_\beta\| < \epsilon_R$ after a new resolution of the evolution problem in $[t_\alpha, t_\beta]$. If the residual criterion is unsatisfied when $t_\beta - t_\alpha < \Delta t$, then we put $t_{P+1} = t_\alpha$ and proceed to an enrichment of the approximation basis, that is performed as follows.

The new reduced approximation basis $\underline{B}^{(m+1)}$ is defined by adding some Krylov's subspaces to the significant information extracted from $\underline{a}_p^{(m)}$, $\forall p \in [1, \dots, P+1]$. We now explain the construction of $\underline{B}^{(m+1)}$.
We define the matrix $\underline{Q}$ containing the reduced vectors $\underline{a}_p^{(m)}$, $\forall p$. We solve the eigenproblem defined by

$$ \underline{Q} \underline{Q}^T \underline{\phi} = \mu \underline{\phi} \qquad (51) $$

whose solution results in $r^{(m)}$ couples $(\underline{\phi}_k, \mu_k)$, where we assume that $\mu_1 \ge \mu_2, \dots, \ge \mu_{r^{(m)}}$. We select the $\tilde{r}^{(m+1)}$ eigenvectors $\underline{\phi}_k$ related to the eigenvalues verifying $\mu_k > 10^{-8}\mu_1$, defining the matrix $\underline{\phi}$ of dimension $r^{(m)} \times \tilde{r}^{(m+1)}$. We can write

$$ \tilde{\underline{B}}^{(m+1)} = \underline{\tilde{B}}^{(m)} \underline{\phi}^{(m)} \qquad (52) $$

Obviously, the change in the reduced approximation basis implies a change in the expression of the reduced vectors $\underline{a}_p^{(m)}, \forall p$. For this purpose we can write

$$ \underline{\tilde{B}}^{(m+1)} \underline{\tilde{a}}_p^{(m+1)} = \underline{\underline{B}}^{(m)} \underline{a}_p^{(m)} \qquad (53) $$

from which we have

$$ (\underline{\tilde{B}}^{(m+1)})^T \underline{\tilde{B}}^{(m+1)} \underline{\tilde{a}}_p^{(m+1)} = (\underline{\tilde{B}}^{(m+1)})^T \underline{\underline{B}}^{(m)} \underline{a}_p^{(m)} \qquad (54) $$

or

$$ \underline{\tilde{a}}_p^{(m+1)} = ((\underline{\tilde{B}}^{(m+1)})^T \underline{\tilde{B}}^{(m+1)})^{-1} (\underline{\tilde{B}}^{(m+1)})^T \underline{\underline{B}}^{(m)} \underline{a}_p^{(m)}, \quad \forall p \qquad (55) $$

Now, we can add to $\underline{\tilde{B}}^{(m+1)}$ some (3 in our simulations) Krylov's subspaces defined at time $t_\beta = t_\alpha$:

$$ \underline{\underline{B}}^{(m+1)} = (\underline{\tilde{B}}^{(m+1)}, \underline{R}_\beta, \underline{H} \underline{R}_\beta, \underline{H}^2 \underline{R}_\beta) \qquad (56) $$

where $\underline{R}_\beta$ is defined by Eq. (50) and the matrix $\underline{H}$ is given by

$$ \underline{H} = \underline{M}^{-1}(\underline{M} + \Delta t \underline{G}) $$

The introduction of the Krylov's subspaces implies a new change in the expression of the reduced vectors $\underline{\tilde{a}}_p^{(m+1)}$. For this purpose we add to those vectors 3 new components that are assumed null because the new approximation functions just introduced operate only for $t > t_{P+1}$. Thus we can write

$$ (\underline{a}_p^{(m+1)})^T = ((\underline{\tilde{a}}_p^{(m+1)})^T, 0, 0, 0), \quad \forall p \qquad (57) $$

The following algorithm summarizes the procedure described above.

$$ t = t_\alpha = 0, t_\beta = t_\alpha + \tau (\tau \gg \Delta t) \\ \text{While } t_\alpha < T $$

$$ \text{While } t_\beta - t_\alpha > \Delta t $$

Compute the evolution of the reduced variables in $[t_\alpha, t_\beta]$

If $\|\underline{R}(t = t_\beta)\| < 0.001$ then:

$$ \delta = t_\beta - t_\alpha $$

$$ t_\alpha = t_\beta $$

$$ t_\beta = t_\alpha + 2\delta $$

If $\|\underline{R}(t = t_\beta)\| > 0.001$ then:

$$ \delta = t_\beta - t_\alpha $$

$$ t_\alpha = t_\beta $$

$$ t_\beta = t_\alpha + \frac{\delta}{2} $$

Karlhunen–Loève expansion in $[0, t_\alpha]$

Enrichment from the Krylov's subspaces related to the residual at $t_\alpha$

$$ t_\beta = t_\alpha + \tau $$

End.

<figure><img src="image_9.png" /><figcaption>Fig. 3. 2D start-up elongation $We = 5, b = 10$: (a) distribution function $t = 0.0$; (b) distribution function $t = 0.1$; (c) distribution function $t = 0.2$; (d) distribution function $t = 2.0$.</figcaption></figure>
## 5. Numerical examples

The first numerical example concerns a one-dimensional start-up elongation problem defined by the velocity field $u = 5 \cdot x$ which corresponds to the elongation rate grad($u$) = $We = 5$. The domain is partitioned using 315 linear 1D finite elements of size $h_1 = 0.01$. The time step is fixed to $\Delta t = 10^{-5}$ and the parameter $b$ is set to 10. The initial solution is given by the equilibrium state defined in Section 2. Basis enrichment operates when $\epsilon_R$ exceed 0.001. Each enrichment is performed using the three first Krylov subspaces. The Karhunen-Loève decomposition allows to select the most significant basis functions according to the criterion $\mu/\mu_{max} > 10^{-8}$ ($\mu$ being the problem eigenvalues). The total time is $T = 2$ for which the steady state is almost reached. Fig. 1(a) shows the evolution of the solution at times $t = 0, 0.1, 0.2, 0.5, 1, 1.5$ and 2. The finite element description uses 316 shape functions. In our case only 6 significant functions are found to represent the evolution of the solution during the

whole time interval (see Fig. 1(c)). The evolution of the coefficients affecting the different approximation functions is depicted in Fig. 1(d). The basis enrichment occurs mainly at the beginning of the simulation due to the fast solution changes taking place. In the same way the number of approximation functions needed to represent the solution evolution increases at the beginning of the simulation (Fig. 1(e)) and becomes constant when the solution evolution becomes smoother (the number of approximation functions is represented at the times at which the approximation basis enrichment takes place). The polymer stress is depicted in Fig. 1(f). The exact steady solution can be obtained in this case to evaluate the accuracy of the computed solution involving only six approximation functions. The difference between the computed steady solution and the exact one is plotted in Fig. 1(b).

Now, we consider the benchmark problem done in [7]. The problem lies in a one-dimensional start-up elongation ($u = 10 \cdot x$) being $We = 10$. The domain has been partitioned into 705

<figure><img src="image_10.png" /><figcaption>Fig. 4. 2D start-up elongation $We = 5, b = 10$: (a) first most significant function; (b) second most significant function; (c) third most significant function; (d) fourth most significant function; (e) evolution of the reduced order approximation coefficients; (f) evolution of the number of approximation functions.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

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

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

<figure><img src="image_14.png" /><figcaption>Fig. 5. 3D start-up elongation We = 5, b = 10: (a) distribution function $t = 0.0$; (b) distribution function $t = 0.1$; (c) distribution function $t = 0.2$; (d) distribution function $t = 2.0$.</figcaption></figure>

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

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

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

<figure><img src="image_18.png" /><figcaption>Fig. 6. 3D start-up elongation We = 5, b = 10: (a) first most significant function; (b) second most significant function; (c) evolution of the reduced order approximation coefficients; (d) evolution of the number of approximation functions.</figcaption></figure>
linear 1D finite elements of size $h_1 = 0.01$, being the time step $\Delta t = 10^{-6}$ and the parameter $b$ is set to 50. Fig. 2(a) shows the evolution of the solution at times $t = 0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5$ and $1$. This simulation has been done using only 11 functions (rather than 706 needed without a model reduction). The five most representative functions are plotted in Fig. 2(c). The evolution of the associated coefficients of the distribution function expressed in the reduced basis is depicted in Fig. 2(d). According to Fig. 2(e) it seems that the basis enrichment is active when the solution approximates the highest gradients (molecular extensions close to $\sqrt{b}$). The number of approximation functions becomes constant for $t > 0.4$. The polymer stress is depicted in Fig. 2(f) in perfect agreement with the results reported in [7]. Finally Fig. 2(b) depicts the error between exact and computed steady solutions.

It can be noticed that the above problems are symmetric, and as a consequence only half of the conformation domain has been considered in the computations. The extension of the problem analysis to the whole domain implies double of the number of degrees of freedom if one proceeds in the context of a standard finite element discretisation. On the contrary, when one proceeds with a model reduction technique the number of significant functions remains unchanged, and thus the size of the problem does not depend on the size of the domain where the model is defined, but only on the solution regularity. Moreover, we will prove in the following 2D and 3D examples that the number of significant approximation functions depends on the solution regularity but it does not depend on the dimension of the configuration space.

We consider the 2D start-up elongation flow characterized by the velocity gradient $\text{grad}(u) = \text{diag}(1, -1)We$. The finite element size is defined by the element length in the x-direction, $h_1 = 0.02$, and in the y-direction $h_2 = 0.05$, the time step being $\Delta t = 10^{-5}$ and the parameter $b = 10$. The evolution of the solution at times $t = 0, 0.1, 0.2$ and $2$ is shown in Fig. 3(a)–(d). The reduced order simulation involves only 10 approximation functions, whose four most significant ones are represented in Fig. 4(a)–(d). The evolution of the associated coefficients are depicted in Fig. 4(e). As previously, the basis enrichment operates mainly at the beginning of the simulation according to the Fig. 4(f).

Finally, a 3D start-up elongation flow is considered. The problem symmetry is considered again. Planes $x_1 = 0, x_2 = 0$ and $x_3 = 0$ are used for graphical representation of the solution evolution and the different basis functions. The flow is now characterized by $\text{grad}(u) = \text{diag}(1, -0.5, -0.5)We$. The element size is defined by $h_1 = 0.03$ and $h_2 = h_3 = 0.1$, being $\Delta t = 10^{-5}$ and $b = 10$. The solution computed at times $t = 0, 0.1, 0.2$ and $2$ is shown in Fig. 5(a)–(d). This simulation involves 11 basis functions, whose two most significant ones are represented in Fig. 6(a) and (b). It must be noticed that the number of approximation functions involved is almost the same as in the 2D case. The evolution of the associated reduced basis coefficients is depicted in Fig. 6(c) and the evolution of the number of approximation functions in Fig. 6(d). In this case the exact steady solution can be obtained, which allows one to conclude about the excellent accuracy of the proposed reduced order technique.

<figure><img src="image_19.png" /><figcaption>Fig. 7. 3D steady state elongation $We = 5, b = 10$: (a) exact steady distribution function; (b) distribution computed by searching the steady solution.</figcaption></figure>

Moreover, we can notice that the model reduction strategy described in the present paper allows the direct computation of the steady solution of the kinetic model, i.e. the steady solution of the Fokker-Planck equation, when it exists, in contrast to usual stochastic techniques. For this purpose, the technique uses an initial prediction that is used in the initial reduced approximation basis, which is enriched at each iteration using the Krylov's subspaces. Fig. 7 compares the exact 3D steady distribution in the elongation flow just considered [2] with the one computed by using the reduced model technique when the steady solution is directly searched without considering the evolution problem. As it is noticed both solutions are in perfect agreement.

## 6. Conclusions

We have proposed in this paper a new model reduction technique allowing accurate and fast simulations of kinetic theory models. The new approximation functions used to define the reduced model are defined in the whole domain in an appropriate manner (the most characteristic functions related to the model solution). Thus, the number of degrees of freedom involved in the solution of the Fokker-Planck equation is significantly reduced. The construction of those new approximation functions is done with an 'a priori' approach, which combines a basis reduction (using the Karhunen-Loève decomposition)
with a basis enrichment based on the use of some Krylov
subspaces. We have illustrated the potentiality of that technique
by solving the FENE viscoelastic model in 1D, 2D and 3D. We
noticed that the simulations involved the same number (around
10) of approximation functions. Thus, we can conclude that the
number of significant approximation functions used to define
the reduced model depends on the solution regularity but does
not depend significantly on the dimension of the space.

The extension of these ideas to multi-bead-spring (MBS)
models is not direct because such models are defined in a multi-
dimensional conformation space, and even if a reduced number
of approximations functions are enough for describing accu-
rately the whole evolution, these function cannot be defined us-
ing a grid support (whose nodes evolves exponentially with the
space dimension). One possibility consists of the use of appro-
priate tensor products which constitutes one of our main works
in progress.

References

[1] A. Ammar, F. Chinesta, A particle strategy for solving Fokker-Planck equation governing the fiber orientation distribution in steady recirculating flows involving short fiber suspensions, Lecture Notes on Computational Science and Engineering. Springer 43 (2004) 1–16.

[2] B.B. Bird, C.F. Curtiss, R.C. Armstrong, O. Hassager, Dynamics of polymeric liquids, in: Kinetic Theory, vol. 2, John Wiley & Sons, 1987.

[3] C.V. Chaubal, A. Srinivasan, O. Egecioglu, L.G. Leal, Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems, J. Non-Newtonian Fluid Mech. 70 (1997) 125–154.

[4] C. Chauviere, A. Lozinski, Simulation of dilute polymer solutions using a Fokker-Planck equation, Comput. Fluids 33 (2004) 687–696.

[5] F. Chinesta, G. Chaidron, A. Poitou, On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions, J. Non-Newtonian Fluid Mech. 113 (2-3) (2004) 97-125.

[6] R. Jendrejack, J. de Pablo, M. Graham, A method for multiscale simulation of flowing complex fluids, J. Non-Newtonian Fluid Mech. 108 (2002) 123–142.

[7] R. Keunings, On the Peterlin approximation for finitely extensible dumb-ells, J. Non-Newtonian Fluid Mech. 68 (1997) 85-100.

[8] R. Keunings, Micro-macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory, in: D.M. Binding, K. Walters (Eds.), Rheology Reviews, British Society of Rheology, 2004, pp. 67–98.

[9] A. Lozinski, C. Chauviere, A fast solver for Fokker-Planck equation applied to viscoelastic flows calculations: 2D FENE model, J. Comput. Phys. 189 (2003) 607–625.

[10] M.A. Martinez, E. Cueto, I. Alfaro, M. Doblar, F. Chinesta, Updated Lagrangian free surface flow simulations with the natural neighbour Galerkin methods, Int. J. Numer. Meth. Eng. 60 (2004) 2105–2129.

[11] M.A. Martinez, E. Cueto, M. Doblar, F. Chinesta, Natural element meshless simulation of injection processes involving short fiber suspensions, J. Non-Newtonian Fluid Mech. 115 (2003) 51–78.

[12] H.C. Öttinger, M. Laso, Smart polymers in finite element calculation, International Congress on Rheology, Brussel, Belgium, 1992.

[13] D. Ryckelynck, A priori hyperreduction method: an adaptive approach, J. Comput. Phys. 202 (2005) 346–366.

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

[15] M. Somasi, B. Khomami, N.J. Woo, J.S. Hur, E.S.G. Shaqfeh, Brownian dynamics simulations of bead-rod and bead-spring chains: numerical algorithms and coarse-graining issues, J. Non-Newtonian Fluid Mech. 108 (1-3) (2002) 227–255.

[16] G. Venkiteswaran, M. Junk, A QMC approach for high dimensional Fokker-Planck equations modelling polymeric liquids, Math. Comput. Simulat. 68 (2005) 43–56.