|
TrioCFD 1.9.8
TrioCFD documentation
|
|
TrioCFD 1.9.8
TrioCFD documentation
|
The Front-Tracking discontinu model aims to solve unsteady problems such as wall bubble growth and instabilities. The existing numerical methods all have defects that make this type of computation difficult or impossible. For this reason, a numerical method adapted to these problems was built: the mixed Front-Tracking/VOF method. At present, it mainly handles axisymmetric or two-dimensional simulations in order to concentrate efforts on the quality of the physics of the simulations, but the three-dimensional method is already being implemented.
The objective of this module is to numerically simulate the interactions of interfaces with small objects (optical probes, for example) and the phenomena related to boiling and the boiling crisis. The interaction of interfaces with small structures is a problem governed by contact-line dynamics and surface tension. During boiling, the following mechanisms may play a fundamental role and must be treatable numerically:
All these problems require accounting for a two-phase flow with phase change, where the surface tension is of the same order of magnitude as the inertia, viscous dissipation or gravity forces. Contact lines play an essential role that the numerical method must account for.
The review presented here is relatively brief, as such an exercise has already been carried out elsewhere (see [26], [43], [131]). Only the characteristics of the available methods are made explicit, in order to justify the choice of the starting method.
Lagrangian tracking methods. The approach that immediately comes to mind consists in describing the interfaces by a moving mesh. This leads to two classes of methods: Lagrangian and mixed Eulerian-Lagrangian tracking methods (called ALE methods, for Arbitrary Lagrangian Eulerian) on the one hand, and front-tracking methods on the other. In the former (see for example [117]), the whole fluid mesh deforms so that the interfaces are at each instant a line of the mesh. With such a formulation, the boundary conditions of mechanical stress, velocity or heat flux are natural, but the changes of interface topology are practically out of reach for reasons of geometric complexity. This method is used when very precise results are sought on relatively simple geometries.
Dual Eulerian/Lagrangian mesh methods. In the front-tracking method [174], [159], [43], a moving surface mesh represents the interfaces while the velocity, pressure, temperature and other volume quantities are discretized on a fixed mesh. In this approach, the boundary conditions at the interfaces are difficult to control. Certain properties accessible in a Lagrangian approach, such as the volume conservation properties or the energy balances on the phases, are very complicated to verify when the interfaces do not coincide with the mesh lines. These methods are therefore characterised by inaccuracies in the balances, compensated by a large number of tricks to make these inaccuracies acceptable. Much less heavy than the Lagrangian methods, the front-tracking methods benefit from all the experience acquired in single-phase flow simulation. The management of the surface mesh of the interfaces remains a major difficulty, in particular during topology changes, but the constraints on this mesh are much less severe than for the volume mesh of the Lagrangian methods.
Purely Eulerian methods. The most widespread and one of the most robust methods to date is the VOF method (Volume Of Fluid [79]). In this method, the density field is discretized on the same fixed mesh as the velocity, pressure and temperature. From this density field, the position of the interfaces in each cell is deduced by geometric algorithms. This position makes it possible to compute the effect of surface tension on the fluid and the density variations in each cell. The accounting for surface tension depends enormously on this reconstruction algorithm, and the most precise methods are just as complex to implement as the front-tracking methods. On the other hand, the tracking of interfaces is performed by a transport algorithm on the Eulerian mesh, which avoids a large part of the algorithmic difficulties related to the moving mesh. Among the other purely Eulerian methods, there are also the so-called "level-set" methods, where the phase indicator is a continuous discretized function and the interfaces are implicitly defined as the 0.5 isovalue of the phase indicator. One can also cite the method based on the second-gradient equations [27], in which the physical properties specific to the interfaces (such as surface tension or latent heat of phase change) are naturally accounted for by the evolution equations of the Eulerian variables.
The complexity/precision/efficiency compromise. The methods have been cited in decreasing order of numerical efficiency, precision and complexity. The Lagrangian methods are naturally the most precise because the discretization of the physical-quantity fields respects the topology of these fields: the mesh lines follow the discontinuities of the physical field. These methods are also the most complicated to implement. The mixed Eulerian-Lagrangian methods constitute an interesting compromise. The purely Eulerian methods are generally the simplest because they avoid the management of a moving mesh, but the discontinuities of the fields are represented by a quantity on the Eulerian mesh, and the localisation is therefore not very precise; these methods often require a finer mesh to obtain the same precision.
Duquennoy ([43], [13]) implemented a front-tracking method in which he more specifically improved the accounting for phase change and contact lines. His work highlighted several weak points of the front-tracking method:
These weak points constitute as many selection criteria for a method better suited to the problems considered.
Parasitic currents. Parasitic currents are the currents observed in a numerical simulation that has reached a stationary equilibrium state even though no energy is injected into the system. These currents result from discretization errors of the surface tension and have serious consequences on the results of the computations. The first consequence is that they make certain computations impossible. According to Duquennoy [43], the intensity of the parasitic currents increases when the value of the Ohnesorge number decreases (the length \(D\) being the characteristic dimension of the drops or bubbles):
\[Oh = \dfrac{\mu_{l}}{\sigma \rho_{l} D} \]
In practice, with water at atmospheric pressure and a characteristic dimension of the order of a centimetre, the parasitic currents make the interfaces completely unstable at long times. The second consequence is that they make anisothermal computations unrealistic where the heat transfer is not dominated by conduction. It is therefore extremely important to reduce the intensity of the parasitic currents at least to a level such that the Péclet number is much smaller than 1; in that case, the parasitic currents no longer have any influence on the temperature field.
Volume conservation. This property is a very important issue for the front-tracking method as for the VOF method. By nature, the VOF method has good volume-conservation properties, but certain steps of the algorithm are approximate and conservation is inexact in most existing algorithms [47]. The front-tracking method implemented by Duquennoy presents an even larger error than the VOF methods on the volume balance of the phases, on the one hand because of the spline interpolation method used to guarantee the regularity of the Lagrangian mesh, and on the other hand because of the transport algorithm of this mesh. The error on the volume balance has important consequences on the simulation of wall bubble growth: in some cases, the bubble volume decreases even though the fluid is superheated everywhere. Consequently, the accuracy of the phase volume balance is an indispensable ingredient of the numerical method.
Contact lines. The formulation used by Duquennoy currently prevents the simulation of angles smaller than about \(30^{\circ}\), and the contact-angle modelling is poorly controlled (the dynamic contact angle in particular is not modelled rigorously).
Coherence of the temperature field. The previous formulation of front-tracking does not exactly ensure the boundary condition \(T = T_{sat}\) at the interfaces. Thus, one can observe that a vapour bubble initially at \(T_{sat}\) heats up in contact with a superheated liquid, which contradicts physical sense. This incoherence, as well as the imprecision of the mass-transfer formulation at the interfaces, is not compatible with the implementation of a precise model of singular heat transfer at the contact lines.
The choice fell on a front-tracking method for the following reasons:
The original front-tracking method nevertheless presents many drawbacks, the most important being perhaps the presence of parasitic currents. This is why many aspects of the initial formulation were deeply modified, starting with the discretization approach of the continuous equations.
Discretization approach of the continuous equations. Fundamentally, this approach changes radically compared to the initial method [174] : the old approach consists in applying a filtering operator to the Navier-Stokes equations to make all the quantities continuous, then discretizing these filtered equations. The disadvantage is that one converges quickly towards the solution of the filtered problem and not towards the solution of the real problem; in this sense, the method is not even "consistent". The new approach is based on balances over control volumes, exactly as in a classical finite-volume method. The Navier-Stokes equations, including the discontinuities at the interfaces, are thus discretized directly. To obtain mesh convergence of the method, a subgrid model of the physical quantities must be introduced in the elements where these quantities are discontinuous. The convergence obtained is of order 1 if a zeroth-order subgrid model is used in the elements containing interfaces.
Summary of the essential contributions of the method. The essential contributions of the work carried out in TrioCFD concern:
The proposed method can be extended to three-dimensional schemes and to other types of Eulerian discretizations (for example the Finite-Element-Volume discretization [49], [73]).
The following notation conventions are used for the discrete quantities:
In two dimensions, the interface is defined by a mesh of nodes connected by segments. By convention, the segments are oriented so that the vapour is on the left. In three dimensions, the interface is a union of triangles whose normal is oriented towards the vapour. For the algorithms to work correctly, the mesh must verify certain topological properties: two interface segments never intersect; either the interfaces are closed or their ends (the nodes with only one connected segment) are located on a boundary of the domain; the interfaces define closed volumes, and all the interfaces defining the boundary of a volume must have their normals oriented in the same direction. These properties ensure the topological coherence of the mesh. The interfaces are considered as a succession of segments; the order of the method is not increased by the use of splines for the computation of the indicator, because the topological properties stated above are easier to verify with segments.
We call \(\Gamma\) the set of interfaces and \(E\) an interface element (segment in 2D and triangle in 3D). We sometimes use the notion of "complete connected interface portion": connectivity means that one can pass from one node to any other of the interface portion by traversing elements one by one; a portion is complete if, for any of its elements, its neighbours are also in it.
The curvature is computed at the nodes of the surface mesh. This choice is imperative so that the resulting surface tension only admits the "constant curvature" solution as the stationary solution of the problem. Two formulations of the discrete curvature can be considered: one based on geometric interpolation, and one based on the differential of the interface energy.
Geometric formulation. The first is the most intuitive in two dimensions: the geometric definition of curvature is used,
\[c(s) \mathrel{\hat{=}} \frac{\partial t}{\partial s} \cdot n \]
where \(s\) is the curvilinear abscissa in normal parametrization. Tracing the circle passing through the interface point and its two neighbours, the two-dimensional curvature is the signed inverse of the radius \(R\) of the circle (curvature positive if the centre of the circle is in the vapour, negative otherwise). In axisymmetric geometry, the curvature in the other direction must be added, which is written \(c_{2} = \sin\theta / x\).
Several steps of the algorithm require passing from a quantity defined at the interface nodes to a quantity defined on the fixed mesh and vice versa. We denote by \(\overline{\mathcal{G}}\) the operator passing from a surface quantity to a volume quantity, and \(\hat{\mathcal{G}}\) the reciprocal operator. For discrete fields \(\hat{f}\) and \(\overline{f}\):
\[\hat{f} = \hat{\mathcal{G}}(\overline{f}), \qquad \overline{f} = \overline{\mathcal{G}}(\hat{f}) \]
Important remark: these operators are not inverses of each other. If one passes from a surface quantity to a volume quantity, then again to a surface quantity, the integral is conserved but the nodal values undergo a numerical diffusion due to the successive interpolations:
\[\hat{\mathcal{G}} \left( \overline{\mathcal{G}} \left( \hat{f} \right) \right) \neq \hat{f} \]
The most important property verified by these operators is the conservation relation (where \(\Gamma\) represents a complete connected portion of interfaces and \(\Omega\) the set of Eulerian mesh elements crossed by \(\Gamma\)):
\[\int_{\Gamma} \tilde{f}(x)\, ds = \int_{\Omega} \overline{f}(\Omega)\, d\Omega \]
The discrete indicator \(\overline{I}\) of a volume element \(\Omega\) is defined as the fraction of the element volume occupied by the gas phase (the mean void fraction in the element). The density \(\overline{\rho}\) of an element is the mean density in that element, \(\overline{\rho} \mathrel{\hat{=}} \frac{1}{\overline{V_{\Omega}}}\int_{\Omega}\rho\,d\Omega\), and the following relations hold:
\[\overline{\rho} = \rho_{v}\, \overline{I} + \rho_{l}\, (1-\overline{I}), \qquad \overline{I} = \dfrac{\rho - \rho_{l}}{\rho_{v} - \rho_{l}} \]
The indicator is computed as a function of the position of the interface nodes by an exact geometric algorithm, identical to that proposed by Popinet [137], except that only the interface segments are used (not a spline interpolation).
Having a discretization of the density (through the indicator), one must choose whether to discretize the velocity \(v\) or the momentum \(\rho v\) — i.e. whether to favour volume conservation or momentum conservation. If the velocity is discretized, the incompressibility condition is very easy to write and translates exactly in discrete variables, \(\overline{\nabla \cdot \overline{v}} = 0 \Leftrightarrow \forall\Omega,\ \int_{\partial\Omega} \overline{v} \cdot \overline{n}\, ds = 0\), but a local mass balance and momentum conservation are much more complicated to obtain. If the momentum \(\overline{\rho v}\) is discretized, its conservation is easy to obtain but the incompressibility condition becomes locally ambiguous in regions where the density varies, since the velocity must be evaluated by a formula of the type \(\overline{v} \mathrel{\hat{=}} \overline{\rho v}/\overline{\rho}\). Such an error would have unfortunate consequences in front-tracking (and would increase with the density ratio): the interfaces would not be transported without deformation even if the solution is a uniform velocity field, modifying the surface energy of the interfaces and leading to parasitic currents or numerical instabilities. The choice retained for now is to use a velocity discretization.
In the current implementation, the velocity and pressure discretization is of "marker and cell" type. The pressure is discretized at the centre of the control volumes \(\Omega_x\) and \(\Omega_y\). The vertical faces carry a horizontal velocity component \(\overline{v_x}\), the horizontal faces a vertical velocity component \(\overline{v_y}\) (and similarly the third component in 3D). The velocity then has the two interpretations:
\[\overline{v_{x}} = \int_{\Gamma_{x}} v \cdot x\, ds, \qquad \overline{v_{x}} = \dfrac{1}{\rho} \int_{\Omega_{x}} \rho v_{x}\, d\Omega \]
The first allows writing the conservative mass balance on \(\Omega\) and defining the fluid incompressibility on the discrete equations, \(\int_{\partial\Omega} v \cdot n\, ds = \int_{\Omega} \nabla \cdot v\, d\Omega = 0\). The second allows writing a discrete conservative momentum balance, provided \(\overline{\rho}\) is constant.
Here again, the discrete quantity representing the internal energy of the fluid must be chosen. A discretization of the energy makes it easy to write a conservative scheme (important for phase change), but computing the Fourier heat flux and the evaporation requires estimating the temperature from the internal energy, which involves dividing by a quantity that tends to zero in the vapour, possibly leading to unbounded temperatures (in violation of the second principle of thermodynamics). If the temperature is discretized instead, the second principle is easy to verify, but energy conservation is more difficult to ensure in regions where \(\rho c_P\) varies. The current formulation uses a temperature discretization and is therefore not energy-conservative:
\[\overline{T} \mathrel{\hat{=}} \dfrac{\int_{\Omega} \rho c_{P} T\, d\Omega}{\int_{\Omega} \rho c_{P}\, d\Omega} \]
and, consistently, the heat capacity of the element is \(\overline{\rho c_{P}} \mathrel{\hat{=}} \frac{1}{\overline{V_{\Omega}}}\int_{\Omega} \rho c_{P}\, d\Omega = \rho_{v} c_{Pv} \overline{I} + \rho_{l} c_{Pl} (1 - \overline{I})\).
Saturation temperature. According to the continuous equations used, the local saturation temperature depends on the curvature of the interfaces and is expressed as a function of the pressures on either side of the interface:
\[T_{sat} = T_{sat}(P), \quad \text{with } P = P_{v} + \dfrac{\rho_{v}}{\rho_{l} - \rho_{v}} \sigma c + \dfrac{1}{2} P_{r} = P_{l} + \dfrac{\rho_{l}}{\rho_{l} - \rho_{v}} \sigma c + \dfrac{1}{2} P_{r} \]
With the half-sum of the two expressions of the pressure:
\[P = \dfrac{1}{2} \left( P_{v} + P_{l} + \dfrac{\rho_{l} + \rho_{v}}{\rho_{l} - \rho_{v}} \sigma c\right) \]
At equilibrium, the discrete pressure field \(\overline{P}\) verifies the equation constructed from the surface-tension source terms, \(\overline{P} = \overline{\kappa}\,\overline{I} - \overline{\rho}\,\overline{\phi} + P_{0}\), with \(\overline{\kappa} = \sigma c + (\rho_{v} - \rho_{l})\phi = \text{constant}\) and \(P_{0} = \text{constant}\). From this property, an expression of the pressure in the liquid and in the vapour on either side of the interface is deduced. Substituting these into the previous relation suggests using the following expression, exact at equilibrium, of the saturation temperature as a function of the discrete fields \(\overline{P}\), \(\overline{I}\) and \(\overline{\kappa}\):
\[\overline{T_{sat}} \mathrel{\hat{=}} T_{sat}(P), \quad \text{with } P \mathrel{\hat{=}} \overline{P} + \overline{\rho}\,\overline{\phi} + (1 - I)\overline{\kappa} \]
If \(\overline{\kappa}\) is constant, the discrete system can be simultaneously at mechanical and thermal equilibrium. Unfortunately, this introduces an extremely strong coupling between the temperature equation and the momentum equation; without a particular treatment (implicit treatment of the saturation-temperature variation), the numerical scheme is very unstable.
Phase-change enthalpy. The volume power of phase change on a control volume \(\Omega\) is the integral of the heat flux from the phases towards the interfaces contained in this volume, \(\overline{h} \mathrel{\hat{=}} \frac{1}{\overline{V_{\Omega}}}\int_{\Gamma\cap\Omega} \dot{q}\, ds\). Rather than interpolating the temperature beyond the interface (which is problematic near the domain boundaries), a formulation closer to the VOF reasoning is preferred: from the discrete temperature \(\overline{T}\) of an element, a subgrid model of the continuous temperature field is reconstructed using the position of the interfaces. The currently implemented subgrid model is simple and not very precise; the heat flux at the interface is modelled as
\[\dot{q} \mathrel{\hat{=}} \dfrac{k(\overline{T}-\overline{T_{sat}})}{\delta_{R}+\delta} \]
where \(k\) is the thermal conductivity of the liquid, \(\delta_{R}\) the equivalent interface resistance thickness and \(\delta\) a characteristic dimension equal to one quarter of the element size. The volume power of phase change is then \(\overline{h} \mathrel{\hat{=}} \dot{q} \cdot \text{surface}(\Gamma \cap \Omega) / \overline{V_{\Omega}}\), and the resulting temperature variation over a time step \(\Delta t\) is \(\Delta \overline{T} = -\overline{h}\, \Delta t / \overline{\rho c_{P}}\). This value must be bounded to ensure the time stability of the scheme: \(|\overline{h}| \leq \overline{\rho c_{P}}\, |\overline{T} - \overline{T_{sat}}|\, \Delta t\).
An even cruder model can be used: an infinite thermal conductivity in the element, so that the elements containing an interface always have a temperature \(T = T_{sat}\). Its spatial precision is not much lower (consistent and first-order in space), but the heat flux is discontinuous in time when the interface enters a new element, causing velocity and pressure jolts.
Mass flux at the interfaces. The mass transfer on the Eulerian mesh is defined as \(\hat{\dot{m}} \mathrel{\hat{=}} \frac{1}{\overline{V_{\Omega}}}\int_{\Omega \cap \Gamma} \dot{m}\, ds\). Since \(\dot{q} = \mathcal{L} \dot{m}\), the relation between \(\overline{h}\) and \(\overline{\dot{m}}\) is \(\overline{\dot{m}} = \overline{h}/\mathcal{L}\). For reasons of numerical stability of the interfacial-tension energy, it is necessary to reduce the irregularities of the \(\overline{h}\) field by a slight conservative spatial filtering of the mass flux (passing from \(\overline{\dot{m}}\) to \(\hat{\dot{m}}\), then back to \(\overline{\dot{m}}\)). The singular contribution \(\hat{\dot{m_{s}}}\) of the contact lines and the contribution \(\overline{\dot{m_{f}}}\) of the evaporation of wall liquid films are added:
\[\hat{\dot{m}} \mathrel{\hat{=}} \dfrac{1}{\mathcal{L}} \hat{\mathcal{G}}(\overline{h}) + \hat{\dot{m_{s}}}, \qquad \overline{\dot{m}} \mathrel{\hat{=}} \overline{\mathcal{G}}(\hat{\dot{m}}) + \overline{\dot{m_{f}}} \]
For an element \(\Omega\) crossed by an interface \(\Gamma\), with \(\Omega_l\) and \(\Omega_v\) the liquid and vapour volumes, the continuous velocity field satisfies
\[\int_{\partial\Omega} v \cdot n\, ds = \int_{\Gamma} (v_{v} - v_{l}) \cdot n\, ds = \int_{\Gamma} \dot{m} \left( \dfrac{1}{\rho_{v}} - \dfrac{1}{\rho_{l}}\right) ds \]
This translates into the definition of the discrete velocity divergence in the element \(\Omega\):
\[\overline{\nabla \cdot \overline{v}}(\Omega) \mathrel{\hat{=}} \overline{\dot{m}} \left( \dfrac{1}{\rho_{v}} - \dfrac{1}{\rho_{l}}\right) \]
The projection scheme then ensures \(\int_{\partial\Omega} \overline{v} \cdot n\, ds = \overline{\nabla \cdot \overline{v}}(\Omega)\, \overline{V_{\Omega}}\) for any volume \(\Omega\) made of Eulerian mesh elements.
Using the indicator relation, the variation of the gas volume in a volume \(\Omega\) is \(\partial_{t}\int_{\Omega} I\, d\Omega = -\frac{1}{\rho_{v}-\rho_{l}}\int_{\partial\Omega}\rho v \cdot n\, ds\). Inspired by this, the discrete quantity \(\overline{V'}\) is defined on the Eulerian mesh elements, with the density \(\rho\) on a face computed as \(\rho_l\) (resp. \(\rho_v\)) if a neighbour element is full of liquid (resp. vapour), and \((\rho_l + \rho_v)/2\) otherwise (the choice in the third case does not affect the conservation properties sought).
The goal of the (somewhat technical) phase-mass-balance derivation is to obtain a geometric condition on the displacement of the interface nodes so that the phase mass balance is respected exactly — an essential and non-trivial ingredient of the method. Two ingredients are essentially used: the exact volume balance on the discrete velocity field in the phases, provided by the projection scheme, and the conservative distribution of the mass fluxes \(\dot{m}\) on the volumes containing the interfaces. The result is semi-local: to ensure phase volume conservation, the geometric volume variation of each connected interface portion over the time step must equal the discrete integral of \(V'\). To this end, a volume variation \(V'_i\) is attributed to each interface node (via the conservative passage relations), and the displacement of each node must contribute to its associated volume variation. The discrete phase volume balance is exact only if the contact lines are located on closed boundaries of the domain. This algorithm provides not only an exact global volume conservation but also a "quasi-local" conservation, and it implicitly accounts for the mass fluxes at the domain inlets and the volume variations of the other interfaces.
1.16.1