Discussion:
[ff3d-users] help!
Augusto Almeida
2007-07-03 21:54:12 UTC
Permalink
I need run a enviroment (a quadratic tube), with a viscous fluid in there,
the aim is steady state and export the vector components (discretization
70x70x70) in a archive with this format :

vx vy vz

The archive contais 343000 registers.

Follow the code (sorry some parts is in portuguese!) :

Thanks too much!

double limite_inf = 0.0;
double limite_inf_r = 0.0;
double limite_sup = 1.0;
double limite_sup_r = 1.0;
double dimensao = 20.0;
double dimensao_refinada = 70.0;

// Definindo a dimensao do solido (deve ser menor para evitar calculos
// desnecessarios pois dentro do solido a velocidade e zero)
vector normal = (dimensao, dimensao, dimensao);

// Definindo a dimensao da malha onde o fluido percorre
vector refinado = (dimensao_refinada, dimensao_refinada, dimensao_refinada);

// Definindo os limites superiores e inferiores da "caixa" onde ocorre a
// simulacao
vertex ini = (limite_inf , limite_inf , limite_inf);
vertex ini_r = (limite_inf_r, limite_inf_r, limite_inf_r);
vertex fin = (limite_sup , limite_sup , limite_sup);
vertex fin_r = (limite_sup_r, limite_sup_r, limite_sup_r);

// Criando as malhas
mesh MALHA = structured(normal , ini, fin);
mesh MALHA_REFINADA = structured(refinado, ini_r, fin_r);

// Carregando o arquivo com as geometrias usadas no experimento
scene placas = pov("tubo_quadrado.pov");

// Definindo o dominio fora (outside) dos solidos
domain DOMINIO = domain(placas, outside(<0,0,1>));

// Definindo solido com refinamento inferior ao da malha
mesh obstacle = surface(DOMINIO, MALHA);

// Funcao de penalizacao para nao realizar calculo dentro do solido
function P = 1E3 * (one(<0,0,1>)) + 1;

// Calculo da area
double area = int[MALHA](1);

function uin = 0.001;
function mu = 0.1;
femfunction u0(MALHA) = 0;
femfunction v0(MALHA) = 0;
femfunction w0(MALHA) = 0;
femfunction p0(MALHA_REFINADA) = 0;

double dt = 0.01;
double pp;
double qq;

double i = 1;

function u = 0;
function v = 0;
function w = 0;

do{
cout << "========== Tubo, Navier-Stokes passo " << i << "\n";

// Calculo da componente x
solve(u) in MALHA
memory(matrix=none),
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 MALHA xmin;
u = 0 on MALHA ymin;
u = 0 on MALHA ymax;
u = 0 on MALHA zmin;
u = 0 on MALHA zmax;
};

// Calculo da componente y
solve(v) in MALHA
memory(matrix=none),
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 MALHA xmin;
v = 0 on MALHA ymin;
v = 0 on MALHA ymax;
v = 0 on MALHA zmin;
v = 0 on MALHA zmax;
};

// Calculo da componente z
solve(w) in MALHA
memory(matrix=none),
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 MALHA xmin;
w = 0 on MALHA ymin;
w = 0 on MALHA ymax;
w = 0 on MALHA zmin;
w = 0 on MALHA zmax;
};

qq = int[MALHA](dx(u) + dy(v) + dz(w)) / area;
solve(q) in MALHA_REFINADA
krylov(type=cg,precond=ichol)
{
pde(q)
- div(dt * grad(q)) = (dx(u) + dy(v) + dz(w) - qq);

q = 0 on MALHA_REFINADA xmax;
};

p0 = p0 - q;
pp = int[MALHA](p0) / area;
p0 = p0 - pp;
u0 = u + (dt * dx(q));
v0 = v + (dt * dy(q));
w0 = w + (dt * dz(q));

i = i + 1;

} while(i <= 100);

// transforma o objeto em uma figura tetraédrica
mesh T = tetrahedrize(DOMINIO, MALHA_REFINADA);

// salva a malha
//save(vtk, "ns_placas",[u,v,w],MALHA_REFINADA);

//save(raw,"veloctuboquadrado00.vec", [u, v], MALHA_REFINADA);
//save(raw,"veloctuboquadrado003d.vec", [u, v, w], MALHA_REFINADA);
//save(medit,"tubo_quadrado_vel",T);
//save(medit,"tubo_quadrado_vel",[u,v,w],T);
//save(medit,"tubo_quadrado_pressao",T);
//save(medit,"tubo_quadrado_pressao",p0,T);


.pov
union{
box {<0.0,0.0,0.0> , <1.0,1.0,0.2>}
box {<0.0,0.0,0.2> , <1.0,0.2,0.8>}
box {<0.0,0.8,0.2> , <1.0,1.0,0.8>}
box {<0.0,0.0,0.8> , <1.0,1.0,1.0>}
pigment{ color rgb <0,0,1>}
}

_________________________________________________________________
Mande torpedos SMS do seu messenger para o celular dos seus amigos
http://mobile.msn.com/
Stephane Del Pino
2007-07-04 06:40:29 UTC
Permalink
Hello.
Post by Augusto Almeida
I need run a enviroment (a quadratic tube), with a viscous fluid in there,
the aim is steady state and export the vector components (discretization
vx vy vz
I am not sure I understand your question. If you want to save the velocity
field writting first all vx values then vy values and then vz values, use the
following command:

save(raw,"velocity.dat",[vx,vy,xz],mesh);

I hope that helps,

Best regards,
Stéphane.

Loading...