# A Fully Deterministic Micro-Macro Simulation of Complex Flows Involving Reversible Network Fluid Models

Bechir Mokdad, Amine Ammar, Magdeleine Normandin, Francisco Chinesta, Jean-Robert Clermont

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

## To cite this version:

Bechir Mokdad, Amine Ammar, Magdeleine Normandin, Francisco Chinesta, Jean-Robert Clermont. A Fully Deterministic Micro-Macro Simulation of Complex Flows Involving Reversible Network Fluid Models. Mathematics and Computers in Simulation, 2010, 80 (9), pp.1936-1961. 10.1016/j.matcom.2010.03.002. hal-01004833

HAL Id: hal-01004833

https://hal.science/hal-01004833

Submitted on 12 Feb 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A fully deterministic micro–macro simulation of complex flows involving reversible network fluid models

B. Mokdad<sup>a</sup>, A. Ammar<sup>a</sup>, M. Normandin<sup>a</sup>, F. Chinesta<sup>b</sup>, J.R. Clermont<sup>a</sup>

Micro–macro models associate the coarse-grained molecular scale of the kinetic theory to the macroscopic scale of continuum mechanics. The conservation equations are solved along with the microscopic equation or the so-called Fokker–Planck equation. In this paper, a micro–macro approach based on the separated representation introduced in [2,3] with the Stream-Tube method [10–12,21,22] is implemented to study the main features of fiber and polymer networks solutions in complex flows. The Fokker–Planck equation, that defines the fluid microstructure, is solved using a separated representation strategy and is coupled to the macroscopic equations through the macroscopic extra-stress tensor evaluated at the microscopic level. Then, the flow kinematics is solved by applying the Stream-Tube method.

*Keywords:* Kinetic theory; Micro–macro approach; Fiber suspensions; Stream-Tube method; Reversible network models

## 1. Introduction

Many natural and synthetic fluids are viscoelastic materials, in the sense that the stress endured by a macroscopic fluid element depends upon the history of the deformation experienced by that element. Notable examples include polymer solutions and melts, liquid crystalline polymers and fiber suspensions. Rheologists thus face a challenging non-linear coupling between flow-induced evolution of microscopic configurations, macroscopic rheological response, flow parameters (such as the geometry and boundary conditions) and final properties. Theoretical modelling and methods of computational rheology have an important role to play in elucidating this coupling.

Atomistic modelling is the most detailed level of description that can be applied today in rheological studies, using techniques of non-equilibrium molecular dynamics. Such calculations require enormous computer resources, and then they are currently limited to flow geometries of molecular dimensions. Consideration of macroscopic flows found in processing applications calls for less detailed mesoscopic models, such as those of kinetic theory.

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

<footnote><sup>b</sup> EADS Corporate Foundatuation International Chair, GeM UMR CNRS - Centrale Nantes, 1 rue de la Noe, BP 92101, F-44321 Nantes Cedex 3, France</footnote>
Models of kinetic theory provide a coarse-grained description of molecular configurations wherein atomistic processes are ignored. They are meant to display in a more or less accurate fashion the important features that govern the flow-induced evolution of configurations. Over the last few years, different models related to dilute polymers have been evaluated in simple flows by means of stochastic simulation or Brownian dynamics methods.

Kinetic theory models can be very complicated mathematical objects. It is usually not easy to compute their rheological response in rheometric flows, and their use in numerical simulations of complex flows has long been thought impossible. The traditional approach has been to derive from a particular kinetic theory model a macroscopic constitutive equation that relates the viscoelastic stress to the deformation history. One then solves the constitutive model together with the conservation laws using a suitable numerical method, to predict velocity and stress fields in complex flows. The majority of constitutive equations used in continuum numerical simulations are indeed derived (or at least very much inspired) from kinetic theory.

Indeed, derivation of a constitutive equation from a model of kinetic theory usually involves closure approximations of a purely mathematical nature such as decoupling or pre-averaging. It is now widely accepted that closure approximations have a significant impact on rheological predictions for dilute polymer, solutions, or fiber suspensions.

In this context, micro–macro methods of computational rheology that couple the coarse-grained molecular scale of kinetic theory to the macroscopic scale of continuum mechanics have an important role to play. In the one hand, This approach is much more demanding in computer resources than more conventional continuum simulations that integrate a constitutive equation to evaluate the viscoelastic contribution of the stress tensor. On the other hand, micro–macro techniques allow the direct use of kinetic theory models and thus avoid the introduction of closure approximations.

Since the early 1990s the field has developed considerably following the introduction of the CONNFFESSIT method by Öttinger and Laso [25]. Being relatively new, micro–macro techniques have been implemented only for models of kinetic theory with few configurational degrees of freedom, such as non-linear dumbbell models of dilute polymer solutions and single-segment tube models of linear entangled polymers.

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

$$ \frac{d\Psi}{dt} + \frac{\partial}{\partial \underline{X}} (\underline{A}\Psi) = \frac{1}{2} \frac{\partial}{\partial \underline{X}} \frac{\partial}{\partial \underline{X}} : (\underline{B}\Psi) \quad (1) $$

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

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

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

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

$$ \left\{ \begin{array}{l} \nabla \cdot \underline{v} = 0 \\ \nabla \cdot (-p\underline{I} + \underline{\tau}_p + \eta_s\underline{D}) = 0 \end{array} \right. \quad (3) $$

where $p$ is the fluid pressure and $\eta_s\underline{D}$ is a purely viscous component ($\underline{D}$ being the strain rate tensor and $\eta_s$ the solvant viscosity). The set of coupled Eqs. (1)–(3), supplemented with suitable initial and boundary conditions in both physical and configuration spaces, is the generic multiscale formulation. Three basic approaches have been adopted for exploiting the generic multiscale model:
(1) The continuum approach wherein a constitutive equation of continuum mechanics that relates the viscoelastic stress to the deformation history is derived from, and replaces altogether, the kinetic theory model (1) and (2). The derivation process usually involves closure approximations. The resulting constitutive model takes the form of a differential, integral or integro-differential equation. It yields molecular information in terms of averaged quantities, such as the second moment of the distribution: $\int_C \underline{X} \otimes \underline{X} \Psi d\underline{X}$.

(2) The Fokker-Planck approach wherein one solves the generic problem (1)-(3) as such, in both configuration and physical spaces. The distribution function is thus computed explicitly as a solution of the Fokker-Planck Eq. (1). The viscoelastic stress is computed from (2).

(3) The Stochastic approach which draws on the mathematical equivalence between the Fokker-Planck Eq. (1) and the following Ito stochastic differential equation:

$$ \underline{X} = \underline{A} \, dt + \underline{b} \, d\underline{W} \qquad (4) $$

where $\underline{B} = \underline{b} \underline{b}^T$ and $\underline{W}$ is a Wiener stochastic process of dimension $N$. In a complex flow, the stochastic differential Eq. (4) applies along individual flow trajectories, the time derivative is thus a material derivation. Instead of solving the deterministic Fokker-Planck Eq. (1), one solves the associated stochastic differential Eq. (4) for a large ensemble of realizations of the stochastic process $\underline{X}$ by means of a suitable numerical technique. The distribution function is not computed explicitly, and the viscoelastic stress (2) is readily obtained as an ensemble average.

For more details concerning the micro-macro approach the reader can refer to the review paper [17] and the references therein.

The continuum approach has been adopted throughout the development of computational rheology. First simulations were obtained in the early 1980s. Two decades later, macroscopic numerical techniques based upon the continuum approach remains under active development.

The control of the statistical noise is a major issue in stochastic micro-macro simulations based on the stochastic approach. Moreover, to reconstruct the distribution one needs to operate with an extremely large number of particles, however, in general, only the moments of such distribution are required, which can be computed using a much more reduced population of particles. These difficulties do not arise at all in the Fokker-Planck approach, however, the Fokker-Planck Eq. (1) must be solved for the distribution function in both physical and configuration spaces. This necessitates a suitable discretization procedure for all relevant variables, namely position $\underline{x}$, configuration $\underline{X}$ and time $t$. Until now, the dimensionality of the problem could be daunting and consideration of molecular models with many configurational degrees of freedom did not appear feasible. This probably explains why relatively few studies based of the Fokker-Planck approach have appeared in the literature until very recently at least. In [9,20] the resolution of the Fokker-Planck equation involving a moderate number of dimensions is considered. An alternative deterministic particle approach has been proposed in [8]. Sometimes 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 in the domain of interest results simply prohibitory. Some attempts exist concerning the treatment of multidimensional problems. The interested reader can refer to [7] 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 conclusion, the problematic lied to models defined in multidimensional spaces remains still open, and new efforts must be paid to reach significant improvements in next years.

Some of the most used kinetic theory models defined from a Fokker-Planck equation have two important particularities: (i) they can be expressed in a separated form (this is the case of multi-bead models) and (ii) in general, realistic molecular models involve springs with finite extension which implies that the distribution function vanishes on the boundary of the domain where the springs conformation is defined. In this case, the separated representation and the definition of tensor product approximation spaces, run perfectly, and allow to circumvent the difficulties related to the multidimensional character of kinetic theory models, as proved in [2]. This technique consists of the use of a separated representation of the molecular conformation distribution, that introduced in the variational formulation of the Fokker-Planck equation leads to an iteration procedure that involves at each step the resolution of a small-size non-linear problem. The resolution of those non-linear problems can be performed by using a standard Newton strategy, alternating directions resolution or more sophisticated strategies. Thus, the number of degrees of freedom involved in
the discretization of the Fokker-Planck equation can be reduced from $(n_n)^N$ (required in the usual grid/mesh based strategies) to $(n_n) \times N$ $(n_n$ being the number of nodes involved in the discretization of each conformation coordinate and $N$ the dimension of the conformation space).

In [2] 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. In the conclusions of that paper was claimed that the resolution of multidimensional transient Fokker-Planck equations could be considered within an incremental time discretization. However, being the time no more than other coordinate, one could expect a coupled space-time resolution. The main difficulties related to a such approach lies in the fact that the initial condition being non-zero the procedure proposed in [2] cannot be applied in a direct manner. Other difficulty lies in the fact that space approximations are built using standard piece-wise functions, and when this kind of approximation is retained also to construct the time interpolation the known inconsistency related to centered differences of time derivatives is encountered. In [3] we propose alternatives to circumvent both difficulties, and then to allow an efficient treatment of multidimensional transient kinetic theory models.

Until now, only simple rheometric flows were addressed. A first attempt in addressing more complex flows was proposed in [27]. The present paper focuses on the solution of complex fluid models defined in complex geometries. For this purpose we are defining a complex fluid model, a bit more complex than the ones addressed in our former works. The hypothetical fluid consists of a fiber suspension in which the fiber can aggregate. Thus two populations coexist, the one related to the active fibers (making a parallel with the concepts used in associative polymers models [21,15,16]) that participate to the different fibers aggregates, and the other composed by free fibers immersed in the fluid with no direct connection with the other fibers or aggregates. Now, high shear rates are expected destroying the aggregates (therefore increasing the fraction of free tubes), and are reestablished as soon as the shear rate becomes lower enough. Obviously, a such model involves two populations, and then two coupled Fokker-Planck equations. Thus, the results here obtained using a hypothetical model (perhaps unrealistic) could be transferred for addressing more realistic models, as the ones related to associative polymers just referred to.

After describing the model equation related to such hypothetical fluid, we will consider its discretization within the separated representation (finite sums decomposition) framework. The last step will consist in extending this separated representation to cover all the model coordinates (space, configuration and time). This extension will be checked on a contraction flow for the just proposed fluid model.

## 2. A representative micro-macro model of a complex fluid flow

### 2.1. Macroscopic equations

The governing equations to be solved concern the mass and momentum conservation laws for incompressible, isothermal steady state viscoelastic flow. Inertia and body forces are ignored. We may write:

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

$$ \nabla p - \nabla \cdot \underline{\tau} = 2\eta_s \nabla \cdot \underline{D} = 0 \qquad (6) $$

where $\mathbf{v} = (u, v, w)$ denotes the velocity field of the fluid, $p$ is the pressure, $\underline{\tau} = \underline{\tau}_p$ is the extra-stress tensor due to the presence of fibers and $\eta_s$ is the viscosity of the solvent. The extra stress tensor $\underline{\tau} = \underline{\tau}_p$ can be obtained by means of a macroscopic, integral or differential, constitutive equation or using a molecular model involving the resolution of a Fokker-Planck equation. In this work, we consider the microscopic model of fiber network.

### 2.2. Microscopic equations of an hypothetical fiber network model

A network model introduced firstly by Tanaka and Edwards [28,29], is composed of segments and nodes (junctions), Fig. 1. Segment is a microstructure joining two successive nodes and the nodes are the points where the interactions are localized. The internal structure has a transient topology and its dynamics is controlled by the rate of attachment (creation) and detachment (destruction) stochastic and reversible process and depends on the applied deformation.

Network models [15,26,30,21,15,16] incorporate common basic features and physical phenomena. Thus, in our case we consider that a fiber can be active or pendant (dangling). However, the transition from one entity to the other
<figure><img src="image_3.png" /><figcaption>Fig. 1. Fiber network model and kinematic phenomena.</figcaption></figure>

is possible due to flow induced aggregation/desaggregation. We assume that the total number of active and pendant fibers is constant. The flow is considered homogeneous at the scale of fibers independently of the population that it belongs to.

The previous assumptions, seem implying the necessity of considering a set of two Fokker–Planck equations, one for each fiber populations governing the evolution of the respective fiber orientation distributions $\Psi(x, p, t)$ and $\Phi(xp, t)$ for the active and pendant populations, respectively. They describe the fraction of active (respectively pendant) fibers that at point $x$ and time $t$ are oriented in the direction defined by the unit vector $p$. Thus, the orientation distribution function evolution related to the active fibers population $\Psi(x, p, t)$ writes:

$$ \frac{d\Psi}{dt} = -\frac{\partial}{\partial p}\left(\frac{dp}{dt}\Psi\right) + D_{r1}\frac{\partial^2\Psi}{\partial p^2} + V_d\Psi - V_c\Phi \quad (7) $$

where $D_{r1}$ is the diffusion coefficient associated with the active fibers population, $V_d$ and $V_c$ represent the kinetic terms characterizing respectively destruction and aggregation velocities, $(dp/dt)$ denotes the flow induced fiber rotary velocity that is assumed given by the Jeffery's equation:

$$ \frac{dp}{dt} = \underline{\underline{\Omega}} p + \lambda \underline{\underline{D}} p - \lambda (p^T \underline{\underline{D}} p) p \quad (8) $$

where $\underline{\underline{\Omega}}$ and $\underline{\underline{D}}$ are respectively the vorticity and the strain rate tensors associated with the fluid flow and $\lambda$ is a scalar depending on the fiber aspect ratio $r$:

$$ \lambda = \frac{r^2 - 1}{r^2 + 1} \quad (9) $$

with $r = L/d$, $L$ and $d$ referring, respectively, the length and the diameter of the fibers.

Eq. (7) is coupled with the evolution equation for the distribution function of pendant fibers, $\Phi(x, p, t)$:

$$ \frac{d\Psi}{dt} = -\frac{\partial}{\partial p}\left(\frac{dp}{dt}\Phi\right) + D_{r2}\frac{\partial^2\Phi}{\partial p^2} - V_d\Psi + V_c\Phi \quad (10) $$

where $D_{r2}$ denotes the diffusion coefficient related to such fiber population.

The conservation equation involving both population distributions writes:

$$ \int_{\Omega_{\underline{p}}} (\Psi(\underline{x}, t, \underline{p}) + \Phi(\underline{x}, t, \underline{p})) d\Omega_{\underline{p}} = 1, \quad \forall \underline{x}, \forall t \quad (11) $$

$\Omega_{\underline{p}}$ being the fiber configuration space where the fiber orientation is defined, and that corresponds in the 3D case to the surface of the unit sphere (sphere of unit radius). Eq. (11) also implies:

$$ \frac{\partial}{\partial t} \left\{ \int_{\Omega_{\underline{p}}} \left( \Psi(\underline{x}, t, \underline{p}) + \Phi(\underline{x}, t, \underline{p}) \right) d\Omega_{\underline{p}} \right\} = 0, \quad \forall \underline{x}, \forall t \quad (12) $$
Introducing the notation:

$$ \left\{ \begin{aligned} \int_{\Omega_p} \Psi(\underline{x}, t, \underline{p}) \, d\Omega_p &= \Theta(\underline{x}, t) \\ \int_{\Omega_p} \Phi(\underline{x}, t, \underline{p}) \, d\Omega_p &= \Upsilon(\underline{x}, t) \end{aligned} \right. \qquad (13) $$

Eq. (11) writes:

$$ \Theta(\underline{x}, t) + \Upsilon(\underline{x}, t) = 1, \quad \forall \underline{x}, \forall t \qquad (14) $$

We can notice that by adding Eqs. (7) and (10) and integrating over the conformation space $\Omega_p$ we obtain the fiber conservation balance (12).

The integration of the equation Eq. (7) in the configuration space $\Omega_p$ writes:

$$ \frac{d}{dt} \int_{\Omega_p} \Psi \, d\Omega_p = V_d \int_{\Omega_p} \Psi \, d\Omega_p - V_c \int_{\Omega_p} \Phi \, d\Omega_p \qquad (15) $$

that using the notation previously introduced results:

$$ \frac{d\Theta}{dt} = V_d\Theta - V_c(1 - \Theta) \qquad (16) $$

The expression of the extra-stress tensor $\tau_{\underline{p}}$ related to the fibers presence could be derived by generalizing the usual constitutive equation applied for short fiber suspensions:

$$ \tau_{\underline{p}}(\underline{x}, t) = 2\eta_s \left( N_{p1} \int_{\Omega_p} (p \otimes p \otimes p \otimes p) \Psi(\underline{x}, \underline{p}, t) \, d\underline{p} + N_{p2} \int_{\Omega_p} (p \otimes p \otimes p \otimes p) \Phi(\underline{x}, \underline{p}, t) \, d\underline{p} \right) : \underline{D} \qquad (17) $$

where “:” is the tensorial product twice contracted, i.e. the dyadic product.

To solve the associated problem, we consider a decoupled procedure between both the microscopic and the macroscopic scales based on a fixed point strategy (in the steady state regime) or on an explicit time integration (in the transient regime).

Both scales dialog through the extra-stress tensor, calculated at the microscopic scale from both populations distribution functions according to Eq. (17). The extra-stress tensor is then used for updating the suspension flow kinematics (at the macroscopic scale) by solving the resulting anisotropic Stokes problem. The resulting velocity field is then injected to update the orientations distributions by solving (at the microscopic scale) the two coupled Fokker-Planck equations, in which the material derivative (advective terms) makes use of the macroscopic velocity field. The procedure continues until reaching convergence, in steady analysis, or the maximum simulation time (in the transient case).

In what follows we are describing both solvers, the one that serves to solve the two coupled Fokker-Planck equations, and then the flow solver that allows updating the flow kinematics.

### 3. Microscopic scale: a separated representation solver

We proposed in some of our former works (referred in the introduction section) a separated representation technique, making use of finite sums decompositions, that allowed efficient solutions of complex fluid models in simple flows. These models only incorporated a single Fokker-Planck equation eventually defined in high dimensional spaces. The interested reader is invited to refer to [2,3] for additional information on the basic concepts and the implementation details of that novel technique.

The two Fokker-Planck Eqs. (7) and (10) have been strongly coupled. They can be written in a compact form by introducing the vector $\Xi(\underline{p}, t, \underline{x})$:

$$ \Xi(\underline{p}, t, \underline{x}) = \begin{pmatrix} \Psi(\underline{p}, t, \underline{x}) \\ \Phi(\underline{p}, t, \underline{x}) \end{pmatrix} \qquad (18) $$
where we assume that both initial orientation distributions are known at the initial time $t = 0$ defining the vector $\Xi_0$ (where for the sake of clarity the dependence of $\Xi_0$ on the $\underline{p}$ and $\underline{x}$ has been voluntary omitted).

Due to the linearity of the resulting kinetic theory model one could apply a variable transformation in order to define a new couple of unknowns fields subjected to a homogeneous initial condition:

$$ \begin{pmatrix} \Psi(\underline{p}, t, \underline{x}) \\ \Phi(\underline{p}, t, \underline{x}) \end{pmatrix} = \begin{pmatrix} \Psi(\underline{p}, t = 0, \underline{x}) \\ \Phi(\underline{p}, t = 0, \underline{x}) \end{pmatrix} + \begin{pmatrix} \psi(\underline{p}, t, \underline{x}) \\ \phi(\underline{p}, t, \underline{x}) \end{pmatrix} = \underline{\Xi}_0(\underline{p}, \underline{x}) + \underline{\Xi}(\underline{p}, t, \underline{x}) \quad (19) $$

The resulting system can be written in the compact form:

$$ \frac{d\underline{\Xi}(\underline{p}, t, \underline{x})}{dt} + \underline{e}_0(\underline{p}, \underline{x})\underline{\Xi}(\underline{p}, t, \underline{x}) + \underline{e}_1(\underline{p}, \underline{x})\frac{\partial\underline{\Xi}(\underline{p}, t, \underline{x})}{\partial\underline{p}} + \underline{e}_2\frac{\partial^2\underline{\Xi}(\underline{p}, t, \underline{x})}{\partial\underline{p}^2} = -\underline{e}_0(\underline{p}, \underline{x})\underline{\Xi}_0(\underline{p}, \underline{x}) \\ - \underline{e}_1(\underline{p}, \underline{x})\frac{\partial\underline{\Xi}_0(\underline{p}, \underline{x})}{\partial\underline{p}} - \underline{e}_2\frac{\partial^2\underline{\Xi}_0(\underline{p}, \underline{x})}{\partial\underline{p}^2} \quad (20) $$

where

$$ \underline{e}_0(\underline{p}, \underline{x}) = \begin{pmatrix} \frac{\partial \dot{\underline{p}}}{\partial \underline{p}}(\underline{p}, \underline{x}) - V_d & V_c \\ V_d & \frac{\partial \dot{\underline{p}}}{\partial \underline{p}}(\underline{p}, \underline{x}) - V_c \end{pmatrix} \quad (21) $$

$$ \underline{e}_1(\underline{p}) = \dot{\underline{p}} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \quad (22) $$

and

$$ \underline{e}_2 = \begin{pmatrix} -D_{r1} & 0 \\ 0 & -D_{r2} \end{pmatrix} \quad (23) $$

Now, we consider the variational formulation related to Eq. (20) with respect to the conformation coordinates $\underline{p}$ where for the sake of clarity the dependence of all variables on the different coordinates has been omitted:

$$ \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \frac{d\underline{\Xi}}{dt} d\Omega_{\underline{p}} + \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \underline{e}_0 \underline{\Xi} d\Omega_{\underline{p}} + \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \underline{e}_1 \frac{\partial \underline{\Xi}}{\partial \underline{p}} d\Omega_{\underline{p}} + \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \underline{e}_2 \frac{\partial^2 \underline{\Xi}}{\partial \underline{p}^2} d\Omega_{\underline{p}} = - \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \underline{e}_0(\underline{p}) \underline{\Xi}_0 d\Omega_{\underline{p}} \\ - \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \underline{e}_1 \frac{\partial \underline{\Xi}_0}{\partial \underline{p}} d\Omega_{\underline{p}} - \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \underline{e}_2 \frac{\partial^2 \underline{\Xi}_0}{\partial \underline{p}^2} d\Omega_{\underline{p}} \quad (24) $$

where $\underline{\Xi}^*$ represents the usual test functions.

Integrating by parts the last term of the first member, we obtain:

$$ \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \frac{d\underline{\Xi}}{dt} d\Omega_{\underline{p}} + \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \underline{e}_0 \underline{\Xi} d\Omega_{\underline{p}} + \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \underline{e}_1 \frac{\partial \underline{\Xi}}{\partial \underline{p}} d\Omega_{\underline{p}} - \int_{\Omega_{\underline{p}}} \frac{\partial \underline{\Xi}^*}{\partial q} \underline{e}_2 \frac{\partial \underline{\Xi}}{\partial \underline{p}} d\Omega_{\underline{p}} = - \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \underline{e}_0 \underline{\Xi}_0 d\Omega_{\underline{p}} \\ - \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \underline{e}_1 \frac{\partial \underline{\Xi}_0}{\partial \underline{p}} d\Omega_{\underline{p}} - \int_{\Omega_{\underline{p}}} \underline{\Xi}^* \underline{e}_2 \frac{\partial^2 \underline{\Xi}_0}{\partial \underline{p}^2} d\Omega_{\underline{p}} \quad (25) $$

Owing to the advective character of Eq. (25) in the conformation space, a SUPG stabilization is applied on the advective term (the one involving first derivatives with respect to the conformation coordinates):
$$ \int_{\Omega_p} \Xi^* \frac{d\Xi}{dt} d\Omega_p + \int_{\Omega_p} \Xi^* e_0 \Xi d\Omega_p + \int_{\Omega_p} \Xi^* e_1 \frac{\partial \Xi}{\partial p} d\Omega_p + \int_{\Omega_p} \frac{\Gamma h_p}{2} \frac{\partial \Xi^*}{\partial p} e_1 \frac{\partial \Xi}{\partial p} d\Omega_p - \int_{\Omega_p} \frac{\partial \Xi^*}{\partial p} e_2 \frac{\partial \Xi}{\partial p} d\Omega_p \\ = - \int_{\Omega_p} \Xi^* e_0 \Xi_0 d\Omega_p - \int_{\Omega_p} \Xi^* e_1 \frac{\partial \Xi_0}{\partial p} d\Omega_p - \int_{\Omega_p} \Xi^* e_2 \frac{\partial^2 \Xi_0}{\partial p^2} d\Omega_p \quad (26) $$

where $\Gamma$ is calculated according to:

$$ \Gamma = \max \left( \left( \coth(P_{e1}) - \frac{1}{P_{e1}} \right), \left( \coth(P_{e2}) - \frac{1}{P_{e2}} \right) \right) \quad (27) $$

and the Peclet numbers $P_{e1}$ and $P_{e2}$ are given by:

$$ \left\{ \begin{array}{l} P_{e1} = \frac{\|e_1(\underline{p}, \underline{x})\| h_p}{2D_{r1}} \\ P_{e2} = \frac{\|e_1(\underline{p}, \underline{x})\| h_p}{2D_{r2}} \end{array} \right. \qquad (28) $$

where $h_p$ represents the characteristic size of the conformation space discretization.

In what follows we are considering the steady state solution of a complex flow. For this purpose we are coupling the separated representation solver just described with the Stream-Tube method (summarized in the next section) for computing the flow kinematics. The Fokker-Planck equation will be integrated along the flow streamlines. For this purpose we use the method of characteristics. In this context a point $\underline{x}$ located on a streamline can be represented as $\underline{x}(t)$. Thus, in what follows we consider that both distribution functions write: $\psi(\underline{p}, t)$ and $\phi(\underline{p}, t)$. Thus, the solution representation in the conformation space is searched in the separated form:

$$ \Xi = \begin{pmatrix} \psi(\underline{p}, t) \\ \phi(\underline{p}, t) \end{pmatrix} = \sum_{j=1}^{n} \begin{pmatrix} \alpha_j E_j(\underline{p}) F_j(t) \\ \beta_j G_j(\underline{p}) H_j(t) \end{pmatrix} \qquad (29) $$

where again, for the sake of clarity, the dependence of variables on the physical space $\underline{x}$ is omitted.

This solution representation is built-up within an iteration scheme. Thus, if at iteration $n$ the solution involves the known functional couples $E_j(\underline{p})F_j(t)$ and $G_j(\underline{p})H_j(t)$ (see Eq. (29)), then we perform a projection, an enrichment step and a convergence check:

(1) **Projection**: By introducing the representation (29) into the stabilized variational formulation one could compute (after discretization) the approximation coefficients $\alpha_i$ and $\beta_i$, $i \in [1, \dots, n]$.

(2) **Enrichment**: With the just computed approximation (29) where both the functional couples and the approximation coefficients are known, we look for the new functional couples: $R(\underline{p})S(t)$ and $V(\underline{p})W(t)$ by introducing the enriched approximation

$$ \Xi = \begin{pmatrix} \psi(\underline{p}, t) \\ \phi(\underline{p}, t) \end{pmatrix} = \sum_{j=1}^{n} \begin{pmatrix} \alpha_j E_j(\underline{p}) F_j(t) \\ \beta_j G_j(\underline{p}) H_j(t) \end{pmatrix} + \begin{pmatrix} R(\underline{p})S(t) \\ V(\underline{p})W(t) \end{pmatrix} \qquad (30) $$

again into the variational formulation. The enrichment stage results (after discretization and numerical integration) in a non-linear problem that is solved using an appropriate iteration procedure. This question is addressed later. After convergence of the non-linear solver we obtain the expression of both functional couples that after
normalization, i.e.

$$ \left\{
\begin{aligned}
E_{n+1}(\underline{p}) &= \frac{R(\underline{p})}{||R(\underline{p})||} \\
F_{n+1}(t) &= \frac{S(t)}{||S(t)||} \\
G_{n+1}(\underline{p}) &= \frac{V(\underline{p})}{||V(\underline{p})||} \\
H_{n+1}(t) &= \frac{W(t)}{||W(t)||}
\end{aligned}
\right.
\qquad (31) $$

leads to the updated separated representation:

$$ \Xi = \begin{pmatrix} \psi(\underline{p}, t) \\ \phi(\underline{p}, t) \end{pmatrix} = \sum_{j=1}^{n+1} \begin{pmatrix} \alpha_j E_j(\underline{p}) F_j(t) \\ \beta_j G_j(\underline{p}) H_j(t) \end{pmatrix} \qquad (32) $$

(3) Checking the convergence: This iteration scheme needs the definition of a stoping criterion. Many possibilities exist, being the simplest ones:

* The iteration process stops as soon as $||R(\underline{p})S(t)||$ and $||V(\underline{p})W(t)||$ become small enough. These norms can be evaluated very fast, but even when both become small enough nothing proves that the solution is reached. Only in some cases we proved that these norms decrease monotonically.

* The other possibility consists in evaluation the residual of the strong formulation of the problem, Eq. (20). When the residual becomes small enough the model solution is attained. However, the computational cost of this residual is higher than the one needed for computing the norms of the enrichment functions just introduced.

In what follows we use the stoping criterion based on the residual of the strong formulation. The iteration at which the convergence is attained defined the number of sums involved in the finite sums decomposition.

The solution procedure needs firstly a discrete representation of all the fields involved in the variational formulation. The associated notation is grouped in the first appendix. These discrete representations allow us to transform the variational formulation (26) into a matrix form as described in the second appendix for the projection step and in the third appendix for the approximation basis enrichment step.

Here, we would like only to address the main benefits that such separated representation introduces:

* The projection step (at iteration $n$) needs the solution of a linear system consisting of $2n$ unknowns (the coefficients $\alpha_i$ and $\beta_i$ in the expression (29). This system is fully populated, but it has a very reduced size (in the order of tens in all the simulation performed until now).

* The enriching step involves the solution of a non-linear problem, because now the searched terms are product of functions. We use an alternating directions scheme (see [3]) in which the functions defined in the conformation space $\underline{p}$ are searched assuming known the ones related to the time coordinate (that were computed at the previous iteration). Thus, from the just updated function defined in the conformation space $\underline{p}$, the time function can be updated. The iteration procedure continues until reaching convergence (other more sophisticated procedures – e.g. Newton-Raphson – exist but its application does not increase significantly the solver efficiency). The solution of the non-linear enrichment problem is attained in few iterations, that will be noted by $P$. Thus, if we are using $N_p$ nodes in the conformation space and $N_t$ time-nodes, and even if we are computing the whole space-time evolution, we must solve linear systems of size $2N_p$ (related to the discretization of the conformation domain $\Omega_p$ – the unit sphere surface) and 1D backward discrete models related to the time evolution.

* This procedure, that seems to be specially appropriate for solving multidimensional models because here the complexity scales linearly with the dimension of the space (instead the exponential scaling of mesh based discretization strategies), allows even in the case just described important CPU time savings. It allows also a simple treatment of
non-linear kinetic theory models, where at iteration n, the non-linear terms could be evaluated at the solution defined at the previous iteration.

### 3.1. Addressing complex flows

The previous analysis also useful for the treatment of transient models with uniform solutions in the physical space $\Omega_x$. These kind of flows are usually encountered in rheometric devices. However, in complex flows, the evolution of the different fields with respect to the physical coordinates $\underline{x}$ cannot be neglected. Thus, we must address the discretization of the material derivative $d\underline{\Xi}/dt$ involved in the weak form of the Fokker-Planck Eq. (26).

In this more general situation, one could define a general separated representation consisting of a finite sum decomposition in which each term of the sum involves the product of three functions: a function of the physical space $\underline{x}$, a function of the conformation space $\underline{p}$ and finally a function of time $t$. This kind of representation needs for an efficient stabilization of the advection term that applies in the physical space $\Omega_x$. This route needs for further analysis and developments. Thus, in what follows we are focussing on steady complex flows that as we are proving can profit of all the analysis introduced above.

If we have a steady state flow, one could integrate the Fokker-Planck equation along the characteristic lines (that coincide in this case with the flow streamlines) from the initial condition assumed known on the inflow boundary. Thus, one could expect that the solution at each point $P$ only depends on curvilinear coordinate $s$ associated with that point $s_P$. We assume that the streamline that passes through point $P$ has its origin $s = 0$ at the point in which the streamline intersects the domain inflow boundary. The stabilization by upwinding used for the treatment of the time discretization is automatically applicable to the discretization on the curvilinear coordinate.

With these ideas in mind, one could write, for the orientation distribution field $\underline{\Xi}$ defined in $\Omega_x \times \Omega_p$ the following separated representation:

$$ \underline{\Xi}_{st} = \begin{pmatrix} \psi_{st}(\underline{p}, s) \\ \phi_{st}(\underline{p}, s) \end{pmatrix} = \sum_{j=1}^{n} \begin{pmatrix} \alpha_j^{st} E_j^{st}(\underline{p}) F_j^{st}(s) \\ \beta_j^{st} G_j^{st}(\underline{p}) H_j^{st}(s) \end{pmatrix} \qquad (33) $$

where the index $st$ refers to the particular streamline along which the integration is performed.

If the integration is performed on a large enough number of streamlines, the solution everywhere can be obtained from the interpolation of solutions defined on the streamlines.

## 4. Macroscopic scale: flow kinematics solver

In order to facilitate the coupling between the microscopic solver and the macroscopic one, we propose the use of the stream tube method (deeply described in [6,10–14,23,24]) for solving the resulting macroscopic anisotropic Stokes problem. This technique is based on a domain transformation that maps the physical domain into a simple reference domain in which the streamlines are parallel to one of the coordinate axis. Thus, instead of looking for the velocity and pressure fields (as in usual mixed formulations) one looks for the domain transformation function and the pressure field.

We are summarizing the main ideas that this technique involves, but the interested reader can refer to the papers just addressed (some o them cited in the reference section) where more details in computational implementation are introduced and where several complex flows are addressed (including some ones involving primary and secondary recirculating areas). However, until now this technique was only applied for different macroscopic non-Newtonian constitutive equations, but never coupled with kinetic theory descriptions of complex fluids.

The Stream-Tube method can be directly and easily coupled with the separated representation of the Fokker-Planck equation on the streamlines, described in Section 3.1.

### 4.1. The Stream-Tube method revisited: basic concepts

We consider here a two-dimensional flow. In Cartesian coordinates, $k$ denotes the mapping function, between the respective reference and physical domains $\Omega^*$ and $\Omega$ (see Fig. 2). Under the assumption of a flow involving only open streamlines (for more complex flow condition the interested reader can refer to [10,23] and the references therein), we
<figure><img src="image_4.png" /><figcaption>Fig. 2. Stream-Tube method: domain transformation.</figcaption></figure>

may write the following relations:

$$ \begin{cases} x = k(X, Z) \\ z = Z \end{cases} \qquad (34) $$

For the sake of clarity, we introduce the notation $(\partial x/\partial X) = k_{,X}$ and $(\partial x/\partial Z) = k_{,Z}$. The Jacobian of the transformation writes:

$$ J = \begin{pmatrix} k_{,X} & k_{,Z} \\ 0 & 1 \end{pmatrix} \qquad (35) $$

being the Jacobian of the inverse transformation

$$ \begin{pmatrix} \frac{\partial X}{\partial x} & \frac{\partial X}{\partial z} \\ \frac{\partial Z}{\partial x} & \frac{\partial Z}{\partial z} \end{pmatrix} = \begin{pmatrix} k_{,X} & k_{,Z} \\ 0 & 1 \end{pmatrix}^{-1} = \begin{pmatrix} \frac{1}{k_{,X}} & -\frac{k_{,Z}}{k_{,X}} \\ 0 & 1 \end{pmatrix} \qquad (36) $$

The differential operators are transformed according to:

$$ \begin{cases} \frac{\partial}{\partial x} = \frac{1}{k_{,X}} \frac{\partial}{\partial X} \\ \frac{\partial}{\partial x} = -\frac{k_{,Z}}{k_{,X}} \frac{\partial}{\partial X} + \frac{\partial}{\partial X} \end{cases} \qquad (37) $$

In the two-dimensional case here addressed the components of the velocity vector $v^T = (u, w)$ are transformed in $v^T = (U, W)$ when the coordinate system $(X, Z)$ is considered instead of the physical one $(x, y)$. The just referred transformation writes:

$$ \begin{pmatrix} U \\ W \end{pmatrix} = \begin{pmatrix} \frac{1}{k_{,X}} & -\frac{k_{,Z}}{k_{,X}} \\ 0 & 1 \end{pmatrix} \begin{pmatrix} u \\ w \end{pmatrix} = \begin{pmatrix} \frac{1}{k_{,X}} u - \frac{k_{,Z}}{k_{,X}} w \\ w \end{pmatrix} \qquad (38) $$
Now, if we enforce

$$u = \frac{\partial x}{\partial Z} w = k, Z w \qquad (39)$$

the X-component of the velocity vector $U$ vanishes, fact that guarantees that the flow streamlines consist of the lines $X = cte$.

Solving the problem by using the Stream-Tube method consists of determining the mapping function that minimizes the residual of the conservation equations, by solving the resulting non-linear model using an appropriate technique (e.g. the Newton-Raphson strategy). In the two-dimensional case here considered, the momentum balance (neglecting mass and inertia terms) writes

$$\nabla \cdot (\underline{\tau} - p\underline{I}) = 0 \qquad (40)$$

with $\tau = 2\eta_s D + \underline{\tau}_p$. Using the just introduced transformation, the momentum balance equations read:

$$\left\{ \begin{aligned} 0 &= \frac{\partial \tau_{xx}}{\partial X} + k, X \frac{\partial \tau_{xz}}{\partial Z} - k, Z \frac{\partial \tau_{xz}}{\partial X} - \frac{\partial p}{\partial X} \\ 0 &= \frac{\partial \tau_{xz}}{\partial X} + k, X \frac{\partial \tau_{zz}}{\partial Z} - k, Z \frac{\partial \tau_{zz}}{\partial X} - k, X \frac{\partial p}{\partial Z} + k, Z \frac{\partial p}{\partial X} \end{aligned} \right. \qquad (41)$$

Finally, we should address the mass conservation balance for an incompressible flow. For this purpose we consider two neighbor streamlines, closer enough to suppose that the velocity between both streamlines does not evolve significantly. Now, we consider the inflow boundary $z = 0$ and another generic section $z$ within the flow domain. These two sections result in the transformed sections $Z = 0$ and $Z$. At the inflow section $dX \equiv dx$, that is, $k, X|_{Z=0} = 1$ and then $dX = dx$. At section $Z$, $dX = k, X dx$, and the mass conservation writes:

$$w_{X,Z=0} dx = w(X, Z) k, X|_{X,Z} dx \rightarrow w(X, Z) = w(X, Z=0) \frac{1}{k, X|_{X,Z}} \qquad (42)$$

The classical “velocity/pressure” formulation is then replaced by “transformation function/pressure” formulation. To solve the related equation, we propose an iterative algorithm that allows to compute the “transformation function/pressure” and the stress resulting from the microscopic calculations within a fixed point strategy. For a given value of microscopic stress at the previous iteration $l-1$, the “transformation function/pressure” derived from the momentum balance equation:

$$-(\nabla p)^l + (2\eta_s \nabla \cdot \underline{D})^l = -(\nabla \cdot \underline{\tau}_p)^{l-1} \qquad (43)$$

where $\eta_s$ is the solvent viscosity and the strain rate tensor $\underline{D}$ is expressed in terms of the unknown transformation function and its derivatives with respect to $X$ and $Z$ using the relationship (37) and the development of the divergence operator given in Eq. (41). The superscripts “l” and “l-1” denote, respectively, the current and the previous iteration.

Since in this equation the elliptic contribution related to the solvent viscosity can be negligible against the stress contribution, the introduction of an arbitrary numerical viscosity must be done. This viscosity is chosen to be in the same order of the apparent viscosity of the microscopic behavior (that decreases for high values of shear rates). A Carreau model has been arbitrarily chosen here to represent this arbitrary viscosity. Thus, Eq. (43) is written as:

$$-(\nabla p)^l + (2(\eta_s + \eta_{carr})\nabla \cdot \underline{D})^l = -(\nabla \cdot \underline{\tau}_p)^{l-1} + (2\eta_{carr}\nabla \cdot \underline{D})^{l-1} \qquad (44)$$

with the arbitrary viscosity given by:

$$\eta_{carr} = \frac{\eta_0}{[1 + (a\dot{\gamma})^2]^n} \qquad (45)$$

where $\eta_0$, $a$ and $n$ are three model parameters that are adjusted in order to describe with a Carreau model the rheological behavior that the considered fluid exhibits a simple shear flow.

The boundary conditions are transformed using all the previous expressions, resulting in different type of conditions on the transformation function $k$. The interested reader can refer to [10,23] and the references therein.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Steady distribution function related to the fiber network model in a 3D transient state shear flow with $D_{r1} = 0.5$, $D_{r2} = 0.2$, $V_c = 0.1$ and $V_d = 0.1$: (top-left) active fibers for $W_e = 1$, (top-right) pendant fibers for $W_e = 1$, (bottom-left) active fibers for $W_e = 10$, (bottom-right) pendant fibers for $W_e = 10$.</figcaption></figure>

# 5. Numerical results

## 5.1. Transient network analysis in a steady simple shear flow

First, we consider the resolution of the Fokker-Planck equation related to the fiber network model. Fig. 3 depicts the distribution function for active and pendant fibers obtained by considering $D_{r1} = 0.5$, $D_{r2} = 0.2$, $2\eta_s N_{p1} = 1$, $2\eta_s N_{p2} = 1$, $V_c = 0.1$ and $V_d = 0.1$. An isotropic orientation of fiber distribution is assumed as initial orientation state. A mesh of 642 nodes of the unit sphere related to the configuration space $\Omega_p$ and 100 nodes on the 1D domain $\Omega_t$ related to the time coordinate. After imposing the shear flow, the fiber orientation concentrates around the direction of the flow. The computed solutions for different Weissenberg numbers ($W_e = \max((\dot{y}/D_{r1}), (\dot{y}/D_{r2}))$) in simple 3D shear flows defined by ($u = \dot{y}y, v = 0$), are compared. Under these conditions, we can notice that the active fibers have a considerable contribution in comparison with the pendant fibers. This contribution increases with the Weissenberg number.

Figs. 4 and 5 depict the most significant approximation functions, related to the conformation space and time, involved in the separated representation of the distribution function for both $W_e = 1$ and $W_e = 10$, respectively. The dimensionless shear stresses for fiber network model is depicted in Fig. 6. It was pointed out that these results confirm a classical profile of the polymer stress translating a rheological behavior of such models.

Fig. 7 depicts the time evolution of the active population in a simple shear flow characterized by $W_e = 0.1$. This Obviously, the pendant population constitutes the complementary part, such that the addition of both becomes constant during the entire simulation. Moreover, we can notice that independently of the initial concentration of active fibres (40%, 60% and 80% in our simulations), the same steady state is reached.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Most significant approximation functions involved in the separated representation of the distribution function related to the fiber network model in a 3D transient shear flow characterized by $W_e = 1$, $D_{r1} = 0.5$, $D_{r2} = 0.2$, $V_c = 0.1$ and $V_d = 0.1$: (left) active fibers, (right) pendant fibers.</figcaption></figure>
<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" /></figure>

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

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

<figure><img src="image_13.png" /><figcaption>Fig. 5. Most significant approximation functions involved in the separated representation of the distribution function related to the fiber network model in a 3D transient shear flow characterized by $W_e = 10$, $D_{r1} = 0.5$, $D_{r2} = 0.2$, $V_c = 0.1$ and $V_d = 0.1$: (left) active fibers, (right) pendant fibers.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 6. Dimensionless microscopic shear stress evolution related to the fiber network model: (top) $D_{r1} = 0.5$ and $D_{r2} = 0.2$; (bottom) $V_c = V_d = 0.1$.</figcaption></figure>

Finally, Fig. 8 compares the evolution of active populations for two different Weissenberg numbers $W_e = 0.1$ and $W_e = 1$. We can notice that the higher is the shear rate ($W_e$) the higher is the time evolution of both populations.

## 5.2. Analysis of a contraction flow

The algorithm consists in repeating until convergence the following three steps:

<figure><img src="image_15.png" /><figcaption>Fig. 7. Time evolution of the active population of fibres in a simple shear flow at $W_e = 0.1$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 8. Time evolution of active populations in a simple shear flow for $W_e = 0.1$ (top) and $W_e = 1$ (bottom).</figcaption></figure>

(1) Compute the microscopic stress from the orientation distributions obtained by solving the Fokker-Planck equation.

(2) Compute the flow kinematics.

(3) Check the convergence (of the fixed point strategy).

This algorithm will be applied on a physical domain involving a contraction. In this section, we present the macroscopic results. Fig. 9 plots the axial velocity component $w$, the stress components $\tau_{xz}$ and $\tau_{zz}$ and the pressure field. It can be noticed the velocity overshoot just downstream of the convergent section, the peaks of the stress near the re-entrant corner and the variation of the total pressure drop with elastic effects.

These macroscopic results reveal a variation of the analyzed field near the re-entrant corner. We can conclude that these expected macroscopic behavior is a consequence of the fiber orientation at the microscopic level. To illustrate this dependence, we depict in Fig. 10, the distribution function along some streamlines. Indeed, at the central streamline the microscopic state shows the existence of only the elongational effects are activated. For the other streamlines the distribution functions become a compromise between elongation and shear.

### 5.3. Convergence analysis

The fluid model considered in this paper, that we solved in a complex geometry in the previous section, does not correspond to any physical model. We proposed this model because it involves two coupled Fokker-Planck equations governing the evolution of two populations in competition. The resulting model is quite similar to kinetic theory models of associative polymers. Numerical solutions of kinetic theory models, making use of the Fokker-Planck formalism and defined in complex geometries, are quite rare. By this reason we cannot compare our results with others previously published.
<figure><img src="image_17.png" /><figcaption>Fig. 9. (Top-left) Dimensionless velocity component $w$, (top-right) dimensionless stress component $\tau_{xz}$, (bottom-left) dimensionless stress component $\tau_{zz}$ and (bottom-right) dimensionless pressure $p$.</figcaption></figure>

However, to be confident with the results reported in the previous sections, a convergence analysis should be performed.

First, a simpler model, as the one that governs the evolution of the orientation of a single population of fibers, was addressed in [27]. In that paper, the Fokker-Planck equation governing the evolution of the orientation of the fibers suspended in a Newtonian fluid was solved and coupled with the flow kinematics computed by using finite elements. The integration of the Fokker-Planck equation along the flow streamlines (as we consider in the present paper) produced the best results, in comparison with the ones obtained by using finite elements (discontinuous Galerkin), in some flows where the exact solutions were known.

In what follows we are discussing first, the convergence of the kinetic-theory solver, then the convergence of the Stream-Tube method, and finally the coupling of both.

Concerning the convergence of the kinetic theory solver we considered 4 different meshes of the unit surface consisting of: MM2 (42 nodes), MM3 (162 nodes), MM4 (642 nodes) and MM5 (2562 nodes) illustrated in Fig. 11. The kinetic theory model was solved in a simple shear flow by using each one of these meshes without modifying the time step. The parameters involved in the kinetic theory model were set to: $D_{r1} = 0.5$, $D_{r2} = 0.1$ and $V_c = V_d = 0.1$. Fig. 12 depicts the computed stress for these four meshes as a function of the Weissenberg number. The convergence of the numerical solution with the mesh refinement can be noticed.

Concerning the convergence of the Stream-Tube method it was analyzed in many former works of some of the authors of the present work. In the numerical simulations made by the Stream-Tube method, mesh convergence tests have been performed in order to validate the computed solutions. Those have been reported in [11]. Moreover, in other works the predictions of the Stream-Tube method were compared with the solutions obtained with other alternative discretization techniques (see [6,13]) and even with commercial simulation codes (Polyflow and Fluent) [14].
<figure><img src="image_18.png" /><figcaption>Fig. 10. Micro-macro coupling for $h = -1$ and $\eta_s = 1$: distribution functions for active, $\Psi$, and pendant, $\Phi$, populations of fibers along the streamlines.</figcaption></figure>

Finally, in order to check the convergence of the coupling strategy we are considering the contraction flow and the two populations kinetic theory model defined by the set of parameters $D_{r1} = 0.5$, $D_{r2} = 0.2$ and $V_c = V_d = 0.1$. The integration of the Fokker-Planck equation along the flow streamlines is performed by decoupling, within the separated representation framework, the conformational coordinates $p$ and a coordinate mapping the streamline. To parameterize the streamline we are using a pseudo-time, as is usual in the method of characteristics, where any point on the streamline can be defined from $\underline{x}(t)$, where $t$ represents the time that a fluid particle spent for reaching position $\underline{x}$. Different time steps were considered to represent the different temporal functions involved in the separated representation: $\Delta t = 0.1, 0.05, 0.02, 0.005$ and $\Delta t = 0.001$. Fig. 13 shows the computed $a_{1212}$ component of the fourth order orientation tensor associated with the first population using the different time steps on a streamline segment $[\underline{x}(t=0), \underline{x}(t=0.5)]$. The numerical convergence can be easily noticed in that figure.

The interest of using a separated $(\underline{x}(t), \underline{p})$ representation, instead a standard backward integration along the streamline, allows for impressive CPU time savings, that can attain some orders of magnitude depending on the considered problem.
<figure><img src="image_19.png" /><figcaption>Fig. 11. Different meshes of the unit sphere surface containing 42, 162, 642 and 2562 nodes.</figcaption></figure>

If one proceeds by applying a standard integration by the method of characteristics, a problem defined on the unit sphere surface should be solved at each time step. In the contraction flow previously addressed, standard integration schemes need the solution of thousands of problems defined on the surface of the unit sphere, along each streamline. However, when one uses separated representations, if the solution consists of a finite sum of *N* terms, and *P* iterations were needed for computing each term of the sum, one must solve *N* × *P* problems defined on the surface of the unit sphere and other *N* × *P* defined on the time axis (or the curvilinear abscise related to the considered streamline). In the example here considered, *N* ≈ 10 and *P* ≈ 10 and then the computing time savings are of one order of magnitude. However, by decreasing the time step (needed when the Weissenberg number increases) the computing time involved in the separated representation remains more or less unchanged because only the 1D integration is affected, but this integration is extremely fast. On the contrary, when one proceeds using an incremental technique the computing time increases linearly with the time step reduction. Finally, we must mention that the stability constraints are much less critical in the case of separated representations that in usual incremental techniques. We have computed solutions using larger time steps that the ones allowed by the stability constraint of incremental procedures, without stability issues. The same time steps induces stability issues when they are used within incremental strategies.
<figure><img src="image_20.png" /><figcaption>Fig. 12. Stress vs. Weissenberg number for four different meshes.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 13. $a_{1212}$ along a streamline segment for different integration time steps.</figcaption></figure>

# 6. Conclusions

In this paper we proposed a fully deterministic micro-macro approach for simulation complex fluid flows. This approach performs, at the microscopic scale, a separated representation of the fields involved in the kinetic theory representation to circumvent the curse of dimensionality that kinetic theory models usually involve.

At the macroscopic level, and for general complex flows, the microscopic information can be easily taken into account by integrating the kinetic theory model along the flow streamlines (in steady flows) or along the nodal pathlines in the transient case.

The Stream-Tube method allows a simple description of the flow streamlines, allowing simple microscopic integrations. For this reason it has been adopted in the present paper when a steady-state geometrically complex flow was considered.

For general transient and geometrically complex flows one must integrate the microscopic model along the flow pathlines, that are automatically defined when a Lagrangian description of the flow kinematics is adopted. Obviously, this kind of descriptions induces some numerical difficulties when one proceeds in the finite element framework because the high element distortion. However, the use of meshless techniques could be an appealing choice for addressing this kind of micro-macro models involving complex fluid models, complex geometries and transient analysis. This analysis constitutes a work in progress.
# Appendix 1. Separated representation solver: notation

Using the notation introduced in Section 3 we start defining the different discrete fields. Thus, functions $E_j(p)$ (resp. $F_j(t)$, $G_j(p)$, $H_j(t)$) and $R(p)$ (respectively $S(t)$, $V(p)$ and $W(t)$) are defined using usual finite element interpolations. The conformation vector $p$ related to a node on the unit sphere $p_i$ is defined as

$$ \underline{p}_i^T = (x_i, y_i, z_i) \qquad (46) $$

where $x_i^2 + y_i^2 + z_i^2 = 1$.

From now on, we denote by $\underline{N}$ (resp. $\underline{M}$) the vector containing the shape functions $N_r(p)$ (resp. $M_s(t)$). Thus, any function depending on $p$, e.g. $E_j(p)$, can be written as:

$$ E_j(\underline{p}) = \sum_{r=1}^{N_r^p} N_r(\underline{p}) \cdot E_j(\underline{p}_r) \qquad (47) $$

being $N_r^p$ the number of nodes distributed on the surface of the unit sphere, and $N_r(p)$ the shape function related to node $p$.

In our simulations the unit sphere is substituted by polyhedra composed by a set of triangular facets whose vertices are the nodes $p$ distributed on the unit surface. In each triangle the shape functions are the usual piecewise linear shape functions used in the context of finite elements.

Finally $\underline{E}_j$, $\underline{F}_j$, $\underline{G}_j$, $\underline{H}_j$, $\underline{R}$, $\underline{S}$, $\underline{V}$ and $\underline{W}$ define the vector containing the respective nodal values.

We define the following matrix:

$$ \underline{N} = \int_{\Omega_p} \underline{N} \underline{N}^T d\Omega_p \underline{M} = \int_{\Omega_t} \underline{M} \underline{M}^T d\Omega_t \qquad (48) $$

$$ \underline{N}_{11} = \int_{\Omega_p} \left\{ (\underline{e}_0)_{1,1} \underline{N} \underline{N}^T + (\underline{e}_1)_{1,1} \underline{N} \underline{d} \underline{N}^T + \frac{\Gamma \underline{h}_p}{2} (\underline{e}_1)_{1,1} \underline{d} \underline{N} \underline{d} \underline{N}^T - (\underline{e}_2)_{1,1} \underline{d} \underline{N} \underline{d} \underline{N}^T \right\} d\Omega_p \qquad (49) $$

$$ \underline{N}_{12} = \int_{\Omega_p} \{ (\underline{e}_0)_{1,2} \underline{N} \underline{N}^T \} d\Omega_p \qquad (50) $$

$$ \underline{N}_{21} = \int_{\Omega_p} \{ (\underline{e}_0)_{2,1} \underline{N} \underline{N}^T \} d\Omega_p \qquad (51) $$

$$ \underline{N}_{22} = \int_{\Omega_p} \left\{ (\underline{e}_0)_{2,2} \underline{N} \underline{N}^T + (\underline{e}_1)_{2,2} \underline{N} \underline{d} \underline{N}^T + \frac{\Gamma \underline{h}_p}{2} (\underline{e}_1)_{2,2} (\underline{p}) \underline{d} \underline{N} \underline{d} \underline{N}^T - (\underline{e}_2)_{2,2} \underline{d} \underline{N} \underline{d} \underline{N}^T \right\} d\Omega_p \qquad (52) $$

and

$$ \underline{M} = \int_{\Omega_t} \underline{M} \underline{d} \underline{M}^T d\Omega_t \qquad (53) $$

# Appendix 2. Separated representation solver: projection step

If we assume the functions $E_j$, $F_j$, $G_j$ and $H_j$ ($\forall j \in [1, \dots, n]$) to be known (verifying the boundary conditions), the coefficients $\alpha_j$ and $\beta_j$ can be computed by introducing the approximation of $\Xi$ into the Galerkin variational formulation associated with Eq. (26).

At the $n$ th iteration, the components of $\Xi$ yield (using the notation introduced in the previous appendix):

$$ \psi(\underline{p}, t) = [\underline{N}^T \underline{E}_1 \underline{M}^T \underline{F}_1 \dots \underline{N}^T \underline{E}_n \underline{M}^T \underline{F}_n] \cdot \begin{bmatrix} \alpha_1 \\ \vdots \\ \alpha_n \end{bmatrix} \qquad (54) $$
and

$$ \phi(p, t) = [\underline{\underline{N}}^T \underline{G}_1 \underline{M}^T \underline{H}_1 \dots \underline{\underline{N}}^T \underline{G}_n \underline{M}^T \underline{H}_n] \cdot \begin{bmatrix} \beta_1 \\ \vdots \\ \beta_n \end{bmatrix} \quad (55) $$

Then, the components of the test field $\underline{\underline{\Xi}}^*$, which is involved in the variational formulation, are given by:

$$ \psi^*(\underline{p}, t) = [\underline{\underline{N}}^T \underline{E}_1 \underline{M}^T \underline{F}_1 \dots \underline{\underline{N}}^T \underline{E}_n \underline{M}^T \underline{F}_n] \cdot \begin{bmatrix} \alpha_1^* \\ \vdots \\ \alpha_n^* \end{bmatrix} \quad (56) $$

and

$$ \phi^*(\underline{p}, t) = [\underline{\underline{N}}^T \underline{G}_1 \underline{M}^T \underline{H}_1 \dots \underline{\underline{N}}^T \underline{G}_n \underline{M}^T \underline{H}_n] \cdot \begin{bmatrix} \beta_1^* \\ \vdots \\ \beta_n^* \end{bmatrix} \quad (57) $$

Now, if all these expressions are introduced into the weak form of the Fokker-Planck equation

$$ \int_{\Omega} \underline{\underline{\Xi}}^* \frac{d\underline{\underline{\Xi}}}{dt} d\Omega + \int_{\Omega} \underline{\underline{\Xi}}^* e_{\underline{0}} \underline{\underline{\Psi}} d\Omega + \int_{\Omega} \underline{\underline{\Xi}}^* e_1 \frac{\partial \underline{\underline{\Xi}}}{\partial \underline{p}} d\Omega + \int_{\Omega} \frac{\Gamma h_p}{2} \frac{\partial \underline{\underline{\Xi}}^*}{\partial \underline{p}} e_1 \frac{\partial \underline{\underline{\Xi}}}{\partial \underline{p}} d\Omega \\ - \int_{\Omega} \frac{\partial \underline{\underline{\Xi}}^*}{\partial \underline{p}} e_2 \frac{\partial \underline{\underline{\Xi}}}{\partial \underline{p}} d\Omega = - \int_{\Omega} \underline{\underline{\Xi}}^* e_{\underline{0}} \underline{\underline{\Xi}}_0 d\Omega - \int_{\Omega} \underline{\underline{\Xi}}^* e_1 \frac{\partial \underline{\underline{\Xi}}_0}{\partial \underline{p}} d\Omega - \int_{\Omega} \underline{\underline{\Xi}}^* e_2 \frac{\partial^2 \underline{\underline{\Xi}}_0}{\partial \underline{p}^2} d\Omega \quad (58) $$

the following linear system, from which we can determine the $\alpha_i$ and $\beta_i$ coefficients $\forall i = 1, \dots, n$, is obtained:

$$ \underline{\underline{K}} \underline{\underline{\alpha}} = \underline{\underline{V}} \quad (59) $$

where

$$ \underline{\underline{\alpha}} = [\alpha_1 \ \dots \ \alpha_n \ \beta_1 \ \dots \ \beta_n]^T \quad (60) $$

and the components of $\underline{\underline{K}}$ and $\underline{\underline{V}}$ are given by:

$$ \underline{\underline{K}}_{i,j} = \left[ \begin{pmatrix} \underline{\underline{E}}_i^T \\ \underline{\underline{G}}_i^T \end{pmatrix} \left( \frac{\underline{\underline{N}}_{11} + \underline{\underline{N}}_{12}}{\underline{\underline{N}}_{21} + \underline{\underline{N}}_{22}} \right) \begin{pmatrix} \underline{\underline{E}}_j \\ \underline{\underline{G}}_j \end{pmatrix} \cdot \begin{pmatrix} \underline{\underline{F}}_i^T \\ \underline{\underline{H}}_i^T \end{pmatrix} \begin{pmatrix} \underline{\underline{M}} \\ \underline{\underline{M}} \end{pmatrix} \begin{pmatrix} \underline{\underline{F}}_j \\ \underline{\underline{H}}_j \end{pmatrix} \right] \\ + \left[ \begin{pmatrix} \underline{\underline{E}}_i^T \\ \underline{\underline{G}}_i^T \end{pmatrix} \begin{pmatrix} \underline{\underline{N}} \\ \underline{\underline{N}} \end{pmatrix} \begin{pmatrix} \underline{\underline{E}}_j \\ \underline{\underline{G}}_j \end{pmatrix} \cdot \begin{pmatrix} \underline{\underline{F}}_i^T \\ \underline{\underline{H}}_i^T \end{pmatrix} \begin{pmatrix} \underline{\underline{M}} \\ \underline{\underline{M}} \end{pmatrix} \begin{pmatrix} \underline{\underline{F}}_j \\ \underline{\underline{H}}_j \end{pmatrix} \right] \quad (61) $$

and

$$ \underline{\underline{V}}_i = - \begin{pmatrix} \underline{\underline{E}}_i^T \\ \underline{\underline{G}}_i^T \end{pmatrix} \int_{\underline{\underline{\Omega}}_p} \left\{ \left( \frac{\underline{\underline{N}}_{11} + \underline{\underline{N}}_{12}}{\underline{\underline{N}}_{21} + \underline{\underline{N}}_{22}} \right) \begin{pmatrix} \underline{\underline{\Psi}}^0 \\ \underline{\underline{\Phi}}^0 \end{pmatrix} \right\} d\underline{\underline{\Omega}}_p \begin{pmatrix} \underline{\underline{F}}_i^T \\ \underline{\underline{H}}_i^T \end{pmatrix} \int_{\underline{\underline{\Omega}}_t} \begin{pmatrix} \underline{\underline{M}} \\ \underline{\underline{M}} \end{pmatrix} d\underline{\underline{\Omega}}_t \quad (62) $$

where

$$ \left\{ \begin{array}{l} \underline{\underline{\Psi}}^0 \equiv \underline{\underline{\Psi}}(\underline{\underline{p}}, t = 0) \\ \underline{\underline{\Phi}}^0 \equiv \underline{\underline{\Phi}}(\underline{\underline{p}}, t = 0) \end{array} \right. \quad (63) $$

## Appendix 3. Separated representation solver: approximation basis enrichment

From the alpha coefficients just computed the approximation basis can be enriched by adding the “best” couples $E_{n+1}(\underline{\underline{p}})F_{n+1}(t)$ and $G_{n+1}(\underline{\underline{p}})H_{n+1}(t)$. To determine the involved functions, we introduce the enriched approximation (30):
$$ \Xi = \begin{pmatrix} \psi(\underline{p}, t) \\ \phi(\underline{p}, t) \end{pmatrix} = \sum_{j=1}^{n} \begin{pmatrix} \alpha_j E_j(\underline{p}) F_j(t) \\ \beta_j G_j(\underline{p}) H_j(t) \end{pmatrix} + \begin{pmatrix} \underline{R}(\underline{p}) S(t) \\ \underline{V}(\underline{p}) W(t) \end{pmatrix} \qquad (64) $$

into the variational formulation (26).

As described in Section 3, the functions $E_{n+1}(\underline{p})$, $F_{n+1}(t)$, $G_{n+1}(\underline{p})$ and $H_{n+1}(t)$ are obtained by normalizing the functions $R(\underline{p})$, $S(t)$, $V(\underline{p})$ and $W(t)$.

In order to proceed with the discretization of the variational formulation, we must define the expression of the weighting function $\Xi^*$, whose separated representation writes:

$$ \Xi^* = \begin{pmatrix} R^*(\underline{p})S(t) + R(\underline{p})S^*(t) \\ V^*(\underline{p})W(t) + V(\underline{p})W^*(t) \end{pmatrix} \qquad (65) $$

Thus, introducing Eqs. (64) and (65) into the weak form (26) we obtain the non-linear discrete system:

$$ \underline{V}_1(\underline{R}, \underline{S}, \underline{V}, \underline{W}) + \underline{K}(\underline{R}, \underline{S}, \underline{V}, \underline{W}) \begin{bmatrix} \underline{R} \\ \underline{S} \\ \underline{V} \\ \underline{W} \end{bmatrix} = \underline{V}_2(\underline{R}, \underline{S}, \underline{V}, \underline{W}) \qquad (66) $$

where

$$ \left\{ \begin{aligned} \underline{K}(\underline{R}, \underline{S}, \underline{V}, \underline{W}) &= \underline{[K_1 \underline{K_2}]} \\ \underline{V}_1(\underline{R}, \underline{S}, \underline{V}, \underline{W}) &= \underline{[V_1^1 \underline{V_2^1}]} \\ \underline{V}_2(\underline{R}, \underline{S}, \underline{V}, \underline{W}) &= \underline{[V_1^2 \underline{V_2^2}]} \end{aligned} \right. \qquad (67) $$

with

$$ \underline{K}_1 = \left( \underline{\frac{N}{N_{22}}} \right) \left[ \left( \underline{\frac{S^T}{W^T}} \right) \left( \underline{\frac{M}{M}} \right) \left( \underline{\frac{S}{W}} \right) \right] + \left( \underline{\frac{N}{N}} \right) \left[ \left( \underline{\frac{S^T}{W^T}} \right) \left( \underline{\frac{M}{M}} \right) \left( \underline{\frac{S}{W}} \right) \right] \qquad (68) $$

$$ \underline{K}_2 = \left( \underline{\frac{M}{M}} \right) \left[ \left( \underline{\frac{R^T}{V^T}} \right) \left( \underline{\frac{N}{N}} \right) \left( \underline{\frac{R}{V}} \right) \right] + \left( \underline{\frac{M}{M}} \right) \left[ \left( \underline{\frac{R^T}{V^T}} \right) \left( \underline{\frac{N}{N_{22}}} \right) \left( \underline{\frac{R}{V}} \right) \right] \qquad (69) $$

$$ \begin{aligned} \underline{V}_1^1 = & \sum_{j=1}^{n} \begin{pmatrix} \alpha_j \\ \beta_j \end{pmatrix} \left[ \left( \underline{\frac{N}{N_{22}}} \right) \left( \underline{\frac{E_j}{G_j}} \right) \right] \left[ \left( \underline{\frac{S^T}{W^T}} \right) \left( \underline{\frac{M}{M}} \right) \left( \underline{\frac{F_j}{H_j}} \right) \right] \\ & + \sum_{j=1}^{n} \begin{pmatrix} \alpha_j \\ \beta_j \end{pmatrix} \left[ \left( \underline{\frac{N}{N}} \right) \left( \underline{\frac{E_j}{G_j}} \right) \right] \left[ \left( \underline{\frac{S^T}{W^T}} \right) \left( \underline{\frac{M}{M}} \right) \left( \underline{\frac{F_j}{H_j}} \right) \right] \\ & + \sum_{j=1}^{n} \begin{pmatrix} \beta_j \\ \alpha_j \end{pmatrix} \left[ \left( \underline{\frac{N}{N_{12}}} \right) \left( \underline{\frac{G_j}{E_j}} \right) \right] \left[ \left( \underline{\frac{S^T}{W^T}} \right) \left( \underline{\frac{M}{M}} \right) \left( \underline{\frac{H_j}{F_j}} \right) \right] \\ & + \left[ \left( \underline{\frac{N}{N_{12}}} \right) \left( \underline{\frac{V}{R}} \right) \right] \left[ \left( \underline{\frac{S^T}{W^T}} \right) \left( \underline{\frac{M}{M}} \right) \left( \underline{\frac{W}{S}} \right) \right] \end{aligned} \qquad (70) $$
$$ \underline{V}_2^1 = \sum_{j=1}^{n} \begin{pmatrix} \alpha_j \\ \beta_j \end{pmatrix} \left[ \begin{bmatrix} \underline{M} \\ \underline{M} \end{bmatrix} \begin{pmatrix} \underline{F}_j \\ \underline{H}_j \end{pmatrix} \right] \left[ \begin{bmatrix} \underline{R}^T \\ \underline{V}^T \end{bmatrix} \begin{pmatrix} \underline{N} \\ \underline{N} \end{pmatrix} \begin{pmatrix} \underline{E}_j \\ \underline{G}_j \end{pmatrix} \right] \right. \\ \left. + \sum_{j=1}^{n} \begin{pmatrix} \alpha_j \\ \beta_j \end{pmatrix} \left[ \begin{bmatrix} \underline{M} \\ \underline{M} \end{bmatrix} \begin{pmatrix} \underline{F}_j \\ \underline{H}_j \end{pmatrix} \right] \left[ \begin{bmatrix} \underline{R}^T \\ \underline{V}^T \end{bmatrix} \begin{pmatrix} \underline{N}_{11} \\ \underline{N}_{22} \end{pmatrix} \begin{pmatrix} \underline{E}_j \\ \underline{G}_j \end{pmatrix} \right] \right. \\ \left. + \sum_{j=1}^{n} \begin{pmatrix} \beta_j \\ \alpha_j \end{pmatrix} \left[ \begin{bmatrix} \underline{M} \\ \underline{M} \end{bmatrix} \begin{pmatrix} \underline{H}_j \\ \underline{F}_j \end{pmatrix} \right] \left[ \begin{bmatrix} \underline{R}^T \\ \underline{V}^T \end{bmatrix} \begin{pmatrix} \underline{N}_{21} \\ \underline{N}_{12} \end{pmatrix} \begin{pmatrix} \underline{G}_j \\ \underline{E}_j \end{pmatrix} \right] \right. \\ \left. + \left[ \begin{bmatrix} \underline{M} \\ \underline{M} \end{bmatrix} \begin{pmatrix} \underline{W} \\ \underline{S} \end{pmatrix} \right] \left[ \begin{bmatrix} \underline{R}^T \\ \underline{V}^T \end{bmatrix} \begin{pmatrix} \underline{N}_{21} \\ \underline{N}_{12} \end{pmatrix} \begin{pmatrix} \underline{V} \\ \underline{R} \end{pmatrix} \right] \right] \quad (71) $$

$$ \underline{V}_1^2 = - \left[ \int_{\Omega_p} \begin{pmatrix} \underline{N}_{11} \\ \underline{N}_{22} \end{pmatrix} \begin{pmatrix} \Psi(p, t=0) \\ \Phi(p, t=0) \end{pmatrix} d\Omega_p \right] \left[ \int_{\Omega_t} \begin{pmatrix} \underline{M}^T \\ \underline{M}^T \end{pmatrix} d\Omega_t \begin{pmatrix} \underline{S} \\ \underline{W} \end{pmatrix} \right] \right. \\ \left. - \left[ \int_{\Omega_p} \begin{pmatrix} \underline{N}_{21} \\ \underline{N}_{12} \end{pmatrix} \begin{pmatrix} \Phi(p, t=0) \\ \Psi(p, t=0) \end{pmatrix} d\Omega_p \right] \left[ \int_{\Omega_t} \begin{pmatrix} \underline{M}^T \\ \underline{M}^T \end{pmatrix} d\Omega_t \begin{pmatrix} \underline{S} \\ \underline{W} \end{pmatrix} \right] \right] \quad (72) $$

and

$$ \underline{V}_2^2 = - \left[ \begin{bmatrix} \underline{R}^T \\ \underline{V}^T \end{bmatrix} \int_{\Omega_p} \begin{pmatrix} \underline{N}_{11} \\ \underline{N}_{22} \end{pmatrix} \begin{pmatrix} \Psi(p, t=0) \\ \Phi(p, t=0) \end{pmatrix} d\Omega_p \right] \left[ \int_{\Omega_t} \begin{pmatrix} \underline{M}^T \\ \underline{M}^T \end{pmatrix} d\Omega_t \right] \right. \\ \left. - \left[ \begin{bmatrix} \underline{R} \\ \underline{V} \end{bmatrix} \int_{\Omega_p} \begin{pmatrix} \underline{N}_{21} \\ \underline{N}_{12} \end{pmatrix} \begin{pmatrix} \Phi(p, t=0) \\ \Psi(p, t=0) \end{pmatrix} d\Omega_p \right] \left[ \int_{\Omega_t} \begin{pmatrix} \underline{M}^T \\ \underline{M}^T \end{pmatrix} d\Omega_t \right] \right] \quad (73) $$

As discussed in Section 3 the solution of this non-linear system is accomplished within an alternating direction scheme. Thus, at each iteration *m* of this non-linear iterative solver, assuming the couple ($\underline{R}^m, \underline{S}^m$) to be known, the other couple ($\underline{V}^{m+1}, \underline{W}^{m+1}$) is computed by solving a system of size $N_p + N_t$. Now, with the just computed couple ($\underline{V}^{m+1}, \underline{W}^{m+1}$), the first one is updated by solving a system of size $N_p + N_t$ which leads to ($\underline{R}^{m+1}, \underline{S}^{m+1}$). The iteration procedure stops after reaching convergence, i.e. when $||\underline{R}^{m+1}, \underline{S}^{m+1} - \underline{R}^m, \underline{S}^m|| < \epsilon$ and $||\underline{V}^{m+1}, \underline{W}^{m+1} - \underline{V}^m, \underline{W}^m|| < \epsilon$, where $\epsilon$ is a sufficiently small parameter.

## References

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

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

[6] Y. Breaux, J.R. Clermont, Numerical simulation of non-Newtonian complex flows using the Stream-Tube method and memory integral constitutive equations, Int. J. Numer. Meth. Fluids 21 (1995) 371-389.

[7] H.J. Bungartz, M. Griebel, Sparse grids, Acta Numer. 13 (2004) 1-123.

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

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

[10] J.R. Clermont, Analysis of incompressible three-dimensional flows using the concept of stream tubes in relation with a transformation of the physical domain, Rheol. Acta 27 (1988) 357-362.

[11] D.R. Grevov, Formulations théoriques et simulations numériques d'écoulements de fluides viscoélastiques entre cylindres excentriques et en conduites, PhD Thesis, INPG-Grenoble, 2000.

[12] D.R. Grevov, J.R. Clermont, Numerical study of flows of complex fluids between eccentric cylinders using transformation functions, Int. J. Numer. Meth. Fluids 40 (2002) 669-695.
[13] D.R. Grecov, M. Normandin, J.R. Clermont, A numerical approach for computing flows by local transformations and domain decomposition using an optimization algorithm, Comput. Meth. Appl. Mech. Eng. 191 (2002) 4401-4419.

[14] D.R. Grecov, J.R. Clermont, Numerical simulations of non-stationary flows of non-Newtonian fluids between concentric and eccentric cylinders by Stream-Tube method and domain decomposition, Rheol. Acta 47 (5-6) (2008) 609-620.

[15] J.G. Hernandez Cifre, Th.M.A.O.M. Barenbrug, J.D. Schieber, B.H.A.A. van der Brule, Brownian dynamics simulation of reversible polymer networks under shear using a non-interacting dumbbell model, J. Non-Newtonian Fluid Mech. 113 (2003) 73-96.

[16] J.G. Hernandez Cifre, R. Pamies, A.-L. Kjoniksen, K.D. Knudsen, B. Nystrom, J. Garcia de la Torre, Brownian dynamics simulation of reversible polymer networks using a non-interacting bead-and-spring chain model, J. Non-Newtonian Fluid Mech. 146 (2007) 3-10.

[17] R. Keunings, Micro-macro methods for the multiscale simulation of 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.

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

[21] G. Marrucci, S. Bhargava, S.L. Copper, Models of shear-thickening behaviour in physically cross-linked networks, Macromolecules 26 (1993) 6483-6488.

[22] B. Mokdad, E. Pruliere, A. Ammar, F. Chinesta, On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach, Appl. Rheol. 17 (2007), 26494-1-26494-14.

[23] M. Normandin, J.R. Clermont, J. Guillet, C. Raveyre, Three-dimensional extrudate swell experimental and numerical study of a polyethylene melt obeying a memory-integral equation, J. Non-Newtonian Fluid Mech. 87 (1999) 1-25.

[24] M. Normandin, D.G. Radu, A. Mahmoud, J.R. Clermont, Finite element and Stream-Tube formulations: two-dimensional applications, Math. Comput. Simul. 60 (2002) 129-134.

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

[26] L. Pellens, K.H. Ahn, S.J. Lee, J. Mewiq, Evaluation of a transient network model for telechelic associative polymers, J. Non-Newtonian Fluid Mech. 121 (2004) 87-100.

[27] 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. Meth. Eng., State of the Art Rev. 16 (2009) 1-30.

[28] E. Rincon, A.E. Chavez, R. Herrera, O. Manero, Rheological modelling of complex fluids: a transient network model with microstates, J. Non-Newtonian Fluid Mech. 131 (2005) 303-334.

[29] A.S. Sarvestani, C.R. Picu, Network model for the viscoelastic behavior of polymer nanocomposites, Polymer 45 (2004) 7779-7790.

[30] F. Tanaka, S.F. Edwards, Viscoelastic properties of physically cross-linked networks: transient network theory, Macromolecules 25 (1992) 1516-1523.