How to compute drag and lift forces on voxelized objects?

[gustavo8942]

Hello, an amazing program you have developed, I appreciate that you published it for free. I'm simulating a high lift airfoil at high Reynolds, which was a nightmare with Ansys fluent. My goal is to obtain the drag and lift coefficients or at least the force at the airfoil, is it possible or implemented in the program? I exported the velocity and force vtk file, but the force vtk when I opened it at OpenFoam it showed nothing. Thanks :)

https://drive.google.com/file/d/1oRb9zjwCpU86CRqfl6DP671Ibxth5gM1/view?usp=sharing (a short clip of the simulation)

[ProjectPhysX]

Yes everything you need is implemented already. Quick guide through the setup:

void main_setup() {
    // make a list of all variables you have in SI units (m, kg, s)
    const float si_x = 2.0f; // 2m wingspan
    const float si_u = 3.0f; // 3m/s wind speed
    const float si_rho = 1.225f; // 1.225kg/m^3 air density
    const float si_nu = 1.48E-5f; // 1.48E-5f m^2/s kinematic shear viscosity of air

    // set velocity in LBM units, should be between 0.01-0.2
    const float lbm_u = 0.1f;
    const float lbm_rho = 1.0f; // density in LBM units always has to be 1

     // set unit conversion factors between SI units and LBM units
    units.set_m_kg_s(lbm_x, lbm_u, lbm_rho, si_x, si_u, si_rho);

    // compute kinematic shear viscosity ion LBM units
    const float lbm_nu = units.nu(si_nu);

    // set grid resolution based on lbm_x
    const uint Nx = to_uint(lbm_x);
    const uint Ny = to_uint(lbm_x);
    const uint Nz = to_uint(lbm_x);

    LBM lbm(Nx, Ny, Nz, lbm_nu ); // create LBM object

    // load geometry from stl file, mark all grid points that belong to your wing geometry with (TYPE_S|TYPE_X)
    lbm.voxelize_stl(get_exe_path()+"../stl/wing.stl", center, size, TYPE_S|TYPE_X);

    // set box boundary conditions
    for(uint n=0u, x=0u, y=0u, z=0u; n<lbm.get_N(); n++, lbm.coordinates(n, x, y, z)) {
        if(!(lbm.flags[n]&TYPE_S)) lbm.u.y[n] = lbm_u; // initial velocity
        if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==0u||z==Nz-1u) lbm.flags[n] = TYPE_E;
        if(z==0u) lbm.flags[n] = TYPE_S;
    }

    // run simulation
    lbm.run(units.t(0.4f)); // run for 0.4 seconds

    // calculate forces on boundaries on GPU, then copy force field to CPU memory
    lbm.calculate_force_on_boundaries();
    lbm.F.read_from_device();

    // optional: export force field as .vtk file
    lbm.F.write_host_to_vtk(); // F is already copied from GPU memory to host memory, no need for ..._write_device_to_vtk() here

    // sum force over all boundary nodes marked with (TYPE_S|TYPE_X)
    const float3 lbm_force = lbm.calculate_force_on_object(TYPE_S|TYPE_X);

    // transition back to SI units
    const float si_force_x = units.si_F(lbm_force.x); // force in Newton
    const float si_force_y = units.si_F(lbm_force.y);
    const float si_force_z = units.si_F(lbm_force.z);

    // print lift/drag force components
    print_info("lift force = "to_string(si_force_z)+" N");
    print_info("drag force = "to_string(si_force_y)+" N");
}

[randomwangran]

@ProjectPhysX

    // transition back to SI units
    const float si_force_x = units.si_F(lbm_force.x); // force in Newton
    const float si_force_y = units.si_F(lbm_force.x);
    const float si_force_z = units.si_F(lbm_force.x);

Is this a typo or something that I don't know?

    // transition back to SI units
    const float si_force_x = units.si_F(lbm_force.x); // force in Newton
    const float si_force_y = units.si_F(lbm_force.y);
    const float si_force_z = units.si_F(lbm_force.z);

[randomwangran]

constg float lbm_nu = units.nu(si_nu);

should be?

const float lbm_nu = units.nu(si_nu);

[ProjectPhysX]

@randomwangran yep was typos, thanks! Fixed it.

[randomwangran]

@ProjectPhysX

Also:

    print_info("lift force = "to_string(si_force_z)+" N");
    print_info("drag force = "to_string(si_force_y)+" N");

should be:

   print_info("lift force = " + std::to_string(si_force_z) + " N");
   print_info("drag force = " + std::to_string(si_force_y) + " N");

[ProjectPhysX]

@randomwangran thanks, the to_string is correct, I'm using my own conversion functions.