Numerical Analyses of Transcritical Evaporation using Diffuse-Interface Method

Liquid fuel placed in a supercritical environment (with respect to the fuel’s critical properties) can undergo different evaporation processes that is determined by the droplet’s surface temperature T s in comparison with the critical mixing temperature T c . In the transcritical regime, the ambient condition is such that at initial time, T s < T c , a sharp liquid-vapor interface exists, and evaporation occurs at subcritical condition. Eventually, the surface temperature approaches the critical mixing temperature, and the evaporation process becomes supercritical, with no distinct liquid-vapor interface. Traditional methods for simulating this complex evaporation behavior employ separate formulations for the liquid and vapor phases during subcritical evaporation and a continuous formulation during supercritical evaporation. Here, we developed a holistic formulation for the simulation of transcritical evaporation that is valid for both sub- and super-critical evaporation regimes using a diffuse-interface method derived from Linear Gradient Theory. The formulation allows for direct modeling of surface tension and interface thickness during the process of interface thickening as T s → T c . The diffuse-interface model was employed in simulations of a nanoscale transcritical planar evaporation conﬁguration. Furthermore, the ability of the diffuse-interface model to resolve both macroscale ﬂuid behavior and nanoscale liquid-vapor interface structure was demonstrated.


Introduction
In modern rockets, gas turbines, and diesel engines, the operating conditions are often such that liquid fuel is injected into an ambient environment at pressures that exceed the critical pressure of the fuel. A liquid fuel droplet at subcritical temperature placed in such a supercritical pressure environment can undergo different evaporation processes that depend on the local thermodynamic state of the interface. If the droplet's surface temperature is less than the critical mixing temperature (T s < T c ) for the entire droplet's lifetime, the evaporation process is characterized as subcritical, where a sharp liquid-vapor interface is always present. If the ambient temperature is sufficiently high so that T s reaches T c during the droplet's lifetime, the evaporation process is transcritical. In this evaporation regime, the surface of the droplet initially exists as a sharp liquid-vapor interface. During the subsequent evaporation process, the liquidvapor interface heats up and becomes diffused as T s approaches T c . Eventually, when T s = T c , the droplet's surface ceases to exist as a sharp liquid-vapor interface and the system transitions to a supercritical diffusive evaporation regime. Compared to subcritical evaporation, there is relatively little understanding and modeling capability of transcritical evaporation configurations. Numerical simulations of droplet evaporation traditionally employ distinct formulations for the liquid and gas phases, along with matching conditions for density, partial density, and energy fluxes across the liquid-vapor interface [1]. This hybrid numerical model assumes an infinitely thin liquid-vapor interface and thus is only valid in subcritical regimes. To extend this hybrid numerical model to transcritical evaporation, Zhu & Aggarwal [2] proposed to switch to a continuous single-phase formulation when T s = T c is achieved to simulate the supercritical regime. Nevertheless, the proposed method by Zhu & Aggarwal is unable to capture the gradual thickening process of the interface. Dahms & Oefelein [3] resolved the steady liquid-vapor interface structure and demonstrated the thickening of the steady interface as the temperature approaches the critical mixing temperature. The interface structure reconstruction was achieved by employing the Linear Gradient Theory (LGT) to model interfacial properties. The use of LGT to model interfacial dynamics has also been employed by Yue et al. [4] in simulations of of viscoelastic droplets in a Newtonian matrix and by Shen et al. [5] in modeling the contact-line behavior of vapor bubble during boiling. Mo & Qiao [6] contributed to the understanding of the liquid-vapor interface behavior during transcritical evaporation using molecular dynamics simulation. Despite these recent achievements in understanding of the liquid-vapor interface structure at transcritical conditions, a holistic continuum formulation for simulation of transcritical evaporation that is able to resolve liquid-vapor interface thickening remains elusive. In this paper, we developed a diffuse-interface method for capturing the unsteady behavior of liquid-vapor interface during transcritical evaporation. Specifically, we employ an unsteady diffused-interface model derived from LGT [7,8] to simulate the nanoscale transcritical evaporation process of Dodecane/Nitrogen. The thermodynamic condition is such that the pressure is supercritical for both species and temperature is supercritical for Nitrogen. Furthermore, the ability of the diffuse-interface model to simulate engineer-scale evaporation processes was demonstrated.

Mathematical Model
Linear Gradient Theory, formulated by van der Waals and later by Cahn & Hilliard [9], models the molecular interaction effects inside the interface by augmenting thermodynamic potentials with a contribution from the interface that is dependent on gradients of partial densities. For example, in a binary mixture, is the augmented internal energy of a binary system, where GD stands for gradient-dependent and e is the specific internal energy of the binary mixture. ρ is the density and ρY i for i = 1, 2 are partial densities of the fuel species (in this case Dodecane), and ambient species (in this case Nitrogen), respectively. The interfacial energy contribution is scaled by the capillary coefficients κ ij , which, similar to surface tension, is dominantly a function of temperature. In this paper, we employ a diffuse-interface model that is derived from LGT [7,8]. The governing equations are Here, D/Dt = ∂/∂t + u · ∇ denotes material derivative, u is the flow velocity, P GD = P − 0.5 i j κ ij ∇ρY i · ∇ρY j is the gradient-dependent form of pressure, E = e GD + 0.5u · u is the total internal energy, τ = η(∇u + ∇u T ) − 2η/3(∇ · u)I is the the viscous stress tensor, K is Korteweg stress resulting from interfacial effects, J i and q are the diffusive species and heat fluxes, and J i and Q are the interfacial species and heat fluxes. For binary mixtures studied in this paper, the diffusive fuel species and heat fluxes J 1 and Q are The augmented interfacial stress and fluxes (K, J i , and Q) can be closed by assuming that the interfacial stress tensor is elastically restoring and considering that the total species and heat fluxes should be irreversible per the second law of thermodynamics. For binary mixture studied in this paper, the interfacial stress and fluxes are In this study, Soret and Dufour effects are neglected. Here, µ i , s i , and h i for i = 1, 2 are the partial chemical potential, entropy, and enthalpy, respectively. W i for i = 1, 2 are molecular weight of each species, R is the universal gas constant, and X i for i = 1, 2 are species molar fractions. β v is the volume expansivity, λ is the thermal conductivity, and D 12 is the binary diffusion coefficient. Furthermore, we compute the interfacial surface tension σ from the net under-pressure inside the interface: The thermodynamic variables are closed using the Peng-Robinson equation of state. The capillary coefficient κ ij are calculated using correlation from experimentally measured surface tension [10]. The method of Chung et al. [11] is used to calculated mixture viscosity and conductivity at elevated pressure, and the Chapman-Enskog theory [12], along with Takahashi's high pressure correction [13], is used to calculate binary diffusivity. The diffuse-interface equations are discretized on a 1-D grid using a central finite-difference scheme. We use a Backward-Differentiation-Formular (BDF) scheme for time integration to handle stiffness in the governing equations due to the interface.

Configurations
In this study, we perform 1-D simulations to study subcritical and transcritical planar evaporation in a nanoscale configuration. We impose symmetric condition at the liquid (left) boundary and Dirichlet boundary conditions for thermodynamic variables at the ambient (right) boundary. For both cases, the initial half-width of the liquid film is x 0 = 5 nm and the ambient condition is fixed at x amb = 100 nm. For both configurations, the operating pressure is supercritical at P amb = 60 bar and the initial liquid phase temperature is T l = 450 K. The ambient heating temperature is T amb = 600 K for the subcritical configuration and T amb = 800 K for the transcritical configuration.
To obtain a steady liquid-vapor interface at initial conditions, we first solve the governing diffuseinterface equations with isothermal condition T = T l . The starting point of this isothermal simulation is a tanh profile for density and partial densities with left and right conditions obtained from Vapor-Liquid Equilibrium (VLE) calculations.

Subcritical Evaporation
To begin the subcritical evaporation process, we apply a heating layer near the ambient boundary to bring the ambient temperature up to T amb = 600 K. This heating layer is imposed as a tanh layer for temperature located at x = 75 nm (see Figure. 1a). Initially, the fluid is in compositional equilibrium across the domain (see Figure. 1b). As a result, no evaporation occurs and the behavior of the liquid-vapor interface is largely dominated by the heating and expansion of the liquid film, causing the interface to move outward (see the upper panels of Figure. 1). As the interface heats up, it maintains quasi-steady state and exhibit different equilibrium VLE composition at its immediate surroundings. This behavior can be observed in the increase of the gas-phase Dodecane mass fraction in Figure. 1b. Thus, eventually, the diffusive species flux caused by compositional inequilibrium drives the subsequent evaporation process, causing the interface to move inward (see the bottom panels of Figure. 1. In conjunction, the heating rate of the Dodecane liquid film slows as some energy is contributed towards the enthalpy of vaporization. Figure. 2a demonstrates this movement of the interface during the evaporation process. Furthermore, Figure. 2b shows that the interface's thickness increases and surface tension decreases as interfacial temperature rises toward the critical mixing temperature. However, the liquid-vapor interface remains in existence throughout the evaporation process because interface temperature is always lower than the critical mixing temperature T c = 650K.

Transcritical Evaporation
To simulate the process of transcritical evaporation, we return to the isothermal steady interface profile and raise the ambient temperature to a supercritical value T amb = 800 K. Similar to subcritical evaporation, the interface initially moves outward due to heating and expansion of the liquid film. Eventually, a compositional inequilibrium develops, which drives the subsequent evaporation process, and the interface moves inward (see Figure. 3).
As the interfacial temperature approaches the critical mixing temperature, the liquid-vapor interface thickness increases and starts to diverge. Concurrently, the interfacial surface tension starts to vanish. At t = 14 ns, the interfacial temperature reaches the critical mixing temperature, the liquid-vapor interface ceases to exist, and the evaporation problem becomes diffusionally dominated. Figure. 5 shows a comparison of surface tension observed during subcritical and transcritical evaporation. In the subcritical evaporation configuration, the interfacial temperature is always below critical mixing temperature. Therefore, the liquid-vapor interface exists throughout the evaporation process and surface tension never vanishes. In the case of transcritical evaporation, the interface achieves critical mixing temperature at t = 14 ns and surface tension vanishes after this point. Furthermore, the agreement in surface tension in both subcritical and transcritical evaporation cases, despite the difference in transient behavior, demonstrates that the liquid-vapor interface remains in quasi-steady state during the evaporation process. The slight discrepancy between simulated values of surface tension and experimental data [14] can be attributed to the elevated pressure condition studied in this paper as well as the difference between the current binary Dodecane/Nitrogen system with experimental monocomponent Dodecane system.

Extension to Macroscale Configuration
To demonstrate the applicability of the diffuse-interface method in simulating multiphase problems at engineering scales, we set up a transcritical planar evaporation case in a macroscale domain. The initial liquid film half-width is x 0 = 0.05 mm and ambient condition is enforced at x amb = 1 mm. The pressure condition is P = 60 bar. The initial liquid temperature is T l = 450 K and the ambient gas temperature is T amb = 1500 K. The initial condition imposes composition equilibrium across the domain, constant temperature in the liquid phase, and constant temperature gradient in the gas phase. In order to resolve the nanoscale liquid-vapor interface structure, we use an initially stretched grid around x 0 = 0.05 mm. This allows us to discretize the diffuse-interface equations with 1200 grid points while achieving good resolution inside the interface, where the minimum grid spacing is 0.18 nm. During the macroscale movement of the interface as part of the evaporation process, we employ an adaptive grid refinement routine to coarsen and refine the grid accordingly so that both the grid resolution at the interface and the total grid size is maintained. The initial process of liquid film heat-up and expansion is simulated and presented in Figure. 6. We observe the same behavior as nanoscale planar evaporation, where the movement of the interface is at first dominated by the expansion of the liquid film. Meanwhile, an evaporation flux is developing as the heating up liquid-vapor interface induces compositional gradients in the liquid and gas phases. Note that by employing a stretched grid discretization, we are able to resolve both the macroscale fluid profile around the interface, where the liquid phase heats up and evaporates (see Figure. 6a), and the nanoscale interface structure, where the interface thickens as its temperature approaches critical mixing temperature (see Figure. 6b).

Conclusion
In this paper, we proposed a holistic formulation using diffuse-interface model for simulating transcritical evaporation processes. The formulation allows for natural modeling of surface tension and interface thickness during the process of interface thickening as T s → T c . The diffuse-interface formulation was first employed to simulate a nanoscale planar evaporation configuration. We demonstrated that the proposed continuum formulation was able to resolve the fluid behavior as previously observed using molecular dynamics, including the initial heatup and expansion of the liquid phase, the build-up of compositional gradient and subsequent evaporation process, as well as the thickening of the liquid-vapor interface and vanishing of the surface tension as interface temperature approached critical mixing state. The ability of the proposed diffuse-interface formulation in resolving both macroscale fluid behaviors around the interface and the nanoscale interface structure was demonstrated by simulation of an engineering-scale planar evaporation configuration. Computational feasibility was achieved while maintaining adequate grid resolution across the interface by employing a stretched grid discretization and an adaptive grid refinement routine. We showed that the proposed diffuse-interface formulation is a promising candidate for simulation of droplet evaporation at transcritical environments. The relative simplicity of this holistic formulation compared to traditional methods allows for the future extension to multicomponent reacting flows and the simulation of transcritical droplet combustion.