CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Fluent UDF and Scheme Programming

Ishii-Zuber drag in Fluent

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By CeesH

Reply
 
LinkBack Thread Tools Display Modes
Old   August 23, 2014, 05:11
Default Ishii-Zuber drag in Fluent
  #1
Senior Member
 
Cees Haringa
Join Date: May 2013
Location: Delft
Posts: 258
Rep Power: 6
CeesH is on a distinguished road
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.

Code:
#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;
}
CeesH is offline   Reply With Quote

Old   August 27, 2014, 09:43
Default
  #2
Sun
Member
 
Join Date: Nov 2010
Posts: 86
Rep Power: 6
Sun is on a distinguished road
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!
Sun is offline   Reply With Quote

Old   September 1, 2014, 05:08
Default
  #3
Senior Member
 
Cees Haringa
Join Date: May 2013
Location: Delft
Posts: 258
Rep Power: 6
CeesH is on a distinguished road
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..
CeesH is offline   Reply With Quote

Old   September 1, 2014, 05:23
Default
  #4
Sun
Member
 
Join Date: Nov 2010
Posts: 86
Rep Power: 6
Sun is on a distinguished road
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!
Sun is offline   Reply With Quote

Old   September 1, 2014, 07:28
Default
  #5
Senior Member
 
Cees Haringa
Join Date: May 2013
Location: Delft
Posts: 258
Rep Power: 6
CeesH is on a distinguished road
Part of the fun indeed
Sun likes this.
CeesH is offline   Reply With Quote

Old   September 2, 2014, 05:09
Default
  #6
Senior Member
 
Cees Haringa
Join Date: May 2013
Location: Delft
Posts: 258
Rep Power: 6
CeesH is on a distinguished road
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.
CeesH is offline   Reply With Quote

Old   September 2, 2014, 07:02
Default
  #7
Sun
Member
 
Join Date: Nov 2010
Posts: 86
Rep Power: 6
Sun is on a distinguished road
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).
Sun is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Calculating drag coefficient from Fluent Richard FLUENT 8 February 17, 2015 08:45
Problem in using parallel process in fluent 14 Tleja FLUENT 3 September 13, 2013 10:54
Drag force calculation in FLUENT pychear FLUENT 0 August 20, 2013 10:35
direction vectors for lift and drag in fluent nauman55 FLUENT 6 September 10, 2012 10:34
How does FLUENT calculate lift and drag? xTamx420 FLUENT 0 May 30, 2011 13:35


All times are GMT -4. The time now is 12:54.