[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 (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 }
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
//  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
    bicg(maxiter=5000,epsilon=1e-5) {

          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",[u1, u2, u3],T);


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


//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)+


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


cout<<"forces per unit length: "<<endl;
cout << "total force = " << a<<endl;
cout << "pressure force = " << b<<endl;
cout << "friction force = " << c<<endl;

reply via email to

[Prev in Thread] Current Thread [Next in Thread]