Numerical simulation of two-phase flows with heat and mass transfer

We present a finite element method for simulating complex free surface flow. 
 The mathematical model and the numerical method take into account two-phase non-isothermal 
 flow of an incompressible liquid and a gas phase, capillary forces at the interface of both 
 fluids, Marangoni effects due to temperature variation of the interface and mass transport 
 across the interface by evaporation/condensation. The method is applied to two examples from 
 microgravity research, for which experimental data are available.

1. Introduction. Mass transport across an interface liquid/gas is a fundamental process in many technical applications. Coupled to fluid dynamics in either phase and then in turn also to an evolution of the interface leads to a very complex system. Thus numerical simulations play an important role in understanding the relative influence of the underlying different physical mechanisms as well as predicting the behavior of the system over time.
The application motivating the present study originates from space technology, more precisely from basic questions of propellants and their management during different flight phases (propelled/ballistic) of launchers.
In particular, knowing the location of the free surface liquid/gas and the heat transfer from the wall towards the fluid is of paramount importance in the design and the optimization of cryogenic upper stage tanks for launchers with ballistic phases, where residual accelerations are smaller by up to four orders of magnitude compared to the gravity acceleration on earth. This changes the driving forces drastically: free surfaces become capillary dominated and natural or free convection is replaced by thermocapillary convection if a non-condensable gas is present.
The corresponding thermocapillary convection (driven by the temperature dependence of the surface tension) creates a velocity field directed away from the hot wall towards the colder liquid and then again back towards the wall. A deformation of the free surface resulting in an apparent contact angle is observed. The thermocapillary flow convects the heat from the wall to the liquid and significantly increases the heat transfer compared to pure conduction.
In this article we present finite element based numerical methods capable of simulating such phenomena for model geometries and scenarios. The present article is a report about the numerical methods developed during a long term collaboration with the Fluid Dynamics Group at ZARM, Bremen (Germany), headed by M. Dreyer. We gratefully acknowledge this stimulating collaboration, providing excellent experimental devices and data.
Our mathematical model and the numerical method take into account two-phase non-isothermal flow of an incompressible liquid and a gas phase, capillary forces at the interface of both fluids, Marangoni effects due to temperature variation of the interface and mass transport across the interface by evaporation/condensation.
Most numerical methods for moving boundary problems can be categorized into two fundamentally distinct classes: interface tracking and interface capturing methods. In interface capturing methods, the interface is described implicitly by an additional scalar function. Prominent representatives of this method are VOF [20,42], phase field approaches [10,24] and level set methods [18,37,45,47]. The latter is the by far most common approach in two phase flow. One of the reasons is that it computes the signed distance function to the interface, which has several appealing features.
Extensions of these methods in order to incorporate evaporation/condensation are treated for instance in [16,46,48,53] for the level-set method, in [43,51] for VOF and for the phase-field method in [25,39].
In contrast to interface capturing methods, in interface tracking the phase boundary is described explicitly in terms of the computational mesh. Consequently, the computational mesh deforms according to the interface motion, see for example [2,7,23]. While usually limited to situations in which only moderate deformations and no topological changes of the interface occur, the explicit representation of the interface allows for a very accurate treatment of surface tension. Interface tracking is thus the method of choice in the present situation, where in the range of the considered experiments the topology of the interface stays fixed. Our method further benefits from interface tracking by a variational formulation of the curvature terms, allowing also for a semi-implicit treatment of the latter, see Section 3 below.
The rest of the article is organized as follows. In Section 2 we sketch the derivation of our mathematical model and then formulate it in dimensionless form. Section 3 gives an overview of the finite element method used to simulate the system. At the heart of the method is a mesh moving/interface tracking method and a variational formulation of the capillary free surface condition. In Section 4 we demonstrate the performance of our method for two experiments that were carried out under microgravity conditions, one in the ballistic phase of a sounding rocket flight, the second one in the drop tower facilities at ZARM, Bremen. The article is concluded by some remarks in Section 5.

2.
Mathematical model of two-phase flows with heat and mass transfer. In this chapter, we briefly revisit the mathematical model to describe the system under consideration. The model is derived from basic balance laws with appropriate simplifications due to the specific application. The modeling is rather classical, see for instance [44]. We follow the presentation in [28], where a more detailed derivation can be found.
2.1. Setting and generic balance equations. In view of the application we have in mind a sketch of the typical geometry is depicted in Fig. 1. The overall domain Ω is divided into a domain Ω l occupied by a liquid phase, a domain Ω g occupied by a gaseous phase Ω g and (optionally) a domain Ω s consisting of a solid wall. The liquid Figure 1. Sketch of the geometry and gaseous phases are assumed to be separated by a smooth sharp interface Γ S which may move over time. n denotes the normal vector on the interface Γ S pointing from Ω l into Ω g . The phase Ω l is considered to consist of a pure substance, while the gaseous phase Ω g may consist of a mixture of non-condensable gas and the liquid's vapor. In what follows, all material parameters are assumed to be constant in each phase for simplicity. Our model is based on the balance equations for mass, momentum and energy. Consider the density of some generic extensive thermodynamical quantity a. Let V ⊂ Ω be a test volume. The rate of change of ρa in V is determined by contributions from fluxes across the boundary and sources and sinks in V : Here, ρ denotes the density of the fluid, u the fluid velocity, q d denotes the diffusive flux of a, f denotes the sources and sinks and n S the outward pointing normal of S := ∂V . Assuming the test volume V to be completely contained in one of the phases Ω g , Ω l , applying Gauss' theorem to Eq. (1) yields Since Eq. (2) holds for arbitrary test volumes V it follows that We now set a = 1, q d = 0 and f = 0 in Eq. (3) to obtain the equation for the balance of mass: (4) Combining the last two equations finally yields the generic balance equation In a similar way, but working with a test volume transported by the flow field, one derives the balance equations for the generic quantity a on the interface Γ S (we refer to [28] for a detailed derivation): Here, [g] := lim ε→0 + g(· + εn) − lim ε→0 − g(· + εn) denotes the jump of a generic function g across Γ S and u Γ denotes the normal velocity of the interface.

2.2.
Balance equations for mass, momentum and energy. Starting from the generic balance equations Eq. (5) and Eq. (6), the corresponding equations for the balance of mass, momentum and energy can be derived.
2.2.1. Conservation of mass for individual species. The equation for the conservation of mass was already derived for the overall mass density of the fluid ρ, see Eq. (4). Furthermore, we assume the pure liquid phase to be incompressible, i.e. the density to be constant in time and space. Eq. (4) then simplifies to In contrast to the liquid, as already mentioned, the gaseous phase may consist of a mixture of non-condensable gas and vapor. Therefore, the conservation equation for mass has to be formulated with respect to the individual species. To this end denote by ω the mass fraction of vapor to non-condensable gas. The diffusive flux for the species in Eq.
2.2.2. Conservation of momentum. Inserting a = u in Eq. (5) yields the equation for the conservation of momentum. According to Cauchy's theorem [19], stresses acting on the surface of a test volume are given by Tn S , where T is the stress tensor and n S the outward pointing normal of S. Therefore, −T defines the diffusive flux of momentum in Eq. (5). The stress tensor T is given by where p denotes the pressure and τ the viscous stresses, which is, assuming the liquid to be Newtonian, given by with µ the dynamic viscosity, η the bulk viscosity and D(u) = ∇u + (∇u) T the rate of strain tensor. In this work, we consider external forces (gravity) of the type with a function g that may depend on time but is constant in space. The conservation of momentum therefore reads 2.2.3. Conservation of energy. The energy density in the system is given by the sum of internal energy density ρe u and kinetic energy density ρ u·u 2 . The stress exerted by the fluid on S is given by Tu. Therefore, −Tu defines the diffusive flux of kinetic energy in Eq. (5). Let q denote the diffusive flux of internal energy. As for the momentum balance, we take into account an external source of energy given by ρu · g and assume to have no source or sink of internal energy. Inserting these considerations in Eq. (5) yields Expanding the terms in Eq. (11) and using relations Eqs. (10) and (7) leads to ρ ∂ t e u + u · ∇e u + ∇ · q = τ : (∇u) According to Fourier's law q = −λ∇ϑ, with the temperature ϑ and λ denoting heat conductivity. In view of ρe u = ρc p ϑ ρc p ∂ t ϑ + u · ∇ϑ − λ∆ϑ = τ : (∇u) (13) holds in the bulk of each phase.

2.2.4.
Scale analysis for small Mach numbers and further simplifications. In the gas phase the ideal gas law, is assumed to hold. Hereby R s denotes the specific gas constant. For ideal gases, the energy is proportional to the thermodynamic pressure p td : where c V denotes heat capacity. Furthermore, the speed of sound is given by γp ρ (γ denoting the adiabatic exponent) see for instance [35]. Since in our case the magnitude of the velocity is far below the speed of sound, the scale analysis of [27] for small Mach numbers applies. From this analysis one infers that the thermodynamic pressure in the gas is constant in space (however, may vary in time) and that the fluid flow is quasi-incompressible (or weakly compressible). Since no chemical reactions are considered and the temperature differences are small, the divergence of the flow field in the gas is a constant in space compression rate accounting for the change of mass in the gas phase by evaporation/condensation in (our) case of a closed system. The compression rate in turn is determined by the flow field in the liquid phase and the mass flux.
The thermodynamic pressure is calculated from an ordinary differential equation involving the compression rate.
With the help of (6) one derives the jump conditions on the interface: where j denotes the mass flux and Λ the latent heat of evaporation. The temperature on the interface is assumed to be the (pressure dependent) saturation temperature and for the mass flux we have Alternatively and physically more sound, one can use the Hertz-Knudsen law derived from kinetic theory, see Eq. (41) below. However, since over/undercooling is very small in our applications, the above relation is a good approximation.
In the following, we assume that stresses acting from the flow in the gas phase can be neglected in comparison with stresses acting from the flow in the liquid phase, i.e. |T g n| |T l n|. As a consequence, and the flow problems in the liquid and gas phase are coupled only by heat and mass transfer. We also neglect the term τ : (∇u) in Eq. (13), which describes the transformation of kinetic energy into heat by friction.
More details are given below in Subsection 2.3.
ρ0,g . The dependent variables are made dimensionless by settinĝ Then the equations for the system in dimensionless form read as follows. In the bulk of the liquid phase Ω l it holds: with e g the unit vector in direction of gravity. On the phase boundary Γ S the following conditions hold:T with p v the partial pressure of the vapor. For the gaseous phase one gets (using the asymptotic results for Ma → 0): Note that the scale analysis yields two different pressures in the gas phase: the thermodynamic pressure p td and the (hydrodynamic) pressure p. The latter is a Lagrange multiplier for the constraint on the divergence of u in Eq. (30). The ordinary differential equation for the scaled pressure reads: and the dimensionless form of the compression rate in case of a closed system is calculated byĉ The boundary conditions for the gas flow at the interface are given by: Note that this implies no slip of the flow field in tangential directions. The saturation temperature in Eq. (28) and Eq. (36) can be calculated from the mass fraction of vapor and the scaled pressure bŷ with Finally, the mass flow through the interface follows the relation As mentioned above, one may replace Eq. (36) and Eq. (28) by modelling the deviation of the temperature on the interface to the saturation temperature as a result of evaporation. In this case, the so called Hertz-Knudsen model is obtained, relating the local temperature ϑ Γ S on the interface and evaporation rate j (see for instance [26] where m denotes the molar mass, k B the Boltzmann constant, p sat (ϑ Γ S ) the saturation pressure corresponding to the saturation temperature ϑ Γ S via the Clausius-Clapeyron equation and α the evaporation coefficient. Let us note that for moderate evaporation rates as they occur in our applications, the difference between the saturation temperature ϑ s and ϑ Γ S is usually small.
3. Numerical methods. In this chapter our numerical approach to solving the two-phase model introduced in Section 2 is presented. In what follows, the main ingredients of the method are introduced and described in detail. The entire approach is then summarized in the last Subsection 3.4 of this chapter. 3.1. ALE formulation of the model problem. In view of the application we have in mind, the interface will not change its topological type. Thus a mesh moving method, i.e. a method where the mesh is always aligned with the interface (or equivalently, the interface is explicitly tracked by the mesh), is appropriate. Besides its simplicity this methods also yields a high numerical accuracy. In order to formulate this mesh moving method, so called arbitrary Lagrangian-Eulerian (ALE) coordinates are used (see also [7,21,23]). To this end, consider a reference domainΩ g ∪Ω l =Ω ⊂ R d and the corresponding reference mapping ξ : R×R d → R d such that ξ(t,Ω) = Ω(t) and ξ(t,Γ S ) = Γ S (t) for the (fixed) reference interfaceΓ S , ruling out the occurrence of topological changes. We assume ξ to be sufficiently smooth and to be a homeomorphism for each time instant t. The domain velocity is defined by For any sufficiently smooth function F : In what follows, occurring time derivatives in the balance equations are replaced by the ALE time derivative Eq. (43). We will comment on the construction of the mapping ξ in Subsection 3.3.
3.2. The CFD solver Navier. Our implementation is based on the finite element CFD solver Navier, whose key features will be briefly described in the following subsection. A detailed description of its components can be found in [22], and an overview of Navier's structure can be found in [1].
3.2.1. Space discretization. The space discretization of Navier is based on the P2/P1 Taylor-Hood finite element on simplicial meshes, i.e. the trial space for the velocity consists of (globally continuous in either phase) piecewise quadratic polynomials while the pressure is approximated by (globally continuous in either phase) piecewise linear polynomials. This pairing yields an LBB-stable discretization of the Stokes equations, see for instance [17]. The interface Γ S is represented by (topologically) fixed edges of the (moving) triangulation of Ω, see Fig. 2 for a sketch. Furthermore, isoparametric, quadratically deformed mesh elements are used to increase accuracy of the interface representation.

3.2.2.
Variational treatment of interface forces. One of the crucial points of both space and time discretization is the treatment of the capillary boundary condition Eq. (26). This condition can be implicitly integrated into the weak formulation of the momentum equation Eq. (23) rendering an explicit evaluation of the curvature of Γ S unnecessary. To this end, following a technique proposed by [9], note that the curvature vector can be represented as where ∆ S denotes the Laplace-Beltrami operator and ∇ S = (Id − nn T )∇ denotes the tangential gradient. In order to simplify notation, we restrict the presentation to the isothermal case and assume Dirichlet boundary conditions for the velocity on ∂Ω l \ Γ S in Eqs. (23), (24).
The variational formulation of the core problems then reads: find u ∈ X such that Hereby, X and Y denote suitable function spaces for velocity and pressure and f = Bo We e g − Ra Re 2 P r ϑe g .
Integrating by parts the stress term yields The stress tensor T is symmetric and since ϕ vanishes on the Dirichlet boundary, inserting boundary condition Eq. (26) into Eq. (48) yields In the non-isothermal case an additional temperature dependent term for the Marangoni stresses (which is treated explicitly, see also Subsection 3.4) is obtained. Further integration by parts yields 1 Here, ∂Γ S denotes the contact line between interface and wall, n S denotes the outward pointing normal of the interface, i.e. the unit vector of the interface's tangent plane perpendicular to the contact line. With τ denoting the unit vector in the wall's tangent plane which is perpendicular to the interface, the condition holds. As a boundary condition for the shape of the surface, the contact angle γ is prescribed. In our applications, γ will be a fixed, static contact angle, but dynamical contact angle conditions may also be treated in the same way. See for instance [14] for an assessment of dynamical contact angles. Furthermore, in order to allow for a movement of the interface, a slip boundary condition for the velocity in a small vicinity around the contact line is prescribed.
The advantage of the variational strategy arises in the right hand side of Eq. (50), yielding a term which is well defined on polygonally bounded domains for which a point-wise evaluation of the curvature does not make any sense.
In terms of time discretization, two dependencies enter the variational formulation of the interface force in Eq. (50): the integration domain Γ S , but also the transformation ξ depend on time. A fully explicit treatment leads to a severe CFL condition of τ ≤ C √ We h 3/2 , where τ denotes time step and h denotes mesh size (see [3]). A fully implicit treatment, on the other hand, would lead to high computational extra costs. In Navier, these terms are treated in a semi-implicit way: in each time step, the transformation ξ is treated implicitly with integration performed on the current geometry. To be more precise, denote by . (n) quantities for the nth time step. For the (n + 1)st time step, one computes Hereby, we made use of the kinematic boundary condition Eq. (27), i.e. the fact that the interface moves according to the fluid velocity. All terms on the right hand side of Eq. (52) can be evaluated on the current geometry at time step n.
The semi-implicit treatment of the curvature only yields additional contributions in the system matrices and right hand sides for the computation of u n+1 (which, due to the time dependence of the integration domain have to be re-assembled in each time step). The main benefit lies in the decoupling of the geometry update and computation of the fluid velocity, enabling for an efficient solution strategy. Unconditional stability of the scheme (in the context of a space-time finite element discretization) was proved in [2].

Time discretization.
In each time step, we will have to solve for the unknowns u n+1 and p n+1 on the domains Ω n l and Ω n g . As the main time discretization scheme, Navier uses the fractional step-θ scheme which was introduced in Bristeau et al. [4]. To this end (and to simplify notation) consider the following flow problem for the velocity u and pressure p on a domain Ω and time interval [t n , t n + τ ]: with corresponding right hand sides f and g plus boundary conditions. In our model, g = 0 in the liquid phase and g = −c in the gaseous phase. Define θ = 1 − 1 √ 2 ≈ 0.29289 and divide the time interval into three subintervals [t n , t n + θτ ], [t n +θτ, t n +(1−θ)τ ] and [t n +(1−θ)τ, t n +τ ]. Let u (n) be the velocity of the previous time step. The flow problem Eq. (53) is discretized in these time subintervals in the following way: In [t n , t n + θτ ] solve for u (n+θ) , p (n+θ) : In [t n + θτ, t n + (1 − θ)τ ] solve for u (n+1−θ) : Finally, in [t n + (1 − θ)τ, t n + τ ] solve for u (n+1) , p (n+1) : Hereby, α = (1 − 2θ)/(1 − θ) and β = 1 − α. This approach decouples the two numerical main difficulties of the Navier-Stokes equations: In the first and last sub-step, the nonlinearity is treated explicitly and a Quasi-Stokes problem has to be solved, see Subsection 3.2.4. In the second substep, the pressure is treated explicitly and the divergence condition for the velocity is neglected. In this step, a non-linear equation for the velocity has to be solved, see Subsection 3.2.5. Further details can be found in [1,4].

3.2.4.
Solution of the Quasi-Stokes system. Let X and Y be suitable function spaces for velocity and pressure and X resp. Y their corresponding dual spaces. The variational formulation of the Quasi-Stokes problems Eq. (54) and Eq. (56) can be written in the following way: find u ∈ X, p ∈ Y such that: Hereby the linear operators A : X → X , B : X → Y and B T : Y → X are defined as The linear forms l ∈ X and m ∈ Y are defined as wheref andg correspond to the right hand sides of Eqs. (54) and Eq. (56). The problem Eq. (57) is equivalent to the Schur complement formulation: For an LBB stable space discretization by a finite element pair (X, Y ) the Schur complement operator BA −1 B T is symmetric and positive definite. Therefore, problem Eq. (63) can be solved efficiently using a (preconditioned) conjugate gradient method. A variant of this method which solves Eqs. (63) and (64) in parallel, is described in [4]. The method can be easily modified to account for the non-homogeneous divergence condition Eq. (30) which occurs for the flow problem in the gaseous phase. The only difference lies in the presence of the additional term m on the right hand side of Eq. (63).

3.2.5.
Solution of the non-linear subproblem. The non-linear subproblem Eq. (55) in the second step of the fractional step-θ can be formulated in the following way: find u ∈ X such that Au + N (u, u) = l.
(65) Hereby the operators A : X → X , N : X × X → X are defined as The linear form l ∈ X is defined by wheref corresponds to the right hand side of equation Eq. (55).
In Navier, problem Eq. (65) is solved by the following fixed point iteration: the kth iterate u (k) is obtained from The inner problem (69) is linear and solved by a GMRES method with restart (see, for instance, [41]). In each GMRES restart, the convection operator is updated using the new velocity iterate. In Navier, the nodes located on the discrete surface discretizing u n Γ are shifted using Eq. (70) to approximate Γ (n+1) S . In order to prevent the computational mesh from degenerating, this movement has to be extended to the whole computational domain. This can be done in different ways (see for instance [50,52] for a study on different extension techniques for problems with rather large deformations), and in Navier we use one of the following two options: • Harmonic extension: for the coordinates ξ of the computational domain, a Laplace problem − ∆ξ = 0 (71) is solved and Dirichlet values for the nodes on the interface, i.e. ξ = x for x ∈ Γ (n+1) S , are prescribed.
• An extension operator based on a non-linear variational approach to obtaining optimal meshes, derived from principles of non-linear elasticity theory. This approach was introduced and analyzed by M. Rumpf in [40]. The main idea is to find minimizers of a certain class of functionals, which measure the quality of mesh deformation. While the first option is preferable, when deformations of the interface and the mesh can be expected to be small, some of our applications (see Subsection 4.2) require the treatment of large deformations. In this case, the non-linear variational approach is used.

3.4.
Numerical method: Overview. Using the ingredients presented in the previous subsections, we are now able to summarize our numerical approach to solving the model presented in Section 2.
Let τ be the time step size, Ω (n) = Ω Recall that due to assumption (21), fluid forces from the gas phase acting on the interface are neglected. This allows us to decouple the flow problem into a free boundary flow problem for the liquid phase (see Steps 4 and 5 below), and a flow problem with a fixed boundary for the gas phase (Step 7 below).
At each time step t (n+1) the following subproblems are solved: 1. From ω (n) and p in Ω l , (73) in Ω g and in Ω s , using the Dirichlet values S . R c,s (R c,g ) denotes the ratio between the solid's (gas) and fluid's heat capacities. 3. On the interface Γ 7. Flow problem in the gaseous phase: solve for u : using boundary condition Eq. (37) on the interface.
Remark 1. Some further remarks on the scheme seem to be in order: • In view of a linear stability analysis which was presented in [31] showing that evaporation has a destabilizing effect on the liquid/gas interface, terms involving evaporation are treated explicitly. In contrast, surface tension, which stabilizes the interface, is treated (semi-) implicitly. • The applications from space industry that we have in mind are mainly concerned with cylindrical tank geometries. Thus, we furthermore assume the flow to be axially symmetric and solve all fluid equations in a rotationally symmetric formulation. The weak formulation is derived with appropriately r-weighted inner products and unknowns to avoid the singularities of cylindrical differential operators, see [49]. • Most of the subproblems are convection dominated, and therefore a streamline diffusion formulation (see for instance [5]) to solve the problems for temperature and mass fraction in Steps 2 and 8 is used.
4. Numerical results. In this section computational results for the method described above are presented for two different applications and experiments.

4.1.
Surface deformation by thermocapillary convection: SOURCE. The dynamic behavior of liquids in tanks subject to reduced gravity (i.e. microgravity conditions characterized by very small Bond numbers), where capillary forces dominate, can be strongly influenced by the contact angle between the interface and the container wall, see for instance [8,14,36]. Despite being a material property in a first approximation, the apparent contact angle, i.e. the observed macroscopic contact angle, can be influenced by thermal conditions at the wall. It is known from experiments performed under microgravity that the apparent contact angle of a cold liquid spreading over a hot wall increases [15]. These experiments and effects were described and simulated numerically in [29,30,32]. In this work, details about the numerical simulation of an experiment which was conducted aboard a sounding rocket, the so called SOURCE experiment (SOUnding Rocket Compere Experiment) are provided. A detailed description of this experiment can be found in [11,12], and in the following we provide an outline of the setting for our numerical simulations. 4.1.1. Description of the experiment. Our numerical simulations were concerned with the evolution of the apparent contact angle and interface position during one stage of the SOURCE experiment. In this stage, which takes place under microgravity conditions, a cylindrical tank is partly filled with the experimental liquid, and the gaseous phase consists of a mixture of inert gas and the liquid's vapor.
The container wall is heated by a ring heat foil mounted on top of the cylindrical experiment cell, creating an axial temperature gradient towards the contact line. An approximation to this profile is depicted in Fig. 3 (left). In view of the heated wall, the temperature dependency of the surface tension leads to a tangential stress which, in turn, creates a capillary (Marangoni) convection, transporting warmer liquid from the wall towards the symmetry axis. This mechanism leads to a change of the apparent contact angle between the liquid interface and the wall.
In our simulations which are described in detail in [13], we were interested in the influence of a characteristic temperature difference between the liquid close to the wall and the wall itself (the driving force for the Marangoni convection). This temperature difference may be used in the definition of the dimensionless numbers from Subsection 2.3. The influence of different temperature differences on the apparent contact angle may then be studied in terms of these dimensionless numbers, see [13] for detailed results.

4.1.2.
Numerical results. Since in the SOURCE experiment, the driving mechanism behind the apparent contact angle change is thermocapillary convection, we simplified the model to take only the liquid phase into account. This excludes mechanisms such as heat and mass transfer between the two phases.
For the dimensionless numbers, a characteristic velocity scale U has to be specified. In the case of thermocapillary convection of boundary layer type, Ostrach proposed the following scaling for the characteristic velocity in [38] for a given characteristic temperature difference Θ, Despite the fact that this scaling was derived under simplified assumptions, our results in [13] show that it is also valid for the situation presented here. For our computations, an initial configuration of the interface under 0g conditions was computed using the Young-Laplace equation, where ξ denotes the coordinate vector and z Γ is an integration constant chosen in such a way that the volume of the domain Ω l corresponds to the liquid volume in the experiment. Furthermore, a boundary condition in terms of the contact angle γ has to be prescribed. To prevent a simulation break down due to degenerating mesh elements, this contact angle has to be bounded somewhat from 0 from below. In the case of the SOURCE experiment, a static contact angle of 5 • (the experimental liquid had a static contact angle of 0 • ) was used. Also, initial values for the temperature distribution in the cylinder wall were modelled according to measurements obtained from the experiment. Fig. 3 (left) shows a three-dimensional representation of the axially symmetric temperature distribution and a zoom of a locally adapted triangulation (right) needed to capture the viscous and thermal boundary layers of the problem.
During the simulation, the change of apparent contact angle and the temperatures defining the characteristic temperature difference were measured. In Figs. 4(a) to 4(c), the temperature distribution at three different time instants is shown, corresponding to the initial state (t = 2 s), after t = 25 s and at the final time

Reorientation behavior of cryogenic liquids.
In the second application, we are interested in the reorientation behavior of cryogenic liquids undergoing a step reduction of gravity in the presence of non-isothermal boundary conditions. The corresponding experiment was performed at the drop tower facilities at ZARM, University of Bremen. A detailed description of the experiment can be found in [34] and in the following only a brief outline is given.
4.2.1. Description of the experiment. As in the SOURCE experiment, a partly filled cylindrical container is considered. The container holds the test liquid, liquid argon (LAr), and is housed within a cryostat to maintain cryogenic conditions. In contrast to the SOURCE experiment, no additional inert gas is present in the gaseous phase, making it a single species system. In this case, the temperature at the liquid/gaseous interface corresponds to the saturation temperature Eq. (28) and is constant (or, if the Hertz-Knudsen model is used, almost constant for typical evaporation rates observed in the experiment). As a consequence, thermocapillary convection is not expected for this experiment. At the beginning of the experiment, the experiment capsule is subject to normal gravity conditions and the interface obtains its 1g equilibrium configuration. During a pre-heating phase before the release of the capsule, axial wall temperature gradients of different magnitude are established in the quartz wall of the container. Once the capsule is released ("dropped") and subject to approximately 4.7 s of microgravity conditions, these non-isothermal conditions influence the capillary driven reorientation behavior of the liquid.
In the experiments, the frequency and dampening of the interface oscillations as well as the evolution of thermodynamic pressure were found to correlate strongly (see [33,34]). Our goal was to verify this correlation numerically.

Numerical simulations.
A detailed description of the numerical setup is given in [33], and our findings are summarized in this subsection. The simulation presented here corresponds to an axial wall gradient of 1.73K/mm in the experiments presented in [34]. The initial 1g position of the interface was computed using the Young-Laplace equation Eq. (80). Although liquid argon is a completely wetting liquid (γ = 0 • ), a static contact angle of γ = 10 • had to be prescribed to prevent the triangulation from degenerating in the vicinity of the contact point.
The simulation was started at rest (u l = u g = 0) shortly before the step reduction of gravity and an initial temperature profile corresponding to measurements from the experiment was imposed. It is worth mentioning that an accurate approximation of the initial temperature distribution in the vicinity of the interface is essential to obtain a physically realistic evaporation rate j, see Eq. (40). Initially, the interface had a constant temperature corresponding to the saturation temperature of ϑ s = 90.8 K defined by the thermodynamic pressure of p td = 1, 044 hPa through equation Eq. (28). Within a layer of thickness δ Γ S parallel to Γ S , the temperature drops to the sub-cooled temperature of the liquid bulk, ϑ = 85.0 K, see Fig. 5(a). An estimate for the thickness δ Γ S can be found in [6], and within this layer we assumed to have a linear temperature profile.
At t = 0 s, gravity is "turned off" by setting Bo = Ra = 0, corresponding to the release of the experiment capsule in the drop tower. The temperature and velocity profiles during the reorientation are depicted in Figs. 5(a) to 5(d) at four different time instances. Subfigure 5(b) depicts the moment of the largest deflection of the z-coordinate of the wall contact line. We make the following observations: the initial distribution in the temperature layer below the interface is immediately distorted by the capillary driven movement of the interface. Furthermore, this movement induces a forced convection in the gaseous phase transporting "cold" gas at saturation temperature from the interface into the gaseous phase. Also a vortex forming underneath the interface is observed, see Fig. 5(d).
We now compare different characteristic quantities of the simulation with experimental data: Fig. 6(a) shows the evolution of the z-coordinates of the center point z cp , i.e. the z-coordinate of the interface at the symmetry axis and the evolution of the z-coordinate of the wall contact point z wp . While perfect quantitative agreement between numerical and experimental results cannot be expected due to our numerical restrictions (for instance, the static contact angle γ = 10 • ), qualitative agreement is clearly visible: the number and frequency of oscillations roughly agrees. Fig. 6(a) also depicts the evolution of the normalized thermodynamic pressure p td /p td (t = 0 s) over time. Again, qualitative agreement between numerics and experiment can be observed. Furthermore, there is a clear correlation between interface motion and pressure evolution.
To explain this behavior, Fig. 6(b) shows the evolution of the mass flux J = Γ S j over time. As can be seen, this mass flux strongly correlates with the position of the interface: as soon as the largest deflection of the interface in z-direction occurs     (at approximately t = 0.9 s), the interface is in contact with the hot wall and surrounding hot gas. Consequently, the evaporation rate in the vicinity of the contact point should be strongly magnified. That this is indeed the case shown by Fig. 6(c). Once the contact point recedes and reaches its minimum z-deflection, the evaporation rate also decreases. We also make the observation that although large local evaporation rates are observed close to the contact line, condensation occurs in a region close to the symmetry axis. Our numerical simulations also suggest that the region separating condensation and evaporation roughly follows the evolution of the forming vortex, compare Figs. 6(d), 5(d).

Conclusion.
In this article we have presented a finite element method for the simulation of non-isothermal two-phase free surface flow including surface tension, Marangoni effects and mass transport across the interface by evaporation/condensation. The method was applied to two examples from microgravity flow, for which experimental data were available. Comparing experimental and numerical results shows good agreement between both.