[Liufeng_ustb]

Hi friends. I need your help eagerly. If you are willing to help, please send your e-mail id to me, I can tell you more info, or I send you to my mesh and my case and my udf.

My e-mail is liufengustb@hotmail.com

I really, really, really need you guys' help.

I am a senior,my undergraduate thesis is to simulate a reduction process. I am stuck for a few days. I am worried about my thesis. If I have no results, there will be a trouble for me to graduate from college. So save me.

Problem Description:
1.The reactor is a cylindrical furnace.
2.The whole bottom is the inlet of the gas.
3.The reduction process can be described as follows: At first, put the solid in the reactor; Then the gas at a low speed goes through the layer of solid which can be called packed bed. The solid is some regular spheres. During the time when gas contacts with solid, the reduction reaction occurs as well as other transfer phenomenon.
4.Gas phase is a mixture which includes 3 species (h2, h2o (vapor), n2).
5.Solid phase is also a mixture which includes 2 species. (fe (s, iron), fe2o3 (s,hematite))
6.The reaction can be expressed: 3h2(g) + fe2o3(s) = 2fe(s) + 3h2o(g)
7.Initially, the mole fraction of h2o in gas phase is zero, the same as fe(s) in the solid phase.
Then, if the calculation runs as my expectation, my udf will be loaded. So, the mole fraction will change.

When I start to calculate my case which disables heterogeneous reaction (i.e I set Reaction Rate Function as none in Phase Interaction box), it ran successfully. However, once I loaded my UDF, it went wrong and a dialog told me it was of divergence.

Here are the details:

General

Models:
I disable the two additional Eulerian models (I have tried the situation of enable and disable, the results told me I should select disabling them). Besides, I imitated the tutorial of FLUENT about Eulerian Granular Model and Heat Transfer.

As for Species Model, I still don’t get clear. I am still not sure whether I should enable Reactions or not. Because in both two situations (enable Reactions or disable), FLUENT will occur divergence when UDF is loaded. So, at present, I temporarily disable Volumetric Reactions.


Materials:
you can check my solid.scm or case where I define my solid material as fluid.

Phases:
I don’t select Packed Bed when Granular is chosen.
I change Diameter into 0.01mm.
Other options are set according to the Tutorial. (Eulerian Granular Model and Heat Transfer)
Packing Limit is set 0.6, and I will set my volume of fraction of secondary phase 0.6 as well.

Boundary Conditions:
I summarize my experience and conclude that the boundary conditions are very important for the stability of solution, especially for Backflow which needs to be set in Pressure Outlet.

Velocity Inlet:
Primary phase (gas phase)
The velocity inlet is 2m/s.

Pressure Outlet
Mixture:
For gas, I only care Thermal Conditions and I ignore Species. I set Backflow Total Temperature as 1200K which is the same as inlet. But I keep Species default.

For Solid, it’s the same as gas phase. Only care Thermal Conditions.

Wall
Set thermal condition for mixture and momentum condition for both phases.

Solution Methods:
I only revise Transient Formulation as Second Order Implicit and keep Scheme as well Spatial Discretization default.

Solution Controls:
Pressure: 0.2; Momentum: 0.2; Species: all 0.95


Solution Initialization:
Compute from in.
As for patch, I patch several options for solid (X Velocity 0; Y Velocity 0; Temperature 800K; Granular Temperature: 1e-05; fe: 0; Volume Fraction: 0.6)


That's basic info. If you need more, please contact with me(send me e-mail).
I look forward your valuable suggestions.

Pleas help me. You do this is to save me.

Thanks you all. Forgive my long article.

All the best.

[RodriguezFatz]

This is a huge amount of information. Can I extract this?

When I start to calculate my case which disables heterogeneous reaction (i.e I set Reaction Rate Function as none in Phase Interaction box), it ran successfully. However, once I loaded my UDF, it went wrong and a dialog told me it was of divergence.

Meaning, your UDF doesn't work?

[Liufeng_ustb]

Originally Posted by RodriguezFatz

This is a huge amount of information. Can I extract this?

When I start to calculate my case which disables heterogeneous reaction (i.e I set Reaction Rate Function as none in Phase Interaction box), it ran successfully. However, once I loaded my UDF, it went wrong and a dialog told me it was of divergence.

Meaning, your UDF doesn't work?

I don't know why. But sometimes when I open contours where I can Heterogeneous Rate and Species, it shows that it reacts which means my udf is truely loaded.

However, it always is of divergence.

Do you have experience of using macro which defines heterogeneous reaction rate? If you do, do you have time to help me check my case and other documents? I really need help. Because this would influence my graduating from college.

My e-mail is liufengustb@hotmail.com. I can send more details to your e-mail...

Best wishes!

[RodriguezFatz]

The idea of a forum is, that others can read this stuff later on and get help for their own cases. If people help you privatly by mail it is pretty kind, but this doesn't get the forum better.

My question was: Do you only have the problem, when you switch on the UDF-reactions?

[Liufeng_ustb]

Originally Posted by RodriguezFatz

The idea of a forum is, that others can read this stuff later on and get help for their own cases. If people help you privatly by mail it is pretty kind, but this doesn't get the forum better.

My question was: Do you only have the problem, when you switch on the UDF-reactions?

Well, well. I am so sorry, I am too worried about my thesis, I just forgot what meaning the forum is. Forgive me.

/*This UDF aims to define the reaction rate requied to simulate the reduction of pellet. 
Here, gas is a primary phase mixture of three species: h2o and h2 and n2. 
solid phase is also a mixture which consists of fe2o3 and its reduction product fe. */

#include "udf.h"
#include "sg_mphase.h"

/*user inputs and some constants*/

#define MAX_SPE_EQNS_PRIM 3 
/*total number of species in primary phase*/

#define index_reduction_primary 0  
/*reduction species index in primary phase*/

#define  prim_index 0  
/*index of primary phase*/

#define  P_OPER 101325 
/*operating pressure equal to GUI value*/

#define diam2 0.01
/*the diameter of the secondary phase particle*/

#define Sc_no 0.7 
/*set Schmidt Number as constant*/

/*end of user inputs*/



/*************************************************************/

/* Computes equilibrium pressure of hydrogen as function of temperature */
/* Equation is taken from THERMODYNAMIC THEORY IN SI, by Gibbs Equilibrium */
/* Returns pressure in PASCALS, given temperature in KELVIN  */
/* fe2o3 is reduced step by step, fe2o3->fe3o4->feo, and the final step is most difficult. 
So the whole reaction rate is determined by the last step, that is 
                feo+h2=fe+h2o                                */
double p_equilibrium_h2(double t_solid, double p_gas)                                 

{
   double delt_g = 10040 - 5.18 * t_solid;
   double interm = -delt_g / 8.314 / t_solid;
   double p_equilibrium_h2 = p_gas / (1+exp(interm));
   return p_equilibrium_h2;
} 


     
DEFINE_HET_RXN_RATE(diff_rate, c, t, hr, mw, yi, rr, rr_t)
{
	Domain **domain_reactant = hr->domain_reactant;
	real *stoich_reactant = hr->stoich_reactant;
	int *reactant = hr->reactant;

	Thread **pt = THREAD_SUB_THREADS(t);
    Thread *tp = pt[0];
    Thread *ts = pt[1];
	/*These two sentences corresponds primary phase and secondary phase 
	which are defined in FLUENT to them defined in UDF*/

    int i; 
    real concentration_reduction_primary, accum = 0., mole_frac_reduction_prim,
         concentration_reduction_particle_surface;
	/*Due to the rate of reduction which is related to the difference of h2 pressure between
	in gas phase and at the surface of particle, we need to find the value of h2 pressure or concentration in 
	these two different places*/

    real T_prim = C_T(c, tp); /*primary phase (gas) temperature*/
    real T_sec = C_T(c, ts);  /*secondary phase (particle) temperature*/

    real D_coeff, Re, Sc, Nu, urel, urelx,urely,urelz=0., tpressure, mass_coeff, area_density, 
        flux_reduction ;



	D_coeff = C_MU_L(c,tp) / C_R(c,tp) / Sc_no;

    if(Data_Valid_P())                          
    {
      urelx = C_U(c,tp) - C_U(c,ts);
      urely = C_V(c,tp) - C_V(c,ts);
      
      urel = sqrt(urelx*urelx + urely*urely);
      /*relative velocity*/
	  
       Re = urel * diam2 * C_R(c,tp) / C_MU_L(c,tp);
	   /*Reynold number*/

       Nu =  2. + 0.6 * pow(Re, 0.5)* pow(Sc_no, 0.333);
		/*Nusselt number*/

       mass_coeff = Nu * D_coeff / diam2 ; 

	   tpressure = P_OPER +  C_P(c, tp);

       for (i=0; i < MAX_SPE_EQNS_PRIM ; i++)
         {
           accum = accum + C_YI(c, tp, i) / mw[i][prim_index];
		   /*The function of variable accum is to acquire the total moles 
		   in primary phase if given a 1kg primary phase*/
         } 

       mole_frac_reduction_prim = C_YI(c, tp, index_reduction_primary) 
                   / mw[index_reduction_primary][prim_index] / accum;
	   /*According the meaning of accum, we can readily know this formula can get the 
	   mole fraction of reduction specie in primary phase if given 1kg primary phase*/
	  

       concentration_reduction_primary = mole_frac_reduction_prim * tpressure 
                   / 8.314 / T_prim ;
	   /*This formula is an application of P=cRT, besides, Dalton's law is also used for 
	   p(h2) in primary phase*/

       concentration_reduction_particle_surface = p_equilibrium_h2(T_sec, tpressure) /8.314 / T_sec ;
	   /*This formula is much easier because at the surface of particle, the h2 pressure is 
	   exactly the equilibrium pressure, but we need to notice the temperature is the temperature of 
	   secondary phase*/

       area_density = 6. * C_VOF(c,ts) / diam2 ;
	   /*area_density is a variable representing the area of the particle of secondary phase in a cell
	   which also includes primary phase. It will be used to get rr*/

       flux_reduction = mass_coeff * sqrt((concentration_reduction_primary - concentration_reduction_particle_surface) * 
		   (concentration_reduction_primary - concentration_reduction_particle_surface));
	   /*flux_reduction is rate of transporting mass, its theoretical base is obvious for us*/

       *rr = area_density * flux_reduction ;        
	   /* *rr is a pointer, its entry is exactly the reaction rate for heterogeneous reaction*/

       }
}

That is my udf. I already know the reaction will cause solution more unstable. But, some people already worked out fluided bed whose heterogeneous reaction rate is defined by udf. I cannot find out why mine doesn't work well.
Because if my udf is not loaded, everything is ok.

All the best.

[Liufeng_ustb]

1.The reactor is a cylindrical furnace. (diameter: 0.05m; height: 0.3m)

2.The whole bottom is the inlet of the gas.

3.The reduction process can be described as follows: At first, put the solid in the reactor; Then the gas at a low speed goes through the layer of solid which can be called packed bed (its height is 0.1m). The solid is some regular spheres. During the time when gas contacts with solid, the reduction reaction occurs as well as other transfer phenomenon.

4.Gas phase is a mixture which includes 3 species (h2, h2o (vapor), n2).
5.Solid phase is also a mixture which includes 2 species. (fe (s, iron), fe2o3 (s,hematite))

6.The reaction can be expressed: 3h2(g) + fe2o3(s) = 2fe(s) + 3h2o(g)

7.Initially, the mole fraction of h2o in gas phase is zero, the same as fe(s) in the solid phase.

Then, if the calculation runs as my expectation, my udf will be loaded. So, the mole fraction will change.

[RodriguezFatz]

Did you decrease the time-step radically and get the same divergence?

[Liufeng_ustb]

Originally Posted by RodriguezFatz

Did you decrease the time-step radically and get the same divergence?

Yeah, I decreased it to 0.0002s. But it still didn't work. Should I continue to decrease it?

I double if my udf is something wrong, but I can not find it

[RodriguezFatz]

1) The appropriate time step size depends on the time scale of your problem. If the whole reaction will occour in 0.001s, a time step of 0.0002s will of course cause problems. Did you try to determine the time scale?

2) For debugging the udf: I would calculate the initial reaction rates by hand and just hardwire these numbers into the UDF, instead of calculate them like you do. If that works, you can try to spot the error in your UDF by stepwise switching on lines / modules of your UDF...

[Liufeng_ustb]

Originally Posted by RodriguezFatz

1) The appropriate time step size depends on the time scale of your problem. If the whole reaction will occour in 0.001s, a time step of 0.0002s will of course cause problems. Did you try to determine the time scale?

2) For debugging the udf: I would calculate the initial reaction rates by hand and just hardwire these numbers into the UDF, instead of calculate them like you do. If that works, you can try to spot the error in your UDF by stepwise switching on lines / modules of your UDF...

I will try it right now.

However, I have another questions. The smaller time step size I set, does it mean the more suitable it is (except long time running)? Is time step size related to or dependent on mesh quality?

Besides, after I get reaction rate by hand, how could I put it into my udf?
I think it(*rr) must be a number, does it matter if I directly set it to a constant number?

Thank you for your advice, I will calculate my reaction rate.

[RodriguezFatz]

The smaller the time step the more stable is your calculation. Time step size is related to the size of your grid cells. I don't think it is strongly related to the quality of the mesh.

I don't know how you set the reactions, I never did that in Fluent. But what you write sounds logical to me. And yes, what I mean is to set it to a constant number instead of letting the UDF do the calculation. Of course this will be only correct for the first time step, but you can at least try to spot the error source.

[thanhndb]

Dear all Fluent Users,

Even I do not success yet in writing a UDF for my problem. I have also asked several users but I did not get any replies so far.

I have just found a nice code that might be helpful for us.
Science is belonged to the world not for any single person.

If you are interested in, you can find the code attached here.

Enjoy your research work on CFD simulation of Multiphase processes!