Documentation related to velocity correction (flow rate prescription at the...

Documentation related to velocity correction (flow rate prescription at the oulet of the box-domain ...)

Docs/velocity_correction.tex

+\documentclass[a4paper,10pt]{article}
+\usepackage[utf8x]{inputenc}
+\usepackage{amsmath}
+% \usepackage{showkeys} % affiche les labels et references
+\usepackage[T1]{fontenc}
+\usepackage{multirow}
+\usepackage{amsthm,amssymb}
+\setlength{\parskip}{12pt} %taille du saut entre deux paragraphes
+\usepackage{graphicx}
+\graphicspath{{./images/}}
+\usepackage{pdfsync} % permet de cliquer sur le pdf pour avoir l endroit ou il y a le code latex
+ \usepackage{cancel}
+\usepackage{xcolor}
+\usepackage[left=3cm,right=3cm,top=2.5cm,bottom=2.5cm]{geometry}
+\overfullrule 5pt %Pour d\'esigner les listes couples d'un *full carre noir la ou l on sort
+
+\usepackage[pdftex]{hyperref}
+\hypersetup{
+  pdftitle={Flow past 3D sphere},
+  pdfauthor={Mimeau Chloe},
+  pdfcreator={PARMES},
+  pdfkeywords={sphere, potential flow, velocity correction, flow rate},
+  pdfsubject={Flow past 3D sphere}
+% pagebackref=false,
+%colorlinks=false,citecolor=red,
+%linkcolor=blue,
+%plainpages=false,
+%breaklinks=true,
+% pdfpagelabels=true
+}
+
+
+\theoremstyle{plain}
+\newtheorem{Theoreme}{Th\'eor\`eme}[section]
+%opening
+\title{Correction de vitesse pour la modélisation d'un écoulement incompressible dans une boîte}
+\author{}
+\date{}
+
+\begin{document}
+
+\maketitle
+Lors de la résolution de l'équation de Poisson $\Delta \mathbf{u} = - \nabla \boldsymbol{\omega} $, le mode 0 des transformées de Fourier rapides de la vitesse est imposé à 0, ce qui entraîne une modification des valeurs moyennes des composantes de la vitesse et donc des composantes de la vitesse elles-mêmes lorsque la transformée de Fourier inverse est appliquée. Une correction de la vitesse est donc nécessaire après la résolution de l'équation de Poisson, à chaque itération.
+
+On introduit la notation suivante pour exprimer la vitesse corrigée :
+%
+\begin{align}
+\bf{u}= \widetilde{u} + \overline{u}\label{ux} 
+\end{align}
+
+où $\widetilde{\bf{u}} $ désigne la vitesse en sortie de FFT et où $\overline{\bf{u}} $ dénote la moyenne en espace de la vitesse, que l'on cherche ici à déterminer. L'écoulement s'effectue dans la direction $x$ et la vitesse est initialisée en posant ${\bf{u}}(t=0)=(u_x, u_y, u_z)=(u_{\infty}, 0, 0)$.
+
+
+$\overline{u_x}$ est un flux constant en espace, par conséquent $\dfrac{\partial \overline{u_x}}{\partial x} = \dfrac{\partial \overline{u_x}}{\partial y} = \dfrac{\partial \overline{u_x}}{\partial z} = 0$. Ainsi le champ de vorticité moyen en espace est donné en 2D par :
+\begin{align}\label{vortMoy}
+\overline{\omega} & = \dfrac{\partial \overline{u_y}}{\partial x}
+\end{align}
+et en 3D par :
+\begin{align}
+\overline{\omega_x} & = \dfrac{\partial \overline{u_z}}{\partial y} - \dfrac{\partial \overline{u_y}}{\partial z} \\
+\overline{\omega_y} & = -\dfrac{\partial \overline{u_z}}{\partial x} \label{wy} \\
+\overline{\omega_z} & = -\dfrac{\partial \overline{u_y}}{\partial x} \label{wz}
+\end{align}
+
+La correction de la vitesse se fait par réajustement du débit à l'entrée du domaine de calcul. Cette correction est développée ci-dessous composante par composante tout d'abord en 2D puis sera étendue au cas 3D.
+
+%%%%%%%%%%% 2D %%%%%%%%%%%%%%
+{\bf{2D}}
+
+\underline{Correction de $u_x$}:\\
+
+On note $B$ le coté correspondant au bord d'entrée du domaine de calcul. En partant de la décomposition \ref{ux} on a alors :
+\begin{align}
+\underbrace{\int\limits_B  u_x \ dy}_{\text{débit souhaité}} = \underbrace{\int\limits_B \widetilde{u_x} \ dy}_{\text{débit calculé}} + \int\limits_B \overline{u_x} \ dy \\
+~ 
+\Longleftrightarrow \underbrace{\int\limits_B  u_x \ dy}_{\text{débit souhaité}} = \underbrace{\int\limits_B \widetilde{u_x} \ dy}_{\text{débit calculé}} + \ L_y \overline{u_x} \\ \label{debSouhaiteUx}
+~ 
+\Longleftrightarrow \ \ \overline{u_x} = \dfrac{\text{débit souhaité}}{L_y} - \dfrac{\text{débit calculé}}{L_y}
+\end{align}
+
+où $L_y$ dénote la longueur du coté $B$. Or le débit que l'on souhaite imposer à l'entrée du domaine est égal à $u_{\infty} L_y$ et donc :
+\begin{equation}
+\boxed{u_x = \widetilde{u_x} + u_{\infty} - \dfrac{\text{débit calculé}}{L_y}}
+\end{equation}
+
+
+
+\underline{Correction de $u_y$}:\\
+
+$\overline{u_y}$ est un flux invariant en $y$. \\
+En effet, la vitesse est a divergence nulle, donc la moyenne en espace de la vitesse l'est  également et ainsi $\partial_x \overline{u_x} + \partial_y \overline{u_y} = 0$. Or $\overline{u_x}$ est un flux constant en espace, par conséquent on a bien $\partial_y \overline{u_y} = 0$ (invariance de $\overline{u_y}$ en $y$).\\
+
+$\overline{u_y}$ est donc un flux variant en $x$ et d'après la relation \ref{vortMoy} on a $\overline{u_y}= \overline{\omega} x + c$, avec $c$ une constante. Ainsi :
+\begin{align}
+\underbrace{\int\limits_B  u_y \ dy}_{\text{débit souhaité}} = \underbrace{\int\limits_B \widetilde{u_y} \ dy}_{\text{débit calculé}} + \int\limits_B \overline{\omega} x \ dy + \int\limits_B c \ dy 
+\end{align}
+
+Or ici le débit souhaité est égal a 0. Donc :
+\begin{align}
+0 = \underbrace{\int\limits_B \widetilde{u_y} \ dy}_{\text{débit calculé}} + \overline{\omega} x_0 L_y + L_y  c  \\
+~ 
+\Longleftrightarrow c = -\dfrac{\text{débit calculé}}{L_y} - \overline{\omega} x_0
+\end{align}
+
+et donc :
+\begin{equation}
+\boxed{u_y = \widetilde{u_y} + \overline{\omega} x - \overline{\omega} x_0 -  \dfrac{\text{débit calculé}}{L_y}}
+\end{equation}
+
+où $x_0$ dénote la coordonnée dans la direction $x$ du bord d'entrée du domaine de calcul.
+
+%%%%%%%%%%% 3D %%%%%%%%%%%%%%
+{\bf{3D}}
+
+\underline{Correction de $u_x$}:\\
+
+On note $S$ la surface correspondant au bord d'entrée du domaine de calcul. En partant de la décomposition \ref{ux} on a alors :
+\begin{align}
+\underbrace{\iint\limits_S  u_x \ dydz}_{\text{débit souhaité}} = \underbrace{\iint\limits_S \widetilde{u_x} \ dydz}_{\text{débit calculé}} + \iint\limits_S \overline{u_x} \ dydz \\
+~ 
+\Longleftrightarrow \ \ \overline{u_x} = \dfrac{\text{débit souhaité}}{L_yL_z} - \dfrac{\text{débit calculé}}{L_yL_z} \label{debSouhaiteUx3D}
+\end{align}
+
+et donc :
+\begin{equation}
+\boxed{u_x = \widetilde{u_x} + u_{\infty} - \dfrac{\text{débit calculé}}{L_yL_z}}
+\end{equation}
+
+\underline{Correction de $u_y$}:\\
+
+D'après la relation \ref{wz} on a $\overline{u_y} = \overline{\omega_z}x +c$ avec $c$ une constante. D'après \ref{ux} on obtient donc que $u_y = \widetilde{u_y} + \overline{\omega_z} x +c$. Ainsi,
+\begin{align}
+\iint\limits_S u_y \ dydz = \iint\limits_S \widetilde{u_y} \ dydz + \iint\limits_S \overline{\omega_z}x \ dydz + \iint\limits_S c \ dydz \\
+~ 
+\Longleftrightarrow \ \underbrace{\iint\limits_S u_y \ dydz}_{\text{débit souhaité}} = \underbrace{\iint\limits_S \widetilde{u_y} \ dydz}_{\text{débit calculé}} + \ \overline{\omega_z}\ x_0L_yLz + cL_yL_z
+\end{align}
+
+où $x_0$ désigne la coordonnée dans la direction $x$ de la surface d'entrée $S$ du domaine de calcul.
+
+Donc $\ \ \ c= \dfrac{\text{débit souhaité}}{L_yL_z} - \dfrac{\text{débit calculé}}{L_yL_z} - \overline{\omega_z}\ x_0$
+
+Or le débit souhaité est égal a 0, par conséquent :
+\begin{equation}
+\boxed{u_y = \widetilde{u_y} + \overline{\omega_z}x - \dfrac{\text{débit calculé}}{L_yL_z} - \overline{\omega_z}\ x_0}
+\end{equation}
+
+\underline{Correction de $u_z$}:\\
+
+Sur le même principe on a d'après la relation \ref{wy} que $\overline{u_z} = -\overline{\omega_y}x +c$ avec $c$ une constante. D'après \ref{ux} on obtient donc que $u_z = \widetilde{u_z} + \overline{\omega_y} x +c$. Ainsi en procédant de même que pour $u_y$ on obtient :
+\begin{equation}
+\boxed{u_z = \widetilde{u_z} - \overline{\omega_y}x - \dfrac{\iint\limits_S \widetilde{u_z} \ dydz}{L_yL_z} + \overline{\omega_y}\ x_0}
+\end{equation}
+
+{\bf{Remarque :}}
+
+Dans le cadre d'un écoulement autour d'un cylindre (2D) ou d'une sphère (3D), on peut également choisir l'écoulement potentiel comme condition initiale de l'écoulement.
+
+Par décomposition de Helmholtz on a:
+\begin{equation}
+\mathbf{u} = \nabla \times \psi + \nabla \phi
+\end{equation}
+où la fonction courant $\psi$ et le potentiel scalaire $\phi$ vérifient les propriétés suivantes :
+\begin{equation}
+\text{div} \ (\psi) = 0, \ \ \ \
+\nabla \times \nabla \phi= 0
+\end{equation}
+Donc si l'écoulement est potentiel, c'est à dire irrotationnel, alors :
+\begin{equation}
+\mathbf{u} = \nabla \phi
+\end{equation}
+et \underline{en 2D}, la fonction courant et le potentiel scalaire sont liés par la relation suivante :
+\begin{equation}
+\dfrac{\partial \phi}{\partial x} = \dfrac{\partial \psi}{\partial y}, \ \ \
+\dfrac{\partial \phi}{\partial y} = -\dfrac{\partial \psi}{\partial x} 
+\end{equation}
+Ainsi $u_x = \dfrac{\partial \psi}{\partial y} $ et $u_y= -\dfrac{\partial \psi}{\partial x}$.
+
+La fonction courant associée à un écoulement potentiel 2D autour d'un cylindre de rayon $R$ est exprimée de la façon suivante :
+\begin{equation}
+\psi = u_{\infty} y \left ( 1 + \dfrac{R^2}{(x-x_C)^2 + (y-y_C)^2} \right )
+\end{equation}
+
+avec $(x_C, y_C)$ les coordonnées du centre du cylindre dans le repère cartésien du domaine de calcul.
+
+Ainsi l'expression du débit souhaité, intervenant dans l'équation \ref{debSouhaiteUx}, est désormais donnée par :
+\begin{equation}
+\underbrace{\int\limits_B  u_x \ dy}_{\text{débit souhaité}} = \int\limits_B  \dfrac{\partial \psi}{\partial y} \ dy = [\psi]_B = \left[ u_{\infty} y \left (1-\dfrac{R^2}{(x_0-x_C)^2-(y-y_C)^2}\right ) \right ]_B
+\end{equation}
+
+\underline{En 3D}, il est pratique d'exprimer le potentiel d'un écoulement autour d'une sphère de rayon $R$ dans un système de coordonnées sphériques $(r, \theta, \varphi)$, dans lequel l'angle dans la direction azimutale $\theta$ est mesuré à partir de la direction de l'écoulement. Dans ce système de coordonnées, le potentiel est indépendant de l'angle méridional $\varphi$ et est donné par :
+\begin{equation}
+\phi = -u_{\infty} r \ \text{cos} \theta\left ( 1 + \dfrac{R^3}{2r^3} \right )
+\end{equation}
+%
+Ainsi, les composantes de la vitesse dans les directions $r$ et $\theta$ sont les suivantes:
+\begin{align}
+u_r = \dfrac{\partial \phi}{\partial r} = -u_{\infty} \text{cos} \theta \left ( 1 - \dfrac{R^3}{r^3} \right ) \\
+u_{\theta} = \dfrac{1}{r}\dfrac{\partial \phi}{\partial \theta} = u_{\infty} \text{sin} \theta \left ( 1 + \dfrac{R^3}{2 r^3} \right )
+\end{align}
+%
+Dans un système de coordonnées cartésiennes, on obtient donc :
+\begin{align}
+u_x &= \dfrac{-u_{\infty} x}{r} \left ( 1 - \dfrac{R^3}{r^3} \right ) \\
+u_y &= \dfrac{u_{\infty} y}{r}\left ( 1 + \dfrac{R^3}{2 r^3} \right ) \\
+u_z &= 0
+\end{align}
+
+où $r=\sqrt{(x-x_S)^2 + (y-y_S)^2 + (z-z_S)^2}$, avec $(x_S, y_S, z_S)$ les coordonnées du centre de la sphère dans le repère cartésien du domaine de calcul.
+
+Ainsi l'expression du débit souhaité, intervenant dans l'équation \ref{debSouhaiteUx3D}, est désormais donnée par:
+\begin{equation}
+\underbrace{\iint\limits_S  u_x \ dydz}_{\text{débit souhaité}} = \iint\limits_S \ \dfrac{-u_{\infty} x}{r} \left ( 1 - \dfrac{R^3}{r^3} \right )  \ dydz = -u_{\infty} x_0 \iint\limits_S \ \dfrac{1}{r} \left ( 1 - \dfrac{R^3}{r^3} \right )  \ dydz 
+\end{equation}
+\end{document}