|
[Sponsors] |
August 23, 2014, 05:11 |
Ishii-Zuber drag in Fluent
|
#1 |
Senior Member
Cees Haringa
Join Date: May 2013
Location: Delft
Posts: 607
Rep Power: 0 |
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; } |
|
August 27, 2014, 09:43 |
|
#2 |
Senior Member
Join Date: Nov 2010
Posts: 103
Rep Power: 15 |
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! |
|
September 1, 2014, 05:08 |
|
#3 |
Senior Member
Cees Haringa
Join Date: May 2013
Location: Delft
Posts: 607
Rep Power: 0 |
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.. |
|
September 1, 2014, 05:23 |
|
#4 |
Senior Member
Join Date: Nov 2010
Posts: 103
Rep Power: 15 |
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! |
|
September 1, 2014, 07:28 |
|
#5 |
Senior Member
Cees Haringa
Join Date: May 2013
Location: Delft
Posts: 607
Rep Power: 0 |
Part of the fun indeed
|
|
September 2, 2014, 05:09 |
|
#6 |
Senior Member
Cees Haringa
Join Date: May 2013
Location: Delft
Posts: 607
Rep Power: 0 |
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.
|
|
September 2, 2014, 07:02 |
|
#7 |
Senior Member
Join Date: Nov 2010
Posts: 103
Rep Power: 15 |
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). |
|
March 16, 2017, 09:52 |
|
#8 |
New Member
Join Date: Jan 2016
Posts: 1
Rep Power: 0 |
CeesH - even though this is a somewhat old thread, do you mind sharing the final code with us, if you're still around
|
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
direction vectors for lift and drag in fluent | nauman55 | FLUENT | 7 | September 23, 2020 14:47 |
Calculating drag coefficient from Fluent | Richard | FLUENT | 12 | April 11, 2018 15:35 |
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 |
How does FLUENT calculate lift and drag? | xTamx420 | FLUENT | 0 | May 30, 2011 13:35 |