# Reduced Numerical Modeling of Flows Involving Liquid-Crystalline Polymeres

Amine Ammar, Etienne Prulière, Francisco Chinesta, Manuel Laso

► To cite this version:

Amine Ammar, Etienne Prulière, Francisco Chinesta, Manuel Laso. Reduced Numerical Modeling of Flows Involving Liquid-Crystalline Polymeres. Journal of Non-Newtonian Fluid Mechanics, 2009, 160 (2 3), pp.140-156. 10.1016/j.jnnfm.2009.03.013. hal-01007243

HAL Id: hal-01007243

https://hal.science/hal-01007243

Submitted on 8 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
Reduced numerical modeling of flows involving liquid-crystalline polymers

A. Ammarª, E. Pruliereª, F. Chinestaᵇ, M. Lasoᶜ

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

ᵇ EADS Corporate Foundation International Chair, Ecole Centrale of Nantes, 1 rue de la Noe, BP 92101, F-44321 Nantes cedex 3, France

ᶜ Laboratory of Non-Metallic Materials, Department of Chemical Engineering and Institute for Optoelectronics and Microsystems (ISOM) ETSII/UPM, Jose Gutierrez Abascal, 2, 28006 Madrid, Spain

Kinetic theory models involving the Fokker-Planck equation can be accurately discretized using a mesh support (Finite Elements, Finite Differences, Finite Volumes, Spectral Techniques, ...). However, these techniques involve a high number of approximation functions. In the finite element framework, widely used in complex flow simulations, each approximation function has only local support and is related to a node that defines the associated degree of freedom. In the technique proposed here, a reduced approximation basis is constructed. The new shape functions have extended support and 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 very 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 paper analyzes the application of model reduction to the simulation of non-linear kinetic theory models involving complex behaviors, such as those coming from stability analysis, complex geometries and coupled models. We apply our model reduction approach to the Doi's classical constitutive equation for viscoelasticity of liquid-crystalline polymer.

Keywords:

Complex fluids

Kinetic theory

Model reduction

Liquid-crystalline polymers

Numerical modeling

# 1. Introduction

The purpose of this work is to present a robust and simplified computational scheme to carry out micro-macro-simulations of molecular behavior within the kinetic theory framework. This framework has been used traditionally for deriving constitutive equations. However, in many cases the explicit derivation of a constitutive equation implies the introduction of some closure approximations. The impact of these approximations on the rheological response can be significant [24], even if some of them retain the main features of the kinetic theory model. To avoid the introduction of such approximations, or for checking the validity of the different closure relations, one can resort, within the kinetic theory framework, to a formalism based on the Fokker-Planck equation. The Fokker-Planck equation describes the spatial and temporal evolution of the configurational distribution function associated with the fluid microstructure. The main difficulty related to the solution of that equation concerns the high-dimensional space in

which the distribution function is defined. This space typically includes the physical and conformation coordinates, and time.

Based on the experience gathered in previous works [1,8] it seems that when diffusion effects can span a large numerical range, continuous approximations using a fixed or moving mesh are suitable numerical approaches. In this case accurate stabilizations are required for dealing with small diffusion effects. However, as just pointed out, since the Fokker-Planck equation is defined in a multidimensional space, additional difficulties appear because of the high dimensionality. Some attempts for solving the Fokker-Planck equation using a fixed mesh discretization have been developed [6,33]. The main difficulties in this approach are related to the high-dimensional character of the problem. As a consequence of this high dimensionality, the linear systems obtained after conventional implicit or semi-implicit space-time discretizations result to be extremely large for practical inversion. On the other hand, explicit discretizations, which do not require matrix inversions, tend to suffer from the stability constraint of too small time steps. Another deterministic particle approach, very close to the smooth particles hydrodynamics proposed in [4] was analyzed in [1] using smooth particles. In [1] it was noticed that the numerical effect of smoothing on the solution can be significant.

Another family of simulation techniques exploits the equivalence between the Fokker-Planck equation and Ito's stochastic
<figure><img src="image_2.png" /><figcaption>Fig. 1. Distribution function at times $t = 10^{-3}$, $t = 0.5$, $t = 1$ and $t = 2$.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Four most representative approximation functions.</figcaption></figure>
equation. In these techniques, instead of directly solving the Fokker-Planck equation, a large number of trajectories (individual instances of a stochastic process associated with the Fokker-Planck equation) are integrated. Relevant moments of the distribution function are then calculated by averaging over the ensemble of trajectories. These moments are needed to compute the stress tensor. The first works on this approach are related to the CONNFFESSIT method [40,30,25] (see also the references therein, in particular those concerning the Brownian Configuration Fields extension). This approach was considered in [46] for treating MBS (Multi Bead Spring) models. A similar technique was proposed in [47] in the context of MBS kinetic theory models, which introduces a change of variable and uses a Monte-Carlo technique for accounting for the diffusion term. In these techniques a high number of particles are introduced in the simulation and a stochastic technique is used to account for the Brownian effects. A multi-scale approach using deterministic particles for treating the advection and a different set of particles for accounting for diffusion effects was considered in [23]. A similar idea was applied to the flow calculations of short fiber suspensions in [8]. In that work, the discretization of the advection-dominated Fokker-Planck equation governing fiber orientation was carried out using a particle technique, where the diffusion term was modeled using random walks. An important feature of this approach was that the ensemble size (number of fibers) required for computing the configurational distribution function and its moments, involved in the expression of the stress, increases significantly with the diffusion coefficient. Consequently, for practical applications the use of the method of particles in the framework of a stochastic simulation is restricted to small diffusion effects.

The purpose of the present work is to propose an efficient and accurate discretization technique able to solve non-linear kinetic theory models involving complex behaviors, such as those coming from stability analysis, complex geometries and coupled models. An appealing candidate for carrying the just referred analysis concerns the kinetic theory description of liquid-crystalline polymers. We have chosen Doi's model because it contains all the numerical difficulties that we want to address, and because there is an extensive bibliography on the solution of this model. There is no other reason concerning the election of that model than the pure numerical motivation of using a well-known and established model, even if other physically more realistic LC models exist [3,36]. Hence, our aim is not to provide new rheological results, but to show that the new methodology is of very general applicability, faithfully reproducing the stress tensor.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Evolution of the coefficients related to the four most representative approximation functions.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Evolution of the number of representative functions included in the reduced approximation basis.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Evolution of the order parameter.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Evolution of the first component of the second-order orientation tensor.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

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

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Dependence of the instability related to the apparition of a nematic phase with the value of the dimensionless interaction potential: (top) $U=7$; (middle) $U=8$; (bottom) $U=12$.</figcaption></figure>

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

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

<figure><img src="image_13.png" /><figcaption>**Fig. 8.** Orientation evolution in a shear flow for different diffusion coefficients: (up) $D_r = 1/10$; (middle) $D_r = 1/20$; (down) $D_r = 1/40$.</figcaption></figure>

duces results already reported in the literature, and leads to very significant CPU time savings.

The technique considered here allows for a substantial reduction in the number of degrees of freedom involved in the model discretization, leading to major reductions of the computing time,

sometimes of the order of tens of thousands [2]. It operates by extracting automatically, and in a way completely transparent for the user, the most relevant information of the unknown solution for constructing the functional approximation from the information just extracted. The new shape functions (the most characteristic
<figure><img src="image_14.png" /><figcaption>Fig. 9. Distribution function at three different times: log rolling phenomena.</figcaption></figure>

functions related to the model solution) have a large (non-local) support, i.e. they are defined in the whole domain, as will be described later. Thus, the number of degrees of freedom involved in the solution of the Fokker-Planck equation is reduced by several orders of magnitude. The construction of those new approxima-

tion functions is done with an 'a priori' approach, which combines a basis reduction technique (using the Karhunen-Loève decomposition) with a basis enrichment based on the use of some Krylov subspaces [43,44].

The efficiency of reduction methods for linear models involving homogeneous solutions has been demonstrated in some previous works [2, 7, 41, 44]. However, several questions remained open whose clarification is crucial to guarantee general applicability of the reduction method. These questions concern the ability of the model reduction technique: (i) to treat non-linear models, where the capacity of the algorithm to build the optimal reduced approximation basis has never been checked; (ii) to describe richer behaviors arising from stability analyses; (iii) to treat models defined in complex geometries involving non-homogenous solutions and (iv) to deal with flows in which kinematics and microstructure evolution are intimately coupled. These four aspects constitute some original attempts never until now analyzed (nor published) in the context of model reduction techniques that could contribute to extend the domain of applicability of model reduction techniques in computational rheology for treating efficiently and accurately very complex models defined in 3D complex geometries.

## 2. Liquid-crystalline polymers: the Doi's model

The numerical calculation of flows of complex fluids consisting of highly anisotropic particles (be it at the macroscopic level, e.g. fibers, at the mesoscopic level, e.g. coarse-grained models for polymeric liquid crystals (LC), or at the microscopic level, e.g. low molecular mass liquid-crystalline materials) are not entirely devoid of difficulties. Most of these difficulties stem from the necessity to capture the rich dynamic fluid behavior, which spans both isotropic and anisotropic fluid microstructures (nematic, chiral nematic, smectic, etc.). For example, models describing single particle orientation in a dilute suspension, i.e. in the absence of interparticle interaction, are not able to describe the behavior of a nematic suspension.

Successful treatment of collective flow-orientation phenomena appearing at high concentrations require that interparticle interaction be considered, either explicitly, like in microscopic LC models or in a mean-field sense. Notable examples of mean field approaches to anisotropic particle interaction described in the kinetic theory framework are the model initially proposed by Hess in 1976 as mentioned by Doi and Edwards in [12], or the physically more realistic LC models proposed by Marrucci and coworkers [36–39] and of Beris and Edwards [3, 13–15]. In spite of the complexities associated with the constitutive equations resulting from solving the kinetic theory, major advances have been made in the coupled solution of the hydrodynamic field equations and the LC constitutive equation as well as in the evolution of LC structure in given velocity fields [9–11, 17–22, 26, 35, 42, 45, 48, 49].

Different closure approximations have been proposed for alleviating its numerical simulation operating in the macroscopic scale [5, 16, 28], but in this work we focus on the direct solution of the Fokker-Planck equation associated with Doi's model, that we introduce in the following.

In order to describe the fluid microstructure we firstly define the unit vector **u** describing the orientation of the liquid crystal (LC) molecule axis. Now, the distribution function which gives the probability of finding, at a certain point **X** and time t, molecules aligned in the direction **u** is represented by $\Psi(\mathbf{X},u,t)$. The dynamics of the evolution of this distribution function is induced by (i) hydrodynamic effects, (ii) Brownian effects and (iii) LC-molecule interaction described by a nematic potential. The corresponding
<figure><img src="image_15.png" /><figcaption>Fig. 10. Distribution function at three different times: wagging phenomena.</figcaption></figure>

Fokker-Planck equation related to Doi's model writes:

$$ \frac{\partial \Psi}{\partial t} = -\frac{\partial}{\partial \mathbf{u}} [(\mathrm{grad} \, \mathbf{v}\mathbf{u}) - ((\mathbf{u} \otimes \mathbf{u})(\mathrm{grad} \, \mathbf{v})\mathbf{u}\Psi)] + \frac{\partial}{\partial \mathbf{u}} D_r \left[ \frac{\partial \Psi}{\partial \mathbf{u}} + \Psi \frac{\partial}{\partial \mathbf{u}} \left( \frac{V(\mathbf{u})}{k_b T} \right) \right] \quad (1) $$

where $\partial\Psi/\partial t$ represents the material derivative and therefore it incorporates all the spatial gradients encountered in non-homogeneous suspensions, $k_b$ is Boltzmann's constant, $D_r$ the diffusion coefficient describing the Brownian effects, $T$ the absolute temperature, $\mathbf{v}$ the fluid velocity field and $\partial/\partial\mathbf{u}$ the gradient operator defined on the unit sphere where the distribution function $\Psi(\mathbf{X}, \mathbf{u}, t)$ is defined. A frequently used form for the nematic

potential is:

$$ V(\mathbf{u}) = -\frac{3}{2} U k_b T (\mathbf{u} \otimes \mathbf{u}) : \mathbf{S} \quad (2) $$

where $U$ is a dimensionless interaction potential and the second-order traceless orientation tensor $\mathbf{S}$ is defined by:

$$ \mathbf{S} = \langle \mathbf{u} \otimes \mathbf{u} \rangle - \frac{\mathbf{I}}{3} \quad (3) $$

where the averaging $\langle \cdot \rangle$ is defined on the unit sphere by:

$$ \langle \mathbf{u} \otimes \mathbf{u} \rangle = \int (\mathbf{u} \otimes \mathbf{u}) \Psi(t, \mathbf{u}) d\mathbf{u} \quad (4) $$
<figure><img src="image_16.png" /><figcaption>Fig. 11. Steady-state solution for a Poiseuille flow by coupling the flow kinematics and the microstructure state: (a) U=0; (b) U=2 and (c) U=4.</figcaption></figure>

scalar order parameter can be associated with the second-order tensor S giving the intensity of the orientation, according to:

$$S = \sqrt{\frac{3}{2} S : S} \qquad (5)$$

### 3. Discretization

A natural discretization of this problem defined on the unit surface is based on the use of the spherical polar and azimuthal angles

($\varphi, \theta$):

$$\mathbf{u} = \begin{pmatrix} \cos(\varphi) \cdot \sin(\theta) \\ \sin(\varphi) \cdot \sin(\theta) \\ \mathrm{con}(\theta) \end{pmatrix} \qquad (6)$$

where both angular coordinates are defined in: $\varphi, \theta \in [0,2\pi \times [0,\pi]$. However singularities are encountered at the ends of the $\theta$ definition interval, because of the term $((D_r/\sin^2\theta)(\partial^2\Psi/\partial\varphi^2))$ in the polar expression of the Fokker-Planck equation. Thus, we prefer to work in Cartesian coordinates, where the surface of the unit sphere is approximated by a set of planar triangular facets. The nodal coordinates of the vertices of these triangles are given by:

$$\mathbf{u}_i = \begin{pmatrix} x_i \\ y_i \\ z_i \end{pmatrix} \qquad (7)$$

which satisfy:

$$\sqrt{x_i^2 + y_i^2 + z_i^2} = 1 \qquad (8)$$

and where the natural periodicity is implicitly verified.

From now on, by the sake of simplicity, we consider homogeneous strain rates (complex flows will be addressed later), and then no discretization in physical space is required. Thus, the Fokker-Planck equation can be rewritten in the following form that only involves the conformation (orientation) coordinates and time:

$$\frac{\partial \Psi}{\partial t} + E_0(\mathbf{u})\Psi + \mathbf{E}_1^T(\mathbf{u})\frac{\partial \Psi}{\partial \mathbf{u}} - D_r \frac{\partial^2 \Psi}{\partial \mathbf{u}^2} + H_0(\mathbf{u}, \mathbf{S})\Psi + \mathbf{H}_1^T(\mathbf{u}, \mathbf{S})\frac{\partial \Psi}{\partial \mathbf{u}} = 0 \qquad (9)$$

where the two scalar functions $E_0(\mathbf{u})$, $H_0(\mathbf{u}, \mathbf{S})$ as well as the two vector fields $\mathbf{E}_1(\mathbf{u})$, $\mathbf{H}_1(\mathbf{u}, \mathbf{S})$ can be easily deduced from Eq. (1). We denote by $\Omega$ the domain in the conformation space where Eq. (9) is defined, namely the two-dimensional manifold spanned by the unit orientation vector $\mathbf{u}$, i.e. the surface of the unit sphere. First, the problem is formulated in the Finite Element framework using a weighting function $\Psi^*$:

$$\int_{\Omega} \Psi^* \frac{\partial \Psi}{\partial t} d\mathbf{u} + \int_{\Omega} \Psi^* E_0(\mathbf{u}) \Psi d\mathbf{u} + \int_{\Omega} \Psi^* \mathbf{E}_1^T(\mathbf{u}) \frac{\partial \Psi}{\partial \mathbf{u}} d\mathbf{u} \\ - \int_{\Omega} \Psi^* D_r \frac{\partial^2 \Psi}{\partial \mathbf{u}^2} d\mathbf{u} + \int_{\Omega} \Psi^* H_0(\mathbf{u}, \mathbf{S}) \Psi d\mathbf{u} + \int_{\Omega} \Psi^* \mathbf{H}_1^T(\mathbf{u}, \mathbf{S}) \frac{\partial \Psi}{\partial \mathbf{u}} d\mathbf{u} = 0 \qquad (10)$$

The computational domain $\Omega$ is partitioned into a collection of non-overlapping finite elements (triangular facets covering the surface of the unit sphere). A linear and continuous interpolation of the distribution function is then built in each triangle, form which we can write:

$$\Psi^e(\mathbf{u}) = \sum_{i=1}^{3} N_i(\mathbf{u}) \Psi_i^e \qquad (11)$$

and

$$\Psi^{e*}(\mathbf{u}) = \sum_{i=1}^{3} N_i(\mathbf{u}) \Psi_i^{e*} \qquad (12)$$

where $\Psi_i^e$ and $\Psi_i^{e*}$ are the values at node i of $\Psi$ and $\Psi^*$ respectively, and $N_i(\mathbf{u})$ is the associated shape function which takes a unit value at the node i vanishing at the other nodal positions.
<figure><img src="image_17.png" /><figcaption>Fig. 12. Start-up flow in the driven cavity flow problem. Flow streamlines for: (a) $U=0, t'=0.25$; (b) $U=0, t'=0.75$; (c) $U=0, t'=1.5$; (d) $U=0, t'=5.5$; (e) $U=4, t'=0.25$; (f) $U=4, t'=0.75$; (g) $U=4, t'=1.5$; (h) $U=4, t'=5.5$.</figcaption></figure>

Integrating Eq. (10) by parts and taking into account that the configuration space for **u** is unbounded, it results:

$$ \int_{\Omega} \Psi^* \frac{\partial \Psi}{\partial t} du + \int_{\Omega} \Psi^* E_0(\mathbf{u}) \Psi du + \int_{\Omega} \Psi^* E_1^T(\mathbf{u}) \frac{\partial \Psi}{\partial \mathbf{u}} du \\ + \int_{\Omega} \left( \frac{\partial \Psi^*}{\partial \mathbf{u}} \right)^T D_r \frac{\partial \Psi}{\partial \mathbf{u}} du + \int_{\Omega} \Psi^* H_0(\mathbf{u}, \mathbf{S}) \Psi du $$

$$ + \int_{\Omega} \Psi^* \mathbf{H}_1^T(\mathbf{u}, \mathbf{S}) \frac{\partial \Psi}{\partial \mathbf{u}} du = 0 \qquad (13) $$

Due to the advection-diffusion character of that equation, an appropriate stabilization of the Finite Element scheme is needed to avoid numerical instabilities induced by the convection term. Stabilization is achieved by a usual upwinding formulation, which modifies the weighting functions related to the advection term $\tilde{\Psi}^*$
<figure><img src="image_18.png" /><figcaption>Fig. 13. Start-up flow in the driven cavity flow problem. Orientation of the director for: (a) $U=0, t^*=0.25$; (b) $U=4, t^*=0.25$; (c) $U=0, t^*=5.5$; (d) $U=4, t^*=5.5$.</figcaption></figure>

as described later. The stabilized variational formulation results in:

$$ 
\begin{aligned}
& \int_{\Omega} \Psi^* \frac{\partial \Psi}{\partial t} du + \int_{\Omega} \Psi^* E_0(\mathbf{u}) \Psi du + \int_{\Omega} \tilde{\Psi}^* \mathbf{E}_1^T(\mathbf{u}) \frac{\partial \Psi}{\partial \mathbf{u}} du \\
& + \int_{\Omega} \left( \frac{\partial \Psi^*}{\partial \mathbf{u}} \right)^T D_r \frac{\partial \Psi}{\partial \mathbf{u}} du + \int_{\Omega} \Psi^* H_0(\mathbf{u}, \mathbf{S}) \Psi du \\
& + \int_{\Omega} \tilde{\Psi}^* \mathbf{H}_1^T(\mathbf{u}, \mathbf{S}) \frac{\partial \Psi}{\partial \mathbf{u}} du = 0
\end{aligned}
\quad (14) $$

where in each element

$$ \tilde{\Psi}^{e*}(\mathbf{u}) = \sum_{i=1}^{3} \tilde{N}_i(\mathbf{u}) \Psi_i^{e*} \quad (15) $$

with

$$ \tilde{N}_i(\mathbf{u}) = N_i(\mathbf{u}) + \frac{\beta h}{2} \frac{\partial N_i(\mathbf{u})}{\partial \mathbf{u}} \quad (16) $$

where $h$ is the characteristic element size and $\beta$ is the upwinding parameter given by

$$ \beta = \coth(Pe) - \frac{1}{Pe} \quad (17) $$

where the Peclet number $Pe$ is given by

$$ Pe = \frac{\max(||\mathbf{E}_1(\mathbf{u})||, ||\mathbf{H}_1(\mathbf{u},\mathbf{S})||)h}{2D_r} \quad (18) $$

In order to construct $E_0(\mathbf{u})$ and $H_0(\mathbf{u},\mathbf{S})$ the differential operators must be projected on the plane tangent to the unit sphere. The required gradient projection is as follows: if we consider a triangular facet defined by its three vertices: $(x_1, y_1, z_1)$, $(x_2, y_2, z_2)$, $(x_3, y_3, z_3)$

reference coordinates $(\xi, \eta)$ can be defined using an isoparametric geometrical interpolation:

$$ 
\begin{aligned}
x &= \sum_{i=1}^{3} N_i(\xi, \eta) x_i \\
y &= \sum_{i=1}^{3} N_i(\xi, \eta) y_i \\
z &= \sum_{i=1}^{3} N_i(\xi, \eta) z_i
\end{aligned}
\quad (19) $$

We denote by $\mathbf{e}_1$ and $\mathbf{e}_2$ two vectors defined in the plane containing the triangle, which can be expressed as:

$$ 
\begin{aligned}
\mathbf{e}_1 &= \begin{pmatrix} \frac{\partial x}{\partial \xi} \\ \frac{\partial y}{\partial \xi} \\ \frac{\partial z}{\partial \xi} \end{pmatrix} = \begin{pmatrix} \frac{\partial}{\partial \sum_{i=1}^3 N_i x_i} \\ \frac{\partial}{\partial \sum_{i=1}^3 N_i y_i} \\ \frac{\partial}{\partial \sum_{i=1}^3 N_i z_i} \end{pmatrix} \\
\mathbf{e}_2 &= \begin{pmatrix} \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \eta} \\ \frac{\partial z}{\partial \eta} \end{pmatrix} = \begin{pmatrix} \frac{\partial}{\partial \sum_{i=1}^3 N_i x_i} \\ \frac{\partial}{\partial \sum_{i=1}^3 N_i y_i} \\ \frac{\partial}{\partial \sum_{i=1}^3 N_i z_i} \end{pmatrix}
\end{aligned}
\quad (20) $$

Now, the relation between the gradient in the Cartesian coordinates and the one expressed on the manifold tangent to the unit
<figure><img src="image_19.png" /><figcaption>Fig. 14. Streamlines related to the flow of Newtonian fluid at very low Reynolds number in the 4:1 contraction.</figcaption></figure>

sphere is given by [31]:

$$ \begin{pmatrix} \frac{\partial \Psi^e}{\partial x} \\ \frac{\partial \Psi^e}{\partial y} \\ \frac{\partial \Psi^e}{\partial z} \end{pmatrix} = (e_1 \quad e_2) \begin{pmatrix} e_1^T e_1 & e_1^T e_2 \\ e_2^T e_1 & e_2^T e_2 \end{pmatrix}^{-1} \begin{pmatrix} \frac{\partial \Psi^e}{\partial \xi} \\ \frac{\partial \Psi^e}{\partial \eta} \end{pmatrix} \quad (21) $$

where

$$ \Psi^e = \sum_{i=1}^{3} N_i(\xi, \eta) \Psi_i^e \quad (22) $$

and the matrix

$$ \begin{pmatrix} e_1^T e_1 & e_1^T e_2 \\ e_2^T e_1 & e_2^T e_2 \end{pmatrix} $$

contains the elements of the metric tensor g corresponding to a given triangular facet. Note that the required properties of differentiability, symmetry, non-singularity and positive-definiteness are satisfied by construction in Cartesian coordinates.

As it is usual in the FEM framework, after numerical integration and assembling of element contributions the following non-linear system is obtained:

$$ M \dot{\Psi} + (G + H(\Psi)) \Psi = 0 \quad (23) $$

where $\Psi$ refers to the vector containing the nodal values of the distribution function.

The non-linearity follows from the dependence of the nematic potential on the average of the distribution function itself on the unit sphere. When one uses an explicit time discretization technique, the non-linear system can be linearized by considering the nematic potential calculated at the previous time step.

## 4. Model reduction

In this section we summarize the reduction technique widely described in [2] that will be applied to the discrete problem resulting from the discretization of the Fokker-Planck equation that governs the evolution of the configuration distribution. We describe in this section the main ideas of the adaptation strategy that we are using. However, in order to be efficient from a computational point of view, the numerical algorithm that we are using is somewhat more complex. Interested readers are addressed to [2,41] and [44] for a fully detailed description.

We start by assuming that $\Psi$ can be expressed from the initial condition according to $\Psi(t) = \Psi(t=0) \times \xi(t) = B\xi(t)$, where the basis transformation matrix B has only a column and vector $\xi$ reduces to a scalar. Now, we compute S iterations of the evolution problem using the reduced model derived from Eq. (23):

$$ B^T M B \dot{\xi} + B^T (G + H(B\xi)) B \xi = 0 \quad (24) $$

without changing the approximation basis. In Eq. (24) the size of the matrix that must be inverted to update vector $\xi$ coincides with the

number of approximation functions defining matrix B. After these S iterations the residual R is evaluated:

$$ R = M B \dot{\xi} + (G + H(B\xi)) B \xi \quad (25) $$

If the norm of the residual is small enough, $||R|| < \epsilon$, with $\epsilon$ a threshold value small enough, we can continue for other S iterations using the same reduced approximation basis B. On the contrary, if the residual norm is too large, $||R|| \ge \epsilon$, we need to enrich the approximation basis and compute again the last S iterations. This enrichment is built using some Krylov's subspaces related to the residual. The enrichment when only the first Krylov's subspace is considered reads: $B \leftarrow (B, R)$.

One could expect that the enrichment process is increasing continuously the size of the reduced approximation basis, but in fact after reaching the convergence a Karhunen-Loève decomposition is performed on the whole past time interval in order to extract the significant information and to define an orthogonal reduced approximation basis allowing to stabilize and keep moderate the size of the discrete reduced model (see [2] for more details).

## 5. Results

### 5.1. LCP in an elongational flow

We first present calculation results for an elongational flow whose kinematics is defined by grad $v = \text{diag}(1, -0.5, -0.5)$, the initial orientational distribution being isotropic: $\psi(u, t=0) = 1/4\pi$. The dimensionless interaction potential is set to $U=7.5$ and the isotropic diffusion coefficient to $D_r = 0.25$. A mesh of 2560 nodes on the unit sphere was used. This size of the mesh is also the size of the linear system that must be solved at each iteration when using a semi-implicit finite element technique, in contrast with the only six characteristic functions that are needed to represent the whole evolution of the distribution function when using the present method as described later.

Some snapshots taken during the simulation calculation are represented in Fig. 1 both on the unit sphere and using a spherical coordinates representation. The solution at time $t=2$ is already very close to the steady-state solution. Here, only 6 functions are required to represent accurately the evolution of the solution, instead the 2560 degrees of freedom necessary in a standard finite element formulation. These six functions are computed and introduced into the approximation basis when their contribution is needed using the adaptive procedure described in the previous section. Thus, the adaptive construction of the reduced approximation basis takes place in a fully automatic and transparent manner. Fig. 2 shows the four most significant functions defined by a colored contour map on the unit sphere. The evolution of the weights associated with these four most representative functions is depicted in Fig. 3. Fig. 4 represents the number of approximation functions used in the problem discretization at the instant when an enrichment of the reduced approximation basis takes place. The evolution of the order parameter is shown in Fig. 5 in agreement with the analytical result that predicts a steady order parameter of 0.86. Fig. 6 illustrates the evolution of the first component of the second-order orientation tensor. This component is the most relevant of the orientation tensor because it indicates the intensity of the fiber alignment along the elongation axis.

### 5.2. Appearance of a nematic phase

The second numerical example emphasizes the ability of the reduced technique just proposed to describe the appearance of the nematic phase, which marks the onset of the instability of the isotropic orientational distribution. This phenomenon has been
<figure><img src="image_20.png" /></figure>

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

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

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

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

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

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

<figure><img src="image_27.png" /><figcaption>**Fig. 15.** Solution along the streamline beginning at point (0, -0.1): (a) velocity gradient components; (b-e) fourth most significant approximation functions; (f) approximation functions weights; (g) reconstructed distribution function and (h) evolution of the number of approximation functions.</figcaption></figure>

analyzed, for example, in [34]. In this numerical test we have
prescribed a very small shear rate ($O(10^{-7}) s^{-1}$), a constant dif-
fusion coefficient $D_r = 0.25$, and a series of increasing values of
the dimensionless interaction potential $U$. For small enough val-
ues of $U$ the suspension remains isotropic. As the critical value is
reached a nematic phase appears. Fig. 7 confirms this expected
result, and proves that the characteristic time related to the

appearance of the nematic phase increases as the value of *U*
decreases.

5.3. LCP in a simple shear flow

In order to carry out a quantitative comparison of the results
obtained with the present method (using reduced approximation
<figure><img src="image_28.png" /></figure>

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

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

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

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

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

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

<figure><img src="image_35.png" /><figcaption>**Fig. 16.** Solution along the streamline beginning at point (0, -0.5): (a) velocity gradient components; (b-e) fourth most significant approximation functions; (f) approximation functions weights; (g) reconstructed distribution function and (h) evolution of the number of approximation functions.</figcaption></figure>

bases) with those reported in [29] we prescribe a unit shear rate, a dimensionless interaction potential $U = 10.67$ and three different diffusion coefficients: $D_r = 1/10, 1/20$ and $1/40$. We consider again an isotropic initial orientation distribution. Fig. 8 shows the expected oscillatory response, without reaching a steady state (at least in a reasonable integration time). Furthermore, it can be

observed that the oscillation amplitude decreases as the diffusion coefficient increases. The value of the oscillation period is in complete agreement with the one reported in [29]. Besides classical solutions presented in [29], the Fokker-Planck equation further exhibits irregular or chaotic orbits, and a rich non-equilibrium phase diagram [27].
### 5.4. Stability analysis

We analyze in the present section the ability of the reduced basis approximation to capture rich microstructure evolutions, as encountered in 3D shear flows, and reported in [29]. We consider again an interaction potential $U=10.67$, a diffusion coefficient $D_r=1$ and a flow kinematic given by: $(u, v, w) = (0, 0, \dot{\gamma}x)$. Now, we consider two points in the stability diagram reported in Fig. 9 of [29], the first one given by an orientation director angle of $-45^\circ$ and a shear rate $\dot{\gamma} = 8$ and the second one by the same orientation of the director but a higher shear rate $\dot{\gamma} = 25$. From the numerical solution we confirm that in the first case the distribution evolves towards a log rolling state (see Fig. 9) whereas in the second case the solution exhibits the expected wagging phenomena (see Fig. 10).

### 5.5. Coupling microstructure and flow kinematics

We address in the present section a first coupling between the flow kinematics and the microstructure evolution. For this purpose we consider two flow problems: first the steady solution of a Poiseuille flow involving a liquid crystal; second, the start-up of the flow confined in a driven cavity. In both cases we must solve the Fokker-Planck Eq. (1) and the momentum equation that neglecting inertia terms can be written as:

$$-\nabla p + 2\eta_s \nabla \mathbf{D} = \nabla \boldsymbol{\tau}_p \qquad (26)$$

where $p$ is the pressure field, $\eta_s$ the solvent viscosity, $\mathbf{D}$ the strain rate tensor (symmetric component of the gradient of velocity tensor) and $\boldsymbol{\tau}_p$ the polymer contribution to the stress. Eq. (26) must be solved under the incompressibility restriction:

$$\nabla \cdot \mathbf{v} = 0 \qquad (27)$$

and after specifying the expression of the polymer contribution to the stress $\boldsymbol{\tau}_p$. In the next simulations we are assuming $\boldsymbol{\tau}_p$ given by:

$$\boldsymbol{\tau}^p = \boldsymbol{\tau}_B^p + \boldsymbol{\tau}_{VE}^p \left\{ \begin{array}{l} -\frac{\nu k_B T}{2} \int_{S(0,1)} \left( (\mathbf{u} \otimes \mathbf{u}) - \frac{\mathbf{I}}{3} \right) \nabla_u \cdot (\nabla_u) \, d\mathbf{u} \\ -\frac{\nu k_B T}{2} \int_{S(0,1)} \left( (\mathbf{u} \otimes \mathbf{u}) - \frac{\mathbf{I}}{3} \right) \nabla_u \cdot \left( \nabla_u \left( \frac{V}{k_B T} \right) \right) \, d\mathbf{u} \end{array} \right. \qquad (28)$$

where $\nabla_u$ denotes the gradient differential operator in the conformation space $\mathbf{u}$, that was refereed by $\partial/\partial\mathbf{u}$ in Section 2, $\mathbf{I}$ is the unit tensor, $v$ represents the number of molecules by unit of volume, $K_B$ is the Boltzmann's constant and $T$ is the temperature.

Before solving the system (26)–(28) we proceed to write these equations in a dimensionless form, taking $1/D_r$ as characteristic time, the channel width $L$ as characteristic length and $\nu k_B T$ as characteristic pressure, and by applying the following relations (where the dimensionless variables are marked with asterisks):

$$\left\{ \begin{array}{l} \mathbf{t}^* = t \cdot D_r \\ \mathbf{x}^* = \frac{\mathbf{x}}{L} \\ \text{grad } \mathbf{v}^* = \frac{\text{grad } \mathbf{v}}{D_r} \\ \boldsymbol{\tau}^* = \frac{\boldsymbol{\tau}}{\nu k_B T} \\ p^* = \frac{p}{\nu k_B T} \\ \eta_s^* = \frac{\eta_s \cdot D_r}{\nu k_B T} \end{array} \right. \qquad (29)$$

We have solved the dimensionless model equations for 3 different dimensionless interaction potentials ($U=0, U=2$ and $U=4$) in the nematic potential defined by Eq. (2). In the first test, the one concerning the steady solution of the Poiseuille flow, we fixed

the gradient of pressure $dp^*/dx^* = -1$ and the dimensionless viscosity $\eta^* = 0.1$. The flow equations (dimensionless form of Eqs. (26)–(28)) were discretized by using a mixed stable finite element discretization. On the other hand, the steady solution of the Fokker-Planck equation was computed on each streamline by integrating it on the conformational space (the orientation distribution does not evolve in the flow direction). The coupling between both solvers was accomplished within a fixed point strategy. In these simulations we look directly for the steady Poiseuille solution starting from a parabolic velocity profile and an isotropic microstructure state. Fig. 11 depicts the initial velocity profile (black curve) and the one obtained after convergence (red curve) that was reached in few iterations. We can notice that for low values of the interaction potential, Fig. 11(a), the diffusion effects are dominant, leading to a microstructure distribution whose director (also illustrated in this figure at different points) remains everywhere misaligned with respect to the flow direction. However, the misalignment magnitude decreases as the domain walls are approached because at those regions the shear rate that tends to align the molecules in the flow direction increases significantly. Due to this misalignment the final steady velocity profile shows a plateau in the neighborhood of the flow symmetry axis. When the interaction potential increases, the alignment is facilitated and consequently the molecules tend to align in the flow direction, as can be observed by comparing the directors depicted in Fig. 11(a) and (c). The increase in interaction potential also causes the velocity profile to approach the parabolic one.

Now we considered a more complex flow: the start-up flow of a LCP confined in a square driven cavity whose lateral walls have a length $L$. The upper wall moves from the left to the right, with a horizontal velocity profile defined by a fourth-order polynomial:

$$\mathbf{v}(x, y, L)^T = (16\nu_{\text{max}}x^2(L-x)^2, 0) \qquad (30)$$

The maximum velocity on the upper wall, $\nu_{\text{max}}$ is fixed in order to have a unit Weissenberg number, that is: $W_e = \nu_{\text{max}}/L = 1$. The flow equations (dimensionless form of Eqs. (26)–(28)) were discretized by using a mixed stable finite element discretization. On the other hand, the transient Fokker-Planck equation was integrated by using the method of characteristics along several flow trajectories, from an isotropic orientation distribution assumed at the initial time, using the adaptive reduced approximation basis and the algorithm previously described. We must notice that because the Fokker-Planck equation does not involve diffusion terms in the physical coordinates, its integration does not require the prescription of boundary condition for the orientation distribution on the cavity walls. The coupling between both solvers was accomplished within a fixed point strategy. Fig. 12 compares the streamlines at four different times $t^* = 0.25, t^* = 0.75, t^* = 1.5$ and $t^* = 5.5$ for the two interaction potentials ($U=0$ and $U=4$). There is a significant impact of the interaction potential on the flow kinematics, in particular on the position of the stagnation point of the recirculation. Fig. 13 depicts the director field at $t^* = 0.25$ and $t^* = 5.5$ for both interaction potentials.

### 5.6. The 4:1 contraction

As the in the driven cavity flow problem, for computing the steady-state solution in a 4:1 contraction we couple the flow kinematics and the Fokker-Planck equation solvers by using a fixed point strategy. Thus, the flow kinematics is searched from the microstructure computed at the previous iteration (or at the previous time step in transient simulations), and from this solution the fluid microstructure can be updated, and so on. As the solution of the flow kinematics does not present particular dif-
<figure><img src="image_36.png" /><figcaption>**Fig. 17.** Solution along the streamline beginning at point (0, -0.9): (a) velocity gradient components; (b-e) fourth most significant approximation functions; (f) approximation functions weights; (g) reconstructed distribution function and (h) evolution of the number of approximation functions.</figcaption></figure>

ficulties, we will focus on the solution of the microstructure evolution.

Obviously in the scenario considered in this section that involves non-homogeneous conformational distributions, the spatial gradi-

ent involved in the material derivative of the Fokker-Planck Eq. (1) must be incorporated by considering an appropriate discretization technique accounting for its advective character. The advection character in the physical coordinates implies that the solution at a
<figure><img src="image_37.png" /></figure>

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

<figure><img src="image_39.png" /><figcaption>**Fig. 18.** Evolution of the order parameter along the three streamlines starting at points: (a) (0, −0.1); (b) (0, −0.5) and (c) (0, −0.9).</figcaption></figure>

certain point of the physical domain only depends on the upstream solution along the associated trajectory. Thus, the Fokker-Planck equation can be easily integrated along an arbitrary number of streamlines. Now, from the solution of the distribution function along those streamlines, the solution can be defined everywhere

<figure><img src="image_40.png" /><figcaption>**Fig. 19.** Steady-state flow streamlines computed after convergence of the coupled model.</figcaption></figure>

using an appropriate interpolation. The integration by the method of characteristics allows the direct use of the reduction strategy previously described.

First, we assume a flow kinematics coming from the solution of the Stokes equations (the inertia term is always neglected in our simulations). Fully developed flow conditions are imposed on the inlet and outlet boundaries. Fig. 14 depicts some streamlines, where only a half of the computational domain is considered due to the flow symmetry. Now, assuming a 3D isotropic orientation distribution on the inflow boundary, Figs. 15–17 depict the evolution of the components of the gradient of velocity tensor along the streamlines starting at points (0, −0.1), (0, −0.5) and (0, −0.9) respectively, the fourth most significant approximation functions, the evolution of the weight of each one of these functions along the streamlines, the final distribution (just at the intersection between the considered streamline and the outflow boundary) and the evolution of the number of functions involved in the reduced approximation. In all cases $D_r = 0.23$ and $U = 7.21$. Fig. 18 compares the evolution of the order parameter along these three streamlines.

From the computed solution of the Fokker-Planck equation along numerous streamlines, the extra-stress tensor can be computed on these streamlines, and from that, extended to any point in the domain by using an appropriate interpolation. Now, the flow kinematics can be recomputed, and then the microstructure updated. The iteration process continues until reaching convergence.

In is well known that short fibers suspensions, even at very low fiber concentrations, exhibit significant non-Newtonian behaviors that is the case of the 4:1 contraction here considered results in the appearance of a noticeable recirculating flow located in the contraction corner and whose size increases with the fiber concentration [32]. Thus, one could expect a similar behavior in the case of a LCP flow in the 4:1 contraction. To elucidate this question we computed different streamline of the steady-state flow obtained after the convergence of the fixed point iteration algorithm. These streamlines are depicted in Fig. 19 from which we can notice the existence of the expected recirculating flow, invalidating any attempt of computing the microstructure by using the Stokes velocity field (fully decoupled technique).

## 6. Conclusion

The model reduction strategy described in this paper has been successfully applied for treating non-linear kinetic theory models involving complex behaviors, as the ones arising from stability analysis, complex geometries and coupled models. While significantly reducing in the number of approximation functions involved in the simulation – and consequently the computational effort – the computed solutions in flows involving elongation or shear, remain very accurate. Major CPU time savings, ranging from thousands to tens of thousands, can be achieved. Comparisons with previous published simulations show that the present method is also able to efficiently and accurately predict the isotropic-nematic transition, complex flow simulation as well as complex microstructure evolutions implying richer out-of-plane orientations.
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, in: Lecture Notes on Computational Science and Engineering, vol. 43, Springer, 2004, pp. 1–16.

[2] A. Ammar, D. Ryckelynck, F. Chinesta, R. Keunings, On the reduction of kinetic theory model related to finitely extensible dumbbells, J. Non-Newtonian Fluid Mech. 134 (2006) 136–147.

[3] A.N. Beris, B.J. Edwards, Thermodynamics of Flowing Systems With Internal Microstructure, Oxford University Press, 1994.

[4] 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.

[5] C.V. Chaubal, L.G. Leal, A closure approximation for liquid-crystalline polymer models based on parametric density estimation, J. Rheol. 42 (1998) 177–202.

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

[7] F. Chinesta, A. Ammar, A. Falco, M. Laso, On the reduction of stochastic kinetic theory models of complex fluids, Model. Simulat. Mater. Sci. Eng. 15 (2007) 639–652.

[8] 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 (2004) 97–125.

[9] Z.L. Cui, et al., On weak plane Couette and Poiseuille flows of rigid rod and platelet ensembles, SIAM J. Appl. Math. 66 (2006) 1227–1260.

[10] L.R.P. de Andrade Lima, A.D. Rey, Poiseuille flow of discotic nematic liquid crystals' onion textures, J. Non-Newtonian Fluid Mech. 119 (2004) 71–81.

[11] L.R.P. de Andrade Lima, A.D. Rey, Poiseuille flow of Leslie-Ericksen discotic liquid crystals: solution multiplicity, multistability, and non-Newtonian rheology, J. Non-Newtonian Fluid Mech. 110 (2003) 103–142.

[12] M. Doi, S.F. Edwards, The Theory of Polymer Dynamic, Oxford University Press, London, 1989.

[13] B.J. Edwards, A.N. Beris, M. Grmela, The dynamic behavior of liquid-crystals—a continuum description through generalized brackets, Mol. Cryst. Liquid Cryst. 201 (1991) 51–86.

[14] B.J. Edwards, A.N. Beris, M. Grmela, Generalized constitutive equation for polymeric liquid-crystals. 1. Model formulation using the Hamiltonian (Poisson Bracket) formulation, J. Non-Newtonian Fluid Mech. 35 (1990) 51–72.

[15] B.J. Edwards, A.N. Beris, M. Grmela, R.G. Larson, Generalized constitutive equation for polymeric liquid-crystals. 2. Nonhomogeneous systems, J. Non-Newtonian Fluid Mech. 36 (1990) 243–254.

[16] J. Feng, C.V. Chaubal, L.G. Leal, Closure approximations for the Doi theory: which to use in simulating complex flows of liquid-crystalline polymers? J. Rheol. 42 (1998) 1095–1119.

[17] J. Feng, L.G. Leal, Pressure-driven channel flows of a model liquid-crystalline polymer, Phys. Fluids 11 (1999) 2821–2835.

[18] J. Feng, L.G. Leal, Simulating complex flows of liquid-crystalline polymers using the Doi theory, J. Rheol. 41 (1997) 1317–1335.

[19] M.G. Forest, R.H. Zhou, Q. Wang, Nano-rod suspension flows: a 2D Smoluchowski–Navier–Stokes solver, Int. J. Numer. Anal. Model. 4 (2007) 478–488.

[20] M.G. Forest, R.H. Zhou, Q. Wang, Kinetic structure simulations of nematic polymers in plane Couette cells. II: In-plane structure transitions, Multiscale Model. Simul. 4 (2005) 1280–1304.

[21] D. Grecov, L.R.P. de Andrade Lima, A.D. Rey, Multiscale simulation of flow-induced texture formation in polymer liquid crystals and carbonaceous mesophases, Mol. Simul. 31 (2005) 185–199.

[22] D. Grecov, A.D. Rey, Texture control strategies for flow-aligning liquid crystal polymers, J. Non-Newtonian Fluid Mech. 139 (2006) 197–208.

[23] 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.

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

[25] 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.

[26] D.H. Klein, L.G. Leal, C.J. García-Cervera, H.D. Ceniceros, Ericksen number and Deborah number cascade predictions of a model for liquid crystalline polymers for simple shear flow, Phys. Fluids 19 (2007) 023101.

[27] M. Kroger, Simple models for complex non-equilibrium fluids, Phys. Rep. 390 (2004) 453–551.

[28] M. Kroger, A. Ammar, F. Chinesta, Consistent closure schemes for statistical models of anisotropic fluids, J. Non-Newtonian Fluid Mech. 149 (2008) 40–55.

[29] R.G. Larson, H.C. Öttinger, Effect of molecular elasticity on out-of-plane orientation in shearing flows of liquid-crystalline polymers, Macromolecules 24 (1991) 6270–6282.

[30] M. Laso, H.C. Öttinger, Calculation of viscoelastic flow using molecular-models—the CONNFESSIT approach, J. Non-Newtonian Fluid Mech. 47 (1993) 1–20.

[31] M. Lipschutz, Differential Geometry, McGraw-Hill, 1970.

[32] G.G. Lipscomb, M.M. Denn, D.U. Hur, D.V. Boger, J. Non-Newtonian Fluid Mech. 26 (1988) 297–325.

[33] 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.

[34] P.L. Maffettone, M. Grosso, C. Friedenberg, G.G. Fuller, Extensional flow of a two-dimensional polymer liquid crystal, Macromolecules 29 (1996) 8473–8478.

[35] D. Marenduzzo, E. Orlandini, J.M. Yeomans, Hydrodynamics and rheology of active liquid crystals: a numerical investigation, Phys. Rev. Lett. 98 (2007) 118102.

[36] G. Marrucci, Rheology of liquid-crystalline polymers, Pure Appl. Chem. 57 (1985) 1545–1552.

[37] G. Marrucci, P.L. Maffettone, Description of the liquid-crystalline phase of rod-like polymers at high shear rates, Macromolecules 22 (1989) 4076–4082.

[38] G. Marrucci, P.L. Maffettone, Nematic phase of rodlike polymers. 1. Prediction of transient-behavior at high shear rates, J. Rheol. 34 (1990) 1217–1230.

[39] G. Marrucci, P.L. Maffettone, Nematic phase of rodlike polymers. 2. Polydomain predictions in the tumbling regime, J. Rheol. 34 (1990) 1231–1244.

[40] H.C. Öttinger, M. Laso, Smart Polymers in Finite Element Calculations, Int Congr. on Rheology, Brussels, Belgium, 1992.

[41] E. Pruliere, A. Ammar, N. El Kissi, F. Chinesta, Recirculating flows involving short fiber suspensions: numerical difficulties and efficient advanced micro-macro solvers, Arch. Comput. Methods Eng. State Art Rev. 16 (2009) 1–30.

[42] A.D. Rey, M.M. Denn, Dynamical phenomena in liquid-crystalline materials, Annu. Rev. Fluid Mech. 34 (2002) 233–266.

[43] D. Ryckelynck, A priori hyper-reduction method: an adaptive approach, J. Comput. Phys. 202 (2005) 346–366.

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

[45] G. Sgalari, L.G. Leal, E. Meiburg, Texture evolution of sheared liquid crystalline polymers: numerical predictions of roll-cells instability, director turbulence, and striped texture with a molecular model, J. Rheol. 47 (2003) 1417–1444.

[46] 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 (2002) 227–255.

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

[48] H.J. Yu, P.W. Zhang, A kinetic-hydrodynamic simulation of microstructure of liquid crystal polymers in plane shear flow, J. Non-Newtonian Fluid Mech. 141 (2007) 116–127.

[49] R.H. Zhou, M.G. Forest, Q. Wang, Kinetic structure simulations of nematic polymers in plane Couette cells. I: The algorithm and benchmarks, Multiscale Model. Simul. 3 (2005) 853–870.