[CeesH]

Hi all,

Since for some reason ANSYS decided not to include Ishii-Zuber drag in FLUENT (eulerian multiphase, in this case); I'm attempting to put it in myself. I basically put the CFX-implementation of the law in an UDF, hooked it up to FLUENT and then... hell broke loose. In other words; divergence.

Of course, there were some bugs but I believe I got all of them now, yet the simulation still fails. I tried hooking from both steady and transient, from t = 0 of starting from a converged case with S-N drag, but nothing works. Of course, I tried fiddeling the under-relaxation factors. I also tried several meshes, mimicking a case that has been reported in literature as being done in FLUENT with I-Z drag, with an equally rough an with finer meshes, but no success.

When I output the values of C_D vs. Re_p outside of FLUENT, the results look exactly as expected: A C_D between 2.67 and 0.42, depending on size, in the Newton regime and a switch to viscous between Re 10 and 100, depending on size.

In most cases, FLUENT does iterate a few times before crashing - when outputting u_slip vs. C_D the results also look as expected. But, uslip goes up dramatically every iteration, hitting light-speed after some 10 iterations. C_D seems to be fine, but the adaption of velocity to it does not.

So, I'm lost. Below is my code. Does anybody see any mistakes? Does anybody have experience implementing drag models and knows the problem? Does anybody have a working implementation of I-Z drag for me to test/compare? All help would be awesome.

#include "udf.h"
#include "stdio.h"
#include "math.h"

#define GRA 9.81
#define ST 0.072

/*note it only works in 3d
gravity is 9.81, surface tension 0.072*/

DEFINE_EXCHANGE_PROPERTY(ishiizuber,c,mix_thread,s_col,f_col)
{
Thread *thread_l, *thread_g;
double C_D, bs, drho, frd, eotvos, mustar, mumix, rd, remix,slipx,slipy,slipz,uslip, cds, cdcap, cdel;

thread_l = THREAD_SUB_THREAD(mix_thread, s_col);/* continuous */
thread_g = THREAD_SUB_THREAD(mix_thread, f_col);/* discontinuo*/
bs = C_PHASE_DIAMETER(c,thread_g);
rd = C_VOF(c,thread_g);

/*calculate slip velocity*/
slipx = C_U(c, thread_g) - C_U(c, thread_l);
slipy = C_V(c, thread_g) - C_V(c, thread_l);
slipz = C_W(c, thread_g) - C_W(c, thread_l);
uslip = sqrt((slipx*slipx) + (slipy*slipy) + (slipz*slipz));

/*caluclate mixture viscosity and reynolds, and eotvos number*/
mustar = (C_MU_L(c,thread_g) + 0.4*C_MU_L(c,thread_l))/(C_MU_L(c,thread_g) + C_MU_L(c,thread_l));
mumix = C_MU_L(c,thread_l)*(pow((1.-rd),(-2.5*mustar)));
remix = ((C_R(c,thread_l)*bs)/mumix)*uslip;
drho = (C_R(c,thread_l) - C_R(c,thread_g));
eotvos = (drho*bs*bs* GRA)/ST;

/*cd for spherical bubbles*/
cds = (24./remix)*(1. + (0.15*pow(remix,0.687));

/*cd for elliptical bubbles*/
frd = (C_MU_L(c,thread_l)/mumix) * pow((1.-rd),0.5);
cdel = ((1. + 17.67*pow(frd,(6./7.)))/(18.67*frd)) * (2./3.) * (pow(eotvos,0.5));

/*cd for the cap regime*/
cdcap = (8./3.)*(1.-rd)*(1.-rd);

/*determine which coefficient to use*/
if(cdel < cdcap)
{
    C_D = cdel;
}
else
{
    C_D = cdcap;
}

if(cds > C_D)
{
    C_D = cds;
}

return C_D;
}

[Sun]

Cees,
I've done some implementations of Cds for stirred tank and I don't think anything is wrong with your code. But what i can think of is you are returning C_D (sounds weird!) whereas you need to return the "fluid-fluid exchange coefficient". Have a look at the theory guide and search for fluid-fluid exchange coefficient (Kpq or Klg if I remember) in the expression of Kpq there is a drag function for which you've computed C_D via I-Z. So, to put it in a nutshell calculate C_D with I-Z then calculate and return Kpq with the computed C_D.
Cheers!

[CeesH]

Hi Sun,

I completely overlooked that, thanks! Now it seems to be running well (fingers crossed - it hasn't converged yet..) without lightspeed velocities or so.

edit: Of course, pretty much as soon as I posted this, divergence. But I guess I'm going in the right direction now..

[Sun]

Sorry for the diverged solution, but I guess it's part of fun! have a look at the values of slip velocity (a few time steps before divergence) and see if it's going crazy high, if it is not the case you'll be sure that you are in the right direction.
And if I were you I would run it with Schiller-Naumann as well to see if the high velocity is only caused by the choice of C_D or other things might contribute.
cheers!

[CeesH]

Part of the fun indeed

[CeesH]

It seems the problem was that in the mixture reynolds number the term (1-VOF(gas)) ends up in the denominator, which leads to a singularity in regions where the volume fraction goes to 1. Using a capping term or avoiding mixture viscosity solves the problem.

[Sun]

Great!
glad to hear you've found/solved the problem.
you can also add an if for the denominator (if (1-rd!=0) cdcap = ...) or maybe use the max function like max((1-rd), 1e-7).

[stuffen17]

CeesH - even though this is a somewhat old thread, do you mind sharing the final code with us, if you're still around