[Cschwatz]

Hello everyone,

I'm trying to write an UDF that implements the Gidaspow drag model in my simulation. The reason I'm coding it instead of using the one provided in Fluent is that I'm not working with a granular phase, meaning that this model becomes unavailable for me to select. I have been using the example given in the Fluent UDF manual as a base for my code. When I compile it in Fluent, it gives back no errors. However, when I try to run my simulation, it gives back a floating point error like this:

---
Divergence detected in AMG solver: pressure correction Stabilizing vof-1 to enhance linear solver robustness.

Divergence detected in AMG solver: vof-1
Error at host: floating point exception

Error at Node 0: floating point exception

Error: floating point exception
Error Object: #f
---

I'll attach what I wrote so far to see if it helps:

---
#include "udf.h"
#include "math.h"

DEFINE_EXCHANGE_PROPERTY(gidaspow_drag_model, c, mix_thread, s_col, f_col)
{
Thread *thread_p, *thread_rbc;
real x_vel_p, x_vel_rbc, y_vel_p, y_vel_rbc, abs_vel, slip_x, slip_y,
rho_p, rho_rbc, mu_p, Re_p, Cd,vof_p, vof_rbc, D_p_rbc,
C_rbc, diam_rbc;

/* Finding the threads for the gas (primary) */
/* and RBCs (secondary phases) */

thread_p = THREAD_SUB_THREAD(mix_thread, s_col);/* Gas phase */
thread_rbc = THREAD_SUB_THREAD(mix_thread, f_col);/* Solid phase*/

/* Finding phase velocities in x and y*/

x_vel_p = C_U(c, thread_p);
y_vel_p = C_V(c, thread_p);

x_vel_rbc = C_U(c, thread_rbc);
y_vel_rbc = C_V(c, thread_rbc);

/*Computing the slip between the phases*/
slip_x = x_vel_p - x_vel_rbc;
slip_y = y_vel_p - y_vel_rbc;

/*Calculating the properties*/
rho_p = C_R(c, thread_p);
rho_rbc = C_R(c, thread_rbc);

mu_p = C_MU_L(c, thread_p);

/*Defining the diameter of the discrete phase*/
diam_rbc = C_PHASE_DIAMETER(c, thread_rbc);

/*Compute slip*/
abs_vel = sqrt(slip_x * slip_x + slip_y * slip_y);

vof_p = C_VOF(c, thread_p);/* Gas volume fraction*/
vof_rbc = C_VOF(c, thread_rbc);/* RBC volume fraction*/

/*Computing Reynold's number*/

Re_p = (vof_p*rho_p*abs_vel*diam_rbc) / mu_p;

/*Computing the drag coefficients based on the RBC volume fraction*/

if (vof_rbc < 0.2) {

Cd = pow(vof_rbc, -1.65)*fmax((24. / Re_p)*(1 + 0.15*pow(Re_p, 0.687)), 0.44);

C_rbc = (3 / 4)*((Cd*vof_rbc*vof_p*rho_p*abs_vel) / diam_rbc);
}

else if (vof_rbc >= 0.2)
{

C_rbc = (150 * (pow((1 - vof_rbc), 2)*mu_p) / (vof_rbc*pow(diam_rbc, 2))) + (1.35*((1 - vof_rbc)*rho_p*abs_vel) / diam_rbc);
}

/*Computing the Drag*/

D_p_rbc = C_rbc * (slip_x + slip_y);

return D_p_rbc;
}

---

Thanks in advance!

[AlexanderZ]

if the problem is not in your code, than it should be model settings. Something with VOF.
Try to simplify your model and code and reach convergence

on the other hand you have

C_rbc = (150 * (pow((1 - vof_rbc), 2)*mu_p) C_rbc = (150 * (pow((1 - vof_rbc), 2)*mu_p) / (vof_rbc*pow(diam_rbc, 2))) + (1.35*((1 - vof_rbc)*rho_p*abs_vel) / diam_rbc); + (1.35*((1 - vof_rbc)*rho_p*abs_vel) / diam_rbc);

where you are dividing on vof_rbc, which could be 0 potencially, check it

best regards

[Cschwatz]

Thank you for the reply Alexander!

I guess the problem was the division by vof_rbc, since fluent initializes a case by setting the secondary volume fraction as being zero. I solved this problem by changing the solid volume fraction (vof_rbc) by (1 - fluid volume fraction), which is essentially the same thing, but it avoids possible divisions by zero.