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
[email protected]
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);
_______________________________________________
ff3d-users mailing list
[email protected]
http://lists.nongnu.org/mailman/listinfo/ff3d-users
