ff3d-users
[Top][All Lists]

## [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
//  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)
-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;

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);
```