[ff3d-users] Force in a pressure-driven Stokes flow

[Semin]

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