# On the Simulation of Kinetic Theory Models of Complex Fluids Using the Fokker-Planck Approach

Béchir Mokdad, Etienne Prulière, Amine Ammar, Francisco Chinesta

► To cite this version:

Béchir Mokdad, Etienne Prulière, Amine Ammar, Francisco Chinesta. On the Simulation of Kinetic Theory Models of Complex Fluids Using the Fokker-Planck Approach. Applied Rheology, 2007, 17/2 (26494), pp.1-14. 10.3933/ApplRheol-17-26494. hal-01007154

HAL Id: hal-01007154

https://hal.science/hal-01007154

Submitted on 10 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# ON THE SIMULATION OF KINETIC THEORY MODELS OF COMPLEX FLUIDS USING THE FOKKER-PLANCK APPROACH

B. MOKDAD, E. PRULIERE, A. AMMAR, F. CHINESTA*

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

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

**ABSTRACT:**

Models of kinetic theory provide a coarse-grained description of molecular configurations wherein atomistic processes are ignored. The Fokker-Planck equation related to the kinetic theory descriptions must be solved for the distribution function in both physical and configuration spaces. When the model involves high dimensional spaces (including physical and conformation spaces and time) standard discretization techniques fail due to excessive computation requirements. In this paper, we revisit some model reduction techniques recently proposed to circumvent those difficulties, exploring other new application areas related to entangled polymer models as well as the use of such reduced models for treating complex flows in which the distribution function involves both the physical and the conformation coordinates.

**ZUSAMMENFASSUNG:**

Modelle der kinetischen Theorie stellen eine vergröberte Beschreibung molekularer Konfigurationen zur Verfügung, bei der atomistische Details vernachlässigt werden. Die zur kinetischen Theorie gehörende Fokker-Planck-Gleichung muss für eine Verteilungsfunktion sowohl im physikalischen als auch im Konfigurationsraum gelöst werden. Wenn das Modell dabei einen hochdimensionalen Raum erfordert (einschlisslich physikalischer Raum, Konfigurationsraum und Zeit) versagen konventionelle Diskretisierungstechniken aufgrund ihrer sehr hohen Anforderungen an die Rechenzeit. Mit diesem Beitrag geben wir einen Überblick über kürzlich vorgeschlagene Modell-Reduktionstechniken, die in der Lage sind, diese Schwierigkeiten zu überwinden, und mithilfe derer neue Anwendungsgebiete, in denen die Verteilungsfunktionen sowohl die physikalische als auch konfigurelle Dimension beschreiben muss, angegangen werden können. Dazu zählen insbesondere verschlaufte Polymersysteme in komplexen Strömungen.

**RÉSUMÉ:**

Les modèles de théorie cinétique définissent une description globale de la configuration moléculaire dans laquelle l'échelle atomistique est ignorée. L'équation de Fokker-Planck relative à la théorie cinétique doit être résolue en donnant la fonction de distribution aussi bien sur l'espace physique que de configuration. Quand le modèle se définit sur un espace de dimension élevée (en superposant l'espace physique, de configuration et temporel) les techniques usuelles de discrétisation se heurtent à la difficulté du traitement numérique excessivement gourmant en temps de calcul et en stockage. Dans ce travail, on présente une revue de quelques techniques de réduction en explorant des nouvelles applications liées à la reptation des polymères fondus et également aux traitements d'écoulements complexes dans lesquels la fonction de distribution s'exprime sur l'espace physique et de configuration.

**KEY WORDS:** kinetic theory, Fokker-Planck equation, model reduction, separated representation, Karhunen-Loève decomposition
# 1 INTRODUCTION

In kinetic theory models the molecular conformation, and more in general the fluid microstructure, is described by a distribution function $\psi(x, t, X)$ that represents the probability of finding, at point $x$ in the physical space and time $t$, the fluid microstructure defined by the conformation vector $X$. The evolution of such distribution function is governed by the advection-diffusion Fokker-Planck equation that establishes the conservation balance of $\psi$. The simplest form of such equation is given by:

$$ \frac{d\psi}{dt} = -\frac{\partial}{\partial \mathbf{X}}(\mathbf{X}\psi) + \frac{\partial}{\partial \mathbf{X}}\left(D_r \frac{\partial\psi}{\partial \mathbf{X}}\right) \quad (1) $$

where $d/dt$ represents the material derivative in the physical space.

Since the early 1990's the use of kinetic theory models has developed considerably following the introduction of the CONNFFESSIT method by Öttinger and Laso [1] that is based on the stochastic solution of the Fokker-Planck equation. Since its introduction, it has been successfully applied for solving numerous models, most of them related to dilute polymer solutions. In recent years, kinetic theory of entangled systems such as concentrated polymer solutions and polymer melts, has known major developments that go well beyond the classical reptation tube model developed by Edwards, de Gennes, and Doi. Micro-macro simulations allow to introduce in the macroscopic scale (the one related to the momentum and mass conservation) a constitutive equation based on a microscopic kinetic theory description usually accomplished in the context of stochastic simulations (see [2] for a deep review on this topic).

As claimed by different authors, the control of the statistical noise is a major issue in stochastic simulations. This difficulty does not arise at all in the Fokker-Planck approach. The difficulty, however, is that the Fokker-Planck equation must be solved for the distribution function in both physical and configuration spaces.

When the model involves high dimensional spaces (including physical and conformation spaces and time) standard discretization techniques fail due to exorbitant computation requirements to perform accurate numerical simulations. This difficulty justifies the reduced

number of works that consider the direct solution of the Fokker-Planck equation using standard deterministic discretization techniques, that have been only applied to solve models involving moderate number of dimensions [5-7].

One appealing strategy that allows to alleviate the computational effort is based on the use of reduced approximation bases within an adaptive procedure, which induces the associated reduction of the number of degrees of freedom and then the size of the linear systems to be solved. This strategy was deeply analyzed in [8] and successfully applied for solving FENE molecular models in [9]. The main ideas of this technique will be briefly summarized in Section 2.1.

An important drawback of such an approach is the fact that the approximation functions are defined in the whole domain, and until now, the simplest and most usual form to represent a function is by giving its values at some points in the domain of interest, being its value defined in any other point by interpolation. Some times, the treated model results highly multidimensional, and in this case the possibility of describing functions from their values at the nodes of a mesh (or a grid) of the domain of interest results exorbitant. Some attempts exist concerning the treatment of multidimensional problems. The interested reader can refer to [10] for a review on sparse grids methods involving sparse tensor product spaces, but despite of its optimality, the interpolation is defined in the whole multidimensional domain, and consequently only problems defined in spaces of dimension of the order of tens can be treated. In [11] multidimensional problems are revisited and deeply analyzed, and for this purpose new mathematical entities are introduced to be applied in the numerical treatment of such problems. Thus, the problematic lies to models defined in multidimensional spaces remains still open, and new efforts must be paid to reach significant improvements in the next years.

Some of the most usual kinetic theory models defined from the Fokker-Planck formalism have two important particularities: (i) they can be expressed in a separated form and (ii) the conformation domain is bounded, with the distribution function vanishing on its boundary. In this case, the separated representation and the definition of tensor product approximation spaces work perfectly, and allow to circumvent the difficulties related to the multidimensional charac-
ter of kinetic theory models, as proved in [12]. In [12] we considered the steady state solution of some classes of multidimensional partial differential equations, and in particular those governing the molecular configuration distribution in kinetic theory models of complex fluids associated with dilute polymer solutions (FENE and multi-bead-spring FENE models). In [13] we proposed an extension of that technique to solve multidimensional transient kinetic theory models. An overview about single and multichain FENE models and approaches to their numerical solution had been presented in [14, 15].

We would like to emphasize that we are not looking for a general numerical procedure for solving multidimensional PDE. In this context, the sparse grids or sparse tensor product bases [10] are excellent candidates, but they are not able to treat highly multidimensional problems ($N \gg 20$, $N$ representing that dimension) as claimed in [16]. The technique proposed in our former work [12] and extended in [13] to transient simulations, is, in our opinion, a suitable choice when we are dealing with highly multidimensional parabolic PDE with homogeneous boundary conditions, and despite of the apparent lack of generality, usual molecular descriptions make use of a configuration distribution function defined in a bounded domain that vanishes on its boundary.

## 1.1 SOME EXAMPLES OF KINETIC THEORY MODELS

As the different numerical approaches will be illustrated through the solution of some kinetic theory models, we start by introducing those models.

### Short fiber suspensions

In the case of a short fiber suspension, the configuration distribution function (also known as the orientation distribution function) gives the probability of finding the fibers in a given direction. Obviously, this function depends on the physical coordinates (space and time) as well as on the configuration coordinates, that taking into account the rigid character of the fibers, are defined on the surface of the unit sphere. Thus, we can write $\psi(x, t, \mathbf{p})$, where $x$ defines the position of the fiber center of mass, $t$ the time and $\mathbf{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 \mathbf{p}}(\psi \dot{\mathbf{p}}) + \frac{\partial}{\partial \mathbf{p}}\left(D_r \frac{\partial \psi}{\partial \mathbf{p}}\right) \quad (2) $$

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

$$ \int_{S(o,1)} \psi(\mathbf{x}, t, \mathbf{p}) d\mathbf{p} = 1, \forall \mathbf{x}, \forall t \quad (3) $$

where $S(o, 1)$ denotes the surface of the unit sphere that in this case defines the conformation space. When the fibers are assumed ellipsoidal and when the suspension is dilute enough, the rotation velocity can be obtained from the Fenyery's equation [17]

$$ \dot{\mathbf{p}} = \Omega \mathbf{p} + k \mathbf{D} \mathbf{p} - k (\mathbf{p}^T \mathbf{D} \mathbf{p}) \mathbf{p} \quad (4) $$

where $\Omega$ and $\mathbf{D}$ are the vorticity and the strain rate tensors respectively, associated with the fluid flow undisturbed by the presence of the fiber, and $k$ is a scalar parameter 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} \quad (5) $$

### Dilute polymer solutions: the FENE model

We consider here the dimensionless form of the FENE model [18]

$$ \frac{d\psi}{dt} = -\frac{\partial}{\partial \mathbf{q}}\left[\left(\text{Grad } \mathbf{v}\mathbf{q} - \frac{1}{2}h(q)\mathbf{q}\right)\psi\right] + \frac{1}{2}\frac{\partial}{\partial \mathbf{q}}\left(\frac{\partial\psi}{\partial \mathbf{q}}\right) \quad (6) $$

where $\psi(x, t, q)$ being $q$ the dimensionless connector vector, $q = \|\mathbf{q}\|$, and the function $h$ is given by

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

where $b$ is the dimensionless square of the maximum molecule stretch. Besides, a normality condition is associated with Eq. 6
$$ \int_C \psi(\mathbf{x}, t, \mathbf{q}) d\mathbf{q} = 1, \forall \mathbf{x}, \forall t \qquad (8) $$

where C defines the conformation space. Finally, the link between statistical distribution of dumb-bell configurations and the polymer stress $T_p$ is provided by the Kramer's expression

$$ T_p(\mathbf{x}, t) = \langle h(q) \mathbf{q} \otimes \mathbf{q} \rangle = \int_C \psi(\mathbf{x}, t, \mathbf{q}) h(q) \mathbf{q} \otimes q dq \qquad (9) $$

**Remark 1.1:** The derivative $d\psi/dt$ refers to the material derivative in general. In problems leading to homogeneous solutions this derivative becomes the partial derivative. This is the case for the models considered in this paper, except when for the complex flow involving the coupling between the physical and the conformation coordinates addressed in the last section.

**Remark 1.2:** Equation 6 defines the time evolution of the molecule distribution function. To be solved, an initial distribution must be prescribed. The equilibrium state is defined for a null velocity gradient. Thus, the distribution function becomes isotropic, it depends on the norm of vector $\mathbf{q} = (||\mathbf{q}|| = q)$, and it verifies

$$ o = -\frac{\partial}{\partial q} \left[ \left( -\frac{1}{2} h(q) q \right) \psi_o \right] + \frac{1}{2} \frac{\partial^2 \psi_o}{\partial q^2} \qquad (10) $$

that results in

$$ \psi_o(q) = \frac{h(q)^{-b/2}}{\int h(q)^{-b/2} dq} \qquad (11) $$

Thus, a natural choice for the initial condition consists in $\psi(\mathbf{x}, t = 0, \mathbf{q}) = \psi_o(q)$.

### Dilute polymer solutions: the MBS model

The multi-bead-spring (MBS) chain consists of S + 1 beads connected by S springs. The bead serves as an interaction point with the solvent and the spring contains the local stiffness information depending on local stretching (see [18] for more details). The associated dimensionless Fokker-Planck writes:

$$ \frac{d\psi}{dt} + \frac{\partial}{\partial \mathbf{q}} \left[ (\text{Grad} \mathbf{vq} - \mathbf{Ag}(\mathbf{q})) \psi \right] - \frac{1}{4} \frac{\partial}{\partial \mathbf{q}} \left( \mathbf{A} \frac{\partial \psi}{\partial \mathbf{q}} \right) = o \qquad (12) $$

where $\mathbf{q}$ is the configuration vector $\mathbf{q} = (\mathbf{q}_1, \mathbf{q}_2, \dots, \mathbf{q}_S)$ and the components of $\mathbf{g}$ are given by:

$$ \mathbf{g}_k = \frac{1}{4} \mathbf{q}_k h(\mathbf{q}_k), \quad \mathbf{q}_k = \|\mathbf{q}_k\| \qquad (13) $$

being the Rouse matrix

$$ \mathbf{A}_{kl} = \begin{cases} 2 \times \mathbf{l} & \text{if } k = l \\ -\mathbf{l} & \text{if } k = l \pm 1 \\ \mathbf{o} \times \mathbf{l} & \text{otherwise} \end{cases} $$

being $\mathbf{l}$ the unit matrix whose dimension coincides with the dimension of the physical space where vectors $\mathbf{q}_k$ are defined.

### Entangled polymers: the Doi-Edwards model

In most polymer processing operations such as injection molding, film blowing and extrusion, the polymers are in the molten state. A widely applied class of molecular-based models for concentrated polymer solutions and melts relies on the reptation notion.

The key idea of this model is the use of the reptation mechanism introduced by de Gennes [19] to a tube (along which the molecule can move) in order to describe the viscoelastic behaviour of entangled polymers. The molecule is described as sliding or reptating through a tube whose contours are defined by the locus of entanglements with neighbouring molecules. The motion of a molecular chain in any other direction than the one defined by the tube axis is strongly restricted, except at both tube ends, where it can move in any possible direction. The tube moves itself due to two mechanisms: (i) by means of the motion of the central chain itself, which partially leaves its original tube, for extending it in other directions, and (ii) by the fluctuation induced by the motions of the neighbour chains defining the tube lateral border. In addition to the reptation mechanism, the Doi-Edwards model [20] assumes affine tube deformation induced by the macroscopic flow, but neglects other phenomena like the stretch of the chain and the Convective Constraint Release. Within this reptation picture and these assumptions, the dynamics of a single segment is given by:

$$ \frac{du}{dt} = (\mathbf{l} - \mathbf{u} \otimes \mathbf{u}) \cdot (\text{Grad} \mathbf{v} \cdot \mathbf{u}) \qquad (14) $$
where **u** is the unit vector describing the orientation of the tube segment, *s* the contour length coordinate ($s=0$ and $s=1$ represent the chain ends), Grad **v** the velocity gradient, and **I** the unit tensor. The distribution function $\psi$ is such that $\psi(\mathbf{x}, t, \mathbf{u}, s)duds$ represents the joint probability that at time $t$ and position $\mathbf{x}$ a tube segment has an orientation in the interval $[\mathbf{u}, \mathbf{u} + d\mathbf{u}]$ and contains the chain segment labelled in the interval $[s, s + ds]$. Thus, the configuration space results $\Omega = S(0, 1) \times [0, 1]$, where $S(0, 1)$ is the surface of the unit sphere centered at the origin. The Fokker-Planck equation related to the Doi-Edwards model is the convection-diffusion equation that represents the conservation balance of the distribution function:

$$ \frac{d\psi}{dt} + \frac{\partial}{\partial \mathbf{u}} \left( \left[ (\mathbf{I} - \mathbf{u} \otimes \mathbf{u}) \cdot (\text{Grad} \mathbf{v} \cdot \mathbf{u}) \right] \psi \right) - \frac{1}{\pi^2 \tau_d} \left( \frac{\partial^2 \psi}{\partial s^2} \right) = \mathbf{0} \quad (15) $$

where $d/dt$ represents the material derivative, $\tau_d$ is the disengagement time, namely the characteristic time for a chain to come out the tube by reptation. We can define the diffusion coefficient related to the $s$-coordinate as $D_s = 1/(\pi^2\tau_d)$. To solve the Fokker-Planck equation, one needs to prescribe appropriate boundary conditions at the tube borders $s=0$ and $s=1$ (the orientation coordinate being defined on the unit sphere does not require any boundary condition). Sometimes, as previously indicated, an isotropic orientation distribution is assumed at both ends, which reads:

$$ \psi(\mathbf{x}, t, \mathbf{u}, s = 0) = \psi(\mathbf{x}, t, \mathbf{u}, s = 1) = \frac{1}{4\pi} \quad (16) $$

However, other boundary conditions can be also applied as described in Sect. 6.2 of [16], where flow-induced alignment of chain ends had been considered, with a strong effect on dynamic viscosities. Knowing the distribution function, the stress can be computed from:

$$ T_p(\mathbf{x}, t) = G \int_0^1 \int_{S(0,1)} \mathbf{u} \otimes \mathbf{u} \psi(\mathbf{u}, s, \mathbf{x}, t) duds \quad (17) $$

where $G$ is an elastic modulus and $T_p$ is the polymer stress tensor.

## Advanced reptation models

Recognizing the large discrepancy between the Doi and Edwards model and the experimental

data many attempts of modifying the original model have been recently made: avoiding independent alignment, introducing double reptation, chain stretching, convective constraint release as well as anisotropic tube cross sections (see [21] and the references therein). In the present section we focus on the model proposed in [21] that will be referred as the Öttinger's model. This model introduces a further variable $\lambda$ to account for the chain stretching, which is defined as the ratio $\lambda = L/L_o$, where $L$ and $L_o$ denotes the contour length and the equilibrium contour length of the primitive chain.

The equation governing the evolution of the contour length, involving stretching and relaxation, reads:

$$ \frac{d\lambda}{dt} = \dot{\lambda} = \dot{\lambda}_{\text{convect}} + \dot{\lambda}_{\text{dissip}} \quad (18) $$

where $\dot{\lambda}_{\text{convect}}$ and $\dot{\lambda}_{\text{dissip}}$ denote the convective and dissipation contributions respectively. The first contribution is defined by:

$$ \dot{\lambda}_{\text{convect}} = \lambda (\text{Grad} \mathbf{v}^T : \bar{\mathbf{t}}) \quad (19) $$

where $\bar{\mathbf{t}}$ is the symmetric second-moment orientation tensor:

$$ \bar{\mathbf{t}} = \int_0^1 \int_{S(0,1)} (\mathbf{u} \otimes \mathbf{u}) \psi(\mathbf{u}, s, \mathbf{x}, t) duds \quad (20) $$

The dissipative contribution reads:

$$ \dot{\lambda}_{\text{dissip}} = - \frac{1}{\tau_s} \frac{c(\lambda)}{3Z} (\lambda - 1) \quad (21) $$

where $Z$ is the number of entanglements per chain at equilibrium, $\tau_s$ is the characteristic stretching time and $c(\lambda)$ is an effective spring coefficient which models the elastic nature of the polymer molecules:

$$ c(\lambda) = \frac{3Z\lambda_{\text{max}}^2 (\lambda + 1)}{\lambda(\lambda_{\text{max}}^2 - \lambda^2)} \quad (22) $$

where $\lambda_{\text{max}}$ is the maximum possible stretching ratio of the contour length. The balance equation for the configurational distribution function $\psi(\mathbf{x}, \mathbf{u}, \mathbf{s})$ takes the form:
$$ \frac{d\psi}{dt} = -\frac{\partial}{\partial \mathbf{u}}\left(\left(\mathbf{I} - \mathbf{u}\otimes\mathbf{u}\right) \cdot \left(\text{Grad}\cdot\mathbf{u}\right)\right)\psi \\ -\frac{\partial}{\partial s}(\dot{S}_{\text{tot}}\psi) - \frac{\dot{\lambda}_{\text{dissip}}}{\lambda}\psi \\ + \frac{1}{\pi^2\tau_d}\frac{\partial^2\psi}{\partial s^2} + \frac{\partial}{\partial \mathbf{u}}\left(D_u\frac{\partial\psi}{\partial \mathbf{u}}\right) \qquad (23) $$

where $1/(\pi^2\tau_d)$ is the diffusion coefficient related to s-coordinate, $D_u$ the diffusion coefficient related to u-coordinate, being the diffusion orientation term defined on the unit surface, I the identity tensor and $\dot{S}_{tot}$ the drift velocity for s which states that there is only a re-scaling of the position label s for the tube segment when the chain relaxes in the tube:

$$ \dot{S}_{tot} = -\frac{1}{\lambda}\left(s - \frac{1}{2}\right)\dot{\lambda}_{dissip} \qquad (24) $$

The creation/destruction term - $(\lambda_{dissip}/\lambda)\psi$ compensates for configurations lost or gained at the boundaries, because the equality

$$ \bar{\mathbf{t}} = \int_{0}^{1} \int_{S(o,1)} \dot{S}_{tot} \psi du ds = \int_{0}^{1} \int_{S(o,1)} \frac{\dot{\lambda}_{dissip}}{\lambda} \psi du ds \qquad (25) $$

The terms involving second-order derivatives are of irreversible nature and express the erratic representation motion along the chain contour (second derivative with respect to the s-coordinate) and the constraint release (second-order derivative with respect to u whose associated diffusion coefficient is noted by $D_u$).

The following expression for $D_u$ is proposed in [21]:

$$ D_u = \frac{1}{6} \left[ \delta_1 \frac{1}{\tau_d} - \delta_2 \frac{\dot{\lambda}_{dissip}}{\lambda} H\left(\frac{-\dot{\lambda}_{dissip}}{\lambda}\right) \right] \qquad (26) $$

Where $H(x)$ is the Heaviside step function ($H(x) = 1$ for $x \ge 0$, $H(x) = 0$ for $x < 0$). The quantities $1/\tau_d$ and $\lambda_{dissip}/\lambda$ determine the rates of constraint release due to the loss of entanglements caused by the reptation motion and the retraction of the surrounding chains, respectively. Hence, the $\delta_1$-term can be interpreted as representing “double reptation”, while the $\delta_2$-term can be regarded as the convective constraint release mechanism. At

the surrounding chain ends, we impose random orientation by specifying the distribution:

$$ \psi(\mathbf{x}, t, \mathbf{u}, s = 0) = \psi(\mathbf{x}, t, \mathbf{u}, s = 1) = \frac{1}{4\pi} \qquad (27) $$

The polymer stress $T_p$ is the obtained from

$$ T_p = G \left( 1 + \lambda_{max}^2 \frac{(\lambda^2 - 1)}{\lambda_{max}^2 - \lambda^2} \right) \bar{\mathbf{t}} \qquad (28) $$

where $G$ is an elastic modulus and $\bar{\mathbf{t}}$ is the symmetric second-moment orientation tensor defined in Eq. 20.

Remark 1.3: Network-based entangled polymer models can be also defined in the kinetic theory framework by using some specific Fokker-Planck equations whose expressions are close to the ones just presented.

## 2 REDUCED ORDER MODELLING

The above introduced Fokker-Planck equations will be solved by using two different reduced model techniques, the one based on the proper orthogonal decomposition (also known as Karhunen-Loève decomposition), which is an appealing choice in the case of kinetic theory models involving a reduced number of conformation coordinates, and the other one based on a separated representation, which is specially addressed to solve problems in highly dimensional spaces. We briefly summarize both strategies in the next sections.

### 2.1 INTRODUCTION: THE KARHUNEN-LOEVE DECOMPOSITION

We assume that the evolution of a certain field $u(\mathbf{x}, t)$ is known (being its evolution governed by a PDE). For practical purposes, this field is expressed in a discrete form, that is, it is defined at the $N_p$ nodes of a spatial mesh $\mathbf{x}_i$ and for some times $t^p$, i.e. $u(\mathbf{x}_i, t^p) \equiv u_i^p$ (we assume $N_p$ nodes and $P$ time steps). The main idea of the Karhunen-Loève (KL) decomposition is how to obtain the most typical or characteristic structure $\phi(\mathbf{x})$ among these $u_i^p(\mathbf{x}_i)$, for all $p$. For this purpose, we look for a function $\phi(\mathbf{x})$ that maximizes $\beta$ defined by

$$ \beta = \frac{\sum_{p=1}^{p=N} \left[ \sum_{i=1}^{i=N_p} \phi(\mathbf{x}_i) u_i^p(\mathbf{x}_i) \right]^2}{\sum_{i=1}^{i=N} (\phi(\mathbf{x}_i))^2} \qquad (29) $$
The maximization ($\delta\beta = 0$) leads to:

$$ \sum_{p=1}^{p=N_p} \left[ \left( \sum_{i=1}^{i=N_p} \tilde{\phi}(\mathbf{x}_i) u^p(\mathbf{x}_i) \right) \left( \sum_{j=1}^{j=N_p} \phi(\mathbf{x}_j) u^p(\mathbf{x}_j) \right) \right] = \beta \sum_{i=1}^{i=N_p} \tilde{\phi}(\mathbf{x}_i) \phi(\mathbf{x}_i); \forall \tilde{\phi} \qquad (30) $$

which can be rewritten in the form

$$ \sum_{i=1}^{i=N_p} \left\{ \sum_{j=1}^{j=N_p} \left[ \sum_{p=1}^{p=N_p} u^p(\mathbf{x}_i) u^n(\mathbf{x}_j) \phi(\mathbf{x}_j) \right] \tilde{\phi}(\mathbf{x}_i) \right\} = \beta \sum_{i=1}^{i=N_p} \tilde{\phi}(\mathbf{x}_i) \phi(\mathbf{x}_i); \forall \tilde{\phi} \qquad (31) $$

Defining the vector $\Phi$ such that its i-component is $\phi(\mathbf{x}_i)$, Eq. (31) takes the following matrix form

$$ \tilde{\Phi}^T \mathbf{k} \Phi = \beta \tilde{\Phi}^T \phi; \forall \tilde{\Phi} \Rightarrow \mathbf{k} \Phi = \beta \Phi \qquad (32) $$

where the two points correlation matrix is given by

$$ k_{ij} = \sum_{p=1}^{p=N_p} u^p(\mathbf{x}_i) u^n(\mathbf{x}_j) \Leftrightarrow \mathbf{k} = \sum_{p=1}^{p=N_p} u^p(\mathbf{u}^p)^T \qquad (33) $$

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

$$ \mathbf{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} \qquad (34) $$

it is direct to verify that the matrix $\mathbf{k}$ in Eqs. 32 and 33 results

$$ \mathbf{k} = \mathbf{Q}\mathbf{Q}^T \qquad (35) $$

Thus, the functions defining the most characteristic structure of $u^p(\mathbf{x})$ are the eigenfunctions $\phi_k(\mathbf{x}) \equiv \Phi_k$ associated with the highest eigenvalues.

### A posteriori reduced modelling

If some direct simulations have been carried out, we can determine $u(\mathbf{x}_i, t^p) \equiv u_i^p$, for all $i \in [1, ..., N_p]$ for all $p \in [1, ..., P]$ and from these simulations the $r$ eigenvectors related to the $r$-highest eigenvalues $\phi_k^T = [\phi_k(\mathbf{x}_i), ..., \phi_k(\mathbf{x}_{N_p})]$ for all $k \in [1, ..., r]$ (with $r \ll N_p$). Now, we can try to use these eigenfunctions for approximating the solution of either a problem slightly different to the one that has served to define $u(\mathbf{x}_i, t^p)$, or the problem under consideration at higher times $t \gg t^p$. For this purpose we need to define the matrix $\mathbf{B}$

$$ \mathbf{B} = \begin{pmatrix} \phi_1(\mathbf{x}_1) & \phi_2(\mathbf{x}_1) & \cdots & \phi_r(\mathbf{x}_1) \\ \phi_1(\mathbf{x}_2) & \phi_2(\mathbf{x}_2) & \cdots & \phi_r(\mathbf{x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_1(\mathbf{x}_{N_p}) & \phi_2(\mathbf{x}_{N_p}) & \cdots & \phi_r(\mathbf{x}_{N_p}) \end{pmatrix} \qquad (36) $$

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

$$ \mathbf{K} \mathbf{U}^p = \mathbf{F}^{p-1} \qquad (37) $$

where $\mathbf{F}^{p-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

$$ \mathbf{U}^p = \sum_{i=1}^{i=r} \Phi_i \alpha_i^p = \mathbf{B} \alpha^p \qquad (38) $$

implying

$$ \mathbf{K} \mathbf{B} \alpha^p = \mathbf{F}^{p-1} \qquad (39) $$

Multiplying both terms by $\mathbf{B}^T$ it results

$$ \mathbf{B}^T \mathbf{K} \mathbf{B} \alpha^p = \mathbf{B}^T \mathbf{F}^{p-1} \qquad (40) $$

which proves that the final linear system is of small size, i.e. the dimensions of $\mathbf{B}^T \mathbf{K} \mathbf{B}$ are $r \times r$, with $r \ll N_p$, and the dimensions of both $\alpha$ and $\mathbf{B}^T \mathbf{F}$ are $r \times 1$.

**Remark 2.1:** Equation 40 can be also derived introducing the approximation (38) into the PDE Galerkin form.

**Remark 2.2:** In the case of non-linear models $\mathbf{K}$ represents the linearized matrix.
## Adaptivity via an "a priori" model reduction

In order to compute reduced model solutions without an a priori knowledge, Ryckelynck proposed in [22] to start with a low order approximation basis, using some simple functions (e.g. the initial condition in transient problems) or using the eigenvectors of a similar problem previously solved. Now, we compute S time steps of the evolution problem using the reduced model (Eq. 40) without changing the approximation basis, and then compute the residual. If the norm of the residual is small enough, then we can compute other S time steps using the same reduced approximation basis that involves a reduced number of degrees of freedom. In the other case, when the residual norm is greater than a threshold value, we proceed to enrich the approximation basis by adding the residual and some Krylov subspaces generated by the residual to the former approximation basis, and then the previous S time steps are recomputed. The enrichment procedure continues until convergence. When the convergence is reached, and before proceeding with the next S time steps, a Karhunen-Loève decomposition is performed in the whole computed time interval in order to extract only the significant approximation functions. Thus, although the enrichment increases the size of the reduced approximation basis, the Karhunen-Loève decomposition reduces it, and in general the number of approximation functions remains approximatively constant during the entire simulation. See [8] for a deep description of this procedure and [9] for a detailed description of its application for solving the FENE model summarized in Section 1.1.

## 2.2 Separated representation and tensor product approximation spaces

The Fokker-Planck equation being a parabolic PDE, we consider in this section, for the sake of simplicity, a similar equation but involving simpler notation, as it is the case of the transient Poisson's problem defined in a N-dimensional space, given by:

$$ \frac{\partial u}{\partial t} - \Delta u = f(x_1, x_2, \dots, x_N) \quad (41) $$

where u is assumed a scalar function depending on space and time, that is $u(x, t)$ or $u(x, t) = u(x_1,$

$x_2, \dots, x_N, t)$. Equation 41 is assumed to be defined in the domain $\Omega = \Omega_x \times \Omega_t = ]-L, +L[^N X]_0, t_{max}]$ and u is assumed vanishing on the space boundary, i.e. $u(x \in d\Omega_x, t) = 0$ as well as at the initial time, i.e. $u(x, t = 0) = 0$.

**Remark 2.3:** Most of kinetic theory models are defined in a bounded domain, with the associated distribution function vanishing on its boundary. Despite the fact that these models are subjected to a non zero initial condition, it is direct to prove that the subtraction of this initial condition to the original distribution function leads to another function, verifying a slightly different Fokker-Planck equation, which vanishes on the domain boundary as well as at the initial time.

The problem solution is assumed in the form:

$$ u(x_1, x_2, \dots, x_N, t) = \sum_{j=1}^{\infty} \alpha_j \left[ \left( \prod_{k=1}^{N} F_{kj}(x_k) \right) F_{(N+1)j}(t) \right] \quad (42) $$

where $F_{kj}$ is the $j^{th}$ basis function which only depends on the $k^{th}$ space coordinate, the time being no more than other coordinate. The construction of such solution (Eq. 42) consists of an iteration procedure involving at each iteration n two steps:

1.) Projection of the solution in a discrete basis. If we assume the functions $F_{kj}$ (for all $j \in [1, \dots, n]$ for all $k \in [1, \dots, N+1]$) known (verifying the boundary and initial conditions), the coefficients $\alpha_j, j \in [1, \dots, n]$, can be computed by introducing the approximation of u (Eq. 42) into the Galerkin variational formulation associated with Eq. 41 which results in a linear system of size $n \times n$.

2.) Enrichment of the approximation basis. From the alpha coefficients just computed the approximation basis can be enriched by adding the new function $\Pi^{k=N}_{k=1} F_{k(n+1)}(x_k) F_{(N+1)(n+1)}(t)$. For this purpose we solve the Galerkin variational formulation related to Eq. 41 using the approximation of given by:

$$ u(x_1, x_2, \dots, x_N, t) = \sum_{j=1}^{n} \alpha_j \left( \prod_{k=1}^{N} F_{kj}(x_k) \right) F_{(N+1)j}(t) + \left( \prod_{k=1}^{N} R_k(x_k) \right) R_{(N+1)}(t) \quad (43) $$

which results in a non-linear problem of size $\sum_{i=1}^{N+1} n_i$ being $n_i$ the number of nodes used to approximate the function $R_i$. Functions $F_{k(n+1)}$ are
<figure><img src="image_2.png" /><figcaption>Figure 1 (left): Evolution of the steady distribution function defined on the unit surface along the chain arc-length for $D_s = 1$ and $We = 1$: Shear flow (left) and longitudinal flow (right).</figcaption></figure>

<caption>Figure 2: Approximation functions: (top-left) most significant approximation function, (top-right) second most significant function, (bottom-left) third most significant function, and (bottom-right) fourth most significant approximation function.</caption>

finally obtained by normalizing the functions $R_1, R_2, ..., R_{N+1}$.

The interested reader can refer to [12] and [13] for a detailed description of this technique for the solution of the MBS-FENE model summarized in Section 1.1 in the steady and transient cases respectively.

### 3 NUMERICAL EXAMPLES OF KINETIC THEORY MODEL REDUCTION

As the application of the described reduction techniques in the framework of dilute polymer flows was considered in some of our former works, as previously indicated, we will focus in this section on the solution of more complex models concerning entangled polymers. The main difficulties encountered in the problems addressed here concern the complexity of the configuration space, the non-linearity of some models as well as the complex flow simulations coupling physical and configuration spaces.

#### 3.1 REDUCED KARHUNEN-LOEVE APPROXIMATION OF THE DOI-EDWRADS MODEL

In order to analyze the capabilities of the technique described in Section 2.1 we consider the 3D Doi-Edwards model (Eq. 15) in some simple

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

rheometric flows in which the distribution function can be expressed in the form $\psi(\mathbf{u}, s, t)$. The reduction technique based on the use of the Karhunen-Loève decomposition allows to define a reduced approximation basis that is able to represent the entire time evolution of the distribution function using only 8 approximation functions. The evolution of the steady-state distribution function along the curvilinear coordinate s is depicted in Figure 1 in a shear and in an elongational flow, both characterized by a Weissenberg number $We = 1$ (the diffusion coefficient being $D_s = 1$). In both cases only 8 degrees of freedom are needed, instead of the 642 approximation functions related to the 642 nodes defining the coarse finite element mesh on the unit surface. Thus, at each time step one needs to perform the inversion of a matrix of size 8 x 8. These four most significant approximation functions are depicted in Figure 2.

#### 3.2 SEPARATED REPRESENTATION OF REPTATION MODELS OF ENTANGLED POLYMERS

In this section we analyze the use of a variable separation methodology that allows to solve the problem in the multidimensional space by decoupling the different subspaces. Moreover, the approximation functions associated with each subspace are computed during the simulation in order to keep a reduced number of degrees of freedom.

First, we consider again the Fokker Planck equation of the Doi-Edwards model (Eq. 15) whose solution in some simple rheometric flows is assumed not depending on the physical coordinates:

$$ \psi(\mathbf{u}, s, t) = \sum_{j=1}^{\infty} \alpha_j F_j(\mathbf{u}) G_j(s) H_j(t) + \Xi \quad (44) $$

where the conformation is defined at time t from: $(\mathbf{u}, s, t) \in \Omega = \Omega_u \times \Omega_s \times \Omega_t = S(0, 1) \times [0, 1] \times [0, t_{max}]$,
<figure><img src="image_4.png" /><figcaption>Figure 3 (above): Comparison of the Doi-Edwards model (left) and the Öttinger one (right) in a steady state 2D shear flow characterized by We = 1 and δ₁ = δ₂ = 1/λ: (top) distribution functions, (middle) reduced approximation bases related to the orientation coordinate, and (bottom) reduced approximation bases related to the contour length coordinate.</figcaption></figure>

<caption>Figure 4: Comparison of the Doi-Edwards model (left) and the Öttinger one (right) in a steady state 2D shear flow characterized by We = 100 and δ₁ = δ₂ = 1/λ: (top) distribution functions; (middle) reduced approximation bases related to the orientation coordinate and (bottom) reduced approximation bases related to the contour length coordinate.</caption>

and the term Ξ allows to enforce homogeneous boundary conditions to the separated representation. Because the boundary conditions

$$ \psi(\mathbf{u}, s = 0) = \psi(\mathbf{u}, s = 1) = \frac{1}{4\pi} \qquad (45) $$

a natural choice consists in Ξ = 1/4π.

In the transient case, as the initial condition $\psi(\mathbf{u}, s, t = 0)$ also verifies Eq. 45, one could consider Ξ = $\psi(\mathbf{u}, s, t = 0)$ that can be expressed in a separated form by using the singular value decomposition or the alternating least squares technique. In the numerical examples that follow, we focuss in the steady state solution and then the choice Ξ = 1/4π is considered.

The same strategy can be applied for solving the advanced reptation model including double reptation as well as convective constraint release just described, by applying an appropriate technique for accounting for the non-linearity. In the simulations that follow, the simplest technique, based on a fixed point strategy, is considered. We consider both the Doi-Edwards and the Öttinger models in simple 2D and 3D shear flows. The computed solutions for two

different Weissenberg numbers (We = 10 and We = 100) are then compared.

The results computed in the 2D case are depicted in Figures 3 and 4. We can notice that the parabolic profile of the distribution function along the contour length coordinate becomes for the Öttinger model the more and more planar as the Weissenberg number increases, proving that the relaxation is localized in the neighborhood of the chain ends.

The 3D case is analyzed. Figure 5 depicts the results obtained by setting $δ₁ = 1/λ$, $λ_{max} ≈ 10$ and $Z = 20 ≫ 1$. In this figure the distribution functions related to the Öttinger model are depicted at the middle chain section $s = 0.5$ for different Weissenberg numbers. Moreover, we included in this figure the profile of the distribution function along the contour length coordinate for different values of the conformation vector $\mathbf{u}$, proving that these profiles, parabolic at low Weissenberg numbers, become highly planar by increasing it. The computed dimensionless shear stresses for different choices of the parameters $δ₁$ and $δ₂$ are in perfect agreement with the ones reported by Fang et al. (see Figure 8 in [14]). Figures 6 and 7
<figure><img src="image_5.png" /><figcaption>**Figure 6:** Approximation functions involved in the separated representation of the distribution function related to the Öttinger model in a 3D steady shear flow characterized by $We = 1$.</figcaption></figure>

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

depict the different approximation functions involved in the separated representation of the distribution function for both $We = 10$ and $We = 100$, respectively.

### 3.3 COMPLEX FLOW SIMULATION: COUPLING PHYSICAL AND CONFORMATION SPACES

In complex flows simulation two difficulties appear: (i) the first difficulty is related to the necessity of solving the macroscopic flow kinematics that is coupled with the microscopic fluid description, and (ii) the microscopic fluid description accomplished in the kinetic theory framework introduces now a distribution function that depends on the physical and the conformation coordinates.

From now on, we will assume, as widely employed, that the first difficulty is circumvented considering a decoupled solution in which the flow kinematics at a certain iteration (or a certain time step in transient simulations) is accomplished by using the microscopic fluid description at the previous iteration (or time step), and so on. The second difficulty is also well known. Now, the distribution function is defined in a multidimensional space including the physical and the conformation coordinates. This difficulty motivated the development of stochastic simulations. However, the separated representation previously described, and its associated reduced tensor product approximation basis, makes possible the solution of the deterministic multidimensional Fokker-Plank equation. This section focusses in this new potentiality.
<figure><img src="image_7.png" /></figure>

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

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

<figure><img src="image_10.png" /><figcaption>Figure 7 (above): Approximation functions involved in the separated representation of the distribution function related to the Öttinger model in a 3D steady shear flow characterized by We = 100.</figcaption></figure>

<caption>Figure 8 (left below): Fiber orientation distribution function along a circular trajectory.</caption>

<figure><img src="image_11.png" /><figcaption>Figure 9 (right below): Comparison of the exact and computed fiber orientation distributions at a point located in the circular trajectory.</figcaption></figure>

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

We consider the kinetic theory model of a short fiber suspension summarized in Section 1.1, defined by the Fokker-Planck equation

$$ \frac{d\psi}{dt} = -\frac{\partial}{\partial \mathbf{p}}(\dot{\mathbf{p}}\psi) + \frac{\partial}{\partial \mathbf{p}}\left(D_r \frac{\partial\psi}{\partial \mathbf{p}}\right) \quad (46) $$

As the equation is an advection equation in the physical domain, the solution at a certain point only depends on its upstream solution. Thus, two possibilities exist, the first one based on the solution of that equation along a certain number of trajectories, from which the solution could be interpolated everywhere; and the other one based on the solution in the whole physical domain using an appropriate stabilized discretization technique (to account for its hyperbolic character). In that follows we investigate the first strategy, the second one being under development.

We consider a 2D steady recirculating shear flow whose kinematics is defined by:

$$ \mathbf{v} = \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} -y\sqrt{x^2 + y^2} \\ x\sqrt{x^2 + y^2} \end{pmatrix} \quad (47) $$

and we look for the steady solution of Eq. 46 along the circular streamline defined by $(x^2 + y^2)^{0.5} = 0.5$, being the diffusion coefficient $D_r = 0.1$. The distribution function can be expressed as

$$ \psi(s, \varphi) = \sum_{j=1}^{\infty} \alpha_j F_j(\varphi) G_j(s) \quad (48) $$

where $\varphi$ is the angular coordinate defining the planar fiber orientation and $s$ the curvilinear coordinate along the circular streamline. Figure 8 depicts the resulting distribution function. If we compare at a certain point on the streamline the computed distribution with the exact one (that can be easily determined [23]), a perfect agreement is noticed as proved in Figure 9.

Now, the same problem is solved assuming a 3D fiber orientation distribution which is
<figure><img src="image_13.png" /><figcaption>Figure 10 (above): Steady 3D fiber orientation distribution at some locations on the circular trajectory.</figcaption></figure>

<caption>Figure 11: 3D orientation approximation functions $F_j(\mathbf{p})$.</caption>

defined on the unit surface. Thus, the separated representation reads:

$$ \psi(s, \mathbf{p}) = \sum_{j=1}^{\infty} \alpha_j F_j(\mathbf{p}) G_j(s) \qquad (49) $$

Figure 10 depicts the computed results and Figure 11 the 9 orientation approximation functions used for computing the distribution evolution along the whole circular trajectory. The computed solution is in good agreement with the ones computed by using the standard finite element technique.

## 4 CONCLUSIONS AND PERSPECTIVES

In this paper we explored the ability of some recent numerical strategies for solving steady or transient multidimensional partial differential equations as the ones encountered in the kinetic theory description of complex fluids. The technique based on the use of reduced approximation bases constructed from the application of the Karhunen-Loève decomposition can be successfully applied to simulate kinetic theory models defined in spaces of moderate dimension, but it fails in highly dimensional spaces. In highly dimensional problems, the technique based on a separated representation of the solution seems to be an excellent candidate for treating this kind of complex models, as proved in this paper throughout several examples.

<figure><img src="image_14.png" /></figure>
The main challenge in the simulation of complex flows is the coupling between the physical and the conformation discretization spaces. The separated representation of the advective stabilization terms requires further developments. On the other hand, fully coupled models, coupling the macroscopic flow kinematics and the microscopic fluid description, deserve additional works.

## ACKNOWLEDGEMENTS

We thank R. Keunings for its participation in our former works concerning the application of model reduction techniques for solving kinetic theory models, as well as for its helpful suggestions.

## REFERENCES

[1] Öttinger HC, Laso M: Smart polymers in finite element calculation, Int. Congr. on Rheology, Brussel, Belguim, 1992.

[2] Keunings R: Micro-macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory, Rheol. Rev., DM Binding and K Walters (Eds.), Brit. Soc. Rheol. (2004) 67-98.

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

[4] Chaubal CV, Srinivasan A, Egecioglu O, Leal LG: Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems, J. Non-Newtonian Fluid Mech. 70 (1997) 125-154.

[5] Lozinski A, Chauviere C, Fang J, Owens RG: Fokker-Planck simulations of fast flows of melts and concentrated polymer solutions in complex geometries, J. Rheol. 47 (2003) 535-561.

[6] Lozinski A, Owens RG, Fang J: A Fokker-Planck-based numerical method for modelling strongly non-homogeneous flows of dilute polymeric solutions, Journal of Non-Newtonian Fluid Mechanics 122 (2004) 273-286.

[7] Ammar A, Chinesta F, Ryckelynck D: Deterministic particle approach of multi-bead-spring polymer models, European Journal of Computational Mechanics, in press.

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

[9] Ammar A, Ryckelynck D, Chinesta F, Keunings R: On the Reduction of Kinetic Theory Models Related to Finitely Extensible Dumbbells, J. Non-Newtonian Fluid Mech. 134 (2006) 136-147.

[10] Bungartz HJ, Griebel M: Sparse grids, Acta Num. 13 (2004) 1-123.

[11] Beylkin G, Mohlenkamp M: Algorithms for numerical analysis in high dimensions, SIAM J. Sci. Comp., 26 (2005) 2133-2159.

[12] Ammar A, Mokdad B, Chinesta F, Keunings R: A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids. Journal of Non-Newtonian Fluid Mechanics, in press.

[13] Ammar A, Mokdad B, Chinesta F: A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids. Part 2: Transient simulation using space-time separated representations. Journal of Non-Newtonian Fluid Mechanics, submitted.

[14] Kröger M: Simple models for complex nonequilibrium fluids, Phys. Rep. 390 (2004) 453-551.

[15] Kröger M: Models for polymeric and anisotropic liquids. Springer, Berlin, 2005.

[16] Achdou Y, Pironneau O: Computational methods for option pricing, Siam Frontiers in Applied Mathematics, 2005.

[17] Jeffery GB: The motion of ellipsoidal particles immersed in a viscous fluid, Proc. R. Soc. London A 102 (1922) 1922.

[18] Bird RB, Curtiss CF, Armstrong RC, Hassager O: Dynamics of polymeric liquids, John Wiley & Sons, 2: Kinetic Theory, 1987.

[19] de Gennes PG: Reptation of a polymer chain in the presence of fixed obstacles, J. Chem. Phys. 55 (1971) 572-579.

[20] Doi M, Edwards SF: The theory of polymer dynamics, Oxford Science Publications, 1978.

[21] Fang J, Kröger M, Öttinger HC: A thermodynamically admissible reptation model for fast flows of entangled polymers. II. Model predictions for shear and extensional flows, J. Rheol. 40 (2000) 1293-1318.

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

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