{'fraction_non_alphanumeric': 0.038862626790507636, 'fraction_numerical': 0.01147350167470366, 'mean_word_length': 4.74795193883124, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Self-similar solution of a supercritical two-phase laminar mixing layer

Previous works for a liquid suddenly contacting a gas at a supercritical pressure show the coexistence of both phases and the generation of diffusion layers on both sides of the liquid-gas interface due to thermodynamic phase equilibrium. A related numerical study of a laminar mixing layer between the liquid and gas streams in the near field of the splitter plate suggests that mass, momentum and thermal diffusion layers evolve in a self-similar manner at very high pressures. In this paper, the high-pressure, two-phase, laminar mixing-layer equations are recast in terms of a similarity variable. A liquid hydrocarbon and gaseous oxygen are considered. Freestream conditions and proper matching conditions at the liquid-gas interface are applied. To solve the system of equations, a real-fluid thermodynamic model based on the Soave-Redlich-Kwong equation of state is selected. A comparison with results obtained by directly solving the laminar mixing-layer equations shows the validity of the similarity approach applied to non-ideal two-phase flows. Even when the gas is hotter than the liquid, condensation can occur at high pressures while heat conducts into the liquid. Finally, a generalized correlation is proposed to represent the evolution of the mixing layer thickness for different problem setups.


Introduction
Combustion chambers used in many engineering applications (e.g., power units, gas turbines or liquidpropellant rockets) are designed to operate at elevated pressures. That is, a thermodynamic regime is sought where less dissociation of the reaction products occurs and a better combustion efficiency and specific energy conversion are obtained. In many situations, the chamber pressure is well above the critical pressure of the liquid fuel that is being injected. Well-known fuels are based on hydrocarbon mixtures (e.g., diesel fuel, Jet A, RP-1) with critical pressures around 20 bar, while operating pressures may range from 25 to 40 bar in diesel engines or gas turbines and 70 to 200 bar in rocket engines.
The performance of the combustion reaction also depends on how fast the liquid fuel vaporizes and mixes with the surrounding oxidizer. Therefore, understanding the physical phenomena involved in this process is necessary for a proper design of the injectors' shape and distribution, combustion chamber size, etc. At subcritical pressures, a clear distinction exists between liquid and gas phases, which allows for extensive experimental studies. However, for near-critical or supercritical pressures, experimental studies show that a thermodynamic transition occurs where the liquid and gas present similar fluid properties near the liquid-gas interface, which is suddenly immersed in a variable-density layer and is rapidly affected by turbulence [1][2][3][4][5][6][7][8][9].
Past works have described this behavior as a very fast transition of the liquid phase to a supercritical state, assuming that a two-phase behavior cannot be sustained under these thermodynamic conditions as suggested by experimental results where a gas-like turbulent structure is observed [10,11]. But evidence of a two-phase behavior at supercritical pressures in multicomponent fluids exists based on a requirement of thermodynamic phase equilibrium at the liquid-gas interface [12][13][14][15][16][17]. At supercritical pressures, phase equilibrium enhances the dissolution of the gas into the liquid phase, thus generating diffusion layers on both sides of the liquid-gas interface where strong variations of fluid properties occur. Mixture critical properties near the interface differ from pure fluid critical properties and, in general, the new critical pressure is higher than the chamber pressure. However, this does not seem to be a requirement to sustain two phases under very high pressures as discussed in Poblador-Ibanez and Sirignano [17]. Some authors suggest caution at pressures near the critical pressure of the mixture, where the liquid-gas interface enters a continuum region of a few nanometers thickness [18,19]. Accordingly, a two-phase behavior cannot exist. Nevertheless, diffusion around the interface occurs rapidly enough to reach diffusion layer thicknesses of the order of micrometers while phase equilibrium is well established [17]. Thus, the interface may still be treated as a discontinuity under phase equilibrium with a clear jump in fluid properties across it, similar to how a compressive shock wave is considered as a discontinuity in compressible fluids. In that case, the nonequilibrium layer thickness for the shock is an order of magnitude or more greater than the layer thickness for phase non-equilibrium.
There are further reasons that support the fact that liquid injection at supercritical pressures is still a two-phase problem. Surface tension forces are reduced at very high pressures due to both phases presenting similar densities at the interface. Therefore, the aerodynamic effects on the liquid breakup process are enhanced [20]. Moreover, mixing causes the liquid viscosity to drop to gas-like values near the interface, reducing viscous damping of surface instabilities [21]. Altogether, fast growing instabilities associated with smaller wavelengths develop at the interface, which may cause a very fast atomization of the liquid jet with a cloud of very small droplets around it. This reasoning is further supported by numerical results of incompressible liquid round jets and planar sheets injected into a gas. For similar liquid and gas densities with reduced surface tension, a fast breakup process is induced with early small droplet formation and enhanced radial development of the two-phase mixture [22][23][24][25][26]. Traditional experimental techniques have difficulties with visual penetration of dense sprays due to the scattering caused by the large amount of small droplets and the intrinsic difficulties of variable-density fluids. Thus, it becomes difficult to distinguish between liquid and gas phases. Some new techniques have been able to penetrate dense sprays and show clear liquid structures (i.e., ligaments, lobes, bigger droplets) emerging from the liquid core, although they have not been tested yet in supercritical environments [27,28].
The present work does not address the hydrodynamic instabilities that occur farther downstream as part of the atomization process. It focuses on the initial laminar mixing between the injected liquid and the gas before substantial growth of surface instabilities or transition to turbulence. Understanding how the mixing between both streams evolves is crucial to understand the beginning of high-pressure atomization and can also help create models to use as initial conditions for more complex numerical studies.
Davis et al. [29] present a two-dimensional numerical study involving the non-ideal two-phase laminar mixing layer equations and provide guidelines on the validity of these relations. Results are presented for a binary mixture where a liquid stream and a gas stream come together at the end of a splitter plate at various ambient pressures, ranging from subcritical to supercritical for the pure liquid. The liquid stream is initially composed by pure n-decane and the gas stream is pure oxygen. The mixing layer thickness grows to a few micrometers in the liquid phase and to tens of micrometers in the gas phase with downstream distance from the end of the splitter plate. At the same time, the interface reaches a near-steady solution sufficiently far away from the splitter plate. Most importantly, each individual configuration (i.e., each set of boundary conditions) shows signs of a self-similar behavior of the mixing layer evolution, even when the defined mixing layer equations involve non-ideal terms and are coupled to a non-ideal thermodynamic model involving a real-gas equation of state, thermodynamic fundamental principles and various high-pressure correlations. A self-similar behavior in this high-pressure environment is also suggested in Poblador-Ibanez and Sirignano [17], although in that work a non-ideal unsteady one-dimensional model is implemented.
Reducing the system of partial differential equations (i.e., mixing layer equations) to a system of ordinary differential equations in terms of a similarity variable is of special interest, not only because a simplified mathematical model describing the same phenomena is obtained, but also because of the fundamental analysis implied in the transformation. Similarity solutions have always been sought in classical fluid mechanics (e.g., the Blasius solution for flow over a flat plate or the Falker-Skan wedge flows), usually limited to incompressible flows. Some transformations are also available for compressible flows or flows with variable fluid properties, although simplifications are made to obtain a generic self-similar model (e.g., use of ideal-gas model). In general, these approaches are well documented in textbooks [30,31]. When analyzing compressible or variable-density laminar shear layers, the focus is usually limited to single-phase configurations under ideal-gas assumptions with momentum and thermal mixing. However, some works also include species mixing in their analysis [32].
In views of the discussed numerical results [17,29], the non-ideal two-phase laminar mixing layer equations are recast in terms of a similarity variable, thus forming a system of ordinary differential equations.
Boundary conditions and liquid-gas interface matching relations are provided, which allow for mass, momentum and energy transfer across the interface to occur. Then, the self-similar model is coupled to a non-ideal thermodynamic model to close the system of equations. The self-similar solution is compared against the results from Davis et al. [29] with different problem configurations to assess the validity of the self-similar model. Finally, a proper non-dimensional scaling is presented, which provides a semi-collapse of the self-similar solutions for different sets of boundary conditions. This allows the development of a correlation based on freestream and interface conditions which is able to estimate mixing layer thicknesses with sufficient accuracy.

Non-ideal two-phase laminar mixing layer equations
The behavior of two parallel streams coming together at the edge of a splitter plate can be modeled under the boundary-layer approximation after a transitional region [30]. For a sufficiently large Reynolds number, the flow begins to develop with the transverse velocity being much smaller than the streamwise velocity, v u, and variations in the x-direction become negligible compared to variations in y (i.e., ∂()/∂x ∂()/∂y).
The transverse pressure gradient, ∂p/∂y, is also negligible. Therefore, the transverse velocity solution is directly obtained from the continuity equation and not from the transverse momentum equation [17,33].
For free-shear, low-Mach number flows at high pressures without confining walls, pressure can be assumed constant (Dp/Dt = ∂p/∂x = ∂p/∂y = 0). Furthermore, viscous dissipation and kinetic energy can be safely neglected in the energy equation. Figure 1: Sketch of the mixing layer problem between a liquid stream and a gas stream. The liquid-gas interface or dividing streamline is assumed to be fixed at y = 0 m as discussed in [29].
Both fluid streams are initially pure components (e.g., a pure liquid hydrocarbon and a pure gas, such as oxygen where ρ, u and v are the mixture density, streamwise velocity and transverse velocity, respectively. h is the mixture specific enthalpy and h 1 and h 2 are the partial specific enthalpies of each species of the binary mixture. Other mixture fluid properties are the mixture dynamic viscosity, µ, the mass diffusion coefficient, D, the thermal conductivity, λ, and the specific heat at constant pressure, c p . Furthermore, the energy equation is written as an enthalpy transport equation by using the relation λ∇T = (λ/c p )∇h − N i=1 (λ/c p )h i ∇Y i . Note that since the problem is diffusion-driven, fluid properties will vary as species and energy diffuse.
These steady-state equations are valid as long as the flow is laminar and surface instabilities are negligible. Davis et al. [29] provide some guidance for oxygen-hydrocarbon mixtures at various ambient pressures.
The quasi-parallel laminar region may exist up to a distance from the splitter plate of the order of O(10 −2 m). The mean flow velocity is fixed at 10 m/s and the velocity difference between both streams is varied in the range of ∆u = 0.3 − 5 m/s to limit the maximum Reynolds number atL = 0.01 m to ReL = 10, 000. ReL is defined as [29] where the subscript G refers to gas freestream values. Analogously, the subscript L will refer to liquid freestream values.
With this analysis, one of the main goals is to show the existence of a distinct two-phase behavior before hydrodynamic instabilities dominate.

Similarity transformation
The system of partial differential equations (PDEs) describing the problem (Eqs. (1), (2), (3) and (4)) can be transformed into a system of ordinary differential equations (ODEs) by using a proper similarity variable, η. However, the transformation requires an intermediate step due to variable density. The mapping follows The first transformation is defined as where y is a dummy variable for integration purposes and z is a density-weighted transverse coordinate following the so-called Howarth-Dorodnitsyn transformation [31]. This type of transformation was independently introduced in 1942 by Dorodnitsyn [34] and later in 1948 by Howarth [35]. Other authors followed with similar work [36,37].
As shown in Figure 1, the interface between the two phases is placed at y = 0. To transform the original system of equations from (x, y) → (x, z), the following relations between partial derivatives are obtained Defining the transformed transverse velocity, w, as the velocity components can be related to the stream function for compressible flow, Ψ(x, z), as Here, w does not have units of velocity and can be considered as a convenient variable in the solution process.
This arrangement yields a system of PDEs in the (x, z) space in terms of the compressible stream function, Eqs. (11)- (14). Note the continuity equation, Eq. (11), is identically satisfied.
After this intermediate step, the self-similar transformation is taken by defining the similarity variable, η, as in the Blasius solution for a flat plate [30], but where x is substituted byx and y is substituted by z (see Eq. (15)).
The transformation from (x, z) → η requires the following relations between partial derivatives wherex =x = x.
A self-similar solution exists if Ψ/ √x , Y 1 and h depend only on η. The following expressions are assumed, where From Eq. (10), it reads that where the ( ) operator applied to any variable a (i.e., a ) means differentiation of a with respect to η or da/dη. Second and third derivatives follow the same notation. For the sake of simplicity, single dependence on η is implicit in the equations to follow.

Evaluation of the transverse velocity
The determination of the transverse velocity, v, across the mixing layer from the self-similar solution in variable-density mixing layers has been given little to no attention in previous works and textbooks.
Generally, the interest is only focused on solving the velocity profile in the main flow direction, together with the mixing profiles of other variables such as density and temperature. Some works addressing this issue are those of Kennedy and Gatski [32] and Libby and Liu [38]. However, relations for v are presented without further development. On the other hand, Pruett [39] offers some insights in the process to develop such relations. Nevertheless, the analysis performed on the boundary layer of an axisymmetric body complicates the final result.
In the present study, a straightforward development of an equation for v based on η and downstream location is presented in a form also shown in Sirignano [40]. The latter work develops a self-similar model for mixing in a variable-density laminar shear layer with imposed counterflow.
Eqs. (8) and (16) provide a differential relation between y and η for a fixed downstream position or constantx which can be integrated to relate the transverse location to the similarity variable as where η is a dummy variable for integration purposes andĨ(η) ≡ η 0 dη /ρ(η ). By combining expressions of the transformed transverse velocity, w, in (x, y) and (x, η) coordinates (Eqs. (9) and (19)), the following relation is obtained which shows that w (and thus v) are self-similar when multiplied by the square root ofx or downstream position. In the previous equation, the equality y 0 (∂ρ/∂x)dy = ∂( y 0 ρdy )/∂x = ∂z/∂x is used. Substitution of ∂z/∂x in Eq. (24) is found by differentiating dy/dx, dz/dx and dz/dx at constant y using Eq. (7).

Comments on the thermodynamic model
Another requirement needed to achieve similarity relates to the thermodynamic model used to evaluate fluid properties and transport coefficients. For these variables to depend only on η, the thermodynamic model must ultimately depend on pressure, temperature (i.e., enthalpy) and mixture composition. The self-similar model presented here relies on constant pressure everywhere in the domain, while temperature and mixture composition depend only on η.
The thermodynamic model implemented in this work is based on a volume-corrected Soave-Redlich-Kwong (SRK) cubic equation of state [41], which is able to represent non-ideal fluid states. The volumetric correction is needed to obtain accurate density predictions of high-density fluids (i.e., liquids) because the original SRK equation of state [42] presents liquid density errors of up to 20% when compared to experimental values [20,43]. The modified SRK equation of state is expressed in terms of the compressibility factor, Z, as where a(T ) is a temperature-dependent cohesive energy parameter, b represents a volumetric parameter related to the space occupied by the molecules and c(T ) is a temperature-dependent volume correction. R and R u are the specific gas constant and the universal gas constant, respectively. The solution of this cubic equation provides the density of the fluid mixture, ρ.
The volume-corrected SRK equation of state is used, together with high-pressure correlations, to evaluate fluid properties and transport coefficients [44][45][46] (e.g., viscosity). Details on the development and implementation of this thermodynamic model can be found in Davis et al. [29].

Boundary conditions and interface matching relations
The governing equations representing the non-ideal two-phase laminar mixing layer have been reduced to a system of ODEs, Eqs. (20)- (22), which depends only on the similarity variable, η. To solve this system of equations, proper boundary conditions must be imposed. These can be divided into freestream conditions for each individual fluid and interface matching relations between both fluid streams.
The self-similar transformation is defined such that the liquid-gas interface, which is denoted by Γ, does not move (V Γ = 0) and is located at y = z = η = 0. The justification for this approximation has been shown by Davis et al. [29]. Positive values of η define the gas stream, while negative values of η represent the liquid stream. Following this definition, freestream boundary conditions for f , Y and h become and At ±∞, both streams are pure components. Thus, if Y 1 is chosen to represent the pure gas species, Y G = 1 and Y L = 0 in Eqs. (32) and (33).
The streamwise velocity component and the tangential stress are continuous across the interface. These conditions are expressed in the original space (x, y) as which converts in the η space to Therefore, f is continuous across the interface while a jump in the second derivative, f , exists for Another condition for f is obtained from the mass balance across the interface. Since the interface location is assumed to remain fixed, the net mass flux,ω, across the interface readṡ sinceĨ = 0 at the interface. This equation shows that f is also continuous across the interface (i.e., f g (0) = f l (0) = f (0)). From this relation, the transformed velocity evaluated at the interface corresponds to the mass Finally, expressions for the species mass balance and energy balance at the interface are needed. Eqs. (40) and (41) show these relations in the (x, y) space, respectively.
Eqs. (40) and (41) Thermodynamic phase equilibrium is used to evaluate fully the interface solution of the system of ODEs.
That is, an equality in chemical potential for each species on either side of the interface is imposed. This condition is expressed through an equality in fugacity for each species [42,44], φ i , as which is a function of mixture temperature, pressure and composition. Under constant pressure across the interface, phase equilibrium can be expressed using the fugacity coefficient, where X i represents the mole fraction of species i.
Note that liquid and gas compositions are only the same at the mixture critical point. Therefore, phase equilibrium shows that Y and h are discontinuous across the interface, as well as other fluid properties (e.g., density).
Since the thickness of the interface is of the order of O(10 −9 m) [18,19] and diffusion layers in non-ideal high-pressure conditions quickly reach thicknesses of the order of O(10 −6 m) [17,29], the interface thickness is neglected and temperature is assumed to be continuous (i.e., T g (0) = T l (0) = T Γ ). The quick adjustment of the interface solution supports the thermodynamic equilibrium assumption [17,29].

Solution method
The system of ODEs, Eqs. (20)- (22), must be solved numerically. To avoid solving a third-order differential equation, Eq. (20) is split into two equations: a first-order differential equation to solve for f and a second-order differential equation to solve for f . Therefore, Eq. (20) is rewritten as where the new variables are g 1 = f and g 2 = f .
Eqs. (21), (22), (46) and (47) are discretized using a second-order finite difference scheme and solved using a tri-diagonal matrix solver. The asymptotic behaviors at +∞ and −∞ are well represented with a numerical range for the similarity variable given by −0.5 ≤ η ≤ 0.5. Since the system of equations is highly coupled, some iterations are needed before reaching a converged solution. A mesh size of ∆η = 1.5625×10 −4 kg/(m 5/2 ) provides a grid-independent solution.
The interface solution is also updated at every iteration. The streamwise momentum matching conditions, Eqs. (36) and (37), are combined numerically to solve for U Γ or f (0). On the other hand, Eqs. (42), (43) and (45) form a closed system that can be solved at every iteration in terms of the interface temperature. The solution of this system of equations provides the interface values for f (0), Y g (0), Y l (0), h g (0) and h l (0). More information on this interface algorithm is provided in Poblador-Ibanez and Sirignano [17].

Results
The system of ODEs has been solved for the analyzed cases from Davis et al. [29], where numerical results for a non-ideal two-phase laminar mixing layer using the system of partial differential equations, Eqs. (1), (2), (3) and (4), are provided.
Oxygen is chosen to be the pure gas species, while the liquid consists of pure n-decane. Mass fraction freestream conditions correspond to those discussed in Section 4, with Y 1 representing the pure gas species  Table 1.  In Section 6.1, the mixing layer solution from the system of ODEs is given for the four pressure cases from Table 1 and compared with the solution from the system of PDEs obtained in Davis et al. [29] for the same configurations. Then, in Section 6.2, the interface equilibrium solution is obtained with both approaches. Lastly, in Section 6.3, the evolution of the mass, momentum and thermal mixing layer thicknesses is determined and a generalization of their evaluation is made to avoid solving either the system of partial differential equations or the self-similar system.

Mixing layer profiles
For visualization purposes, the results are non-dimensionalized by the liquid freestream conditions to provide a consistent comparison between different cases. The non-dimensional similarity variable and f function are, respectively, The non-dimensional second derivative of f (or velocity gradient) follows that Any other variable, φ, is made dimensionless by dividing by the liquid freestream value, φ L (e.g., h * = h/h L or f * = f /u L ). Note that Y is already non-dimensional. Furthermore, for continuous variables across the interface, such as f and T , the following non-dimensional definitions are made to obtain distributions ranging from 0 to 1, being 1 the largest value of the variable, The transformed transverse velocity, w, is not self-similar since it depends onx and η (see Eq. (19)).
However, as stated in Section 3.1, w becomes a function of η only when multiplied by the square root ofx.
To obtain a non-dimensionalŵ for representation purposes, the results are scaled asŵ * =ŵ/ √ ρ L µ L u L . Figure 2 presents the non-dimensional solution of f * , f * and f * . Figure 2b shows how the interface velocity tends to be very close to the freestream liquid velocity. That is, the liquid phase is much more dense and viscous than the gas phase and it becomes hard for the slower gas stream to slow the liquid stream down.
The solutions of Y and h * are plotted in Figures 3 and 4.
As pressure increases well above the liquid critical pressure, the dissolution of the lighter gas species into the liquid phase is enhanced [16,17] (see Figure 3a). This generates sharper variations in the fluid properties within the liquid phase, which can be responsible for the change in behavior seen in the liquid velocity gradient (see Figure 5). As seen in Figure  bar. A discussion based on a thermodynamic analysis is provided in Poblador-Ibanez and Sirignano [17].
Considering a mass element containing the whole mixing layer (liquid and gas), the mixing process is just an internal phenomenon. Negligible net heat flux, Q, crosses the boundaries of the mass element. Therefore for a constant pressure process, ∆H = Q| p = 0. To compensate for the increase in enthalpy, H, in the liquid phase, the gas phase enthalpy must decrease accordingly. At 10 bar, this requirement cannot be solely  ity profiles is provided in Davis et al. [29]. Figure 6 shows the non-dimensional continuous distributions of streamwise velocity, θ u , and temperature, θ T . As pressure increases and the gas dissolves more easily into the liquid phase, the interface velocity and temperature deviate further from the freestream liquid values.
That is, as the liquid and gas phases come closer together, the liquid loses inertia and the gas stream can slow it down more. Furthermore, the enhanced dissolution of the gas into the liquid increases the heat flux by diffusion into the liquid, raising its temperature. Figure 7:ŵ * profiles from the solution of the self-similar system of ordinary differential equations.
Distributions of the non-dimensional modified transverse velocity,ŵ * , are plotted in Figure 7. Information on the phase change behavior at the interface can be extracted from these results. At the interface, w = ρv from Eq. (9). Therefore, w > 0 orŵ > 0 at the interface implies that the transverse velocity is also positive, v > 0. That is, net vaporization at the interface is occurring under the assumption that V Γ = 0. On the other hand, w < 0 orŵ < 0 implies v < 0 or net condensation. As inferred from Eq. (19), the strength of the mass flux across the interface decreases with downstream distance as the 1/ √x or 1/ √ x. At 10 and 50 bar, net vaporization is occurring, while condensation dominates at higher pressures. This interface behavior and its implications on the liquid jet breakup are further discussed in [17,21].
The agreement between the solution of the self-similar system of ODEs  at Re x ≈ 60 in the gas mixing layer and Re x ≈ 240 in the liquid mixing layer, as defined by Eq. (5).
This corresponds to a streamwise distance from the splitter plate of 0.06 mm and 0.24 mm, respectively.
Following this restriction, the analyzed locations are x = 1 mm, x = 5 mm and x = 10 mm.
The profiles shown in Figures 11-15 tend to collapse onto each other although some small deviations are seen depending on the variable and the fluid phase analyzed. The streamwise velocity profiles agree almost perfectly in both phases, but the mass fraction and density profiles deviate slightly. The greatest errors appear when comparing the temperature profile in the liquid phase, although the gas phase seems to match perfectly. The temperature range across the liquid thermal mixing layer is very small compared to the temperature difference between the two pure fluids, which can magnify errors on the liquid side plot while hiding them on the gas side plot.

Interface equilibrium solution
This section compares the interface equilibrium solution obtained with the self-similar model with the results shown in Davis et al. [29] obtained with the mixing layer equations. To assess the performance of the self-similar system of equations, the error of a given variable Θ between the solution of the ODE system and the PDE system is defined as representing a relative deviation error from the PDE solution. size. That is, diffusion is almost negligible in the liquid phase at very low pressures and it could affect the evaluation of enthalpy and concentration gradients needed to solve the interface matching conditions.
Larger errors are observed when comparing the mass flux across the interface at a given location (see Table 6). For 50, 100 and 150 bar, the errors are about 5%; but again the 10 bar case presents a larger error, in line with the comparison of other interface properties. Nevertheless, the transition from net vaporization to net condensation is still well captured.
Sensitivity issues might arise when evaluatingω. Even though the relative error in interface temperature, composition, etc. is well below 1%, the magnitude ofω itself is very small and easily perturbed by a small change in the other variables. This issue is then reflected in the transverse velocity profiles, as discussed in Section 6.1 and seen in Figures 10 and 15, and might be one of the main causes of the small mismatch between the ODE and the PDE results when comparing mixing profiles across the mixing layer.

Mixing layer thickness
Discussion in prior subsections addressed the similarity of solutions at differing downstream positions for a given problem configuration. The similarity between cases with differing constraints is discussed in this section.
The self-similar solution can provide information on the exact evolution of the mixing layer thickness.  Table 2: Comparison of the interface velocity, U Γ , and equilibrium temperature, T Γ , between the solutions obtained using the system of ordinary differential equations (ODE) and the system of partial differential equations (PDE) [29].  Table 3: Comparison of the interface equilibrium composition on each phase, Y g,O 2 and Y l,O 2 , between the solutions obtained using the system of ordinary differential equations (ODE) and the system of partial differential equations (PDE) [29].  Table 4: Comparison of the interface equilibrium mixture densities on each phase, ρ g and ρ l , between the solutions obtained using the system of ordinary differential equations (ODE) and the system of partial differential equations (PDE) [29]. This is true for simpler analyses, such as the self-similar solution of incompressible flow over a flat plate (a.k.a., the Blasius equation) [30]. However, semi-empirical formulas are usually provided to estimate boundary layer or mixing layer thicknesses for more complex flows, such as compressible flows, with different boundary conditions [30]. These equations rely on many flow parameters and are not as accurate as obtaining the mixing layer thickness directly from the similarity solution of the specific constraints.
A more accurate approach is presented in this paper, which covers variable density flows under non-ideal conditions. It relies on solving the self-similar system of equations coupled with a non-ideal thermodynamic  Table 5: Comparison of the interface equilibrium mixture enthalpies on each phase, h g and h l , between the solutions obtained using the system of ordinary differential equations (ODE) and the system of partial differential equations (PDE) [29].
p (bar)ω ODE (kg/ms 2 )ω PDE (kg/ms 2 ) Eω (%)  Table 6: Comparison of the net mass flux across the interface,ω, at the downstream location x =x = 0.01 m between the solutions obtained using the system of ordinary differential equations (ODE) and the system of partial differential equations (PDE) [29].
model. It is not as simple to implement as other semi-empirical models, but as discussed later, some generalization is possible. To evaluate the layer thickness on each side of the interface for mass, momentum and energy mixing regions, the following non-dimensional parameters are defined where all profiles vary from 0 to 1 going from the interface to the freestream condition. The mixture composition is used to analyze the mass mixing layer, while the streamwise velocity represents the momentum mixing layer. For the energy mixing layer, temperature has been used since it is a more representative variable than enthalpy. That is, the thermal mixing layer is represented, separately from composition influence on the enthalpy.

Case
Mixture p (bar) u G (m/s) u L (m/s) T G (K) T L (K)  Davis et al. [29] report the non-similarity across different problem constraints. However, when the proper fluid properties are used to non-dimensionalize the similarity variable and the mixing profiles, the thickness of the mixing layers almost collapse onto each other for different configurations. The analyzed cases are defined in Table 7. This is shown in Figures 16, 17 and 18. The profiles may also be nearly self-similar, although, in general, the variations in the behavior of the selected thermodynamic model creates different profiles as the problem configuration changes. The non-dimensional similarity variables used to plot the mass mixing layer are defined from freestream conditions for each phase as and the non-dimensional similarity variables for the momentum mixing layer and the thermal mixing layer are defined, respectively, as and The mixing layer thickness, δ, is obtained by using Eq. (23). δ = √ 2xĨ(η) just represents the transverse coordinate evaluated at the η value defining the limit of the mixing layer. This limit is obtained from A polynomial of degree 9 is fitted using a least-squares method among the data defining the distributions of mass fraction, velocity and temperature within the mixing layer for each case. That is, the following function is used to represent the diffusion profiles where p(η * ) represents θ * Y,l (η * Y,l ), θ * Y,g (η * Y,g ), θ * u,l (η * u,l ), θ * u,g (η * u,g ), θ * T,l (η * T,l ) or θ * T,g (η * T,g ). The polynomial reasonably predicts the diffusion profiles with the exception of the gas temperature profile and the liquid streamwise velocity profile, where larger variations are observed (see . The coefficients of Eq. (58) to predict each distribution are given in Table 8.  Table 8: Coefficients of the polynomial fit of θ * Y,l (η * Y,l ), θ * Y,g (η * Y,g ), θ * u,l (η * u,l ), θ * u,g (η * u,g ), θ * T,l (η * T,l ) and θ * T,g (η * T,g ).
A representative value of η * defining the mixing layer thickness for all scenarios is obtained from the polynomial fitting of each profile (i.e., where each polynomial profile reaches 99% of the freestream value).   values. This situation also happens for the liquid temperature profile. However, the velocity profiles in the liquid phase showed greater variations between each other. Nevertheless, the interface streamwise velocity tends to be very close to the liquid freestream velocity for all analyzed cases. Therefore, the magnitude errors between the exact solution and the polynomial estimate are very small even though θ * u,l profiles do not collapse onto each other.
The difference in the evaluation of the integral from Eq. (23) orĨ(η), i.e., whether the exact solution of the ODE system or the approximate method is used, causes different mixing layer thicknesses. Figure 19 shows the value of this integral where 99% of the freestream value is obtained using both approaches. Overall, estimated values are in good agreement with the ODE solution. The relative deviation error of the mixing layer thickness value from the exact ODE solution is shown in Figure 20. In general, errors are small and below 2%. However, some estimates present higher errors, but still below 5%.
Information on the mixing layer thickness can be directly obtained from Figure 19. As previously men-    tioned, the magnitude ofĨ(η) is used to estimate the mixing layer thickness at any downstream location as δ = √ 2xĨ(η). It is interesting to note that while the thermal mixing layer is always larger in the gas phase, the momentum mixing layer dominates by a larger factor in the liquid phase. In both phases, the mass mixing layer is the smallest.
The effects of increasing pressure on the mixing layer thickness are shown through cases A-D. The gas becomes more dense as pressure increases, reducing diffusivity. Thus, the gas mixing layer becomes thinner at larger pressures. On the other hand, the increased solubility of the gas into the liquid when pressure  increases widens the momentum and thermal mixing layers in the liquid phase. Figure 21 presents the evolution of the mass, momentum and thermal gas mixing layers for case D (i.e., a very high pressure case). The mixing layer thickness is computed using Eq. (23) and a comparison between the solution obtained with the exact ODE solution and the polynomial estimation is shown. Furthermore, data points from the PDE solution obtained in [29] are overlapped. The mixing layer thicknesses based on the PDE solution are obtained as the transverse location, y, where θ Y,g , θ u,g and θ T,g become 0.99. Figure 21: Mass, momentum and thermal gas mixing layer thickness evolution for case D (see Table 7). A comparison between the ODE solution, the polynomial approximation and the PDE solution is presented.
The evolution of the different mixing layer thicknesses is in accordance for the three methods used in this comparison. Errors between the ODE solution and the estimated solution based on the polynomial fit have already been discussed (see Figure 20). For both of them, the mixing layer must grow with downstream distance as √ x, which is corroborated with the PDE results from [29]. However, deviations between the self-similar model and the mixing layer equations are seen, especially for the evolution of the streamwise velocity and the temperature mixing layers, with errors of the same order as the relative deviations between the ODE solution and the approximate solution (< 5%).
Overall, the polynomial fitting approach is valid for the type of mixture and range of thermodynamic states evaluated in this work. It simplifies the evaluation of the mixing layer thickness with reasonable accuracy and only relies on the freestream conditions and the interface solution, which could be tabulated from previous simulations. Thus, it is not necessary to solve the system of ODEs or the system of PDEs for all cases and an estimate of the mixing layer thicknesses can be readily obtained. It is expected that other problem configurations not used in the development of the correlations could show larger errors than the ones reported for cases A-G.

Summary and conclusions
The non-ideal, two-phase, laminar mixing-layer equations have been reduced to a system of ordinary differential equations in terms of a similarity variable, η. The similarity transformation follows classical techniques used in compressible flows [31] and has been generalized to be implemented with any non-ideal thermodynamic model for the equation of state and phase equilibrium. In this work, the high-pressure thermodynamic model is based on a volume-corrected Soave-Redlich-Kwong cubic equation of state and other fundamental thermodynamic principles, coupled with high-pressure correlations to evaluate transport coefficients and phase equilibrium [29].
Good agreement between the self-similar solution and the solution of the system of partial differential equations [29] is obtained, proving the validity of the self-similar approach. This is an important step towards reducing the complexity of the analysis of supercritical two-phase flows, while still capturing the main physics involved (e.g., enhanced diffusion in the liquid phase, phase change reversal at the interface).
This type of approach can be helpful when used to implement realistic initial conditions to more complex flow simulations (i.e., high-pressure atomization). The results show the existence of two phases at pressures above the critical pressure of any chemical component. It is seen that condensation can occur even though the gas is hotter and heat is conducting into the liquid, as reported in previous works [17,21].
A correlation using a ninth order polynomial is used to fit data from various problem setups into a generalized approach to represent the mixing layer evolution. This reduces the set of ordinary differential equations to a set of functions which depend only on the freestream conditions, the interface equilibrium solution and a non-ideal equation of state. Good estimates of the evolution of the mixing layer thickness on both sides of the interface are obtained from this simplified model, with errors well below 5% in most cases.

Conflict of interest
The authors declared that there is no conflict of interest.