Hello.
Le mercredi 23 mai 2007, Dominique Deloison a écrit :
Hi Stéphane,
...
2) I should change something else (P,p0 or eps or dt ...). I have to say
I am not at all familiar with this resolution scheme and that I do not
really understand it. It is clear for me that the dt should be small
because it defines the time discretisation. The role of eps, P and the
definition of p0 are far from being obvious to me.
In fact this example is not very good. It is just here to have an idea of how
ff3d can be used to write projection algorithms...
The main difficulty with this example is the treatment of the pressure. The
pressure is treated using a penalty trick to avoid spurious modes. If I am
right (but this should be checked) this epsilon must be of the form C*h,
where C is a constant and h the mesh step size.
P is important : it is the penalty term that forces the velocity to be 0 in
the cylinder.
I modified again the algorithm using a stable element : so called Q2-isoQ1/Q1
which means that the pressure grid is twice coarser than the velocity grid.
This does not require anymore 'eps'. I join you a modified version where I
also initialized p0 to 0 (probably more reasonable).
Just tell me if you get better results. This is not the best projection
algorithm. You can for example have a look to the one proposed by Jean-Luc
Guermond.
Best regards,
Stéphane.
------------------------------------------------------------------------
// -*- 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.6 2007/05/21 21:23:17 delpinux Exp $
// Navier-Stokes flow around a cylinder
// Uses a projection scheme over the velocity divergence free space.
// Fictitious domain method uses volume penalty method
// Only 5 iterations are performed
vector n2 = (25, 5, 5);
vector n1 = 2*n2;
vertex A = (0,0,0);
vertex B = (5,1,1);
mesh M = structured(n1,A,B);
mesh M2= structured(n2,A,B);
scene S = pov("navier-stokes.pov");
domain Omega = domain(S,outside(<1,0,0>));
mesh obstacle = surface(Omega,M);
function P = 1E3*one(<1,0,0>)+1;
// Computes the area of Omega
double area = int[M](1);
function uin = 12.5;
function mu = 0.1;
femfunction u0(M) = 0;
femfunction v0(M) = 0;
femfunction w0(M) = 0;
femfunction p0(M2) = 0;
double dt = 0.01;
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 M
krylov(type=cg,precond=ichol)
{
pde(u)
(P/dt)*u - div(mu*grad(u))
= ((1/dt)*convect([u0,v0,w0],dt,u0) - dx(p0));
u = uin on M xmin;
u = 0 on M ymin;
u = 0 on M ymax;
u = 0 on M zmin;
u = 0 on M zmax;
};
solve(v) in M
krylov(type=cg,precond=ichol)
{
pde(v)
(P/dt)*v - div(mu*grad(v))
= ((1/dt)*convect([u0,v0,w0],dt,v0) - dy(p0));
v = 0 on M xmin;
v = 0 on M ymin;
v = 0 on M ymax;
v = 0 on M zmin;
v = 0 on M zmax;
};
solve(w) in M
krylov(type=cg,precond=ichol)
{
pde(w)
(P/dt)*w - div(mu*grad(w))
= ((1/dt)*convect([u0,v0,w0],dt,w0) - dz(p0));
w = 0 on M xmin;
w = 0 on M ymin;
w = 0 on M ymax;
w = 0 on M zmin;
w = 0 on M zmax;
};
qq = int[M](dx(u)+dy(v)+dz(w))/area;
solve(q) in M2
krylov(type=cg,precond=ichol)
{
pde(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 <= 10);
// Final results.
mesh T = tetrahedrize(Omega,M);
save(vtk,"ns",[u,v,w],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
address@hidden
http://lists.nongnu.org/mailman/listinfo/ff3d-users