Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Propeller rotation axis off-axis during voxelization #163

Closed
wjsjtu123 opened this issue Mar 15, 2024 · 1 comment
Closed

Propeller rotation axis off-axis during voxelization #163

wjsjtu123 opened this issue Mar 15, 2024 · 1 comment
Labels
setup question question for a specific setup

Comments

@wjsjtu123
Copy link

wjsjtu123 commented Mar 15, 2024

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?
image
image
image
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
Copy link
Owner

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:
image

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
setup question question for a specific setup
Projects
None yet
Development

No branches or pull requests

2 participants