Hi,
Dear all,
I am simulating the flow and heat transfer in a cavity containing an obstacle(for simplifying),
| |
| |
| |
|_____________________| LX/4
| obstacle|
|_____________________|
The code that I use, buonce back conditions are applyed just in the hydrodynamic field, but not applyed in the temperature field,
For that, I add the BB conditions in the function of collision in the T field,
Unfortunately, the streamlines that I find are not good, especialy near the obstacle and in it,
You find below This function, would you mind helping me to correct this instructions, please?
[code=“cpp”]
//##############################################################################
//
// Copyright ©, 2005, Michael Sukop and Danny Thorne
//
// collide.c
//
//##############################################################################
//
// void collide( lattice_ptr lattice)
//##############################################################################
//void collide( lattice_ptr lattice)
//{
//double *feq;
//double *f;
// double *ftemp;
//#if ZHANG_AND_CHEN_ENERGY_TRANSPORT
// double force;
//#endif / ZHANG_AND_CHEN_ENERGY_TRANSPORT */
void collide( lattice_ptr lattice)
{
double *feq;
double *f;
double *ftemp;
#if ZHANG_AND_CHEN_ENERGY_TRANSPORT
double force;
#endif / ZHANG_AND_CHEN_ENERGY_TRANSPORT */
double omega;
int bc_type;
int n, a;
int subs;
double fl;
int i, j;
int ip, jp,
in, jn;
int ni = lattice->param.LX,
nj = lattice->param.LY;
for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++)
{
for( n=0; nNumNodes; n++)
{
feq = lattice->pdf[subs][n].feq;
f = lattice->pdf[subs][n].f;
ftemp = lattice->pdf[subs][n].ftemp;
bc_type = lattice->bc[subs][n].bc_type;
#if INAMURO_SIGMA_COMPONENT
if( ( get_bc_sigma_walls(lattice) && subs==1)
|| !( bc_type & BC_SOLID_NODE))
#else /* !(INAMURO_SIGMA_COMPONENT) /
if( !( bc_type & BC_SOLID_NODE))
#endif / (INAMURO_SIGMA_COMPONENT) */
{
// C O L L I D E
//getchar();
// f = ftemp - (1/tau[subs])( ftemp - feq)
for( a=0; a<=8; a++)
{
#if 1
f[a] = ftemp[a] - ( ( ftemp[a] / lattice->param.tau[subs] )
- ( feq[a] / lattice->param.tau[subs] ) );
#else
f[a] = ftemp[a] - ( ( ftemp[a] )
- ( feq[a] ) ) / lattice->param.tau[subs];
#endif
} /* for( a=0; a<=8; a++) */
#if PUKE_NEGATIVE_DENSITIES
for( a=0; a<=8; a++)
{
if( f < 0.)
{
printf("\n");
printf(
"collide() – Node %d (%d,%d), subs %d, "
“has negative density %20.17f "
“in direction %d "
“at timestep %d. Exiting!\n”,
n, n%lattice->param.LX,
n/lattice->param.LX,
subs,
f[a], a,
lattice->time );
printf(”\n”);
exit(1);
}
} / for( a=0; a<=8; a++) /
#endif / PUKE_NEGATIVE_DENSITIES */
} /* if( !( bc_type & BC_SOLID_NODE)) */
else // bc_type & BC_SOLID_NODE
{
f[1] = ftemp[1];
f[2] = ftemp[2];
f[3] = ftemp[3];
f[4] = ftemp[4];
f[5] = ftemp[5];
f[6] = ftemp[6];
f[7] = ftemp[7];
f[8] = ftemp[8];
// B O U N C E B A C K
if( lattice->param.bc_slip_north
&& n >= lattice->NumNodes - lattice->param.LX)
{
// Slip condition on north boundary.
f[1] = ftemp[1];
f[2] = ftemp[4];
f[3] = ftemp[3];
f[4] = ftemp[2];
f[5] = ftemp[8];
f[6] = ftemp[7];
f[7] = ftemp[6];
f[8] = ftemp[5];
} /* if( lattice->param.bc_slip_north && ... ) */
else
{
if( subs==0) // Hydrodynamic field
{
// Usual non-slip bounce-back condition.
f[1] = ftemp[3];
f[2] = ftemp[4];
f[3] = ftemp[1];
f[4] = ftemp[2];
f[5] = ftemp[7];
f[6] = ftemp[8];
f[7] = ftemp[5];
f[8] = ftemp[6];
} /* if( subs==0) */
#if NUM_FLUID_COMPONENTS==2
else // subs==1 // T field
{
#if INAMURO_SIGMA_COMPONENT
if( lattice->param.bc_sigma_slip)
{
//
// Slip BC for solute on side walls.
// Will this make a difference on Taylor dispersion?
//
if( lattice->FlowDir == /Vertical/2)
{
if( /*west*/(n )%lattice->param.LX == 0
|| /*east*/(n+1)%lattice->param.LX == 0)
{
// Slip condition on east/west boundary
f[1] = ftemp[3];
f[2] = ftemp[2];
f[3] = ftemp[1];
f[4] = ftemp[4];
f[5] = ftemp[6];
f[6] = ftemp[5];
f[7] = ftemp[8];
f[8] = ftemp[7];
}
}
else if( lattice->FlowDir == /*Horizontal*/1)
{
if( /*north*/ n >= lattice->NumNodes - lattice->param.LX
|| /*south*/ n < lattice->param.LX )
{
// Slip condition on north/south boundary.
f[1] = ftemp[1];
f[2] = ftemp[4];
f[3] = ftemp[3];
f[4] = ftemp[2];
f[5] = ftemp[8];
f[6] = ftemp[7];
f[7] = ftemp[6];
f[8] = ftemp[5];
}
else
{
// ERROR: Solid exists somewhere other than as side walls.
printf("%s (%d) >> "
"ERROR: "
"bc_sigma_slip is on. "
"FlowDir is determined to be horizontal. "
"Encountered solid node somewhere other than side walls. "
"That situation is not supported. "
"Exiting!", __FILE__, __LINE__);
exit(1);
}
}
else
{
printf("%s (%d) >> "
"FlowDir is indeterminate. "
"Cannot apply slip BC (bc_sigma_slip). "
"Exiting!", __FILE__, __LINE__);
exit(1);
}
} /* if( lattice->param.bc_sigma_slip) */
else if((n > lattice->param.LX *(lattice->param.LX/4 -1)+1)&& (n < lattice->param.LX *lattice->param.LX/4))
{ // The Buonce back instructions that I add for T field (subs==1)
#endif /* INAMURO_SIGMA_COMPONENT /
// Usual non-slip bounce-back conditions
f[4] = ftemp[2];
f[7] = ftemp[5];
f[8] = ftemp[6];
}
else {
f[1] = ftemp[3];
f[2] = ftemp[4];
f[3] = ftemp[1];
f[4] = ftemp[2];
f[5] = ftemp[7];
f[6] = ftemp[8];
f[7] = ftemp[5];
f[8] = ftemp[6];
#if INAMURO_SIGMA_COMPONENT
} / if( lattice->param.bc_sigma_slip) else /
#endif / INAMURO_SIGMA_COMPONENT */
} /* if( subs==0) else*/
#endif /* NUM_FLUID_COMPONENTS==2 */
} /* if( lattice->param.bc_slip_north && ... ) else */
} /* if( !( bc_type & BC_SOLID_NODE)) else */
} /* for( n=0; n<lattice_NumNodes; n++) */
if( subs==0)
{
//printf(“test_inamuro_collllid-PC_effetporus_CL %d\n”,subs);
// Compute the solid density term for fluid component.
for( n=0; nNumNodes; n++)
{
//printf(“debut boucle %f\n”);
//if(lattice->fs1[n].fs>0.0)
//{
ftemp = lattice->pdf[subs][n].ftemp;
bc_type = lattice->bc [subs][n].bc_type;
i = n%ni;
j = n/ni;
jp = ( j<nj-1)?( j+1):( 0 );
jn = ( j>0 )?( j-1):( nj-1);
ip = ( i<ni-1)?( i+1):( 0 );
in = ( i>0 )?( i-1):( ni-1);
//printf("n= %d,fl= %f\n"n,lattice->fs1[n].fs);
//getchar();
if( !( bc_type & BC_SOLID_NODE))
{
//printf("test_inamuro_collllid-PC_effetporus_CL2 %d %d\n",n,subs);
// Variable solid density.
fl =(1-(lattice->fs1[n].fs));
//printf("n= %d,fl= %f\n",n,fl);
if((fl<1)&(fl>0))
{
fl=fl*fl*fl/((1.0-fl)*(1.0-fl)+1.E-2)/1.E+12; //-2
}
//printf("fl2=%d %f %f%\n",n,fl,lattice->fs1[n].fs);
/* 0 */ ftemp++;
/* 1 */ ftemp++ = fl( lattice->pdf[subs][ j *ni + ip].f[3]
- lattice->pdf[subs][ j ni + i ].f[1]);
/ 2 / ftemp++ = fl( lattice->pdf[subs][ jpni + i ].f[4]
- lattice->pdf[subs][ j ni + i ].f[2]);
/ 3 */ ftemp++ = fl( lattice->pdf[subs][ j *ni + in].f[1]
- lattice->pdf[subs][ j ni + i ].f[3]);
/ 4 / ftemp++ = fl( lattice->pdf[subs][ jnni + i ].f[2]
- lattice->pdf[subs][ j ni + i ].f[4]);
/ 5 / ftemp++ = fl( lattice->pdf[subs][ jpni + ip].f[7]
- lattice->pdf[subs][ j ni + i ].f[5]);
/ 6 / ftemp++ = fl( lattice->pdf[subs][ jpni + in].f[8]
- lattice->pdf[subs][ j ni + i ].f[6]);
/ 7 / ftemp++ = fl( lattice->pdf[subs][ jnni + in].f[5]
- lattice->pdf[subs][ j ni + i ].f[7]);
/ 8 / ftemp++ = fl( lattice->pdf[subs][ jnni + ip].f[6]
- lattice->pdf[subs][ j ni + i ].f[8]);
} / if( !( *bc_type++ & BC_SOLID_NODE)) /
} / for( n=0; n<lattice_NumNodes; n++) */
for( n=0; n<lattice->NumNodes; n++)
{
f = lattice->pdf[subs][n].f;
bc_type = lattice->bc [subs][n].bc_type;
if( !( bc_type & BC_SOLID_NODE))
{
f++;
for( a=1; a<9; a++, f++)
{
*f += *(f+9);
} /* for( a=1; a<9; a++) */
} /* if( !( bc_type & BC_SOLID_NODE)) */
// } // if (lattice->fs1[n].fs>0)
} /* for( n=0; n<lattice->NumNodes; n++, f+=18) */
} /* if( subs==0) */
} /* for( subs=0; subs<NUM_FLUID_COMPONENTS; subs++) */
#if SAY_HI
printf(“collide() – Bye!\n”);
#endif /* SAY_HI */
} /* void collide( lattice_ptr lattice) */
Best regards,
Bolt
[edit yann@flowkit] Pleas try to use the "formatted code" tag when posting cod in this forum.