Řešíme stacionární laminární proudění nestlačitelné tekutiny s konstantní hustotou popsané systémem Navierových-Stokesových rovnic \begin{eqnarray*} \nabla\cdot\left(\vec{u}\otimes\vec{u}\right) - \nabla\cdot(\nu\nabla \vec{u}) &=& - \nabla p, \\ \nabla\cdot\vec{u} &=& 0, \end{eqnarray*}
$\vec{u}$ je rychlost, $p$ je podíl tlaku a hustoty (někdy nazývaný kinematický tlak) a $\nu$ je (konstantní) kinematická vazkost.
$f$ probíhá přes všechny stěny buňky $\Omega_c$, $\vec{u}_f$ je aproximace $\vec{u}$ v těžišti stěny $f$. Označíme $$ \phi_f := \vec{u}_f\cdot\vec{S}_f. $$
potom aproximace rovnice kontinuity v buňce $c$ je $$ \nabla\cdot\vec{u} \approx \frac{1}{|\Omega_c|}\sum_f \phi_f. $$
Na kartézské síti s krokem $h=\Delta x = \Delta y$ označme sousedy buňky $C$ indexy $N$ - soused nahoře, $W$ - soused vlevo, $S$ - soused dole a $E$ soused vpravo. Potom je (pro kartézskou síť) ve vnitřních buňkách sítě $$ \nabla\cdot(\nu \nabla\vec{u}) \approx \nu \frac{\vec{u}_N + \vec{u}_W + \vec{u}_S + \vec{u}_E - 4\vec{u}_C}{h^2}. $$
Diskrétní verzi operátoru $\nu\nabla\cdot\cdot$ je tedy možné reprezentovat (na kartézské síti a při správném očíslování neznámých) pětidiagonální neostře diagonálně dominantní maticí. Při vhodně implementovaných okrajových odmínkách je tato matice symetrická.
Konvektivní člen je nelineární a proto je třeba provést linearizaci. Označme horním indexem $n$ novéhodnoty a horním indexem $o$ stávající známe hodnoty rychlosti. Linearizaci lze provést různým způsobem, pro algoritmus SIMPLE se používá tzv. Picardova aproximace
$$ \nabla\cdot(\vec{u}^n \otimes \vec{u}^n) \approx \nabla\cdot(\vec{u}^o \otimes \vec{u}^n) $$Potom aproximace tohoto členu pomocí MKO je
$$ \iint_{\Omega_C} \nabla\cdot(\vec{u}^o\otimes\vec{u}^n)\,dV = \sum_{f} \int_{\Gamma_f} \vec{n}\cdot\vec{u}^o\otimes\vec{u}^n\,dS \approx \sum_{f} \phi^o_f \vec{u}^n_f. $$Tedy $$ \nabla\cdot(\vec{u}^o\otimes\vec{u}^n) \approx \frac{1}{|\Omega_C|}\sum_{f} \phi^o_f \vec{u}^n_f. $$
Použijeme-li schéma upwind (např. pro tok mezi buňkami $C$ a $N$) $$ \phi_f \vec{u}_f = (\phi_f)^+ \vec{u}_C + (\phi_f)^- \vec{u}_N, $$ dostáváme
$$ \nabla\cdot(\vec{u}^o\otimes\vec{u}^n) \approx \frac{1}{|\Omega_C|}\left( \sum_{f} (\phi^o_f)^+ \vec{u}^n_C + \sum_{f} (\phi^o_f)^- \vec{u}^n_F\right), $$kde $\vec{u}^n_F$ je hodnota řešení v buňce sousedící s $C$ stěnou $f$.
Diskrétní verze konvektivního členu je tedy reprezentovaná maticí s diagonálním prvkem $$ \frac{1}{\Omega_C}\sum_{f} (\phi^o_f)^+. $$ Mimo diagonálu jsou na odpovídajících místech (sloupce $N$, $E$, $S$ a $W$) příslušné hodnoty $(\phi_f)^-$.
Pokud je splněna rovnice kontinuity, tak je výsledná matice (při správně implementovaných okrajových podmínkách) nesymetrická a neostře diagonálně dominantní.
Algoritmus existuje v několika variantách (např. řešení rovnice pro tlak nebo pouze pro korekci tlaku, detaily diskretizace, ...). Na tomto místě popíšeme variantu používanou v sofwarovém balíku OpenFOAM.
Algoritmus je iterační a používá jako hlavní proměnné $\vec{u}$ a $p$ definované ve středech buňek a jako pomocnou proměnnou $\phi$ definovanou na stěnách.
Na začátku výpočtu se určí $\phi$ na stěně $f$ jako $\phi_f = \vec{u}_f\cdot\vec{S}_f$ kde $\vec{u}_f = (\vec{u}_C + \vec{u}_F)/2$ je lineární interpolací rychlosti do středu stěny. V dalších krocích se $\phi$ počítá jinak.
Nejprve diskretizujeme rovnici pro rychlost (gradient tlaku zatím ponecháme formálně zapsany jako $\nabla p$): $$ a_C^o \vec{u}_C^* = \sum_{f} a_F^o \vec{u}^*_F + Q_C - \nabla p_C^o = H(\vec{u}^*)_C - \nabla p_C^o. $$ Protože byl použit starý tlak $p^o$, nedostaneme řešením této rovnice $\vec{u}^n$, ale pouze odhad $\vec{u}^*$. Ten však nemusí splňovat rovnici kontinuity. Vyjádříme proto $$ \vec{u}_C^* =\frac{1}{a_C^o} H(\vec{u}^*)_C - \frac{1}{a_C^o}\nabla p_C^o = \hat{u}_C^* - \frac{1}{a_C^o}\nabla p_C^o. $$ kde $\hat{u}_C^* := \frac{1}{a_C^o} H(\vec{u}^*)_C$.
Tento odhad rychlost interpolujeme na stěnu $f$ $$ \vec{u}_f^* = \hat{u}_f^* - \frac{1}{a_f^o}\nabla p_f^o. $$
Pro $\vec{u}^n$ pak požadujeme splnění rovnice (pozor na horní indexy $n$, $o$ a $*$)> $$ \vec{u}_f^n = \hat{u}_f^* - \frac{1}{a_f^o}\nabla p_f^n. $$
Tato rychlost by však měla splňovat rovnici kontinuity, takže $$ 0 = \sum_f \phi_f^n = \sum_f \vec{u}_f^n \cdot \vec{S}_f = \sum_f \hat{u}_f^* \cdot \vec{S}_f - \sum_f \frac{1}{a_f^o}\nabla p_f^n \cdot \vec{S}_f. $$
Diskretizovaná rovnice pro tlak má tedy tvar $$ \sum_f \frac{1}{a_f^o}\nabla p_f^n \cdot \vec{S}_f = \sum_f \hat{u}_f^* \cdot \vec{S}_f, $$ což je diskrétní verze rovnice $\nabla\cdot(\frac{1}{a}\nabla p) = \nabla\cdot \hat{u}^*$. Z nové hodnoty tlaku $p^n$ určíme novou hodnotu rychlosti a objemového toku $\phi$ \begin{eqnarray} \vec{u}^n_C &=& \hat{u}_C^* - \frac{1}{a_C^o}\nabla p_C^n, \\ \phi_f^n &=& \hat{u}_f^*\cdot \vec{S}_f - \frac{1}{a_f^o}\nabla p_f^n\cdot\vec{S}_f. \end{eqnarray}
Jedna iterace algoritmu SIMPLE jedna iterace algoritmu SIMPLE zahrnující relaxaci v rovnice pro rychlost a explicitní relaxaci tlaku tedy vypadá následovně:
1) ze známých hodnot $\vec{u}^o$ a $\phi^o$ sestavíme aproximaci levé strany rovnice pro rychlost (tj. bez vlivu gradientu tlaku) $$ \nabla\cdot\left(\vec{u}^o\otimes\vec{u}^*\right) - \nabla\cdot(\nu\nabla \vec{u}^*) \approx a^o_C \vec{u}_C - H(\vec{u}^*)_C, $$
2) provedeme relaxaci s koeficientem $\alpha$
3) sestavíme a vyřešíme rovnici pro odhad rychlosti $\vec{u}^*$ se gradientem tlaku počítaným z $p^o$ $$ a^o_C \vec{u}^*_C = H(\vec{u}^*)_C - \nabla p^o_C, $$
4) pomocí lineární interpolace určíme $\hat{u}_f$ ($F$ je index sousední buňky) $$ \hat{u}_f^* = \text{interpolate}\left(\frac{1}{a^o_C}H(\vec{u}^*)_C, \frac{1}{a^o_F}H(\vec{u}^*)_F\right). $$
5) sestavíme a vyřešíme rovnici pro tlak $p^*$ $$ \sum_f \frac{1}{a_f^o}\nabla p_f^* \cdot \vec{S}_f = \sum_f \hat{u}_f^* \cdot \vec{S}_f, $$
6) nový tlak určíme z $p^o$ a $p^*$ pomocí relaxace s koeficientem $\beta$ $$ p^n = p^o + \beta(p^* - p^o) $$
7) určíme nové hodnoty rychlosti a objemového toku \begin{eqnarray} \vec{u}^n_C &=& \hat{u}_C^* - \frac{1}{a_C^o}\nabla p_C^n, \\ \phi_f^n &=& \hat{u}_f^*\cdot \vec{S}_f - \frac{1}{a_f^o}\nabla p_f^n\cdot\vec{S}_f. \end{eqnarray}
Výše uvedený postup opakujeme dokud proces nezkonverguje ke stacionárnímu stavu.
In [1]:
using PyPlot;
In [2]:
include("mesh.jl");
include("fields.jl");
include("operators.jl");
Vazké členy jsme již implementovali v minulé přednášce. Zbývá pouze konvektivní člen. Ten zapíšeme jako div(phi, U)
In [3]:
function div(ϕ::ScalarList , U::Field{T}) where {T}
@assert length(ϕ) == length(U.mesh.owner)
A = spzeros(length(U.values),length(U.values))
b = zeros(T, length(U.values))
mesh = U.mesh
for f in internal_faces(mesh)
o = mesh.owner[f]
n = mesh.neighbor[f]
α = max(ϕ[f], 0.0)
β = min(ϕ[f], 0.0)
A[o,o] += α / mesh.volume[o]
A[o,n] += β / mesh.volume[o]
A[n,o] -= α / mesh.volume[n]
A[n,n] -= β / mesh.volume[n]
end
for p in boundary_patches(mesh)
name = mesh.patch[p].name
bc = U.boundaries[name]
for f in patch_faces(mesh, p)
o = mesh.owner[f]
c1, c2 = boundary_coeffs(bc, f)
α = max(ϕ[f], 0.0)
β = min(ϕ[f], 0.0)
A[o,o] += (α + β*c2) / mesh.volume[o]
b[o] += β*c1 / mesh.volume[o]
end
end
return Equation{T}(A, U.values, b)
end
Out[3]:
Při spuštění výpočtu musíme vypočítat $\phi_f$ pomocí lineární interpolace
In [4]:
function create_ϕ_by_interpolation(U::VectorField)
mesh = U.mesh
ϕ = zeros(Scalar, length(all_faces(mesh)))
for f in internal_faces(mesh)
o = mesh.owner[f]
n = mesh.neighbor[f]
ϕ[f] = dot(mesh.surface[f], (U[o] + U[n]) / 2)
end
for p in boundary_patches(mesh)
name = mesh.patch[p].name
bc = U.boundaries[name]
for f in patch_faces(mesh, p)
o = mesh.owner[f]
Uf = boundary_value(bc, f)
ϕ[f] = dot(mesh.surface[f], Uf)
end
end
return ϕ
end
Out[4]:
V průběhu výpočtu budeme $\phi_f$ počítat z $\hat{u}$, $p$ a $1/a_f$. Máme $$ \phi^n_f = \hat{u}_f\cdot\vec{S}_f - \frac{1}{a_f} \nabla p_f\cdot\vec{S}_f = \hat{u}_f\cdot\vec{S}_f - \frac{1}{a_f} \frac{\partial p}{\partial n} ||\vec{S}_f||. $$
In [5]:
function calculate_ϕ!(ϕ::ScalarList, Ubar::VectorField, p::ScalarField, ra::ScalarList)
mesh = Ubar.mesh
for f in internal_faces(mesh)
o = mesh.owner[f]
n = mesh.neighbor[f]
S = mesh.surface[f]
δ = norm(mesh.centre[n] - mesh.centre[o])
pn = (p[n] - p[o]) / δ
ϕ[f] = dot(S, (Ubar[o]+Ubar[n])/2.) - (ra[o]+ra[n])/2 * pn * norm(S)
end
for pa in boundary_patches(mesh)
name = mesh.patch[pa].name
bc = Ubar.boundaries[name]
pbc = p.boundaries[name]
for f in patch_faces(mesh, pa)
o = mesh.owner[f]
S = mesh.surface[f]
Ub = boundary_value(bc, f)
pb = boundary_value(pbc, f)
pn = (pb - p[o]) / norm(mesh.centre[o] - mesh.facecentre[f])
ϕ[f] = dot(S, Ub) - ra[o] * pn
end
end
end
Out[5]:
In [6]:
function SIMPLE(U, p, ν; fix_pressure=true, iters=50)
mesh = U.mesh
α = 0.7
β = 0.3
ϕ = create_ϕ_by_interpolation(U)
for iter = 0:iters
UOld, pOld = copy(U.values), copy(p.values)
UEqn = div(ϕ,U) - Δ(ν,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);
if fix_pressure
pEqn.A[1,1] -= 1/ν #length(mesh.centre)
end
solve!(pEqn)
p ← β*p + (1-β)*pOld
U ← Ubar - ra .* grad(p)
calculate_ϕ!(ϕ, Ubar, p, ra)
if rem(iter,5)==0
nxny = length(mesh.centre)
pRez = norm(pOld - p.values) / nxny
URez = norm(UOld - U.values) / nxny
println(iter, "\t", pRez, "\t", URez)
end
end
return U,p
end
Out[6]:
In [7]:
function create_fields(mesh)
U = VectorField(mesh)
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 = ScalarField(mesh);
for name ∈ ["left", "right", "bottom", "top"]
set_neumann_patch!(p, name)
end
return (U,p)
end;
In [8]:
function cavity_ns(ν, mesh)
U,p = create_fields(mesh);
U,p = SIMPLE(U, p, ν)
return U,p
end
Out[8]:
In [9]:
ν = 0.01;
msh25 = cartesian_mesh(25,25);
U,p = cavity_ns(ν, msh25);
In [10]:
plot_contourf(p; cmap="jet"); colorbar();
plot_arrows(U);
axis("equal");
In [11]:
ν = 1.e-3;
msh50 = cartesian_mesh(50,50);
In [12]:
U,p = cavity_ns(ν, msh50);
In [13]:
plot_contourf(p; cmap="jet"); colorbar();
plot_arrows(U);
axis("equal");
In [14]:
component(f::VectorField, c) = [ v[c] for v in f.values ]
n=50
uu = reshape(component(U,1), (n,n));
vv = reshape(component(U,2), (n,n));
pp = reshape(p.values, (n,n));
In [15]:
# Ψ_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,1]+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
#Ψ[end,end] = (Ψ[end-1,end]+Ψ[end,end-1]) / 2
return Ψ
end
Out[15]:
In [16]:
Ψ = streamfunction(uu,vv,1/50,1/50);
pos = range(0, maximum(Ψ), length=10)
contour(Ψ', pos, colors="black");
neg = range(minimum(Ψ), 0, length=5)
contourf(pp'; cmap="jet"); colorbar();
contour(Ψ', neg, colors="red");
axis("equal");
Řešíme problém ustáleného obtékání eliptického tělesa v kanále nevazkou nestlačitelnou tekutinou popsanou systémem Navierových-Stokesových rovnic. Při řešení uvažujeme následující okrajové podmínky:
obrázek oblasti viz níže
In [17]:
include("gmsh_mesh.jl")
Out[17]:
In [18]:
gmsh = gmsh_mesh("domain.msh");
In [19]:
plot_mesh(gmsh)
axis("equal");
In [20]:
patch_names(gmsh)
Out[20]:
In [21]:
U = VectorField(gmsh)
set_dirichlet_patch!(U, "INLET", Vector(1,0));
set_dirichlet_patch!(U, "TOP", Vector(0,0));
set_dirichlet_patch!(U, "BOTTOM", Vector(0,0));
set_dirichlet_patch!(U, "PROFILE", Vector(0,0));
set_neumann_patch!(U, "OUTLET");
U.values[:] = [Vector(1,0) for c in cells(gmsh)]
p = ScalarField(gmsh);
for name ∈ ["INLET", "PROFILE", "TOP", "BOTTOM"]
set_neumann_patch!(p, name)
end
set_dirichlet_patch!(p, "OUTLET", 0.0);
In [22]:
U,p = SIMPLE(U, p, 1.e-3, fix_pressure=false, iters=100);
In [23]:
plot_contourf(U; cmap="jet"); colorbar();
plot_arrows(U);
axis("equal");
In [24]:
gmf = gmsh_mesh("domain_bl.msh");
In [25]:
plot_mesh(gmf);
axis("equal");
xlim(-1.5,0); ylim(-0.1,1.1);
In [26]:
length(cells(gmf))
Out[26]:
In [27]:
patch_names(gmf)
Out[27]:
In [28]:
Uf = VectorField(gmf)
set_dirichlet_patch!(Uf, "INLET", Vector(1,0));
set_dirichlet_patch!(Uf, "TOP", Vector(0,0));
set_dirichlet_patch!(Uf, "BOTTOM", Vector(0,0));
set_dirichlet_patch!(Uf, "PROFILE", Vector(0,0));
set_neumann_patch!(Uf, "OUTLET");
Uf ← [Vector(1,0) for c in cells(gmf)]
pf = ScalarField(gmf);
for name ∈ ["INLET", "PROFILE", "TOP", "BOTTOM"]
set_neumann_patch!(pf, name)
end
set_dirichlet_patch!(pf, "OUTLET", 0.0);
In [29]:
Uf,pf = SIMPLE(Uf, pf, 1.e-2, fix_pressure=false, iters=10);
In [30]:
Uf,pf = SIMPLE(Uf, pf, 1.e-3, fix_pressure=false, iters=50);
In [31]:
plot_contourf(Uf; cmap="jet"); colorbar();
#plot_arrows(Uf);
axis("equal");
In [32]:
plot_contourf(Uf; cmap="jet"); colorbar();
plot_arrows(Uf; units="xy", scale=10);
axis("equal");
xlim(-1.5,1.5); ylim(-2.5,2.5);
In [ ]: