Evaporation Boundary Conditions for the Linear R13 Equations Based on the Onsager Theory

Due to the failure of the continuum hypothesis for higher Knudsen numbers, rarefied gases and microflows of gases are particularly difficult to model. Macroscopic transport equations compete with particle methods, such as the Direct Simulation Monte Carlo method (DSMC), to find accurate solutions in the rarefied gas regime. Due to growing interest in micro flow applications, such as micro fuel cells, it is important to model and understand evaporation in this flow regime. Here, evaporation boundary conditions for the R13 equations, which are macroscopic transport equations with applicability in the rarefied gas regime, are derived. The new equations utilize Onsager relations, linear relations between thermodynamic fluxes and forces, with constant coefficients, that need to be determined. For this, the boundary conditions are fitted to DSMC data and compared to other R13 boundary conditions from kinetic theory and Navier–Stokes–Fourier (NSF) solutions for two one-dimensional steady-state problems. Overall, the suggested fittings of the new phenomenological boundary conditions show better agreement with DSMC than the alternative kinetic theory evaporation boundary conditions for R13. Furthermore, the new evaporation boundary conditions for R13 are implemented in a code for the numerical solution of complex, two-dimensional geometries and compared to NSF solutions. Different flow patterns between R13 and NSF for higher Knudsen numbers are observed.


Introduction
For modelling ideal gas flow, there are in general two approaches, the microscopic and the macroscopic approach. In the microscopic approach, the Boltzmann equation [1,2] is solved, e.g., with the Direct Simulation Monte Carlo method (DSMC) [3]. However, tracking particles is computationally expensive, and for engineering applications, determining the macroscopic quantities is often sufficient. In the macroscopic approach, microscopic information is condensed into quantities such as mass density, bulk velocity, temperature, heat flux and stress. Macroscopic transport equations reduce the number of variables and when simplified allow for analytical solutions. The advantage of faster calculations is associated with the restriction to certain flow regimes. Flow regimes can be characterized by the Knudsen number, which is the ratio of the mean free path, i.e., the average distance a molecule travels between two subsequent collisions, and a characteristic length, e.g., the diameter of a pipe. For Knudsen numbers larger than Kn ≈ 4 × 10 −2 [4], the classical Navier-Stokes-Fourier (NSF) equations start to fail [4,5]. Applications for Knudsen numbers in

The R13 Equations
In the following, all equations are non-dimensionalized and linearized around an equilibrium state defined by a reference density for the vapour ρ 0 and reference temperature T 0 . The equilibrium saturation pressure for both liquid and vapour is defined as p 0 = p sat (T 0 ). We shall consider small deviations from equilibrium, caused by pressure or temperature gradients, to drive evaporation or condensation. Non-dimensionalizing allows one to introduce meaningful coefficients into the equations, e.g., Prandtl or Knudsen numbers. The connection between variables denoting non-dimensional deviation to an equilibrium state (with hat) and the regular variables with dimension is: Here, T is temperature, ρ mass density, p pressure, v k the velocity vector, q k the heat flux vector, σ ik the stress tensor, h enthalpy, u internal energy, η = ρs entropy density, x k the position vector and t time. From now on, the hats are not shown. The governing macroscopic equations that describe the gas are given by the conservation laws for mass, momentum and energy, which in linearized and dimensionless form read: ∂v i ∂t 3 2

∂T ∂t
Here, F i is a body force, e.g., gravitational force. One has five equations for the five unknowns ρ, v i and T. An algebraic equation for p is found in the ideal gas law p = ρRT, which assumes for the non-dimensional and linear case the form p = ρ + T, with all variables describing the deviation to the equilibrium state.
It is necessary to find equations for the heat flux vector q k and stress tensor σ ik , which beyond the hydrodynamic regime become full balance equations. By means of the order of magnitude method, Struchtrup and Torrilhon derived the following (here linearized and non-dimensionalized) balance equations from the Boltzmann equation, known as the regularized 13 moment equations [9], ∂σ ij ∂t ∂q i ∂t The higher moments are defined over the relations: By using the Chapman-Enskog expansion, while considering low Knudsen numbers, Equations (5) and (6) reduce to the laws of Navier-Stokes and Fourier, i.e., the left-hand sides become zero [5]. The balance laws (5) and (6) use the higher moments ∆, R ik and m ijk . Here, Pr = µc p k denotes the Prandtl number, with µ as the shear viscosity. For a monatomic gas, one has c p = 5 2 R as the isobaric specific heat and k = 15 4 µ as the thermal conductivity. The Knudsen number is Kn = µ √ RT pL , with L as the characteristic length, e.g., the diameter of a pipe. Here, θ 2 , θ 4 , w 2 and w 3 are coefficients for different collision models, such as Maxwell, HS and BGK models. In the following sections, only Maxwell molecules are considered; nevertheless, the corresponding coefficients for Maxwell, HS or BGK models for stress tensor, heat flux vector and higher moments can be found in Table 1 [12].

Macroscopic Evaporation Boundary Conditions for Maxwell Molecules
For the case that a vapour molecule hitting the liquid interface is reflected back to the vapour and not being absorbed, Maxwell proposed an accommodation model, which is based on the assumption that the fraction χ of the vapour molecules hitting the liquid surface are diffusively reflected, i.e., with momentum and energy exchange, and the remaining fraction (1 − χ) is specularly reflected, without energy exchange [7].
Based on microscopic evaporation boundary conditions of the Boltzmann equation, which are derived from a Maxwell model for the interface, Struchtrup et al. derived Macroscopic evaporation Boundary Conditions (MBC) for the R13 equations [11]. In these, interface effects are described through the accommodation coefficient χ and the evaporation coefficient ϑ. The evaporation coefficient equals the condensation coefficient, which is the probability that a vapour particle hitting the liquid interface will condense [16].
After non-dimensionalization and linearization around an equilibrium state, the MBC for evaporation [11] read: Here, the index n refers to the direction normal to the interface. The Einstein notation, i.e., A jj = 3 ∑ j=1 A jj , is not applicable for the index n. The variables are tensor components, where the overbar denotes the normal-tangential and the tilde the tangential-tangential parts; see Appendix A. Note that all variables describe the deviation to an equilibrium state.

Evaporation Boundary Conditions for Linear R13 Based on the Second Law of Thermodynamics
The MBC have the major drawback of stability problems; see [17]. Therefore, we aim to derive stable Phenomenological Boundary Conditions (PBC) for the regularized R13 equations for a liquid-gas interface. The approach follows [12], in which a reduced entropy balance is used to derive boundary conditions for a wall-gas interface. The entropy balance for a fluid with dimensionless entropy density η, entropy flux Ψ k and entropy generation rate Σ gen reads:

∂ η ∂t
Equation (16) shall be integrated over a small volume of area ∆A and height ∆z across the liquid-vapour interface. By using Gauss' theorem, the integrated entropy balance becomes: For ∆z → 0, the first term vanishes, and (17) reduces to the entropy balance for the interface, Hence, the entropy generation rate Σ sur f ace = 1 dA ∆A∆z Σ gen dV is equal to the difference in entropy fluxes entering and leaving the interface. In the following, all variables on the liquid side are denoted with l and all variables on the vapour side with g. A linear combination of manipulated mass, energy and entropy balances (Appendix B) leads to the (linearized and non-dimensional) entropy flux on the liquid side as: Here, T, ρ and v are deviations from an equilibrium state defined by T 0 , ρ 0 and p 0 = p sat (T 0 ). For the linear R13 equations and the vapour side, the linearized and dimensionless entropy flux (Appendix B) is: Furthermore, the (linearized and non-dimensional) balance laws for mass, momentum and energy, integrated around the interface similar to (18), become: p l n i + σ l ik n k = p g n i + σ g ik n k , The variables v l k and v g k are the velocities on the liquid and vapour sides from the perspective of an observer resting on the interface.
The entropy fluxes (19) and (20) are plugged into the integrated entropy balance (18). Equations (21)-(23) are used to eliminate the variables v l k , σ l ik and q l k . All variables describe the deviation to equilibrium, are dimensionless and linearized. After applying the appropriate coefficients for Maxwell molecules, according to Table 1, using the Clausius-Clapeyron equation [18] (linearized and dimensionless) in the form p sat T l = h 0 gl RT 0 T l and by considering ρ l ρ 0 , one may write (18) as: where k n k and the corresponding ideal gas law, given as ρ g = p g − T g , was used. To accomplish a proper entropy balance for the linearized equations, terms up to second order are kept [19].
Next, the entropy balance is split into contributions from normal and tangential components (see Appendix A); all matrices and higher moments are symmetric and trace free, As before, the overbar denotes normal-tangential, and the tilde denotes tangential-tangential components. In the case that the mass flow J g n vanishes, Equation (25) simplifies to the entropy generation at a wall-gas-interface; see [12].
The entropy generation may be written as a superposition of thermodynamic fluxes J i and forces X i [13,14]: Here, moments with odd degree in the normal direction n are identified as fluxes, i.e., J n , q n , m nnn , σ nk , R nk and m nij , while moments with even degree in n are identified as the corresponding forces, i.e., p g , T g , T l , σ nn , R nn , ∆, V k , q k , m nnk and σ ij . Note that p g , T g , T l , σ nn , R nn , ∆, J n , q n and m nnn are scalars, V k , q k , m nnk , σ nk and R nk are vectors and σ ij and m nij are tensors. Furthermore, a linear force-flux relation is stated within the Onsager theory, to satisfy Equation (26): Here, L ij is a positive-definite matrix of Onsager coefficients with the Onsager reciprocity relation, requiring symmetry of L ij . Only equations of the same tensor rank are coupled over the reciprocity relation (Curie principle [20]). This means that all force terms of the same tensor rank superimpose on each other and impact all fluxes of the same tensor rank; hence: Scalar fluxes: Vector fluxes: Tensor fluxes: For λ 0 = λ 1 = λ 2 = 0, one obtains the full set of phenomenological boundary conditions for a wall-gas interface, which are independent of evaporation as in [12]. The interface conditions (29) and (30), which consist of first order tensors (vectors) and second order tensors (matrices), respectively, have been fitted for a wall-gas interface in [12]. The fitting of (28) for evaporation at liquid-vapour interfaces shall be discussed in Section 3. In the following, the new evaporation boundary conditions (28)-(30) shall be referred to as PBC.

Comparison to Previous Macroscopic Boundary Conditions
The structure of PBC and MBC is very similar; the main difference lies in the values of the coefficients. As a first step for determining the Onsager coefficients of the PBC (28)-(30), we aim to use the coefficients of the MBC in a way that all terms, except those where higher order moments, i.e., ∆, R ij , m ijk , occur, agree with the MBC. This is justified due to the fact that the MBC predict effects in the Navier-Stokes regime very well. In the rarefied gas regime, however, their application seems to be more limited [11]. Since the higher moments are responsible for predicting a simplified Knudsen layer and also for rarefaction effects, a difference between PBC and MBC in these terms is desired. For a liquid-gas interface, the matrix of Onsager coefficients of those boundary conditions with variables of zero tensor rank (28) assumes the dimension 3 × 3, in contrast to the wall-gas interface, where the matrix reads 2 × 2 [12]. Based on these thoughts, the following Onsager coefficients are suggested: with: To leave the coefficients adjustable, the factors a... f have been introduced. For a = b = ... = f = 1, the PBC differ from the MBC, only in the higher order terms; see Appendix C. The boundary conditions (29) and (30) have been fitted for a wall-gas interface in [12] and shall not further be investigated here. To determine the coefficients a... f by fitting to a DSMC solution, two evaporation problems will be discussed, for which analytical solutions for R13 with PBC can be obtained.

Simplification of R13 for 1D Problems
As can be expected, the present PBC, just like the MBC, give less accurate results than methods that solve the full Boltzmann equation. The R13 equations and their corresponding interface and boundary conditions are simplifications of the Boltzmann equation and carry less information. The adjustable coefficients a... f in (31)-(36) leave six degrees of freedom to determine the Onsager coefficients. It is of interest whether the simplification of R13 to the Boltzmann equation can be partly corrected by adjusting the Onsager coefficients. In this context, we simplify the linear R13 equations for one-dimensional and steady systems and solve them for two problems, previously discussed in [11]. Then, the new solutions are fitted to DSMC data.
All variables depend only on the location x. For the equilibrium rest state, the saturation pressure of the liquid interface is set to p sat (T 0 ) = p 0 . We assume that the liquid temperature at the interface is controlled. Small pressure or temperature changes are sufficient to drive evaporation or condensation. All equations are linear and dimensionless and describe the deviation to their equilibrium state. The simplified balance equations for mass, momentum and energy read: After, simple integration follows: Hence, velocity and conductive heat flux are constant in the vapour phase. The normal components of the linear and non-dimensional constitutive equations for (7)-(9) obtain the form: with data to adjust between the molecule models from Table 1. The linear and non-dimensional equations for normal stress σ and conductive heat flux q o become: Integration yields: with A, B, K as constants of integration. There are six unknowns (V 0 , P 0 , Q 0 , A, B, K), that must be determined for finding the solution. For evaporating interfaces and by taking ∆ = R = 0 (39) into account, the normal boundary conditions (28) simplify to: with V o = n k V k and q o = q k n k .

Problem I: Vapour Layer between Two Liquid Reservoirs
In the first problem for fitting the coefficients a... f and also for gaining insight into the Knudsen layers, we consider one-dimensional, steady-state heat and mass transfer within a vapour phase in between two liquid reservoirs with controlled temperature on the liquid side of the liquid-vapour interfaces. The configuration shown in Figure 1 has been discussed in [11] and shall be outlined only briefly here. The interfaces are located at x = ± 1 2 with the normal vector n pointing from liquid into vapour and the superscripts 0 for x = − 1 2 and 1 for x = 1 2 , i.e., V 0 0 = −V 1 0 = V 0 . The driving force for evaporation and condensation is the temperature difference between T 0 l and T 1 l . The required six equations are found by evaluating the boundary conditions (28) at both interfaces. For evaluation of the equations, it is convenient to take both the sums and the differences at both interfaces. For the three sums, it follows: Stress profile Equation (42) and temperature profile Equation (43) follow as: The three differences of the normal boundary conditions form a linear system for V 0 , Q 0 and A as: 1 Kn Here, A is the amplitude of the Knudsen layer. We refrain from showing the solution, and will only show results from the inversion in the figures. For the linear NSF-Onsager boundary conditions (see Appendix D), one finds: The given solution for NSF is a simplification for χ = ϑ = 1; see Appendix D. For the NSF-Onsager coefficients r 11 , r 12 and r 22 , the Onsager matrix (A30) or the corrected Onsager matrix (A31) can be used. The solution of the MBC for this system can be found in [11]. Results shall be compared in Sections 3.5 and 3.6.

Problem II: Evaporation in the Half-Space Problem
In the half space problem, a liquid interface evaporates into the equilibrium state, as discussed previously in [11]. The driving force is the prescribed pressure p ∞ far away from the interface; see Figure 2. The six unknowns are found by considering evaporation boundary conditions on one side and constant velocity v ∞ = V 0 , pressure p ∞ = P 0 and temperature T ∞ far away from the interface. For reaching constant pressure p ∞ and due to the momentum balance (38), it is necessary to set the normal stress far away from the interface to σ ∞ = 0. Moreover, conductive heat flux q 0 is set to zero, as well. With T ∞ prescribed, one finds the constant K. For (50) and (51), it follows: Evaluating the boundary conditions (28) at the interface between liquid and vapour leads to: For the Navier-Stokes-Fourier equation out of Equation (A29), it follows: With prescribed pressure p ∞ and by setting p sat (T l ) − p ∞ = ∆p and T l − T ∞ = ∆T, there are three unknowns v ∞ , T ∞ and A, which can be calculated with (59)-(61) for PBC and (62) and (63) for NSF. The solution for the MBC can again be found in [11]. Note that for NSF, A is zero, and the two given equations are sufficient.
Ytrehus, who discussed the half space problem in [15], proposed dimensionless ratios in which the prescribed pressure p ∞ is eliminated. The ratios that make it easy to compare different models, e.g., Maxwell molecules, BGK, Navier-Stokes-Fourier, etc., read: Note that (59)-(63) and therefore also (64) and (65) are independent of the Knudsen number.

Fitting of the Onsager Coefficients: Standard Temperature Profile
The ratios (64) and (65) from Problem II together with DSMC data for Problem I shall be used to fit the coefficients a... f in (31)-(36). The temperatures and saturation pressures at the liquid boundaries are given as T 0 l = p sat (T 0 l ) = 1.05 and T 1 l = p sat (T 1 l ) = 0.95. All results in the following are based on full evaporation and fully-diffusive reflection, by setting the evaporation and accommodation coefficients ϑ = χ = 1. Maxwell molecules are considered, and their data are taken out of Table 1. In Table 2, factors for the Onsager coefficients, used in Equations (31)-(36), which have been found by trial and error, are suggested to adjust the PBC, Equations (28), for the best fit. The results of the new PBC are compared with the previously derived evaporation boundary conditions (MBC) and also with Navier-Stokes-Fourier solutions. NSF is based on Onsager boundary conditions, as well, and uses the Onsager matrix (A30) or the corrected Onsager matrix (A31). Ytrehus used a moment method to solve the half space problem with high precision [15] and his results are used here as a reference. Ytrehus' ratios α p , α θ (64) and (65) have been calculated for PBC, MBC, NSF and corrected NSF. Together with the percentual deviation to Ytrehus' solution, they are given in Table 3. By trial and error fitting of the Onsager coefficients, it was not possible to achieve superior agreement between PBC and DSMC for Problem I (Section 3.3) and proper results for Ytrehus' ratios (64) and (65) at the same time. Forcing good agreement between Ytrehus' solution of the half space problem and PBC regarding the dimensionless ratios showed a significant decrease in agreement between PBC and DSMC for Problem I. The fittings that are chosen here are compromises between Problem I and Problem II, but with strong emphasis on achieving proper results for Problem I, which means proper agreement with DSMC results. Figure 3 shows temperature and normal stress profiles for Kn = 0.078. R13 with PBC (solid, purple) and MBC (solid, red) are in good agreement with DSMC (green, dashed). The amplitude of the Knudsen layer A is zero for NSF (black, dashed) and corrected NSF (blue, dashed). As a result, both NSF solutions slightly deviate from DSMC close to the boundaries. A = 0 removes the last term in (51) and therefore leads to a linear function. In Problem I, NSF is not able to predict normal stress at all; see Equations (55) and (56).
In Figure 4, temperature and normal stress profiles are illustrated for Kn = 0.235. Both sets of boundary conditions for R13 reconstruct the DSMC results well, but slightly underpredict the Knudsen layers both for temperature and normal stress. For the temperature profile, they are in better agreement with DSMC than the two NSF solutions. For both Kn = 0.078 and Kn = 0.235, one notes the significant temperature jumps at the boundaries.
In addition to temperature and normal stress profiles, we seek to gain insight into the three integration constants velocity V 0 , heat conduction q 0 and Knudsen Layer amplitude A, depending on the Knudsen number. The three variables are plotted over Kn = {0, ..., 1} in Figure 5.
The signs of velocity V 0 and heat conduction q 0 are positive. That is, mass and conductive heat flux are transferred from warm to cold, which means they are transported at x = − 1 2 into the system via evaporation, and due to the steady state, the same amount of mass and conductive heat is transported at x = 1 2 out of the system into the colder reservoir via condensation.   The purple, large, dashed line represents R13 with PBC for a = b... = f = 1; see Appendix C. Although there are differences in the higher order terms between PBC and MBC, if the adjustable coefficients are set to unity, the order of magnitude of the maximum deviation between the two models is with ±10 −7 very small, i.e., at first glance, both plots appear to be identical.
R13 with PBC shows very good agreement with DSMC for V 0 and q 0 for all Knudsen numbers. The PBC results for normal stress are better than those of MBC for Kn < 0.3. For higher Knudsen numbers, both PBC and MBC fail to predict σ in precise agreement with DSMC. Again, normal stress cannot be predicted by NSF.
Interestingly, for this PBC fit, Ytrehus' ratios are similar to those of the MBC, i.e., 1.4% (PBC) and 0.74% (MBC) deviation for α p and 10.02% (PBC) and 10.44% (MBC) for α θ ; see Table 3. Corrected NSF is under 1% deviation for both ratios. Uncorrected NSF shows zero deviation for α θ and 6.18% for α p . For Knudsen numbers larger than Kn = 0.235, the deviation between DSMC and PBC becomes slightly larger for the temperature profile and stays similar for the normal stress profile. The temperature jump at the boundaries increases with increasing Knudsen number.

Fitting of the Onsager Coefficients: Inverted Temperature Profile
By adjusting the values for ∆T and ∆p, it can be shown that the sign of the conductive heat flux q 0 switches. This leads to an inverted temperature profile as depicted below. The negative sign of q 0 indicates conductive heat transport from x = 1 2 to x = − 1 2 ; see Figure 1. However, the second law is not violated, since the overall heat transport is given with Q = ρV 0 h + q 0 , and the advective term ρV 0 h is dominant. Hence, the overall heat Q is transported from hot to cold as expected. One notes that due to the reversed sign of the conductive heat flux, the necessary vapourization enthalpy is partly provided by the colder boundary. The liquid temperatures at the boundaries are set to T 0 l = 1.01 and T 1 l = 0.99 and the respective saturation pressures to p sat (T 0 l ) = 1.0752 and p sat (T 1 l ) = 0.9248. Therefore, the evaporating material of the system is different from the one considered for the standard temperature profile. The small temperature difference between hot and cold boundaries and the large difference between the saturation pressures allow for a temperature jump large enough to reverse the sign of the conductive heat flux.
By fitting with trial and error, it was not possible to achieve good fits for the standard and inverted temperature profiles at the same time. We believe that this is due to the evaporating material being different between the standard and inverted cases, since the saturation pressures are different. Therefore, we present a fitting for the adjustable factors within the PBC for the inverted case, which is given in Table 4. The ratios α p ,α θ , as well as the percentual deviation to Ytrehus' solution are presented in Table 5. The temperature and stress profiles for Kn = 0.078 are given in Figure 6. As a comparison to the new fitting, a PBC solution, which uses the previous coefficients, is given, as well (purple, dashed). R13 with PBC and MBC both overpredict the Knudsen layer at the interfaces. This inaccuracy of Knudsen layer modelling is due to the small number of moments, used in the R13 equations; see [21]. For the temperature profile, corrected NSF shows the best agreement with DSMC here. Normal stress is predicted well for PBC and MBC and is again zero for NSF.
For Kn = 0.235, the overprediction of the R13 boundary conditions becomes so large that the profiles are no longer inverted, as shown in Figure 7. Note that it is possible to "turn" the PBC temperature profile to match the DSMC results; however, this leads to worse results for other plots.
In this case, MBC shows slightly better results for temperature and normal stress profiles than PBC.   For evaporation velocity V 0 and conductive heat flux q 0 , R13 with PBC is in very good agreement with DSMC. In comparison to the standard temperature profile, the normal boundary stress of the PBC starts to differ from DSMC earlier, i.e., for Kn > 0.1. Corrected NSF is in surprisingly good agreement with DSMC for Kn < 0.3, but fails to predict normal boundary stress. Except for temperature and normal stress profiles for Kn = 0.235, R13 with PBC shows the best agreement with DSMC compared to all discussed models here. One notes that for this PBC fitting, the deviations of 5.11% in α θ and 0.46% in α p to Ytrehus' solution become smaller than for the standard profile.

Impact of Evaporation and Accommodation Coefficients
To gain a better understanding of the impact of evaporation and accommodation coefficients, the PBC shall be tested for the standard temperature profile of the previously discussed problem and a variety of ϑ, χ. Figure 9 illustrates solutions of the PBC for Problem I (Section 3.3) together with the fitting from Table 2 and Kn = 0.078. The plots are based on χ = 0.1 (green), χ = 0.5 (red), χ = 1 (blue), ϑ = 0.1 (solid), ϑ = 0.5 (dashed) and ϑ = 1 (large dashed). For ϑ = 1, the solutions are independent of χ. Since the evaporation coefficient is defined through the condensation coefficient, this may be explained due to the fact that for the condensation coefficient being unity, no reflection occurs, and all vapour molecules hitting the liquid interface are condensed. The largest temperature jump between gas and the boundary is found for ϑ = 0.1 and χ = 0.1 and the smallest for χ = 1.
The stress profile seems to be dependent mainly on the evaporation coefficient. The accommodation coefficient has only a small impact for ϑ = 0.5. The largest stress can be found for ϑ = 1. Evaporation velocity V 0 , conductive heat flux q 0 and boundary normal stress σ for various values of ϑ and χ are depicted in Figure 10.

Notes on the Meaning of the Individual Onsager Coefficients of the Normal Fluxes
The fittings used in the Tables 2 and 4 are based on a trial and error procedure, in which the factors a... f within the Onsager coefficients (31)-(36) are individually adjusted. Due to symmetry of the Onsager matrix, six independent parameters need to be determined. The tuning of the Onsager coefficients one by one gives an insight into their respective impact. However, one notes that due to the coupling within the Onsager matrix in Equation (28), the individual Onsager coefficient impacts multiple fluxes. The following is an attempt to highlight some trends, which were observed during the fitting procedure.
Since λ 0 appears only in the equation for the normal velocity, it has a strong impact on V 0 and no impact on the conductive heat flux q 0 . Apparently, it has no impact on the boundary normal stress σ. Temperature and stress profiles appear to be independent of λ 0 as well. The coefficient λ 1 has a big impact on V 0 and q 0 and a small impact on σ. It has a major impact on the temperature profile and a smaller impact on the stress profile. λ 2 strongly influences V 0 and σ and very slightly q 0 . Since λ 2 does not appear in the equation for q 0 , this is expected. It has an impact on temperature and stress profiles, but with clear emphasis on the stress profile. The coefficient λ 3 seems to play a key role in the fitting. Even though it appears only in the equation for q 0 , it has not only a strong impact on the magnitude and slope of q 0 , but also on those of V 0 and σ. Regarding the profiles, λ 3 seems to impact mainly the temperature and only very slightly the stress. The Onsager coefficient λ 4 mainly impacts σ, but also V 0 , q 0 and both profiles, with stronger impact on the stress profile, as expected. λ 5 appears only in the equation for the normal component of the higher moment m nnn . The coefficient has a strong impact on σ, a medium impact on V 0 and no impact on q 0 . It influences the stress profile significantly and the temperature profile slightly.
After these dependencies were established, several rounds of fitting were done, until a reasonable fitting was obtained.

R13 with Onsager Boundary Conditions in Numerical Simulation
It shall be shown that the applicability of R13 with PBC (Phenomenological Boundary Conditions) is not limited to one-dimensional systems. The code of Torrilhon and Sarna [22], written in C++, is used in this section to solve the R13 equations with PBC for evaporation. As a comparison, simplified NSF (Navier-Stokes-Fourier) is solved with the same program. Torrilhon and Sarna's code allows for generic implementation of macroscopic transport equations. The numerical solver relies on a Discontinuous Galerkin (DG) method, which utilizes finite elements to discretize the system. Here, the code is extended by implementing the evaporation boundary conditions previously derived in Section 3 and also simplified Onsager boundary conditions for NSF.
The PBC for R13, given in Equations (28)-(30), are adjusted by using data for Maxwell molecules out of Table 1. The liquid phase is not solved and therefore can be treated in the same manner as a wall, which allows for mass transfer. Adjustment of the Onsager coefficients allows one to derive other boundary conditions, such as the wall with energy transfer or inflow/outflow. Table 6 gives an overview of these modifications. For an adiabatic wall (fully specular reflective), all Onsager coefficients are set to zero, which leads to v g n = q g n = m nnn = σ g nk = R nk = m nij = 0. The Onsager coefficients for a wall with energy transfer are taken from [12]. The adjustable coefficients within the Onsager coefficients for the different boundaries were already implemented in Table 6.
Note: Compared to Section 3.1, a slightly different fitting is used here. Additionally, the coefficients used in λ 0 , ..., λ 5 are based on adjustments as in Problem I (Section 3.3); however, different definitions of the Knudsen number between DSMC and R13 were used. Therefore, a small error is introduced here.
The coefficients in ζ 0 , ..., ζ 2 and κ 0 are not fitted and set to unity. The adjustable coefficients for a wall with energy transfer λ 3 , ..., λ 5 and ζ 0 , ..., ζ 2 are taken from [12], and κ 0 is set to unity here. Depending on the boundary, different pressures and temperatures are assumed, as depicted in Table 7. Table 7. Overview of input parameters for the boundary conditions.

E Vapouration/Condensation W All with Energy Transfer I Inflow/Outflow
For a detailed description of the numerical solution, see [22].

Navier-Stokes-Fourier with Onsager Boundary Conditions in Numerical Simulation
For obtaining a comparison to the R13 solutions for two-dimensional systems, the Navier-Stokes-Fourier equations together with Onsager boundary conditions for evaporation/condensation are used here. For χ = ϑ = 1 and considering one-dimensional geometry, evaporation boundary conditions for NSF are given in Appendix D; see (A29). For two-and three-dimensional geometries, an additional boundary condition is found in [11] and reads: Note that Equations (A29) are simplified equations for 1D geometry. Again, by considering χ = ϑ = 1 and after full linearization and non-dimensionalization, Equation (66) becomes:

Numerical Solutions for Two-Dimensional Channel-Flow with four Evaporating Cylinders
The system of interest for the two-dimensional, steady-state simulation is a channel with four evaporating cylinders, which is discretized as depicted in Figure 11. The left boundary is the inlet of the channel flow, and the right boundary is the outlet. The top and bottom are walls, which allow energy transfer. The cylinder walls use evaporation boundary conditions given by (28)-(30) with Table 6 for R13 and (67), (A29) and (A31) for NSF.
The input parameters, which are given in Table 8, are non-dimensional and describe the deviation to equilibrium. They are chosen in a way that evaporation at the cylinders can be observed clearly. Table 8. Input parameters for two-dimensional channel flow with four evaporating cylinders.

E Vapouration/Condensation W All with Energy Transfer I Inflow/Outflow
The plots in Figure 12 show pressure contours, superimposed by velocity streamlines, for R13 and NSF, for the three Knudsen numbers: Kn = {0.1, 0.5, 1}. For Kn = 0.1, the velocity streamlines are similar between R13 and NSF. The inflow of the left boundary collides with the evaporating flow, which leaves the two cylinders on the left-hand side. The largest flow velocity is observed in between the two cylinders on the right-hand side.
For Kn = 0.5, the evaporation overcomes the inflow and leaves the system at the inlet of the channel. This interesting effect is observed for R13 and NSF, but with different flow behaviour. For R13, the streamlines, which leave the inlet, have their origin mainly in the left bottom cylinder. The dominance of the left cylinder of R13 becomes even more apparent for Kn = 1. The NSF velocity streamlines at the inlet for Kn = {0.5, 1} come almost equally from both cylinders on the left-hand side.
For Kn = 0.1, the pressure contours of R13 and NSF show very similar behaviour. With increasing Kn, the R13-pressure contours on the right-hand side of the diagrams disconnect from each other and become almost vertical for Kn = 1.
Furthermore, for Kn = 1, significant differences between R13 and NSF are found for the temperature profiles, which are depicted in Figure 13.
The overall temperature around the four evaporating cylinders is much lower for NSF than for R13. As can be seen by the conductive heat flux streamlines, the enthalpy of vapourization is provided by the boundaries, as in the previous simulations. The magnitude of the R13 heat flux shows interesting peaks in between the two cylinders on the right-hand side for Kn = {0.5, 1}. The large differences between R13 and NSF for Kn = {0.5, 1} are likely due to rarefaction effects, which cannot be captured by NSF. It has to be taken into account, as mentioned in Section 4.2, that simplified NSF boundary conditions are used here. Note that R13 is limited to flow regimes below Kn = 1 and can only describe a tendency here. For validation of the R13 results, a reliable reference, such as from a DSMC simulation, is necessary, which might be part of future work.

Conclusions
Based on the Onsager theory, which utilizes the second law of thermodynamics, evaporation boundary conditions (PBC) for the R13 equations are derived. The Onsager coefficients have been determined by following a process consisting of three steps: In the first step (Section 3.1), the boundary conditions are compared with previously discussed boundary conditions for evaporation (MBC), which represent an alternative approach for deriving boundary conditions for R13. Under the assumption of proper results for MBC in the Navier-Stokes-Fourier (NSF) regime and by keeping in mind that higher moments develop a significant impact only for higher Knudsen numbers, coefficients are taken over from MBC to PBC so that the differences between the sets of boundary conditions lie only in the terms with higher moments [12]. The idea is to find boundary conditions that are just as reliable as MBC in the NSF regime and more accurate in the rarefied gas regime. In the next step, adjustable coefficients are suggested for the PBC. These coefficients are fitted by trial and error to DSMC data for the analytical solution of a finite, one-dimensional system (Section 3.3). In the third step for finding meaningful Onsager coefficients, the half space problem (Section 3.4) is solved analytically, and ratios suggested by Ytrehus [15] are used to fine-tune the coefficients. The overall agreement between PBC and DSMC (Section 3.5 and 3.6) has been shown to be better than for MBC/NSF and DSMC. Even though there are differences in the higher order terms, when setting the adjustable coefficients a... f of the PBC to unity, the maximum deviation to the MBC, for the boundary values of the finite problem, is in the order of magnitude of ±10 −7 , only.
For a general approach to convert MBC to PBC, with differences in the higher order terms only, see [17]. Kinetic boundary conditions, such are used in [6,11,21,23,24], might lead to violation of the Onsager symmetry relations. Furthermore, due to the approximative nature of the models, there can be small inaccuracies in the results, e.g., due to the details of the Knudsen layers that cannot be fully described [21]. The present approach uses fitting of coefficients to recover Onsager symmetry and also to improve the accuracy of the results by small adjustments of the kinetic coefficients.
The impact of the evaporation and accommodation coefficients is discussed in Section 3.7. In Section 3.8, it is explained how the trial and error fitting gives an insight into the meaning of the individual Onsager coefficients.
Due to lack of a mathematical approach for the fitting, i.e., an optimization algorithm, it is uncertain if significantly better fittings for the presented problems are possible. This may be part of a future analysis. Even though NSF fails to predict normal stress for the presented systems, it shows surprisingly good results for low to moderate Knudsen numbers. The advantage of R13 with PBC compared to NSF might be shown even more clearly in numerical simulations for complex geometries. The Onsager coefficients appear to be dependent on the evaporating material, which in the practical application becomes problematic. Therefore, we recommend an investigation considering the fitting of Onsager coefficients as a function of the enthalpy of vapourization, which defines the material.
In Section 4, the new evaporation/condensation boundary conditions are implemented in a code for the numerical solution of two-dimensional, steady-state problems. Results for Knudsen numbers of Kn = {0.1, 0.5, 1.0} are obtained and compared to simplified Navier-Stokes-Fourier solutions. It is observed that with increasing Knudsen number, R13 shows different flow behaviour than NSF.
It is necessary to compare these results to a reliable reference, such as a DSMC solution, which shall be a future effort. Additionally, it might be of interest to compare the numerical R13 results to those of a 26-moment method; see [25]. . (A30) The solutions based on (A30) are referred to as uncorrected NSF. A correction can be found in kinetic theory, which yields [11,26]: