[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[ff3d-users] Force in a pressure-driven Stokes flow
From: |
Semin |
Subject: |
[ff3d-users] Force in a pressure-driven Stokes flow |
Date: |
Mon, 18 Feb 2008 14:07:05 +0100 |
User-agent: |
Thunderbird 2.0.0.9 (Windows/20071031) |
Dear Stéphane Del Pino,
I would like to compute the hydrodynamic force acting on a motionless
cylinder, which is set in the centre of a rectangular pipe where a
steady Stokes flow is established. As a first step, I performed
computation of the hydrodynamic (pressure + friction forces) forces (per
unit length) acting on an infinite cylinder; in that case the problem is
2D and I successfully used FreeFem++.
In the 3D case and in order to smooth the shape of the cylinder, I added
two half spheres at both ends of the cylinder. I modified the Stokes
file provided in one of your example file (see attached files).
Unfortunately, I get the following problems:
1. I set the pressure to range between -10 and 10: yet, some values of
the pressure are out of this range (it goes from -12 to 12).
2. The pressure gradient and the velocity are in the same direction,
what is not physically correct. Yet, the order of magnitude of the
velocity seems to be ok.
3. The pressure force (even its sign) is consistent with the 2D
computation [the length of the cylinder is approximately 6 ; this helps
to compare the two simulations]. Yet, it is not consistent with the
pressure field. [Please note that I am not sure that my computation of
the surface integral is good.]
4. The order of magnitude and the sign of the friction force are not
correct; its sign is consistent with the velocity field.
5. I try to use the force computation method you proposed in a
ff3d-users mail (see “other evaluation of the total hydrodynamic
force”). It does not work. I do not understand this method, and do not
know if it is relevant here. I would like to know when it is relevant,
and to have a reference about this method.
I use the ff3d-win32-20071108.zip
<http://www.freefem.org/ff3d/binaries/ff3d-win32-20071108.zip> release.
Thanks a lot for helping a young experimenter.
Best regards
Benoît Semin
union {sphere{
<2.5,2,0.5>,0.3 }
cylinder{<2.5,2,0.5>,<2.5,8,0.5>,0.3}
sphere{<2.5,8,0.5>,0.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: stokes.ff,v 1.5 2005/11/28 00:54:32 delpinux Exp $
// Stokes equations
vector n = (40,40,40);
vertex A = (0,0,0);
vertex B = (5,10,1);
double I=0;double J=0; double K=0;
scene S = pov("cylinder1.pov");
mesh M = structured(n,A,B);
domain Omega = domain(S,outside(<1,0,0>));
mesh obstacle = surface(Omega,M);
function ki=one(<1,0,0>);
function teta0 = 0;
double eps = 1e-5;
solve(u1,u2,u3,p) in Omega by M
//memory(matrix=none),
krylov(precond=diagonal,type=bicg),
bicg(maxiter=5000,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)
-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 = 0 on M ymax;
u1 = 0 on M zmin;
u1 = 0 on M zmax;
u1 = 0 on <1,0,0>;
u2 = 0 on M xmin;
u2 = 0 on M xmax;
u2 = 0 on M zmin;
u2 = 0 on M zmax;
u2 = 0 on <1,0,0>;
u3 = 0 on M;
u3 = 0 on <1,0,0>;
p = -10 on M ymin;
p = 10 on M ymax;
};
// Post-processing
mesh T = tetrahedrize(Omega,M);
save(medit,"cylinder1-u",T);
save(medit,"cylinder1-u",[u1, u2, u3],T);
save(medit,"cylinder1-p",T);
save(medit,"cylinder1-p",p,T);
I=int[obstacle](-p*ny+(dx(u2)+dy(u1))*nx+2*dy(u2)*ny+(dz(u2)+dy(u3))*nz);
J=int[obstacle](-p*ny);
K=int[obstacle]((dx(u2)+dy(u1))*nx+2*dy(u2)*ny+(dz(u2)+dy(u3))*nz);
double fy = int[M](ki/eps*u2);
cout<<"forces along the y direction \n";
cout <<"total hydrodynamic force "<<I << "\n" ;
cout <<"pressure force "<<J << "\n" ;
cout <<"friction force "<<K << "\n" ;
cout<<"other estimation of the total hydrodynamic force "<<fy<<"\n";
cout<<"p(2.5,1,0.5)="<<p(2.5,1,0.5)<<"\n";
cout<<"u2(2.5,1,0.5)="<<u2(2.5,1,0.5)<<"\n";
cout<<"u2(3.5,0.5,0.5)="<<u2(3.5,0.5,0.5)<<"\n";
//Computation of the force acting
//on an infinite cylinder with FreeFem++
real n=0.5,L=2.5;
real r=0.3,p=50;
real a=0,b=0,c=0;
func f=2; //pressure gradient
border C0(t=0,2*pi){x=r*cos(t);y=-r*sin(t);}
border C11(t=0,1){x=-L+2*L*t;y=-n;}
border C12(t=0,1){x=L;y=-n+2*n*t;}
border C13(t=0,1){x=L-2*L*t;y=n;}
border C14(t=0,1){x=-L;y=n-2*n*t;}
mesh Th=buildmesh(C0(35*r*p)+C11(L*p)+
C12(n*p)+C13(L*p)+C14(n*p));
//plot(Th,wait=true);
fespace Vh(Th,P2); Vh phi,v;
solve Laplace(phi,v)=int2d(Th)(dx(phi)*dx(v) + dy(phi)*dy(v))
-int2d(Th)(f*v) +
on(C0,C11,C13,phi=0)+on(C12,C14,phi=0);
a=(-int1d(Th,C0)(dx(phi)*N.x+dy(phi)*N.y)
+f*pi*(r)*(r));
b=(f*pi*(r)*(r));
c=(-int1d(Th,C0)(dx(phi)*N.x+dy(phi)*N.y));
cout<<"forces per unit length: "<<endl;
cout << "total force = " << a<<endl;
cout << "pressure force = " << b<<endl;
cout << "friction force = " << c<<endl;
plot(phi);
- [ff3d-users] Force in a pressure-driven Stokes flow,
Semin <=