Propeller rotation axis off-axis during voxelization

[wjsjtu123]

When performing a propeller rotational motion, I define a rotational motion around the x-axis with a rotational speed of 4145r/min, but in the voxelized motion, instead of following a fixed x-axis, the propeller will rotate periodically off-center, what is the cause of this, and how can I fix it? My code is as follows

void main_setup() { // NACA X-57 lift propeller; required extensions in defines.hpp: FP16C, EQUILIBRIUM_BOUNDARIES, MOVING_BOUNDARIES, SUBGRID, INTERACTIVE_GRAPHICS or GRAPHICS
// ################################################################## define simulation box size, viscosity and volume force ###################################################################
    const uint3 lbm_N = resolution(float3(2.0f, 1.0f, 1.0f), 4000u); // input: simulation box aspect ratio and VRAM occupation in MB, output: grid resolution
    const float lbm_u = 0.11f;
    const float lbm_length = 0.5f*(float)lbm_N.x;
    const float si_n =4145.0f; //RPM of propeller
    const float si_n_t = 20.0f;
    const float si_T = (60.0f/si_n)si_n_t; // 20 revolutions of the propeller
    const uint lbm_dt = 30u; // revoxelize rotor every dt time steps
    const float si_length=0.984f;
    const float si_J = 0.6f; // Advance ratio
    const float si_D = 0.57592f; //propeller radius
    const float si_u = si_Jsi_D*(si_n/60.0f); //airspeed
    const float si_w_u = (si_npif/30.0f)(si_D0.5f); //propeller linear velocity
    const float si_nu=1.48E-5f, si_rho=1.225f;
    units.set_m_kg_s(lbm_length, lbm_u, 1.0f, si_length, si_u, si_rho);
    print_info("lbm_u="+to_string(units.u(si_u)));
    LBM lbm(lbm_N, 1u, 1u, 1u, units.nu(si_nu));
    // ###################################################################################### define geometry ######################################################################################
    Mesh fuselage = read_stl(get_exe_path()+"../stl/propeller_fuselage_semi.stl"); //
    Mesh* propeller = read_stl(get_exe_path()+"../stl/propeller_semi.stl"); //
    const float scale = lbm_length/fuselage->get_bounding_box_size().x; // scale body and rotors to simulation box size
    fuselage->scale(scale);
    propeller->scale(scale);
    const float3 offset = lbm.center()-fuselage->get_bounding_box_center(); // move body and rotors to simulation box center
    fuselage->translate(offset);
    propeller->translate(offset);
    fuselage->set_center(fuselage->get_bounding_box_center()); // set center of meshes to their bounding box center
    propeller->set_center(propeller->get_bounding_box_center());
    print_info("diameter="+to_string(units.x(si_D)));
    print_info("propeller_max_size="+to_string(propeller->get_max_size()));
    const float propeller_radius=0.5funits.x(si_D), propeller_omega=(units.si_u(si_w_u))/(propeller_radius), propeller_domega=propeller_omega(float)lbm_dt;
    lbm.voxelize_mesh_on_device(fuselage);
    //lbm.voxelize_mesh_on_device(propeller, TYPE_S, propeller->get_center(), float3(0.0f), float3(propeller_omega, 0.0f, 0.0f)); // revoxelize mesh on GPU
    const uint Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); parallel_for(lbm.get_N(), [&](ulong n) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
    if(lbm.flags[n]!=TYPE_S) lbm.u.x[n] = lbm_u;
    if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==0u||z==Nz-1u) lbm.flags[n] = TYPE_E; // all non periodic
    }); // ####################################################################### run simulation, export images and data ##########################################################################
    lbm.graphics.visualization_modes = VIS_FLAG_SURFACE|VIS_Q_CRITERION;
    lbm.run(0u); // initialize simulation
    while(lbm.get_t()<=units.t(si_T)) { // main simulation loop
        lbm.voxelize_mesh_on_device(propeller, TYPE_S, propeller->get_center(), float3(0.0f), float3(propeller_omega, 0.0f, 0.0f)); // revoxelize mesh on GPU
        lbm.run(lbm_dt); // run dt time steps
        propeller->rotate(float3x3(float3(1.0f, 0.0f, 0.0f), propeller_domega)); // rotate mesh
#if defined(GRAPHICS) && !defined(INTERACTIVE_GRAPHICS)
        if(lbm.graphics.next_frame(units.t(si_T), 10.0f)) {
            lbm.graphics.set_camera_free(float3(0.528513f*(float)Nx, 0.102095f*(float)Ny, 1.302283f*(float)Nz), 16.0f, 47.0f, 96.0f);
            lbm.graphics.write_frame(get_exe_path()+"export/a/");
            lbm.graphics.set_camera_free(float3(0.0f*(float)Nx, -0.114244f*(float)Ny, 0.543265f*(float)Nz), 90.0f+degrees((float)lbm.get_t()/(float)lbm_dtmain_domega), 36.0f, 120.0f);
            lbm.graphics.write_frame(get_exe_path()+"export/b/");
            lbm.graphics.set_camera_free(float3(0.557719f(float)Nx, -0.503388f*(float)Ny, -0.591976f*(float)Nz), -43.0f, -21.0f, 75.0f);
            lbm.graphics.write_frame(get_exe_path()+"export/c/");
            lbm.graphics.set_camera_centered(58.0f, 9.0f, 88.0f, 1.648722f);
            lbm.graphics.write_frame(get_exe_path()+"export/d/");
            lbm.graphics.set_camera_centered(0.0f, 90.0f, 100.0f, 1.100000f);
            lbm.graphics.write_frame(get_exe_path()+"export/e/");
            lbm.graphics.set_camera_free(float3(0.001612f*(float)Nx, 0.523852f*(float)Ny, 0.992613f*(float)Nz), 90.0f, 37.0f, 94.0f);
            lbm.graphics.write_frame(get_exe_path()+"export/f/");
        }
#endif // GRAPHICS && !INTERACTIVE_GRAPHICS
    }
}

[ProjectPhysX]

Hi @wjsjtu123,

when you load an .stl mesh, by default the (rotation) center of the mesh is set to the center of its bounding box. This is correct for 2- or 4-bladed propellers, but fails for 3- or 5-bladed propellers; their rotation center is offset from the bounding box center. Here an illustration:

You can figure out this offset with some geometry.

The solution is to manually correct the position of the mesh center, with

float3 offset = float3(...);
propeller->set_center(propeller->get_bounding_box_center()+offset);

The specified mesh center is used as the rotation center during lbm.voxelize_mesh_on_device(propeller, ...) and propeller->rotate(...).

Kind regards, Moritz