Experimental validation of numerical simulations of free-surface flow within casting mould cavities

Abstract

Purpose

The paper targets on providing new experimental data for validation of the well-established mathematical models within the framework of the lattice Boltzmann method (LBM), which are applied to problems of casting processes in complex mould cavities.

Design/methodology/approach

An experimental campaign aiming at the free-surface flow within a system of narrow channels is designed and executed under well-controlled laboratory conditions. An in-house lattice Boltzmann solver is implemented. Its algorithm is described in detail and its performance is tested thoroughly using both the newly recorded experimental data and well-known analytical benchmark tests.

Findings

The benchmark tests prove the ability of the implemented algorithm to provide a reliable solution when the surface tension effects become dominant. The convergence of the implemented method is assessed. The two new experimentally studied problems are resolved well by simulations using a coarse computational grid.

Originality/value

A detailed set of original experimental data for validation of computational schemes for simulations of free-surface gravity-driven flow within a system of narrow channels is presented.

Citation

Bublík, O., Lobovský, L., Heidler, V., Mandys, T. and Vimmr, J. (2021), "Experimental validation of numerical simulations of free-surface flow within casting mould cavities", Engineering Computations, Vol. 38 No. 10, pp. 4024-4046. https://doi.org/10.1108/EC-08-2020-0458

1. Introduction

Both experimental methods and computational simulations are used in the pre-production stage of many industrial applications, which deal with free-surface flow problems, e.g. liquid transport, sloshing, mixing of fluids, spray formation, etc. This study is a part of a project on material casting and rapid prototyping of large thin structures of complex three-dimensional geometries. In order to computationally resolve the flow within a thin casting mould cavity, a fine computational grid with a large number of grid elements is often required. Thus, the numerical simulation becomes expensive in terms of computational time and memory demands.

Aiming at minimisation of the computational time costs, the lattice Boltzmann method (LBM) is preferred within this study, as it is well suited for massive parallelisation on various architectures while its algorithm is relatively simple and efficient, Körner et al. (2005), Sano et al. (2008), Schornbaum and Rüde (2016), Calore et al. (2016). The LBM utilises a mesoscopic approach to fluid flow modelling as it models the macroscopic variables (e.g. pressure, velocity) by evaluating interactions of virtual fluid particles (expressed by means of particle distribution function), Higuera et al. (1989), Benzi et al. (1992). In the macroscopic limit, the LBM equations yield the Navier–Stokes equations, Chen and Doolen (1998), Wolf-Gladrow (2005).

In order to model the incompressible fluid flow, the LBM method does not require a solution of the Poisson's equation. This decreases its computational demands. Although two-phase flow LBM models of the casting processes exist, Szucki et al. (2017), Kharmiani and Passandideh-Fard (2018), a single phase model with free-surface boundary, Ginzburg and Steiner (2003), Zhang et al. (2019), is considered within this work in order to keep the computational simplicity and cost efficiency.

Similarly to other Eulerian grid-based computational methods, the LBM requires an algorithm for free-surface tracking in order to resolve the free-surface flow. This may be based on the marker-and-cell method, Harlow and Welch (1965), McKee et al. (2008), or the volume of fluid (VOF) method, Hirt and Nichols (1981), Torrey et al. (1987). Due to the mesoscopic description of continuum, the LBM includes the mass transfer calculation, which can be well utilised within the VOF scheme without a need of solving another equation, Thürey (2007), Janßen et al. (2013).

The in-house-implemented LBM solver also focuses on the surface tension and wetting effects at the liquid–solid interface, Bogner et al. (2016). Details on the applied algorithm are described hereafter.

Several benchmark tests with a well-known analytical solution are employed in order to assess the performance of the implemented algorithm. The key part of the paper provides a newly designed set of experiments, which target on free-surface flow propagation within a casting mould. The experiments are analogue to gravity casting of liquified material into a rigid mould with absence of thermodynamic or chemical effects, DeGarmo et al. (2003), Cleary et al. (2014). In order to enhance the repeatibility of the tests and to avoid extensive statistical analysis, the experiments target on laminar flow of Newtonian fluids. In addition, physical properties of the applied fluid samples (viscosity, surface tension and density) are identified during each experimental run.

While some published casting flow benchmarks address flow within large three-dimensional moulds of considerable bulk, e.g. Hirt and Harper (1988), or within large planar moulds, e.g. Sirrell et al. (1996), the presented experiments focus on propagation of casting flow in long narrow channels. Prediction of such flow is essential for applications involving complex three-dimensional channel systems, e.g. multiple-gate runner designs or complex mould structures as in Mun et al. (2015). An attention is paid to channels which length is excessively large in comparison to their cross section as well as to channels which cross section varies along the channel length.

Since the experiments are carried out under well-controlled laboratory conditions, it is aimed at providing a complete set of original experimental data (including rheological measurements) for assessment of performance and accuracy of numerical algorithms for casting flow simulations. The presented experimental results are available in an online repository “https://www.kme.zcu.cz/hydrolab/gravitycasting-newtonianfluid”.

2. Lattice Boltzmann method

The LBM evolved from lattice gas cellular automata (LGCA) of late 1980s, Frisch et al. (1986). First, the Boltzmann equation was used to calculate viscosity and transport coefficients within the LGCA frame, Wolfram (1986), Frisch et al. (1987). Consequently, the LBM was founded as a stand-alone computational scheme that uses continuous distribution functions rather than Boolean fields of LGCA, McNamara and Zanetti (1988). As a result, the LBM lacks the statistical noise of the LGCA method. The original LBM scheme was simplified by introducing the linearised collision operator, Higuera and Jiménez (1989), and further by application of the single time relaxation, Koelman (1991), Chen et al. (1992), which is based on Bhatnagar–Gross–Krook approximation, Bhatnagar et al. (1954). Stability of the method was further improved by the multiple relaxation time (MRT) scheme, d'Humières et al. (2002). An extensive review on the LBM and its substantial development over the last three decades can be found, e.g. in Wolf-Gladrow (2005), Perumal and Dass (2015). Stability of the method is addressed, e.g. in Sterling and Chen (1996), Worthing et al. (1997).

2.1 Discretisation

The LBM describes the continuum dynamics by simulating an interaction of virtual particles within an equidistant lattice grid, Higuera et al. (1989), Chen et al. (1992). These virtual particles statistically represent a swarm of real material particles in the phase space.

In the phase space, evolution of an arbitrary distribution function f=fx,u,t in time t can be described by the Boltzmann equation

(1)∂f∂t+u⋅∂f∂x=Ωf

where x=x1,x2,x3T is the position vector and u=[u1,u2,u3]T is the velocity vector. Assuming that the distribution function f = f(x, u, t) represents the statistical density of particles with a given velocity u at a given point x in space, the term Ω(f) on the right hand side of equation (1) represents the collisions between the particles.

In order to discretise equation (1), the computational domain is covered by an equidistant grid of points x_i and a finite set of microscopic velocity direction vectors e_α is introduced. The velocity vectors e_α are chosen such that the convective scaling Δx = e_αΔt is satisfied for all α. The D3Q19 microscopic velocity model is adopted within this study, i.e. the velocity vectors e_α are defined in 19 directions within a three-dimensional lattice grid, see Figure 1 (left).

First, equation (1) is discretised in the velocity space, which yields a system of equations for all velocity vectors e_α.

(2)∂fα∂t+eα⋅∂fα∂x=Ω(fα)

where fα is a function of the spatial variable x and time t only, i.e. fα = fα(x, t). In the second step, equation (2) is discretised in space and time as

(3)fα(xi+eαΔt,t+Δt)=fα(xi,t)−M−1S(τ)M(fα−fαeq)

where the collision operator Ω(f_α) on the right hand side is approximated by multiple relaxation time operator (MRT), d'Humières et al. (2002). The transformation matrix M represents a linear mapping between the discrete velocity space f and the discrete moment space m=[m0,m1,m2,…,m18]T. The collision matrix S is a diagonal matrix of relaxation rates, which are either constant or depend on the relaxation time τ. The structure of the matrices M, S and m and the applied D3Q19 microscopic velocity model are described in detail in d'Humières et al. (2002).

The MRT collision operator improves the stability of the method for the relaxation time τ → 0.5⁺ and for large values of τ > 1. The extreme values of the relaxation time τ may be problematic when a single relaxation time, Bhatnagar et al. (1954), or a dual relaxation time, Ginzburg et al. (2008), collision operator is applied.

The function fαeq represents an equilibrium function, which is derived from the Maxwell distribution for the given velocity model, Succi (2001). In the case of fluid dynamics, the equilibrium function has the form

(4)fαeq(Ρ,u)=wαΡ1+3(eα⋅u)+92(eα⋅u)2−32(u⋅u)

where w_α are the velocity vector weights, which values are given according to the applied microscopic velocity model, Succi (2001), d'Humières et al. (2002), and Ρ is the density. Using the Chapman–Enskog expansion, the incompressible Navier–Stokes equations can be recovered, Chen et al. (1992).

For the calculation of the new time level according to the LBM scheme (3), it is advantageous to split the calculation into the collision and streaming steps. In the collision step, values of the post-collision distribution functions fαc are computed as follows:

(5)fαc(xi,t)=fα(xi,t)−M−1S(τ)M(fα(x,t)−fαeq).

In the propagation step, the precomputed post-collision distribution functions fαc are propagated into the neighbouring cells in the directions of vectors e_α

(6)fα(xi+eαΔt,t+Δt)=fαc(xi,t).

The collision step as well as the propagation step are evaluated at the given point x_i independently of its neighbours, therefore both of these steps are ideal for parallel processing.

The external body force can be introduced into the collision step as follows:

(7)fαc(xi,t)=fα(xi,t)−M−1S(τ)M(fα(xi,t)−fαeq)+fαg(xi,t)

where fαg(xi,t)=3wαΡ(eα⋅g) and g is the body force acceleration.

In order to enforce the no-slip boundary conditions at rigid wall boundaries, the bounce back scheme of Ladd (1994) is utilised.

2.2 Free-surface modelling

An interface between two immiscible fluids can be treated by a multi-phase flow model, Gunstensen et al. (1991), Swift et al. (1996), Sbragaglia et al. (2007), or in case when a ratio between fluid densities is large (e.g. between liquid and gas), the influence of the low density fluid may be neglected and a free-surface flow model applied Ginzburg and Steiner (2003). Although there are multi-phase flow algorithms capable of resolving high density ratio interfaces, e.g. Bao and Schaefer (2013), significant unphysical fluctuations, so called spurious currents Connington and Lee (2012), may occur at liquid–gas interfaces (e.g. density ratio 1000:1). These currents may deteriorate the resulting velocity field and cause instability of the calculation.

In order to keep the LBM algorithm simple and computational demands low, the free-surface flow algorithm is employed below to track the liquid–gas interface. The implemented LBM free-surface flow model is based on the VOF approach, Thürey (2007), while the free-surface is tracked directly using the distribution functions f_α (which correspond to the mass transfer in the direction of the microscopic velocity vectors e_α).

Three types of computational grid cells are distinguished – liquid, gas and interface cells, see Figure 1 (right). The liquid cell is fully filled with liquid and the standard LBM algorithm apply. The gas cell is void and discarded from calculations. The interface cell is partially filled with the liquid, and the computational algorithm is modified. In order to track the free surface, the mass function m and the fluid fraction χ = m/Ρ are introduced. The value of function χ is equal to 0 for gas, 1 for liquid and χ ∈ (0, 1) for the liquid–gas interface. In order to capture the movement of the free surface, the mass flux between the interface cell and its neighbouring cells is computed as

(8)Δmα(x,t+Δt)=(fα¯(x+eαΔt,t)−fα(x,t))χ(x+eαΔt,t)+χ(x,t)2

where α¯ denotes the opposite direction to α. The total mass change at the given cell is

(9)m(x,t+Δt)=m(x,t)+∑α=0nαΔmα(x,t+Δt).

where n_α = 18 for the D3Q19 model.

When the total mass change is computed, the collision and propagation steps are performed and the macroscopic variables are evaluated. During the mass change, the interface cells can be filled up or emptied and the cell type has to be changed according to the following prescription:

(10)m(x,t+Δt)≥Ρ(x,t+Δt)+ϵ→liquidcell,m(x,t+Δt)≤−ϵ→gascell,

where ϵ is a suitable small constant, which helps to prevent spurious oscillations between the cell types. Its value needs to be tuned for each specific case. Based on the numerical testing, the value of ϵ = 0.05 has been determined for all applications presented below.

When the interface cell changes its state and it gets filled up (emptied respectively), the new fluid cell (the new gas cell respectively) gets right next to at least one gas cell (fluid cell respectively). The gas cells which are adjacent to the filled up cells are changed to interface cells. The same applies to the fluid cells adjacent to the emptied cells. In the end of the step, there are no fluid and gas cells adjacent to each other.

To preserve the local and the global mass conservation during the cell conversion, the excessive or missing mass of liquid within the interface cell

(11)merr(x,t+Δt)=m(x,t+Δt)−1fortransformationtoliquid,m(x,t+Δt)fortransformationtogas,

must be redistributed among the neighbouring interface cells x_b so that

(12)m(xb,t+Δt)=m(xb,t+Δt)+merr(x,t+Δt)Nb

where N_b is a number of the neighbouring interface cells.

The unknown distribution functions fα¯(x,t+Δt) coming from the gas cell to the interface cell can be reconstructed using the equilibrium function fαeq(Ρe,u) as follows:

(13)fα¯(x,t+Δt)=fαeq(Ρe,u)+fα¯eq(Ρe,u)−fα(x,t)

where u is the interface velocity, and the density Ρ_e depends on the external pressure p_e and on the Laplace pressure caused by the surface tension, i.e. Ρ_e = 3(p_e + 2κσ) where σ is the surface tension, and κ denotes the free-surface curvature.

2.2.1 Normals and curvature computation

For computation of normal vectors and curvature of the free surface (the liquid–gas interface), an approximation of the spatial derivatives needs to be introduced, Bogner et al. (2016). When an arbitrary function ϕ is defined on the grid of cells, then

(14)∇ϕ(x)=1∑α=0nαwα∑α=0nαϕ(x+eα)eα.

When the value of function ϕ in the direction e_α would have to be evaluated outside the computational domain, the following linear interpolation is applied:

(15)ϕαoutside=2ϕ(x)−ϕ(x+eα).

Formula (14) can be used to calculate the gradient of the volume fraction χ. The resulting vector ∇χ is parallel to the normal vector at the liquid–gas interface. Thus, the normal vector is computed by normalising the vector ∇χ to unit length.

(16)n=∇χ∣∇χ∣.

A value of the normal vector n at an interface cell adjacent to the solid wall is corrected according to the contact angle θ at the liquid–gas–solid interface. The corrected normal vector is computed using the following relation:

(17)ncorrect=nw cos(θ)+tw sin(θ)

where n_w is the normal vector of the solid wall. The tangent vector t_w of the wall is defined as

(18)tw=n−nwn⋅nw.

In order to precisely resolve the liquid–gas–solid interface, the dynamic contact angle (which reflects the motion of the interface along the solid surface) should be considered, Yokoi et al. (2009). However, for simplicity the static contact angle θ (which is a constant determined for non-moving interface) is employed.

The normal vector n is used to compute the free-surface curvature κ as follows:

(19)κ=∇⋅n

where the divergence of the normal vector field is calculated similarly to (14). In the vicinity of the computational domain boundary, the normal vector n is corrected in the sense of equation (15). The free-surface curvature is calculated at the interface cells only. However, in order to calculate the curvature, the normal vectors must be defined at all cells in the interface cell neighbourhood.

In order to minimise the magnitude of spurious fluctuations resulting from errors in the finite difference approximation of the divergence of the normal vector field, the gradient of the volume fraction respectively, the normal vectors are evaluated based on the smoothed volume fractions χ^ which are calculated as

(20)χ^=1∑α=0nαwα∑α=0nαwαχ(x+eα).

when the value of χ in the direction e_α is outside the domain, and the linear interpolation (15) is utilised. In order to further smooth the volume fraction χ^, the calculation of (20) can be repeated.

2.3 Units conversion

The LBM scheme uses different length and time scales than physical dimensions. For the correct setting of computational parameters, e.g. the relaxation time, and for obtaining physical velocity and pressure field, the unit conversion is necessary. Let index LB denote the variable in lattice Boltzmann units and index phys address the variable in physical units. The unit conversion of arbitrary variable φ from LBM units to physical units is defined using the converter C_φ so that

(21)φphys=CφφLB.

In total, two elementary converters, the space converter C_H and the velocity converter C_U, are defined. If H is an arbitrary physical characteristic dimension of the computational domain and n is the number of corresponding lattice cells related to H, then the space converter can be defined as

(22)CH=Hn.

Similarly, when umaxphys is the expected maximal value of the physical velocity and umaxLB is the corresponding maximal velocity in LBM units, the velocity converter is defined as

(23)CU=umaxphysumaxLB.

For stability reasons, the value of the maximal velocity in LBM units is set to umaxLB=0.1. In addition, the density converter can be considered in the following form:

(24)CΡ=ΡphysΡLB.

Since the density in LBM units Ρ_LB = 1, the density converter C_Ρ is equal to the physical density Ρ_phys. Converters of other physical variables can be derived using the elementary converters C_H and C_U.

(25)CT=CHCU,CP=13CU2,Cν=CH2CT,CG=CHCT2,CST=CH3CT2

where C_T is the time converter, C_P is the pressure converter, C_ν is the kinematic viscosity converter, C_G is the gravity acceleration converter and C_ST is the surface tension converter.

The relaxation time τ in LBM simulations is related to the physical kinematic viscosity ν_phys = C_νν_LB by the relation

(26)τ=12+3νLB.

The physical value of the time step can be expressed as Δt_phys = C_TΔt_LB. The value of Δt_LB = 1 is a common choice.

Using the distributions functions f_α, the macroscopic flow quantities in LBM units could be expressed as

(27)ΡLB=∑α=0nαfα,ΡLBuLB=∑α=0nαfαeα,pLB=13ΡLB.

2.3.1 Algorithm

The LBM scheme with the free-surface model includes several computational steps. To produce efficient and valid LBM code, the following sequence of computational steps is applied:

1. Compute the mass exchange between cells and find filled up/emptied cells (8).

2. Compute the smoothed function χ^ according to formula (20).

3. Compute the normal vectors and set the wall normals according (16) and (17).

4. Compute the curvature according to equation (19).

5. Perform the propagation step (6).

6. Compute the macroscopic variables (27) and perform the collision step (5).

7. Rewrite a flag for filled/emptied cells and mark new interface cells (10).

8. Redistribute mass from the cell conversion, see equation (12).

3. Benchmarks with theoretical solution

The implemented LBM code is first tested using four benchmarks, which address its capability to resolve viscous fluid flow and surface tension effects.

3.1 Poiseuille flow

The Poiseuille flow benchmark focuses on a laminar viscous flow in a fully flooded channel between two parallel infinite plates, Sutera and Skalak (1993). The plates are parallel with the xy plane, and the flow is driven in x direction by a body force g, Figure 2 (left). An incompressible viscous fluid with the kinematic viscosity ν and density Ρ is considered. The fluid velocity in x direction can be determined analytically as

(28)u(z)=Ρg2νz(h−z)

where h is the channel breadth. The velocity u is a function of z coordinate only, and the maximal velocity umax=Ρgh28ν is obtained in the centre of the channel, i.e. for z=h2.

In order to compare the analytical solution with the LBM results, the LBM simulation is performed on a cubic computational domain of h = 0.1 m, while periodic boundary conditions are applied in x and y direction. The numerical results are computed for the flow driven by the body force g = 1 m s⁻², the fluid density Ρ = 1 kg m⁻³ and the visocity ν = 0.001 m² s⁻¹. The theoretically predicted Reynolds number Re=umaxhν is equal to 125. The computed velocity profile for 40 grid cells across the channel breadth is provided in Figure 2 (right).

In order to assess the accuracy of the LBM solution for various number of lattice cells across the channel breadth n, the relative error is evaluated as

(29)errn=|umaxLB−umax|umax

where umaxLB is the maximal velocity computed by the LBM scheme. As displayed in Table 1, a reasonable accuracy is achieved already for a coarse lattice grid. With only five lattice cells across the channel breadth, the relative error is about 4%. The second-order convergence of the LBM results to the analytical solution of the Poiseuille flow is confirmed for large values of n, Figure 3.

3.2 Advection

While the Poiseuille flow test does not involve free-surface interfaces, the second benchmark targets on convergence of the LBM free-surface flow algorithm. A one-dimensional advection of a bulk of fluid is explored below. The problem is solved on a three-dimensional computational domain with absence of rigid boundaries and surface tension effects.

The cuboidal computational domain of dimensions H × H × 2H, where H = 0.01 m, is employed. Periodic boundary conditions are applied in x and y direction. Initially, a part of the computational domain (3/16 of its volume) is filled with liquid of density Ρ = 1000 kg m⁻³ and kinematic viscosity ν = 10^–4 m² s⁻¹, while the rest of the domain is void, Figure 4 (left).

The initial velocity u = 0.1 m s⁻¹ is imposed on a square block of fluid in z direction. Since no body force field (such as gravity) is considered, the analytic solution for displacement of the bulk of fluid at time t = 0.1 s is L = 0.01 m, Figure 4 (right). The LBM solution is computed for several lattice grids with different spatial resolution expressed by means of the number of lattice cells n across the channel breadth H.

The results presented in Table 2 indicate that the solution of the advection problem is first-order convergent. This is due to employment of the free-surface flow model, which utilises the first-order accurate discretisation of the VOF scheme, equation (13).

3.3 Laplace pressure

The Laplace pressure denotes the difference in pressure across a curved interface between static fluids (e.g. water droplet and air) caused by the surface tension, Butt et al. (2003). The Laplace pressure may be predicted theoretically by using the Young–Laplace equation

(30)Δp=σ1R1+1R2

where σ is the surface tension at the fluids interface and R₁, R₂ are the principal radii of the interface curvature (for a spherical droplet R₁ = R₂).

In order to test the implemented LBM surface tension model, a spherical droplet formation from an initially cubic bulk of liquid (with absence of gravity) is simulated, Figure 5. The calculations are performed for liquids of the density Ρ = 1,000 kg m⁻³, the kinematic visocity ν = 0.001 m² s⁻¹ and several values of the surface tension σ = {0.02, 0.04, 0.08} N m⁻¹.

Within the LBM simulation, a cubic computational domain with sidelength of 0.1 m is discretised by 40 × 40 × 40 grid cells. Initially, there is a cube of 30 × 30 × 30 liquid cells in the centre of the domain, which is surrounded by void cells representing the air. During the simulation, a spherical liquid droplet is created solely due to liquid's surface tension. Since the liquid viscosity is relatively high, any considerable droplet shape oscillations are avoided.

When the fluid motion vanishes and the spherical droplet remains still, the resulting Laplace pressure Δp_LB computed by LBM is compared to its theoretically predicted value Δp. The theoretical and computational results are summarised in Table 3. In all cases, the relative error

(31)errσ=|ΔpLB−Δp|Δp

is about 2% and the computed droplet radius is approximately 0.0465 m, which corresponds to theoretically predicted radius 3434π1/3L.

Convergence of the LBM solution of spherical droplet radius to its theoretically predicted value is studied for the surface tension σ = 0.02 N m⁻¹ and provided in Table 4.

3.4 Capillary rise

The third benchmark targets on the LBM capability of resolving the capillary effects caused by the surface tension σ and by the adhesion between the liquid and capillary walls (which leads to wetting the walls under the contact angle θ), Figure 6. When the distance between the capillary walls is sufficiently small, the air–liquid interface curvature is approximately uniform and a significant Laplace pressure is detected, Batchelor (2000). The theoretically predicted height of the capillary rise H is given by balancing the Laplace pressure (30) and the hydrostatic pressure exerted by the column of liquid

(32)H=2σ cosθΡgw

where w is the capillary breadth and g is the gravitational acceleration. Equation (32) is valid for very narrow capillaries only, where a uniform curvature of the air–liquid interface Rθ=wcosθ is expected. In larger capillaries, a correction term must be included, Bullard and Garboczi (2009),

(33)Hc=2σ cosθΡgw−wf(θ)g(γ)e−4.48γ1/8

where f(θ) = 9.24 cos θ + 2.13 cos³θ, g(γ)=0.834γ−0.024γ and γ=Ρgw2σ.

The LBM solver is tested for a capillary rise in a narrow vertical channel between two parallel infinitely wide plates partially immersed in an infinitely large pool of liquid. The test fluid has the density Ρ = 1,000 kg m⁻³, the surface tension σ = 1 N m⁻¹ and the dynamic viscosity η = 0.0002 Pa s. The infinitely large pool is represented by a uniformly discretised computational domain with periodic boundaries. The grid cell size is dx = 0.0005 m.

In total, two sets of numerical tests are carried out. First, an accuracy of LBM results for the capillary breadth w = 0.0095 m and three different values of contact angle θ = {60°, 45°, 30°} is examined. The LBM computed capillary rise H_LB is compared to analytically predicted values in Table 5). The deviation of H_LB from both uncorrected and corrected theoretically predicted values H and H_c is evaluated as

(34)errH=|HLB−H|H,

(35)errHc=|HLB−Hc|Hc.

Since the capillary breadth w is relatively large, the deviation from the corrected capillary rise errHc should be preferred. Nevertheless, both err_H and errHc grow with decreasing value of the contact angle. This can be attributed to the free-surface tracking described in section 2.2, which assumes that each interface cell has at least one liquid and one gas cell in its neighbourhood defined according to D3Q19 velocity model. This assumption is met for contact angle θ ≥ 45° (when the value of errHc is close to 2% or smaller). However, the accuracy of simulations for low contact angle values, such as θ = 30°, is limited.

The second set of LBM tests is performed for the contact angle θ = 45° and three values of the capillary breadth w = {0.005, 0.0095, 0.0145} m. A qualitative comparison of the LBM results is provided in Figure 7, which shows a cropped image of the same part of all capillaries next to each other after the LBM computation has converged. The liquid cells are marked in dark blue, the interface cells are cyan and the void (gas) cells are gray. The quantitative results are presented in Table 6 together with the corresponding analytical data.

As expected, when assuming the uniform curvature of the free surface across the capillary breadth, the difference between the LBM results and the theoretically predicted values by equation (32), expressed in terms of err_H, grows with the increasing capillary breadth w. When the corrected relation (33) is considered, the relative error errHc almost vanishes for w = 0.005 m and for w = 0.0145 m. However, the resulting errHc does not show any general tendency with regards to value of w.

4. Experimental measurements

The study targets on applications in casting industry, problems of free-surface flow within a system of narrow rigid channels and thin cavities in particular. An experimental campaign on gravity driven casting flow is designed in laboratory conditions, which enable a straightforward control of each experimental run. A total of two experimental setups are utilised.

Each test geometry is divided into two parts, a reservoir and a test section, by a removable dam. Before each experimental run, the reservoir holds the prescribed amount of liquid while the test section is empty and dry. When the gate is removed, the liquid flows into the test section driven solely by gravity. The test geometries are made of transparent polymethyl methacrylate, which allows for optical measurements of the fluid flow by a system of digital cameras with native resolution of 4 Mpx @ 30 Hz and 2 Mpx @ 60 Hz.

Chemically stable moderately viscous aqueous solutions are employed such that the resulting casting flow is laminar. The dynamic viscosity η, the surface tension σ and the density Ρ of these Newtonian test fluids are summarised in Table 7. The presented data are based on measurements of fresh fluid samples before every experimental run as well as of already used samples during the casting experiments, Lobovský and Mandys (2019).

4.1 Horizontal channel

Within the first setup, the test section is made of a single horizontal channel of constant rectangular cross section, which does not branch, Figure 8. It has a single inlet and a single outlet. The channel winds in a spiral so that its length is large in comparison to its cross section. This enables to study evolution of the flow while the ratio between the force driving the flow and the resistance of the system decreases significantly as the fluid front propagates along the channel and the liquid level in the reservoir drops. As a result, any inaccuracies get pronounced along the channel length.

The horizontal channel is 3 mm high, 10 mm wide and 1,160 mm long. At its inlet, the channel is connected to a vertical reservoir, which is initially filled up with 90 ml of liquid held by a vertical gate. Details on the test geometry dimensions are displayed in Figure 8. The small height of the channel helps to minimise the gate removal time, which is less than 1/60 s. This estimate is limited by a time resolution of the applied camera system.

Once the gate is removed, the liquid flows rapidly into the test section under gravity. The fluid motion is recorded by the cameras at 20 Hz framerate. The fluid front propagation along the channel centreline (blue line in Figure 9 (left)) is analysed, and a comparison of ten experimental runs with LBM simulations is provided in Figure 9 (right).

The quantitative analysis of the fluid front motion within the channel is given in Table 8, which presents the arrival times of the fluid front at a location specified by the distance from the gate. The average duration of the experimental run, i.e. the time until the entire channel gets flooded, is about 40.38 s with the standard deviation of about 1.11 s. The minimal and maximal duration of the test run is 38.20 and 42.05 s, respectively. Differences between the tests may be attributed to manual removal of the gate, which causes variations in the gate removal times. These variations cannot be quantified due to the temporal resolution of the applied data acquisition system.

In order to address the LBM ability to quickly estimate the studied flow, LBM simulation at coarse resolution is executed. A uniform computational grid of 242,000 cells provides only three grid cells across the channel height. In the initial stage of the simulated flow, the LBM data correspond well to the fastest experimental runs, Figure 9 (right). This is caused by a computational representation of the gate removal, which is instantaneous. However at later stages of the LBM simulation, the fluid front propagation slows down and the complete channel flood time agrees well with the longest experimental runs. This may be observed also in Figure 10, which presents the snapshots of the flow propagation within the horizontal channel as recorded from a randomly selected experimental run and the LBM simulation.

The slower propagation of the LBM fluid front is affected by the applied computational grid resolution. As discussed in section 3.1, for a low number of grid cells across the channel section, the computed flow velocity gets underestimated by several percents. In addition, the flow is also influenced by the viscosity model. Although the Newtonian viscosity is applied within the LBM model, the dynamic viscosity of the experimentally examined fluid is not constant, but slightly decreases with decreasing shear rate, Lobovský and Mandys (2019). Thus, the computationally predicted flux may be slightly underestimated. Despite these limits, the computational model performs well and is capable of resolving the general trends of the experimentally measured flow.

4.2 Vertical labyrinth

The second experimental setup provides data on the fluid flow within a vertically oriented test section, Figure 11. This section consists of three straight vertical channels with constant width and a trident of vertical channels, in which width varies along their length. All channels have a rectangular cross section and a constant breadth 10 mm. There is a free outlet for air at the end of each channel. Similarly to the previous setup, the vertical gate obstructs a 3 mm high horizontal channel, which is a single inlet into the test section. The vertical reservoir is initially filled-up with 72 ml of the test fluid.

The fluid flow is driven solely by gravity. Its evolution at the very beginning of a randomly selected test is shown by four subsequent video frames recorded at 60 Hz in Figure 12. The gate removal time is less than 1/60 s. Propagation of the fluid front within the vertical labyrinth is then recorded at framerate of 30 Hz and analysed along the six pathlines depicted as coloured lines in Figure 13.

Results of ten experimental test runs during the first 6 s of each run are presented in Figure 14. The distinctive change in a slope of the plotted curves corresponds to locations where the pathline changes its direction from horizontal to vertical or vice versa. Similarly to the horizontal setup, variations in between experimental runs (which get pronounced especially along the longer pathlines) are mainly related to nuances in the manual gate removal. The position P_color is defined as a distance from the gate along the pathline of given color. Averaged values of P_color at selected time instants (based on data of all test runs) are given in Table 9.

The LBM simulations are performed for a uniform lattice grid of 600,000 cells with six cells across the height of the horizontal inlet into the labyrinth. As displayed in Figure 14, the numerical and the experimental results are in a good agreement, though there is an evident deviation of the results when the fluid front reaches the base of the channel trident and later along the horizontal channel parts. This is obvious also from qualitative comparison in Figure 15. The delay of the LBM results may be well observed along the magenta pathline where it gets pronounced with increasing time when the force driving the flow decreases and the viscous forces become significant as the test section gets flooded. This difference in the LBM results may be attributed to the spatial resolution of the computational grid as discussed in the section above.

5. Conclusion

The performed experimental study provides original data on problems of gravity-driven free-surface flow analogous to industrial casting processes. Simplified problems are presented and allow for a good control of execution of each experimental run, a repeatibility of the tests, a straightforward analysis of the measured data and relatively simple description of the experimental problem. The newly designed experimental test setups are suitable for benchmarking the computational solvers for problems of free-surface viscous flow within a system of narrow long channels. The experiments target on complex three-dimensional channel structures, which are present in multiple-gate runner designs and geometrically complex casting moulds. All casting flow experimental data are supported by the rheological measurements.

Based on the recorded experimental data, it is shown that the implemented LBM code is capable of resolving the flow propagation with reasonable accuracy for relatively coarse computational grids. This decreases the computational demands of the numerical simulations and makes the presented LBM approach attractive for real-world applications. Results of the LBM simulations for the adopted benchmark tests agree well with their theoretical solution. The method shows the second-order convergence for viscous flows with absence of free-surface boundaries. When the free-surface boundaries apply, the LBM scheme becomes first-order convergent. This is caused by the first-order accuracy of the discretisation applied within the VOF scheme, which is used to track the liquid–gas interface.

Since the LBM algorithm does not require solution of the non-linear Poisson equation for the pressure field in the incompressible fluid, the LBM solution is competitive in terms of computational efficiency as long as transient flow problems are of interest (an explicit time integration scheme is applied). The performance of the presented LBM scheme may be further improved by application of adaptive meshing techniques, which are convenient for discretisation of geometrically complex computational domains with large variations in channel cross section. Nevertheless, the presented LBM algorithm may be efficiently applied in a variety of manufacturing processes.

All presented calculations are performed on a single CPU using up to four cores.

Figures

Figure 1

Lattice Boltzmann method D3Q19 velocity model (left), LBM grid at liquid-gas interface (right)

Figure 2

Poiseuille flow: channel geometry and contour plot of LBM results for velocity u (left), LBM vs analytical solution of velocity profile along channel cross section (right)

Figure 3

Poiseuille flow: convergence of LBM results to analytical solution

Figure 4

Advection: position of the bulk of fluid (blue cells) within the computational domain at time t = 0 and 0.1 s

Figure 5

Laplace pressure: LBM simulation of liquid droplet formation, liquid cells (blue) at 0, 5,000 and 20,000 iterations

Figure 6

Capillary rise: scheme of liquid level rise in between capillary walls

Figure 7

Capillary rise: LBM results for contact angle 45°, capillary breadth w equals to 0.005, 0.0095, 0.0145 m, respectively

Figure 8

Horizontal channel: test geometry with dimensions in [mm]

Figure 9

Horizontal channel: channel centreline (top), experimental and LBM results for fluid front propagation (bottom)

Figure 10

Horizontal channel: comparison of randomly selected experimental run (top) and LBM data (bottom) at time 5, 15, 25 and 35 s

Figure 11

Vertical labyrinth: test geometry with horizontally oriented inlet, dimensions in [mm]

Figure 12

Vertical labyrinth: gate removal and free-surface flow at time −1/60, 0/60, 1/60 and 2/60 s

Figure 13

Vertical labyrinth: six channel pathlines depicted in six different colours

Figure 14

Vertical labyrinth: experimental (solid line) and LBM (circles) results for fluid flow propagation along the six pathlines, coloured according to the pathline they refer to

Figure 15

Vertical labyrinth: comparison of randomly selected experimental run (top) and LBM data (bottom) at time 0.75, 2.4, 5.4 and 15 s

Table 1

Poiseuille flow: maximal velocity as computed by LBM vs number of lattice cells across channel breadth

Number of cells	Parameter τ	LBM umaxLB[ms−1]	Rel. error errn [%]	Order of convergence log2(errnerr2n)
5	0.515	1.1992	4.06	–
10	0.53	1.2247	2.03	1.00
20	0.56	1.2436	0.51	1.99
40	0.62	1.2484	0.132	1.95
80	0.74	1.2495	0.034	1.96

Table 2

Advection: convergence of LBM results to theoretical solution L = 0.01 m

Number of cells n	LBM LLB [mm]	Relative error |LLBM−L|L[%]	Order of convergence log2(errnerr2n)
16	10.312	3.12	–
32	10.156	1.56	0.99
64	10.078	0.78	1.00

Table 3

Laplace pressure: comparison of theoretical and LBM solution of Laplace pressure in spherical droplet for several surface tension values

Surface tension σ[N m−1]	Theory Δp[Pa]	LBM ΔpLB[Pa]	Relative error errσ[%]
0.02	0.860	0.877	1.98
0.04	1.720	1.756	2.09
0.08	3.439	3.517	2.27

Table 4

Laplace pressure: convergence of LBM solution of Laplace pressure to theoretical solution for surface tension σ = 0.02 N m⁻¹

Number of cells n	LBM ΔpLB[Pa]	Relative error errσ[%]	Order of convergence log2(errσnerrσ2n)
16	0.9376	9.06	–
32	0.8833	2.75	1.72
64	0.8720	1.43	0.94

Table 5

Capillary rise: comparison between results of LBM simulations and analytic solution for capillary breadth w = 0.0095 m

θ [°]	Capillary w [m]	Theory H [m]	Theory Hc [m]	LBM HLB [m]	Rel. error errH [%]	Rel. error errHc[%]
60	0.0095	0.0107	0.0103	0.0102	5.17	1.22
45	0.0095	0.0152	0.0145	0.0142	6.26	2.13
30	0.0095	0.0186	0.0178	0.0166	10.68	6.53

Table 6

Capillary rise: comparison between results of LBM simulations and analytic solution for contact angle θ = 45°

θ [°]	Capillary w [m]	Theory H [m]	Theory Hc [m]	LBM HLB [m]	Rel. error errH [%]	Rel. error errHc[%]
45	0.0050	0.0288	0.0285	0.0285	1.12	0.08
45	0.0095	0.0152	0.0145	0.0142	6.26	2.13
45	0.0145	0.0099	0.0090	0.0090	9.33	0.30

Table 7

Parameters of test fluids used in experimental campaigns

Experiment test fluid	Horizontal channel sugar solution I	Vertical labyrinth sugar solution II
Ρ [kg m−3]	1,244	1,240
η [mPa s]	41.6	45.4
σ [mN m−1]	60	64

Table 8

Horizontal channel: arrival times at selected locations and statistics of all executed experimental runs

Distance [mm]	Min(t) [s]	Mean(t) [s]	Max(t) [s]
85	0.45	0.48	0.51
215	1.52	1.59	1.67
345	3.33	3.51	3.69
475	6.01	6.34	6.66
585	9.00	9.52	9.97
695	12.78	13.49	14.08
785	16.42	17.34	18.11
875	20.61	21.78	22.74
945	24.31	25.69	26.77
1,015	28.41	30.01	31.26
1,065	31.61	33.38	34.76
1,115	34.94	36.96	38.50
1,145	37.09	39.21	40.81
1,160	38.20	40.38	42.05

Table 9

Vertical labyrinth: propagation of the fluid front along selected pathlines (averaged data from all experimental runs)

Time [s]	Pred [mm]	Pgreen [mm]	Pblue [mm]	Pmagenta [mm]	Pcyan [mm]	Pblack [mm]
0.067	19.12	19.12	19.12	19.12	19.12	19.12
0.133	41.49	46.03	46.03	46.03	46.03	46.03
0.233	49.91	56.79	63.44	63.44	63.44	63.44
0.333	56.16	61.48	68.39	73.14	73.14	73.14
0.4	59.22	63.72	70.90	78.14	78.14	78.14
0.6	66.37	70.50	77.72	93.77	93.77	93.77
0.8	72.20	76.58	84.10	109.17	109.17	109.17
1.0	77.32	82.13	90.02	123.21	123.21	123.21
1.2	81.89	87.02	94.83	133.89	133.89	133.89
1.3	83.97	89.16	96.88	139.48	139.48	139.48
1.4	85.89	91.10	98.73	142.11	142.11	142.11
1.7	90.94	96.22	103.73	146.35	146.35	146.35
1.8	92.41	97.70	105.21	148.07	148.07	148.07
1.9	93.80	99.17	106.93	153.49	153.49	153.49
2.0	95.19	100.77	109.09	160.16	159.45	160.16
2.1	96.63	102.43	111.41	164.27	162.29	164.27
2.2	98.03	104.16	113.57	167.14	163.05	167.14
2.4	100.81	107.52	117.19	173.44	164.11	173.44
2.6	103.51	110.52	120.21	179.72	165.27	179.72
2.7	104.78	111.87	121.54	182.69	165.71	182.69
2.8	105.97	113.14	122.77	186.34	166.25	184.16
2.9	107.13	114.32	123.93	189.60	166.95	185.48
3.2	110.13	117.43	127.04	200.35	169.63	188.75
3.5	112.69	120.07	129.77	212.16	172.47	192.01
3.6	113.49	120.86	130.64	215.13	173.64	193.38
3.7	114.25	121.66	131.48	217.31	174.97	194.89
3.8	114.97	122.42	132.33	218.58	176.29	196.19
4.0	116.32	123.89	133.95	221.56	178.59	198.58
4.2	117.59	125.26	135.52	224.52	180.52	201.07
4.4	118.78	126.54	136.98	227.29	181.89	203.07
4.6	119.89	127.72	138.24	229.66	182.83	204.35
4.7	120.40	128.26	138.79	232.11	182.87	204.85
4.8	120.89	128.77	139.33	233.61	183.10	205.30
4.9	121.38	129.25	139.77	235.93	183.47	205.69
5.0	121.81	129.68	140.24	237.96	183.84	206.08
5.1	122.21	130.11	140.66	240.18	183.73	206.53
5.4	123.32	131.22	141.80	247.64	185.15	207.63
5.7	124.21	132.14	142.75	255.13	186.48	208.80
6.0	124.96	132.91	143.60	262.02	188.08	210.07