The problem with calculating the lift of the propellers

[nbvtim2008]

Hello, you have created a wonderful program, but could you tell me the code for calculating the lift and torque of the propeller, when I try, the thrust is either 0, or NaN, or Inf, and also please tell me how to set a certain rotation speed, for example 1000 rad/s

[Dick-Tator]

Hi, I have a similar problem. I created a setup which should simualate a rotor and return its thrust. I can simulate the turbulence just fine but when I add the part, where the thrust should be calculate the simulation looks realy strange with white lines.

I also can't find the results of the Force Field operation, altough I copied that line from the example with the Ahmed body.
If this is helpful in some way let me know. And if somebody can help me, please tell me.

void main_setup() { // Propeller simulation; required extensions in defines.hpp: FP16S, EQUILIBRIUM_BOUNDARIES, MOVING_BOUNDARIES, SUBGRID, INTERACTIVE_GRAPHICS or GRAPHICS, Force Field 
	// ################################################################## define simulation box size, viscosity and volume force ###################################################################
	const uint memory = 5000u;
	const uint3 lbm_N = resolution(float3(1.0f, 1.5f, 1.0f), 5000u); // input: simulation box aspect ratio and VRAM occupation in MB, output: grid resolution
	const float lbm_Re = 1000000.0f;
	const float lbm_u = 0.1f;
	const uint lbm_T = 180000u;
	const uint lbm_dt = 4u;
	LBM lbm(lbm_N, units.nu_from_Re(lbm_Re, (float)lbm_N.x, lbm_u));
	// ###################################################################################### define geometry ######################################################################################
	const float3 center = lbm.center();
	const float3x3 rotation = float3x3(float3(0, 0, 1), radians(180.0f));
	Mesh* rotor = read_stl(get_exe_path()+"../stl/Standart3BPropv1.stl", 1.0f, rotation); // https://www.thingiverse.com/thing:3014759/files
	const float scale = 0.98f * rotor->get_scale_for_box_fit(lbm.size()); //rotor to simulation box size
	rotor->scale(scale*0.7);
	rotor->translate(lbm.center()-rotor->get_bounding_box_center()-float3(0.0f, 0.41f*rotor->get_max_size(), 0.0f));
	rotor->set_center(rotor->get_bounding_box_center());
	const float lbm_radius=0.5f*rotor->get_max_size(), omega=lbm_u/lbm_radius, domega=omega*(float)lbm_dt;
	//lbm.voxelize_mesh_on_device(stator, TYPE_S, center);
	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]==0u) lbm.u.y[n] = 0.3f*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
	}); 
	const string path = get_exe_path() + "FP16S/" + to_string(memory) + "MB/";//copied from Ahmed Body
	// ####################################################################### run simulation, export images and data ##########################################################################
	lbm.graphics.visualization_modes = VIS_FLAG_LATTICE|VIS_FLAG_SURFACE|VIS_Q_CRITERION;
	lbm.run(0u); // initialize simulation
	while(lbm.get_t()<lbm_T) { // main simulation loop
		lbm.voxelize_mesh_on_device(rotor, TYPE_S|TYPE_X, center, float3(0.0f), float3(0.0f, omega, 0.0f));
		lbm.run(lbm_dt);
		rotor->rotate(float3x3(float3(0.0f, 1.0f, 0.0f), -domega)); // rotate mesh
		// Force calculation
		lbm.calculate_force_on_boundaries();
		lbm.F.read_from_device();
		const float3 lbm_force = lbm.calculate_force_on_object(TYPE_S | TYPE_X);
		const float si_force_x = units.si_F(lbm_force.x);
		write_line(path+"si_force_x.dat", to_string(lbm.get_t())+"\t"+to_string(si_force_x, 3u)+"\n");//copied from Ahmed Body

#if defined(GRAPHICS) && !defined(INTERACTIVE_GRAPHICS)
		if(lbm.graphics.next_frame(lbm_T, 30.0f)) {
			lbm.graphics.set_camera_centered(-70.0f+100.0f*(float)lbm.get_t()/(float)lbm_T, 2.0f, 60.0f, 1.284025f);
			lbm.graphics.write_frame();
		}
#endif // GRAPHICS && !INTERACTIVE_GRAPHICS
	}
}

[ProjectPhysX]

Hi @nbvtim2008 and @Dick-Tator,

I don't think you can get usable forces on revoxelized moving geometry. The revoxelization with velocity boundary conditions puts nearby cells far off equilibrium, and their DDFs relax only in the few time steps before the next revoxelization. The interval between revoxelizations is not a physical parameter (has to be chosen somewhat small though, 1-10 time steps), and the forces through momentum exchange algorithm will depend on this interval and where within the interval you read them.

For thrust measurement, it's better to not measure forces on the propeller, but average flux (from velocity field) through an area behind the propeller.

Kind regards,
Moritz