TP2 - Problème de Stokes - JULIA


Ce petit "notebook" a pour but de montrer quelques propriétés de résolution d'un systeme linéaire de type "Elliptique contraint". Le langage utilisé est Julia, l'interface utilisée -le notebook- est Ipython avec le package IJulia. Les graphiques sont réalisés à l'aide du package PyPlot.

Les différents tests (benchmarks) seront menés sur un système d'équation provenant de la discrétisation différence finie (schéma de MAC) du problème se Stockes suivant :

$$ \begin{cases} -\mu \Delta \vec{u} + \nabla p = \vec{f} &\text{ dans } \Omega\\ \text{div}(\vec{u}) = 0&\text{ dans } \Omega \\ \vec{u}=\vec{g}& \text{ sur } \partial \Omega \quad \text{ avec } \int_{\partial \Omega} \vec{g}.\vec{n} d\sigma =0 \end{cases} $$


1. Assemblage du système

On se donne une grille uniforme sur $]0,1[^2$ avec un pas $h=\frac{1}{n+1}$. Dans un premier temps prenons ne nous soucions pas des conditions de Dirichlet ($\vec{g}$).

On note

  • $M_{i,j}=(ih,jh)$, $i,j=0,..,n+1$ les noeuds de la grille.
  • $C_{i+\frac{1}{2},j+\frac{1}{2}}=\left(\left(i+1/2\right)h,\left(j+1/2\right)h\right)$, $i,j=0,..,n$ les centres des cellules de la grille.

On cherche

  • les composantes de la vitesse aux noeuds de la grille $u^1_{i,j}\simeq u^1(M_{i,j})$, $u^2_{i,j}\simeq u^2(M_{i,j})$ $i,j=0,..,n+1$.
  • la pression au centre des cellules $p_{i+\frac{1}{2},j+\frac{1}{2}}\simeq p(C_{i+\frac{1}{2},j+\frac{1}{2}})$ $i,j=0,..,n$

Le schéma de MAC donne le système suivant : $$ \begin{cases} - \mu \Delta_h u^1_{i,j} + \frac{1}{2h}\left(p_{i+\frac{1}{2},j+\frac{1}{2}}+p_{i+\frac{1}{2},j-\frac{1}{2}}-p_{i-\frac{1}{2},j+\frac{1}{2}}-p_{i-\frac{1}{2},j-\frac{1}{2}} \right) = f^1_{i,j} & i,j=0,..,n+1\\- \mu \Delta_h u^2_{i,j} + \frac{1}{2h}\left(p_{i+\frac{1}{2},j+\frac{1}{2}}+p_{i-\frac{1}{2},j+\frac{1}{2}}-p_{i+\frac{1}{2},j-\frac{1}{2}}-p_{i-\frac{1}{2},j-\frac{1}{2}} \right) = f^2_{i,j} & i,j=0,..,n+1 \\\frac{1}{2h}\left(u^1_{i,j}+u^1_{i,j+1}-u^1_{i+1,j}-u^1_{i+1,j+1} \right) + \frac{1}{2h}\left(u^2_{i,j}+u^2_{i+1,j}-u^2_{i,j+1}-u^2_{i+1,j+1} \right)=0 & i,j=0,..,n \end{cases}$$

avec $-\Delta_h u_{i,j}=\frac{1}{h^2}\left(-u_{i-1,j} - u_{i,j-1} +4u_{i,j} -u_{i+1,j} -u_{i,j+1}\right)$ et les valeurs $p_{-\frac{1}{2},..}=p_{..,-\frac{1}{2}}=p_{n+\frac{3}{2},..}=p_{..,n+\frac{3}{2}}=0$.

Ceci donne un système matriciel de la forme : $$ \left(\begin{array}{c c c}A_1 & 0 &B_1^T\\0 & A_2 &B_2^T \\B_1& B_2& 0 \end{array}\right)\left(\begin{array}{c}u^1\\u^2\\p\end{array}\right)=\left(\begin{array}{c}f_1\\f_2\\0\end{array}\right) $$

avec

  • $A_1$, $A_2$ matrices de taille $(n+2)^2\times(n+2)^2$
  • $B_1$; $B_2$ matrices de taille $(n+2)^2\times(n+1)^2$

In [1]:
function Laplace1d(n)
    h=1./(n+1)
    d1=[-1/h^2 for i=1:n+1];
    d2=[2/h^2 for i=1:n+2];
    d3=[-1/h^2 for i=1:n+1];
    spdiagm((d1,d2,d3),(-1,0,1));   
end
function Laplace2d(n,m)
    kron(speye(m+2),Laplace1d(n))+kron(Laplace1d(m),speye(n+2))
end
# différentiation 1d (forward)
function Dx1d(n)
    h=1./(n+1)
    d1=[-1/h for i=1:n+1];
    d2=[1/h for i=1:n+1];
    spdiagm((d1,d2),(-1,0))';
end
# différentiation 2d
function Dx2d(n,m)
    d1=ones(m+1);
    d2=ones(m+1);
    tmp=spdiagm((d1/2,d2/2),(-1,0))
    kron(tmp,Dx1d(n))
end
function Dy2d(n,m)
    tmp=spdiagm((ones(n+1)/2,ones(n+1)/2),(-1,0))
    kron(Dx1d(m),tmp)
end


Out[1]:
Dy2d (generic function with 1 method)

2. Conditions aux limites

Pour l'instant nous n'avons pas pris en compte les conditions aux limites autrement dit il n'est pas tenu compte des éléments faisant parti de la frontière.

Supposons que l'on ai numéroté les sommets de tel sorte que le système s'écrive $$\left(\begin{array}{rr}A_{11} & A_{12}\\A_{21} & A_{22}\end{array}\right) \left(\begin{array}{l}V\\V_D\end{array}\right)=\left(\begin{array}{l}b\\b_D\end{array}\right)$$ avec $V_D$ les valeurs de $V$ sur la frontière donc connus ! Ce ne sont pas des inconnues on peut écrire formellement $V_D=g$

Atrement dit on doit résoudre le système : $$ A_{11} V=b-A_{12}g$$

Une autre solution est de travailler par pénalisation, on garde le premiès système et on ajoute 2 termes : $$\left(\begin{array}{rr}A_{11} & A_{12}\\A_{21} & A_{22}+\frac{1}{\varepsilon}Id\end{array}\right) \left(\begin{array}{l}V\\V_D\end{array}\right)=\left(\begin{array}{l}b\\b_D+\frac{1}{\varepsilon}g\end{array}\right)$$

En prenant $\varepsilon<<1$ on fixe numériquement la valeur de $V_D$ car on a de manière équivallente $$\begin{cases} A_{11}V + A_{12}V_D = b\\ \left(Id+{\varepsilon}A_{22}\right)V_D=g+ {\varepsilon}\left(A_{21}V+b_D\right)\end{cases}$$

Et donc $V_D=g+o(\varepsilon)$, le choix de $\varepsilon$ va dépendre du rayon spectral de $A_{11}$.

En prenant $\varepsilon=\rho(A_{11})^{-1}\times10^{-20}$ donne $\rho\left(Id+{\varepsilon}A_{22}\right)\leq 1+10^{-20}$...



In [2]:
n=3; #
mu=1e-2;
A1=mu*Laplace2d(n,n) #Matrice de diffusion pour la première composante de la vitesse 
A2=mu*Laplace2d(n,n) #Matrice de diffusion pour la seconde composante de la vitesse 
B1=Dx2d(n,n) 
B2=Dy2d(n,n)
N=(n+2)^2;
A=[A1,spzeros(N,N),B1;spzeros(N,N), A2 , B2 ; B1' , B2' , spzeros((n+1)^2,(n+1)^2)]


syntax: unexpected semicolon in array expression
while loading In[2], in expression starting on line 8

Problème de la Cavité entrainée :

Dans cet exemple, on impose une vitesse nulle partout sur le bord du carré $]0,1[^2$ sauf sur le bord supérieur $y = 1$ où la vitesse est horizontale.

Plus précisément, on choisit

  • $u =\left(\begin{array}{c}0\\0\end{array}\right)$ sur les bords $x = 0$; $x = 1$ et $y = 0$.
  • $u =\left(\begin{array}{c}15\\0\end{array}\right)$ sur le bord $y = 1$.
  • Avec une force de gravité verticale $f =\left(\begin{array}{c}0\\-30\end{array}\right)$
  • La viscosité $\mu = 10^{-􀀀2}$

In [8]:
full(Dx1d(N))


Out[8]:
26x27 Array{Float64,2}:
 26.0  -26.0    0.0    0.0    0.0  …    0.0    0.0    0.0    0.0    0.0
  0.0   26.0  -26.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0   26.0  -26.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0   26.0  -26.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0   26.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0  …    0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0  …    0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0  …    0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0  …    0.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0     -26.0    0.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0      26.0  -26.0    0.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0   26.0  -26.0    0.0    0.0
  0.0    0.0    0.0    0.0    0.0       0.0    0.0   26.0  -26.0    0.0
  0.0    0.0    0.0    0.0    0.0  …    0.0    0.0    0.0   26.0  -26.0

Problème de poiseuil


In [4]:
full(B1[2,1:12])


Out[4]:
1x12 Array{Float64,2}:
 0.0  2.0  -2.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

In [10]:
size(B2)


Out[10]:
(20,20)

In [58]:
size(A1)


Out[58]:
(25,25)

In [ ]: