Height Function Method Research Articles

Highly scalable DNS solver for turbulent bubble-laden channel flow

We present an efficient and highly scalable solver for direct numerical simulation (DNS) of dispersed gas-liquid flow, containing a large number of deformable bubbles. We apply this to O(104) bubbles in a turbulent flow. This was accomplished by a well-considered combination of state-of-the-art numerical methods, as well as fast and scalable numerical algorithms that have their origin in single-phase and two-phase flows. The features key-elements of the algorithm of curvature computation, bubble collision and variable-coefficient Poisson equation solver are presented. Resolution requirements for an accurate advection of a rising bubble have been established by comparison with the (Hysing, 2009) benchmark. The Generalized Height Function (GHF) method, was adopted for the curvature computations. We observed agreement with theoretical convergence rates, as well as a significant reduction of spurious velocities when employing GHF, compared to a finite difference approach. The main interaction mechanisms between bubbles and walls of the domain were analyzed, showing second order convergence of the underlying numerical methods.A detailed analysis of the parallel performance of the Navier–Stokes (NS) solver and the solver for the gas volume fraction was carried out. The analysis revealed close to linear scaling up to ≈ 18000 cores on a computational grid with 1 billion cells for the NS solver, and an ideal scaling of the gas volume fraction solver up to O(103) bubbles. Beyond that, acceptable overhead for up to O(104) bubbles was found.A simulation of a downflow configuration of a turbulent channel loaded with a total of 10000 bubbles illustrates the computational capabilities. First and second order statistics of the velocity field were computed as well as the profiles of the average gas volume fraction field in the statistically steady state. In the case considered, a 47% increase of the wall shear stress was observed, brought about by turbulence modification arising from the embedded bubbles. The interaction between turbulence and bubbles at high volume fraction resulted in a strong attenuation of the rms of the velocity fluctuations in a wide region of the core of the channel.

A physical modeling and predictive simulation of the laser cladding process

A 3D finite element simulation model of the laser cladding process has been developed taking into account heat transfer, fluid flow, surface tension and free surface movement. All input parameters and data, which are independent of the process parameters but depend only on the material and machine properties, have been obtained from measurements. Thereby the melt pool and the resulting surface contour can be simulated without compromising assumptions or calibration, because the machine parameters are the only variable input parameters of the model. Thus, the model can easily be transferred to other material combinations or other machines. For the surface contour calculation a modified height function method is applied. The model surface follows this contour as an arbitrary Lagrangian Eulerian (ALE) method is used allowing for mesh deformations. The model was implemented using the commercial finite element software COMSOL Multiphysics and validated by comparing the simulation results with caloric measurements of the effective heat input and metallographic cross sections from experiments, where the nickel-base alloy MetcoClad® 625 in powder form was deposited on structural steel S235JRC + C and the process parameters of laser power, feed speed, laser beam spot size and powder mass flow were varied within a range of at least 50% of their mean value each. The maximum deviation of the simulation results compared to the experimental data regarding track geometry is 14% for the parameter sets without weld defects so that these parameter sets could be industrially applied, whereas the average deviation of track width and height is below 5.1%. Moreover, the maximum difference of the simulated absorptivity compared to the measurement results is 6.9%, while the average difference is only 1.7%, meaning the model shows a high predictive capability especially regarding the heat input.

Importance of curvature evaluation scale for predictive simulations of dynamic gas–liquid interfaces

The effect of the scale used to compute the interfacial curvature on the prediction of dynamic gas–liquid interfaces is investigated. A new interface curvature calculation methodology referred to herein as the Adjustable Curvature Evaluation Scale (ACES) is proposed. ACES leverages a weighted least squares regression to fit a polynomial through points computed on the volume-of-fluid representation of the gas–liquid interface. The interface curvature is evaluated from this polynomial. Varying the least squares weight with distance from the location where the curvature is being computed, adjusts the scale the curvature is evaluated on. ACES is verified using canonical static test cases and compared against second- and fourth-order height function methods.Simulations of dynamic interfaces, including a standing wave and oscillating droplet, are performed to assess the impact of the curvature evaluation scale for predicting interface motions. ACES and the height function methods are combined with two different unsplit geometric volume-of-fluid (VoF) schemes that define the interface on meshes with different levels of refinement. We find that the results depend significantly on curvature evaluation scale. Particularly, the ACES scheme with a properly chosen weight function is accurate, but fails when the scale is too small or large. Surprisingly, the second-order height function method is more accurate than the fourth-order variant for the dynamic tests even though the fourth-order method performs better for static interfaces. Comparing the curvature evaluation scale of the second- and fourth-order height function methods, we find the second-order method is closer to the optimum scale identified with ACES. This result suggests that the curvature scale is driving the accuracy of the dynamics. This work highlights the importance of studying numerical methods with realistic (dynamic) test cases and that the interactions of the various discretizations is as important as the accuracy of one part of the discretization.

The evolution of hyperboloidal data with the dual foliation formalism: mathematical analysis and wave equation tests

A long-standing problem in numerical relativity is the satisfactory treatment of future null-infinity. We propose an approach for the evolution of hyperboloidal initial data in which the outer boundary of the computational domain is placed at infinity. The main idea is to apply the ‘dual foliation’ formalism in combination with hyperboloidal coordinates and the generalized harmonic gauge formulation. The strength of the present approach is that, following the ideas of Zenginoğlu, a hyperboloidal layer can be naturally attached to a central region using standard coordinates of numerical relativity applications. Employing a generalization of the standard hyperboloidal slices, developed by Calabrese et al, we find that all formally singular terms take a trivial limit as we head to null-infinity. A byproduct is a numerical approach for hyperboloidal evolution of nonlinear wave equations violating the null-condition. The height-function method, used often for fixed background spacetimes, is generalized in such a way that the slices can be dynamically ‘waggled’ to maintain the desired outgoing coordinate lightspeed precisely. This is achieved by dynamically solving the eikonal equation. As a first numerical test of the new approach we solve the 3D flat space scalar wave equation. The simulations, performed with the pseudospectral bamps code, show that outgoing waves are cleanly absorbed at null-infinity and that errors converge away rapidly as resolution is increased.

Computing interface curvature from volume fractions: A hybrid approach

The Volume of Fluid method is extensively used for the multiphase flows simulations in which the interface between two fluids is represented by a discrete and abruptly-varying volume fractions field. The Heaviside nature of the volume fractions field presents an immense challenge for the accurate computation of the interface curvature and induces the spurious velocities in the flows with surface-tension effects. A 3D hybrid approach is presented combining the Convolution and Generalized Height Function method for the curvature computation. The volumetric surface tension forces are computed using the balanced-force continuum surface force model. It provides a high degree of robustness at lower grid resolutions with first-order convergence and high accuracy at higher grid resolutions with second-order convergence. The present method is validated for several test cases including a stationary droplet, an oscillating droplet and the buoyant rise of gas bubbles over a wide range of Eötvös (Eo) and Morton (Mo) numbers. Our computational results show an excellent agreement with analytical/experimental results with desired convergence behavior.

Simulation of microchannel flows using a 3D height function formulation for surface tension modelling

Simulations of multiphase flow and phase change heat transfer in microchannels require accurate calculation of the surface tension force to provide accurate interface locations and to avoid spurious velocities that distort the flow at the interface. Building on previous work, where we implemented a 2D height function technique, the 3D height function method has been modified and implemented in ANSYS Fluent to enable such calculations to be performed on structured meshes. Firstly, simulations are presented for test cases that demonstrate the correct and accurate calculation of the surface tension force. Then the model is used to investigate both Taylor bubble formation in a square section channel and the development of an elongated vapour bubble from a small spherical nucleated bubble in a heated channel.

Estimation of curvature from volume fractions using parabolic reconstruction on two-dimensional unstructured meshes

This paper proposes a method to estimate the curvature of an interface represented implicitly by discrete volume fractions on an unstructured two-dimensional mesh. The method relies on the computation of local parabolic reconstructions of the interface. The parabolic reconstruction of the interface in a given computational cell is obtained by solving a local non-linear minimisation problem, and only requires additional information from two neighbouring cells. This compactness ensures a robust behaviour on poorly-resolved interfaces. The proposed method is proven to be analogous to the height-function method for Cartesian configurations with consistent heights, and can be interpreted as a generalisation of the height-function method to meshes of any type. Tests are conducted on a range of interfaces with known curvature. The method is shown to converge with mesh refinement with the same order of accuracy as the height-function method for all three types of meshes tested, i.e. Cartesian, triangular, and polygonal.

HFES: A Height Function Method with Explicit Input and Signed Output for High-Order Estimations of Curvature and unit Vectors of Planar Curves

The author proposes a height function (HF) method for estimating curvature and unit vectors of planar curves. The new method is named HFES since its signature consists of explicit input parameters (instead of implicit volume fractions) and signed output estimates. Due to this signature and its formal description, HFES has a number of appealing advantages. First, the convergence rate of HFES can be 4, 6, or a greater even integer. Second, local arcs can be rotated to further reduce leading truncation errors. Third, the sampling size of estimation is completely independent from the Eulerian grid size of the main flow. Fourth, signs of the signed output can be conveniently chosen by orientations of input Jordan curves, e.g., surface tension of a tracked material can be determined in a simple and worry-free manner. Last, HFES is applicable both to numerical simulation of multiphase flows and shape analysis in computer imaging. Furthermore, the author analyzes the effect of input perturbations upon output estimates, derives expressions of the best attainable accuracy afforded by a given input on uniform grids, and explains how to achieve this best accuracy with a concrete example. Results of numerical experiments with both exact and inexact input demonstrate that HFES can be much more accurate than previous HF methods. For the vortex-shear test of a circular disk, the eighth-order HFES coupled with the recent cubic iPAM method is millions of times more accurate than a fourth-order HF method coupled with volume-of-fluid methods! This drastic accuracy improvement is mainly due to the purely explicit spline representation of the curve and the rotation of local arcs.

REMOVED: An adaptive mesh refinement based flow simulation for free-surfaces in thermal inkjet technology

Abstract The microscopic multiphase flows of thermal inkjet (TIJ) devices usually involve a constantly changing liquid–gas interface. Often the curvature of the interface can become too high at some localized positions to be accurately resolved using a conventional static computational grid. A good solution to such challenging problem is to use the Adaptive-Mesh-Refinement (AMR) method, since AMR technique can concentrate the computational effort in the regions where it is needed most. In recent years, an open-source AMR code, Gerris , has been extensively used and benchmarked for a large variety of surface-tension dominated flows. However, the code cannot be directly used to model the TIJ process due to the lack of bubble dynamics and contact angle model for solid boundaries in the code. Therefore in the paper, we propose the following new approaches to model the TIJ flows and implement these algorithms in Gerris : (1) the vapor bubble is treated as an empty cavity so that its complex behavior of the vapor bubble after the bubble nucleation can be modeled using the adiabatic gas expansion and compression; (2) a height-function method based contact angle model is developed and implemented on the embedded solid surfaces. The new code has been verified and validated through a series of numerical examples and experimental tests. The proposed AMR method and the simulation code can be used in a broad range of applications in TIJ technology as well as multiphase flows of microfluidics.

Implementation of a height function method to alleviate spurious currents in CFD modelling of annular flow in microchannels

The use of the continuum surface force (CSF) model in surface tension modelling using the volume of fluid (VOF) method is known to introduce unphysical velocities, called parasitic or spurious currents arising mostly from poor calculation of interface curvature. In this paper, a modified 2-D axisymmetric height function (HF) method is proposed and implemented in parallelised computations. The default CSF method in the commercial software ANSYS Fluent 15.0 is replaced using User-Defined Functions (UDFs) of ANSYS Fluent and the improvements obtained are evaluated in test cases with stationary and moving interfaces. The application of the HF method to annular flow modelling has shown massive improvements in interface modelling and potential for research on annular flow hydrodynamics and heat transfer.

Using Extrapolation Techniques in VOF Methodology to Model Expanding Bubbles

A numerical method is presented to study cavity and bubble dynamics. The liquid phase is assumed to be inviscid and incompressible and separated from the gas or vacuum phase by a free surface. On the free surface the stress tensor reduces to a spatially constant pressure. The flow in the bulk of the liquid is computed using a second-order-in-time projection method. The interface is advected and reconstructed using a Volume-of-Fluid (VOF) method. Setting the pressure on the free surface to the prescribed value involves a modified stencil on nodes close to the interface. This modifed stencil contains interpolated pressures on branches that are cut by the interface. Capillary e_ects are taken into account by adding the Laplace law pressure increment to these prescribed pressures. The curvature that appears in the Laplace law is computed using the height-function method. The VOF advection and momentum advection schemes both require an extension of the velocity in a two-layer wide ghost cell region on the grid across the free surface. This ghost layer is computed in two stages. In the preliminary stage a first-order velocity extrapolation of the liquid velocity field to the ghost layers is performed. In the second stage the ghost layer velocities are projected on the space of divergence free velocities using an auxiliary projection step. The whole procedure is implemented in a free code developed with the help of Gretar Tryggvason and Yue (Stanley) Ling and is available at http://parissimulator.sf.net. Tests are perfomed on radial flows with spherical symmetry except for boundary conditions far from the bubble in a cubic box. In such a geometry the flow is predicted by solutions of the Rayleigh-Plesset equation. Good comparison to the Rayleigh-Plesset solution for a single bubble with low and moderate amplitude oscillations is shown. Perspectives for parallel simulations involving very large numbers of bubbles are given.

A mesh-decoupled height function method for computing interface curvature

In this paper, a mesh-decoupled height function method is proposed and tested. The method is based on computing height functions within columns that are not aligned with the underlying mesh and have variable dimensions. Because they are decoupled from the computational mesh, the columns can be aligned with the interface normal vector, which is found to improve the curvature calculation for under-resolved interfaces where the standard height function method often fails. A computational geometry toolbox is used to compute the heights in the complex geometry that is formed at the intersection of the computational mesh and the columns. The toolbox reduces the complexity of the problem to a series of straightforward geometric operations using simplices. The proposed scheme is shown to compute more accurate curvatures than the standard height function method on coarse meshes. A combined method that uses the standard height function where it is well defined and the proposed scheme in under-resolved regions is tested. This approach achieves accurate and robust curvatures for under-resolved interface features and second-order converging curvatures for well-resolved interfaces.

A high-precision calculation method for interface normal and curvature on an unstructured grid

In the volume-of-fluid algorithm, the calculations of the interface normal and curvature are crucially important for accurately simulating interfacial flows. However, few methods have been proposed for the high-precision interface calculation on an unstructured grid. In this paper, the authors develop a height function method that works appropriately on an unstructured grid. In the process, the definition of the height function is discussed, and the high-precision calculation method of the interface normal is developed to meet the necessary condition for a second-order method. This new method has highly reduced computational cost compared with a conventional high-precision method because the interface normal calculation is completed by solving relatively simple algebraic equations. The curvature calculation method is also discussed and the approximated quadric curve of an interface is employed to calculate the curvature. Following a basic verification, the developed height function method is shown to successfully provide superior calculation accuracy and highly reduced computational cost compared with conventional calculation methods in terms of the interface normal and curvature. In addition, the height function method succeeds in calculating accurately the slotted-disk revolution problem and the oscillating drop on unstructured grids. Therefore, the developed height function method is confirmed to be an efficient technique for the high-precision numerical simulation of interfacial flows on an unstructured grid.

A critical comparison of surface tension models for the volume of fluid method

In many different fields of research, the interactions between two immiscible fluids are of importance. To study these flows in industrial equipment, a multi-scale modeling approach is often used. In this approach, the smallest scale models apply detailed information in the form of closure equations for the larger scale models, which can model complete industrial equipment. This paper will focus on the improvement of the smallest scale model; direct numerical simulations employing the Volume of Fluid model. In this model, mass is inherently conserved because of the surface treatment, but this treatment also poses a challenge in calculating the surface properties like the surface tension. In this paper, three different surface tension models for the Volume of Fluid were tested: the generally used Continuum Surface Force (CSF) model, the height function model and the novel tensile force method. From the verification tests, it was concluded that both the height function model and the tensile force method are an improvement of the CSF model. The single bubble simulations showed that the height function method works best for small bubble (Eo<1). This is due to problems with connectivity for the tensile force method. While for the larger bubbles (Eo>10), the tensile force method is the best functioning surface tension model, because the calculation of the curvature in the height function method uses a stencil in which the distance between two interfaces in the direction of the normal should at least be four grid cells. In all the other tested cases, the height function model and the tensile force method perform equally well. The Morton number changes the ranges for the region of use of the surface tension models slightly when logMo≤−7 (the height function model can only be used when Eo≤2, while the tensile force method can only be used at Eo≥2) and logMo≥1 (the height function model can in this region also be used when Eo>10).

Numerical investigation of the influence of leading and sequential bubbles on slug flow boiling within a microchannel

Multiphase CFD simulations are presently employed to investigate the flow boiling of multiple sequential elongated bubbles in a horizontal microchannel. Most of the computational studies published so far explored the features of boiling flows within microchannels by simulating the fluid-dynamics of a single evaporating bubble, but the present work shows that multiple bubble simulations are necessary to capture the essential features of the heat transfer process of a slug flow. In particular, it is shown that leading and sequential bubbles interact thermally and hydrodynamically due to the evaporation process, thus possessing different growth rates, velocities and thicknesses of the thin liquid films trapped between the bubbles interfaces and the channel wall. The evaporation of this thin liquid film is the dominant heat transfer mechanism in the vapor bubble region and the transit of trailing bubbles strongly enhances the time-averaged heat transfer coefficient of the bubble-liquid slug unit, by as much as 60% higher relative to the leading bubble under the operating conditions presently set. Furthermore, the presence of a recirculating vortex just after the tail of the bubble in the liquid slug trapped between the bubbles was found in the simulations, significantly improving the heat transfer between the wall and the bulk liquid, thus maintaining the heat transfer coefficient much higher than otherwise expected in the liquid slug region as well. Finally, a new multiple bubble heat transfer model is proposed to predict the local variation of the heat transfer coefficient, which might prove to be useful to improve the current boiling heat transfer methods, such as the three-zone model of Thome et al. [1,2]. The numerical framework employed to perform this study was the commercial CFD solver ANSYS Fluent 12 with a Volume Of Fluid interface capturing method, which was improved here by implementing external functions, in particular a Height Function method to better estimate the surface tension force and an evaporation model to compute the phase change.

Numerical investigation of hydrodynamics and heat transfer of elongated bubbles during flow boiling in a microchannel

Flow boiling within microchannels has been explored intensively in the last decade due to their capability to remove high heat fluxes from microelectronic devices. However, the contribution of experiments to the understanding of the local features of the flow is still severely limited by the small scales involved. Instead, multiphase CFD simulations with appropriate modeling of interfacial effects overcome the current limitations in experimental techniques. Presently, numerical simulations of single elongated bubbles in flow boiling conditions within circular microchannels were performed. The numerical framework is the commercial CFD code ANSYS Fluent 12 with a Volume Of Fluid interface capturing method, which was improved here by implementing, as external functions, a Height Function method to better estimate the local capillary effects and an evaporation model to compute the local rates of mass and energy exchange at the interface. A detailed insight on bubble dynamics and local patterns enhancing the wall heat transfer is achievable utilizing this improved solver. The numerical results show that, under operating conditions typical for flow boiling experiments in microchannels, the bubble accelerates downstream following an exponential time-law, in good agreement with theoretical models. Thin-film evaporation is proved to be the dominant heat transfer mechanism in the liquid film region between the wall and the elongated bubble, while transient heat convection is found to strongly enhance the heat transfer performance in the bubble wake in the liquid slug between two bubbles. A transient-heat-conduction-based boiling heat transfer model for the liquid film region, which is an extension of a widely quoted mechanistic model, is proposed here. It provides estimations of the local heat transfer coefficient that are in excellent agreement with simulations and it might be included in next-generation predictive methods.

On the properties and limitations of the height function method in two-dimensional Cartesian geometry

In this study we define the continuous height function to investigate the approximation of an interface line and its geometrical properties with the height function method. We show that in each mixed cell the piecewise linear interface reconstruction and the approximation of the derivatives and curvature based on three consecutive height function values are second-order accurate. We also discuss the quadratic reconstruction and fourth-order accurate expressions of the normal and curvature. We present a hierarchical algorithm to compute the normal vector and curvature of an interface line with the height function method that switches automatically between second- and fourth-order approximations and that can be applied also when the local radius of curvature is of the order of the grid spacing.

Interface curvature via volume fractions, heights, and mean values on nonuniform rectangular grids

Estimating the local curvature of an interface involves the local determination of normals to the interface, and the rates that they turn along the interface. This is challenging in volume-of-fluid type methods since the interface between materials is specified by the relative amount it cuts off from the computational cells that it crosses (also referred to as volume fraction data) rather than by a discrete set of points lying on the interface itself. In this work, we generalize the height function method to nonuniform rectangular grids. We demonstrate analytically and numerically that—using three successive (adjacent) integral mean values (or “column” heights)—interface curvature can be estimated to second-order accuracy, the first derivative (interface normal) to third-order accuracy, and the curve location to fourth-order accuracy—each at its own special points. We also show that there are special points where the curvature can be estimated to fourth-order accuracy when using five successive mean values instead. Underlying all this is a result about the accuracy of the jth-derivative of the kth-degree polynomial that interpolates a function F at k + 1 stencil points placed irregularly in an interval of width h. Namely, for all smooth enough functions F and for k fixed and h getting small, there are k + 1 - j special points in the interval at which the error in the jth derivative is of order O ( h k + 2 - j ) . (This is one order higher than the usual O ( h k + 1 - j ) error holding over the whole interval for that derivative.) The special points in the interval are the k + 1 - j (real) zeroes of certain F-independent polynomials, of degree k + 1 - j , with coefficients depending on the interval’s stencil points. In our case, let F ( x ) be an indefinite integral of the unknown interfacial curve f ( x ) . Then F ( x i ) at k + 1 successive stencil points x i is calculated using cumulative sums of the k successive integral mean values of f, weighted by the successive interval sizes. The k + 1 points ( x i , F ( x i ) ) are now interpolated by a kth degree polynomial ( PF ) ( x ) . It’s first derivative ( PF ) ( 1 ) approximates the unknown curve f; while ( PF ) ( 2 ) and ( PF ) ( 3 ) , respectively, approximate f ( 1 ) and f ( 2 ) . Thus, the results above using three successive mean values correspond to k = 3 and j = 3 , 2 , 1 . The curvature result using five successive mean values correspond to k = 5 and j = 3 . For this last case, since the curvature κ = f ( 2 ) / ( 1 + ( f ( 1 ) ) 2 ) 3 / 2 , we use the facts that f ( 1 ) = F ( 2 ) = ( PF ) ( 2 ) + O ( h 4 ) on the whole stencil interval, not just at the three special points having O ( h 4 ) accuracy for f ( 2 ) .

Development of a conservative 3D free surface code

The paper describes the salient features of a three-dimensional finite difference code at present under development for Large Eddy Simulation of free-surface turbulence. It is essential in this application that the numerical approximations conserve mass, momentum, and kinetic + potential energy while allowing for large movements of the free surface and these are presented in detail for cells cut by the surface. The surface advancement is by a combined Volume-of-Fluid and height function method which conserves mass exactly over finite time steps, but does not permit breaking waves. Momentum is advanced by control volume techniques, and the potential singularity caused by very small surface cells is avoided by a new split-merge technique which is numerically stable. Results are presented for a two-dimensional inviscid large amplitude wave in a periodic box.