ff3d-users
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [ff3d-users] Navier-Stokes problem


From: Thomas B
Subject: Re: [ff3d-users] Navier-Stokes problem
Date: Fri, 08 Jun 2007 14:17:04 +0200
User-agent: Icedove 1.5.0.10 (X11/20070329)

Hello,

thank you very much for your help.

The simulation (basically the script from Dominique) calculated now some forces, but the results were very strange values...

The force is i.e. always very small (<10^-5) and in directions where it should be positive, it is circa -10^9!

My current goal is to simulate a simple plate in the flow channel. The plate has different angles (10°-80°) and the flow different velocities (2-12), but as I said the dependencies are at the moment very unlogical.

I attached the current ff-Script. Maybe you spot a stupid error?
In "final.pov" you'll find something like
box { <0,0,0>, <1.5,0.3,0.08>
pigment {color rgb<1,0,0>}
rotate<0,-10,0>
translate<2,0.35,1.66074>}


Regards,
Thomas

DELPINO Stephane wrote:
Hello. I just reminded that I made a mistake.

Stephane Del Pino a écrit :

Yes. The force on the object can be approximated by the integral over the object of 1/epsilon*U :
    double fx = int[M](1/epsilon*ux);
    double fy = int[M](1/epsilon*uy);
    double fz = int[M](1/epsilon*uz);
where U = (ux,uy,uz) and the epsilon is used for the penalty term P in the example for Dominique, 1/epsilon = 10^3.
One should have read

double fx = int[M](ki/epsilon*ux);

where ki is a function that is 0 outside the obstacle and 1 inside. For instance:
   function ki = one(<1,0,0>);
where <1,0,0> is the reference of the object.

Regards,
Stephane.



_______________________________________________
ff3d-users mailing list
address@hidden
http://lists.nongnu.org/mailman/listinfo/ff3d-users


// -*- c++ -*-


function uin = -4;
double itsteps = 100;

vector n2 = (60,15,30);
vector n1 = 2*n2;

vertex A = (0,0,0);
vertex B = (4,1,2);
mesh M = structured(n1,A,B);
mesh M2= structured(n2,A,B);
scene S = pov("final.pov");
domain Omega = domain(S,outside(<1,0,0>));

mesh obstacle = surface(Omega,M);

double eps = 1000;
function P = eps*one(<1,0,0>)+1;

double area = int[M](1);

function mu = 0.1;
femfunction u0(M) = 0;
femfunction v0(M) = 0;
femfunction w0(M) = 0;
femfunction p0(M2) = 0;

function ki = one(<1,0,0>); 

double dt = 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 = 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 = 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 = 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 M2
      krylov(type=cg,precond=ichol)
    {

        pde(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;


         // Force
         double fx = int[M](ki/eps*u); 
         double fy = int[M](ki/eps*v);
         double fz = int[M](ki/eps*w);

         cout << "<#i=" << i << "\n";
         cout << "<#Fx=" << fx << "\n";
         cout << "<#Fy=" << fy << "\n";
         cout << "<#Fz=" << fz << "\n";
         
} while(i <= itsteps);


mesh T = tetrahedrize(Omega,M);

save(medit,"velocity",T);
save(medit,"velocity",[u,v,w],T);
save(medit,"pressure",T);
save(medit,"pressure",p0,T);


reply via email to

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