Learning Forum

[vineethcns]

Hello, we are trying to simulate the experimental work related to phase change of CO2 in a converging diverging nozzle mentioned in the paper "Characterization of Non-equilibrium Condensation of Supercritical Carbon Dioxide in a de Laval Nozzle" by Lettieri et al.

Can anybody suggest how to get convergence with respect to liquid mass fraction in the droplet based non-equilibrium compressible phase change solver in Ansys CFX?

The Settings are:

Residual target : 1e-4

Conservation target : 1e-3

Liquid Condensation rate relaxation factor : 0.3

Nucleation bulk tension factor : 1.0

Surface tension coefficient: 0.01

[DrAmine]

Can you add more details? Have you check the NES tutorial?

[vineethcns]

Thank you very much dear Abenhadj for such a quick reply. Here are some of the details related to the CFX case

---

(1) We are using an axisymmetric domain for the nozzle (Nozzle boundary data from Lettieri et al. paper) with structured curvelinear mesh.

---

 
(2) We are starting from 55 bar initial pressure through the nozzle and slowly (in 2 milli second) increasing the pressure at inlet to 80 bar and decreasing the pressure at outlet to 28 bar.

Note: This way we are able to perfectly get convergence for the equilibrium solver case.

---

(3) The materials are defined using REFPROP based RGP tables. These tables work perfectly with the equilibrium solver.

 

---

(4) Transient solver setting details
    ================================

  ANALYSIS TYPE
    Option = Transient
    EXTERNAL SOLVER COUPLING
      Option = None
   
    INITIAL TIME
      Option = Automatic with Value
      Time = 0
   
    TIME DURATION
      Option = Total Time
      Total Time = 0.01
   

    TIME STEPS
      First Update Time = 0.0
      Initial Timestep = 1e-08
      Option = Adaptive
      Timestep Update Frequency = 1
      TIMESTEP ADAPTION
        Courant Number = 0.1
        Maximum Timestep = 1e-07
        Minimum Timestep = 1e-9
        Option = MAX Courant Number
     

---

  
(5) Boundary conditions
    =================== 

    DOMAIN Default Domain

    Domain Type = Fluid
    Location = fluid 1

    BOUNDARY inlet
   

---

      Boundary Type = INLET
      Location = inlet

      BOUNDARY CONDITIONS
        FLOW DIRECTION
          Option = Normal to Boundary Condition
       
        FLOW REGIME
          Option = Subsonic
       
        HEAT TRANSFER
          Option = Fluid Dependent
       
        MASS AND MOMENTUM
          Option = Total Pressure
          Relative Pressure = inletPressure
       
        TURBULENCE
          Option = Medium Intensity and Eddy Viscosity Ratio
       
     
      FLUID co2gas
        BOUNDARY CONDITIONS
          HEAT TRANSFER
            Option = Total Temperature
            Total Temperature = 311 [K]
         
          VOLUME FRACTION
            Option = Value
            Volume Fraction = 1.0
         
       
      FLUID co2liq
        BOUNDARY CONDITIONS
          DROPLET NUMBER
            Droplet Number = 0.0 [m^-3]
            Option = Specified Number
         
          VOLUME FRACTION
            Option = Value
            Volume Fraction = 0.0
         
       
  
   
    BOUNDARY outlet
   

---

      Boundary Type = OUTLET
      Coord Frame = Coord 0
      Location = outlet
      BOUNDARY CONDITIONS
        FLOW REGIME
          Option = Subsonic
       
        MASS AND MOMENTUM
          Option = Average Static Pressure
          Pressure Profile Bl = 0.05
          Relative Pressure = outletPressure
       
        PRESSURE AVERAGING
          Option = Average Over Whole Outlet
       
     
   
    BOUNDARY walls
   

---

      Boundary Type = WALL
      Location = walls
      BOUNDARY CONDITIONS
        HEAT TRANSFER
          Option = Adiabatic
       
        MASS AND MOMENTUM
          Option = No Slip Wall
       
        WALL ROUGHNESS
          Option = Smooth Wall
       
     
   
    BOUNDARY wedge0
   

---

      Boundary Type = SYMMETRY
      Location = wedge0
   
    BOUNDARY wedge1
   

---

      Boundary Type = SYMMETRY
      Location = wedge1
   
     
      REFERENCE PRESSURE
        Reference Pressure = 0 [bar]

---

(6) Multiphase solver settings
    ==========================


    MULTIPHASE MODELS
      Homogeneous Model = On
      FREE SURFACE MODEL
        Option = None
     
   
    FLUID DEFINITION co2gas
      Material = CO2v
      Option = Material Library
      MORPHOLOGY
        Option = Continuous Fluid
     
   
    FLUID DEFINITION co2liq
      Material = CO2l
      Option = Material Library
      MORPHOLOGY
        Condensation Rate Relaxation Factor = 0.3
        Option = Droplets with Phase Change
     
  
      FLUID co2gas
        HEAT TRANSFER MODEL
          Include Viscous Work Term = On
          Option = Total Energy
       
     
      FLUID co2liq
        HEAT TRANSFER MODEL
          Option = Small Droplet Temperature
       
        NUCLEATION MODEL
          Nucleation Bulk Tension Factor = 1.0
          Option = Homogeneous
       
     
      HEAT TRANSFER MODEL
        Homogeneous Model = Off
        Option = Fluid Depent
     
      THERMAL RADIATION MODEL
        Option = None
     
      TURBULENCE MODEL
        Option = SST
     
      TURBULENT WALL FUNCTIONS
        Option = Automatic
     
   
    FLUID PAIR co2gas | co2liq
      Surface Tension Coefficient = 0.01 [N m^-1]
      INTERPHASE HEAT TRANSFER
        Option = Small Droplets
     
      INTERPHASE TRANSFER MODEL
        Option = Particle Model
     
      MASS TRANSFER
        Option = Phase Change
        PHASE CHANGE MODEL
          Option = Small Droplets
       

---

 
(7) Initial conditions
    ==================
   
    INITIALISATION
      Coord Frame = Coord 0
      Option = Automatic

      FLUID co2gas
        INITIAL CONDITIONS
          TEMPERATURE
            Option = Automatic with Value
            Temperature = 311 [K]
         
          VOLUME FRACTION
            Option = Automatic with Value
            Volume Fraction = 1.0
         
      FLUID co2liq
        INITIAL CONDITIONS
          DROPLET NUMBER
            Droplet Number = 0 [m^-3]
            Option = Automatic with Value
         
          VOLUME FRACTION
            Option = Automatic with Value
            Volume Fraction = 0.0
         
      INITIAL CONDITIONS
        Velocity Type = Cartesian
        CARTESIAN VELOCITY COMPONENTS
          Option = Automatic with Value
          U = 0 [m s^-1]
          V = 0 [m s^-1]
          W = 0 [m s^-1]
       
        STATIC PRESSURE
          Option = Automatic with Value
          Relative Pressure = 55 [bar]
       
        TURBULENCE INITIAL CONDITIONS
          Option = Medium Intensity and Eddy Viscosity Ratio
       

---

   

( Miscellaneous settings
    ======================
   
  SOLVER CONTROL
    Turbulence Numerics = First Order
    ADVECTION SCHEME
      Option = High Resolution
   
    COMPRESSIBILITY CONTROL
      Clip Pressure for Properties = On
      High Speed Numerics = On
      Minimum Pressure for Properties = 6 [bar]
      Total Pressure Option = Automatic
   
    CONVERGENCE CONTROL
      Maximum Number of Coefficient Loops = 100
      Minimum Number of Coefficient Loops = 1
      Timescale Control = Coefficient Loops
   
    CONVERGENCE CRITERIA
      Conservation Target = 0.001
      Residual Target = 1e-4
      Residual Type = RMS
   
    INTERPOLATION SCHEME
      Pressure Interpolation Type = Trilinear
      Shape Function Option = Geometric
      Velocity Interpolation Type = Trilinear
   
    MULTIPHASE CONTROL
      Volume Fraction Coupling = Coupled
   
    TEMPERATURE DAMPING
      Option = Automatic
   
    TRANSIENT SCHEME
      Option = First Order Backward Euler
   
    VELOCITY PRESSURE COUPLING:
      Rhie Chow Option = High Resolution
   

---

   

Once again, thanks Abenhadj for the help. Please let me know in case more details are required.

[DrAmine]

Quite hard to provide help here: You said you are able to get a converged result with the equilibrium model. Was that run unsteady too? Can you run the NEQ for a fixed operating point using the steady state solver to get a good initial guess?

[vineethcns]

The equilibrium solver was also run under unsteady analysis. The results achieved almost steady state.

The steady state data will provide initial condition for pressure, velocity, liquid and vapour mass fraction but not for droplet number, droplet diameter etc.

Anyway I will try that...I suspect that may cause divergence as the initial condition for droplet number and diameter are unknown. We will take them to be zero and that may be far from steady state value...

Thanks again Sir...We will come back to you...

[DrAmine]

That is why to start with neq steady for operating point. You can assume at entrance no droplet like in the tutorial.

[vineethcns]

 

Dear Abenhadj,  

We tried the procedure as per your suggestions and we would like to thank you as we are getting better non-equilibrium condensation solver behaviour. However, there are still some issues which we would like you to guide us further. Please refer to the attached image:

1) We started with transient, equilibrium solver and then shifted to steady, equilibrium solver and got perfect convergence.  

2) Then we used the results from the steady, equilibrium solver as the initial conditions for the non-equilibrium solver.  

Things are converging but there are two issues:  

(a) We occasionally get some spikes in the residual for the liquid mass fraction. The Maximum coefficient loops for the transient solver was 100 here.  

My question here is that "Is this behaviour normal for this solver?"

(b) We decreased the Maximum coefficient loops to 40 to speed up the simulation but now the frequency of spikes increased.  

My question is that how many coefficient loops we have to take? Are there any guidelines?  

Thanks.

 

[DrAmine]

I won't play with the number of the coefficient loops in transient run. Better to reduce the time step size and use default coefficient loop number.

[vineethcns]

Dear Abenhadj

First of all thanks again for your continued support.

I have few more questions for you which I hope you will answer.

(1) For the non-equilibrium droplet based solver, We tried the adaptive time stepping control based on 'Number of Coefficient Loops' which were kept as minimum=1 and maximum=10. However, the time step taken by the solver was of the order of 1 nano second. It is really difficult to finish the simulation under such
a low time step.

 My question is: Is there a way to increase the time step? Will using second order Backward Euler for transient solution help in this case?

(2) I was going through similar work on CO2 condensation by Alireza Ameli et al. "Non-equilibrium condensation of supercritical carbon dioxide in a converging-diverging nozzle". They used a custom RGP file which had Spinodal limits built in for metastable regions. My question is that where I can get such an RGP file for CO2? Our University has Ansys Licence. Will the Ansys Customer support be able to help in this regard?

(3) My next question is: Does the non-equilibrium solver require custom RGP files with Spinodal limits built in for metastable regions?   Can we manage with just superheated and saturated data in the RGP table?

(4) My last question: In this formum once you discussed with somebody else about 'Subcooled Liquid' feature.
    Do we require the custom RGP file which had Spinodal limits built in for metastable regions for this feature to work?


Lastly: How can we acknowledge your help? Is there a way we can cite your help or Ansys?

Thank you very much again.

Regards

[DrAmine]

1/Probably: You might increase the courant number to a stable state. Regarding 2nd order for some quantities it is bounded scheme.

2/You can generate your RGP files using System Coupling (Earlier under AFD). Check the Fluent Beta Feature 2019R3.

3/Not obligatory

4/If not used, then CFX will assume saturated properties for the liquid even if the RGP contains both properties.

 

 

[ssyadav]

Dear Abenhadj

Vineeth CNS was my student, thank you very much for guiding us so far. We are able to perform some simulations now but the solver diverges suddenly.

We found that gradually increasing the pressure at inlet and decreasing at the outlet from the mean initial pressure is the best way to proceed. The solver diverges when phase change starts some where in the nozzle.

I was going through the paper "Non-equilibrium condensation of supercritical carbon dioxide in a converging-diverging nozzle" by Alireza Ameli which mentions that

"The Gibbs free energy change of the vapor near the critical point is very small which leads the equation
to have an infinite value and causes error in the solver. In order to prevent this numerical instability, the
solver clips the critical radius to a certain amount, as a default 1 mm. This extremely large default
clipping value decreases the accuracy of the simulation and results in non-realistic droplet sizes. The
problem is that it is not possible to change or modify the formula in the Ansys CFX by the user. To
sidestep this problem, the critical radius was calculated according to the equation (1) in the post
processing stage. It was observed that the critical radius, see Figure 4, changes very moderately also in
the nucleation zone. Therefore, the radius of one nano meter has been estimated to be realistic for the
critical droplet size in this study and in the calculations the critical radius is clipped, by a CEL, to one
nano meter."

We are also trying to clip the critical radius but do not know the variable name in for Liquid CO2 critical radius.

Please help us again regarding how we can limit the critical radius of the liquid phase to say 1 nano meter using CEL expression.

Thanks again for the help.

Regards

[ssyadav]

 

Dear Abenhadj 

The error message in the image above may give a clue regarding what may be happening. 

---

It says: FATAL Bounds error detected

             Variable: Liquid CO2 specific heat capacity at constant pressure

---

We are using "subcooled liquid" in the material definition by enabling beta features and the RGP file 

has the data corresponding to the subcooled region.

[DrAmine]

I need to read the post and update later. I think there is way but want to double check

[DrAmine]

And remember to contact the author of the paper for technical information

[DrAmine]

Provide a university mail address if you want to know how to limit the radius as it is supported.

[DrAmine]

Answered with expert parameter.

[ssyadav]

Dear Amine

I hope you are safe from the ongoing pandemic and pray that you remain so.

 I am sorry I am bugging you again and again.

I am stuck at another step, this time related to how to incorporate spinodal data into the RGP tables.

I simulated two different fluid flows inside nozzles.

Case 1) Flow of water in Moses  and Stein nozzle, Real Gas Data: IAPWS equation for water

Case 2) Flow of CO2 in another nozzle by Nakagwa, Real Gas Data: (a) RGP file for CO2 based on REFPROP, (b) Redlich Kwong EOS

I am following Ansys CFX guidelines which say that first switch off the nucleation process so that vapor super-cooling develops and then allow nucleation for the phase change process to occur.

The issue is: I am able to see significant water vapor super-cooling even after phase change in the nozzle with IAPWS data for water. Whereas, the super-cooling vanishes for CO2 once phase change starts. These happens with RGP files as well as with Redlich Kwong equation of state for CO2.  The figures attached above show things more clearly.

Now the CFX manual says that IAPWS library extends to metastable regions.

I am writing my own RGP file generator based on REFPROP.

I have two questions:

1) How to force CFX solver to read metastable data for CO2 when using Redlich Kwong equation of state. Any expert parameter for it?

2) Where I can find more information on the spinodal data to be put in RGP files?

Thanks once again for all the help till now.

Regards

Shyam

 

[ssyadav]

Dear Amine

I found this paragraph in the Ansys CFX solver modelling guide. Any help in this regard?

[DrAmine]

You should activate beta features and use use subcooled liquid properties for RGP files. Next week available again for any feedback

[DrAmine]

and refer to the documentation for creating RGP files out of REFPROP

[ssyadav]

Thanks Amine for the response...

The option "RGP Liquid Properties = Subcooled" applies to RGP Tables.

I could not find this option for Real Gas Equation of state like the Redlich Kwong EOS.

I think there must be an expert parameter for this option while using Real Gas Equation.

Regards

[DrAmine]

Yes that can be.

[DrAmine]

I cannot share the parameter here ( I know the parameter) and sincerely we think that a good RGP is enough to keep track of amlost everything. Cubic EOS are prone to errors for liquid calculations.

[ssyadav]

Thanks dear Amine

The problem I am facing with RGP file is that the solver diverges randomly even with time step as low as 1.0e-9  with errors mostly on 'fatal overflow in linear solver' and 'out of bound message for Cp'

It behaves well with Cubic equation of state and time step is also nice: 1.0e-7

That's why I wanted to know the parameter.

I am sorry I had to delete the post where I posted my mail id.  I got a threat message from some hacker...

Best regards

[DrAmine]

I would recommend to write the RGP file with the described method in the documentation. Regarding the parameter I will share per mail but please that will be the second and last exception...

[Leila Esfahanizadeh]

I am engaged in a project, which involved simulating the discharge and dispersion of supercritical CO2 through the nozzles of a cutting tool. 

The geometry comprises a nozzle with an approximate discharge diameter of 0.2 mm, an inlet boundary with supercritical CO2 at 90 bar and 311K, and an outlet at atmospheric conditions, where CO2 disperses into the air. The importance lies in the downstream conditions upon expansion, with a specific emphasis on solid CO2 volume fraction, CO2 mass fraction, and temperature distribution. Our simulation has thus far achieved convergence for equilibrium phase change modeling of CO2 discharge from the nozzle. We have utilized RGP tables for liquid CO2, vapor CO2, and saturation properties in their homogeneous binary mixtures. Furthermore, we've enabled Beta features in our simulations to account for subcooled liquid properties.

In the next step, for non-equilibrium phase change simulation, I also need to use the expert parameters discussed here. Could you please guide me on how I can manually control the critical radius of sub-cooled CO2 in the nucleation model, and an expert parameter for applying the RK equation of state to account for metastable states in the liquid phase? 

I also appreciate any help regarding the modeling structure in which, I can include Air in addition to the CO2 equilibrium binary mixture in the phase change simulation.

Best regards,

Leila