[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [ff3d-users] Stationary Navier-Stokes equations / Flow past cylinder
|
From: |
Dominique Deloison |
|
Subject: |
Re: [ff3d-users] Stationary Navier-Stokes equations / Flow past cylinder problem |
|
Date: |
Wed, 23 May 2007 00:45:54 +0200 |
|
User-agent: |
Mozilla Thunderbird 1.0.6 (Windows/20050716) |
Hi Stéphane,
I thank you again for your help. Unfortunately I still have problems.
I tried to scale down the "navier-stokes" example from 5 m to 5 cm. (see
attached file).
I modified thus:
- the dimensions (fictitious domain and cylinder),
- the inlet velocity and the viscosity so that the fluid goes past the
cylinder within the 5 iterations, while keeping Reynolds number at the
same value.
Unfortunately the results are very different of the original results:
there is no more deviation of the flow by the cylinder.
There are two options, I believe:
1) I 've done a really stupid mistake,
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.
If you could spare some time to give me some advice, I would be very
grateful.
Dominique.
Stephane Del Pino a écrit :
Hi Dominique.
Le lundi 21 mai 2007, Dominique Deloison a écrit :
Hi Stéphane,
Thank you for your quick answer. I have however still some problems.
To answer your question, I found the out-of-date" navier-stokes"
example here: http://www.freefem.org/ff3d/examples.html
Ooops. These are really obsolete. The same applies for the online
documentation...
I have just updated the files.
Concerning the files you sent me, I have the following comments:
1) - during the resolution, the following warning appears:
"implementation not finished! use of dofposition required". Is is a
problem?
This is not important. It is just a reminder for fictitious domain
development...
2) - the size of each element (0.1) is the same than the
radius of the cylinder, so it is very difficult to recognize anything and
to be able to decide whether the solution looks reasonable or not. I guess
I have to increase the number of elements.
Yes. This is just an example of use, it might not be anyway physical...
By the ways I don't know why the cylinder size changed...
I just have been looking to it and it seems to me that it is not a good
example. So I slightly changed the example. You can download it from the web
site...
3) - the syntax of the "convect" operator is different than the
one used in the documentation.
According to the documentation, the syntax is:
"convect(phi,ux,uy,uz)".
In the file there is: "convect([u0,v0,w0],dt,u0)" which is
identical to the freefem++ syntax.
Yes. The online manual was also obsolete. I just put a newer version, sorry
for the confusion...
4) - I have difficulties to figure out the mathematical description
of the problem described in the example file. For me it looks like a
transient problem. I would be interested in a stationnary solution as
described by the equations I wrote in my previous mail. Is there a non
time dependent version of the convect operator ?
Yes. This is a transient resolution...
There is no stationary convect operator, but, you can write a fix point
algorithm. I join you an example.
5) I tried to postprocess the streamlines with medit (I wrote an .iso
file to describe the seeding points). I was unsuccessful (segmentation
fault). have you ever tried this with this example?
The reason is that medit is not very good for hexahedra.
In the joined file I use a tetrahedrization of the mesh and an alternative VTK
output.
I hope that helps. Thanks for your reports.
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: stokes.ff,v 1.5 2005/11/28 00:54:32 delpinux Exp $
// Stokes equations
vector n = (10,10,10);
vertex A = (0,0,0);
vertex B = (1,1,1);
mesh M = structured(n,A,B);
function teta0 = 0;
double eps = 1e-2;
function u1n=0;
function u2n=0;
function u3n=0;
function u1=0;
function u2=0;
function u3=0;
function p =0;
double i=0;
do {
solve(u1,u2,u3,p) in M
krylov(precond=diagonal,type=bicg),
cg(epsilon=1e-5) {
test(U1,U2,U3,P)
int(grad(u1)*grad(U1) - p*dx(U1)
+grad(u2)*grad(U2) - p*dy(U2)
+grad(u3)*grad(U3) - p*dz(U3)
+u1n*dx(u1)*U1+u2n*dx(u2)*U2+u3n*dx(u3)*U3
-eps*grad(p)*grad(P) - eps*p*P
-dx(u1)*P - dy(u2)*P - dz(u3)*P) = 0;
u1 = 0 on M xmin;
u1 = 0 on M xmax;
u1 = 0 on M ymin;
u1 = 1 on M ymax;
u1 = 0 on M zmin;
u1 = 0 on M zmax;
u2 = 0 on M;
u3 = 0 on M;
};
u1n=u1;
u2n=u2;
u3n=u3;
i++;
} while (i<=3);
// Post-processing
mesh T = tetrahedrize(M);
save(medit,"velocity",T);
save(medit,"velocity",[u1, u2, u3],T);
save(medit,"pressure",T);
save(medit,"pressure",p,T);
save(vtk,"ns",{[u1,u2,u3],p},T);
------------------------------------------------------------------------
_______________________________________________
ff3d-users mailing list
address@hidden
http://lists.nongnu.org/mailman/listinfo/ff3d-users
// 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.3 2007/05/21 21:23:17 delpinux Exp $
cylinder {
<5e-3,5e-3,-1>, <5e-3,5e-3>, 2e-3
pigment { color rgb <1,0,0> }
}
// -*- 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 n = (50,10,10);
vertex A = (0,0,0);
vertex B = (5e-2,1e-2,1e-2);
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);
function P = 1E3*one(<1,0,0>)+1;
// Computes the area of Omega
double area = int[M](1);
function uin = 1;
function mu = 1;
femfunction u0(M) = uin;
femfunction v0(M) = 0;
femfunction w0(M) = 0;
femfunction p0(M) = 2 - x;
double dt = 0.01;
double eps = 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 = convect([u0,v0,w0],dt,u0) on M xmax;
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 = convect([u0,v0,w0],dt,v0) on M xmax;
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 = convect([u0,v0,w0],dt,w0) on M xmax;
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 M
krylov(type=cg,precond=ichol)
{
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.
mesh T = tetrahedrize(Omega,M);
save(vtk,"ns2",[u,v,w],M);
save(medit,"velocity",T);
save(medit,"velocity",[u,v,w],T);
save(medit,"pressure",T);
save(medit,"pressure",p0,T);