[dinh]

Dear Foamer,
Taking a look to ErgunWenYu drag model for the lagrangian tracking, OF 2.3.0 with alphac > 0.8, I think something wrong there, where the factor (1 - alphac) disappears.
Furthermore, can someone explain why alphac appears in the denominator?

Best,

[cmigueis]

dinh, I'm stucked in the very same question, did you have sucess about this thread? The Ergun part seems correct, but I agree with you that something looks bad in the Wen and Yu part.

[lilinmin]

Originally Posted by dinh

Dear Foamer,
Taking a look to ErgunWenYu drag model for the lagrangian tracking, OF 2.3.0 with alphac > 0.8, I think something wrong there, where the factor (1 - alphac) disappears.
Furthermore, can someone explain why alphac appears in the denominator?

Best,

I think that is right because it returns the force not the coefficient, but I think it should be :
(mass/p.rho())*(150.0*(1.0 - alphac)+1.75*Re)*muc/(alphac*sqr(p.d()))
because alphac is divided with (alphac*sqr(p.d())).

does anyone find the same question?

[neiht]

I believe that the formula should be:
(mass/p.rho())*(150.0*(1.0 - alphac)/alphac+1.75*Re)*muc/(sqr(p.d()))

Because Sp = V*beta/(1.0 - alphac)

where beta = F*18*(1.0 - alphac)*alphac*muc/sqr(q.d())

Then there is no /alphac in the end in Ergun approximation.

When I run Single particle sedimentation in openfoam, the velocity of particle is smaller than the result by numerical analysis.

[Shah Akib Sarwar]

Originally Posted by neiht

I believe that the formula should be:
(mass/p.rho())*(150.0*(1.0 - alphac)/alphac+1.75*Re)*muc/(sqr(p.d()))

Because Sp = V*beta/(1.0 - alphac)

where beta = F*18*(1.0 - alphac)*alphac*muc/sqr(q.d())

Then there is no /alphac in the end in Ergun approximation.

When I run Single particle sedimentation in openfoam, the velocity of particle is smaller than the result by numerical analysis.

Can you explain where in the source code you found the definition of Sp and beta?

[joshwilliams]

1-\alpha_c should not appear in the Lagrangian formulation of the model (only the two-fluid model version of the drag law).
If 1-\alpha_c is there, then the drag force will decay to 0 as the solid volume fraction approaches zero (e.g. a single small isolated sphere). However, it should not happen. The drag force should decay to stokes drag.


I attached this figure showing drag force with varying \alpha_p and Re_p.


Script to reproduce is below

import numpy as np 
from numpy import pi
import matplotlib.pyplot as plt

def coeff_drag(Re):
# if Re > 1000.0:
#     return 0.44*Re
# else:
   return24.0*(1.0+0.15*Re**0.687)

def dragmodel_current(alpha_p, Re, d_p=10e-6, rho_p=1207, muc=1.5e-5):
    alpha_c =1-alpha_p
    vol_p =4./3. * pi * (d_p/2)**3
   return vol_p *0.75* coeff_drag(alpha_c*Re)*muc*pow(alpha_c, -2.65) / (alpha_c * np.square(d_p))

def dragmodel_literature(alpha_p, Re, d_p=10e-6, rho_p=1207, muc=1.5e-5):
    alpha_c =1-alpha_p
    vol_p =4./3. * pi * (d_p/2)**3
   return vol_p *0.75* alpha_p * coeff_drag(alpha_c*Re)*muc*pow(alpha_c, -2.65) / (alpha_c * np.square(d_p))

def dragmodel_oneway(alpha_p, Re, d_p=10e-6, rho_p=1207, muc=1.5e-5):
# alpha_p is unused
    num_vals = alpha_p.size
    alpha_p =None
    vol_p =4./3. * pi * (d_p/2)**3
   return vol_p *0.75* coeff_drag(Re)*muc / np.square(d_p) * np.ones(num_vals)

def stokes_drag(alpha_p, Re, d_p=10e-6, rho_p=1207, muc=1.5e-5):
# alpha_p is unused
    num_vals = alpha_p.size
    alpha_p =None
    Re =None
    vol_p =4./3. * pi * (d_p/2)**3
# return vol_p * 0.75 * coeff_drag(Re)*muc / np.square(d_p) * np.ones(alpha_p.size)
    Cd =24
   return vol_p *0.75* Cd * muc / np.square(d_p) * np.ones(num_vals)

alpha_p = np.array([0] + (10**np.linspace(-4, -1, 6)).tolist() + np.linspace(0.1, 0.6, 10).tolist())

Re =10
# print(dragmodel_current(0.1, Re=10), dragmodel_literature(0.1, Re=10))
ls_list = ["-", "--", "dotted", "-."]

fig, ax = plt.subplots()
ax.scatter(alpha_p, stokes_drag(alpha_p, Re), label=f"Stokes drag", c="black")
# plot drag force models at varying Reynolds numbers
for i, Re inenumerate([0, 0.1, 100]):
    ax.plot(alpha_p, dragmodel_current(alpha_p, Re), label=f"Wen-Yu (OF), Re = {Re}", ls=ls_list[i], c="orange", lw=2)
for i, Re inenumerate([0, 0.1, 100]):
    ax.plot(alpha_p, dragmodel_literature(alpha_p, Re), label=fr"Wen-Yu (with 1-$\alpha_f$), Re = {Re}", ls=ls_list[i], c="blue", lw=2)
for i, Re inenumerate([0, 0.1, 100]):
    ax.plot(alpha_p, dragmodel_oneway(alpha_p, Re), label=f"Schiller-Nauman, Re = {Re}", ls=ls_list[i], c="green", lw=2)
ax.set_yscale("log")
ax.set_xscale("log")
ax.set_ylabel(r"$\mathbf{f}_d$ [kg m / $s^2$]")
ax.set_xlabel(r"$\alpha_p$ [-]")
ax.set_title(r"Drag models at $d_p = 10^{-6} m$, $\rho_p = 1000$")

# format plot for legend outside of plot
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width *0.6, box.height])

# tidy up
for spine_i in ["top", "right"]:
    ax.spines[spine_i].set_visible(False)
ax.legend(frameon=False,loc="center left", bbox_to_anchor=(1.,0.5))
plt.savefig("drag_forces.png", dpi=300)
plt.show()