Boat in water simulation setup

[antoniopucciarelli]

Dear Community / @ProjectPhysX ,

I am struggling with the setup of a boat moving through water in FluidX3D.

The issue i have is the following:

As you can see, the flow at the outlet of the domain is flowing in at the inlet. Do you know how to solve this problem?

The setup I use is the following:

void main_setup() { // stable boat

	const uint3 lbm_N = resolution(float3(1.0f, 2.0f, 1.0f), 1000u); // input: simulation box aspect ratio and VRAM occupation in MB, output: grid resolution
	
	const float si_u       = 15.0f;
	const float si_length  = 2.4f;
	const float si_nu      = 1.48E-5f;
	const float si_rho     = 1.225f;
	const float lbm_length = 0.65f * (float)lbm_N.y;
	const float lbm_u      = 0.05f;
	
	units.set_m_kg_s(lbm_length, lbm_u, 1.0f, si_length, si_u, si_rho);
	
	print_info("Re = " + to_string(to_uint(units.si_Re(si_length, si_u, si_nu))));
	
	// ################################################################## define simulation box size, viscosity and volume force ###################################################################
	LBM lbm(lbm_N, units.nu(si_nu));

	// ###################################################################################### define geometry ######################################################################################
	const uint  Nx = lbm.get_Nx();
	const uint  Ny = lbm.get_Ny();
	const uint  Nz = lbm.get_Nz();
	const float R  = 20.0f;
	const float H  = 90.0f;
	
	// const float lbm_length = 0.4f * (float)Nx;
	
	// mesh properties
	string pathName = "??/FluidX3D/stl/boat.stl";
	Mesh* boatMesh  = read_stl(pathName);

	const float scale = lbm_length / boatMesh->get_bounding_box_size().y / 3.5f; 
	boatMesh->scale(scale);
	
	const float3 offset = lbm.center() - boatMesh->get_bounding_box_center() - float3(0.0f, boatMesh->get_bounding_box_center().y / 2.0f, boatMesh->get_bounding_box_center().z * -0.15f);
	boatMesh->translate(offset);
	boatMesh->set_center(boatMesh->get_bounding_box_center());

	lbm.voxelize_mesh_on_device(boatMesh, TYPE_S|TYPE_X, offset, float3(0.0f), float3(0.0f, 0.0f, 0.0f));
 
	parallel_for(lbm.get_N(), [&](ulong n) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
		if(z == 0) 
		{
			lbm.flags[n] = TYPE_S;
			lbm.u.y[n]   = lbm_u;
		}
		else if(lbm.flags[n] != (TYPE_S|TYPE_X) && z < H) 
		{
			lbm.flags[n] = TYPE_F;
			lbm.u.y[n]   = lbm_u;
		}
		
		if(x==0u||x==Nx-1u) 
		{
			lbm.flags[n] = TYPE_S;
		}
	}); 
	// ####################################################################### run simulation, export images and data ##########################################################################
	lbm.graphics.visualization_modes = VIS_FLAG_LATTICE | (lbm.get_D() == 1u ? VIS_PHI_RAYTRACE : VIS_PHI_RASTERIZE);
	lbm.run(0u); // initialize simulation
	
	while(true) 
	{ 
		lbm.u.read_from_device();
	
		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

		// sum force over all boundary
		const float3 lbm_force = lbm.calculate_force_on_object(TYPE_S|TYPE_X);
		const float3 lbm_torque = lbm.calculate_torque_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);
		const float si_torque_x = units.si_M(lbm_torque.x); // torque in Newton * meters
		const float si_torque_y = units.si_M(lbm_torque.y);
		const float si_torque_z = units.si_M(lbm_torque.z);

		// print lift/drag force components
		print_info(">>> FORCE");
		print_info("\tlift force    = " + to_string(si_force_z) + " N");
		print_info("\tdrag force    = " + to_string(si_force_y) + " N");
		print_info("\tlateral force = " + to_string(si_force_x) + " N");
		print_info(">>> TORQUE");
		print_info("\troll torque = " + to_string(si_torque_x) + " Nm");
		print_info("\tpitch torque = " + to_string(si_torque_y) + " Nm");
		print_info("\tyaw torque = " + to_string(si_torque_z) + " Nm");

		// moving sphere
		lbm.u.write_to_device();
		lbm.run(100u);
	}
}

Attached there is the .stl geometry.
boat.zip

My final aim is to simulate the motion of the boat inside water taking into account forces and moments.

If you have some scripts dealing with conservation of momentum and angular momentum, feel free to share it.

Thank you for the support.

[Freekbaars]

hi
i'm currently running through the simulation. (now i myself only started with FluidX3D today and not that much C++ experience) but now i have a little problem it seems like out of nowhere all of a sudden some kind of ball enters the domain. this causes it to stop working completely.

[CbsTau]

Hello @antoniopucciarelli,

I have also struggled with the same thing, and managed to figure out a working solution for your model.

Firstly, by default, all boundaries are periodic unless you specifically set them to be something else like TYPE_S or TYPE_E.

I made a few changes to you setup as there were a couple of mistakes (like using the density and viscosity of air instead of water), but the important changes are to get correct boundary conditions when using free-surface flow. This is not intuitive, but using the "hydraulic jump" example, I was able to apply a similar approach to your model.

Here is the modified code which produces great results. I have tried to add comments to explain. Your Force extraction code at the end seems to work fine, and is not included below to keep things more clear.

void main_setup() { 
	const uint3 lbm_N = resolution(float3(1.0f, 2.0f, 1.0f), 1000u); // input: simulation box aspect ratio and VRAM 
	const float si_u = 15.0f;			// Speed of the water = 15m/s
	const float si_length = 2.4f;		// Length of the boat = 2.4m
	const float si_nu = 1.0E-6f;		// Viscosity of water
	const float si_rho = 1000.0f;		// Density of water
	const float si_g = 9.81f;			// Gravity
	const float si_T = 1.0f;			// Simulation Time = 1s

	const float lbm_length = 0.3f * (float)lbm_N.y;		// equivalent boat length in lbm units (30% of domain size in y)
	const float lbm_h = 0.5f * (float)lbm_N.z;			// height of the water level as % of Domain height in z
	const float lbm_u = 0.1f;							// equivalent water speed in lbm units - Set to 0.075f to 0.1f
	units.set_m_kg_s(lbm_length, lbm_u, 1.0f, si_length, si_u, si_rho);
	const float lbm_f = units.f(si_rho, si_g);			// Generate body force to account for gravity
	print_info("Re = " + to_string(to_uint(units.si_Re(si_length, si_u, si_nu))));

	// ################################################################## define simulation box size, viscosity and volume force ###################################################################
	//LBM lbm(lbm_N, units.nu(si_nu), 0.0f, 0.0f, -lbm_f);	// Using actual si_nu seems to cause numerical instabilities 
	LBM lbm(lbm_N, 0.007f, 0.0f, 0.0f, -lbm_f);				// lbm_nu = 0.007f seems the stable choice for most general water simulations

	// ###################################################################################### define geometry ######################################################################################

	// mesh properties
	string pathName = (get_exe_path() + "../stl/boat.stl");
	Mesh* boatMesh = read_stl(pathName);

	const float scale = lbm_length / boatMesh->get_bounding_box_size().y;
	boatMesh->scale(scale);

	const float3 offset = lbm.center() - boatMesh->get_bounding_box_center() - float3(0.0f, boatMesh->get_bounding_box_center().y / 2.0f, boatMesh->get_bounding_box_center().z * -0.15f);
	boatMesh->translate(offset);
	boatMesh->set_center(boatMesh->get_bounding_box_center());

	lbm.voxelize_mesh_on_device(boatMesh, TYPE_S | TYPE_X);

	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);

	// The boudary conditions for this simualtion are adapted from the "hydraulic jump" example in setup.cpp
	// But this approach seems to make more physical sense.

	const uint  P1 = 20u;

	if (z < lbm_h) {
		lbm.flags[n] = TYPE_F;		// Initialize the Fluid cells below the free surface level (lbm_h)
		lbm.rho[n] = units.rho_hydrostatic(lbm_f, (float)z, lbm_h);			// Initialize the water with hydrostatic pressure gradient
	}

	lbm.u.y[n] = lbm_u;		// initialize all fluid cells at the desired free-stream velocity (lbm_u)

	if (y < P1 && z >= lbm_h) lbm.flags[n] = TYPE_S;	// creates a cutout in the domain to control the water flowing in at the desired height
	//
	if (z == 0u) lbm.flags[n] = TYPE_S; // set bottom to Solid Wall
	if (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == Nz - 1u) lbm.flags[n] = TYPE_E; // all other walls set to Environment (TYPE_E)

	//  The 1st face of cells at y=0u are set to the desired inflow velocity
	if (y == 0u && z < lbm_h) {
		lbm.u.y[n] = lbm_u;
	}
		});

	// ####################################################################### run simulation, export images and data ##########################################################################
	lbm.graphics.visualization_modes = VIS_FLAG_SURFACE | VIS_PHI_RAYTRACE;
	lbm.run(0u); // initialize simulation

	lbm.run();
}

I have also opened a related issue #243 regarding the Raytracing with an STL object which you might like to follow.

I hope this helps.

Kind Regards,
Clinton.