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

HAL
open science

# Solving Parametric Complex Fluids Models in Rheometric Flows

Amine Ammar, Magdeleine Normandin, Francisco Chinesta

► To cite this version:

Amine Ammar, Magdeleine Normandin, Francisco Chinesta. Solving Parametric Complex Fluids Models in Rheometric Flows. Journal of Non-Newtonian Fluid Mechanics, 2010, 165 (23-24), pp.1588-1601. 10.1016/j.jnnfm.2010.08.006. hal-01007311

HAL Id: hal-01007311

https://hal.science/hal-01007311

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

Public Domain
# Solving parametric complex fluids models in rheometric flows

A. Ammarª, M. Normandinª, F. Chinestaᵇ

ª Laboratoire de Rhéologie, 1301 rue de la piscine, BP 53 Domaine universitaire, F-38041 Grenoble Cedex 9, France

ᵇ EADS Corporate Foundation International Chair, GEM: CNRS-Centrale Nantes, 1 rue de la Noe, BP 92101, F-44321 Nantes Cedex 3, France

Inverse identification of complex fluid behaviors is a tricky task because sometimes there are several rheological parameters and the identification procedure itself is quite expensive from the computational time viewpoint. Standard inverse identification procedures solve the model for a choice of the model parameters, and then parameters are updated trying to minimize the gap between the model predictions and some available experimental measures. Thus, the model has to be evaluated for each trial set of the model parameters. When models involve a great number of degrees of freedom the identification procedure becomes a computationally expensive task. In this paper we propose a new procedure able to solve once the model for any value of the model parameters. For this purpose, all the model parameters are considered as extra-coordinates of the model. Thus, the model results finally defined in a multidimensional space including the physical space x, the time t and a number of extra-coordinates related to the model parameters. The solution of such model needs for circumventing the curse of dimensionality illness that suffer multidimensional models.

* Complex fluids
* Inverse identification
* Parametric models
* Curse of dimensionality
* Proper Generalized Decomposition

## 1. Introduction

The main aim of this paper is the analysis of the possibility of solving models in which the parameters are considered as extra-coordinates. This procedure allows solving only once the partial differential equations defining the model to compute the solution for each value of the model parameters, but the price to be paid is the increase of the model dimensionality. The resulting model being defined in a highly dimensional space, its solution needs advanced techniques able to circumvent the redoubtable curse of dimensionality.

In this work we are firstly addressing the curse of dimensionality issue in the context of computational mechanics. Then, we illustrate the application of the Proper Generalized Decomposition (PGD) for solving such kind of multidimensional models suffering the so-called curse of dimensionality. Using the PGD the numerical complexity of the model solution scales linearly with the dimension of the space, instead of the exponential growing characteristic of mesh-based discretization techniques.

Finally, two simple rheological models will be considered and the parameters that they involve will be introduced into the models as additional extra-coordinates. The general solution will be then

computed by applying the PGD to the corresponding multidimensional models.

We would like to emphasize from the beginning of this work that the aim of this paper is only the introduction of a new methodology for solving multidimensional parametric models that could be applied in inverse identification of realistic models encountered in computational rheology. Thus, if we compute only once the solution of such multidimensional model we have access to the solution at each coordinate (space, time and model parameters). Now, as soon as the model parameters are fixed, the solution can be found by particularizing the general solution of such parameter choice. Then one could evaluate the cost function, update parameters and particularize again the same general solution for the new set of parameters. We need only one calculation instead the many calculations needed when using standard optimization or inverse identification strategies. Moreover, because the general solution includes the parameters as extra-coordinates, we could derive an exact expression for the sensibility of the solution with respect to such parameters. This is another advantage with respect to classical strategies where the calculation of sensibilities is a difficult task.

The introduction of parameters as extra-coordinates has enormous advantages with respect to the alternative of solving the model for different values of the parameters. We have just indicated that in the case of introducing the parameters as extra-coordinates we have direct access to the sensibilities needed in the vast majority of optimization strategies. On the other hand, when the number of parameters increases and they can take values in large intervals, the construction of the solution for each possible choice of param-
eters becomes simply prohibitive. Sometimes there are millions of possible choices. Finally, it is very important to notice, as indicated at the end of Section 1.1, that as the original model does not involve derivatives with respect to the parameters, the computation of the functions related to these parameters only implies algebraic calculations, and therefore the impact associated with the introduction of this kind of extra-coordinates on the computing time is negligible. Thus, if a technique for solving multidimensional models is available, the best alternative in general consists of introducing all the parameters as extra-coordinates.

In this paper we are focusing in the definition of parametric rheological models and in its solution by applying the PGD. Thus, we are not considering here neither the optimum searching procedure nor the treatment of advanced rheological models. The answer to the first point lies in procedures able to capture local and global minimum of given multidimensional functions. Concerning the second question, the procedure here introduced can be easily applied to any kind of model, linear or non-linear, defined in the frameworks of kinetic theory but also to models involving conformation tensors.

## 1.1. Efficient treatment of parametric models

Imagine for example that you are interested in solving the heat equation but that you do not know the material thermal conductivity, because it has a stochastic nature or simply because prior to solve the thermal model you should measure it. You have three possibilities: (i) you wait to know the conductivity before solving the heat equation (a conservative solution!); (ii) you solve the equation for many values of the conductivity (a sort of Monte Carlo) and then the work is done (the brute force approach!); or (iii) you solve the heat equation only once for any value of the conductivity (the cleverest alternative!). Obviously the third alternative is the most exciting one. To compute this "magic" solution it suffices to introduce the conductivity as an extra coordinate, playing the same role that the standard space and time coordinates, even if there are no derivatives concerning this extra-coordinate. This procedure runs, very well, and can be extended for introducing many other extra-coordinates: the source term, initial condition... This is the basic idea of spectral approaches in stochastic analysis, which are receiving a growing interest in computational science (see [21] and the references therein for a recent state of the art). It is easy to understand that after performing this type of calculations, a posteriori inverse identification or optimization can be easily handled. But the dream is not finished, many exciting models are waiting for new reformulations: in shape optimization the geometrical parameters could be introduced as new coordinates allowing to compute the solution for all possible geometries, and from that, identify the "optimum" (in many cases the global one instead of a local one reached by standard strategies).

Everything seems exciting, but the main question needs an answer: *How circumventing the curse of dimensionality that those enriched models suffer?*

Different techniques have been proposed for circumventing the curse of dimensionality, being Monte Carlo simulations the most widely used. Their main drawback is the statistical noise, when other than the moments of the distribution functions are computed. Other possibility lies in the use of sparse grids [9], within the deterministic framework, but they suffer also when the dimension of the space increases beyond a certain value (of about 20 dimensions as argued in [1]). In our knowledge there are few precedents of deterministic techniques able to circumvent efficiently the curse of dimensionality in highly multidimensional spaces. An appealing alternative consists of expressing the unknown field as a finite sum of functional products, i.e. expressing a generic multidimensional

function $\Psi(x_1, \cdots, x_N)$ as

$$ \Psi(x_1, \cdots, x_N) \approx \sum_{i=1}^{i=Q} \frac{1}{i!} F_i(x_1) \cdots \frac{N}{i!} F_i(x_N) \quad (1) $$

**Remark:** In this expression the coordinates $x_i$ denote any coordinate, scalar or vectorial, involving the physical space, the time or any other conformation coordinate (e.g. the conductivity in the example previously discussed).

Thus, if $M$ nodes are used to discretize each coordinate, the total number of unknowns involved in the solution is $Q \times N \times M$ instead of the $M^N$ degrees of freedom involved in mesh based discretizations. We must recall that these functions are not "a priori" known, they are computed by introducing the approximation separated representation (1) into the model weak form and then solving the resulting non-linear problem. The interested reader can refer to [12] for a detailed description of the numerical and algorithmic aspects. The construction of such approximation is called Proper Generalized Decomposition because this decomposition is not orthogonal but in many cases the number of terms in the finite sum is very close to the optimal decomposition obtained by applying the Proper Orthogonal Decomposition (POD) (or the Singular Value Decomposition, SVD) on the model solution. As it can be noticed in the expression of the approximation separated representation (1) the complexity scales linearly with the dimension of the space in which the model is defined, instead of the exponential growing characteristic of mesh based discretization strategies. In general, for many models, the number of terms $Q$ in the finite sum is quite reduced (few tens) and in all cases the approximation converges towards the solution associated with a fully tensorial product of approximation bases considered in each space $x_i$. Thus, we can conclude about the generality of the separated representation, but its optimality depends on the solution features. When the solution of physical models can be represented up to certain precision by a reduced number of functional products, the separated representation is an optimal representation, but if we consider a strictly non-separable function, the PGD based solver proceeds to enrich the approximation until including all the elements of the functional space, i.e. the $M^N$ functions involved in the full tensor product of the approximation bases considered in each space $x_i$.

When we consider multidimensional models involving the physical space, the time, and a number of "exotic" extra-coordinates, the verdict concerning the CPU time savings when using the PGD is implacable: the PGD allows solving models never until now solved, suffering the so-called curse of dimensionality, and that were qualified many times as irresolvable! The PGD allows solving them, in some minutes, using a simple laptop! Paradoxes of the history!

The application of the PGD to the solution of physically multi-dimensional models encountered in quantum chemistry and the kinetic theory description of complex fluids were deeply described in some of our former works [2,3,11], but also in [19]. It is important to recall that the use of separated representations for solving multidimensional models is not new, similar approximations were applied in the context of quantum chemistry within the so-called post-Hartree-Fock approaches (the interested reader can refer to the excellent Handbook on that topic edited by Le Bris [17], and in particular the general overview [10], and the references therein).

For stochastic analysis, the PGD has first been introduced by Nouy [21,22] (in this context, the PGD was initially called Generalized Spectral Decomposition). Finally, in the context of efficient non-linear solvers the space-time PGD (radial approximation) was successfully applied in [5] as well as within the LATIN framework in [15,16,20]. The multiscale approaches in space and time were addressed in [6,13], respectively, and in [14] we considered the separated representation of models defined in general domains
with non-homogeneous essential boundary conditions. The first mathematical issues were addressed in [7,18].

The present paper focuses in another issue, the consideration of extra-coordinates related to parametric deterministic models, and in particular models encountered in the description of non-Newtonian fluids. In the next section we illustrate the PGD procedure for solving a simple parametric thermal model: the parametric heat equation where the thermal conductivity, assumed unknown, is introduced as an additional extra-coordinates, like the space and the time. Then, this procedure will be applied to some simple models coming from computational rheology.

## 2. Illustrating the Proper Generalized Decomposition through an academic example

In what follows we are illustrating the construction of the Proper Generalized Decomposition by considering a quite simple problem, the parametric heat transfer equation:

$$ \frac{\partial u}{\partial t} - k \Delta u - f = 0 \qquad (2) $$

where $(\mathbf{x}, t, k) \in \Omega \times I \times \mathfrak{S}$ and for the sake of simplicity the source term is assumed constant, i.e. $f=cte$. Because the conductivity is considered unknown, it is assumed as a new coordinate defined in the interval $\mathfrak{S}$. Thus, instead of solving the thermal model for different values of the conductivity parameter we prefer introducing it as a new coordinate. The price to be paid is the increase of the model dimensionality; however, as the complexity of PGD scales linearly with the space dimension the consideration of the conductivity as a new coordinate allows faster and cheaper solutions.

The solution of Eq. (2) is searched under the form:

$$ u(\mathbf{x}, t, k) \approx \sum_{i=1}^{i=N} X_i(\mathbf{x}) \cdot T_i(t) \cdot K_i(k) \qquad (3) $$

In what follows we are assuming that the approximation at iteration $n$ is already done:

$$ u^n(\mathbf{x}, t, k) = \sum_{i=1}^{i=n} X_i(\mathbf{x}) \cdot T_i(t) \cdot K_i(k) \qquad (4) $$

and that at present iteration we look for the next functional product $X_{n+1}(\mathbf{x}) \cdot T_{n+1}(t) \cdot K_{n+1}(k)$ that for alleviating the notation will be denoted by $R(\mathbf{x}) \cdot S(t) \cdot W(k)$. Prior to solve the resulting nonlinear model related to the calculation of these three functions a model linearization is compulsory. The simplest choice consists in using an alternating directions fixed point algorithm. It proceeds by assuming $S(t)$ and $W(k)$ given at the previous iteration of the non-linear solver and then computing $R(\mathbf{x})$. From the just updated $R(\mathbf{x})$ and $W(k)$ we can update $S(t)$, and finally from the just computed $R(\mathbf{x})$ and $S(t)$ we compute $W(k)$. The procedure continues until reaching convergence. The converged functions $R(\mathbf{x})$, $S(t)$ and $W(k)$ allow defining the searched functions: $X_{n+1}(\mathbf{x}) = R(\mathbf{x})$, $T_{n+1}(t) = S(t)$ and $K_{n+1}(k) = W(k)$.

We are illustrating each one of the just referred steps:

I - Computing $R(\mathbf{x})$ from $S(t)$ and $W(k)$:

We consider the global weak form of Eq. (2):

$$ \int_{\Omega \times I \times \mathfrak{S}} u^* \left( \frac{\partial u}{\partial t} - k \Delta u - f \right) dx dtdk = 0 \qquad (5) $$

where the trial and test functions write, respectively:

$$ u(\mathbf{x}, t, k) = \sum_{i=1}^{i=n} X_i(\mathbf{x}) \cdot T_i(t) \cdot K_i(k) + R(\mathbf{x}) \cdot S(t) \cdot W(k) \qquad (6) $$

and

$$ u^*(\mathbf{x}, t, k) = R^*(\mathbf{x}) \cdot S(t) \cdot W(k) \qquad (7) $$

Introducing (6) and (7) into (5) it results

$$ \int_{\Omega \times I \times \mathfrak{S}} R^* \cdot S \cdot W \cdot \left( R \cdot \frac{\partial S}{\partial t} \cdot W - k \cdot \Delta R \cdot S \cdot W \right) dx dtdk \\ = - \int_{\Omega \times I \times \mathfrak{S}} R^* \cdot S \cdot W \cdot \left( \sum_{i=1}^{i=n} X_i \cdot \frac{\partial T_i}{\partial t} \cdot K_i - \sum_{i=1}^{i=n} k \cdot \Delta X_i \cdot T_i \cdot K_i - f \right) dx dtdk \qquad (8) $$

Now, being known all the functions involving the time and the parametric coordinate, we can integrate Eq. (8) in their respective domains $I \times \mathfrak{S}$. Integrating in $I \times \mathfrak{S}$ and taking into account the notation

$$ \left[ \begin{aligned} w_1 &= \int_{\mathfrak{S}} W^2 dk & s_1 &= \int_I S^2 dt & r_1 &= \int_{\Omega} R^2 d\mathbf{x} \\ w_2 &= \int_{\mathfrak{S}} kW^2 dk & s_2 &= \int_I S \cdot \frac{dS}{dt} dt & r_2 &= \int_{\Omega} R \cdot \Delta R d\mathbf{x} \\ w_3 &= \int_{\mathfrak{S}} W dk & s_3 &= \int_I S dt & r_3 &= \int_{\Omega} R d\mathbf{x} \\ w_4^i &= \int_{\mathfrak{S}} W \cdot K_i dk & s_4^i &= \int_I S \cdot \frac{dT_i}{dt} dt & r_4^i &= \int_{\Omega} R \cdot \Delta X_i d\mathbf{x} \\ w_5^i &= \int_{\mathfrak{S}} kW \cdot K_i dk & s_5^i &= \int_I S \cdot T_i dt & r_5^i &= \int_{\Omega} R \cdot X_i d\mathbf{x} \end{aligned} \right] \qquad (9) $$

Eq. (8) reduces to

$$ \int_{\Omega} R^* \cdot (w_1 \cdot s_2 \cdot R - w_2 \cdot s_1 \cdot \Delta R) d\mathbf{x} \\ = - \int_{\Omega} R^* \cdot \left( \sum_{i=1}^{i=n} w_4^i \cdot s_4^i \cdot X_i - \sum_{i=1}^{i=n} w_5^i \cdot s_5^i \cdot \Delta X_i - w_3 \cdot s_3 \cdot f \right) d\mathbf{x} \qquad (10) $$

Eq. (10) defines an elliptic steady state boundary value problem that can be solved by using any discretization technique operating on the model weak form (finite elements, finite volumes...). Another possibility consists in coming back to the strong form of Eq. (10):

$$ \begin{aligned} & w_1 \cdot s_2 \cdot R - w_2 \cdot s_1 \cdot \Delta R \\ = & - \left( \sum_{i=1}^{i=n} w_4^i \cdot s_4^i \cdot X_i - \sum_{i=1}^{i=n} w_5^i \cdot s_5^i \cdot \Delta X_i - w_3 \cdot s_3 \cdot f \right) \end{aligned} \qquad (11) $$

that could be solved by using any collocation technique (finite differences, SPH...).

II - Computing $S(t)$ from $R(\{\mathbf{x}\})$ and $W(k)$:

In the present case the test function writes:

$$ u^*(\mathbf{x}, t, k) = S^*(t) \cdot R(\mathbf{x}) \cdot W(k) \qquad (12) $$

Now, the weak form reads

$$ \int_{\Omega \times I \times \mathfrak{S}} S^* \cdot R \cdot W \cdot \left( R \cdot \frac{\partial S}{\partial t} \cdot W - k \cdot \Delta R \cdot S \cdot W \right) dx dtdk \\ = - \int_{\Omega \times I \times \mathfrak{S}} S^* \cdot R \cdot W \cdot \left( \sum_{i=1}^{i=n} X_i \cdot \frac{\partial T_i}{\partial t} \cdot K_i - \sum_{i=1}^{i=n} k \cdot \Delta X_i \cdot T_i \cdot K_i - f \right) dx dtdk \qquad (13) $$
that integrated in the domain $\Omega \times \mathfrak{S}$ and taking into account the notation (9) results:

$$ 
\begin{aligned}
& \int_I S^* \cdot \left( w_1 \cdot r_1 \cdot \frac{dS}{dt} - w_2 \cdot r_2 \cdot S \right) dt \\
& = - \int_I S^* \cdot \left( \sum_{i=1}^{i=n} w_4^i \cdot r_5^i \cdot \frac{dT_i}{dt} - \sum_{i=1}^{i=n} w_5^i \cdot r_4^i \cdot T_i - w_3 \cdot r_3 \cdot f \right) dt \quad (14)
\end{aligned}
$$

Eq. (14) represents the weak form of the ODE defining the time evolution of the field S that can be solved by using any stabilized discretization technique (SU, Discontinuous Galerkin, ...). The strong form of Eq. (14) reads:

$$ 
\begin{aligned}
& w_1 \cdot r_1 \cdot \frac{dS}{dt} - w_2 \cdot r_2 \cdot S \\
& = - \left( \sum_{i=1}^{i=n} w_4^i \cdot r_5^i \cdot \frac{dT_i}{dt} - \sum_{i=1}^{i=n} w_5^i \cdot r_4^i \cdot T_i - w_3 \cdot r_3 \cdot f \right) \quad (15)
\end{aligned}
$$

than can be solved by using backward finite differences, or higher order Runge-Kutta schemes, among many other possibilities.

### III - Computing $W(k)$ from $R(\mathbf{x})$ and $S(t)$:

In the present case the test function writes:

$$ u^*(\mathbf{x}, t, k) = W^*(k) \cdot R(\mathbf{x}) \cdot S(t) \quad (16) $$

Now, the weak form reads:

$$ 
\begin{aligned}
& \int_{\Omega \times I \times \mathfrak{S}} W^* \cdot R \cdot S \cdot \left( R \cdot \frac{\partial S}{\partial t} - W - k \cdot \Delta R \cdot S \cdot W \right) \mathrm{d}x \mathrm{d}t \mathrm{d}k \\
& = - \int_{\Omega \times I \times \mathfrak{S}} W^* \cdot R \cdot S \cdot \\
& \qquad \left( \sum_{i=1}^{i=n} X_i \cdot \frac{\partial T_i}{\partial t} \cdot K_i - \sum_{i=1}^{i=n} k \cdot \Delta X_i \cdot T_i \cdot K_i - f \right) \mathrm{d}x \mathrm{d}t \mathrm{d}k \quad (17)
\end{aligned}
$$

that integrated in the domain $\Omega \times I$ and taking into account the notation (9) results:

$$ 
\begin{aligned}
& \int_{\mathfrak{S}} W^* \cdot (r_1 \cdot s_2 \cdot W - r_2 \cdot s_1 \cdot k \cdot W) \mathrm{d}k \\
& = - \int_{\mathfrak{S}} W^* \cdot \left( \sum_{i=1}^{i=n} r_5^i \cdot s_4^i \cdot K_i - \sum_{i=1}^{i=n} r_4^i \cdot s_5^i \cdot k \cdot K_i - r_3 \cdot s_3 \cdot f \right) \mathrm{d}k \quad (18)
\end{aligned}
$$

Eq. (18) does not involve any differential operator. The strong form of Eq. (18) reads:

$$ (r_1 \cdot s_2 - r_2 \cdot s_1 \cdot k) \cdot W = - \left( \sum_{i=1}^{i=n} (r_5^i \cdot s_4^i - r_4^i \cdot s_5^i \cdot k) \cdot K_i - r_3 \cdot s_3 \cdot f \right) \quad (19) $$

that represents an algebraic equation. Thus, the introduction of parameters as additional model coordinates has not a noticeable effect in the computational cost, because the original equation does not contain derivatives with respect to those parameters.

There are other minimization strategies more robust and exhibiting faster convergence for building-up the PGD [12].

### 3. Tensor form

In this section we introduce the tensor form related to the procedure described in the previous section. This tensor formulation allows the discretization of any partial differential equation as we

described in [12]. Parametric non-Newtonian fluid models considered later in this paper make use of the notation introduced in this section.

For illustrating the procedure we are considering again the equation considered in the previous section:

$$ \frac{\partial u}{\partial t} - k \Delta u - f = 0 \quad (20) $$

whose weak form writes:

$$ \int_{\Omega \times I \times \mathfrak{S}} u^* \left( \frac{\partial u}{\partial t} - k \Delta u - f \right) \mathrm{d}x \mathrm{d}t \mathrm{d}k = 0 \quad (21) $$

is discretized using a standard finite element approximation in each model coordinate. For this purpose we define the vectors containing the nodal values of the different functions:

$$ u \approx \sum_{i=1}^{i=N} \mathbf{X}_i \otimes \mathbf{T}_i \otimes \mathbf{K}_i \quad (22) $$

and we introduce $\mathbf{N}$, $\mathbf{M}$ and $\mathbf{Q}$ that represent the vectors containing the shape functions associated, respectively with each space coordinate $\mathbf{x}$, $t$ and $k$.

The discrete weak form needs the introduction of the following matrices:

$$ 
\begin{aligned}
A_1^1 &= \int_{\Omega} \mathbf{N} \mathbf{N}^T d\Omega \\
A_2^1 &= \int_{I} \mathbf{M} \frac{d\mathbf{M}^T}{dt} dl \\
A_3^1 &= \int_{\mathfrak{S}} \mathbf{Q} \mathbf{Q}^T d\mathfrak{S} \\
A_1^2 &= \int_{\Omega} \nabla \mathbf{N} \cdot (\nabla \mathbf{N})^T d\Omega \\
A_2^2 &= \int_{I} \mathbf{M} \mathbf{M}^T dl \\
A_3^2 &= \int_{\mathfrak{S}} k \mathbf{Q} \mathbf{Q}^T d\mathfrak{S} \\
\mathbf{B}_1 &= \int_{\Omega} f \cdot \mathbf{N} \cdot d\Omega \\
\mathbf{B}_2 &= \int_{I} \mathbf{M} dl \\
\mathbf{B}_3 &= \int_{\mathfrak{S}} \mathbf{Q} d\mathfrak{S}
\end{aligned}
$$

Then, the discrete form results:

$$ A u = B \quad (23) $$

with

$$ A = \sum_{j=1}^{2} \frac{j}{A_1} \otimes \frac{j}{A_2} \otimes \frac{j}{A_3} \quad (24) $$

and

$$ B = \mathbf{B}_1 \otimes \mathbf{B}_2 \otimes \mathbf{B}_3 \quad (25) $$

Note that this equation is a tensorial equation. The unknown field $u$ contains the tensor product of the three functions given in 22. Thus, the resulting problem as in the previous section becomes non-linear and an appropriate linearization technique must be applied, as was the case in the description addressed in the previous section.
**Remarks:**

* No special continuity requirements are needed on the coordinate k. Thus, we could consider a piecewise constant approximation of functions defined in such coordinate that results in a unit matrix of $A_1^1$ being $A_2^2$ the diagonal matrix containing the values of this coordinate at the nodal positions.

* We can find in this new tensorial formulation all the operators introduced in the previous section. For example the term

$$\int_{\Omega} R^* \cdot (w_1 \cdot s_2 \cdot R - w_2 \cdot s_1 \cdot \Delta R) d\mathbf{x} \qquad (26)$$

writes in a discrete form as

$$\begin{aligned} & (\mathbf{R}^{*T} \mathbf{A}_1^1 \mathbf{R}) \cdot (\mathbf{S}^T \mathbf{A}_2^1 \mathbf{S}) \cdot (\mathbf{W}^T \mathbf{A}_3^1 \mathbf{W}) \\ & + (\mathbf{R}^{*T} \mathbf{A}_1^2 \mathbf{R}) \cdot (\mathbf{S}^T \mathbf{A}_2^2 \mathbf{S}) \cdot (\mathbf{W}^T \mathbf{A}_3^2 \mathbf{W}) \end{aligned} \qquad (27)$$

where **R**, **S**, **W** are the nodal descriptions of functions $R(x)$, $S(t)$ and $W(k)$, respectively.

* The solution allows computing, in the example here addressed, the temperature field at all positions, at any time and for any value of the thermal conductivity by solving once a higher dimensional model that includes the thermal conductivity as an extra-coordinate. In computational rheology one such approach could be applied by considering as extra coordinate the model parameters or the strain rate applied in simple rheometric flows.

## 4. Examples in computational rheology

In this section we are considering two simple rheological models defined in the kinetic theory framework. We recall that we are mainly interested in describing the treatment of parametric models, more than in addressing rheological issues. Thus, we are considering the model describing the flow induced orientation of fibers immersed in a viscous Newtonian fluid rheometric flow, assuming the solution dilute enough and focusing on the solution of the Fokker-Planck equation related to the fiber orientation distribution (thus the flow kinematics - fiber orientation evolution coupling is neglected). Then, we will focus on the Doi-Edwards model whose parametric form will be solved by introducing one or two extra-coordinates. First, both models will be summarized, but for more details the interested reader can refer to [8].

### 4.1. Parametric model of short fiber suspensions

Many processes involving complex fluids concerns the flow induced anisotropy. The flow kinematics is highly dependant on this induced microstructure as well as the final properties of the conformed parts.

Let **p** be the unit vector aligned in the fiber axis direction, with $\Psi(\mathbf{x}, \mathbf{p}, t)$ the orientation distribution function representing the fraction of fibers that at point $\mathbf{x}$, and time $t$ are aligned in the direction **p**. In the case of homogeneous flows (as is usually the case in rheometric devices) we can ignore the spatial dependence, writing $\Psi(\mathbf{p}, t)$, distribution that satisfies the normality condition

$$\int \Psi(\mathbf{p}, t) d\mathbf{p} = 1, \forall t \qquad (28)$$

We can define the second order orientation tensor as

$$\mathbf{a}_{ij}(t) = \int p_i \cdot p_j \Psi(\mathbf{p}, t) d\mathbf{p} \qquad (29)$$

The isotropic orientation state is defined by the uniform distribution function

$$\Psi(\mathbf{p}) = \frac{1}{4\pi} \qquad (30)$$

For the sake of simplicity we consider fibers with an infinite aspect ratio. In that case the orientation evolution given by the Jeffery's equation reduces to

$$\frac{d\mathbf{p}}{dt} = \nabla\mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla\mathbf{v} \cdot \mathbf{p})\mathbf{p} \qquad (31)$$

where $\nabla\mathbf{v}$ is the gradient of velocity tensor.

The evolution of the fiber orientation distribution $\Psi$ is given by the Fokker-Planck equation:

$$\frac{d\Psi(\mathbf{p}, t)}{dt} + \frac{\partial}{\partial \mathbf{p}} \{ \Psi(\mathbf{p}, t) \frac{d\mathbf{p}}{dt} \} - \frac{\partial}{\partial \mathbf{p}} \{ D_r \frac{\partial \Psi(\mathbf{p}, t)}{\partial \mathbf{p}} \} = 0 \qquad (32)$$

where the diffusion coefficient $D_r$ is the inverse of a characteristic relaxation time.

In this equation we are using the standard notation, but note that $(\partial/\partial\mathbf{p})$ denotes the derivative on the tangent plan to the unit sphere. In practice, we are discretizing this equation using Cartesian coordinates as described in [4].

Thus, if one considers steady state solutions and the dimensionless gradient of velocity obtained by dividing the velocity by that diffusion coefficient, the resulting Fokker-Planck equation can be rewritten as

$$\frac{\partial}{\partial \mathbf{p}} \{ \Psi(\mathbf{p}, t) \frac{d\mathbf{p}}{dt} \} - \frac{\partial}{\partial \mathbf{p}} \left\{ \frac{\partial \Psi(\mathbf{p}, t)}{\partial \mathbf{p}} \right\} = 0 \qquad (33)$$

where now the advective term $d\mathbf{p}/dt$ is computed from the dimensionless strain rate.

The Fokker-Planck equation can be rewritten in the compact form:

$$\mathcal{K}(\Psi(\mathbf{p}, t)) + \mathcal{L}(\Psi(\mathbf{p}, t)) = 0 \qquad (34)$$

expression that involves the following two operators:

$$\mathcal{K}(\Psi(\mathbf{p}, t)) = \frac{\partial}{\partial \mathbf{p}} (\Psi(\mathbf{p}, t) \cdot [\nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p})\mathbf{p}]) \qquad (35)$$

and

$$\mathcal{L}(\Psi(\mathbf{p}, t)) = - \frac{\partial}{\partial \mathbf{p}} \left( \frac{\partial \Psi(\mathbf{p}, t)}{\partial \mathbf{p}} \right) \qquad (36)$$

If we are interested in computing the steady solution of the model for any value of the dimensionless strain rate in a simple shear flow characterized by the dimensionless shear rate $G$ it suffices introducing this parameter as a model extra-coordinate and then look for a separated representation of the distribution function:

$$\Psi(\mathbf{p}, G) \approx \sum_{i=1}^{N} F_1^i(\mathbf{p}) \cdot F_2^i(G) \qquad (37)$$

Thus the model must be solved once instead of the solution of it for any value of $G$. Moreover, as the model does not contain derivatives with respect to this extra-coordinate, the impact of the introduction of $G$ as an extra-coordinate is quite reduced.

The introduction of this extra-coordinate implies the following form of the differential operators:

$$\begin{aligned} \mathcal{K}(\Psi(\mathbf{p}, G)) &= \frac{\partial}{\partial \mathbf{p}} (\Psi(\mathbf{p}, G) \cdot [\nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p})\mathbf{p}]) \\ &= G \frac{\partial}{\partial \mathbf{p}} \left( \Psi(\mathbf{p}, G) \cdot \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \cdot \mathbf{p} \right. \\ & \qquad \left. - \begin{bmatrix} \mathbf{p}^T \cdot \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \cdot \mathbf{p} \end{bmatrix} \mathbf{p} \right) \end{aligned} \qquad (38)$$
Then we can define the following discrete operators:

$$
\begin{align*}
\mathbf{A}_1^1 &= \int_{\Omega} \mathbf{N}\mathbf{K}|_{G=1}(\mathbf{N}^T) d\Omega \\
\mathbf{A}_2^1 &= \begin{pmatrix}
G_1 & 0 & \cdots & 0 \\
0 & G_2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & G_r
\end{pmatrix}
\end{align*}
$$

$$
\begin{gather*}
\mathbf{A}_1^2 = \int_{\Omega} \mathbf{N}\mathcal{L}(\mathbf{N}^T) d\Omega \\
\mathbf{A}_2^2 = \mathbf{I}_{(r,r)}
\end{gather*}
$$

where $\mathbf{I}_{(r,r)}$ denotes the unit matrix of size $r \times r$ (in the previous expressions $r$ denotes the number of nodes considered in the discrete representation of functions involving the extra-coordinate $G$) and $\mathbf{N}$ is the vector containing the shape function defined on the unit surface that corresponds to the conformation space $\Omega$ related to $\mathbf{p}$.

The integrals in $A_1^1$ and $A_2^1$ are computed after applying advective stabilization and integration by parts, respectively.

Finally, the discrete form of the operators defining the problem to solve is given by

$$
\mathcal{A} = \sum_{j=1}^{2} \mathbf{A}_{1}^{j} \otimes \mathbf{A}_{2}^{j} \qquad (39)
$$

$$
\mathcal{B} = 0 \tag{40}
$$

In order to enforce the normality condition we introduce a first set of functions verifying the normality constraint

$$
F_1^1(\mathbf{p}) = \frac{1}{4\pi}
$$

$$
\begin{gather*}
F_1^1(G) = 1 \\
\int_{\Omega} F_1^1(\mathbf{p}) d\Omega = 0, \forall i \ge 2 \tag{41}
\end{gather*}
$$

4.2. Doi–Edwards model

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

The key idea of these models is the use of the repetition mecha-
nism introduced by de Gennes to a tube (along which the molecule
can move) in order to describe the viscoelastic behavior of entan-
gled polymers. The molecule is described as sliding or repeating
through a tube whose contours are defined by the locus of entan-
glements with neighbouring molecules.

We are considering in this work as a second example and with-
out loss of generality the parametric Doi–Edwards model that is the
simplest model within the repetition picture. In addition to the rep-
etition mechanism, this model assumes an affine tube deformation
induced by the macroscopic flow, but neglects other phenomena
like the stretch of the chain and the Convective Constraint Release.
Within this repetition picture and these assumptions, the dynamics
of a single segment is given by

$$
\dot{\mathbf{u}} = \nabla \mathbf{v} \cdot \mathbf{u} - (\mathbf{u}^T \cdot \nabla \mathbf{v} \cdot \mathbf{u})\mathbf{u} \quad (42)
$$

where **u** is the unit vector describing the orientation of the tube seg-
ment and **v** is the velocity gradient. We also introduce the contour
length coordinate s, where s = 0 and s = 1 represent the chain ends.

As previously we assume homogeneous flows allowing to ignore
any space dependence. The distribution function Ψ is such that
Ψduds represents the joint probability that at time t a tube seg-
ment has an orientation in the interval [u, u + du] and contains the
chain segment labelled in the interval [s, s + ds]. Thus, the configu-
ration space results Ω × [0, 1], where Ω is the surface of the unit
sphere centered at the origin. The Fokker-Planck equation related
to the Doi-Edwards model is the convection-diffusion equation
that represents the distribution function balance:

$$
\frac{d\Psi}{dt} = -\frac{\partial}{\partial \mathbf{u}}(\Psi \dot{\mathbf{u}}) + \frac{\partial}{\partial s}\left(\frac{\partial \Psi}{\partial s}\right) \quad (43)
$$

that has been written as previously in a dimensionless form.

To solve the Fokker-Planck equation one needs to prescribe
appropriate boundary conditions at both tube ends s=0 and s=1 (the
orientation coordinate being defined on the unit sphere does not
require any boundary condition). Sometimes, as previously indi-
cated, an isotropic orientation distribution is assumed at both ends,
which reads:

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

Knowing the distribution function, the dimensionless stress can
be computed from:

$$
\tau_{ij} = \int_{0}^{1} \left( \int_{\Omega} \Psi(\mathbf{u}, s) (\mathbf{u}_i \cdot \mathbf{u}_j) d\Omega \right) ds \quad (45)
$$

If we focus in the steady state the Fokker-Planck equation writes

$$
\mathcal{K}(\Psi(\mathbf{u}, s)) + \mathcal{L}(\Psi(\mathbf{u}, s)) = 0 \qquad (46)
$$

where

$$
\mathcal{K}(\Psi(\mathbf{u}, s)) = \frac{\partial}{\partial \mathbf{u}} (\Psi(\mathbf{u}, s) \cdot [\nabla \mathbf{v} \cdot \mathbf{u} - (\mathbf{u}^T \cdot \nabla \mathbf{v} \cdot \mathbf{u})\mathbf{u}]) \quad (47)
$$

and

$$
\mathcal{L}(\Psi(\mathbf{u}, s)) = -\frac{\partial}{\partial s} \left( \frac{\partial \Psi(\mathbf{u}, s)}{\partial s} \right) \quad (48)
$$

In what follows we are considering the solution of the non-
parametric and parametric Doi–Edwards (DE) model.

4.2.1. Non-parametric DE model

The solution of this equation in Ω × [0, 1] is searched under the
form:

$$
\Psi(\mathbf{u}, s) \approx \sum_{i=1}^{N} F_{1}^{i}(\mathbf{u}) \cdot F_{2}^{i}(s) \qquad (49)
$$

Let **M** be the vector containing the shape functions in the space
[0,1]. The tensorial form of the problem writes:

$$
(A_1^1 \otimes A_2^1 + A_1^2 \otimes A_2^2) \left[ \sum_j \mathbf{F}_1^j \otimes \mathbf{F}_2^j \right] = 0 \quad (50)
$$

with

$$
\mathbf{A}_1^1 = \int_{\Omega} \mathbf{N} \mathbf{K} (\mathbf{N}^T) d\Omega
$$

$$
\mathbf{A}_2^1 = \int_{[0,1]} \mathbf{M}\mathbf{M}^T ds
$$
$$A_1^2 = \int_{\Omega} NN^T d\Omega$$

and

$$A_2^2 = \int_{[0,1]} ML(M^T)ds$$

In order to enforce the normality condition we consider a first couple of functions:

$$\begin{aligned} F_1^1(\mathbf{u}) &= \frac{1}{4\pi} \\ F_2^1(s) &= 1 \end{aligned} \qquad (51)$$

and we enforce at each iteration of the algorithm that

$$\int_{\Omega} F_1^i(\mathbf{u}) d\Omega = 0, \forall i \ge 2 \qquad (52)$$

Moreover, for satisfying the boundary conditions we enforce

$$F_2^i(0) = F_2^i(1) = 0, \forall i \ge 2 \qquad (53)$$

### 4.2.2. Parametric model involving a single parameter

In this section we are considering a simple shear flow and the dimensionless shear rate G as a new extra-coordinate of the model.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Evolution of a₁₂ as a function of the dimensionless shear rate.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Functions involved in the four first terms of the finite sums decomposition: $F_1^1(p), \dots, F_1^4(p)$ (top) and $F_2^1(G), \dots, F_2^4(G)$ (bottom).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Convergence and CPU time as a function of the number of terms involved in the finite sums decomposition. The error is defined as the residual L² norm.</figcaption></figure>

The separated representation of the solution writes:

$$ \Psi(\mathbf{u}, s, G) \approx \sum_{i=1}^{N} F_{1}^{i}(\mathbf{u}) \cdot F_{2}^{i}(s) \cdot F_{3}^{i}(G) \quad (54) $$

The first differential operator involved in the Fokker-Planck equation writes in the present case:

$$ \begin{aligned} \mathcal{K}(\Psi(\mathbf{u}, s, G)) &= \frac{\partial}{\partial \mathbf{u}} (\Psi(\mathbf{u}, s, G) \cdot [\nabla \mathbf{v} \cdot \mathbf{u} - (\mathbf{u}^T \cdot \nabla \mathbf{v} \cdot \mathbf{u})\mathbf{u}]) \\ &= G \frac{\partial}{\partial \mathbf{u}} \left( \Psi(\mathbf{u}, s, G) \cdot \left[ \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \cdot \mathbf{u} \right. \right. \\ & \quad \left. \left. - \left( \mathbf{u}^T \cdot \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \cdot \mathbf{u} \right) \mathbf{u} \right] \right) \quad (55) \end{aligned} $$

whose tensor form writes

$$ (A_1^1 \otimes A_2^1 \otimes A_3^1 + A_1^2 \otimes A_2^2 \otimes A_3^2) \left[ \sum_j \mathbf{F}_1^j \otimes \mathbf{F}_2^j \otimes \mathbf{F}_3^j \right] = 0 \quad (56) $$

with

$$ A_1^1 = \int_{\Omega} \mathbf{N} \cdot \mathcal{K}|_{G=1}(\mathbf{N}^T) d\Omega $$

$$ A_2^1 = \int_{[0,1]} \mathbf{M}\mathbf{M}^T ds $$

$$ A_3^1 = \begin{pmatrix} G_1 & 0 & \cdots & 0 \\ 0 & G_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & G_r \end{pmatrix} $$

$$ A_1^2 = \int_{\Omega} \mathbf{N}\mathbf{N}^T d\Omega $$

$$ A_2^2 = \int_{[0,1]} \mathbf{M}\mathcal{L}(\mathbf{M}^T)ds $$

$$ A_3^2 = \mathbf{I}_{(r,r)} $$

In order to enforce the normality condition the first term of the sum involve the functions:

$$ \begin{aligned} F_1^1(\mathbf{u}) &= \frac{1}{4\pi} \\ F_2^1(s) &= 1 \\ F_3^1(G) &= 1 \end{aligned} \quad (57) $$

<figure><img src="image_5.png" /><figcaption>Fig. 4. Distribution function $\Psi(\mathbf{u}, s=0.5)$ (left) and $\Psi(\mathbf{u}_i, s)$, $i=1, 2, \dots$ (right).</figcaption></figure>
and we enforce at the subsequent enrichments that

$$ \int_{\Omega} F_{1}^{i}(\mathbf{u}) d\Omega = 0, \forall i \geq 2 \qquad (58) $$

For ensuring the verification of the boundary conditions we enforce

$$ F_{2}^{i}(0) = F_{2}^{i}(1) = 0, \forall i \geq 2 \qquad (59) $$

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

<figure><img src="image_7.png" /><figcaption>Fig. 5. Functions involved in the first four terms of the finite sums decomposition of $\Psi(\mathbf{u}, s)$.</figcaption></figure>

### 4.2.3. Parametric model involving two parameters

In this section, and for considering different flow conditions, we assume the dimensionless velocity gradient given by

$$ \nabla \mathbf{v} = \begin{pmatrix} -E/2 & G & 0 \\ 0 & E & 0 \\ 0 & 0 & -E/2 \end{pmatrix} \qquad (60) $$

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 6. Evolution of the error - residual L² norm - (on the left) and the computing time in seconds (on the right) as a function of the terms involved in the finite sums decomposition.</figcaption></figure>

considering both G and the ratio e = E/G as two new model extra-coordinates. The r' nodes allowing the discrete description of this last extra-coordinate (e) are given by e<sub>r'</sub>.

Now, the differential operator associated with the orientation
writes:

$$
A_2^2 = \int_{[0,1]} \mathbf{MM}^T ds
$$

$$
\begin{equation}
\begin{aligned}
\kappa(\Psi(\mathbf{u}, s, G, e)) ={}& G \cdot \underbrace{\frac{\partial}{\partial \mathbf{u}} \left( \Psi(\mathbf{u}, s, G, e) \cdot \left[ \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \cdot \mathbf{u} - \left( \mathbf{u}^T \cdot \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \cdot \mathbf{u} \right) \mathbf{u} \right] \right)}_{\kappa_1} \\
& + G \cdot e \cdot \underbrace{\frac{\partial}{\partial \mathbf{u}} \left( \Psi \cdot \left[ \begin{pmatrix} -1/2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -1/2 \end{pmatrix} \cdot \mathbf{u} - \left( \mathbf{u}^T \cdot \begin{pmatrix} -1/2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -1/2 \end{pmatrix} \cdot \mathbf{u} \right) \mathbf{u} \right] \right)}_{\kappa_2}
\end{aligned}
\tag{61}
\end{equation}
$$

In the present case it is easy to verify that the tensor form writes:

$$
(A_1^1 \otimes A_2^1 \otimes A_3^1 \otimes A_4^1 + A_1^2 \otimes A_2^2 \otimes A_3^2 \otimes A_4^2 + A_1^3 \otimes A_2^3 \otimes A_3^3 \otimes A_4^3) \times \\
\left[ \sum_j \mathbf{F}_1^j \otimes \mathbf{F}_2^j \otimes \mathbf{F}_3^j \otimes \mathbf{F}_4^j \right] = 0 \qquad (62)
$$

with

$$
A_1^1 = \int_{\Omega} \mathbf{N} \kappa_1 (\mathbf{N}^T) d\Omega
$$

$$
A_2^1 = \int_{[0,1]} \mathbf{MM}^T ds
$$

$$
A_3^1 = \begin{pmatrix} G_1 & 0 & \cdots & 0 \\ 0 & G_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & G_r \end{pmatrix}
$$

$$
A_4^1 = \mathbf{I}_{(r', r')}
$$

$$
A_1^2 = \int_{\Omega} \mathbf{N} \kappa_2 (\mathbf{N}^T) d\Omega
$$

$$
A_3^2 = \begin{pmatrix} G_1 & 0 & \cdots & 0 \\ 0 & G_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & G_r \end{pmatrix}
$$

<figure><img src="image_10.png" /><figcaption>Fig. 7. Evolution of the shear stress as a function of the dimensionless shear rate G.</figcaption></figure>
$$ \mathcal{A}_4^2 = \begin{pmatrix} e_1 & 0 & \cdots & 0 \\ 0 & e_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & e_{r'} \end{pmatrix} $$

$$ \mathcal{A}_1^3 = \int_{\Omega} \mathbf{N} \mathbf{N}^T d\Omega $$

$$ \mathcal{A}_2^3 = \int_{[0,1]} \mathbf{M} \mathcal{L}(\mathbf{M}^T) \mathrm{d}\mathbf{s} $$

$$ \mathcal{A}_3^3 = \mathbf{I}_{(r,r)} $$

$$ \mathcal{A}_4^3 = \mathbf{I}_{(r',r')} $$

<figure><img src="image_11.png" /><figcaption>Fig. 8. Functions involved in the first four sums of the separated representation.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 9. Evolution of the error – residual L² norm – (on the left) and the computing time in seconds (on the right) as a function of the terms involved in the finite sums decomposition of the single parameter DE model.</figcaption></figure>

As in the previous case the model is solved under the normality conditions enforced as described before.

## 5. Numerical results

In this section we are illustrating the capabilities of the parametric solution procedure through some quite simple academic examples involving short fiber suspensions and the Doi–Edwards parametric models.

All the models here considered were solved by using the Proper Generalized Decomposition, that is, a separated representation of the unknown field making use of functions defined in the configurational spaces and others defined in the parametric spaces. The different solutions where compared with the ones obtained by fixing some particular values of the parameters and solving the resulting low dimensional models by using the FEM. In all the cases the deviations were lower than 0.1 per cent for the converged solutions.

### 5.1. Short fiber suspension parametric model

In this section we consider the kinetic theory model related to suspensions of short fibers deeply described in section 4.1 in a simple shear flow, where the dimensionless shear rate G is assumed as an extra-coordinate.

In this case, the solution of Fokker-Planck equation allows computing the distribution function Ψ(p, G) given the fraction of fibers that are aligned in a direction given by the unit vector p when a dimensionless shear rate G is applied. Thus, instead of solving a problem defined on the unit sphere for each value of G, we should solve a single problem defined in a higher dimensional space that now includes the unit sphere as well as the dimensionless shear rate G. Because of the application of the Proper Generalized Decomposition (PGD) the increase of the model dimensionality has a negligible impact on the computational efforts.

Fig. 1 depicts the evolution of the components $a_{12}$ of the second order orientation tensor as a function of the dimensionless shear rate G. The components of the second order orientation tensor are defined as

$$ \mathbf{a}_{ij}(G) = \int_{\Omega} (\mathbf{p}_i \cdot \mathbf{p}_j) \Psi(\mathbf{p}, G) d\Omega \quad (63) $$

where Ω represents the surface of the unit sphere.

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

Fig. 2 depicts the four most relevant functions depending on the orientation $F_1^1(p), \dots, F_1^4(p)$ as well as the four first functions depending on the parametric coordinate $F_2^1(G), \dots, F_2^4(G)$.

Finally, Fig. 3 shows the evolution of the error (on the left) and the computing time in seconds (on the right) as a function of the number of terms involved in the finite sums decomposition.

### 5.2. Doi–Edwards model

In this section we are considering the Doi–Edwards model that will be analyzed in three scenarios: (i) the standard non-parametric model; (ii) the model defined in a simple shear flow that includes the dimensionless shear rate as a model extra-coordinate; and (iii) the model defined in a prescribed flow exhibiting shear and extension rates, both introduced as extra-coordinates in the DE model.

#### 5.2.1. Non-parametric Doi–Edwards model

In the first scenario, the one concerning the non-parametric DE model, the kinetic theory model involves the distribution function Ψ(u, s) defined in Ω × [0, 1]. The solution of that model in a prescribed simple shear flow characterized by G = 100 allows to compute the orientation distribution at the middle of the chain s = 0.5, that is, Ψ(u, s = 0.5), depicted in Fig. 4(left). Moreover, we could represent for different given orientations uᵢ, i = 1, 2, ..., the

<figure><img src="image_14.png" /><figcaption>Fig. 10. First normal stress difference (τ₁₁ – τ₂₂) as a function of G and e.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 11. Functions involved in the first four sums of the separated representation.</figcaption></figure>
evolution of the distribution function along the chain, that is, $\Psi(\mathbf{u}_i, s)$, $i=1, 2, \dots$. The resulting curves are represented in the same figure, at the right.

Fig. 5 depicts the functions involved in the four first terms of the finite sums decomposition, that is, functions $F_1^1(\mathbf{u})$, $\dots, F_1^4(\mathbf{u})$ (top) and $F_2^1(s)$, $\dots, F_2^4(s)$ (bottom).

Finally, Fig. 6 shows the evolution of the error (on the left) and the computing time in seconds (on the right) as a function of the terms involved in the finite sums decomposition.

### 5.2.2. Parametric Doi–Edwards model involving a single parameter

In this second scenario, the kinetic theory model involves the distribution function $\Psi(\mathbf{u}, s, G)$ defined in $\Omega \times [0, 1] \times [G_{min}, G_{max}]$. The solution of that model in a prescribed simple shear flow characterized by the dimensionless shear rate $G \in [G_{min}, G_{max}]$ allows computing the shear stress as a function of G as depicted in Fig. 7.

Fig. 8 depicts the functions involved in the four first terms of the finite sums decomposition:

$$ \Psi(\mathbf{u}, s, G) \approx \sum_{i=1}^{N} F_{1}^{i}(\mathbf{u}) \cdot F_{2}^{i}(s) \cdot F_{3}^{i}(G) $$

that is, functions $F_1^1(\mathbf{u})$, $\dots, F_1^4(\mathbf{u})$ (top), $F_2^1(s)$, $\dots, F_2^4(s)$ (bottom-left) and $F_3^1(G)$, $\dots, F_3^4(G)$ (bottom-right).

Finally, Fig. 9 shows the evolution of the error (on the left) and the computing time in seconds (on the right) as a function of the terms involved in the finite sums decomposition. As previously, we can notice that the first approximation basis enrichments lead to the higher error reduction. This behavior was deeply discussed in some of our former works, in which we proposed that after some basis enrichments, when the slope of the error curves decreases, one could proceed to refine the different meshes, allowing to come back to the highest slope of the error evolution. In the simulation performed here, we don't include this strategy and the meshes remain unchanged during the whole simulation. Thus, Fig. 9 represents the most unfavorable scenario.

### 5.2.3. Parametric Doi–Edwards containing two parameters

Finally, the kinetic theory model involves the distribution function $\Psi(\mathbf{u}, s, G, e)$ defined in $\Omega \times [0, 1] \times [G_{min}, G_{max}] \times [e_{min}, e_{max}]$. The solution of that model in a prescribed flow combining shear and elongation, both characterized by the dimensionless shear and the elongation rates, $G \in [G_{min}, G_{max}]$ and $e \in [e_{min}, e_{max}]$, respectively, allows computing for example the first normal stress difference for any value of the parameters G and e as depicted in Fig. 10.

Finally, Fig. 11 depicts the functions involved in the four first terms of the finite sums decomposition:

$$ \Psi(\mathbf{u}, s, G, e) \approx \sum_{i=1}^{N} F_{1}^{i}(\mathbf{u}) \cdot F_{2}^{i}(s) \cdot F_{3}^{i}(G) \cdot F_{4}^{i}(e) $$

that is, functions $F_1^1(\mathbf{u})$, $\dots, F_1^4(\mathbf{u})$ (top), $F_2^1(s)$, $\dots, F_2^4(s)$ (middle-left), $F_3^1(G)$, $\dots, F_3^4(G)$ (middle-right) and $F_4^1(e)$, $\dots, F_4^4(e)$ (bottom).

## 6. Conclusions

In this paper we have addressed the treatment of parametric models encountered in computational rheology. Those models involve sometimes some parameters that are a priori unknown and that should be identified by using an appropriate inverse technique. However, in general such inverse strategies runs by assuming a tentative set of parameters and comparing the numerical predictions with the experimental data. Thus, a direct problem must be

solved for each trail set of parameters, with the associate impact on the computing time. In this paper we presented an alternative methodology consisting in the introduction of all the model parameters as extra-coordinates of the model. The price to be paid for avoiding the numerous solutions that standard strategies imply, is the solution of a single problem but now defined in a higher multidimensional space. Mesh based discretization techniques fail when the dimension of the space increases, however, the use of a separated representation allows for computing efficiently models defined in highly multidimensional spaces involving hundreds coordinates as we proved in some of our former works [3]. Thus, the separated representation of parametric models could constitute an appealing alternative to standard methodologies in inverse rheology. The examples analyzed in this works confirm it, even if some deeper analysis should be performed.

## References

[1] Y. Achdou, O. Pironneau, Computational Methods for Option Pricing, Siam Froniers in Applied Mathematics, 2005.

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

[4] A. Ammar, E. Pruliere, F. Chinesta, M. Laso, Reduced numerical modeling of flows involving liquid-crystalline polymers, J. Non-Newtonian Fluid Mech. 160 (2009) 140–156.

[5] A. Ammar, M. Normandin, F. Daim, D. Gonzalez, E. Cueto, F. Chinesta, Non-incremental strategies based on separated representations: applications in computational rheology. Commun. Math. Sci. 8 (3) (2010) 671–695.

[6] A. Ammar, F. Chinesta, E. Cueto, Coupling finite elements and separated representations. Int. J. Multiscale Comput. Eng., in press.

[7] A. Ammar, F. Chinesta, A. Falco, On the convergence of a greedy rank-one update algorithm for a class of linear systems. Arch. Comput. Methods Eng., in press.

[8] R.B. Bird, C.F. Crutiss, R.C. Armstrong, O. Hassager, Dynamic of Polymeric Liquid, Volume 2: Kinetic Theory, John Wiley and Sons, 1987.

[9] H.J. Bungartz, M. Griebel, Sparse grids, Acta Numerica 13 (2004) 1–123.

[10] E. Cancés, M. Defranceschi, W. Kutzelnigg, C. Le Bris, Y. Maday, Computational Quantum Chemistry: A Primer, Handbook of Numerical Analysis, Vol. X, Elsevier, 2003, 3–270.

[11] F. Chinesta, A. Ammar, P. Joyot, The nanometric and micrometric scales of the structure and mechanics of materials revisited: an introduction to the challenges of fully deterministic numerical descriptions, Int. J. Multiscale Comput. Eng. 6/3 (2008) 191–213.

[12] F. Chinesta, A. Ammar, E. Cueto, Recent advances in the use of the proper generalized decomposition for solving multidimensional models. Arch. Comput. Methods Eng., in press.

[13] F. Chinesta, A. Ammar, E. Cueto, Proper generalized decomposition of multiscale models. Int. J. Num. Methods Eng., in press.

[14] D. Gonzalez, A. Ammar, F. Chinesta, E. Cueto, Recent advances in the use of separated representations. Int. J. Num Methods Eng. 81 (5) (2010) 637–659.

[15] P. Ladeveze, Nonlinear Computational Structural Mechanics, Springer, NY, 1999.

[16] P. Ladeveze, J.-C. Passieux, D. Neron, The LATIN multiscale computational method and the proper orthogonal decomposition, Comput. Methods Appl. Mech. Eng. 199 (21–22) (2010) 1287–1296.

[17] C. Le, Bris (Eds.), Handbook of Numerical Analysis, Vol. X: Computational Chemistry, Elsevier, 2003.

[18] C. Le Bris, T. Lelièvre, Y. Maday, Results and questions on a nonlinear approximation approach for solving high-dimensional partial differential equations. Constructive Approximation, DOI 10.1007/s00365-009-9071-1.

[19] G.M. Leonenko, T.N. Phillips, On the solution of the Fokker–Planck equation using a high-order reduced basis approximation. Comput. Methods Appl. Mech. Eng., in press.

[20] A. Nouy, P. Ladeveze, Multiscale computational strategy with time and space homogenization: a radial-type approximation technique for solving micro-problems, Int. J. Multiscale Comput. Eng. 170/2 (2004).

[21] A. Nouy, Recent developments in spectral stochastic methods for the solution of stochastic partial differential equations, Arch. Comput. Methods Eng. 16/3 (2009) 251–285.

[22] A. Nouy, Proper generalized decompositions and separated representations for the numerical solution of high-dimensional stochastic problems. Arch. Comput. Methods Eng. 199 (1–4) (2009) 158–168.