ff3d-users
[Top][All Lists]

## 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
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)
= ((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)
= ((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)
= ((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