Hi.
where did you get the Navier-Stokes example ? It seems that the one in the cvs
tree is correct. I attached to you the files. Let me known if you experience
more troubles.
Best regards,
Stephane.
------------------------------------------------------------------------
// -*- c++ -*-
// This file is part of ff3d - http://www.freefem.org/ff3d
// Copyright (C) 2003 Laboratoire J.-L. Lions UPMC Paris
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software Foundation,
// Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
// $Id: navier-stokes.ff,v 1.5 2007/02/14 21:34:23 delpinux Exp $
// Navier-Stokes flow around a cylinder
// Uses a projection scheme over the velocity divergence free space.
// Fictitious domain method uses the poor elimination method to run
// faster in this example.
// Only 5 iterations are performed
vector n = (50,10,10);
vertex A = (0,0,0);
vertex B = (5,1,1);
mesh M = structured(n,A,B);
scene S = pov("navier-stokes.pov");
domain Omega = domain(S,outside(<1,0,0>));
mesh obstacle = surface(Omega,M);
femfunction Ki(M) = 1 - one(<1,0,0>);
// Computes the area of Omega
double area = int[M](Ki);
// Cylinder diameter is 0.2 => Re = 400
function uin = 12.5;
function mu = 0.005;
femfunction u0(M) = uin;
femfunction v0(M) = 0;
femfunction w0(M) = 0;
femfunction p0(M) = 2 - x;
double dt = 0.01;
double eps = 0.001;
double pp;
double qq;
double i = 1;
function u=0;
function v=0;
function w=0;
do{
cout << "========== Navier-Stokes step " << i << "\n";
solve(u) in Omega by M memory(matrix=none),method(type=eliminate),
krylov(type=cg,precond=diagonal){
pde(u)
(1/dt)*u - div(mu*grad(u))
= ((1/dt)*convect([u0,v0,w0],dt,u0) - dx(p0));
u = 0 on obstacle;
u = uin on M xmin;
u = convect([u0,v0,w0],dt,u0) on M xmax;
};
solve(v) in Omega by M memory(matrix=none),method(type=eliminate),
krylov(type=cg,precond= diagonal){
pde(v)
(1/dt)*v - div(mu*grad(v))
= ((1/dt)*convect([u0,v0,w0],dt,v0) - dy(p0));
v = 0 on obstacle;
v = 0 on M xmin;
v = convect([u0,v0,w0],dt,v0) on M xmax;
v = 0 on M ymin;
v = 0 on M ymax;
};
solve(w) in Omega by M memory(matrix=none),method(type=eliminate),
krylov(type=cg,precond=diagonal){
pde(w)
(1/dt)*w - div(mu*grad(w))
= ((1/dt)*convect([u0,v0,w0],dt,w0) - dz(p0));
w = 0 on obstacle;
w = 0 on M xmin;
w = convect([u0,v0,w0],dt,w0) on M xmax;
w = 0 on M zmin;
w = 0 on M zmax;
};
qq = int[M](dx(u)+dy(v)+dz(w))/area;
solve(q) in Omega by M memory(matrix=none),method(type=eliminate),
krylov(type=cg,precond=diagonal){
pde(q)
(eps*dt)*q - div(dt*grad(q)) = (dx(u)+dy(v)+dz(w) - qq);
q = 0 on M xmax;
};
p0 = p0 - q;
pp = int[M](p0)/area;
p0 = p0 - pp;
u0 = u + (dt*dx(q));
v0 = v + (dt*dy(q));
w0 = w + (dt*dz(q));
i = i+1;
} while(i <= 5);
// Final results.
save(medit,"velocity",M);
save(medit,"velocity",[u,v,w],M);
save(medit,"pressure",M);
save(medit,"pressure",p0,M);
------------------------------------------------------------------------
// This file is part of ff3d - http://www.freefem.org/ff3d
// Copyright (C) 2003 Laboratoire J.-L. Lions UPMC Paris
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software Foundation,
// Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
// $Id: navier-stokes.pov,v 1.2 2004/09/01 21:22:58 delpinux Exp $
cylinder {
<0.5,0.5,0>, <0.5,0.5,1>, 0.1
pigment { color rgb <1,0,0> }
}
------------------------------------------------------------------------
_______________________________________________
ff3d-users mailing list
address@hidden
http://lists.nongnu.org/mailman/listinfo/ff3d-users