Stokesův problém modeluje stacionární velmi pomalé proudění vazké nestlačitelné tekutiny. Je posán následujícími rovnicemi
\begin{eqnarray} \nabla\cdot U &=& 0, \\ \nabla p &=& \nabla \cdot(\mu \nabla U). \end{eqnarray}kde $p$ je tlak, $U$ je rychlost (vektor) a $\mu$ je dynamická viskozita.
Pozn. pro $\mu=konst$ lze ve druhé rovnici nahradit tlak součinem $\mu p$ a rovnice se pak zjednoduší na $\nabla p = \Delta U$.
In [1]:
using PyPlot;
In [2]:
include("mesh.jl");
include("fields.jl");
include("operators.jl");
In [3]:
n = 25
m = cartesian_mesh(n,n);
In [4]:
plot_mesh(m)
axis("equal");
In [5]:
p = ScalarField(m);
for name ∈ ["left", "right", "bottom", "top"]
set_neumann_patch!(p, name)
end
In [6]:
u = ScalarField(m)
set_dirichlet_patch!(u, "left", 0.0);
set_dirichlet_patch!(u, "right", 0.0);
set_dirichlet_patch!(u, "bottom", 0.0);
set_dirichlet_patch!(u, "top", 1.0);
v = ScalarField(m)
set_dirichlet_patch!(v, "left", 0.0);
set_dirichlet_patch!(v, "right", 0.0);
set_dirichlet_patch!(v, "bottom", 0.0);
set_dirichlet_patch!(v, "top", 0.0);
Soustava rovnic pro řešení Stokesova problému je $$ \begin{bmatrix} L & 0 & G_x \\ 0 & L & G_y \\ D_x & D_y & 0 \end{bmatrix}
\begin{bmatrix} u \\ v \\ p \end{bmatrix}
=
\begin{bmatrix} b_u \\ b_v \\ b_0 \end{bmatrix}.
$$ kde $L$ je matice vzniklá aproximací $-\Delta$, $G_x$ je aproximace $\partial_x$ a $G_y$ je aproximace $\partial_y$. Podobně i $D_x$ a $D_y$.
In [7]:
μ = 1.e-4;
In [8]:
eq1 = -Δ(μ, u);
eq2 = -Δ(μ, v);
eq3 = ∂x(p);
eq4 = ∂y(p);
eq5 = ∂x(u);
eq6 = ∂y(v);
In [9]:
L = eq1.A;
Gx = eq3.A;
Gy = eq4.A;
Dx = eq5.A;
Dy = eq6.A;
O = spzeros(length(cells(m)),length(cells(m)));
M = [L O Gx;
O L Gy;
Dx Dy O];
bu = -eq1.b - eq3.b;
bv = -eq2.b - eq4.b;
b0 = -eq5.b - eq6.b;
b = [bu; bv; b0];
Matice $M$ je ale singulární!
Lze to ověřit například výpočtem čísla podmíněnosti $$ \kappa_A = ||A|| ||A^{-1}||. $$
To by sice pro singulární matice nemělo být definované (viz $A^{-1}$ ve vzorci), nicméně pro téměř singulární matice je $\kappa$ veliké. Při odhadu $\kappa$ tak může dojít buď k chybe výpočtu z důvodu singularity matice a nebo díky zaokrouhlovací chybě výpočet proběhne a vypočtený odhad $\kappa$ bude veliký.
In [10]:
# Odhad κ
cond(Array(M))
Out[10]:
Přidám proto k soustavě 1 rovnici p[end]=0. Touto rovnicí fixuji hodnotu tlaku v jednom bodě sítě. Singularita soustavy je totiž způsobena tím, že tlak je ve Stokesově problému určen až na konstantu.
In [11]:
M[end,end] = n*n;
cond(Array(M))
Out[11]:
In [12]:
sol = M \ b;
Řešení v poli sol teď obsahuje složky [u, v, p]. Tyto složky načteme z příslušných částí sol.
In [13]:
uu = sol[1:n*n];
vv = sol[n*n+1:2*n*n];
pp = reshape(sol[2*n*n+1:end], (n,n))';
In [14]:
cx = [v[1] for v in m.centre]
cy = [v[2] for v in m.centre]
contourf(reshape(cx,(n,n))', reshape(cy,(n,n))', pp)
colorbar(); axis("equal")
quiver(cx,cy,uu,vv);
Na jednotkovém čtevrci [0,1]x[0,1] řešme soustavu rovnic \begin{align*}
s okrajovými podmínkami $$ U(x,y) = \left\{ \begin{array}{ll} (1,0) & \text{pro } x \in [0,1] \land y=1,\\ (0,0) & \text{pro } x=0 \lor x=1 \lor y=0. \end{array} \right. $$ a $$ \frac{\partial p}{\partial n} = 0 \,\text{pro } x=0 \lor x=1 \lor y=0 \lor y=1. $$
Diskretizace: $$ a_C U_C = \sum_{f} a_F U_F + Q_C - \nabla p_C = H(U)_C - \nabla p_C, $$
a tedy
$$ U_C = \frac{1}{a_C} H(U)_C - \frac{1}{a_C} \nabla p_C. $$Definujme $$ \hat{U}_C = \frac{1}{a_C} H(U)_C. $$
Rovnice kontinuity
$$ 0 = \nabla \cdot U_C = \nabla \cdot (\hat{U}_C - \frac{1}{a_C} \nabla p_C) $$a tedy
$$ \nabla\cdot(\frac{1}{a_C} \nabla p_C) = \nabla \cdot \hat{U}_C. $$Výsledná rovnice pro tlak však stále obsahuje neznámé rychlosti. Algoritmus SIMPLE řeší danou soustavu rovnic iteračně s tím, že se řeší zvlášť rovnice pro rychlosti a zvlášť rovnice pro tlak.
V této variantě algoritmu SIMPLE vycházíme ze stavu $U^n$ a $p^n$ a pokusíme se vypočítat novou aproximaci $U^{n+1}$ a $p^{n+1}$.
Tato varianta je velmi podobná předchozímu algoritmu. Narozdíl od něj však nepočítá s korekcemi tlaku ale přímo s tlakem.
In [15]:
U = VectorField(m)
set_dirichlet_patch!(U, "left", Vector(0,0));
set_dirichlet_patch!(U, "right", Vector(0,0));
set_dirichlet_patch!(U, "bottom", Vector(0,0));
set_dirichlet_patch!(U, "top", Vector(1,0));
p ← 0.0;
Ručně provedeme 1 iteraci algoritmu SIMPLE:
In [16]:
Ueqn = - Δ(μ, U);
In [17]:
solve!(Ueqn + grad(p));
In [18]:
ra = 1 ./ Ac(Ueqn);
In [19]:
Ubar = VectorField(ra .* H(Ueqn), U.mesh, U.boundaries);
In [20]:
ra .* H(Ueqn);
In [21]:
pEqn = Δ(ra, p) - div(Ubar);
Diskretizovaná rovnice pro tlak s Neumannovskými podmínkami má však singulární matici. Proto přidáme k rovnici v bodě 1 $K$-násobek rovnice $p_1=0$. To se projeví přidáním $K$ k prvku $A[1,1]$. Hodnotu $K$ přitom volíme tak, aby (zhruba) odpovídala velikosti $A[1,1]$.
Z rovnice pro rychlost vidíme, že na kartézské síti ve 2D je ve vnitřním bodě sítě $$ ra = \frac{h^2}{4\mu} $$
Dosazením do rovnice pro tlak je v tomto bodě $$ A[1,1] = - \frac{h^2}{4\mu} \frac{4}{h^2} = - \frac{1}{\mu}. $$
In [22]:
# diagonalni prvek pro bunku odpovidajici cca stredu oblasti
pEqn.A[312,312]
Out[22]:
In [23]:
-1/μ
Out[23]:
In [24]:
# diagonalni prvek v bode 1 (ten je ovlivnen okrajovymi podminkami)
pEqn.A[1,1]
Out[24]:
In [25]:
solve!(pEqn)
In [26]:
plot_contourf(p)
axis("equal"); colorbar();
In [27]:
U ← Ubar - ra .* grad(p);
In [28]:
component(f::VectorField, c) = [ v[c] for v in f.values ]
Out[28]:
In [29]:
plot_contourf(p); colorbar();
plot_arrows(U);
axis("equal")
Out[29]:
Výše uvedený postup bychom měli opakovat dokud řešení nezkonverguje ke stacionárnímu stavu.
Zkusíme provést několik iterací:
In [30]:
p ← 0.0
U ← Vector(0,0);
for iter = 0:50
U_old, p_old = copy(U.values), copy(p.values)
Ueqn = -Δ(μ,U)
solve!(Ueqn + grad(p))
ra = 1 ./ Ac(Ueqn);
Ubar = VectorField(ra .* H(Ueqn), U.mesh, U.boundaries);
pEqn = Δ(ra, p) - div(Ubar);
pEqn.A[1,1] -= 1/μ
solve!(pEqn)
U ← Ubar - ra .* grad(p)
if rem(iter,5)==0
nxny = n*n
pRez = norm(p_old - p.values) / nxny
URez = norm(U_old - U.values) / nxny
println(iter, "\t", pRez, "\t", URez)
end
end
Je vidět, že algoritmus DIVERGUJE
Významnou roli v algoritmu SIMPLE hraje řešení soustav rovnic pro $U$ a $p$. Zatímco soustava rovnic pro $p$ vznikne diskretizací Poissonovy rovnice (a je tedy při vhodné diskretizaci symetrická a pozitivně definitní), soustava pro $U$ už symetrická být nemusí (díky konvektivním členům). Navíc v rovnici pro $U'$ zanedbáváme člen $H(U')/a_C$ a chtěli bychom, aby tento člen byl vzhledem k $U'$ malý. Proto je nutné v rovnici pro rychlost použí relaxaci.
Mějme soustavu lineárních rovnic $Ax=b$. Označme $x^*$ přesné řešení této soustavy a $x^n$ $n$-tou iteraci níže popsané relaxační metody. Tu získáme tak, že k řešíme místo původní soustavy soustavu vzniklou součtem $A x^{n+1} = b$ s vhodným násobek $Dx^{n+1}=Dx^n$, kde $D=diag(A)$, tedy \begin{equation} A x^{n+1} + \frac{1-\alpha}{\alpha} D x^{n+1} = b + \frac{1-\alpha}{\alpha} D x^n, \end{equation} neboli \begin{equation} \left(\frac{D}{\alpha} + (A-D)\right) x^{n+1} = b + \frac{1-\alpha}{\alpha} D x^n. \end{equation} Při vhodné volbě $\alpha \in (0,1)$ je matice relaxované soustavy ODD.
To má význam i pro odhad velikosti $H(U')$. Je-li $U' \approx const.$, pak $H(U')_C/a_C \approx \alpha U'_C$ a tedy lze tento člen v algoritmu SIMPLE zanedbat.
Dalším místem, kde v algoritmu SIMPLE zanedbáváme člen $H(U')$, je vztah pro $U'_C$. Provedeme-li relaxaci rovnice pro rychlost, pak lepší odhad korekce rychlosti bude $$ U'_C = \frac{1}{a_C} H(U')_C - \frac{1}{a_C} \nabla p'_C \approx \alpha U'_C - \frac{1}{a_C} \nabla p'_C, $$ tedy pro $\beta = 1 - \alpha$ $$ U'_C = - \frac{1}{\beta a_C} \nabla p'_C. $$ Označme $p" = p'_C / \beta$. Potom $U'_C = - \frac{1}{a_C} \nabla p"_C$ a $p"$ najdeme řešením původní rovnice \begin{equation} \nabla\cdot(\frac{1}{a_C}\nabla p") = \nabla\cdot U^*. \end{equation} Tlak $p^{n+1}$ nakonec upravíme dle vztahu $p^{n+1} = p^n + p' = p^n + \beta p"$.
Ve variantě s korekcí rychlosti tedy
Vhodné relxační parametry jsou $\alpha + \beta \approx 1$, často $\alpha=0.7$ a $\beta = 0.3$.
Celý algoritmus SIMPLE včetně relaxace bude tedy ve variantě s výpočtem tlaku vypadat následovně:
In [31]:
p ← 0.0
U ← Vector(0,0);
α = 0.7
β = 0.3
for iter = 0:50
U_old, p_old = copy(U.values), copy(p.values)
Ueqn = -Δ(μ,U)
relax!(Ueqn, α)
solve!(Ueqn + grad(p))
ra = 1 ./ Ac(Ueqn);
Ubar = VectorField(ra .* H(Ueqn), U.mesh, U.boundaries);
pEqn = Δ(ra, p) - div(Ubar);
pEqn.A[1,1] -= α/μ
solve!(pEqn)
p ← β*p + (1-β)*p_old
U ← Ubar - ra .* grad(p)
if rem(iter,5)==0
nxny = n*n
pRez = norm(p_old - p.values) / nxny
URez = norm(U_old - U.values) / nxny
println(iter, "\t", pRez, "\t", URez)
end
end
In [32]:
plot_contourf(p); colorbar();
plot_arrows(U);
axis("equal")
Out[32]:
In [33]:
uu = reshape(component(U,1), (n,n));
vv = reshape(component(U,2), (n,n));
pp = reshape(p.values, (n,n));
In [34]:
# Ψ_x = v, Ψ_y = - u
function streamfunction(u, v, Δx=1, Δy=1)
Ψ1 = zero(u)
Ψ1[1,1] = 0.0
for i=2:size(u,1)
Ψ1[i,1] = Ψ1[i-1,1] + Δx * (v[i,1]+v[i-1,1])/2
end
for j=2:size(u,2)
Ψ1[:,j] = Ψ1[:,j-1] - Δy * (u[:,j]+u[:,j-1])/2
end
Ψ2 = zero(Ψ1)
Ψ2[1,1] = 0.0
for j=2:size(u,2)
Ψ2[1,j] = Ψ2[1,j-1] - Δy * (u[1,j]+u[1,j-1])/2
end
for i=2:size(u,1)
Ψ2[i,:] = Ψ2[i-1,:] + Δx * (v[i,:]+v[i-1,:])/2
end
Ψ = (Ψ1 + Ψ2) / 2
return Ψ
end
Out[34]:
In [35]:
Ψ = streamfunction(uu,vv,1/25,1/25);
In [36]:
contour(reshape(cx,(n,n)), reshape(cy,(n,n)), reshape(Ψ,(n,n)), range(0,maximum(Ψ),length=10), colors="black");
contour(reshape(cx,(n,n)), reshape(cy,(n,n)), reshape(Ψ,(n,n)), range(minimum(Ψ),0,length=10), colors="red");
plot_contourf(p); colorbar();
axis("equal");
Totéž na jemnější síti
In [37]:
nf=50;
mf = cartesian_mesh(nf,nf);
In [38]:
pf = ScalarField(mf);
for name ∈ ["left", "right", "bottom", "top"]
set_neumann_patch!(pf, name)
end
In [39]:
Uf = VectorField(mf)
set_dirichlet_patch!(Uf, "left", Vector(0,0));
set_dirichlet_patch!(Uf, "right", Vector(0,0));
set_dirichlet_patch!(Uf, "bottom", Vector(0,0));
set_dirichlet_patch!(Uf, "top", Vector(1,0));
In [40]:
α = 0.7
β = 0.3
for iter = 0:200
U_old, p_old = copy(Uf.values), copy(pf.values)
Ueqn = -Δ(μ,Uf)
relax!(Ueqn, α)
solve!(Ueqn + grad(pf))
ra = 1 ./ Ac(Ueqn);
Ubar = VectorField(ra .* H(Ueqn), Uf.mesh, Uf.boundaries);
pEqn = Δ(ra, pf) - div(Ubar);
pEqn.A[1,1] -= α/μ
solve!(pEqn)
pf ← β*pf + (1-β)*p_old
Uf ← Ubar - ra .* grad(pf)
if rem(iter,10)==0
nxny = nf*nf
pRez = norm(p_old - pf.values) / nxny
URez = norm(U_old - Uf.values) / nxny
println(iter, "\t", pRez, "\t", URez)
end
end
In [41]:
uuf = reshape(component(Uf,1), (nf,nf));
vvf = reshape(component(Uf,2), (nf,nf));
ppf = reshape(pf.values, (nf,nf));
In [42]:
Ψf = streamfunction(uuf,vvf,1/nf,1/nf);
In [43]:
cxf = [v[1] for v in mf.centre];
cyf = [v[2] for v in mf.centre];
In [44]:
contour(reshape(cxf,(nf,nf)), reshape(cyf,(nf,nf)), reshape(Ψf,(nf,nf)), range(0,maximum(Ψf),length=10), colors="black");
contour(reshape(cxf,(nf,nf)), reshape(cyf,(nf,nf)), reshape(Ψf,(nf,nf)), range(minimum(Ψf),0,length=10), colors="red");
plot_contourf(pf); colorbar();
axis("equal");
In [ ]: