General Guidelines for Simulating Evaporation of an Exotic Fluid

General Guidelines for Simulating Evaporation of an Exotic Fluid

Director of Engineering

Overview

While implementing a Volume Of Fluid (VOF) approach, a customer simulating the evaporation of an exotic fluid with highly non-linear properties reported a floating-point error, as illustrated in Figure 1.
      
Floating point error messageFigure 1: Floating point error message

A floating-point error (FPE) like this one can occur for many different reasons, as discussed in the below list, making it difficult to pinpoint the exact source of instability. It is important to understand and control CFD simulation methodologies to enhance stability and reach a converged solution. This article discusses a couple sources of instability and the associated modeling recommendations to mitigate instabilities and enhance simulation speed-up when simulating evaporation, as discussed in the below list.
  1. Domain considerations: Reducing the simulation domain as much as possible without simplifying the problem statement is critical when troubleshooting a FPE and establishing simulation methodology because this decreases turnaround time, increasing simulation speed-up, and reducing computation cost.
  2. Large gradients: In this case, evaporation (without boiling) is a slow process where the heating is due to the rising temperature of the boundaries/surface walls with high temperature gradients. If gradients are too large, the solver may not be able to resolve them properly. Therefore, mesh refinement, a sufficiently small timestep, and appropriate boundary/initial conditions are all critical considerations. 
    1. Boundary/initial conditions: Unrealistic boundary conditions can lead to physical results and instability in the simulation.
    2. Mesh refinement: Simulation accuracy and stability depends on how well the liquid/gas interface is resolved as there can be large gradients in the concentration of vapor phase. Larger temperature gradients at the wall to adjacent fluid can be another source of mesh refinement related instability.
  3. Physics selection: VOF is adopted as the multiphase model. It solves a single set of transport equations together with volume fraction for the phases and can be thought of as a single phase plus volume fraction. Multiphase interaction models are required to model evaporation/condensation, interface momentum dissipation, and multiphase material. The Antoine equation is adopted as empirical formula that correlates the vapor pressure of a liquid with its temperature using a customer field function.
  4. Numerical scheme: Obtaining a solution in a CFD simulation is an iterative process. The segregated solver was adopted for this simulation. Under-relaxation factors (URF) are the most important parameters for this type of solver. Adjusting these parameters can enhance simulation stability. Note, implicit multi-step for VOF is discussed in this article for simulation speed-up.
    1. Time-step size: The time-step must be sufficiently small enough to capture a sharp liquid/gas interface, but large enough to minimize total calculation time – by maximizing the time-step. Basically, the time-step must be limited so that the volume created or destroyed in a cell is never larger than a fraction of the cell volume.

Background

Evaporation is a form of vaporization, where the liquid phase changes to gas at the liquid-gas interface that undergoes simultaneous heat and mass transfer processes. Liquid molecules with sufficiently high kinetic energy diffuse through the liquid-gas interface into the unsaturated gas. During evaporation, the gas directly above the surface of the liquid is saturated with moisture via molecular diffusion; in the absence of gas movement, it spreads further into the bulk gas by diffusion. Note, condensation is the reverse process of vaporization (phase change from gas to liquid).

The case discussed in this article considers evaporation via natural convection where the mass transfer occurs due to movement of air caused by the difference in density of gas at the liquid surface and the gas above it. The exotic liquid phase has a highly non-linear and temperature dependent dynamic viscosity, as defined by a table. Saturation pressure is defined using the Antoine equation. Constant material properties are adopted for the multi-component gas.

Computational Domain Considerations

A simplified test case consisting of a three-dimensional cylindrical flask with a diameter of 0.05 m and height of 0.1 m is filled with an exotic liquid (Eulerian phase) to a height of 0.02 m is adopted in this example, as shown in Figure 2. The remainder of the flask volume is considered as a multi-component gas (Eulerian phase) consisting of air and exotic vapor. There are two issues with this initial approach – 3D and size. Instead, simplifying the domain to a 2D problem with a reduction of size in both axial and radial directions will significantly speed-up the simulation, as shown in Figure 3.

Initial computational domain
Figure 2: Initial domain, 3D, D = 0.05 m, H = 0.1 m, and l_h = 0.02 m

Simplified computational domain
Figure 3: Simplified computational domain, 2D, D = 0.025 m, H = 0.05 m, and l_h = 0.01 m

Boundary and Initial Conditions

The original model specified a constant elevated temperature thermal boundary condition for the vertical and bottom walls, as shown in Figure 4. This is problematic because it causes simulation instability from sharp temperature gradient between the wall surfaces and exotic liquid region. We can enhance stability by replacing the constant elevated temperature with a slow stepwise temperature ramp over an extended time, as shown in Figure 5. A stepwise approach will allow the energy to propagate into the bulk fluid region over a period, thereby reducing the temperature gradient between the heated wall and adjacent fluid. Also, changing the vertical walls to adiabatic and applying the temperature ramp to the bottom surface will remove any sharp thermal gradients between the vertical wall surface and the liquid-gas interface. 
Wall thermal boundary condition
Figure 4: Wall thermal boundary condition – constant specified elevated temperature
Thermal boundary conditions
Figure 5: Thermal boundary conditions - adiabatic vertical wall and bottom surface set to slow stepwise temperature ramp over an extended time

Mass Fraction

When the partial pressure of vapor in the liquid equals the saturation pressure of the liquid at the mixture's temperature, the mixture is termed saturated. Saturated air constitutes a blend of dry air and saturated liquid vapor. The quantity of liquid vapor within moist air varies, ranging from absent in dry air to a maximum level, contingent on the pressure and temperature, when the mixture attains saturation. Alternatively, it can be characterized concerning relative humidity, which is defined as the ratio of the mole fraction of liquid vapor in a given air sample to the mole fraction in a saturated moist air sample under identical temperature and pressure conditions. To facilitate applications involving evaporation and condensation, it is important to convert relative humidity (RH) into mass fraction at the appropriate inlet and outlet boundaries.

For a given RH, reference absolute pressure, and temperature calculate the following:

  1. Calculate the liquid saturation pressure based on the Antoine equation, as described in the section called Liquid Phase Saturation Pressure Definition- 
    Liquid phase saturation pressure
  2. Calculate liquid vapor partial pressure- 
    Liquid vapor partial pressure
  3. If you assume pressure is constant, then the component partial pressure is the same as the mole fraction from the ideal gas law, PV=nRT.  Calculate liquid vapor mole fraction- 
    liquid vapor mole fraction
  4. Calculate mixture molar weight- 
    mixture molar weight
  5. Calculate liquid vapor mass fraction-
    liquid vapor mass fraction
  6. Calculate dry air mass fraction
    dry air mass fraction

Enter the liquid vapor and dry air mass fractions into the appropriate inlet and outlet boundaries.

Mesh Considerations

Mesh refinement and quality are major considerations when dealing with highly unstable simulations. They can be dependent on the application and selected multiphase models. In general, the relative error of the evaporation rate in quasi-steady state simulation is proportional to the inverse of the magnitude of the gradient of volume fraction (or, interface thickness), as stated by equation (2632). So, “the accuracy depends on how well the location of the interface is determined in comparison with the boundary layer”. In the core mesh region, a sufficiently fine mesh is required to resolve a sharp gas-liquid interface, as shown in Figure 6 and 7. This can be achieved by using surface remeshing to improve mesh quality. Basically, the sharper the interface, the more stable the results.

If you decide to implement prism layers to help stabilize your simulation and/or set the URFs to a small value (i.e., 0.1), as discussed in the next section. Then it is a possible indication that the mixing encountered at the near-wall turbulent boundary layer results to multiphase structures (bubbles/droplets) which are not adequately resolved; the smeared/diffused interface of these structures may result to overestimation of the evaporation/condensation rates. Also, the y+ treatment models will have an influence on the approximation of the turbulent properties near the wall, and considering also that in typical engineering grids the near-wall bubbles/droplets are under resolved - then this modelling approach (high y+) may be beneficial, as shown in Figure 8. Basically, a coarser prism layer configuration with high y+ value wall treatment may work better in such cases to avoid divergence as opposed to creating a significant number of prism layers with a slow growth rate. 


Initial Mesh
Figure 6: Gas Volume Fraction, Initial coarse mesh

Refined mesh
Figure 7: Gas Volume Fraction, refined mesh

Coarse prism layers
Figure 8: Gas Volume Fraction, refined mesh with coarse prism layers (2) and high y+ wall treatment

Physics Selection

VOF is adopted as the multiphase model for this problem. It is well suited for simulating flows with immiscible fluids and capable of resolving the interface between the phases of the mixture. 

Flow Solver and Convection Scheme

It is important to understand and control the solver settings to ensure best convergence and enhance simulation stability. STAR-CCM+ has two flow solvers – segregated and coupled. The segregated solver is based on the fluid equations in a decoupled or semi-decoupled way, which means not solving all equations at the same time. In this example, the Segregated Flow solver is automatically selected when running simulations with VOF and multiphase model, as shown in Figure 9.

VOF Evaporation Model Selection
Figure 9: Segregated Flow Solver selected for VOF with multiphase model.

To obtain a sharp interface between liquid and vapor phases, the convection scheme is set to High-Resolution Interface Capturing (HRIC) with a Face Density Reconstruction set to 2nd order. Note, 1st-order convection scheme is more stable but less accurate. HRIC and Modified HRIC are more accurate, but less stable. HRIC is considered more accurate than modified HRIC.
  1. HRIC - requires fine mesh CFL < 1 and small timestep (more accurate than modified HRIC).
  2. Modified HRIC - ok if CFL > 1, not tied to small timestep or fine mesh.

Multiphase Interactions

Phase interaction models must be considered to calculate the impact of one phase on another phase. Evaporation/Condensation is selected as a required model and Interface Momentum Dissipation is selected as an optional model, as shown in Figure 10. Evaporation/Condensation model specifically addresses evaporation without boiling and condensation occurring at the interface between a liquid and gas phase, as part of the VOF method. Interface Momentum Dissipation option introduces additional momentum dissipation in the proximity of the free surface to dissipate parasitic currents. The momentum dissipation is proportional to the velocity gradient and the interface artificial viscosity. 

Multiphase interactions
Figure 10 Multiphase Interaction, Evaporation/Condensation and Interface Momentum Dissipation

Liquid Phase - Saturation Pressure Definition

Evaporation/Condensation introduces a dynamic interface between liquid and vapor phases. Accurate prediction of the phase change rate heavily relies on correct estimation of saturation pressure. Failing to incorporate such effects can lead to erroneous predictions of flow patterns, heat transfer rates, and species transport, impacting the overall simulation reliability and predictive capability. 

The Antoine equation is a widely used empirical formula that correlates the vapor pressure of a liquid with its temperature. It is important to incorporate the Antoine equation into your CFD model because it enables the simulation to reflect changes in saturation vapor pressure accurately as the temperature varies, capturing the interplay between heat transfer, mass transfer, and phase change phenomena.

Antoine Equation:
Antoine equation
Where,
P_sat = liquid saturation pressure (bar)
P_atm = atmospheric (absolute) pressure (bar)
T = temperature (K)
A, B, C = component-specific constants

Coefficients are typically tabulated for various materials with the temperature, T, being an absolute temperature in K and pressure in bars. You may need to modify your values to correspond to these units or convert thereafter using a custom field function (assign saturation pressure to the custom field function). 
Note, you may implement the Antoine equation using the built-in Antoine equation method or a custom field function. If you are having stability problems with the built-in Antoine equation method, then try using a custom field function.

Solver Settings

Under-relaxation factors (URF)

The most important parameters for the segregated solve are the under-relaxation factors (URF). The URF values can range between 0 (slow) and 1 (fast) and represent the rate at which the simulation will update from iteration to iteration. URF is a very useful device for nonlinear problems. It is often adopted to enhance simulation stability in the iterative solution of strongly nonlinear problems. The default values represent a good starting point, but these values may need to be modified if you have strong nonlinear behavior in your problem. Basically, you can lower to URF to slow the change from iteration to iteration which enhances stability at the expense of slower solution convergence rate. 

Lower URF to slow the solver progress for the following under solvers: segregated flow, VOF, and energy as well as evaporation under Continua > Models > Multiphase Interaction >Phase Interaction 1 > Models > Evaporation/Condensation. Decrease the values to 0.3 or 0.4. Further URF troubleshooting may be approached by identifying which re which residual equation is diverging, then reduce the associated URF.

VOF implicit multi-step

Implicit Multi-step is a VOF solution strategy that decouples the flow timestep from the volume fraction timestep. Thus, enabling two timesteps within the simulations. The flow timestep may be larger because it is not limited to a strict CFL condition. Volume fraction timestep is constrained by the CFL to achieve a sharp phase interface. Therefore, the volume fraction timestep is sub-stepped within the larger flow timestep without losing accuracy. This facilitates simulation speed-up. Increase VoF implicit multistep under Solvers > Segregated Multi-Step > Number of Steps. You may select between 1 to 20 sub-steps with a default value of 3. This example increases the number of sub-steps to 4, as shown in Figure 11.
      
Implicit Multi-step
Figure 11: VOF Implicit Multi-Step

This allows you to decouple fluid timestep from volume fraction timestep, and sub-step volume fraction timestep within fluid timestep. This facilitates simulation speedup and was previously implemented in the last model provided.

Time-Step Size Selection

Further Considerations

In multiphase simulations, it is very important to ensure conservation of mass for each time step and phase. Monitoring the mass conservation error allows you to assess the quality of the solution, identify the sources or error in mass conservation, and define stopping criteria based on a specific conservation target. This is in addition to the other simulation convergence checks, i.e., residuals and engineering parameters of interest.  Follow these step-by-step instructions to monitor phase mass conservation:
  1. Right click Reports. Select New > Flow/Energy > Mass Conservation Error
    Mass Conservation Error
  2. Rename the Report to “Mass Conservation Error Vapor”.
    Rename report
  3. In the properties window select Vapor phase.
    Vapor phase mass conservation error
  4. Right-click the report and select Create Monitor and Plot from Report. 
    Create monitor and plot from report
  5. Right-click Reports. Select New > User > Statistics Report
    statistics report
  6. Rename the report to Cumulative Mass Conservation Error Gas.
    cumulative mass conservation error vapor
  7. In the properties window select Units - kg, Statistics – Sum, Sample Collection policy – All Samples, Monitor – Mass Conservation Error Vapor Monitor.
    Monitor properties
  8. Right-click the report and select Create Monitor and Plot from Report.
  9. Repeat for all phases.