À un instant `t`, l'état du fluide est connu en tout point de l'espace. Les équations euleriennes vont permettre de calculer l'état du fluide à l'instant `t"+"dt`.
Les vitesses du courant du fluide forment un champ de vecteurs. Les lignes de ce champ sont appelées lignes de courants. C'est l'approche usuelle, selon le point de vue d'Euler. Remarquez que les trajectoires des particules du fluide ne sont en général pas égales aux lignes de courants sauf si ces lignes de courants sont constantes.
Un point dans cet espace-temps à deux dimensions d'espace est caractérisé par trois composantes réels `(x,y,t)`. Les composantes `x` et `y` représentent la position dans l'espace, tandis que la composante `t` représente la position dans le temps.
On définie le vecteur position `vec r = ((x),(y))`
Cela produit l'élément différentiel de déplacement suivant : `dvec r = ((dx),(dy))`
On définie le vecteur position-temps : `vec X = ((x),(y),(t))`
Cela produit l'élément différentiel suivant : `dvec X = ((dx),(dy),(dt))`
On commence par étudier les fluides incompressibles c'est à dire de densité constante, et de viscosité constante. Ainsi le fluide possède deux caractéristiques importantes que sont sa densité `rho` et sa viscosité dynamique `mu`.
L'état du fluide est décrit par un champ de pression `P` et un champ de vitesse `vecv`, et si on ne précise pas les coordonnées du points où est mesuré le champ, par défaut ce sera à la position `(x,y,t)`. Ainsi nous avons :
`P = P(x,y,t)`
`vecv=vecv(x,y,t)``vecv=((v_x),(v_y))=((v_x(x,y,t)),(v_y(x,y,t)))`
La résolution numérique se fait en subdivisant l'espaces et le temps en une grille de petites mailles. On considère que dans chaque maille l'état du fluide, pression et vitesse, varient très peu. On choisie une taille de maille `(dx,dy,dt)`. On choisie un modèle particulièrement simple de maillage où il n'y a pas vraiment de condition aux limites, qui est topologiquement équivalent à la surface d'un tore. Le maillage est un tableau de `N"×"N` mailles, et on raccorde les bords du tableau pour former un tore, faisant qu'une maille `(x,y)` est toujours entourée de mailles comme suit, en effectuant le calcul des coordonnées selon l'axe des `x` modulo `Ndx` et selon l'axe des `y` modulo `Ndy` :
`( ( "("x"-"dx","y"+"dy")" , "("x","y"+"dy")", "("x"+"dx","y"+"dy")" ) , ( "("x"-"dx","y")" , "("x","y")" , "("x"+"dx","y")" ) , ( "("x"-"dx","y"-"dy")" , "("x","y"-"dy")" , "("x"+"dx","y"-"dy")" ) )`
On a choisie un axe des `x` horizontal allant vers la droite, appelé abscisse. Et on a choisie un axe des `y` vertical allant vers le haut, appelé ordonnée. C'est une convention courante en mathématique. Cependant, il faut noter que sur un ordinateur, les commandes d'affichages standards différent en choississant une abscisse verticale n'allant pas vers le haut mais vers le bas.
Les vecteurs sont représentés verticalement. Leur transposés sont représentés horizontalement. La transposé d'un vecteur `vecv` se note `vecv^TT`. Exemples :
`vec v = ((v_x),(v_y))`
`vecv^TT = ((v_x, v_y))`
`(vecv^TT)^TT= vecv`
On introduit cette distinction pour pouvoir définir un produit matriciel s'appliquant aussi bien aux matrices qu'aux vecteurs. Les vecteurs se généralisent en des matrices et l'opérateur de transposé se généralise aux matrices. Voyez par exemples :
`M =((a,b),(c,d))`
`M^TT = ((a,c),(b,d))`
`(M^TT)^TT= M`
Le produit matriciel entre deux matrices calcule la sommes des produits termes à termes d'une ligne de la première matrice avec une colonne de la seconde matrice et la place dans une matrice résultat au même numéros de ligne et colonne utilisés. Aussi il n'est définie que si le nombre de termes dans chaque ligne de la première matrice est égale au nombre de termes dans chaque colonne de la seconde matrice, autrement dit, si le nombre de colonnes de la première matrice est égale au nombre de lignes de la seconde matrice. Voici quelques exemples :
`(x,y)((a),(b)) = xa+yb`
`((x),(y))(a,b) = ((xa, xb),(ya,yb))`
`((a,b),(c,d))((x),(y)) = ((ax"+"by),(cx"+"dy))`
`(x,y)((a,b),(c,d) )= ((xa"+"yc, xb"+"yd))`
`((a,b),(c,d))((p,q),(r,s)) = ((ap"+"br,aq"+"bs),(cp"+"dr,cq"+"ds))`
Le produit scalaire de deux vecteurs correspond à la somme des produits termes à termes des deux vecteurs. Aussi il n'est définie que si les vecteurs sont de même dimension. Voici des exemples :
`((x),(y))"⋅"((a),(b)) = xa + yb` `vecu"⋅"vecv = u_xv_x + u_yv_y`
Ainsi nous avons :
`vec u "⋅" vec v = vec u^TT vecv`
On peut aussi concevoir des vecteurs de vecteurs, transposés ou non, toutes les configurations sont autorisées. Et lorsque l'on a à faire à un vecteur de vecteurs transposés, ou inversement, un vecteur transposé de vecteurs, de même dimension, alors on l'identifie à une matrice. Par exemples :
`((vecu) , (vecv)) = (( ((u_x),(u_y)) ), (((v_x),(v_y)) ))`
`( (vecu^TT),(vecv^TT) ) = (( ((u_x, u_y)) ), (((v_x, v_y)) )) = (( (u_x, u_y) ), ((v_x, v_y) )) `
`((vecu , vecv)) = ( ((u_x),(u_y))" "((v_x),(v_y)) ) = (( (u_x, v_x) ), ((u_y, v_y) )) `
Les équation locales s'appliquent à un élément différentiel d'espace-temps `dx dy dt` positionné au point `(x,y,t)`. Elles utilisent les opérateurs de dérivée partielle selon chaque axe dont voici la liste :
`(del)/(delx), (del)/(dely), (del)/(delt)`
A partir des ces opérateurs de base, un certain nombre d'autres opérateurs plus complexes sont définis ; la dérivée ou le jacobien `J`, le gradient `vec"grad"` ou le nabla `vecgrad`, la divergence `"div"`, le laplacien `Delta`.
Étant donné un champ scalaire quelconque `varphi`. L'opérateur dérivée partielle selon un axe, par exemple selon `x`, est définie comme suit :
`(del)/(delx)varphi = (varphi(x"+"dx,y,t)-varphi(x,y,t))/dx`
Et les deux notations suivantes sont permisent :
`(del)/(delx)varphi = (delvarphi )/(delx)`
Le champ `varphi`, s'il ne précise pas le point où il est mesuré, alors il est mesuré en sa position par défaut `(x,y,t)`, autrement dit :
`varphi = varphi(x,y,t)`
De même sa dérivée partielle par exemple selon `t` constitue une variable, nommée `(del varphi)/(del t)`, et constitue un champ, et la position par défaut où il est mesuré, est le point `(x,y,t)`. C'est pourquoi les trois écritures suivantes sont valables :
`(del varphi)/(delt) = (del varphi(x,y,t))/(delt) = (del varphi)/(delt)(x,y,t)`
Cela nous permet de calculer le champ au point `(x"+"dx,y"+"dy,t"+"dt)` comme suit par cette égalité que l'on peut démontrer au premier ordre :
`varphi(x"+"dx,y"+"dy,t"+"dt) = varphi + (del varphi)/(delx)dx + (del varphi)/(dely)dy + (del varphi)/(delt)dt`
La position par défaut s'applique également aux composantes manquantes de telle sorte que nous pouvons écrire simplement `varphi(a,b)` pour désigner `varphi(a,b,t)`. Et de même nous pouvons écrire simplement `varphi(a)` pour désigner `varphi(a,y,t)`. Et pour lever l'ambiguité entre produit et application, on note parfois en vert le point où est mesuré le champ comme par exemples `varphi``(x)``, varphi``(x,y)``, varphi``(x,y,t)`.
Cela nous permet de calculer le champ au point `(x"+"dx,y"+"dy)`, sous entendu au point `(x"+"dx,y"+"dy,t)`, comme suit :
`varphi``(x"+"dx,y"+"dy)`` = varphi + (del varphi)/(delx)dx + (del varphi)/(dely)dy`
Et de même, cela nous permet de calculer le champ au point `x"+"dx`, sous entendu au point `(x"+"dx,y,t)`, comme suit :
`varphi``(x"+"dx)`` = varphi + (del varphi)/(delx)dx`
Puis, étant donné un vecteur `vecu` à deux composantes réels quelconques, et un réel `tau` quelconque, on souhaite pouvoir noter `varphi` ou `varphi``(vec u)` ou `varphi``(vec u, tau)` pour désigner `varphi``(u_x, u_y, tau)`. Cela se fait en créant un nouveau contexte par l'instruction suivante :
Contexte : `varphi``(vec u, tau)`` = varphi``(u_x, u_y, tau)`
Le contexte permet non seulement de changer le point de mesure par défaut mais également la forme d'entrée des arguments , ici comme couple du vecteur à deux composantes `x` et `y` et du réel `t`. Aprés cette définition, une expression telle que `varphi``(vec u)` est autorisée et désignera `varphi``(u_x, u_y, tau)`. La position par défaut s'applique dans le cadre de ce nouveau contexte aux listes d'arguments incomplètes. Nous pouvons écrire à partir d'un vecteur `vec w` à deux composantes réels quelconques, l'expression du champ mesuré au point `vec w` et à l'instant `tau`, simplement en notant `varphi``(vec w)` ou de façon complète `varphi``(vec w, tau)`. Et par défaut nous aurons `varphi = varphi``(vec u, tau)`.
Les même considérations sont à faire pour un champ vectoriel `vecv`. La dérivée partielle du champ vectoriel `vecv` selon un axe, par exemple selon `t` est définie par :
`(delvecv)/(delt) = (((delv_x)/(delt)),((delv_y)/(delt)))`
Et nous pouvons calculer le champ vectoriel au point `(x"+"dx,y"+"dy)`, sous-entendu à l'instant `t`, comme suit :
`vecv``(x"+"dx,y"+"dy)`` = vecv + (del vecv)/(delx)dx + (del vecv)/(dely)dy`
Cette équation se met sous une forme utilisant le produit matriciel ; le produit du jacobien par l'élément différentiel de déplacement :
`vecv``(x"+"dx,y"+"dy)`` = vecv + ( ((del vecv)/(delx),(del vecv)/(dely)) )((dx),(dy))`
On définie le jacobien `J` ou la dérivée vectorielle par rapport aux axes d'espace que l'on note comme suit :
`J = (delvecv)/(delvecr) = ( ( (delvecv)/(delx), (delvecv)/(dely) ) ) = ( ( (delv_x)/(delx), (delv_x)/(dely) ),( (delv_y)/(delx), (delv_y)/(dely) ) )`
Puis nous pouvons calculer le champ vectoriel au point d'espace `vecr "+" dr`, sous-entendu à l'instant `t`, comme suit :
Contexte : `vecv``(vec r, t)`` = vecv``(x, y, t)`
`vecv``(vec r "+" dvecr)`` = vecv + (delvecv)/(delvecr) dvecr`
Puis nous pouvons calculer le champ vectoriel au point d'espace-temps `vecX"+" dvecX`, comme suit :
`vecv``(vec X "+" dvecX)`` = vecv + (delvecv)/(delvecX) dvecX`
`(delvecv)/(delvecX) = ( ( (delvecv)/(delx), (delvecv)/(dely), (delvecv)/(delt)) ) = ( ( (delv_x)/(delx), (delv_x)/(dely), (delv_x)/(delt) ),( (delv_y)/(delx), (delv_y)/(dely), (delv_y)/(delt) ) )`
Les opérateurs tels que le jacobien `J`, le nabla `vecnabla` ou gradient `vec"grad"`, la divergence `"div"`, le laplacien `Delta`, sont relatifs aux axes d'espace `x` et `y`. Le nabla est synonyme du gradient. Cet opérateur est défini vectoriellement comme suit :
`vec"grad" = vecnabla = ((del/(delx)),(del/(dely)))`
Le nabla appliqué à un champ scalaire `varphi`, calcul son gradient qui est un champ vectoriel qui indique la plus grande montée :
`vec"grad" varphi = vecgrad varphi = (((delvarphi)/(delx)),((delvarphi)/(dely)))`
Le produit scalaire du nabla avec un champ vectoriel `vecv`, calcule sa divergence qui est un champ scalaire qui exprime la concentration des lignes de courant du champ :
`"div" vecv = vecgrad"⋅"vecv = ((del/(delx)),(del/(dely)))"⋅"((v_x),(v_y)) = (delv_x)/(delx) + (delv_y)/(dely)`
La divergence du gradient d'un champ scalaire `varphi` calcule le laplacien qui est un champ scalaire qui exprime la concentration des lignes de gradient du champ scalaire :
`Delta varphi = "div"(vec"grad" varphi) = vecgrad"⋅"(vecgrad varphi) = (vecgrad"⋅"vecgrad)varphi = vecgrad^2 varphi`
`Delta varphi = ((del/(delx)),(del/(dely)))"⋅"(((delvarphi)/(delx)),((delvarphi)/(dely)))`
`Delta varphi = (del^2varphi)/(delx^2) + (del^2varphi)/(dely^2)`
Le laplacien se définie à partir du nabla comme suit :
`Delta = vecgrad^2 = vecgrad"⋅"vecgrad = ((del/(delx)),(del/(dely)))"⋅"((del/(delx)),(del/(dely)))`
`Delta = (del^2)/(delx^2) + (del^2)/(dely^2)`
Le nabla peut s'appliquer à un champ vectorielle `vecv` :
`vecgrad vecv = (((delvecv)/(delx)),((delvecv)/(dely))) = (( (((delv_x)/(delx)),((delv_y)/(delx))) ),( (((delv_x)/(dely)),((delv_y)/(dely))) ))`
Si le vecteur `vecv` est transposé, on obtient alors le jacobien transposé :
`vecgrad (vecv^TT) = (((delvecv^TT)/(delx)),((delvecv^TT)/(dely))) = (( (((delv_x)/(delx),(delv_y)/(delx))) ),( (((delv_x)/(dely),(delv_y)/(dely))) )) = (( (delv_x)/(delx),(delv_y)/(delx) ),( (delv_x)/(dely),(delv_y)/(dely) ))`
`vecgrad (vecv^TT) = J^TT = ((delvecv)/(del vec r))^TT`
`vecv``(vecr + dvecr)` `= vecv``(vecr)`` + ((delvecv)/(del vec r))d vec r`
Pour lever l'ambiguité entre produit et application, on note parfois en vert le point où est mesuré le champ.
Si on calcul la divergence du nabla appliquée à un champ vectoriel `vecv`, on obtient le laplacien vectoriel, noté par le même symbole delta `Delta` :
`Delta vecv = vecgrad "⋅" (vecgrad vecv) = ((del/(delx)),(del/(dely))) "⋅" (((delvecv)/(delx)),((delvecv)/(dely)))`
`Delta vecv = del/(delx)(delvecv)/(delx)+del/(dely)(delvecv)/(dely) = (del^2vecv)/(delx^2) + (del^2vecv)/(dely^2)`
`Delta vecv = vecgrad^2 vecv = (del^2vecv)/(delx^2) + (del^2vecv)/(dely^2)`
`Delta vecv = ( ((del^2v_x)/(delx^2)),((del^2v_y)/(dely^2)) ) + ( ((del^2v_x)/(delx^2)),((del^2v_y)/(dely^2)) ) = (((del^2v_x)/(delx^2) + (del^2v_x)/(dely^2)),((del^2v_y)/(delx^2) + (del^2v_y)/(dely^2)))`
`Delta vecv = ((Delta v_x),(Delta v_y))`
Le produit scalaire d'un champ vectoriel `vecv` et du nabla, engendre un nouvel opérateur `vecv "⋅"vecnabla` qui peut s'appliquer à un autre champ `vecu` comme suit :
`(vecv "⋅" vecnabla)vecu = (v_xdel/(delx) + v_ydel/(dely))vecu = v_x(delvecu)/(delx) + v_y(delvecu)/(dely)`
Le principe de la résolution passe par un maillage le plus fin possible de l'espace-temps afin que dans chaque maille la vitesse et la pression du fluide soient à peu près identique en tout point de la maille. Puis on transcrit les équations locales en équations aux mailles.
La maille n'étant plus un élément différentiel, si on utilise l'équation locale à une maille telle quelle, cela produira des biais importants toujours dans le même sens. On utilise alors plusieurs fois l'équation locale pour trouver une équation locale qui respecte au moins les symétries sur chaque axe, en `x`, en `y`, et en `t`, et que nous appellons équations aux mailles.
Lorsque `dx` est un élément différentiel et que `x` n'est ni nul ni différentiel, alors nous avons au premier ordre les égalités suivantes :
`x = x+dx`
`varphi = varphi(x+dx) = varphi(x+2dx) = varphi(x+3dx)`
`(del varphi)/(delx) = (varphi-varphi(x"-"dx))/(dx) = (varphi(x"+"dx)-varphi)/(dx)`
A partir de cette règle on déduit des expressions symétriques des dérivés et dérivés secondes qui vont nous permettre de traduire toute équation locale en équation locale symétrique et finalement en équation au mailles :
`(del varphi)/(delx) = (varphi(x"+"dx)-varphi(x"-"dx))/(2dx)`
`(del^2 varphi)/(delx^2) = (varphi(x"+"dx)-2varphi+varphi(x"-"dx))/(dx^2)`
Pour simplifier la notation des coordonnées, on pose que le point par défaut `(x,y,t)` correspond au centre d'une maille, et on crée un système de coordonnées entières associées `(bbx,bby,bbt)` définie par :
`x = bbx dx + dx"/"2`
`y = bby dy + dy"/"2`
`t = bbt dt + dt"/"2`
Mais les éléments `dx, dy, dt` ne sont pas des éléments différentielles. ce sont des éléments de maillage. Nous sommes donc à l'échelle de la maille et non de l'infiniment petit du premier ordre. Le système de coordonnées entières correspond aux indices entiers utilisés par l'ordinateur pour adresser la maille correspondante. Le point de coordonnées `(x,y,z)` correspond au centre d'une maille de coordonnées entières `(bbx,bby,bbz)`.
Néanmoins, on ne souhaite pas utiliser de nouvelles variables pour exprimer les coordonnées entières associées. Pour cela on utilise un artifice, l'usage de crochets spécifiques `(: :)` qui sont par ailleurs de nature ensembliste. Les variables `x,y,t` placées entre ces crochets désigneront leur homologue `bbx,bby,bbt`. Ainsi nous avons créée un second système de coordonnées qui est entier `(:x,y,t:)` et qui porte le même nom que `(x,y,z)` mais qui est différenciée grace au crochet `(: :)`. Et par principe de construction nous avons `(:x,y,t:) = (x,y,t)`. Entre les crochet `(: :)`, toutes les coordonnées ne sont pas obligatoirement mentionnées. Celles non mentionnées ont leur valeur par défaut parmis `(:x,y,t:)`. L'ordre des coordonnées dans la liste n'a pas d'importance. On impose pour permettre cela, une convention d'écriture : Entre les crochet `(: :)`, on commençe l'expression de chaque coordonnée mentionnée toujours par le nom de la coordonnée par défaut, ce qui permet de déterminer de quelle coordonnée il s'agit.
Par exemple, si une seule coordonnée est différente de la valeur par défaut `(:x,y,t:)` on la note de façon condensée en omettant les autres coordonnées et en commençant l'expression par le nom de la coordonnée entière par défaut qui est modifiée et qui sert de valeur de base au calcul :
`varphi(:x"+"1:) = varphi(x"+"dx,y,t)`
`varphi(:y"+"1:) = varphi(x,y"+"dy,t)`
`varphi(:t"+"1:) = varphi(x,y,t"+"dt)`
Autre exemple mettant en exergue la nature ensembliste des crochets : `varphi(:t"+"1, x"-"1:) = varphi(x"-"dx,y,t"+"dt)`
Avec cette convention d'écriture, les versions symétrisées des dérivées partielles se réécrivent comme suit pour produire les équations aux mailles :
`(del varphi)/(delx) = (varphi(:x"+"1:)-varphi(:x"-"1:))/(2dx)`
`(del^2 varphi)/(delx^2) = (varphi(:x"+"1:)-2varphi+varphi(:x"-"1:))/(dx^2)`
Notre maillage en forme topologique de tore, de `N"×"N` mailles, est entièrement défini par la descritpion du voisinage d'une maille `(: :) = (x,y,t)` quelconque, en effectuant le calcul des coordonnées entières selon l'axe des `x` modulo `N` et selon l'axe des `y` modulo `N` :
`( ( (:x"-"1","y"+"1:) , (:y"+"1:), (:x"+"1","y"+"1:) ) , ( (:x"-"1:) , (: :) , (:x"+"1:) ) , ( (:x"-"1","y"-"1:) , (:y"-"1:) , (:x"+"1","y"-"1:) ) )`
La loi de conservation de la masse dite "équation de continuité" ou "équation de bilan des masses" ou encore dans ce cas précis "équation d'imcompressibilité" est :
`vecgrad "⋅" vecv = 0`
La loi de conservation de la quantité de mouvement est de la forme suivante :
Masse `×` Accéleration `=` Force extérieur `+` Force de pression `+` Force de viscosité
`rho(delvecv)/(delt) + rho(vecv*vecnabla)vecv = vecf - vecnablaP+muDeltavecv`
L'évolution de `vecv` est donc régie par cette équation locale :
`rho (delvecv)/(delt) = vecf- vecnablaP + mu Deltavecv - rho (vecv"⋅"vecnabla)vecv`
Si on connait `vecv` et `P` alors on connaitra la dérivé de `vecv` par rapport à `t` grace à la loi de conservation de la quantité de mouvement, et donc on connaitra `vecv(:t+1:)` c'est à dire l'évolution de `vecv` :
| `rho (vecv(:t"+"1:)-vecv)/(dt) = vecf- vecnablaP + mu((vecv(:x"+"1:)-2vecv+vecv(:x"-"1:))/(dx^2) + (vecv(:y"+"1:)-2vecv+vecv(:y"-"1:))/(dy^2)) - rho (v_x(vecv(:x"+"1:)-vecv(:x"-"1:))/(2dx) + v_y(vecv(:y"+"1:)-vecv(:y"-"1:))/(2dy))` |
| `rho (v_x(:t"+"1:)-v_x)/(dt) = f_x - (P(:x"+"1:)-P(:x"-"1:))/(2dx) + mu((v_x(:x"+"1:)-2v_x+v_x(:x"-"1:))/(dx^2) + (v_x(:y"+"1:)-2v_x+v_x(:y"-"1:))/(dy^2)) - rho (v_x(v_x(:x"+"1:)-v_x(:x"-"1:))/(2dx) + v_y(v_x(:y"+"1:)-v_x(:y"-"1:))/(2dy))` `rho (v_y(:t"+"1:)-v_y)/(dt) = f_y - (P(:y"+"1:)-P(:y"-"1:))/(2dy) + mu((v_y(:x"+"1:)-2v_y+v_y(:x"-"1:))/(dx^2) + (v_y(:y"+"1:)-2v_y+v_y(:y"-"1:))/(dy^2)) - rho (v_x(v_y(:x"+"1:)-v_y(:x"-"1:))/(2dx) + v_y(v_y(:y"+"1:)-v_y(:y"-"1:))/(2dy))` |
Mais il nous manque encore l'évolution de `P`. On revient sur les équations locales pour extraire directement une troisième équation locale qui nous permettrait de calculer l'évolution du champ de pression `P`.
`vecgrad "⋅" vecv = 0`
`(del)/(delt) vecgrad "⋅" vecv = 0`
`vecgrad "⋅" (delvecv )/(delt) = 0`
D'autre part
`rho (delvecv)/(delt) = vecf- vecnablaP+mu Deltavecv - rho (vecv"⋅"vecnabla)vecv`
Donc
`vecgrad "⋅" (rho (delvecv )/(delt)) = 0`
`vecgrad "⋅" (vecf - vecnablaP+mu Delta vecv - rho(vecv"⋅"vecnabla)vecv) = 0`
`vecgrad "⋅" vecf - DeltaP + mu vecgrad * Deltavecv - rho vecgrad "⋅" (vecv"⋅"vecnabla)vecv = 0`
`DeltaP = vecgrad "⋅" vecf + mu vecgrad * Deltavecv - rho vecgrad "⋅" (vecv"⋅"vecnabla)vecv`
Si on connait le laplacien `DeltaP` en tout point de l'espace et si on connait une valeur de `P` en un point, alors on connait `P` en tout point de l'espace, par prolongement analytique. C'est une propriété des champs analytiques c'est à dire indéfiniment dérivables. Or nous n'avons pas évoqué la moindre hypothèse de discontinuité sur quelque plan que ce soit. Nos champs sont, de façon sous-entendu, tous analytiques. Le laplacien de `P` étant connue, la détermination de `P` peut s'obtenir par prolongement analytique. Et la méthode de résolution se fait par un calcul itératif résolvant de façon approchée l'équation aux mailles correspondant au laplacien de `P` :
`L = vecgrad "⋅" vecf + mu vecgrad * Deltavecv - rho vecgrad "⋅" (vecv"⋅"vecnabla)vecv` `DeltaP = L`
`(del^2P)/(delx^2) + (del^2P)/(dely^2) = L`
`(P(:x"+"1)-2P+P(:x"-"1:))/(dx^2) + (P(:y"+"1:)-2P+P(:y"-"1:))/(dy^2) = L`
A partir des pressions connues précédemment c'est à dire de `P(:t"-"1:)`. Sur chaque point, on répartie la différence d'inégalité sur les 6 termes comme on ferait une retropropagation d'erreur dans un réseau de neurones, où la maille joue le rôle de neurone.
`L = vecgrad "⋅" vecf + mu vecgrad * Deltavecv - rho vecgrad "⋅" (vecv"⋅"vecnabla)vecv`
`vecgrad "⋅" vecf = (f_x(:x"+"1:)-f_x(:x"-"1:))/(2dx) + (f_y(:y"+"1:)-f_y(:y"-"1:))/(2dy)`
`vecgrad * Deltavecv = (delDeltav_x)/(delx) + (delDeltav_y)/(dely) = (del((del^2v_x)/(delx^2) + (del^2v_x)/(dely^2)))/(delx) + (del( (del^2v_y)/(delx^2) + (del^2v_y)/(dely^2)))/(dely)`
`vecgrad * Deltavecv = del/(delx)(del^2v_x)/(delx^2) + del/(delx)(del^2v_x)/(dely^2) + del/(dely)(del^2v_y)/(delx^2) + del/(dely)(del^2v_y)/(dely^2)`
`vecgrad * Deltavecv = ubrace( del/(delx)((v_x(:x"+"1:)-2v_x+v_x(:x"-"1:))/(dx^2)))_(K_1) + ubrace( del/(delx)((v_x(:y"+"1:)-2v_x+v_x(:y"-"1:))/(dy^2)))_(K_2) + ubrace( del/(dely)((v_y(:x"+"1:)-2v_y+v_y(:x"-"1:))/(dx^2)))_(K_3) + ubrace( del/(dely)((v_y(:y"+"1:)-2v_y+v_y(:y"-"1:))/(dy^2)))_(K_4)`
`K_1 = (v_x(:x"+"2:)-2v_x(:x"+"1:)+v_x)/(2dx^3) - (v_x-2v_x(:x"-"1:)+v_x(:x"-"2:))/(2dx^3)`
`K_2 = (v_x(:x"+"1,y"+"1:)-2v_x(:x"+"1:)+v_x(:x"+"1,y"-"1:))/(2dxdy^2) - (v_x(:x"-"1,y"+"1:)-2v_x(:x"-"1:)+v_x(:x"-"1,y"-"1:))/(2dxdy^2)`
`K_3 = (v_y(:x"+"1,y"+"1:) - 2v_y(:y"+"1:)+v_y(:x"-"1,y"+"1:))/(2dx^2dy) - (v_y(:x"+"1,y"-"1:) - 2v_y(:y"-"1:)+v_y(:x"-"1,y"-"1:))/(2dx^2dy)`
`K_4 = (v_y(:y"+"2:) - 2v_y(:y"+"1:)+v_y)/(2dy^3) - (v_y - 2v_y(:y"-"1:)+v_y(:y"-"2:))/(2dy^3)`
`vecgrad "⋅" (vecv"⋅"vecnabla)vecv = (del( v_x(delv_x)/(delx) + v_y(delv_x)/(dely) ))/(delx) + (del( v_x(delv_y)/(delx) + v_y(delv_y)/(dely) ))/(dely) `
`vecgrad "⋅" (vecv"⋅"vecnabla)vecv = ubrace( del/(delx)(v_x(delv_x)/(delx)))_(Q_1) + ubrace( del/(delx)(v_y(delv_x)/(dely)))_(Q_2)+ ubrace( del/(dely)(v_x(delv_y)/(delx)))_(Q_3) + ubrace( del/(dely)(v_y(delv_y)/(dely) ))_(Q_4)`
`Q_1 = (v_x(:x"+"1:)(v_x(:x"+"2:)-v_x))/(4dx^2) - (v_x(:x"-"1:)(v_x-v_x(:x"-"2:)))/(4dx^2)`
`Q_2 = (v_y(:x"+"1:)(v_x(:x"+"1,y"+"1:)-v_x(:x"+"1,y"-"1:)))/(4dxdy) - (v_y(:x"-"1:)(v_x(:x"-"1,y"+"1:)-v_x(:x"-"1,y"-"1:)))/(4dxdy)`
`Q_3 = (v_x(:y"+"1:)(v_y(:x"+"1,y"+"1:)-v_y(:x"-"1,y"+"1:)))/(4dxdy) - (v_x(:y"-"1:)(v_y(:x"+"1,y"-"1:)-v_y(:x"-"1,y"-"1:)))/(4dxdy)`
`Q_4 = (v_y(:y"+"1:)(v_y(:y"+"2:)-v_y))/(4dy^2) - (v_y(:y"-"1:)(v_y-v_y(:y"-"2:)))/(4dy^2)`
----- 24 septembre 2018 -----
Mais il y a un condition supplémentaire. Le champ de vitesse du courant `vecv` doit vérifier la loi de conservation de la masse, ce qui n'est pas le cas si on pose des valeurs de pression `P` arbitraires.
A priori un fluide incompressible transmet la pression à la vitesse du son. Dans notre modèle théorique, on suppose que la pression se transmet instantanément. Un calcul est donc nécessaire pour déterminer cette pression. Ce champ de pression `P` est tel que la loi de conservation de la quantité de mouvement produise un champ `vecv(:t+1:)` vérifiant justement l'équation de conservation des masses. C'est pourquoi on suppose inconnue ce champ de pression `P`, et que le calcul du champ `vecv(:t+1:)` se fait en même temps le calcul du champ de pressions `P`.
Le nombre de mailles va vite épuisé la capacité de calcule de l'ordinateur. Or le fluide entretient des zones où l'état du fluide est quasi-constant et ne nécessite que de quelques grandes mailles pour être modélisé tandis que dans d'autres zones, le fluide génère des tourbillons ou de forts gradients de vitesse ce qui nécessite un grand nombres de petites maille pour être correctement modélisé.
On commence par étudier les fluides incompressibles c'est à dire de densité `rho` constante, et visqueux c'est à dire d'inertie négligeable devant la viscosité, autrement dit, dont le nombre de Reynolds, le rapport entre les forces d'inertie et les forces visqueuses, est très petit devant `1` :
`(rhoVL)/mu -< 1`
`rho` : Densité du fluide.
`mu` : Viscosité dynamique du fluide.
`V` : Vitesse du fluide.
`L` : Largeur du système.
Le fluide est décrit par un champ de pression `P` et un champ de vitesse `vecv`. Si on ne précise pas les coordonnées du points où est mesurée le champ, par défaut ce sera en position `(x,y,t)`. Ainsi nous avons :
`P = P(x,y,t)`
`vecv=vecv(x,y,t)`
On connait `vecv(:t-1:), P(:t-1:), vecv` et on calcul `P, vecv(:t+1:)`. Il y a `3N^2` équations et `3N^2` inconnus. On résoud le système d'équations de façon itérative jusqu'à obtenir la précision voulue. On formule pour chaque maille la contrainte devant être respectée sous forme d'une fonction `f` d'un certain nombres de pressions `P_1, P_2, P_3,..., P_n` des mailles autours, et devant être nulle. On part du champ de pression `P:=P(:t-1:)`. L'erreur est `epsilon`.
`f(P_1,P_2,P_3,...,P_n) = epsilon`
`df = (delf)/(dP_1)dP_1 + (delf)/(dP_2)dP_2 + (delf)/(dP_3)dP_3 + ... + (delf)/(dP_n)dP_n`
On recherche les variations `dP_n` telque `df = -epsilon`. On choisie la solution équilibrée :
`(delf)/(dP_i)dP_i = epsilon/n`
`dP_i = epsilon/(n(delf)/(dP_i))`
`P_i:= P_i+dP_i`
Puis après avoir parcourue toutes les mailles, on réitère le même calcul jusqu'à ce que tous les `epsilon` soit inférieurs à la précision demandée.