Hello.
I finally had the time to check the navier-stokes projection algorithm from
J-L Guermond. I joined you the file. In my opinion, it gives much better
results.
Best regards,
Stéphane.
// -*- c++ -*-
function uin = -4;
double itsteps = 100;
vector n2 = (60,15,30);
vector n1 = 2*n2;
vertex A = (0,0,0);
vertex B = (4,1,2);
mesh M = structured(n1,A,B);
mesh M2= structured(n2,A,B);
scene S = pov("final.pov");
domain Omega = domain(S,outside(<1,0,0>));
mesh obstacle = surface(Omega,M);
save (vtk,"o",obstacle);
double eps = 1E-3;
function P = 1/eps*one(<1,0,0>)+1;
double area = int[M](1);
function mu = 0.1;
femfunction u0(M) = 0;
femfunction v0(M) = 0;
femfunction w0(M) = 0;
femfunction p0(M2) = 0;
function ki = one(<1,0,0>);
double dt = 0.01;
double pp;
double qq;
double i = 1;
function u = 0;
function v = 0;
function w = 0;
function p1 = 0;
do{
function deltaP = 2*p0-p1;
cout << "========== Navier-Stokes step " << i << "\n";
// Reinitialise pour le gradient conjugue
u = 0;
v = 0;
w = 0;
solve(u) in M
krylov(type=cg,precond=diagonal),
memory(matrix=none)
{
test(U)
int(P/dt*u*U)
+int(mu*grad(u)*grad(U))
=
int(1/dt*convect([u0,v0,w0],dt,u0)*U)
-int(dx(deltaP)*U)
;
u = uin on 0;// M xmin;
u = 0 on 2,3,4,5; // other boundaries but xmax
};
solve(v) in M
krylov(type=cg,precond=diagonal),
memory(matrix=none)
{
test(V)
int(P/dt*v*V)
+int(mu*grad(v)*grad(V))
=
int(1/dt*convect([u0,v0,w0],dt,v0)*V)
-int(dy(deltaP)*V)
;
v = 0 on 0,2,3,4,5;
};
solve(w) in M
krylov(type=cg,precond=diagonal),
memory(matrix=none)
{
test(W)
int(P/dt*w*W)
+int(mu*grad(w)*grad(W))
=
int(1/dt*convect([u0,v0,w0],dt,w0)*W)
-int(dz(deltaP)*W)
;
w = 0 on 0,2,3,4,5;
};
solve(p) in M2
krylov(precond=ichol)
{
test(q)
int(grad(p)*grad(q))
= int(grad(p0)*grad(q))
- int(1/dt*dx(u)*q)
- int(1/dt*dy(v)*q)
- int(1/dt*dz(w)*q);
p = 0 on 1;
}
u0 = u;
v0 = v;
w0 = w;
p1 = p0;
p0 = p;
if (i<10) {
save(vtk,"ns.0".i,{[u,v,w],p},M);
} else {
save(vtk,"ns.".i,{[u,v,w],p},M);
}
i = i+1;
double fx = int[M](ki/eps*u);
double fy = int[M](ki/eps*v);
double fz = int[M](ki/eps*w);
cout << "<#i=" << i << "\n";
cout << "<#Fx=" << fx << "\n";
cout << "<#Fy=" << fy << "\n";
cout << "<#Fz=" << fz << "\n";
} while(i <= itsteps);
mesh T = tetrahedrize(Omega,M);
save(medit,"velocity",T);
save(medit,"velocity",[u,v,w],T);
save(medit,"pressure",T);
save(medit,"pressure",p0,T);
_______________________________________________
ff3d-users mailing list
[email protected]
http://lists.nongnu.org/mailman/listinfo/ff3d-users
