Exhaust Manifold Boundary Conditions

[CbsTau]

Hello Dr. Lehmann and Community,

I continue to be amazed at FluidX3D! As an engineer who has worked with traditional CFD for over 25 years, I continue to be blown away by the speed and overall functionality of your tool! It is amazing what you have accomplished.

I am trying to simulate an engine exhaust manifold of a 4-cylinder engine. I am able to get the basic simulation to run with each cylinder firing in sequence, and repeating in a loop. The problem is that after a cylinder has fired, and injected its gas into the given branch of the manifold, that inlet needs to be "sealed" off. All four inlets need to be TYPE_E to allow the inflow of gas at the appropriate time, but when any branch is not actively firing, it needs to be closed with a reflective boundary TYPE_S.

You will see if you run the setup below, that as the 1st and 2nd cylinder fire in sequence, much of the flow goes out of cylinder 3 and 4, instead of out of the exhaust, because the inactive cylinders are effectively open, even though I tried to force their velocity = 0.0f.

Please can you explain if / how it is possible to change the Type of boundaries during the simulation from TYPE_E to TYPE_S, or how else to close off the three inlets which are not active for any given moment.

Many thanks in advance,
Clinton.

//  Exhaust Manifold - Four cylinders firing in sequence loop
//  
void main_setup() { // Flow thru an Exhaust Manifold; required extensions in defines.hpp: FP16S, EQUILIBRIUM_BOUNDARIES, SUBGRID, INTERACTIVE_GRAPHICS or GRAPHICS
	// ################################################################## define simulation box size, viscosity and volume force ###################################################################
	const uint3 lbm_N = resolution(float3(0.50f, 0.35f, 0.25f), 4000u); // input: simulation box aspect ratio and VRAM occupation in MB, output: grid resolution
	const float si_u = 2.0f;		// 5 m/s speed (y)
	const float si_length = 0.5f;	// 500mm length (x)
	const float si_T = 0.1f;
	const float si_nu = 1.48E-5f, si_rho = 1.225f;
	const float lbm_length = 1.0f * (float)lbm_N.x;		// set the lbm_length to 100% of box width (x)
	const float lbm_u = 0.12f;							// set lbm_u speed to 0.075 for solver stability
	units.set_m_kg_s(lbm_length, lbm_u, 1.0f, si_length, si_u, si_rho);
	const float lbm_nu = units.nu(si_nu);
	const ulong lbm_T = units.t(si_T);
	print_info("Re = " + to_string(to_uint(units.si_Re(si_length, si_u, si_nu))));
	LBM lbm(lbm_N, lbm_nu);
	// ###################################################################################### define geometry ######################################################################################
	Mesh* manifold = read_stl(get_exe_path() + "../stl/manifold_lbm.stl", lbm.size(), lbm.center(), lbm_length); //
	lbm.voxelize_mesh_on_device(manifold);

	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 (x == 0u || x == Nx - 1u || z == 0u || y == Ny - 1u) lbm.flags[n] = TYPE_S; // All closed box walls
	if (z == Nz - 1u) lbm.flags[n] = TYPE_E; // Pipe inlet

	if (y == 0u) lbm.flags[n] = TYPE_E; //  inlet velocity fixed

		}); // ####################################################################### run simulation, export images and data ##########################################################################
	lbm.graphics.visualization_modes = VIS_FLAG_LATTICE | VIS_FIELD | VIS_STREAMLINES;
	lbm.run(0u); // initialize simulation
	lbm.graphics.set_camera_centered(-60.0f, 40.0f, 100.0f, 1.000000f);
	while (true) { // main simulation loop
		lbm.u.read_from_device();
		const float cycle_time = 0.1f;        // Total cycle time in seconds for all four cylinders to fire (repeats every 0.1 s)
		const float active_duration = 0.025f; // Each cylinder is active for 0.025 seconds
		const float current_time = fmodf((float)lbm.get_t(), units.t(cycle_time)); // Time within 0.1s cycle

		// Initialize velocities for each cylinder as zero
		float uz1 = 0.0f, uz2 = 0.0f, uz3 = 0.0f, uz4 = 0.0f;
		
		// Activate the appropriate cylinder based on the current time within the cycle
		if (current_time < units.t(active_duration)) {
			uz1 = lbm_u;  // Cylinder 1 fires
		}
		else if (current_time < 2 * units.t(active_duration)) {
			uz2 = lbm_u;  // Cylinder 2 fires
		}
		else if (current_time < 3 * units.t(active_duration)) {
			uz3 = lbm_u;  // Cylinder 3 fires
		}
		else if (current_time < 4 * units.t(active_duration)) {
			uz4 = lbm_u;  // Cylinder 4 fires
		}

		// Define the bundary faces for each of the cylinders, to apply the inflow velocity
		// The four inflow holes are positioned in the middle of each quarter width in x
		uint z = Nz - 1u;

		for (uint y = 1u; y < Ny - 1u; y++) {
			for (uint x = 1u; x < 0.25f * Nx; x++) {
				const uint n = x + (y + z * Ny) * Nx;
				lbm.u.z[n] = -uz1;
			}
		}
		for (uint y = 1u; y < Ny - 1u; y++) {
			for (uint x = 0.25f * Nx; x < 0.5f * Nx; x++) {
				const uint n = x + (y + z * Ny) * Nx;
				lbm.u.z[n] = -uz2;
			}
		}
		for (uint y = 1u; y < Ny - 1u; y++) {
			for (uint x = 0.5f * Nx; x < 0.75f * Nx; x++) {
				const uint n = x + (y + z * Ny) * Nx;
				lbm.u.z[n] = -uz3;
			}
		}
		for (uint y = 1u; y < Ny - 1u; y++) {
			for (uint x = 0.75f * Nx; x < Nx - 1u; x++) {
				const uint n = x + (y + z * Ny) * Nx;
				lbm.u.z[n] = -uz4;
			}
		}

		lbm.u.write_to_device();
		lbm.run(100u);		// run for 100 steps before updating the veloctiy 
	}

#if defined(GRAPHICS) && !defined(INTERACTIVE_GRAPHICS)
	lbm.graphics.visualization_modes = VIS_FLAG_LATTICE | VIS_FIELD | VIS_STREAMLINES;
	lbm.run(0u); // initialize simulation
	lbm.graphics.set_camera_centered(-60.0f, 40.0f, 100.0f, 1.000000f);
#else // GRAPHICS && !INTERACTIVE_GRAPHICS
	lbm.run();
#endif // GRAPHICS && !INTERACTIVE_GRAPHICS
} /*

manifold_lbm.zip