# Non Incremental Transient Solution of the Rayleigh-Bénard Convection Model Using the PGD

Mohammad Saeid Aghighi, Amine Ammar, Christel Metivier, Magdeleine Normandin, Francisco Chinesta

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

## To cite this version:

Mohammad Saeid Aghighi, Amine Ammar, Christel Metivier, Magdeleine Normandin, Francisco Chinesta. Non Incremental Transient Solution of the Rayleigh-Bénard Convection Model Using the PGD. Journal of Non-Newtonian Fluid Mechanics, 2013, 200, pp.65-78. 10.1016/j.jnnfm.2012.11.007. hal-01007106

HAL Id: hal-01007106

https://hal.science/hal-01007106

Submitted on 6 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Non-incremental transient solution of the Rayleigh-Bénard convection model by using the PGD

M.S. Aghighi<sup>a</sup>, A. Ammar<sup>a</sup>, C. Metivier<sup>b</sup>, M. Normandin<sup>b</sup>, F. Chinesta<sup>c,*</sup>

This paper focuses on the non-incremental solution of transient coupled non-linear models, in particular the one related to the Rayleigh-Bénard flow problem that models natural thermal convection. For this purpose we are applying the so-called Proper Generalized Decomposition that proceeds by performing space-time separated representations of the different unknown fields involved by the flow model. This non-incremental solution strategy allows significant computational time savings and opens new perspectives for introducing some flow and/or fluid parameters as extra-coordinates.

**Keywords:** Rayleigh-Bénard model, Power-law fluids, Model reduction, Non-incremental solution, Proper Generalized Decomposition PGD

## 1. Introduction

This paper focuses on the non-incremental solution of transient coupled non-linear models, in particular the one related to the Rayleigh-Bénard flow problem that models natural thermal convection. Several systems and industrial processes are based on natural convection, justifying the impressive volume of work devoted to its understanding and efficient solution during more than one century. This model, quite simple in appearance, deserves many surprises related to its intricate nature and many issues concerning its numerical solution, mainly in the case of non-Newtonian fluids and/or when the Rayleigh number is large enough to induce the transition to the turbulence.

First studies on this problem were motivated by the Bénard's experiments around 1900 [7] who considered the stability of a fluid layer heated from its basis. A linear stability analysis was proposed in 1916 by Lord Rayleigh [25] underlying the buoyancy-driven source of instability. The first chapters of Chandrasekhar's book [9] present the linear theory within the Boussinesq approximation. Non-linear approaches were reviewed in [8,19].

In the case of Newtonian fluids, linear and non-linear stability analyses of the two-dimensional Rayleigh-Bénard model in arbitrary finite domains were performed by Park and Heo [23]. Vedantam and co-workers [26] performed computational fluid dynamics simulations of the Rayleigh-Bénard flow model for low Prandtl

numbers by using the FLUENT software. In 2008, steady state two-dimensional solutions of the Rayleigh-Bénard problem were obtained numerically by Ouertatani et al. [21] by using finite volumes discretizations. For non-Newtonian fluids Park and Park [22] considered the linear hydrodynamic stability problem of viscoelastic fluids in arbitrary domains and the effects of yield stress on the Rayleigh-Bénard instability was analyzed in [28]. Vikhansky [27] considered the effect of yield stress on the Rayleigh-Bénard model concerning a visco-plastic fluid.

The main contribution here addressed concerns the solution of the non-linear Navier-Stokes equation with a temperature dependent density describing the thermally induced flow. On the other hand, the temperature field is governed by the advection-diffusion heat equation, coupled to the momentum and mass balances. Thus the resulting time-dependent model becomes non-linear (because the inertia term and the eventual non-linear constitutive equations) and strongly coupled. When using standard discretizations one must be careful with respect to many numerical choices concerning the simultaneous solution (monolithic) of all the balance equations, the use of accurate discretizations, adaptive time steps, robust stabilizations of both the advective terms and the mixed formulation, and an adequate treatment of non-linearities.

Many time-dependent models involve a large spectrum of characteristic times that makes difficult their solution by considering basic incremental time discretization techniques. In such cases, the time step is extremely small as a consequence of numerical stability requirements. Thus, simulations over the much larger time interval of interest, which typically requires at least the solution of a large linear algebraic system at each time step, becomes too

<footnote><sup>a</sup>Arts et Métiers ParisTech, 2 Boulevard du Ronceray, BP 93525, F-49035 Angers cedex 01, France</footnote>

<footnote><sup>b</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>c</sup>EADS Corporate Foundation International Chair, GEM, UMR CNRS-Centrale Nantes, 1 rue de la Noe, BP 92101, F-44321 Nantes cedex 3, France</footnote>
expensive in many cases, mainly when the model must be solved many times because we are concerned by optimization or inverse identification issues. Moreover, models coming from the physics of materials and processes are in general non-linear and strongly coupled.

It was in this scenario that Pierre Ladeveze proposed in the 1980s a new powerful simulation paradigm, the LATIN method [15] that combines two key ingredients: (i) an efficient non-linear treatment and (ii) a space-time separated representation. The Ladeveze's group accomplished remarkable progresses in the solution of non-linear models within a multi-scale and multi-physics framework during the last decades [20,16]. An exhaustive review can be found in [13] and the references therein.

In [1] we generalized the Ladeveze's space-time separated representation for addressing models involving many coordinates, as the one encountered in the kinetic theory descriptions of materials and processes that involve many configurational or conformational coordinates describing the rich microstructures. In [2] the time was also included as a new coordinate in the separated representation allowing for non-incremental simulations of high-dimensional models and then in [18] general separated representations were applied for solving non-linear high-dimensional models, where the non-linearities were treated by applying quite standard techniques (e.g. Newton, fixed point ...). Techniques using separated representations, where the functions involved in such approximation are calculated on-the-fly were called Proper Generalized Decompositions - PGD -. Coupled and multi-scale models, where addressed in [20,16,11,6]. Some recent works on the application of the PGD in computational rheology can be found in [24,3]. The interested reader can also refer to [11-13] for some recent reviews on the PGD methodology.

In our knowledge space-time separated representations for computing the non-incremental solution of transient non-linear and coupled flow models have never been explored. This work is a first step in this direction, and for this reason we restrict our analysis to quite low Rayleigh numbers, in order to ensure that the resulting thermal induced flow remains laminar. More complex scenarios will be considered in future works, by introducing visco-elastic constitutive equations, material or process parameters as extra-coordinates as illustrated in [5] or by increasing the Rayleigh number to move beyond the laminar-turbulent transition.

In what follows we start by introducing in Section 2 the main ideas of the PGD. In Section 3 we summarize the space-time separated representation constructor applied in Section 4 for discretizing the Rayleigh-Bénard flow model. Section 5 presents some numerical results concerning both Newtonian and power-law fluids.

## 2. The Proper Generalized Decomposition at a glance

Consider a problem defined in a space of dimension $d$ for the unknown field $u(x_1, ..., x_d)$. Here, the coordinates $x_i$ denote any usual coordinate (scalar or vectorial) related to physical space, time, or conformation space, for example, but they could also include problem parameters such as boundary conditions or material parameters. We seek a solution for $(x_1, ..., x_d) \in \Omega_1 \times \dots \times \Omega_d$.

The PGD yields an approximate solution in the separated form:

$$u(x_1, \dots, x_d) \approx \sum_{i=1}^{N} F_i^1(x_1) \times \dots \times F_i^d(x_d). \quad (1)$$

The PGD approximation is thus a sum of $N$ functional products involving each a number $d$ of functions $F_i^j(x_j)$ that are unknown a priori. It is constructed by successive enrichment, whereby each functional product is determined in sequence. At a particular enrichment step $n+1$, the functions $F_i^j(x_j)$ are known for $i \le n$ from the previous steps, and one must compute the new product involving

the $d$ unknown functions $F_{n+1}^j(x_j), j = 1, \dots, d$. This is achieved by invoking the weak form of the problem under consideration. The resulting discrete system is non-linear, which implies that iterations are needed at each enrichment step. A low-dimensional problem can thus be defined in $\Omega_j$ for each of the $d$ functions $F_{n+1}^j(x_j), j = 1, \dots, d$.

If $M$ nodes are used to discretize each coordinate, the total number of PGD unknowns is $N \times M \times d$ instead of the $M^d$ degrees of freedom involved in standard mesh-based discretizations. Moreover, all numerical experiments carried out to date with the PGD show that the number of terms $N$ required to obtain an accurate solution is not a function of the problem dimension $d$, but it rather depends on the regularity of the exact solution. The PGD thus avoids the exponential complexity with respect to the problem dimension.

In many applications studied to date, $N$ is found to be as small as a few tens, and in all cases the approximation converges towards the solution associated with the complete tensor product of the approximation bases considered in each $\Omega_j$. Thus, we can be confident about the generality of the separated representation (1), but its optimality depends on the solution regularity, the specificities of the differential operator involved and the separated representation constructor.

## 3. Non-incremental solutions of transient models within the PGD framework

In this section we are illustrating the discretization of time dependent partial differential equations using a space-time separated representation (radial approximation in the Ladeveze's terminology) of the unknown field.

Let us consider the advection-diffusion equation

$$\frac{\partial u}{\partial t} - a \cdot \Delta u + \mathbf{v} \cdot \nabla u = f(\mathbf{x}, t) \quad \text{in } \Omega \times (0, t_{\text{max}}] \quad (2)$$

with the following initial and boundary conditions

$$\left\{ \begin{array}{l} u(\mathbf{x}, 0) = u^0 \quad \mathbf{x} \in \Omega, \\ u(\mathbf{x}, t) = u_g \quad (\mathbf{x}, t) \in \partial\Omega \times (0, t_{\text{max}}] \end{array} \right. \quad (3)$$

where $a$ is the diffusion coefficient and $\mathbf{v}$ the velocity field, $\Omega \subset \mathbb{R}^d$, $1 \le d \le 3$, $t_{\text{max}} > 0$. The aim of the separated representation method is to compute $N$ couples of functions $\{(X_i, \Theta_i)\}_{i=1,...,N}$ such that $\{X_i\}_{i=1,...,N}$ and $\{\Theta_i\}_{i=1,...,N}$ are defined respectively in $\bar{\Omega}$ and $[0, t_{\text{max}}]$ and the solution $u$ of this problem can be written in the separate form

$$u(\mathbf{x}, t) \approx \sum_{i=1}^{N} \Theta_i(t) \cdot X_i(\mathbf{x}) \quad (4)$$

The weak form of problem (2) yields:

Find $u(\mathbf{x}, t)$ verifying the boundary conditions (3) such that

$$\int_0^{t_{\text{max}}} \int_{\Omega} u^*(\frac{\partial u}{\partial t} - a \cdot \Delta u + \mathbf{v} \cdot \nabla u - f(\mathbf{x}, t)) \, d\mathbf{x} \, dt = 0 \quad (5)$$

for all the functions $u^*(x, t)$ in an appropriate functional space.

We compute now the functions involved in the sum (4). We suppose that the set of functional couples $\{(X_i, \Theta_i)\}_{i=1,...,n}$ with $0 \le n < N$ are already known (they have been computed at the previous iterations) and that at the present iteration we search the enrichment couple $(R(\mathbf{x}), S(t))$ by applying an alternating directions fixed point algorithm that after convergence will constitute the next functional couple $(X_{n+1}, \Theta_{n+1})$. Hence, at the present iteration, $n$, we assume the separated representation

$$u(\mathbf{x}, t) \approx \sum_{i=1}^{n} \Theta_i(t) \cdot X_i(\mathbf{x}) + S(t) \cdot R(\mathbf{x}) \quad (6)$$
The weighting function $u^\star$ is then assumed as

$$u^\star = S \cdot R^\star + R \cdot S^\star \qquad (7)$$

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

$$ \begin{aligned} & \int_0^{t_{\text{max}}} \int_{\Omega} (S \cdot R^\star + R \cdot S^\star) \cdot \left( R \cdot \frac{\partial S}{\partial t} - a \cdot \Delta R \cdot S + (\mathbf{v} \cdot \nabla R) \cdot S \right) \mathrm{d}\mathbf{x} \, \mathrm{d}t \\ & = \int_0^{t_{\text{max}}} \int_{\Omega} (S \cdot R^\star + R \cdot S^\star) \cdot f(\mathbf{x}, t) - \sum_{i=1}^n X_i \cdot \frac{\partial \Theta_i}{\partial t} + a \cdot \sum_{i=1}^n \Delta X_i \\ & \qquad \cdot \Theta_i - \sum_{i=1}^n (\mathbf{v} \cdot \nabla X_i) \cdot \Theta_i \right) \mathrm{d}\mathbf{x} \, \mathrm{d}t \end{aligned} \qquad (8) $$

We apply an alternating directions fixed point algorithm to compute the couple of functions $(R, S)$:

• Computing the function $R(\mathbf{x})$.

First, we suppose that $S$ is known, implying that $S^\star$ vanishes in (7). Thus, Eq. (8) writes

$$ \begin{aligned} & \int_{\Omega} R^\star \cdot (\alpha_t \cdot R - a \cdot \beta_t \cdot \Delta R + \beta_t \cdot \mathbf{v} \cdot \nabla R) \, \mathrm{d}\mathbf{x} = \\ & \int_{\Omega} R^\star \cdot \gamma_t(\mathbf{x}) - \sum_{i=1}^n \alpha_t^i \cdot X_i + a \cdot \sum_{i=1}^n \beta_t^i \cdot \Delta X_i - \sum_{i=1}^n \beta_t^i \cdot \mathbf{v} \cdot \nabla X_i \right) \, \mathrm{d}\mathbf{x} \end{aligned} \qquad (9) $$

where

$$ \left\{ \begin{aligned} \alpha_t &= \int_0^{t_{\text{max}}} S(t) \cdot \frac{\partial S}{\partial t}(t) \, \mathrm{d}t \\ \alpha_t^i &= \int_0^{t_{\text{max}}} S(t) \cdot \frac{\partial \Theta_i}{\partial t}(t) \, \mathrm{d}t \\ \beta_t &= \int_0^{t_{\text{max}}} S^2(t) \, \mathrm{d}t \\ \beta_t^i &= \int_0^{t_{\text{max}}} S(t) \cdot \Theta_i(t) \, \mathrm{d}t \\ \gamma_t(\mathbf{x}) &= \int_0^{t_{\text{max}}} S(t) \cdot f(\mathbf{x}, t) \, \mathrm{d}t; \quad \forall \mathbf{x} \in \Omega \end{aligned} \right. \qquad (10) $$

The weak formulation (9) is satisfied for all $S^\star$, therefore we could come back to the associated strong formulation

$$ \alpha_t \cdot R - a \cdot \beta_t \cdot \Delta R + \beta_t \cdot \mathbf{v} \cdot \nabla R = \gamma_t(\mathbf{x}) - \sum_{i=1}^n \alpha_t^i \cdot X_i + a \cdot \sum_{i=1}^n \beta_t^i \cdot \Delta X_i - \sum_{i=1}^n \beta_t^i \cdot \mathbf{v} \cdot \nabla X_i \qquad (11) $$

that one could solve by using any appropriate discretization technique.

• Computing the function $R(t)$.

From the function $R(\mathbf{x})$ just computed, we search $S(t)$. In this case $R^\star$ vanishes in (7) and (8) reduces to

$$ \begin{aligned} & \int_0^{t_{\text{max}}} \int_{\Omega} (R \cdot S^\star) \cdot \left( R \cdot \frac{\partial S}{\partial t} - a \cdot \Delta R \cdot S + (\mathbf{v} \cdot \nabla R) \cdot S \right) \mathrm{d}\mathbf{x} \, \mathrm{d}t \\ & = \int_0^{t_{\text{max}}} \int_{\Omega} (R \cdot S^\star) \cdot f(\mathbf{x}, t) - \sum_{i=1}^n X_i \cdot \frac{\partial \Theta_i}{\partial t} + a \cdot \sum_{i=1}^n \Delta X_i \cdot \Theta_i \\ & \qquad - \sum_{i=1}^n (\mathbf{v} \cdot \nabla X_i) \cdot \Theta_i \right) \mathrm{d}\mathbf{x} \, \mathrm{d}t \end{aligned} \qquad (12) $$

where all the spatial functions can be integrated in $\Omega$. Thus, by using the following notations

$$ \left\{ \begin{aligned} \alpha_x &= \int_{\Omega} R(\mathbf{x}) \cdot \Delta R(\mathbf{x}) \, \mathrm{d}\mathbf{x} \\ \alpha_x^i &= \int_{\Omega} R(\mathbf{x}) \cdot \Delta X_i(\mathbf{x}) \, \mathrm{d}\mathbf{x} \\ \beta_x &= \int_{\Omega} R^2(\mathbf{x}) \, \mathrm{d}\mathbf{x} \\ \beta_x^i &= \int_{\Omega} R(\mathbf{x}) \cdot X_i(\mathbf{x}) \, \mathrm{d}\mathbf{x} \\ \lambda_x &= \int_{\Omega} R(\mathbf{x}) \cdot (\mathbf{v} \cdot \nabla R(\mathbf{x})) \, \mathrm{d}\mathbf{x} \\ \lambda_x^i &= \int_{\Omega} R(\mathbf{x}) \cdot (\mathbf{v} \cdot \nabla X_i(\mathbf{x})) \, \mathrm{d}\mathbf{x}, \\ \gamma_x^*(t) &= \int_{\Omega} R(\mathbf{x}) \cdot f(\mathbf{x}, t) \, \mathrm{d}\mathbf{x}; \quad \forall t \end{aligned} \right. \qquad (13) $$

Eq. (12) reads

$$ \int_0^{t_{\text{max}}} S^\star \cdot \beta_x \cdot \frac{\partial S}{\partial t} + (\lambda_x - a \cdot \alpha_x) \cdot S - \gamma_x(t) + \sum_{i=1}^n \beta_x^i \cdot \frac{\partial \Theta_i}{\partial t} + \sum_{i=1}^n (\lambda_x^i - a \cdot \alpha_x^i) \cdot \Theta_i \right) \mathrm{d}t = 0 \qquad (14) $$

As Eq. (14) holds for all $S^\star$, we could come back to the strong formulation

$$ \beta_x \cdot \frac{\partial S}{\partial t} + (\lambda_x - a \cdot \alpha_x) \cdot S = \gamma_x(t) - \sum_{i=1}^n \beta_x^i \cdot \frac{\partial \Theta_i}{\partial t} - \sum_{i=1}^n (\lambda_x^i - a \cdot \alpha_x^i) \cdot \Theta_i \qquad (15) $$

which is a first order ordinary differential equation that can be solved easily (even for extremely small time steps) from its initial condition.

These two steps must be repeated until convergence, that is, until verifying that both functions reach a fixed point. If we denote by $S^{(q)}(t)$ and $S^{(q-1)}(t)$ the computed functions $S(t)$ at the present and previous iteration respectively, and the same for the space functions: $R^{(q)}(\mathbf{x})$ and $R^{(q-1)}(\mathbf{x})$, the stoping criterion used in this work writes:

$$ e = \|S^{(q)}(t) \cdot R^{(q)}(\mathbf{x}) - S^{(q-1)}(t) \cdot R^{(q-1)}(\mathbf{x})\|_2 < 10^{-8} \qquad (16) $$

where $10^{-8}$ represents the square root of the machine precision.

We denote by $Q_{n+1}$ the number of iterations for solving this nonlinear problem to determine the enrichment couple of functions $X_{n+1}(\mathbf{x})$ and $\Theta_{n+1}(t)$. After reaching convergence we write $X_{n+1}(\mathbf{x}) = R(\mathbf{x})$ and $\Theta_{n+1}(t) = S(t)$. The enrichment procedure must continue until reaching the convergence of the enrichment global procedure at iteration $N$, when the separated representation of the unknown field writes:

$$ u(\mathbf{x}, t) \approx \sum_{i=1}^{N} X_i(\mathbf{x}) \cdot \Theta_i(t) \qquad (17) $$

The more usual global stopping criteria are:

• For models whose exact solution $u^{\text{ref}}$ is known:

$$ E = \frac{\|u - u^{\text{ref}}\|_2}{\|u^{\text{ref}}\|_2} < \epsilon \qquad (18) $$

• For models whose exact solution is not known:

$$ E = \frac{\|\frac{\partial u}{\partial t} - a \cdot \Delta u + \mathbf{v} \cdot \nabla u - f(\mathbf{x}, t)\|_2}{\|f(\mathbf{x}, t)\|_2} < \epsilon \qquad (19) $$

with $\epsilon$ a small enough parameter ($\epsilon = 10^{-8}$ in our simulations and the $L_2$-norm applies in the whole space-time domain).

For the procedure of enforcing boundary conditions the interested reader can refer to [14]. For alternative more efficient separated constructors of the separated representation for non-symmetric differential operators the interested reader can refer
to [10]. Finally, error estimators based on quantities of interest allowing the definition of more efficient stopping criteria can be found in [4,17].

The just proposed strategy needs for the solution of about $N \times Q$ space and time problems (with $Q = (Q_1 + \cdots + Q_N)/N$ and $N$ the number of functional couples needed to approximate, up to the desired precision, the searched solution). Thus, one must compute $N \times Q$ d-D problems, $d = 1, 2, 3$, whose complexity depends on the spatial mesh considered, and also $N \times Q$ 1D problems (defined in the time interval $\mathcal{T} = (0, t_{max}]$) that only need the solution of an ordinary differential equation from its initial condition. Obviously, even for extremely small time steps, the solution of these transient 1D problems does not introduce major difficulties.

If instead of the separated representation just discussed, one performs a standard incremental solution, $P$ dD models, $d = 1, 2, 3$, must be solved ($P$ being the number of time steps, i.e. $P = t_{max}/\Delta t$, where the time step $\Delta t$ must be chosen for ensuring the stability conditions).

In all the analyzed cases $N$ and $Q$ are of the order of tens that implies the solution of about hundred $d$-dimensional problems defined in $\Omega$, instead of the thousands (or even millions) needed for solving those models using standard incremental solvers.

## 4. Separated representation of the Rayleigh-Bénard model solution

We consider the dimensionless form of the Rayleigh-Bénard model (see Appendix A for more details concerning the derivation of this model):

$$ egin{cases} 
\nabla \cdot \mathbf{v} = 0 \\ 
\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \theta \cdot \mathbf{j} \\ 
\frac{\partial \theta}{\partial t} + \mathbf{v} \cdot \nabla \theta = (P_r \cdot R_a)^{-\frac{1}{2}} \nabla^2 \theta \\ 
\boldsymbol{\tau} = P_r^{\frac{n}{2}} \cdot R_a^{\frac{n}{2}-1} \cdot (D_{eq})^{n-1} \cdot \mathbf{D} 
\end{cases} 
\quad (20) $$

where $\mathbf{v}$ is the dimensionless velocity, $p$ is the dimensionless pressure, $\theta$ the dimensionless temperature, $\mathbf{D}$ the dimensionless strain rate tensor (symmetric component of the dimensionless velocity gradient), $D_{eq}$ the equivalent strain rate that depends on the second invariant of the strain rate tensor, $P_r$ and $R_a$ the Prandtl and Rayleigh dimensionless numbers respectively (their expressions are given in the appendix A), $\boldsymbol{\tau}$ the deviatoric part of the dimensionless Cauchy's stress tensor and $\mathbf{j}$ the unit vector defining the $y$-direction along which gravity force applies.

The model is defined in the unit square domain $\Omega = (0,1) \times (0,1)$ (related to the square domain of size $H$ depicted in Fig. 1) fulfilled by a fluid (modeled by a power-law constitutive equation characterized by a power index $n$) initially at rest. A dimensionless temperature $\theta(x,y=0) = \theta_H = 0.5$ is enforced at the bottom boundary $y=0$, whereas a low dimensionless temperature $\theta(x,y=1) = \theta_C = -0.5$ is applied on the upper boundary. The heat flux is assumed vanishing on the left and right domain boundaries. The initial temperature distribution compatible with a thermal conduction regime in a fluid at rest, evolves linearly in the $y$-direction:

$$ \theta(x, y, t = 0) = \theta_H - y \quad (21) $$

Eq. (20) defines a mixed formulation involving as unknowns the dimensionless velocity, pressure and temperature fields. When we proceed to the discretization of the weak form related to Eq. (20) some stability conditions must be ensured. One of them concerns the so-called LBB condition that restricts the free choice of pressure and velocity approximations. Because in this work we are considering separated representations of the unknown fields and the question related to the expression of

<figure><img src="image_3.png" /><figcaption>Fig. 1. Square domain and boundary conditions.</figcaption></figure>

LBB conditions within a separated representation framework is not today fully understood, we are considering a penalty formulation of the incompressibility constraint. Thus, the mass balance we are considering writes:

$$ \nabla \cdot \mathbf{v} + \frac{1}{\lambda} \cdot p = 0 \quad (22) $$

with $\lambda$ a large enough constant. This expression implies:

$$ p = -\lambda \cdot (\nabla \cdot \mathbf{v}) \quad (23) $$

that substituted in the momentum equation simplifies the model formulation

$$ egin{cases} 
\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} = \lambda \cdot \nabla (\nabla \cdot \mathbf{v}) + \nabla \cdot \boldsymbol{\tau} + \theta \cdot \mathbf{j} \\ 
\frac{\partial \theta}{\partial t} + \mathbf{v} \cdot \nabla \theta = (P_r \cdot R_a)^{-\frac{1}{2}} \nabla^2 \theta \\ 
\boldsymbol{\tau} = P_r^{\frac{n}{2}} \cdot R_a^{\frac{n}{2}-1} \cdot (D_{eq})^{n-1} \cdot \mathbf{D} 
\end{cases} 
\quad (24) $$

which only implies the dimensionless velocity and temperature fields. Note that the penalty term is not introduced for simplifying the model, but only for circumventing the lack of knowledge concerning the expression of LBB conditions in the framework of separated representations.

The weighted residual form related to Eq. (24) reads:

$$ egin{cases} 
\int_{\Omega \times \mathcal{S}} \mathbf{v}^* \cdot \left\{ \frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} - \lambda \cdot \nabla (\nabla \cdot \mathbf{v}) - \nabla \cdot \boldsymbol{\tau} - \theta \cdot \mathbf{j} \right\} \cdot d\mathbf{x} \cdot dt = 0 \\ 
\int_{\Omega \times \mathcal{S}} \theta^* \cdot \left\{ \frac{\partial \theta}{\partial t} + \mathbf{v} \cdot \nabla \theta - (P_r \cdot R_a)^{-\frac{1}{2}} \nabla^2 \theta \right\} \cdot d\mathbf{x} \cdot dt = 0 
\end{cases} 
\quad (25) $$

where $\boldsymbol{\tau} = P_r^{\frac{n}{2}} \cdot R_a^{\frac{n}{2}-1} \cdot (D_{eq})^{n-1} \cdot \mathbf{D}$ and $\mathcal{S}$ represents the dimensionless time interval.

In what follows we distinguish two cases, the one involving a newtonian fluid characterized by a unit power index, i.e. $n = 1$ and the second one related to non-Newtonian fluids characterized by a non unit power index, i.e. $n \neq 1$. The rheo-thinning fluids are related to $n < 1$, whereas $n > 1$ in the case of rheo-thickening fluids.

### 4.1. Newtonian fluids

In the case of Newtonian fluids $n = 1$ and the constitutive equation reduces to $\boldsymbol{\tau} = \tilde{\boldsymbol{\eta}} \cdot \mathbf{D}$, with $\tilde{\boldsymbol{\eta}} = (P_r(n=1))^{\frac{1}{2}} \cdot (R_a(n=1))^{-\frac{1}{2}}$.

Thus, the strong and weak form of the Rayleigh-Bénard model read:

$$ egin{cases} 
\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} = \lambda \cdot \nabla (\nabla \cdot \mathbf{v}) + \tilde{\eta} \cdot \nabla \cdot \mathbf{D} + \theta \cdot \mathbf{j} \\ 
\frac{\partial \theta}{\partial t} + \mathbf{v} \cdot \nabla \theta = \tilde{\alpha} \cdot \nabla^2 \theta 
\end{cases} 
\quad (26) $$

where $\tilde{\alpha} = (P_r(n=1) \cdot R_a(n=1))^{-\frac{1}{2}}$, and
$$ \begin{cases} \int_{\Omega \times \mathcal{F}} \mathbf{v}^* \cdot \left\{ \frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} - \lambda \cdot \nabla (\nabla \cdot \mathbf{v}) - \tilde{\eta} \cdot \nabla \cdot \mathbf{D} - \theta \cdot \mathbf{j} \right\} \cdot \mathbf{d}\mathbf{x} \cdot \mathbf{d}t = 0 \\ \int_{\Omega \times \mathcal{F}} \theta^* \cdot \left\{ \frac{\partial \theta}{\partial t} + \mathbf{v} \cdot \nabla \theta - \tilde{\alpha} \cdot \nabla^2 \theta \right\} \cdot \mathbf{d}\mathbf{x} \cdot \mathbf{d}t = 0 \end{cases} \quad (27) $$

## 4.2. Power-law fluids

In the case of power-law fluids both formulations, strong and weak, read:

$$ \left\{ \begin{array}{l} \frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} = \lambda \cdot \nabla (\nabla \cdot \mathbf{v}) + \nabla \cdot (\tilde{\eta} \cdot \mathbf{D}) + \theta \cdot \mathbf{j} \\ \frac{\partial \theta}{\partial t} + \mathbf{v} \cdot \nabla \theta = \tilde{\alpha} \cdot \nabla^2 \theta \end{array} \right. \quad (28) $$

where $\tilde{\alpha} = (P_r(n) \cdot R_a(n))^{-\frac{1}{2}}$ and $\tilde{\eta} = (P_r(n))^\frac{1}{2} \cdot (R_a(n))^\frac{1}{2} \cdot (D_{eq})^{n-1}$, and

$$ \left\{ \begin{array}{l} \int_{\Omega \times \mathcal{F}} \mathbf{v}^* \cdot \left\{ \frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} - \lambda \cdot \nabla (\nabla \cdot \mathbf{v}) - \nabla \cdot (\tilde{\eta} \cdot \mathbf{D}) - \theta \cdot \mathbf{j} \right\} \cdot \mathbf{d}\mathbf{x} \cdot \mathbf{d}t = 0 \\ \int_{\Omega \times \mathcal{F}} \theta^* \cdot \left\{ \frac{\partial \theta}{\partial t} + \mathbf{v} \cdot \nabla \theta - \tilde{\alpha} \cdot \nabla^2 \theta \right\} \cdot \mathbf{d}\mathbf{x} \cdot \mathbf{d}t = 0 \end{array} \right. \quad (29) $$

## 4.3. Separated representation of velocity and temperature fields

With the components of the velocity field $\mathbf{v}$ denoted by $(u, v)$, the separated representation of the different unknown fields read:

$$ u(x, y, t) \approx \sum_{i=1}^{i=N} X_i^u(x, y) \cdot \Theta_i^u(t) \quad (30) $$

$$ v(x, y, t) \approx \sum_{i=1}^{i=N} X_i^v(x, y) \cdot \Theta_i^v(t) \quad (31) $$

and

$$ \theta(x, y, t) \approx \sum_{i=1}^{i=N} X_i^\theta(x, y) \cdot \Theta_i^\theta(t) \quad (32) $$

Expression (30) and (31) can be written in the compact separated representation vector form

$$ \mathbf{v}(x, y, t) \equiv \begin{pmatrix} u(x, y, t) \\ v(x, y, t) \end{pmatrix} \approx \begin{pmatrix} \sum_{i=1}^{i=N} X_i^u(x, y) \cdot \Theta_i^u(t) \\ \sum_{i=1}^{i=N} X_i^v(x, y) \cdot \Theta_i^v(t) \end{pmatrix} \equiv \sum_{i=1}^{i=N} \mathbf{X}_i(x, y) \circ \mathbf{\Theta}_i(t) \quad (33) $$

where the symbol “◦” denotes the so-called entry-wise, Hadamard or Schur multiplication for vectors. Thus for two generic vectors $\mathbf{a}$ and $\mathbf{b}$, the i-component of the entry-wise product $(\mathbf{a} \circ \mathbf{b})_i$ is given by $(\mathbf{a} \circ \mathbf{b})_i = \mathbf{a}_i \cdot \mathbf{b}_i$.

This separated representation is built as described in Section 3, by computing a term at each iteration of the PGD constructor. Thus, if we assume at iteration $m$ the solution $\mathbf{v}^m$ and $\theta^m$ given by:

$$ \mathbf{v}^m(x, y, t) = \sum_{i=1}^{i=m} \mathbf{X}_i(x, y) \circ \mathbf{\Theta}_i(t) \quad (34) $$

and

$$ \theta^m(x, y, t) = \sum_{i=1}^{i=m} X_i^\theta(x, y) \cdot \Theta_i^\theta(t) \quad (35) $$

at iteration $m+1$ we look for the next functional products:

$$ \begin{aligned} \mathbf{v}^{m+1}(x, y, t) &= \sum_{i=1}^{i=m} \mathbf{X}_i(x, y) \circ \mathbf{\Theta}_i(t) + \mathbf{R}(x, y) \circ \mathbf{S}(t) \\ &= \mathbf{v}^m(x, y, t) + \mathbf{R}(x, y) \circ \mathbf{S}(t) \end{aligned} \quad (36) $$

and

$$ \begin{aligned} \theta^{m+1}(x, y, t) &= \sum_{i=1}^{i=m} X_i^\theta(x, y) \cdot \Theta_i^\theta(t) + R^\theta(x, y) \cdot S^\theta(t) \\ &= \theta^m(x, y, t) + R^\theta(x, y) \cdot S^\theta(t) \end{aligned} \quad (37) $$

## 4.4. Linearization

The previous models are always non-linear because the advective terms $\mathbf{v} \cdot \nabla \mathbf{v}$ and $\mathbf{v} \cdot \nabla \theta$, and coupled. Moreover, in the case of power-law fluids there is a second non-linearity coming from the fluid constitutive behavior. As we compute the transient solution in a non-incremental way, i.e. we compute simultaneously all the time history, the simplest linearization consists of linearizing at iteration $m+1$ the non-linear terms around the solution at the previous iteration $m$. Thus, we can write at iteration $m+1$:

$$ \mathbf{v}^{m+1} \cdot \nabla \mathbf{v}^{m+1} \approx \mathbf{v}^m \cdot \nabla \mathbf{v}^{m+1} \quad (38) $$

and

$$ \mathbf{v}^{m+1} \cdot \nabla \theta^{m+1} \approx \mathbf{v}^m \cdot \nabla \theta^{m+1} \quad (39) $$

and in the case of power-law fluids:

$$ \tilde{\eta}(\mathbf{v}^{m+1}) \approx (P_r(n))^\frac{1}{2} \cdot (R_a(n))^\frac{1}{2} \cdot (D_{eq}(\mathbf{v}^m))^{n-1} \quad (40) $$

This simple linearization is not optimal as explained in [3], because the number of terms in the decomposition depends on the convergence rate of the fixed point algorithm and not only on the separability of the approximated solution. An enhanced linearization consists of considering a better approximation of $\mathbf{v}^{m+1}$, $\tilde{\mathbf{v}}^{m+1}$, appearing in the non-linear terms:

$$ \mathbf{v}^{m+1} \cdot \nabla \mathbf{v}^{m+1} \approx \tilde{\mathbf{v}}^{m+1} \cdot \nabla \mathbf{v}^{m+1} = (\mathbf{v}^m + \mathbf{R}^{k-1} \circ \mathbf{S}^{k-1}) \cdot \nabla (\mathbf{v}^m + \mathbf{R}^k \circ \mathbf{S}^k) \quad (41) $$

where $k$ refers to the iteration of the non-linear solver that applies at each iteration of the separated representation constructor. Similarly we have:

$$ \mathbf{v}^{m+1} \cdot \nabla \theta^{m+1} \approx \tilde{\mathbf{v}}^{m+1} \cdot \nabla \theta^{m+1} \quad (42) $$

and in the case of power-law fluids:

$$ \tilde{\eta}(\mathbf{v}^{m+1}) \approx (P_r(n))^\frac{1}{2} \cdot (R_a(n))^\frac{1}{2} \cdot (D_{eq}(\tilde{\mathbf{v}}^{m+1}))^{n-1} \quad (43) $$

This enhanced formulation allows to computed decompositions involving less terms because as it can be noticed, see [3] for a more detailed discussion, the formulation is much less sensitive to the non-linear solver convergence rate. In what follows we are using the simplest strategy, the one considering explicitly $\mathbf{v}^m$ in the non-linear terms.

As just discussed, in general, the separated representation of the solution is not optimal. In order to reduce the number of terms we can apply a sort of post-treatment based again on the PGD. Thus, if the computed solution at convergence writes

$$ \left\{ \begin{aligned} \mathbf{v}^{\text{PGD}}(x, y, t) &\approx \sum_{i=1}^{N} \mathbf{X}_i(x, y) \circ \mathbf{\Theta}_i(t) \\ \theta^{\text{PGD}}(x, y, t) &\approx \sum_{i=1}^{N} X_i^\theta(x, y) \cdot \Theta_i^\theta(t) \end{aligned} \right. \quad (44) $$

we can try to reduce the number of terms involved in the decomposition by looking for $\mathbf{v}$ and $\theta$ in a separated form such that:

$$ \left\{ \begin{aligned} \int_{\Omega \times \mathcal{F}} \mathbf{v}^* \cdot (\mathbf{v} - \mathbf{v}^{\text{PGD}}(x, y, t)) \cdot \mathbf{d}\mathbf{x} \cdot \mathbf{d}t &= 0 \\ \int_{\Omega \times \mathcal{F}} \theta^* \cdot (\theta - \theta^{\text{PGD}}(x, y, t)) \cdot \mathbf{d}\mathbf{x} \cdot \mathbf{d}t &= 0 \end{aligned} \right. \quad (45) $$
The resulting solutions involves in general less terms that the original one, because (45) corresponds with the optimal one related to a standard SVD decomposition.

### 4.5. Separated representation based weak form

Thus, by considering the simplest linearization and the approximations at iteration m:

$$ egin{cases} 
\mathbf{v}^m(x, y, t) \approx \sum_{i=1}^{m} \mathbf{X}_i(x, y) \circ \mathbf{\Theta}_i(t) \\ 
\theta^m(x, y, t) \approx \sum_{i=1}^{m} \mathbf{X}_i^\theta(x, y) \circ \mathbf{\Theta}_i^\theta(t) 
\end{cases} 
\quad (46) $$

<figure><img src="image_4.png" /><figcaption>Fig. 2. Comparison of $u_{max}(t)$ (left) and $v_{max}(t)$ (right) obtained from the PGD and FEM solvers ($R_a = 10^4$)</figcaption></figure>

at iteration m + 1 the searched solution writes:

$$ egin{cases} 
\mathbf{v}^{m+1}(x, y, t) \approx \sum_{i=1}^{m} \mathbf{X}_i \circ \mathbf{\Theta}_i + \mathbf{R} \circ \mathbf{S} \\ 
\theta^{m+1}(x, y, t) \approx \sum_{i=1}^{m} \mathbf{X}_i^\theta \circ \mathbf{\Theta}_i^\theta + \mathbf{R}^\theta \circ \mathbf{S}^\theta 
\end{cases} 
\quad (47) $$

where for the sake of clarity we do not specify the dependence of functions $\mathbf{X}_i$, $\mathbf{X}_i^\theta$, $\mathbf{R}$ and $\mathbf{R}^\theta$ on $(x, y)$ and $\mathbf{\Theta}, \mathbf{\Theta}^\theta, \mathbf{S}$ and $\mathbf{S}$ on $t$. The test functions related to (47) write:

$$ egin{cases} 
\mathbf{v}^*(x, y, t) = \mathbf{R}^* \circ \mathbf{S} + \mathbf{R} \circ \mathbf{S}^* \\ 
\theta^*(x, y, t) = (\mathbf{R}^\theta)^* \circ \mathbf{S}^\theta + \mathbf{R}^\theta \circ (\mathbf{S}^\theta)^* 
\end{cases} 
\quad (48) $$

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

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

<figure><img src="image_7.png" /><figcaption>Fig. 3. $u(x = 0.5, y)$ (left) and $\theta(x = 0.5, y)$ (right) at different times ($R_a = 10^4$)</figcaption></figure>

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

<figure><img src="image_9.png" /><figcaption>Fig. 4. $u(x = 0.5, y)$ (left) and $\theta(x = 0.5, y)$ (right) at different times ($R_a = 10^5$)</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 5. $u(x=0.5,y)$ (left) and $\theta(x=0.5,y)$ (right) at different times ($R_a = 10^6$)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 6. Steady state velocity (left) and temperature (right) fields for different $R_a$ numbers: $R_a = 10^4$ (top), $R_a = 10^5$ (center) and $R_a = 10^6$ (bottom)</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 7. Velocity and temperature fields associated with $R_a = 10^6$ at three different times instants: $t = 0.3$ (top), $t = 18$ (center) and $t = 55$ (bottom) at which the steady state is almost reached.</figcaption></figure>

By introducing Eqs. (47) and (48) into the general weak form (29) it results:

$$ \begin{aligned} & \int_{\Omega \times \mathcal{J}} (\mathbf{R}^+ \circ \mathbf{S} + \mathbf{R} \circ \mathbf{S}^+) \cdot \left\{ \sum_{i=1}^{i=m} \mathbf{X}_i \circ \frac{d\mathbf{\Theta}_i}{dt} + \mathbf{R} \circ \frac{d\mathbf{S}}{dt} + \sum_{i=1}^{i=m} \mathbf{X}_i \circ \mathbf{\Theta}_i \right\} \\ & \qquad \cdot \left( \sum_{i=1}^{i=m} \nabla(\mathbf{X}_i \circ \mathbf{\Theta}_i) + \nabla(\mathbf{R} \circ \mathbf{S}) \right) - \lambda \cdot \sum_{i=1}^{i=m} \nabla(\nabla \circ (\mathbf{X}_i \circ \mathbf{\Theta}_i)) \\ & \qquad - \lambda \cdot \nabla(\nabla \circ (\mathbf{R} \circ \mathbf{S})) - \nabla \tilde{\eta} \cdot \left( \sum_{i=1}^{i=m} \mathbf{D}_i + \mathbf{D} \right) - \tilde{\eta} \cdot \sum_{i=1}^{i=m} \nabla \cdot \mathbf{D}_i + \nabla \cdot \mathbf{D} \\ & \qquad - \sum_{i=1}^{i=m} \mathbf{X}_i^\theta \circ \mathbf{\Theta}_i^\theta + \mathbf{R}^\theta \cdot \mathbf{S}^\theta \right) \cdot \mathbf{j} \right\} \cdot d\mathbf{x} \cdot dt = 0 \quad (49) \end{aligned} $$

and

$$ \begin{aligned} & \int_{\Omega \times \mathcal{J}} ((\mathbf{R}^\theta)^+ \cdot \mathbf{S}^\theta + \mathbf{R}^\theta \cdot (\mathbf{S}^\theta)^+) \cdot \left\{ \sum_{i=1}^{i=m} \mathbf{X}_i^\theta \circ \frac{d\mathbf{\Theta}_i^\theta}{dt} + \mathbf{R}^\theta \cdot \frac{d\mathbf{S}^\theta}{dt} + \sum_{i=1}^{i=m} \mathbf{X}_i \circ \mathbf{\Theta}_i \right\} \\ & \qquad \cdot \left( \sum_{i=1}^{i=m} \nabla \mathbf{X}_i^\theta \circ \mathbf{\Theta}_i^\theta + \nabla \mathbf{R}^\theta \cdot \mathbf{S}^\theta \right) - \tilde{\alpha} \cdot \sum_{i=1}^{i=m} \nabla^2 \mathbf{X}_i^\theta \circ \mathbf{\Theta}_i^\theta + \nabla^2 \mathbf{R}^\theta \cdot \mathbf{S}^\theta \right\} \cdot d\mathbf{x} \\ & \qquad \cdot dt = 0 \end{aligned} \quad (50) $$

where

$$ \left\{ \begin{aligned} \mathbf{D}_i &= \frac{1}{2} (\nabla (\mathbf{X}_i \circ \mathbf{\Theta}_i) + (\nabla (\mathbf{X}_i \circ \mathbf{\Theta}_i))^T) \\ \mathbf{D} &= \frac{1}{2} (\nabla (\mathbf{R} \circ \mathbf{S}) + (\nabla (\mathbf{R} \circ \mathbf{S}))^T) \end{aligned} \right. \quad (51) $$

If we consider a vector **A** with components ($A_1, A_2$) depending on $(x, y)$ and a vector **B**, with components ($B_1, B_2$) depending on $t$, the gradient differential operator applying on the Hadamard product **A** ◦ **B** writes:
<figure><img src="image_13.png" /><figcaption>Fig. 8. Local Nusselt's number at different times and for different values of $R_a$: $R_a = 10^4$ (top), $R_a = 10^5$ (center) and $R_a = 10^6$ (bottom)</figcaption></figure>

$$ \nabla(\mathbf{A} \circ \mathbf{B}) = \frac{\partial A_1}{\partial x} \cdot B_1 \frac{\partial A_1}{\partial y} \cdot B_1 = \left( \frac{\partial A_1}{\partial x} \frac{\partial A_1}{\partial y} \right) \circ \begin{pmatrix} B_1 & B_1 \\ B_2 & B_2 \end{pmatrix} = \nabla \mathbf{A} \circ \hat{\mathbf{B}} \quad (52) $$

It is easy to verify that $(\nabla(\mathbf{A} \circ \mathbf{B}))^T = (\nabla \mathbf{A})^T \circ \hat{\mathbf{B}}^T$.

Finally when the divergence operator applies on tensors $\mathbf{D}$ and $\mathbf{D}_i$ it results:

$$ 2 \cdot \nabla \cdot \mathbf{D} = \nabla \cdot (\nabla(\mathbf{R} \circ \mathbf{S}) + (\nabla(\mathbf{R} \circ \mathbf{S}))^T) = \nabla^2 \mathbf{R} \circ \hat{\mathbf{S}} + \nabla \cdot (\nabla(\mathbf{R} \circ \mathbf{S}))^T \quad (53) $$

and

$$ 2 \cdot \nabla \cdot \mathbf{D}_i = \nabla \cdot (\nabla(\mathbf{X}_i \circ \mathbf{\Theta}_i + (\nabla(\mathbf{X}_i \circ \mathbf{\Theta}_i))^T) = \nabla^2 \mathbf{X}_i \circ \hat{\mathbf{\Theta}}_i + \nabla \cdot (\nabla(\mathbf{X}_i \circ \mathbf{\Theta}_i))^T \quad (54) $$

The last terms in the right hand members of Eqs. (53) and (54) cannot be written in a compact form, however, the material incompressibility implies

$$ \begin{aligned} & \nabla \cdot \left( \sum_{i=1}^{i=m} (\nabla(\mathbf{X}_i \circ \mathbf{\Theta}_i))^T + (\nabla(\mathbf{R} \circ \mathbf{S}))^T \right) = \nabla \cdot \nabla \mathbf{v}^{m+1} \\ & = \nabla (\nabla \cdot \mathbf{v}^{m+1}) = \mathbf{0} \end{aligned} \quad (55) $$

an consequently term $\sum_{i=1}^{i=m} \nabla \cdot \mathbf{D}_i + \nabla \cdot \mathbf{D}$ in Eq. (49) reduced to:

$$ \sum_{i=1}^{i=m} \nabla \cdot \mathbf{D}_i + \nabla \cdot \mathbf{D} = \frac{1}{2} \sum_{i=1}^{i=m} \nabla^2 \mathbf{X}_i \circ \hat{\mathbf{\Theta}}_i + \nabla^2 \mathbf{R} \circ \hat{\mathbf{S}} \quad (56) $$

### 4.6. Fixed point alternated direction linearization

Now for computing functions $\mathbf{R}$, $R^\theta$, $\mathbf{S}$ and $S^\theta$, we proceed as described in Section 3, by applying a fixed point alternated direction strategy, that starting from an arbitrary $\mathbf{S}$ and $S^\theta$ computes $\mathbf{R}$ and $R^\theta$ from Eqs. (49) and (50). Then from the just calculated couple of functions we can update functions $\mathbf{S}$ and $S^\theta$. Both steps repeat until reaching the fixed point, i.e. until the just computed functions are close enough to the previous ones.

It must be highlighted that when functions $\mathbf{S}$ and $S^\theta$ are known, time integrals in Eqs. (49) and (50) can be performed leading to a weak form that only involved the fields $\mathbf{R}$, $R^\theta$ and theirs derivatives. Obviously, because the original problem involved second order space derivatives of the velocity components and the temperature field, the resulting weak form will involve second order space derivatives on both, the components of $\mathbf{R}$ and $R^\theta$. Obviously, integration by parts can by then applied in order to reduce the derivatives order and proceed to the discretization of the weak form by using standard continuous finite element interpolations for example, but other choices exist. A strong form could be derived, as shown in Section 3, and then solved by using any collocation technique applying on it.

On the other hand, when functions $\mathbf{R}$ and $R^\theta$ are assumed known, space integrals in Eqs. (49) and (50) can be performed, leading to a weak form that only concerns time functions $\mathbf{S}$, $S^\theta$ and their first time derivatives. Then two possibilities exist: (i) solving the weak form by using some stabilized discretization (e.g. discontinuous Galerkin) or (ii) coming back to the associated strong form, as described in Section 3, and then apply any forward discretization technique on the resulting ODE (e.g. Euleur or Runge-Kutta schemes).

<figure><img src="image_14.png" /><figcaption>Fig. 9. $u(x = 0.5, y)$ (left) and $\theta(x = 0.5, y)$ (right) for different values of the power index $n$ ($R_a = 10^4$)</figcaption></figure>
**5. Numerical results**

In what follows we are applying the technique just described to the solution of the Rayleigh–Bénard problem for both Newtonian and non-Newtonian fluids.

For approximating functions depending on the physical space $(x,y): \mathbf{X}_i, \mathbf{R}, \mathbf{X}_i^\theta$ and $\mathcal{R}^\theta$ we are using standard 8-nodes quadrilateral $\mathcal{R}^\theta$ finite elements. The considered computational meshes of the square cavity involve 1160 and 2296 nodes. The finer mesh was used to conclude on the convergence of the results computed on the coarser one. The time step considered to discretize the dimensionless time interval will be 0.1. Convergence in time was also checked by considering finer time steps. Convergence is assumed reached when the L2 norm of the error attains the value $\epsilon = 10^{-6}$.

At the initial time fluids filling the cavity were assumed at rest and the temperature distribution varies linearly with the y-coordinate, as previosuly mentioned: $\theta(x,y,t=0) = 0.5 - y$.

**5.1. Newtonian fluid**

First we consider the air filling the whole cavity, choice that implies $P_r = 0.71$. The Rayleigh–Bénard problem is then solved for different increasing values of $R_a$: $10^4$, $10^5$ and $10^6$. The PGD solution involves 20 terms, that is, $N = 20$ in the separated approximation of the velocity components $(u, v)$ and the temperature $\theta$ given by Eqs. (30)–(32) sufficed for attaining the desired accuracy.

In order to validate the obtained solution, we compared our results after reaching the flow steady state with the ones reported in [21] in the same conditions. The differences were of order $10^{-4}$ for the three values of the Rayleigh number $R_a$, proving the accuracy of the PGD solution. It is important to mention that each transient solution was computed in 15 min when using the PGD solver, instead the 2 days computation (2880 min) when using the finite

elements with equivalent meshes and time step. Obviously, by using better time integrations and more powerful computational resources, finite element solutions can be obtained in reasonable computing times, much lower than the 2 days just indicated. This comparison was only to justify the impressive reduction accomplished when using non-incremental discretizations because in the solutions reported below the whole transient solution was obtained by solving about 20 two-dimensional problems ($N = 20$) and transient solutions computed by applying standard incremental discretizations must solve a two-dimensional problem at each time step. Thus, one must be careful when applying standard techniques in order to reduce as much as possible the number of time steps, by using higher order discretizations, implicit integrations schemes, etc. However, when using space-time separated representations the simplest numerical choices allows computing time savings of many order of magnitude.

Fig. 2 compares the time evolution of the maximum value of both velocity components, $u_{max}(t)$ and $v_{max}(t)$ computed by using both the PGD and standard finite elements. We can notice that both evolutions agree in minute (the maximum difference being lower than 1%).

Figs. 3–5 depict $u(x=0.5,y)$ (left) and $\theta(x=0.5,y)$ (right) at different times and for different values of the $R_a$ number ($R_a = 10^4$ in Fig. 3, $R_a = 10^5$ in Fig. 4 and $R_a = 10^6$ in Fig. 5). We can see the gradual evolution from the initial state: $u(x=0.5,y,t=0)=0$ and $\theta(x=0.5,y)=0.5-y$; to the almost steady state (when the different curves almost superpose).

From the comparison of these results we can conclude that the velocity field becomes more complex as $R_a$ increases, and that this increase induces thermal boundary layers in the vicinity of the upper and lower boundaries. The steady state velocity and temperature fields are depicted in Fig. 6 ($R_a = 10^4$ (top), $R_a = 10^5$ (center) and $R_a = 10^6$ (bottom)).

<figure><img src="image_15.png" /><figcaption>**Fig. 10.** $u(x=0.5,y)$ (left) and $\theta(x=0.5,y)$ (right) for different values of the power index $n$ ($R_a = 10^5$)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 11.** $u(x=0.5,y)$ (left) and $\theta(x=0.5,y)$ (right) for different values of the power index $n$ ($R_a = 10^6$)</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 12.** Local Nusselt's number at different times and for different values of $R_a$: $R_a = 10^4$ (top), $R_a = 10^5$ (center) and $R_a = 10^6$ (bottom)</figcaption></figure>

In order to appreciate the time evolution of those fields we depict in Fig. 7 the velocity and temperature fields associated with $R_a = 10^6$ at three different times instants: $t = 0.3$ (top), $t = 18$

(center) and $t = 55$ (bottom) at which the steady state is almost reached.

Finally, in order to quantify the thermal efficiency, we define the local Nusselt's number on the bottom boundary $N_u(x,y = 0, t) = -\frac{\partial\theta}{\partial y}|_{y=0}$. In Fig. 8 we depict $N_u(x,y = 0)$ for the three values of the $R_a$ number: $R_a = 10^4$ (top), $R_a = 10^5$ (center) and $R_a = 10^6$ (bottom); each one at different time instants. From the comparison of these results we can conclude the increase of the thermal transfer by increasing the value of $R_a$, because the mean value of the Nusselt's number on the bottom wall increases from the conductive regime at $t = 0$ characterized by a constant unit value of the Nusselt's number, i.e. $N_u(-x,y = 0, t = 0) = 1$, to a mean value greater than one and increasing with $R_a$. Moreover, this mean value increases in time until reaching a maximum value and then it decreases lightly until reaching its steady state value.

## 5.2. Power-law fluids

In this section we consider power-law fluids characterized by different values of the power index, from $n = 0.7$ to $n = 1.8$. We consider $P_r = 7$ and three different values of $R_a$: $R_a = 10^4$, $R_a = 10^5$ and $R_a = 10^6$. In Fig. 9 we depict the steady state velocity $u(x = 0.5,y)$ and temperature $\theta(x = 0.5,y)$ for the different values of the power index in the case of $R_a = 10^4$. For the largest values of the power index (rheo-thickening fluid with $n > 1$) the fluid remains almost at rest and the temperature profile is very close to the initial one $\theta(x,y,t = 0) = 0.5 - y$. When the power index decreases $n < 1$, characterizing rheo-thinning fluids, the variation of both the velocity and the temperature fields is quite significant. Figs. 10 and 11 depict similar results for $R_a = 10^5$ and $R_a = 10^6$ respectively. We van notice that with the increase of $R_a$ the perturbation is the more and more significant and that both, the velocity and temperature fields localize the more and more in the upper and bottom walls neighborhood.

Fig. 12 depicts the steady state $N_u(x,y = 0)$ for the three values of $R_a$: $R_a = 10^4$ (top), $R_a = 10^5$ (center) and $R_a = 10^6$ (bottom); each one for different values of the power index. The mean value of $N_u(x,y = 0)$ increases with $R_a$ and with the decrease of the power index. Thus, rheo-thinning behavior increase the thermal transfer due to the enhanced convection effects.

Fig. 13 summarizes these facts by showing the evolution of the maximum velocity $u_{max}$ at the steady state (left) and the steady state mean value of the Nusselt's number on the bottom wall for different values of $R_a$ and the power index $n$.

In order to appreciate the time evolution of both $u_{max}(t)$ and the mean value of the Nusselt number $\overline{N_u}(t)$ we depicts in Fig. 14 these evolutions for $R_a = 10^5$ and three different values

<figure><img src="image_18.png" /><figcaption>**Fig. 13.** $u_{max}$ at the steady state (left) and the steady state mean value of the Nusselt's number on the bottom wall for different values of $R_a$ and the power index $n$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 14. $u_{max}(t)$ (left) and $\overline{N}_u(t)$ (right) for $R_a = 10^5$ and different values of the power index: $n = 0.8$ (top), $n = 1$ (center) and $n = 1.2$ (bottom).</figcaption></figure>

of the power index: $n = 0.8$ (top), $n = 1$ (center) and $n = 1.2$ (bottom). We can notice that the mean Nusselt's number increases rapidly, to reach a maximum and then it decreases until reaching its steady state.

Finally Fig. 15 depicts the steady state velocity (left) and temperature (right) fields for $R_a = 10^5$ and three different values of the power index $n$: $n = 0.8$ (top), $n = 1$ (center) and $n = 1.2$ (bottom). It can be noticed that the flow and thermal perturbations increase as the power index decreases (rheo-thinning behavior).

## 6. Conclusions

In this work we considered successfully the transient solution of non-linear coupled models related to the Rayleigh-Bénard flow model of both Newtonian and non-Newtonian fluids, by using the Proper Generalized Decomposition - PGD -. PGD proceeds by building-up a space-time separated representation of the different unknown fields, in our case the two components of the velocity field and the temperature. About 20 terms were needed to represent these fields, number that implied the necessity of solving

few tens of two-dimensional problems, instead the thousands required when using standard incremental discretizations. Thus, significant computing time savings were noticed.

This work constitutes a first attempt of using PGD decompositions for addressing complex flows of complex fluids, opening exciting perspectives concerning the solution of the Rayleigh-Bénard parametric models, involving pseudo-doplastic but also viscoelastic fluids and the increase of the Rayleigh numbers in order to approach and even move beyond the laminar-turbulent transition.

## Appendix A. Dimensionless Rayleigh-Bénard problem

Let's be the mass, momentum and energy balances:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} = -\nabla p + \nabla \cdot \tau + \rho g \cdot \beta \cdot (T - T_r) \cdot \mathbf{j} \\ \\ \displaystyle \frac{\partial \tau}{\partial t} + \mathbf{v} \cdot \nabla T = \alpha \cdot \nabla^2 T \end{array} \right. \qquad (A.1) $$

where $\mathbf{v}$ is the velocity field, $p$ is the pressure, $T$ the temperature, $\mathbf{D}$ the strain rate tensor (symmetric component of the velocity
<figure><img src="image_20.png" /><figcaption>Fig. 15. Steady state velocity (left) and temperature (right) fields for $R_a = 10^5$ and three different values of the power index $n:n = 0.8$ (top), $n = 1$ (center) and $n = 1.2$ (bottom).</figcaption></figure>

gradient), $\alpha$ the thermal diffusivity ($\alpha = \frac{k}{\rho C_p}$, k being the thermal conductivity – assumed isotropic –, $\rho$ the density and $C_p$ the specific heat), $\tau$ the deviatoric part of the Cauchy's stress tensor, $g$ the gravity's acceleration, $\beta$ the expansion coefficient, $T_r$ a reference temperature and $\mathbf{j}$ the unit vector defining the y-direction along which the gravity applies.

The constitutive equation is assumed given by the power law:

$$ \tau = K \cdot D_{eq}^{n-1} \cdot \mathbf{D} \qquad (A.2) $$

where K is the consistency index, n the power index characterizing the fluid behavior ($n = 1$ for Newtonian fluids and $n \neq 1$ in the non-Newtonian case) and $D_{eq}$ the equivalent deformation related to the second invariant of the rate of strain tensor.

The dimensionless form of these equations is performed by considering the following relations, in which the star superscript refers to the dimensionless variables:

$$ \left\{ \begin{aligned} \mathbf{x} &= \mathbf{x}^* \cdot \mathbf{H} \\ t &= t^* \cdot \frac{H}{(g \cdot \beta H \cdot \Delta T)^{1/2}} \\ \mathbf{v} &= \mathbf{v}^* \cdot (g \cdot \beta H \cdot \Delta T)^{1/2} \\ p &= p^* \cdot (g \cdot \beta H \cdot \Delta T) \cdot \rho \\ \tau &= \tau^* \cdot (g \cdot \beta H \cdot \Delta T) \cdot \rho \\ T &= \theta \cdot \Delta T + T_r \end{aligned} \right. \qquad (A.3) $$

where H is the length of the square cavity in which the flow takes place, $\Delta T$ the temperature difference between the upper and the bottom boundaries having temperatures $T_C$ and $T_H$ ($T_H > T_C$) respectively. In this work we considered $T_r = \frac{T_C+T_H}{2}$.

By introducing relations (A.3) into the balance Eqs. (A.1) it results

$$ \left\{ \begin{aligned} \nabla \cdot \mathbf{v}^* &= 0 \\ \frac{\partial \mathbf{v}^*}{\partial t^*} + \mathbf{v}^* \cdot \nabla \mathbf{v}^* &= -\nabla p^* + \nabla \cdot \mathbf{\tau}^* + \theta \cdot \mathbf{j} \\ \frac{\partial \theta}{\partial t^*} + \mathbf{v}^* \cdot \nabla \theta &= \frac{\alpha}{H \cdot (g \cdot \beta H \cdot \Delta T)^{1/2}} \cdot \nabla^2 \theta \end{aligned} \right. \qquad (A.4) $$
where in the above equations the differential operator $\nabla$ applies with respect to the dimensionless space coordinates $\mathbf{x}^*$.

Now, making use of the expressions of the dimensionless Prandtl ($P_r$) and Rayleigh ($R_a$) numbers, the first defined from the ratio of momentum diffusivity to thermal diffusivity, and the second one related with buoyancy driven flows

$$ \begin{cases} P_r = \frac{K H^{2-n}}{\rho \alpha^2 \bar{n}} \\ R_a = \frac{\rho g \beta H^{2n+1} \Delta T}{K \alpha^n} \end{cases} \qquad (A.5) $$

Eq. (A.4) reduces to:

$$ \begin{cases} \nabla \cdot \mathbf{v}^* = 0 \\ \frac{\partial \mathbf{v}^*}{\partial t} + \mathbf{v}^* \cdot \nabla \mathbf{v}^* = -\nabla p^* + \nabla \cdot \boldsymbol{\tau}^* + \theta \cdot \mathbf{j} \\ \frac{\partial \theta}{\partial t} + \mathbf{v}^* \cdot \nabla \theta = (P_r \cdot R_a)^{-\frac{1}{2}} \nabla^2 \theta \end{cases} \qquad (A.6) $$

When we consider the constitutive Eq. (A.2), its dimensionless form results

$$ (g \cdot \beta H \cdot \Delta T) \cdot \rho \cdot \boldsymbol{\tau}^* = K \cdot \left( \frac{(g \cdot \beta H \cdot \Delta T)^{1/2}}{H} \right)^{n-1} \cdot (D_{eq}^*)^{n-1} \cdot \mathbf{D}^* \qquad (A.7) $$

that can be written as:

$$ \boldsymbol{\tau}^* = P_r^\beta \cdot R_a^{\beta-1} \cdot (D_{eq}^*)^{n-1} \cdot \mathbf{D}^* \qquad (A.8) $$

Thus, by omitting star superscripts, the dimensionless Rayleigh-Bénard model writes

$$ \begin{cases} \nabla \cdot \mathbf{v} = 0 \\ \frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \theta \cdot \mathbf{j} \\ \frac{\partial \theta}{\partial t} + \mathbf{v} \cdot \nabla \theta = (P_r \cdot R_a)^{-\frac{1}{2}} \nabla^2 \theta \\ \boldsymbol{\tau} = P_r^\beta \cdot R_a^{\beta-1} \cdot (D_{eq}^*)^{n-1} \cdot \mathbf{D} \end{cases} \qquad (A.9) $$

## References

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

[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. Part II: Transient simulation using space-time separated representations, J. Non-Newtonian Fluid Mech. 144 (2007) 98-121.

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

[4] A. Ammar, F. Chinesta, P. Diez, A. Huerta, An error estimator for separated representations of highly multidimensional models, Comput. Methods Appl. Mech. Eng. 199 (2010) 1872-1880.

[5] A. Ammar, M. Normandin, F. Chinesta, Solving parametric complex fluids models in rheometric flows, J. Non-Newtonian Fluid Mech. 165 (2010) 1588-1601.

[6] A. Ammar, F. Chinesta, E. Cueto, M. Doblare. Proper Generalized Decomposition of time-multiscale models. Int. J. Numer. Methods Eng. 90(5) (2012) 569-596.

[7] H. Bénard, Les tourbillons cellulaires dans une nappe liquide, Rev. Gén. Sci. Pures Appl. 11 (1900) 1261-1271.

[8] F.H. Busse, Non-linear properties of thermal convection, Rep. Prog. Phys. 41 (1978) 1929-1967.

[9] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Clarendon press, Oxford, 1961.

[10] F. Chinesta, A. Ammar, E. Cueto, Recent advances and new challenges in the use of the Proper Generalized Decomposition for solving multidimensional models, Arch. Comput. Methods Eng. 17/4 (2010) 327-350.

[11] F. Chinesta, A. Ammar, E. Cueto, Proper Generalized Decomposition of multiscale models, Int. J. Numer. Methods Eng. 83/8-9 (2010) 1114-1132.

[12] F. Chinesta, A. Ammar, A. Leygue, R. Keunings, An overview of the Proper Generalized Decomposition with applications in computational rheology, J. Non-Newtonian Fluid Mech. 166 (2011) 578-592.

[13] F. Chinesta, P. Ladeveze, E. Cueto, A short review in model order reduction based on Proper Generalized Decomposition, Arch. Comput. Methods Eng. 18 (2011) 395-404.

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

[15] P. Ladeveze, Non-linear Computational Structural Mechanics, Springer, NY, 1999.

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

[17] P. Ladeveze, L. Chamoin, On the verification of model reduction methods based on the proper generalized decomposition, Comput. Methods Appl. Mech. Eng. 200 (2011) 2032-2047.

[18] 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/2 26494 (2007) 1-14.

[19] A.C. Newell, Th. Passot, J. Lega, Order parameter equations for patterns, Ann. Rev. Fluid Mech. 25 (1993) 399-453.

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

[21] N. Ouertatani, N. Ben Cheikh, B. Ben Beya, T. Lili, Numerical simulation of two-dimensional Rayleigh-Bénard convection in an enclosure, Comptes Rendus Mécanique 336/5 (2008) 464-470.

[22] H.M. Park, K.S. Park, Rayleigh-Bénard convection of viscoelastic fluids in arbitrary finite domains, Int. J. Heat Mass Transfer 47/10-11 (2004) 2251-2259.

[23] H.M. Park, Y.M. Heo, Rayleigh-Bénard convection in two-dimensional arbitrary finite domains, Int. J. Thermal Sci. 45/7 (2006) 697-705.

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

[25] Lord Rayleigh, On the convective currents in a horizontal layer of fluid when the higher temperature is on the underside, Philos. Mag. 32 (1916) 529-546.

[26] S. Vedantam, M.T. Dhotre, J.B. Joshi, Three-dimensional CFD simulation of Rayleigh-Bénard convection for low Prandtl number fluids, Chem. Eng. Res. Des. 84/1 (2006) 29-37.

[27] A. Vikhansky, Thermal convection of a viscoplastic liquid with high Rayleigh and Bingham numbers, Phys. Fluids 21 (2009) 103103.

[28] J. Zhang, D. Vola, I.A. Frigaard, Yield stress effects on Rayleigh-Bénard convection, J. Fluid Mech. 566 (2006) 389-419.