Questions about wall bouncing

[Darcher007]

I found that when running a model with a parallel phase interface, the liquid phase immediately below the wall would bounce up

Image at time step 0

And when looking at the velocity map, I noticed abnormal y-direction velocities at the upper and lower walls (I set the upper and lower boundaries to be walls, and the left and right cycle boundaries)

my XML
<Geometry nx="300" ny="150">
		<MRT>
			<Box />
		</MRT>
		<LIQUID name="blobb">
			 ...
		</LIQUID>
		<Wall mask="ALL" name="border">
			<Box ny="1" />
			<Box dy="-1" />
		</Wall> 
</Geometry>

According to the lattice Boltzmann bounce boundary and code, there should not be a velocity in the y-direction, and this velocity looks as if the upper and lower boundaries have periodicity. If my guess is true, can you tell me how to lift this periodicity?

[TravisMitchell]

Hi @Darcher007 - are you able to provide the specific model you are using and full xml file?

[llaniewski]

@Darcher007 "the upper and lower boundaries have periodicity", yes, in TCLB the domain is always periodic in all directions. But this should not be a problem if you added the upper and lower walls with

<Wall mask="ALL" name="border">
	<Box ny="1" />
	<Box dy="-1" />
</Wall> 

Nevertheless it would be useful to see what case (xml) are you using.

[Darcher007]

My XML file is similar to the following, It can be run with the "d2q9_ShanChen" model. I use TCLB-6.6.1

<?xml version="1.0"?>
<CLBConfig output="output/" permissive="true">
	<Geometry nx="128" ny="128">
		<MRT>
			<Box/>
		</MRT>
                <None name="blobb">
			<Box dx="0" nx="128" dy="0" ny="28"/>
		</None>
		<Wall mask="ALL" name="border">
			<Box ny="1"/>
			<Box dy="-1"/>
		</Wall>
	</Geometry>
	<Model>
		<Param name="Density" value="0.056"/>
		<Param name="Density" value="2.659" zone="blobb"/>
		<Param name="G_ff" value="-6.0"/>
		<Param name="VelocityX" value="0.1" zone="blobb" />
		<Param name="viscosity" value="0.166"/>
		<Params GravitationX="1e-5"/>
	</Model>
	<Solve Iterations="100" output="output/">
		<VTK Iterations="1"/>
	</Solve>
	<Solve Iterations="1000" output="output/">
		<VTK Iterations="10"/>
	</Solve>
	<Failcheck Iterations="1000"/>
	<Solve Iterations="10000" output="output/">
		<VTK Iterations="500"/>
	</Solve>
</CLBConfig>

The specific case is two-phase flow, with the liquid phase flowing on the lower wall . The model I am using is a multiple relaxation model. The specific problem is that when I write the program in C++ I am also applying a bounce boundary on the lower wall, at which point the liquid phase does not bounce off the wall; however, when I use TCLB's wall bounce boundary, the liquid phase bounces off the wall boundary.

I found that even when I used "mask=all", there seems to be a periodicity between the upper and lower wall surfaces, which is reflected in the appearance of the gas phase at the lower wall (see my last statement for the exact image)

Do you think the cause of the gas phase at the lower wall is periodicity or something else?

[Darcher007]

Also if this is periodicity causing a gas phase to appear underneath, why would adding a wall boundary still make the upper and lower boundaries periodic

[TravisMitchell]

The code assumes periodicity, but the wall condition will over-write the values (i.e. densities) streaming in from the other side of the domain.

I think the issue here is that you want to prescribe density for the walls, and the upper and lower walls will have different densities - something like the below should stop the bouncing up behaviour I believe. Your previous set up, when you masked all removed the additionals flag name "blobb" so the initial bottom wall density was set to 0.056.

`

<Geometry nx="128" ny="128">
	<MRT>
		<Box/>
	</MRT>
    <None name="blobb">
		<Box dx="0" nx="128" dy="0" ny="28"/>
	</None>
	<Wall mask="ALL" name="upper">
		<Box dy="-1"/>
	</Wall>
	<Wall mask="ALL" name="lower">
		<Box ny="1"/>
	</Wall>
</Geometry>
<Model>
	<Param name="Density" value="0.056"/>
	<Param name="Density" value="0.056" zone="upper"/>
	<Param name="Density" value="2.659" zone="blobb"/>
	<Param name="Density" value="2.659" zone="lower"/>
	<Param name="G_ff" value="-6.0"/>
	<Param name="VelocityX" value="0.1" zone="upper" />
	<Param name="viscosity" value="0.166"/>
	<Params GravitationX="1e-5"/>
</Model>

`

[Darcher007]

@TravisMitchell Thanks for the reminder that the liquid phase is now no longer bouncing, but the periodicity at the left and right boundaries still affects the fluid flow. I wrote a kind of inlet boundary and outlet boundary, and the fluid works fine when adding them individually to the XML, but when adding them together to the XML it reports the error "nan"

<?xml version="1.0"?>
<CLBConfig output="output/" permissive="true">
	<Geometry nx="300" ny="200">
		<MRT>
			<Box/>
		</MRT>
                <Win>
		        <Box dx="10" nx="1" />
		</Win>
		<Eout >
			<Box dx="-10" nx="1"/>
		</Eout> 
                <None name="blobb">
			<Box dx="0" nx="300" dy="0" ny="40"/>
		</None>
		<Wall mask="ALL" name="upper">
		         <Box dy="-1"/>
	         </Wall>
	         <Wall mask="ALL" name="lower">
		          <Box ny="1"/>
	         </Wall>
	</Geometry>
	<Model>
		<Param name="Density" value="0.056"/>
		<Param name="Density" value="0.056" zone="upper"/>
	        <Param name="Density" value="2.659" zone="blobb"/>
	        <Param name="Density" value="2.659" zone="lower"/>
		<Param name="G_ff" value="-6.0"/>
		<Param name="VelocityX" value="0.1" zone="blobb" />
		<Param name="viscosity" value="0.166"/>
		<Params GravitationX="1e-5"/>
	</Model>
	<Solve Iterations="100" output="output/">
		<VTK Iterations="1"/>
	</Solve>
	<Solve Iterations="1000" output="output/">
		<VTK Iterations="10"/>
	</Solve>
	<Failcheck Iterations="1000"/>
	<Solve Iterations="10000" output="output/">
		<VTK Iterations="500"/>
	</Solve>
</CLBConfig>

only Win（Non-equilibrium extrapolation boundaries）

only Eout（Neumann Boundary Condition）

both（about 200+ time step）

What's causing this to happen?

[TravisMitchell]

@Darcher007 this is hard for me to say, as I am unsure what you have implemented for Win and Eout. The main issue I can see with the xml file would be that you are offsetting the inlet and outlet away from the boundaries of the domain (dx=10 and dx=-10), is there a reason for this?

[Darcher007]

@TravisMitchell
Thanks for your previous reply, now I'm trying a new approach

Dynamics.R

AddDensity( name="f[0]", dx= 0, dy= 0, group="f")
AddDensity( name="f[1]", dx= 0, dy= 0, group="f")
AddDensity( name="f[2]", dx= 0, dy= 0, group="f")
AddDensity( name="f[3]", dx= 0, dy= 0, group="f")
AddDensity( name="f[4]", dx= 0, dy= 0, group="f")
AddDensity( name="f[5]", dx= 0, dy= 0, group="f")
AddDensity( name="f[6]", dx= 0, dy= 0, group="f")
AddDensity( name="f[7]", dx= 0, dy= 0, group="f")
AddDensity( name="f[8]", dx= 0, dy= 0, group="f")

AddField(name="q0_stream", stencil2d=1, group="streaming")
AddField(name="q1_stream", stencil2d=1, group="streaming")
AddField(name="q2_stream", stencil2d=1, group="streaming")
AddField(name="q3_stream", stencil2d=1, group="streaming")
AddField(name="q4_stream", stencil2d=1, group="streaming")
AddField(name="q5_stream", stencil2d=1, group="streaming")
AddField(name="q6_stream", stencil2d=1, group="streaming")
AddField(name="q7_stream", stencil2d=1, group="streaming")
AddField(name="q8_stream", stencil2d=1, group="streaming")

AddStage("BaseIteration", "Run"     ,  save=Fields$group %in% c("f""streaming", "neighbour_type_group","density") , load=DensityAll$group %in% c("f","neighbour_type_group","density"))
AddStage("PsiIteration" , "calcPsi" ,  save=Fields$name=="psi", load=DensityAll$group %in% c("f","density"))
AddStage("StreamIteration" , "Stream" ,save=Fields$group %in% c("f","density"), load=DensityAll$group %in% c("streaming","density")) #

AddAction("Init"     , c("BaseInit",      "PsiIteration","rhoIteration","uIteration"))
AddAction("Iteration", c("BaseIteration" ,"StreamIteration","PsiIteration","rhoIteration","uIteration"))

Dynamics.c

CudaDeviceFunction void Run() {
    rho = calcRho();
    if (NodeType & NODE_MRT) 
    {
        CollisionBGK();
    }
     neighbour_type = neighbour_type(0,0);
}
CudaDeviceFunction void Stream() {
 if (NodeType & NODE_MRT) 
    {
        Streaming();
    }
    
    switch (NodeType & NODE_BOUNDARY) {
    case NODE_Solid:
    case NODE_Wall:
        break;
    case NODE_Sboundary:
        Sboundary();
        break;
    case NODE_Nboundary:
        Nboundary();
        break;
    case NODE_Win:
        Win();
        break;  
    case NODE_Eout:
        Eout();
        break;    
    }
   

}
CudaDeviceFunction void Streaming() {
    f[0] = q0_stream(0,0);
    f[1] = q1_stream(-1,0);
    f[2] = q2_stream(0,-1);
    f[3] = q3_stream(1,0);
    f[4] = q4_stream(0,1);
    f[5] = q5_stream(-1,-1);
    f[6] = q6_stream(1,-1);
    f[7] = q7_stream(1,1);
    f[8] = q8_stream(-1,1);
}

It makes it easier for me to understand the code. What do you think of that?

[llaniewski]

Hi @Darcher007, when are you writing to the q?_stream fields?

If you want to copy fields from previous iteration (that what AddDensity(..., dx=0, dy=0, ...) does) you can do it just with f1=f1(0,0), without doubling the memory.

I don't exactly understand what the model is supposed to do. If you want to meet for a call, please contact me with an email.

[Darcher007]

CudaDeviceFunction void CollisionBGK() {
    f[?] = ....

    q0_stream = f[0];
    q1_stream = f[1];
    q2_stream = f[2];
    q3_stream = f[3];
    ....
}

The reason I created a new model like this is that I can't understand the model that comes with TCLB very well. In LBM theory, the initialization is done first, and then there is a three-step cycle of collision, migration, and boundary conditions, My model is written according to this theory, with q?_stream being set for migration

The most difficult thing to understand in TCLB is the data reading. In the manual, there is a paragraph

5 model infoemation/Dynamic.R

f_100 == f_100(-1,0,0) # the value is read from the left (`x=-1` direction) <==> it is streamed to the right (`x=1` direction).

It confuses me a lot. For example

d2q9_shanchen dynamic.c 
CudaDeviceFunction void CollisionBGK() {
    // Here we perform a single relaxation time collision operation.
    // We save memory here by using a single dummy variable
    real_t u_eq[2], f_temp[9];
  
    real_t d = getRho();

    vector_t u_eq_vect = getUeq(d);
    u_eq[0] = u_eq_vect.x; 
    u_eq[1] = u_eq_vect.y;

    f_temp[0] = f[0];
    f_temp[1] = f[1];
    f_temp[2] = f[2];
    f_temp[3] = f[3];
    f_temp[4] = f[4];
    f_temp[5] = f[5];
    f_temp[6] = f[6];
    f_temp[7] = f[7];
    f_temp[8] = f[8];
    SetEquilibrium(d, u_eq); //stores equilibrium distribution in f[0]-f[8]

    f[0] = f_temp[0] - omega*(f_temp[0]-f[0]);	
    f[1] = f_temp[1] - omega*(f_temp[1]-f[1]);
    f[2] = f_temp[2] - omega*(f_temp[2]-f[2]);
    f[3] = f_temp[3] - omega*(f_temp[3]-f[3]);	
    f[4] = f_temp[4] - omega*(f_temp[4]-f[4]);
    f[5] = f_temp[5] - omega*(f_temp[5]-f[5]);
    f[6] = f_temp[6] - omega*(f_temp[6]-f[6]);	
    f[7] = f_temp[7] - omega*(f_temp[7]-f[7]);
    f[8] = f_temp[8] - omega*(f_temp[8]-f[8]);
}

The f[1] in f_temp[1] = f[1]; here is f[1]==f[1](-1,0) or f[1]==f[1](0,0)
The same situation occurs at the setting of the boundary conditions

CudaDeviceFunction void EPressure()
{
        real_t ru, ux0;
	real_t rho = Density;
	ux0 = -1. + ( f[0] + f[2] + f[4] + 2.*(f[1] + f[5] + f[8]) ) / rho;
	ru = rho * ux0;

	f[3] = f[1] - (2./3.) * ru;
	f[7] = f[5] - (1./6.) * ru + (1./2.)*(f[2] - f[4]);
	f[6] = f[8] - (1./6.) * ru + (1./2.)*(f[4] - f[2]);
}

The f[1] in f[3] = f[1] - (2./3.) * ru; here is f[1]==f[1](-1,0) or f[1]==f[1](0,0)
Since the boundaries in my case are easily dispersed, and in this case I don't know if it's due to my wrong approach or the setup of TCLB, so I created the new model that came with my last reply

My e-mail address is [email protected] ,I hope you can remind me if I'm wrong about anything.