My case of dam break flow-how to use movable boundaries

[syy-mudong]

Hi, First of all, thank you very much for allowing me to use this great program for free. As someone who has just come into contact with the LBM field. I made some modifications to the dam break flow and wanted to obtain a simulation result within 5 seconds, but I found some problems after outputting it. I guess there may be an issue with my LBM unit settings.

May I ask for my set_ What are the issues with up? I will be extremely grateful. And thank you for your work! Shen

setup.cpp

void main_setup() { // dam break
    // ######################################################### define simulation box size, viscosity and volume force ############################################################################

    float const si_nu = 1E-6f; // kinematic shear viscosity [m^2/s] 
    const float si_rho = 997.f; // fluid density [kg/m^3] 
    const float si_sigma = 73.81E-3f; // fluid surface tension [kg/s^2] 
    const float si_u = 100.; //velocity [m/s]
    const float si_x = 20;
    const float si_g = 9.8f;

    const uint size = 100;
    const int box_size = size * 2u;
    const float u = 0.1; //velocity
    const float rho = 1.0f; // density
    units.set_m_kg_s((float)box_size, u, rho, si_x, si_u, si_rho);

    const float nu = units.nu(si_nu); // calculate values for remaining parameters in simulation units
    const float sigma = units.sigma(si_sigma);
    const float si_f = units.f_from_g(si_g, si_rho);
    const float f = units.f(si_f);

    LBM lbm(size , size * 2u, size * 2u, nu, 0.0f, 0.0f, -f, sigma);
    // #############################################################################################################################################################################################
    const ulong N=lbm.get_N(); const uint Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); for(ulong n=0ull; n<N; n++) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
        // ########################################################################### define geometry #############################################################################################
        if(z<Nz*6u/8u && y<Ny/8u) lbm.flags[n] = TYPE_F;
        if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==0u||z==Nz-1u) lbm.flags[n] = TYPE_S; // all non periodic
    }   // #########################################################################################################################################################################################
    lbm.run(0u);
    while (lbm.get_t() < units.t(5.0f))
    {
        lbm.phi.write_device_to_vtk();
        lbm.run(20u);
    }
    print_info("m: " + to_string(units.si_x(1u)) + " kg: " + to_string(units.si_M(1.0f)) + " s: " + to_string(units.si_t(1.0)) + " Force: " + to_string(units.si_F(1.f)));
} /**/

The results were perfect before 709

Afterwards, the problem arose

[syy-mudong]

I probably found my problem. In my case, the Reynolds number is very large, but I did not uncomment SUBGRID. When I uncomment this item, it can run normally.

[syy-mudong]

Now I have a new problem. When I try to add a rising gate to simulate a dam break flow with movable boundaries, my boundaries have not changed. What are some suggestions for using movable boundaries.

void main_setup() { // dam break
    // ######################################################### define simulation box size, viscosity and volume force ############################################################################

    float const si_nu = 1E-6f; // kinematic shear viscosity [m^2/s] 
    const float si_rho = 997.f; // fluid density [kg/m^3] 
    const float si_sigma = 73.81E-3f; // fluid surface tension [kg/s^2] 
    const float si_u = 6.; //velocity [m/s]
    const float si_x = 0.2;
    const float si_g = 9.8f;

    const uint size = 40;
    const float u = 0.1; //velocity
    const float rho = 1.0f; // density
    units.set_m_kg_s((float)size, u, rho, si_x, si_u, si_rho);

    const float nu = units.nu(si_nu); // calculate values for remaining parameters in simulation units
    const float sigma = units.sigma(si_sigma);
    const float si_f = units.si_f_from_si_g(si_g, si_rho);
    const float f = units.f(si_f);

    LBM lbm(size  , size * 4u, size * 3u, nu, 0.0f, 0.0f, -f, sigma);
    // #############################################################################################################################################################################################
    const ulong N=lbm.get_N(); const uint Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); for(ulong n=0ull; n<N; n++) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
        // ########################################################################### define geometry #############################################################################################
        if(z<Nz * 3u/4u && y<Ny/4u) lbm.flags[n] = TYPE_F;
        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 (y == Ny / 4u ) {
            lbm.flags[n] = TYPE_S;
            lbm.u.z[n] = units.u(si_u);
        }
    }   // #########################################################################################################################################################################################
    lbm.run();
    /*while (lbm.get_t() < units.t(10.0f))
    {
        lbm.phi.write_device_to_vtk();
        lbm.run(units.t(0.02f));
    }*/
    print_info( "m: " + to_string(units.si_x(1u)) + " kg: " + to_string(units.si_M(1.0f)) + " s: " + to_string(units.si_t(1.0)) + " rho: " + to_string(units.si_rho(1.f)));
} /**/

[ProjectPhysX]

Hi @syy-mudong,

yes at too large Reynolds numbers it will become unstable, and then you need to enable the SUBGRID model extension.

The "moving boundaries" refer to momentarily stationary boundaries with an instantaneous velocity. For example a rotating wheel has a tangential velocity at the boundaries, but the boundaries themselves remain in place. When the velocity is not tangential, like for a rising gate, you have to re-voxelize, or change the gate position in the LBM flags every few time steps:

1. Clear the old gate flags.
2. Set the flags at the new, slightly moved position. The movement in every step should not be more than 1 cell.
3. Update the modified flag lattice in VRAM with lbm.flags.write_to_device().
4. Initialize the new moving boundaries with lbm.update_moving_boundaries();.
5. Run 1-10 time steps with lbm.run(10u).

Regards, Moritz

[syy-mudong]

Thank you for your suggestion. I will try to create a new set_ Up to achieve my goal. Finally, this is really a great program,.As a student, it deepened my understanding of LBM in the process of learning it.

Thank you for your work! Shen

[syy-mudong]

I'm sorry to bother you @ProjectPhysX . Because the powerful computing power of this program fascinates me, I want to use it to help me complete my work. I would like to ask if it is possible to simulate the situation where a small ball collides and enters the water surface. Previously, opening the gate required manual updating of the voxel, but now I want to simulate fluid structure coupling. Can I do this? In addition, I want to output the pressure exerted by the dam break flow on the boundary, but through F.write_ device_ to_ vtk . The output is not pressure(Pa), may I ask if I want to output pressure? How to make my set_up?

Shen

[syy-mudong]

Hi @ProjectPhysX ,

Through reading the code during this period, I have roughly understood that F represents the volumetric force at the LBM point, so it cannot to get the pressure (Pa) by F.write_ device_ to_ vtk. Because of my poor coding ability, there may be a misunderstanding in my understanding of the code.If you could answer my question in your free time, I would greatly appreciate it.

1---

void LBM::calculate_force_on_boundaries() { // calculate forces from fluid on TYPE_S nodes
    for(uint d=0u; d<get_D(); d++) lbm[d]->enqueue_calculate_force_on_boundaries();
    for(uint d=0u; d<get_D(); d++) lbm[d]->finish_queue();
}
float3 LBM::calculate_force_on_object(const uchar flag_marker) { // add up force for all nodes flagged with flag_marker
    double3 force(0.0, 0.0, 0.0);
    for(ulong n=0ull; n<get_N(); n++) {
        if(flags[n]==flag_marker) {
            force.x += (double)F.x[n];
            force.y += (double)F.y[n];
            force.z += (double)F.z[n];
        }
    }
    return float3(force.x, force.y, force.z);
}

If I want to calculate the pressure of a small region I am interested in, I need to call these two functions to first calculate the force of the region I am interested in, and then I use float3/dx^2 to obtain the pressure (Pa).

2--- Secondly, in my understanding, nu is related to the relaxation time. Perhaps due to my carelessness, I did not see how nu works in the code. Now I have an idea if it is possible to change nu after running for a period of time to achieve the purpose of simulating non-Newtonian fluid.

Best regards, Shen