Computer Science
Computer Science-2021.gz:143199
Computer Science
cb9d7361e432fc0a45393e13ceecad91841fbf0a
Experimental validation of numerical simulations of free-surface flow within casting mould cavities

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.


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 Experimental validation of free-surface flow architectures while its algorithm is relatively simple and efficient, K€ orner et al. (2005), Sano et al. (2008), Schornbaum and R€ ude (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), , 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 twophase 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€ urey (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 threedimensional 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".

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 EC 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 enez (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 eres 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).

Discretisation
The LBM describes the continuum dynamics by simulating an interaction of virtual particles within an equidistant lattice grid, , 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 ¼ f x; u; t ð Þin time t can be described by the Boltzmann equation T is the position vector and u ¼ ½u 1 ; u 2 ; u 3 T is the velocity vector. Assuming that the distribution function f 5 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 5 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 α.
vf where fα is a function of the spatial variable x and time t only, i.e. fα 5 fα(x, t). In the second step, equation (2) is discretised in space and time as where the collision operator Ω(f α ) on the right hand side is approximated by multiple relaxation time operator ( 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 where w α are the velocity vector weights, which values are given according to the applied microscopic velocity model, Succi (2001), d'Humi eres 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: f c α ðx i ; tÞ ¼ f α ðx i ; tÞ À M −1 SðτÞ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 α f α ðx i þ e α Δt; t þ ΔtÞ ¼ f c α ðx i ; 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: where f g α ðx i ; 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.

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 EC 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 freesurface 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€ urey (2007), while the freesurface 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 distinguishedliquid, 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 χ 5 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 where α denotes the opposite direction to α. The total mass change at the given cell is Δm α ðx; t þ ΔtÞ: where n α 5 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: mðx; t þ ΔtÞ ≥ ϱðx; t þ ΔtÞ þ e → liquid cell; mðx; t þ ΔtÞ ≤ À e → gas cell; where e 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 e 5 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 m err ðx; t þ ΔtÞ ¼ mðx; t þ ΔtÞ À 1 for transformation to liquid; mðx; t þ ΔtÞ for transformation to gas; & (11) must be redistributed among the neighbouring interface cells x b so that where N b is a number of the neighbouring interface cells.

Experimental validation of free-surface flow
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: 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 5 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 w is defined on the grid of cells, then When the value of function w in the direction e α would have to be evaluated outside the computational domain, the following linear interpolation is applied: 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.
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: n correct ¼ n w cosðθÞ þ t w sinðθÞ (17) where n w is the normal vector of the solid wall. The tangent vector t w of the wall is defined as t w ¼ n À n w n$n w ð Þ: 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: 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 b χ which are calculated as 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 b χ, the calculation of (20) can be repeated.

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 f from LBM units to physical units is defined using the converter C f so that 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 Similarly, when u max phys is the expected maximal value of the physical velocity and u max LB is the corresponding maximal velocity in LBM units, the velocity converter is defined as For stability reasons, the value of the maximal velocity in LBM units is set to u max LB ¼ 0:1. In addition, the density converter can be considered in the following form: Since the density in LBM units ϱ LB 5 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 .
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 5 C ν ν LB by the relation The physical value of the time step can be expressed as Δt phys 5 C T Δt LB . The value of Δt LB 5 1 is a common choice.

Experimental validation of free-surface flow
Using the distributions functions f α , the macroscopic flow quantities in LBM units could be expressed as 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).
(7) Rewrite a flag for filled/emptied cells and mark new interface cells (10).

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.

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 where h is the channel breadth. The velocity u is a function of z coordinate only, and the maximal velocity u max ¼ ϱgh 2 8ν is obtained in the centre of the channel, i.e. for z ¼ h 2 . In order to compare the analytical solution with the LBM results, the LBM simulation is performed on a cubic computational domain of h 5 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 5 1 m s À2 , the fluid density ϱ 5 1 kg m À3 and the visocity ν 5 0.001 m 2 s À1 . 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 EC where u max LB 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.

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 3 H 3 2H, where H 5 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 ϱ 5 1000 kg m À3 and kinematic viscosity ν 5 10 -4 m 2 s À1 , while the rest of the domain is void, Figure 4 (left).
The initial velocity u 5 0.1 m s À1 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 5 0.1 s is L 5 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 firstorder 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).

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 where σ is the surface tension at the fluids interface and R 1 , R 2 are the principal radii of the interface curvature (for a spherical droplet R 1 5 R 2 ). 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.    Within the LBM simulation, a cubic computational domain with sidelength of 0.1 m is discretised by 40 3 40 3 40 grid cells. Initially, there is a cube of 30 3 30 3 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 is about 2% and the computed droplet radius is approximately 0.0465 m, which corresponds to theoretically predicted radius 3 4 3 4π À Á 1=3 L.
Convergence of the LBM solution of spherical droplet radius to its theoretically predicted value is studied for the surface tension σ 5 0.02 N m À1 and provided in Table 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 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 θ ¼ w cos θ is expected. In larger capillaries, a correction term must be included, Bullard and Garboczi (2009), where f(θ) 5 9.24 cos θ þ 2.13 cos 3 θ, gðγÞ ¼ 0:834 ffiffi ffi γ p − 0:024γ and γ ¼ ϱgw 2 σ . 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 ϱ 5 1,000 kg m À3 , the surface tension σ 5 1 N m À1 and the dynamic viscosity η 5 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 5 0.0005 m.
In total, two sets of numerical tests are carried out. First, an accuracy of LBM results for the capillary breadth w 5 0.0095 m and three different values of contact angle θ 5 {608, 458, 308} 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 Figure 6. Capillary rise: scheme of liquid level rise in between capillary walls EC Since the capillary breadth w is relatively large, the deviation from the corrected capillary rise err H c should be preferred. Nevertheless, both err H and err H c 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 θ ≥ 458 (when the value of err H c is close to 2% or smaller). However, the accuracy of simulations for low contact angle values, such as θ 5 308, is limited. The second set of LBM tests is performed for the contact angle θ 5 458 and three values of the capillary breadth w 5 {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 err H c almost vanishes for w 5 0.005 m and for w 5 0.0145 m. However, the resulting err Hc does not show any general tendency with regards to value of w.

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 y and Mandys (2019).

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 EC 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 y 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.

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 Distance [mm] Min ( 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 Experimental validation of free-surface flow 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  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 EC 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.

Conclusion
The performed experimental study provides original data on problems of gravity-driven freesurface flow analogous to industrial casting processes. Simplified problems are presented  Table 9. Vertical labyrinth: propagation of the fluid front along selected pathlines (averaged data from all experimental runs) Experimental validation of free-surface flow 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.