Rayleigh-Benard convection blows up

[rodionstepanov]

The simulation with (main_setup code below) seems to blows up between 77000 and 78000 steps. I have changed just the box size and boundary condition. In defines.hpp #define TEMPERATURE and #define SUBGRID are enable. By the way it is not clear where additional parameters for convection, like, thermal diffusivity set.

void main_setup() { // Rayleigh-Benard convection
// ######################################################### define simulation box size, viscosity and volume force ############################################################################
LBM lbm(256u, 128u, 128u, 0.02f, 0.0f, 0.0f, -0.001f, 0.0f, 1.0f, 1.0f);
// #############################################################################################################################################################################################
const uint N=lbm.get_N(), Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); for(uint n=0u, x=0u, y=0u, z=0u; n<N; n++, lbm.coordinates(n, x, y, z)) {
// ########################################################################### define geometry #############################################################################################
lbm.u.x[n] = random_symmetric(0.015f);
lbm.u.y[n] = random_symmetric(0.015f);
lbm.u.z[n] = random_symmetric(0.015f);
if(z==1u) {
lbm.T[n] = 1.75f;
lbm.flags[n] = TYPE_T;
} else if(z==Nz-2u) {
lbm.T[n] = 0.25f;
lbm.flags[n] = TYPE_T;
}
if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==0u||z==Nz-1u) lbm.flags[n] = TYPE_S; // all non periodic
//if(z==0u||z==Nz-1u) lbm.flags[n] = TYPE_S;
} // #########################################################################################################################################################################################
key_4 = true;
Clock clock;
lbm.run(0u);
while(lbm.get_t()<1000000u) {
lbm.graphics.set_camera_free(float3(1.0f*(float)Nx, -0.4f*(float)Ny, 2.0f*(float)Nz), -33.0f, 42.0f, 68.0f);
lbm.graphics.write_frame_png(get_exe_path()+"export/t/");
lbm.graphics.set_camera_free(float3(0.5f*(float)Nx, -0.35f*(float)Ny, -0.7f*(float)Nz), -33.0f, -40.0f, 100.0f);
lbm.graphics.write_frame_png(get_exe_path()+"export/b/");
lbm.graphics.set_camera_free(float3(0.0f*(float)Nx, 0.51f*(float)Ny, 0.75f*(float)Nz), 90.0f, 28.0f, 80.0f);
lbm.graphics.write_frame_png(get_exe_path()+"export/f/");
lbm.graphics.set_camera_free(float3(0.7f*(float)Nx, -0.15f*(float)Ny, 0.06f*(float)Nz), 0.0f, 0.0f, 100.0f);
lbm.graphics.write_frame_png(get_exe_path()+"export/s/");
lbm.run(1000u);
}
write_file(get_exe_path()+"time.txt", print_time(clock.stop()));
lbm.run();
} /**/

[MarcoAurelioFerrari]

The additional parameters are set on lbm(256u, 128u, 128u, 0.02f, 0.0f, 0.0f, -0.001f, 0.0f, 1.0f, 1.0f);
For your case, you set the domain as 256x128x128. With a viscosity of 0.02, body force term of -0.001 in the z-direction, alpha =1, and beta =1

It is important to notice that when you increased the domain in the z direction from 64 to 128, you effectively also increased the Rayleigh number by 8x. Since you didn't change any other property, the maximum fluid velocity will also increase, and if the velocity is above the compressibility limit, the simulation will diverge.

It could also be a boundary condition problem since you change it too. If you make the domain periodic again, does the simulation diverges?

[rodionstepanov]

@MarcoAurelioFerrari thanks, I see your point but LES model is supposed to stabilize simulation

[ProjectPhysX]

@rodionstepanov I can confirm that this blows up :)
Problem here is that the temperature gradient in LBM units is very large (1.75-0.25), volume force is large (-0.001f), and box height is large (128), so velocity will get larger than the LBM speed of sound (0.57) and then it blows up. Finding suitable parameters is an art in itself. Note that you don't need SUBGRID here.

Thermal diffusion and expansion coefficients are the alpha and beta in the LBM constructor. These are in LBM units, and you'll have to do the unit conversion between SI and LBM units manually.