Stacionarni inkompresibilni Navier-Stokes
Podsjetimo slaba (mješovita) formulacija za Navier-Stokes glasi.
$$ \text{Pronađi } (u,p)\in u_0 +V $$ $$ \displaystyle \forall (v,q)\in V,\quad \int_{\Omega} \sigma(u,p):\varepsilon(v)+Du u \cdot v - q \textrm{div} u\,\mathrm{d}x = \int_{\Gamma_\text{out}} h_0\cdot v\,\mathrm{d}S $$gdje je $\sigma(u,p)=2\mu\varepsilon(u)-pI$, a prostori su:
$V=\{(u,p)\in \mathbf{H}^1(\Omega,\mathbb{R}^d)\times \mathbf{L}^2(\Omega)\;:\; u_{|\Gamma_{\text{in}}\cup\Gamma}=0\},$
$u_0+V=\{(u,p)\in \mathbf{H}^1(\Omega,\mathbb{R}^d)\times \mathbf{L}^2(\Omega)\;:\; u_{|\Gamma_{\text{in}}}=u_0, \: u_{|\Gamma}=0\}.$
Odabir konačnih elemenata u rješavanju inkompresiblnih modela fluida.
Koristimo $\mathbb{P}1-\mathbb{P}1$ sa stabilizacijskim članom:
$$ \delta\, h^2\, \int_{\Omega} \nabla u \cdot \nabla v \,\mathrm{d}x $$Alternativno je potrebno koristi Taylor-Hood elemente koji zadovoljavaju inf-sup uvjet ili Ladyzhenskaye-Babuška-Brezzi uvjet, tipa $\mathbb{P}2-\mathbb{P}1$. Ovo je korisno ako nam je potrebna dodatna regularnost brzine kod korištenja derivacije oblika dan preko integrala po rubu od $\Omega$.
Za tretiranje tlaka (obično tlak $p$ pripada prostoru $\mathbf{L}_0^2(\Omega)=\{ p\in\mathbf{L}^2(\Omega)\: :\:\int_\Omega p =0 \}$), jedan od jednostavnih rješenja je kroz penalizirajući element $$ -\epsilon \int_\Omega p q , \qquad \epsilon>0. $$
// stabilizacija penaliziranjem za P1/P1 element kod Stokesa
real delta = 0.01; // za P1/P1 elemente, ako se koriste P2/P1 stavi delta=0
real epsln = 1.0e-6;
problem Stokes(ux,uy,p, vx,vy,q)=
int2d(Omega)(2*mu*(dx(ux)*dx(vx) + dy(uy)*dy(vy) + 0.5*(dx(uy) + dy(ux))*(dx(vy) + dy(vx)))
-p*(dx(vx)+dy(vy))
-q*(dx(ux)+dy(uy))
- delta * hTriangle * hTriangle * (dx(p)*dx(q) + dy(p)*dy(q))
- p*q*epsln
)
+on(1, ux=u0x, uy=u0y)
+on(0, ux=0.0, uy=0.0)
+on(3, ux=0.0, uy=0.0);
problem LinearNavierStokes(dux,duy,dp, vx,vy,q)=
int2d(Omega)(2*mu*(dx(dux)*dx(vx) + dy(duy)*dy(vy) + 0.5*(dx(duy) + dy(dux))*(dx(vy) + dy(vx)))
+ux*dx(dux)*vx + ux*dx(duy)*vy + uy*dy(dux)*vx +uy*dy(duy)*vy
+dx(ux)*dux*vx + dx(uy)*dux*vy + dy(ux)*duy*vx + dy(uy)*duy*vy
-dp*(dx(vx)+dy(vy))
-q*(dx(dux)+dy(duy))
- delta * hTriangle * hTriangle * (dx(dp)*dx(q) + dy(dp)*dy(q))
- epsln*dp*q
)
+int2d(Omega)(2*mu*(dx(ux)*dx(vx) + dy(uy)*dy(vy) + 0.5*(dx(uy) + dy(ux))*(dx(vy) + dy(vx)))
+ux*dx(ux)*vx + ux*dx(uy)*vy + uy*dy(ux)*vx +uy*dy(uy)*vy
-p*(dx(vx)+dy(vy))
-q*(dx(ux)+dy(uy))
- delta * hTriangle * hTriangle * (dx(p)*dx(q) + dy(p)*dy(q))
- p*q*epsln
)
+on(1, dux=0.0, duy=0.0)
+on(0, dux=0.0, duy=0.0)
+on(3, dux=0.0, duy=0.0);
...
for (int iter=0; iterMaxIter; ++iter)
{
...
Stokes;
for (int i=0; i<20; i++){
LinearNavierStokes;
ux=ux+dux;
uy=uy+duy;
p=p+dp;
real Error=sqrt(int2d(Omega)((dux)^2))/sqrt(int2d(Omega)((ux)^2))
+sqrt(int2d(Omega)((duy)^2))/sqrt(int2d(Omega)((uy)^2))
+sqrt(int2d(Omega)((dp)^2))/sqrt(int2d(Omega)((p)^2));
if (Error<1e-5) break;
...
}
Proširena Lagrangeova metoda
Problem tipično sadrži neki uvjet na volumen ili perimetar u obliku $G(\Omega)=0$. Umjesto problema uvjetne optimizacije: $$ \min_{\Omega,\: G(\Omega)=0} J(\Omega),$$ numerički tretiramo bezuvjetni problem: $$ \min_{\Omega} J(\Omega)+\lambda G(\Omega)+\frac{a}{2}G(\Omega)^2, $$ gdje je $\lambda$ Lagrangeov multiplikator, a $a>0$ penalizirajući faktor.
Pravilo generiranja $(\lambda^{k+1},a^{k+1})$ za dani $(\Omega^k,\lambda^k,a^k)$: $$ \lambda^{k+1}=\lambda^k+a^k G(\Omega^k), $$ $$ a^{k+1} = \left\{ \begin{array}{lr} \alpha a^{k}, & \text{$k$ manji od $k_0$} \\a^k, &\text{ inače. }\end{array} \right. $$
G=Omega.area-vol0;
if (iter==0) {lagr=30; a=1;}
lagr=lagr+10*a*G;
a=(iter<100)?1.005*a:a;
a=(a<1e3)?a:1e3;
Rezultati simulacija
