Adding drag force from the particles in NS equation for two-way coupling

[Abishek Sarkar]

Dear all,

I am trying to simulate a particle-laden, isotropic turbulent flow using Dedalus 3. The one-way coupled simulations have given me good results till now. Now I trying to incorporate two-way coupling, by considering the drag forces from the particles as body force in the fluid momentum equation. But the issue I am facing is that, when this additional body force term is added into the momentum equation, I am getting an abrupt flow like behavior towards the negative x-direction, which destroys the isotropy in the flow and affects the dissipation rate, which causes the other turbulent parameters to get modified. I will share the lines of code relevant to this. Kindly let me know, what changes are required to tackle this persistent issue.

particle_force_field = dist.VectorField(coords, name='particle_force', bases=(xbasis, ybasis, zbasis))

# Problem Setup
problem = d3.IVP([u, p, tau_p], namespace=locals())
problem.add_equation("dt(u) + (1/rho_f)*grad(p) - nu*lap(u) = - u@grad(u) - iBody*(1/rho_f)*particle_force_field + iForcing*Fw")
problem.add_equation("div(u) + tau_p = 0")
problem.add_equation("integ(p) = 0")  # Pressure gauge

def compute_force_field_delta(....):

    # Fine and coarse grid dimensions
    dx_d, dy_d, dz_d = Lx / Nx_d, Ly / Ny_d, Lz / Nz_d
    V_cell_d = dx_d * dy_d * dz_d

    # Reset fine-grid deposition field (high-res drag force array)
    particle_force_field_deposit.fill(0)

    # Particle loop: deposit drag force directly into fine grid
    for i, (x, y, z) in enumerate(particle_positions):
        fluid_velocity_x = U_x[i]   #same in  y and z direction

        particle_velocity_x = particle_velocities[i, 0]    # same in y and z direction

        drag_force_x = (1 / V_cell_d) * (3 * np.pi * rho_f * nu * d_p) * (fluid_velocity_x - particle_velocity_x) * body_force_switch    # same in y and z direction

        # Map particles to fine-grid cell indices using robust symmetric rounding
        xg = x / dx_d    #  same in y and z direction

        ix = int(np.floor(xg + 0.5)) % Nx_d    #  same in y and z direction


        # Deposit drag force into fine grid using np.add.at
        np.add.at(particle_force_field_deposit[0], (ix, iy, iz), drag_force_x)   #  same in y and z direction

    particle_force_field['g'][:] = particle_force_field_deposit

Thank you

SIncerely

Abishek

[Keaton Burns]

Hi Abishek,

There’s not really enough detail here to see what’s going on, but this isn’t really a built-in feature of Dedalus so I think you’ll have to mostly debug it yourself. E.g. try to build the force field you want, and then plot it to see what might be going wrong. This will be especially tricky in parallel when the grid is distributed, etc.  If you have any specific questions about Dedalus data layouts, etc., please feel free to follow up.

Best,

-Keaton

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/d5196199-71dc-4aff-ad14-1b7147456c7dn%40googlegroups.com.

[Abishek Sarkar]

Thanks Keaton for your speedy response.

Right now I am trying to deposit the drag forces from the particles, as delta deposition by approximating the indices of the particle's location and assigning the forces to the approximated grid point in the dealiased grid of particle_force_field['g'] array (which has of (3, 1.5*Nx, 1.5*Ny, 1.5*Nz)). I was hoping that after depositing the numerical values on the grid points,  Dedalus will convert (or interpolate) back to the particle_force_field vector which has a dimension of (3, Nx.Ny,Nz).

But it seems that this deposition of the drag forces on the dealiased grid is causing some numerical errors, resulting in a net flow like behavior in the negative x-direction. May be dedalus is not being able to generate the particle_force_field vector from the deposited values of the particle drag forces in the particle_force_field['g'] array.

Therefore, I kindly request you (and others as well), if there is a safe way to deposit the numerically calculated values of the particle drag forces on the particle_force_field['g'] array, so that it is able to correctly capture the delta like deposition of the particle drag forces. Or if there is a way to directly deposit the numerical values of drag forces on the particle_force_field vector.

Thank you

Sincerely

Abishek

Abishek Sarkar

PhD

Mechanical Engineering Department

IIT Gandhinagar

You received this message because you are subscribed to a topic in the Google Groups "Dedalus Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dedalus-users/aYk10aUmx4k/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/CADZXxBjUxzpphzuXQyFyAGAUM%2BXC3hig5rCgoyOKjSgzkRdSFg%40mail.gmail.com.

[Keaton Burns]

Hi Abishek,

Yes that is how Dedalus should interpolate if you explicitly change the scales on field. Again I suggest you try to debug this by examining the numpy arrays your code is producing to set the grid values, and then look at how the field transforms are modifying that with dealiasing, etc.

It’s probably also a good idea to put in some explicit smoothing, i.e. use narrow gaussians instead of delta functions for the particle impulses.

Best,

-Keaton

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/CAL1Oi9JBP%3Dymn1Udc89rf_BE6Qx9eP6ahoy5xbJP2M-7MeBvBw%40mail.gmail.com.