Calculating a jet striking a plate

[PetterTh]

Hi

Thank you for a very impressive tool!

I'm trying to setup a calculation for a water jet striking a plate like shown below. However, I'm not able to get the inlet/outlet working. How do I setup the solver so that there is a circular inlet boundary on one side where the water-jet can enter and strike the plate on the other side before the water exits on the side of the domain?

Best Regards Petter

[PetterTh]

I was able to make it work like this: `

void main_setup() { // Jet hitting wall

// ######################################################### define simulation box size, viscosity and volume force ############################################################################
const uint size = 96;
uint r=0;
uint ncells = 0;

const int box_diameter = size*3u; 
float jet_diameter = box_diameter / 4;

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_sos = 1400.f; // Speed of sound in water [m/s]

const float si_D =0.05; // jet diameter [m] 
const float si_u = 100.; // jet velocity [m/s] 

const float u = 1/sqrt(3)/si_sos*si_u; // jet velocity
const float rho = 1.0f; // density
units.set_m_kg_s((float)jet_diameter, u, rho, si_D, si_u, si_rho); // calculate 3 independent conversion factors (m, kg, s)
print_info("m" + to_string(units.get_m())+" kg "+ to_string(units.get_kg()) + " s " + to_string(units.get_s())+" Force "+to_string(units.si_F(1.f)));
float vel = u;

const float nu = units.nu(si_nu); // calculate values for remaining parameters in simulation units
const float sigma = units.sigma(si_sigma);

LBM lbm(size * 3u, size * 3u, size * 1u, nu, 0.0f, 0.0f, 0.0f, sigma);

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)) {
    const uint D = max(Nx, Nz);
     r = jet_diameter/2;

    if ((z <= 5u) && sq(x - Nx / 2) + sq(y - Nx / 2) < sq(r)) { // Initial waterjet
        lbm.flags[n] =  TYPE_F; 
        lbm.u.z[n] = vel;;
        lbm.u.x[n] = 0.0f;;
        lbm.u.y[n] = 0.0f;;
        lbm.rho[n] = 1.0f;;
        if (z == 1u) {
            ncells += 1;
        }
    }
    if (z == Nz - 1u) lbm.flags[n] = TYPE_S; // Walls on the top plane

    if (z == Nz - 1u && sq(x - Nx / 2) + sq(y - Nx / 2) < sq(D / 3)) lbm.flags[n] = TYPE_X+TYPE_S; // To Calculate Force away from the boundaryconditions

    if ((z == 3u)) { // Inlet condition

        if (sq(x - Nx / 2) + sq(y - Nx / 2) < sq(r)) {
            lbm.flags[n] = TYPE_E;
            lbm.u.z[n] = vel;;
            lbm.u.x[n] = 0.0f;;
            lbm.u.y[n] = 0.0f;;

        }
    }
    if (z <= 10u && sq(x - Nx / 2) + sq(y - Nx / 2) > sq(r)){ // Floor
            lbm.flags[n] = TYPE_S;
    }
    else if (x <= 1u ) { // Outlet
        lbm.flags[n] = TYPE_E;
        lbm.u.x[n] = -vel;
    }
    else if (x >= Nx-2u) { // Outlet
        lbm.flags[n] = TYPE_E;
        lbm.u.x[n] = vel;
    }
    else if (y <= 1u) { // Outlet
        lbm.flags[n] = TYPE_E;
        lbm.u.y[n] = -vel;
    }
    else if (y >= Nx - 2u) { // Outlet
        lbm.flags[n] = TYPE_E;
        lbm.u.y[n] = vel;
    }
}

lbm.run(0u);
print_info("Inlet Cells:" + to_string(ncells));
int c = 0;
float forceAve = 0;

int nrun = 50;
while (running) {

    lbm.run(nrun);

    c++;
    if (c > int(Nz/vel/nrun)) { // Don't collect data until the jet has hit the back plate
        lbm.calculate_force_on_boundaries();
        lbm.F.read_from_device();
        const float3 force = lbm.calculate_force_on_object(TYPE_X + TYPE_S);
        forceAve += force.z;

        if (c % 25 == 0) {
            lbm.u_write_device_to_vtk("u.vtk");
            lbm.rho_write_device_to_vtk("rho.vtk");
            print_info("AverageForce " + to_string(forceAve / (25.)) + " N"+ " Estimated Force "+to_string(ncells*sq(vel)*1.0));
            print_info("SI_Units" + to_string(units.si_F(forceAve) / (25.)));
            forceAve = 0.;
        }
    } 
}

 }`

However the calculated force acting on the plate is not calculated according to the analytical solution. I should be $F= \rho A v^2$, however the calculated force does not seem to scale with the velocity squared. However, if I change the function calculate_force_on_boundaries() from

    F[                 n] = 2.0f*fx*Fb; // 2 times because fi are reflected on solid boundary nodes (bounced-back)
    F[    def_N+(ulong)n] = 2.0f*fy*Fb;
    F[2ul * def_N + (ulong)n] = 2.0f * fz* Fb; 

to

    F[                 n] = 4.0f*sq(fx)*Fb; // 2 times because fi are reflected on solid boundary nodes (bounced-back)
    F[    def_N+(ulong)n] = 4.0f*sq(fy)*Fb;
    F[2ul * def_N + (ulong)n] = 4.0f * sq(fz)* Fb; 

I get a result which is quite close to the analytical result. Do you think this change makes sense? I've never used LBM before so I might be totally wrong, but should it not be a pressure component in the calculated force aswell?

[ProjectPhysX]

Hi @PetterTh! Sorry for the delayed response; I just tested your setup. The force implementation with the "2" is correct, changing it to "4" makes it wrong.

There is some potential issues:

• Make sure the verlocity in LBM units is <<0.2. This does not have to be bound to the SI speed of sound in any way, you can really choose it freely, and the viscosity and surface tension in LBM units get adjusted automatically to retain Reynolds and Capillary number.
• The subgrid model may make forces inaccurate here. Quite often there is only a very narrow corridor for suitable parametrization without subgrid model; I'm not sure if this could work here at all. Maybe test with lower jet velocity if you have corresponding experimental data.
• There is no wetting model implemented. Forces only work accurately when the solid boundary is entirely wetted. Shouldn't be too much of an issue here as the main part of the striking force is in the middle.