# Résolution numérique des équations de transfert par une méthode de réduction a priori

Nicolas Verdon, Cyrille Allery, Aziz Hamdouni

► To cite this version:

Nicolas Verdon, Cyrille Allery, Aziz Hamdouni. Résolution numérique des équations de transfert par une méthode de réduction a priori. CFM 2007 - 18ème Congrès Français de Mécanique, Aug 2007, Grenoble, France. hal-03361697

HAL Id: hal-03361697

https://hal.science/hal-03361697v1

Submitted on 1 Oct 2021

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Résolution numérique des équations de transfert par une méthode de réduction a priori

Nicolas Verdon, Cyrille Allery & Aziz Hamdouni

Université de La Rochelle
Laboratoire d'Etude des Phénomènes de Transfert Appliqués au Bâtiment
Avenue Michel Crépeau F-17042 La Rochelle cedex 1
nverdon@univ-lr.fr

## Résumé :

Ce papier décrit l'application de la méthode de réduction a priori (APR) à plusieurs cas académiques. Cet algorithme, qui consiste a déterminer sans échantillonage préalable une base d'un écoulement, utilise la décomposition de Karhunen-Loève et des sous-espaces de Krylov, et a été appliquée dans le cas de l'équation de convection-diffusion 2D, ainsi que dans le cas de la cavité entraînée. Couplée à un système dynamique, cette méthode présente l'avantage d'être extrémement rapide et permet un réajustement dynamique de la base lors de modifications de paramètres contrôlant l'écoulement.

## Abstract :

This paper describes the application of the A Priori Reduction method (APR) in several academic test cases. This algorithm, which consists in determining a basis of a flow without prerequisite knowledge, uses the Karhunen-Loève Decomposition and Krylov subspaces, and has been applied to the convection-diffusion equation and also for the driven cavity configuration. Coupled with a dynamical system, this method has the important advantage to be computationally very fast and allows a dynamical update of the basis when control parameters are changing.

## Mots-clefs :

### Modèles d'Ordre Réduit ; Navier-Sokes ; Système Dynamique

## 1 Introduction

La résolution numérique des équations de transfert a une importance grandissante du point de vue fondamental mais aussi à un niveau industriel. Cependant, certaines méthodes requièrent à l'heure actuel un effort de calcul trop important pour prétendre à des applications industrielles. Par exemple, le contrôle actif nécessite des temps de calcul très courts pour permettre des adaptations du système mis en jeu en temps réel. Dans ce contexte on trouve bon nombre de méthodes de réduction de modèles. Dans le cas de la mécanique des fluides, qui présente des applications variées, la méthode intensivement utilisée est la décomposition orthogonale aux valeurs propres (POD) qui possède l'avantage de construire une base optimale du point de vue énergétique. Plus de propriétés et de détails théoriques sur cette méthode peuvent être trouvés dans Holmes et al. (1996). Associée à un système dynamique obtenu par une méthode Galerkin, cette approche permet de représenter l'évolution d'un écoulement sur des temps bien plus longs que ceux nécessaires à la construction de la dite base. Cependant, il est connu que le système dynamique n'est valable que pour un certain temps, voir Sirisup et Karniadakis (2004), avant que la solution ne diverge. De même, un changement dans les paramètres de contrôle de l'écoulement, à savoir le nombre de Reynolds ou la modification d'une condition aux limites, n'est pas bien représenté par le système dynamique. La base calculée n'est effectivement va-
lable que pour une gamme étroite de ces paramètres.

Pour s'affranchir des problèmes liés à la POD, nous proposons une méthode itérative de construction de la base. Cette approche présente l'avantage de *réajuster* la base lorsqu'elle n'est plus valable (temps trop longs, changement de paramètres de contrôle, etc...).

Le papier s'organise de la façon suivante : dans un premier temps la méthode est décrite pour une équation de transfert quelconque et dans un deuxième temps les résultats obtenus pour la convection-diffusion ainsi que pour les équations de Navier-Stokes dans la configuration de la cavité entraînée seront présentés.

## 2 Présentation de la méthode APR

La méthode APR étant itérative, il convient d'initialiser le processus de résolution, c'est-à-dire ici de s'imposer une base initiale. Comme il a déjà été mis en évidence par Ryckelynck (2006), le choix de cette base n'est pas un élément très important pour la qualité de la solution, mais détermine la rapidité et le nombre d'itérations que l'algorithme devra effectuer. En effet, choisir une base plus physique requiert moins d'itérations de calcul APR pour faire converger la solution. On appelle cette base initiale $\psi^0$. La méthode étant expliquée en détails dans Verdon et al. (2006), seules ses grandes lignes seront expliquées dans les paragraphes suivants.

### 2.1 Construction du système réduit

On considère donc le problème physique<fnref n="1" /> que l'on doit résoudre, que l'on note :

$$ \mathcal{L}(\mathbf{u}, p) = b \qquad (1) $$

La discrétisation de l'équation (1) conduit, après linéarisation si nécessaire, au système d'équations linéaires suivant :

$$ A\zeta = B \qquad (2) $$

où A est la matrice issue de la discrétisation de l'équation (1), de taille N, ζ est le vecteur inconnues et enfin B le second membre, incluant les conditions aux limites et les données de la solution calculée au pas de temps antérieur.

Notons maintenant n le nombre de vecteurs de base de $\psi$ ($n \ll N$). Le système réduit associé au modèle complet (2) est obtenu en écrivant l'inconnue $\zeta$ sur la base $\psi$

$$ \zeta(\mathbf{x}, t) = \sum_{i=1}^{n} \psi_i(\mathbf{x})a_i(t) \qquad (3) $$

et en projetant le système obtenu sur cette même base. Le système réduit dont les nouvelles inconnues sont les coefficients temporels $a_i(t)$ s'écrit de la façon suivante :

$$ \mathcal{F}(a_i(t)) = 0 \quad \forall t \in T \quad \forall i = 1, \dots, n \qquad (4) $$

où T est l'intervalle de temps sur lequel on désire construire la base de la solution. Il faut donc résoudre ce système réduit pour tous les pas de temps de l'intervalle T. La solution complète est alors reconstruite en utilisant l'expression de sa décomposition (3) sur la base $\psi$. La qualité de la solution est alors évaluée sur tout l'intervalle en considérant le résidu du modèle complet $r = A\zeta - B$. Si la solution n'est pas satisfaisante, la base est alors adaptée comme expliqué dans le paragraphe suivant.

<footnote n="1">¹équation de convection-diffusion, équations de Navier-Stokes, etc...</footnote>
## 2.2 Adaptation et enrichissement de la base

Lors de l'itération courante, on a résolu l'équation réduite sur tout l'intervalle de temps, c'est-à-dire que l'on a déterminé les coefficients temporels $a(t)$, $\forall t \in T$. Pour la phase d'amélioration de la base, on réalise une décomposition de Karhunen-Loève sur cet ensemble. On est donc ammené à résoudre le problème aux valeurs propres suivant :

$$C\Phi = \lambda\Phi \quad (5)$$

où C est la matrice de corrélation définie de la façon suivante :

$$C_{ij} = \sum_{l=1}^{N_t} a_i(t_l) a_j(t_l) \quad (6)$$

avec $N_t$ le nombre de pas de temps. On obtient ainsi $n$ vecteurs propres correspondant aux $n$ valeurs propres solutions de cette équation. On considère alors les $\eta$ plus grandes valeurs propres, définies telles que $\lambda_{\eta} > \mu_{sel} \cdot \lambda_1^2$ afin de ne garder que les modes propres dont la contribution dans l'écoulement est significative de la dynamique observée sur l'intervalle de temps. La matrice d'amélioration $V$ est alors formée de la façon suivante : ses vecteurs colonnes sont les $\eta$ vecteurs propres conservés, c'est-à-dire $V = [\Phi_1 \ \Phi_2 \ \dots \ \Phi_{\eta}]$. La nouvelle base est alors améliorée à l'aide de cette matrice $V$ :

$$\tilde{\psi} = \psi \cdot V \quad (7)$$

Ensuite, il faut également tenir compte de l'erreur que l'on a comise à l'itération $k$, c'est pourquoi on conserve le premier résidu du modèle complet $r = A\zeta - B$ qui ne satisfait pas le critère de convergence et on construit l'espace de Krylov<fnref n="3" /> d'ordre $m$ associé à ce résidu :

$$K_m = \{r, Ar, A^2 r, \dots, A^{m-1} r\} \quad (8)$$

A la base améliorée $\tilde{\psi}$, on rajoute les vecteurs du sous-espace de Krylov que l'on vient de calculer, c'est ce qu'on appelle la phase d'enrichissement de la base, c'est-à-dire le moment où l'on prend en considération les erreurs de l'itération précédente. Dans la suite, on prend $m = 1$. La base ainsi mise à jour, c'est-à-dire améliorée par la décomposition de Karhunen-Loève et enrichie par les vecteurs de Krylov, est utilisée pour l'itération suivante. Si elle permet d'obtenir une solution suffisamment précise, elle sera conservée par la suite pour la construction du système dynamique. Dans le cas contraire, elle sera mise à jour jusqu'à ce que la qualité de la solution calculée soit jugée satisfaisante.

## 3 Système dynamique

Une fois que l'on a obtenu une base qui représente correctement la solution du problème considéré sur l'intervalle de temps étudié, on espère simuler la solution sur un intervalle de temps beaucoup plus long que celui nécessaire au calcul de la base. Pour ce faire, on construit un système dynamique en effectuant une projection de Galerkin de l'équation sur la base de l'écoulement.

En utilisant l'expression (3) de la décomposition de la vitesse sur la base $\psi$, et par projection, on obtient le système dynamique suivant :

$$\frac{da_i}{dt} = \sum_{j=1}^{n} B_{ij} a_j + \sum_{j=1}^{n} \sum_{k=1}^{n} C_{ijk} a_j a_k + D_i \quad \text{pour } i = 1, \dots, n \quad (9)$$

<footnote n="2">²où les valeurs propres sont rangées dans l'ordre décroissant et $\mu_{sel}$ est le coefficient de sélection que l'on se fixe</footnote>

<footnote n="3">³Pour plus de détails sur les méthodes de Krylov, le lecteur pourra se référer à Saad (1996)</footnote>
où l'on a pour l'équation de convection-diffusion :

$$B_{ij} = (\psi_n, -\alpha\psi_m + \nu\Delta\psi_m) \qquad C_{ijk} = 0 \quad \text{et} \quad D_i = 0 \qquad (10)$$

et pour les équations de Navier-Stokes :

$$B_{ij} = (\psi_n, \nu\Delta\psi_m) \qquad C_{ijk} = -(\psi_n, \nabla\psi_m) \quad \text{et} \quad D_i = \int_{\partial\Omega} P \cdot \psi_n \, n \, dS \qquad (11)$$

La résolution du système d'équations différentielles (9) par un algorithme de Runge-Kutta est très rapide en comparaison avec la résolution des équations pendant la précédente phase. La solution complète est alors calculée à l'aide de (3) sur des temps bien plus long que le temps nécessaire pour la méthode APR.

# 4 Résultats

Cette partie présente les deux cas tests que nous avons étudiés, à savoir l'équation de convection diffusion et les équations de Navier-Stokes pour le cas d'une cavité entraînée.

## 4.1 Cas 1 : équation de convection-diffusion

L'équation de convection-diffusion a été étudiée dans le domaine carré $0 \le x, y \le 2$ avec la connaissance de la solution analytique donnée dans Tian et Ge (2007)<fnref n="4" /> :

$$\phi(x, y, t) = \frac{1}{4t+1} \exp \left[ -\frac{(x - \alpha t - 0.5)^2}{\Gamma(4t+1)} - \frac{(y - \alpha t - 0.5)^2}{\Gamma(4t+1)} \right] \qquad (12)$$

Les conditions aux limites de type Dirichlet et les conditions initiales ont été obtenues à partir de cette solution analytique. Par ailleurs, la base initiale de notre algorithme correspond aux conditions initiales du problème.

Plusieurs expérimentations numériques ont été effectuées sur cet exemple.

Dans un premier temps, nous avons testé la faculté de réadaptation de la base APR lorsque le système dynamique diverge, et dans une deuxième approche, on a étudié l'adaptation rapide de la base pour un changement de paramètre de contrôle, à savoir lorsque $\Gamma$ varie.

Dans tout ce qui suit, l'intervalle de calcul est noté $T$ et correspond à une durée $T = 0.0125s$. Dans un premier temps les calculs on été menés sur un intervalle $T$ à partir de $t = 0s$ afin de déterminer une première base APR. A partir de cette base, comportant 5 vecteurs, on construit un système dynamique que l'on résout sur $10T$. On constate alors qu'à partir de $t = 6T$, la base APR ne représente plus correctement la dynamique du système et l'erreur calculée par rapport à la solution analytique dépasse les 5%<fnref n="5" /> pour atteindre 22% à $t = 10T$. On effectue alors un nouveau calcul APR sur un intervalle $T$ à partir de $t = 0.075s$ et avec la nouvelle base obtenue, on construit le nouveau système dynamique correspondant qui est résolu jusqu'à $t = 10T$ c'est-à-dire l'intervalle que l'on voulait étudier. L'erreur à la fin du calcul est alors aux alentours de 5% (voir fig. 1), ce qui est fortement amélioré par rapport aux performances de la première base et le temps de calcul nécessaire a été divisé par deux par rapport au modèle complet.

Dans une deuxième approche, on a testé la sensibilité de la base à un changement de paramètre, ici $\Gamma$. On repart de la première base calculée sur l'intervalle $T$ et on construit le système dynamique correspondant, mais cette fois-ci avec $\Gamma = 0.02$ alors que l'itération d'APR a été effectuée pour $\Gamma = 0.01$. On se rend alors compte que même sur l'intervalle $T$ où a été calculée la

<footnote n="4">⁴avec $\alpha = 0.8$ et $\Gamma = 0.01$</footnote>

<footnote n="5">⁵on juge arbitrairement que la base n'est plus valable lorsque l'erreur dépasse 5%</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1 – Résultats du cas 1 : isolignes de $\phi$ pour $t = 0.0125s$ et $t = 0.125s$</figcaption></figure>

base, le système dynamique n’est pas capable de reproduire fidèlement la solution puisque l’on trouve une erreur égale à 22.16% pour $t = T$. L’idée est alors de relancer la méthode APR, en considérant comme base initiale la première base calculée. En seulement 3 itérations APR, on arrive à une base, comprenant 4 vecteurs, qui représente correctement la solution, c’est-à-dire que le système dynamique formé à partir de cette base ne donne plus qu’une erreur d’environ 1.9% à $t = 0.0125s$.

Par ailleurs, la modification de la base n’a nécessité que 6s de temps de calcul alors que la construction de la première base en a requis environ 20s, ce qui présente un intérêt en vue d’applications en contrôle actif.

## 4.2 Cas 2 : équations de Navier-Stokes

Le cas étudié ici est la cavité entraînée unitaire ($0 \le x, y \le 1$) représentée sur la figure 2 pour un nombre de Reynolds $Re = 1$ avec, initialement, le fluide au repos. Ce cas a été traité

<figure><img src="image_3.png" /><figcaption>FIG. 2 – Configuration de la cavité entraînée 2D</figcaption></figure>

sans le système dynamique mais nous a tout de même permis de valider la méthode APR pour les équations de Navier-Stokes. On constate que la méthode APR permet de retrouver la solution stationnaire (fig. 3) pour un nombre de Reynolds valant $Re = 1$ en environ deux fois moins de temps que le modèle complet, c’est-à-dire 413.3s de calcul contre 781,8s pour le modèle complet, et une erreur relativement faible, comme présenté sur la figure 3(c).
<figure><img src="image_4.png" /><figcaption>FIG. 3 – Résultats du cas 2 : isolignes de la vitesse pour $\mathcal{R}_e = 1$</figcaption></figure>

Cette étude est une étude préliminaire à un cas instationnaire qui sera présenté lors de la communication orale.

## 5 Conclusions

Dans ce travail, la méthode APR a été appliquée avec succès à deux configurations : l'équation de convection-diffusion instationnaire et les équations de Navier-Stokes dans le cas de la cavité entraînée à $\mathcal{R}_e = 1$. Les performances qui sont ici mises en avant sont d'une part, le gain de temps comparé à un calcul avec un modèle complet, puis une facilité et une rapidité d'adaptation de la base APR, d'une part lorsque le système dynamique diverge, et d'autre part lorsqu'un paramètre de contrôle change. La prochaine étape du travail consiste donc à appliquer cette méthode à des écoulements turbulents.

## Références

* Holmes P., Lumley J.L., and Berkooz G., Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge University Press, Cambridge, ISBN 0 521 55142 0, 1996.

* Saad Y., Iterative Methods for Sparse Linear Systems, 2nd Edition, 2000.

* Ryckelynck D., A priori hyperreduction method : an adaptive approach, *Journal of Computational Physics*, Volume 202, Issue 1, pp. 346-366, 1er janvier 2005.

* Verdon N., Allery C., Ryckelynck D., Hamdouni A., An adaptive ROM approach for solving transfer equations, *Revue européenne de mécanique numérique*, Volume 15, N. 5, 2006

* Sirisup S., Karniadakis G. E., A spectral viscosity method for correcting the long term behavior of POD models, *Journal of Computational Physics*, Volume 194, p. 92-116, 2004

* Tian Z. F., Ge Y. B. 2007 A fourth-order compact ADI method for solving two-dimensional unsteady convection-diffusion problems, *Journal of Computational and Applied Mathematics* Volume 198, pp. 268-286