Net flow in the -ve x-direction due to the body force from particles

[Abishek Sarkar]

Dear all,

I am trying to implement two-way coupling to track particles in a isotropic turbulent flow. To do that I am defining a function that evaluates the drag forces from each and every particle, and add this contribution of the drag forces from the particles as a body force in the fluid momentum equation. 

But I have observed that, as soon as I implement this two-way coupling, a flow is generated in the -ve x-direction, primarily in the form a jet like structure from the right end of the domain, and generating a flow towards the left direction. This net flow in the -ve x-direction is giving unphysical results, primarily affecting the settling velocity of the particles.

def compute_force_field_gaussian(particle_positions, particle_velocities, U_x, U_y, U_z, d_p, nu, particle_force_field, sigma, body_force_factor):
    forces = np.zeros((3, Nx, Ny, Nz))

    particle_force_field['g'].fill(0)

     for i, (x, y, z) in enumerate(particle_positions):

        ix = int(np.floor(x / dx)) % Nx
        iy = int(np.floor(y / dy)) % Ny
        iz = int(np.floor(z / dz)) % Nz

        # Particle-fluid slip velocity
        slip_x = U_x[i] - particle_velocities[i, 0]
        slip_y = U_y[i] - particle_velocities[i, 1]
        slip_z = U_z[i] - particle_velocities[i, 2]

        fx = -drag_coeff * slip_x
        fy = -drag_coeff * slip_y
        fz = -drag_coeff * slip_z

        #  Gaussian weights
        r2 = (x - xg)**2 + (y - yg)**2 + (z - zg)**2
        weights = np.exp(-r2 / (2 * sigma**2))
        weights /= np.sum(weights)  # Normalize per particle

        # Apply weighted deposition
        np.add.at(particle_force_field['g'][0], (gx, gy, gz), fx * weights)
        np.add.at(particle_force_field['g'][1], (gx, gy, gz), fy * weights)
        np.add.at(particle_force_field['g'][2], (gx, gy, gz), fz * weights)

        Downscale using Gaussian smoothing followed by interpolation

Kindly look into the issue, if you know what is happening, or this is a numerical error caused by approximating the indices of the forces array by using the np.floor function, to get the nearest index to the particle location (which may not be grid point).

Thank you in advance.

Sincerely

Abishek