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 [email protected] http://lists.nongnu.org/mailman/listinfo/ff3d-users
