1 Introduction
1.1 Jet Stripping of Liquid Coatings
We present here a numerical study of the liquid metal coating process. First, liquid film formation on a vertically climbing wall is simulated. Subsequently – in most cases in the same simulation – we simulate wiping of the created film by a planar air jet. These processes are of major industrial significance e.g. in metallurgy (Takeishi et al., 1995), photography, painting and manufacturing of materials (Bajpai, 2018), where the need arises to control the thickness of the deposit. One of the means to establish this control is the use of a airflow, for example with flat planar jets known as “airknives”. These, employed above the coat reservoir, will act by thinning the film deposed onto the product in a controlled manner. However, their effect is far less predictable once airflow issuing from them becomes turbulent, especially around product edges. Similarly, significant kinetic energy of the incoming turbulent airflow may cause unwanted coat atomization, forcing the operators to lower injected air velocity below certain thresholds which are in practice found empirically. Thus, opportunity arises to optimize the industrial process – at the very least, there is a sustained need for studies of such a configuration.
The process of film formation, which is the basis of the coat formation procedure, has been studied both experimentally and analytically by many authors starting with now classical results of (Landau and Levich, 1942). Analytic solutions were found e.g. by (Groenveld, 1970) who focused on withdrawal with “appreciable” inertial forces (relatively high Re flows) or (Spiers et al., 1973) who has modified the withdrawal theory of Landau and Levich, obtaining improved predictions for film thickness that were also confirmed experimentally. Later, (Snoeier et al., 2008) investigated extensively the film formation regimes in which bulges are formed, focusing on the transition between zeroflux and LLtype films.
In the process of coating, liquid is drawn from a reservoir onto a retracting sheet, forming a coat characterized by phenomena such as longitudinal thickness variation (in 3D) or waves akin to that predicted by Kapitza & Kapitza (Cheng, 1994) (visible in two dimensions as well). While the industry standard configuration for Zinc coating is marked by coexistence of medium Capillary number (Ca=0.03) and film Reynolds number , we present also parametric studies in order to establish if our numerical method influences the film regimes obtained in the target configuration. Note that metallurgical effects (solidification) are neglected, as they don’t play a role in the initial stages of film formation (Hocking et al., 2011).
Once a stable film is formed on the retracted sheet, it can be further thinned by striping/wiping with airflow. The latter, in most cases, will be a turbulent flow, as the high Re in the gas are required to exert sufficient pressure on the liquid coat. Although the airflow effects on the coat can be studied using time averaging (Myrillas et al., 2013), certain instantaneous effects, such as forming of bulges and/or edge effects will not be accounted for. Thus, numerical simulations are a promising tool to supplement experimental studies in this field. One of the first systematic accounts of the jet stripping of liquid coatings comes from (Ellen and Tu, 1984) who have shown analytically that not only pressure gradient acting on the film, but also surface shear stress term plays an important role in the coat thickness modification. (Tuck, 1983) derived analytical expressions for a dependency between jet airflow velocity and resulting film thickness – assuming only pressure gradients plays role in film deformation – and adopting the lubrication approximation for the film flow. The work (Takeishi et al., 1995) provided numerical solutions for velocity and shear profiles at the filmair interface during wiping (using a glycerine solution as the coating liquid).
In 2017 the authors of (Hocking et al., 2011) has analysed the problem numerically using a simplified model (including empirically determined shape functions) and a method of lines to study the modified equations of (Tuck, 1983). They concluded, for example, that disturbances to the coating (including bulges and dimples) above the impact zone will persist more likely for thinner coats, as thick ones ’compensate’ for that with surface tension and solidification intensity.
In this work, we follow the DNS (Tryggvason et al., 2011) approach, i.e. we solve a complete set of Navier Stokes equations describing the flow in both phases (in the onefluid formulation (Delhaye, 1973)) with proper boundary conditions, if permitted by the computational code used. A similar approach has previously been adapted e.g. by (Lacanette et al., 2006), however their 2006 paper was limited to twodimensional Large Eddy Simulation approach. Still, they were able to recover pressure profiles of an impinging jet, as well as give some rudimentary prediction of the splashing which takes place below the impingement area. The authors of (Myrillas et al., 2013) performed a study very similar to Lacanette et al. (2006) – but using parameters of dipropylene glycol as a coating liquid – yielding e.g. profiles of the film in the impingement zone. An even more basic 2D study using the VOF method was published in (Yu et al., 2014), yielding information e.g. about droplet trajectories after impact. In this paper, we continue such a numerical approach, this time applying a threedimensional code with very high spatiotemporal resolutions and grid adaptivity.
1.2 Problem Specification
The investigated configuration is visible in Figure 1. As we can see in the sideview, the nozzleband distance is measured at 10mm in industrial configuration. Nozzle diameter is mm. The proportions in the twodimensional schematic are forgone for presentation purposes, hence the vertical character is slightly more visible in the 3D rendering. Liquid is drawn from the reservoir C at the bottom, which coats the moving band A. Subsequently, air injected from the nozzle(s) B collides with the coated band A and leaves the flow domain below and above the nozzle(s); outlets are drawn in Figure 1 (left) with grayed lines.
Gravity is taken into account, and upward band velocity is in most cases taken at m/s. Regardless of the choice of liquid contained in the reservoir C, band A will be coated, although characteristics of the resulting films will depend strongly on the liquid characteristics. Except where noted, we have decided to choose liquid zinc as the coating liquids. Properties of Zn are assumed, that is surface tension , density and viscosity Properties of the surrounding gas  which in all cases is air  are density and viscosity
As explained below, we introduce multiple sets of boundary conditions in three dimensions. To concisely refer to them, we introduce the following nomenclature to designate the investigated configurations. Three major geometries considered will be termed with
If present, the second lower index may be used to designate the grid resolutions used. This index will equal the power of two corresponding to the maximum refinement used by the Basilisk code described further. And so, for example, stands for the first configuration at equivalent refinement level. Most of the distinguishing features of the three geometries have been delineated in Table 1. In case other quantities (such as injection velocity ) are varied between configurations, it will be designated in parenthesis (example: ) Using above terminology, we can now revisit Figure 1: the configuration presented on the lefthandside is recognized as in 2D, while the rhs of Fig. 1 depicts
Conf.  nozzles  

Our departure point is the full ”industrial” configuration , visible in Fig. 1 on the right. As sketched in Figure 1, we orient the geometry so that is the vertical direction, and air injection takes place along axis with nozzles extended in the directions. As visible in Table 1 this configuration involves both ”airknife” nozzles; additionally there are outlet areas at the and domain walls. Split boundary conditions are used to ensure that fluid outflow takes place e.g. only above liquid bath level. As shown in the table, the thickness of the coated band is kept at m, and the centered wall moves up (m/s). Due to the fact that the extent (depth) of the coated all is smaller than the nozzle depth, the configuration allows the air issues from both nozzles to collide.
Two additional configurations are rendered in Figure 2. As with Figure 1, note that rendering is not fully uptoscale: dimensions used in actual simulations are given in Table 1. The configuration has been created from by including only half of the latter and a symmetry boundary condition at the direction. The width (extent) of the coated wall has also been slightly decreased (from 15 to 5 centimeters) to limit computational cost. Still in the configuration the film is formed gravitationally and the airknifeliquid interaction is maintained. Since the coated wall is placed exactly at only half of its width is included in the configuration, which makes less suited for studies e.g. of the edge effects of the coated band. Instead, more computational resources can be directed at studying the airliquid interactions.
The third introduced configuration, is shown at the bottom of Figure 2. It is a ”synthetic” subproblem, designed for a fully academic investigation of the airmetal impact phenomenon, and the initial stages of the twophase flow postimpact. This is made possible by further reducing the domain size, which now is limited to a cube, encompassing a m deep (extent) fragment of the flat nozzle, a nozzlefilm gap and the coat. The coated wall is not present except as a boundary condition on the direction. In this configuration, film is predefined (at thickness as explained below) and no gravitational coating is present. A combination of outflow and symmetry boundary conditions are used on all domain walls, with the exception of a partial inlet at the . Using the configuration, we further reduce the associated cost of simulating the innozzle flow as well as the gasliquid impact.
2 Description of the Flow
2.1 Governing Equations
In all cases presented here, full NavierStokes equations:
(1) 
are solved, assuming also incompressibility of the flow:
(2) 
In (1), stands for the velocity and signifies pressure. Liquid properties (which vary with the phase) are designated and for viscosity and density, respectively. Symbols and
represent unitary and rate of strain tensors, respectively, with
defined asBody force (gravity) is taken into account and represented by We will occasionally refer to the directions “up” and “down” which in both two and threedimensional simulations are to be associated with axis. Surface tension is taken into account into the presented simulations, and represented in (1) by where is a coefficient, is the curvature of the interface , while is Dirac distribution centered on it. We assume a onefluid approach (Delhaye, 1973), in which density and viscosity can change at and a pressure jump is possible there in case of nonzero surface tension. Below, we will occasionally denote fluid properties with suffixes and (liquid/gas) to denote phases.
2.2 Gravitational Film Formation
Regarding the film formation on a moving wall similar to A in Figure 1, we may assume, after (Groenveld, 1970) that for a film with locally constant thickness equations (1) simplify to
(3) 
Integrating (3) one obtains a parabolic profile inside the vertically moving film
(4) 
with arbitrary Using this profile, one can define a dimensionless flux
(5) 
and dimensionless film thickness
(6) 
subsequently establishing a following dependency between the two:
(7) 
Note that knowing the upwardmoving wall velocity and liquid properties, finding is possible from (6) given that has been precomputed and the flow regime, governed mainly by film Reynolds number
, is applicable. We will employ this to estimate the Groenveld’s thickness (
) of the film when studying its formation in Section 4.2.2.3 Liquidair interaction
Once a stable film is formed on the substrate, it is acted upon by airknives, a process we will briefly discuss below. In general terms (especially for situations when velocity profile inside the film can be assumed known), the approach to modeling airliquid interaction is to write the film equation for . This equation needs to include terms representing gravity, surface tension, as well as pressure gradient and the shear stress Comparing magnitudes of pressure and shear stress with that of surface tension often results in dropping the latter from the model. Subsequently, film equations are solved (steady state solutions are sought) and thickness predictions given, dependant upon and The problem is formulated in such a way that the boundary condition imposed on the moving wall () is . At the interface , author (Hocking et al., 2011) expect
as a velocity condition^{1}^{1}1At this stage, can be seen as a ordercorrect prediction of film thickness, Hocking uses a value of one micrometer., and
for the pressure. This is supplemented by a transport condition for thickness
(8) 
With these assumptions, Hocking et al. Hocking et al. (2011) use linearization and a thin film assumption (treating film thickness as an infinitesimal) to simplify the governing NavierStokes equations to:
(9) 
In the above, is the Stokes number – in the cited work,
and is of minusfourth order – while and are dimensionless shear stress and pressure distribution functions. For this model, interesting empirically established forms of and are presented by (Tuu and Wood, 1996) for pressure:
(10) 
and by (Elsaadawy et al., 2007) for the shear stress:
(11) 
with Naturally, should approximately correspond to the height of air impact zone, while and should be calibrated by supplying correct values associated with gas velocity. (One observes distribution (10) to be similar to simulated pressure bell curves e.g. in Fig. 19a). As far as (9) is considered, (Hocking et al., 2011) apply and as discussed above. Note that the expression in parenthesis in (9) is the flux of coating material; within a thinfilm approximation we might request at certain critical points, as well as thus finding
A slightly different way of obtaining the constraint is to first assume a known, e.g. Poiseuille velocity profile inside the film. This could be written as
(12) 
We integrate (12) to get the flux equation:
(13) 
Note that in (13) the zeroflux thickness can be obtained by zeroing pressure gradient and shear strain terms. Similarly to (Hocking et al., 2011) we can now request the flux (13) to have zero derivative with respect to leading to
yielding a solution for thickness we denote
(14) 
where We rephrase the above result by introducing a pressure to shear strain ratio defined as
with this, (14) becomes
(15) 
Further on, by approximating pressure and strain, for example by
(16) 
we can represent as
with and being pressure and shear strainrelated dimensionless constants. Even if we – somewhat optimistically – calculated using as reference length, in the discussed applications we still could have which would reduce (15) to
This simple result relates thinlayerapproximated thickness to shearstress; the actual thickness approximations for using that formulation will be given below. In the general, nonlaminar case, the relation between and cannot be established beforehand. Still, as a working measure, we can define the average velocity in the reference frame of the moving wall, such that
Using one may examine the balance of air and wall stresses, pressure and shear strain, viscosity and gravity in the form
(17) 
with being a wall friction dimensionless coefficient. By once again zeroing the flux derivative one obtains a result similar to (15) with a related correction
(18) 
where
(19) 
While it is reasonable to expect the actual estimates are nontrivial to obtain; we can however conclude that films thicker than in the laminar case are possible in this regime. Nevertheless, postimpact film thickness values of order of microns should be expected, which pose a significant challenges to computational simulations. Additionally, it must be noted unsteady solutions, as presented in two dimensions by (Hocking et al., 2011) involve wavy structures pushed away from the impact zone; while we don’t present a quantitative description of such dips and depressions, their appearance is expected. This further complicates the task of establishing effective coating thickness above the impact zone. Our estimates of will be given below (see Section 4.3), as well as summarized for the industrial parameters in Table 3.
3 Computational methods
In the research presented here we have applied the “Basilisk” computational code (Popinet, 2015), which is an inhouse, GPLlicensed code whose main developer is one of the present authors (SP). It is a descendant of the “Gerris” code (Popinet, 2009) and as the latter, it enables local adaptive mesh refinement (AMR) (Puckett and Saltzman, 1992) using and quad/octtree type mesh (regular, structured cubic meshes without refinement are also possible). The code is optimised for speed (which differs it from Gerris) and capable of both OpenMP (single node) parallelism and MPItype (multinode) operation. Most recently, Basilisk has been applied e.g. to model compressible flows connected to bubble dynamics (Fuster and Popinet, 2018), propose singlecolumn models in meteorological simulations (van Hooft et al., 2018), or simulate turbidity currents (Yang et al., 2018). We conclude our description of the code by briefly remarking about two features that make it stand out: firstly, it is a multiequation solver, i.e. a broad framework that allows choosing equations to be solved, making it defacto a multiphysics code. Secondly, it contains a builtin parser/lexer providing “targeted”, minimal recompilations for the configuration currently used.
NavierStokes equations are solved using a well known projection method (Tryggvason et al., 2011) with a procedure similar to that applied in Gerris (Popinet, 2003, 2009)
. Centered discretization is used for all scalar and vector fields, with additional facecentered values defined for
which are used e.g. to ensure divergencefree condition during mesh refinement. For consistency reasons, advection term of (1) is defined and calculated on cell faces as is Advective fluxes are obtained using the BellCollelaGlaz advection scheme (Bell et al., 1989). Discretizations are generally finitedifferencing up to second order unless noted otherwise. The RungeKutta scheme is used for time advancement, and a certain optimisation of Poisson equation’s solution is given by implementing the Multigrid (MG) method (Brandt, 1984).The Volume of Fluid (VOF) method (Tryggvason et al., 2011) is used to track the interface using geometric interface reconstruction (Aniszewski et al., 2014). In this method, fraction function (equal to one or zero in either phases ) is passively advected with the flow. Grid cells with fractional values are those in which interface is geometrically reconstructed and represented by a line/plane (in two and three dimensions, respectively). Note that and are usually local functions of Interface curvature is also computed from using the HeightFunctions method (Afkhami and Bussmann, 2008; Popinet, 2009) taking into account proper treatment close to solution boundaries.
Basilisk’s procedure for local mesh adaptation employs a wavelet transform of a given scalar field to assess the latter’s discretization error. If the error is above the userspecified threshold, the grid is locally refined by subdividing it onto four (quadtree) and eight (octree) subcells and performing a prolongation of the coursermesh scalar onto children cells to obtain their initial values (the inverse process is termed restriction). For the simulations presented herein, we use and fields’ error as the refinement criteria with and error thresholds, respectively.
3.1 Ensuring Momentum Conservation in TwoPhase Flow
The momentumconserving methods (Vaudor et al., 2017) derive from a variant of VOF (Hirth and Nichols, 1979) method originally proposed in (Rudman, 1998) to treat twophase flows with considerable density ratios. General idea is that instead of a simple incompressibility assumption
(20) 
we instead write the mass transport equation in full, as is done in compressible formulation (Pilliod, 1996):
(21) 
using also the conservative form of the NavierStokes equations (not shown) (Vaudor et al., 2017), which contain the momentum term
(22) 
Subsequently, in implementation, we calculate density from the fraction function definition:
(23) 
instead of the other way around as it is done traditionally (Hirth and Nichols, 1979). The way in which the momentumconserving methods stand out from traditional twophase NavierStokes equation models using VOF is that subsequently, the products found in both (21) and (22) are calculated consistently in the same control volumes. This can be nontrivial if staggered grid (Tyliszczak, 2014) discretizations are used, and can be solved either by gridcell subdivisions (Rudman, 1998) or using subfluxes of fraction function (Vaudor et al., 2017). Thus consistency between transports of mass and momentum are ensured numerically, resulting in a far more robust computation.
3.2 Implementation of embedded solids
Problem geometry illustrated in Figure 1 is nontrivial, due to the fact that flow is expected to take place around walls of the coated band, as well as the edges and corners defining the flat nozzle, i.e. space containing embedded (or immersed) solids, and the computational code used must allow for this. We use a rudimentary technique of locally modifying the velocity field for this purpose. Local modification of scalar fields is a relatively simple technique used when simulating largescale systems involving solids (LinLin et al., 2016). It is a strongly simplified variant of the Immersed Boundary Method (IBM) of Peskin (Peskin, 2002), which does not modify the grid data structure and is thus compatible with MPI protocol. If we denote the interior of the solid contained by boundary by we can note:
(24) 
that is, all velocity components are set to zero within the solid. As long as no provisions are needed for the crude approximation provided by (24) yields satisfactory results (LinLin et al., 2016). A moving wall can be prescribed by using a nonzero () righthand side in (24). Note however, that pressure is not modified in any way inside the solid which, in principle, may result in its incorrect values especially at boundary This could be addressed for by locally modifying pressure gradients, which in a physical sense is equivalent to defining a certain force which would only be nonzero at the boundary (Gibou and Min, 2012). This however complicates the technique to a degree comparable with implementation of domain reshaping, as optimally, it should be followed by removal of the interior points from the grid.
Instead, we note that for geometries presented – even the most complicated setup – the domain interior is merely a sum of cuboids: it contains no inclined nor curved surfaces. The noslip condition at the surface of the substrate wall moving with velocity can be reasonably approximated using (24). Thus, for the current calculations we adopt this simple technique.
3.3 Spatially Restricted Refinement
To limit the associated CPU cost of the grid refinement, we have employed additional technique of spatially restricted refinement (for short, we will use the abbreviation ’SRR’ below). Using SRR is straightforward. The quad/octtree data structure in Basilisk results in subdivisions of cells into four/eight subcells as the grid is refined in two or three dimensions respectively. The entire domain is a 0level (parent/root) cell with four/eight 1level subcells and so on. The subdivision is performed based on criteria stemming from error estimation on chosen scalar fields performed in wavelet space (usually, components of and/or VOF fraction function). If refinement criterion is locally fulfilled, Basilisk will keep refining the grid up until reaching the maximum allowed level. The SRR technique locally limits this maximum grid level using a spatial criterion. This means larger discretization errors are intentionally allowed far from regions of interest. Thę latter regions have to be predefined before the simulation. Then, dynamic grid refinement will act as usual, the only difference being that refinement to maximum level will take place only in chosen domain subareas while outside of them lower maximum level is forced. This tactic of refinement situates the presented simulation between the blockbased (Lakehal, 2010) and pointbased (Popinet, 2009) mesh refinement.
4 Results
4.1 Planning of the Simulations
The full airknife configuration poses numerous challenges for computational simulation. Firstly, it comprises of subprocesses – such as interaction of turbulent structures with planar interface – which are very demanding on their own. This is either for reasons pertaining code stability (Vaudor et al., 2017), reliability of results (Tyliszczak et al., 2008) or CPUcost (Aniszewski, 2011). Secondly, the geometry of the problem, as presented in Figures 1 and 2 results in a complicated flow. The latter includes a range of scales – from domain size to liquid sheet thickness – that require very fine resolution. However resolution could be limited only in region of interest, which amounts to a relatively small part of the simulation domain. For this reason, we have implemented a broad campaign of simulations focused on individual subproblems. For reasons of brevity we will only present here a subset of the obtained results, namely:

A film formation study: (w/o the air injection nozzles) and in 2D;

A brief, 2D validation on the dynamics of the jet impinging on flat plate ( w/o the liquid phase);

Studies of film formation and airknifeliquid interaction with ”relaxed” and industrial parameters.

Simulations of the full configuration using and geometries with varying spatial resolutions and injection velocities.
In the above the “relaxed” parameters simulations assume a decreased and as a means of preparation, converging to the final solution with increasing dimensionless numbers. For reference, Table 2 contains parameters for both industrial and relaxed parametrisations of the considered problem. Most important differences include an order of magnitude lower liquid density and higher : both of these contribute to sway the balance between gravity and liquid uptake towards the latter. This subsequently leads to a thicker film formed, thus decreasing associated CPU cost needed to perform simulation. (For the same reasons, in gas phase, velocity is decreased twofold in relaxed parametrisation.) This results for example in the zeroflux thickness of the film in relaxed parameters being fourteen times that of its value in industrial parameters.
Case  

Unit  (kg/m)  (kg/m)  (Pas)  (Pas)  (m/s)  (m)  (m)  (m/s) 
Relaxed  
Industrial 
Additional difference between the relaxed and industrial configurations is the coated plate thickness, it is held at mm for the relaxed variant and mm in industrial. Nozzle wall thickness is configured analogically. Both changes facilitate the implementation of simulation geometry in the relaxed case, meaning that coarser grids suffice to implement (24) formulation as more gridpoints end up contained in the region.
4.2 Film formation studies
Film formation studies focus on the steel band emerging from the zinc reservoir. As said above, this is implemented using the domain reshaping technique of Basilisk. Such studies allow for observation of e.g. edge effects at the stage well before the initial coat is modified with airknives (Ellen and Tu, 1984). Even by studying this problem in two dimensions a lot can be learned e.g. about the flow inside the film. We can define the film Reynolds number, describing internal liquid flow as
(25) 
where is liquid density, is the upwardmoving band (wall) velocity, and is a zeroflux film thickness (Groenveld, 1970), i.e. that at which liquid fluxes associated with wall movement (upwards) and gravity (downwards) balance out. Thickness can be found by assuming parabolic velocity profile and comparing dimensionless flux and film thickness, leading to
(26) 
Using (26) and calculating from (25) we arrive at values of and for industrial parameters. Indeed, one could say simulations prove that industrialclass metal coating is a man made system on the edge of criticality, as this is very close to critical values delineating laminar and turbulent film formation regimes. A slightly more delicate interpretation is suggested once we modify our expectations towards film thickness as follows.
Assuming that the withdrawal is dominated by inertial forces, one can employ Groenveld’s analysis mentioned in the context of equation (6). In (Groenveld, 1970), values of and characterizing ”industrial” parameters place our case in the high regime, for which Groenveld proposes values of and at and , respectively. Using these with (5) and (6) one can estimate the associated thickness We will use this value below as a rough estimate of the expected film thickness for gravitational withdrawal simulated in this work. Using thus calculated thickness value we may modify (25) like so:
(27) 
While this value is three times smaller than , one could expect the film to be at least in the intermittent regime.
Figure 3 presents the interface geometry at simulated time values s for the industrial parameters simulation. (This simulation is configured in such way that the upwardmoving band is defined as a boundary condition, so no solid embedding technique is needed.) It is visible that in this case the film head penetrated upwards somewhat faster than would suggest; we could attribute it to the boundary condition for the function (Afkhami et al., 2018). Wavy character of the film is easily observable, especially for the final curve corresponding to s. This is emphasized in Fig. 3b showcasing the very same curves with logarithmic scale used for the axis. Moreover, dashed line in Fig. 3b, representing Groenveld’s prediction using high theory is reasonably approximated by our result, save for the aforementioned wavy film character. More specifically, the recovered Groenveld’s thickness is i.e. less than ten percent of the thickness discussed in (Myrillas et al., 2013) in the context of dipropylene glycol, and requires a substantial grid resolution to resolve. The result visible in Fig. 3 has been obtained with levels of refinement. Domain size has been m (only a part is visible in Figure 3). Thus, an individual grid element has the size of
resulting in approximately four grid elements per film thickness at its thinnest point.
Figure 4
shows the creation of a boundary layer for the various moments in time (in the
s range) of the same flow. The profiles have been sampled at or above the reservoir. The velocity profile remains parabolic, however it clearly becomes steeper for s with an apparent plateau extending for suggesting a detachment of the layer adjacent to the plate (Snoeier et al., 2008). In Fig. 4 we additionally compare the profile for s with analytical expression (4) (dots) using Consistency is visible especially closest to the wall, suggesting that the final profiles lend themselves well to those assumed in (Groenveld, 1970), as hinted previously by Figure 3. This serves as a convincing argument that the film evolution is reasonably well described by the high theory.Moreover, in Fig. 4 profiles are sampled only for (i.e. inside the liquid film). Thus, for each of the lines, the abscissa of its righthand endpoint corresponds to the film thickness at m. As one can observe for we have m whereas for the thickness drops, suggesting a bulge has passed over the point and retracted.
Using the equivalent grid, we have performed a threedimensional simulation to study film formation. Its results are presented in Figure 5, which could be seen as a 3D analog of the interface geometry presented above in Fig. 3a. Similar time instance, s is chosen in Figure 5. A heavily ”rugged” film surface is easily recognizable in Fig. 5a, in which it has been colored by the vertical velocity component As we can see, distinct liquid boundary layer develops directly adjacent to the wall, traveling with velocity . This is fully consistent with liquid velocity profiles presented in Fig. 4 for s. As we get further from the boundary layer, velocity at which the film is climbing drops sharply; Fig. 5a indicates also that surface material crumbles back into the bath (blue areas close to the reservoir height). We have included, as an inset in Fig. 5b, an isosurface for the zero vertical velocity , rendered in turquoise against the gray interfacial surface. (Note that occurs as well in the gas far from the coated wall. For this reason, parts of the isosurface were removed from Figure 5b artificially to not obscure the view of the coated wall area.) In this way, we are able to approximate the stagnation height for s as m e.g. m above the bath level. Above this height, all flow is upwards. The interface formations visible throughout the height of the film surface seem sufficiently resolved and not numerically induced. For example, halfway through the film height in Fig. 5 film thickness is of order m (or eighty times the grid size at levels of grid refinement).
Even using a slightly less refined grid ( levels of refinement, or equivalent), we still observe a wrinkling of the interface as well, mostly in the configuration which involves the coated edge. The evolution during the phase of fully developed film is visible in Figure 6, prepared using a equivalent grid. In this phase the band is fully coated, while some “dimples” appear close to the reservoir surface once zinc is drained. Only a very thin layer of zinc is deposited close to the band edges, as can be seen by the surface color which corresponds to . The surface of the film undergoes progressive distortion starting from the side of coated band. This applies especially to the coated and walls, in which wrinkling appears progressively further from the band edges. The turbulent nature of the film is evident, along with fully threedimensional character of the wrinkles/waves. To our knowledge, this is the first published result of a 3D coating simulation including the edge, and is far higher than the previously published 2D results (Myrillas et al., 2013; Lacanette et al., 2006). We observe the wrinkles similarly in the film thickness profiles presented for all configurations further on (e.g. Figure 15).
We continue our look at the physics of the threedimensional film formation with Figure 7, which contains velocity profiles for the vertical component () and the transverse component () along the wall height – only the range of is included, as all values included are smaller than s. At that time, the liquid reaches roughly to m, consistent with Fig. 3. Note that profiles are averaged, and include data only for m, thus the measurement window supplying data for Figure 7 includes only the closest proximity of the coated wall. In Fig. 7a, we observe a transition from a rather smooth profile at s to a much more varied, at final pictured stages. Notably, we observe a stagnation region forming close to the bath level (itself drawn with a dashed line) which is consistent with interface geometry observed in Figure 6. It is expected that velocities are present in this region further from the wall – this however has not been captured with the profile measurement window. Average values are consistent with Figure 6 as well (note that gas velocity is also taken into account in Fig. 7). We now focus our attention on the curve for at s (red color in Fig. 7a). This curve, although calculated using a threedimensional simulation, is comparable with Fig. 4 (curve for s). If, using the latter of the mentioned curves, one calculates a mean value (for ) of , it is equal to m/s. This value should be at least comparable with Fig. 7a taken for s and m; in fact, we find which is within ten percent of the twodimensional simulation. The slight discrepancy might be attributed to the averaging in three dimensions; e.g. presence of the coated band edge, as well as wrinkles pictured in Fig. 6.
Further evidence of the strictly threedimensional character of the film is found in Figure 7b, showcasing the profiles of the transverse velocity component . While close to the beginning of the flow at s (blue squares) this component is nearly zero, oscillates with increasing amplitude in the entrainment zone as time progresses, and remains negative everywhere below the bath level. That is to say the net flow of the liquid layers contacting the coated wall is from the coated edge towards the symmetry plane (at ). As the film forms and its top edge moves further from the bath, transverse net flow is positive, i.e. towards the coated edge, which is consistent with Fig. 6 and explains the rugged surface of the film in the edge neighborhood. Summarizing, it is obvious that at this values, threedimensional effects are strongly present and decisive in determining the liquid flow character.
Finally, note also that Figure 7 features the most resolved of the 3D simulations presented in the paper, at which is locally equivalent to a grid of points ^{2}^{2}2Due to CPU time and memory restrictions we have not continued this simulation into the injection stages..
4.3 Singlephase Impinging Jet Study
A brief study has been performed on the velocity and pressure profiles obtained in boundary layers of a (twodimensional), singlephase jet impinging on the flat plate. This is motivated by the need to ”calibrate” the analytic predictions for the result of jet interaction with film. More precisely: velocity, pressure and shear stress profiles can be applied to extract coefficient when calculating the film thickness above the impact zone (Tuu and Wood, 1996). Note that in this 2D simulation, nozzle walls have been defined as slightly thicker (mm instead of mm used in the ”industrial” configuration); this should have no effect on the pressure distribution in the airwall impact zone.
To begin with, the velocity profile in the wall jet region (Gauntner et al., 1970) has been extracted from the simulation using ms We have applied a combination of time and ensembleaveraging in order to ensure smoothness of profiles (in most situations 15 simulations were ensembleaveraged). Timespans used for temporal averages were relatively large: of the order needed for largest vortices to leave the domain. In Figure 8a an example of wall velocity profile is presented.
Figure 8a presents the profile of the velocity component parallel to the wall (coated band in full simulation) normalized by the mean value taken at distance from the stagnation point. The curve is accompanied by a fit to the analytic prediction presented by (Ozdemir and Whitelaw, 1992), namely:
(28) 
where we assume m/s as per problem specification, while coefficients and are taken and respectively. Value of needs to be set such that it corresponds to the position in the outer layer where velocity is half of the maximum recorded between the two layers. Clearly, we observe in Figure 8a a boundary layer whose thickness is about mm.
Figure 8b displays the spatial distribution of mean pressure normalized by the dynamic flow pressure, or
for the impingement simulation. One can observe correct symmetry in the distribution, as well as easily recoverable stagnation point directly opposing the nozzle outlet. Zooming into Figure 8b reveals that pressure changes sign very close to the wall, which is consistent with velocity curves predicted by (28) and the existence of boundary layer. According to (Tuu and Wood, 1996), peak pressure evolves with the distance from the nozzle as
(29) 
Values presented in Fig. 8b are about half of predicted by (29), meaning that the potential core is not resolved well enough for This warrants an increase in refinement, and is one of the reasons for using at least levels of refinement in majority of the simulations. Profiles presenting raw pressure values are shown below, e.g. in Figure 19a.
By fitting the simplified and curves, constructed using (16) to the results obtained in this section, we were able to establish the values of the and coefficients mentioned in 2.3. For the industrial parameters, they amount to This results in the correction of (15) equal to . The abovenozzle film thickness is thus between and for this configuration with thin layer assumptions. By using an iterative computation with the ColebrookWhite equation (Menon, 2015) as is traditionally done in Moody diagram applications (Moody, 1944), we have estimated the wall friction coefficient for (19) at resulting in a value for at which is a correction associated with turbulent film regime; this increases the expected thickness (18) above the injector by no more than compared to the thinlayer approximated prediction. This small difference is in accordance with predictions of (Tharmalingnam and Wilkinson, 1978) who have investigated a similar case of a rotating roll coating. They show that including or disregarding inertial effects changes the resulting thickness prediction only very slightly, especially for low capillary numbers () which is the case in our work ().
4.4 Threedimensional Wiping Simulations. Configuration
In this subsection, we present the first results of the threedimensional simulations performed. This includes geometry specification given in Table 1 and Figure 1a. This first presented simulation ran at refinement level; with , translating to gridcell size of
(30) 
or 8 cells in nozzle diameter . This is not sufficient to resolve turbulent flow within the nozzle, nor the thickness, however Groenveld’s thickness of m is in reach. Simulations of this (and higher) resolutions are possible using the SRR technique described above; refinement level is decreased by as much as  levels (or times) far from the region of interest, e.g. close to outer domain borders. At the equivalent resolution, it is reasonable to expect liquid bulges and wrinkles whose thickness surpasses these should be captured in current simulation.
To further examine the influence of the airknives, we will proceed to threedimensional, macroscopic visualization. An example is contained in Figure 9a which displays interface geometry at s. Nozzle locations are drawn using shading and black outlines in this view – this is done in postprocessing and only for orientation purposes. Additionally, a cutplane is positioned in the backdrop (parallel to coordinate) colored by vorticity. At a relatively wide impact zone is already visible, with individual droplets ejected from the film, as well as rich wrinkling. The structure of the air trace is threedimensional: even if its character is homogeneous above the nozzle, below it we see two zones with larger traces. Additionally, edge area is visibly atomized. Figure 9b emphasizes the consequences of certain geometrical differences between the industrial and relaxed parameter sets. Thicker coated plate is visible; the coat on the plate edge is seemingly not disturbed except in the impact area where it interacts directly with the turbulent structures resulting from collision of air that emanates from opposing nozzles. Moreover, some liquid deposits on the nozzle walls (the nozzles are not rendered in Fig. 9b) partly obscuring the view. Large amounts of the coating material crumble down below the impact zone, resembling the ”peeling” effect observed in (Myrillas et al., 2013).
Another comparison of relaxed and industrial parameters (as defined in Table 2) is visible in Figure 10. Due to a difference in timescale, picture visible in Fig. 10a has been chosen based on the width of the impact zone.
The main focus of the comparison is however the degree of film atomization which, in case of relaxed parameters, is visibly higher. This can be attributed to smaller density of the liquid phase, facilitating the momentum transfer from the gas. Also, film deposit is far thicker to the left ( direction), which, at the same resolution, means that zeroflux film is better represented in Fig. 10a than in b which is consistent with the Re being much lower in the relaxed case: (a) versus in (b) (“optimistic” estimation based on nozzle diameter ) or (a) versus (b) for the film (based on and ).
Continuing the flow analysis for the industrial case, we turn our attention to Figure 11 which displays film geometry at approximately s. By this time, the lower bulge (A) starts forming (below the nozzle level) leading to the onset of backflow into the reservoir. It is visible at the side of the film as a distinct line in Fig. 11a. Basing the measurement on the profile of the fraction function , we estimate the width of the impact zone at cm. By ”impact zone” we understand an area (AB in Fig. 11) with minimum not the entire area of gas/liquid interaction which, as visible in Figure 11 is much larger: film shearing takes place in a nearly mwide area above and below the airknives.
No significant edge effects are visible yet (C in Fig. 11, the wavy edge film structure is an early effect of film formation). However, inspecting Figure 11b we may conclude that the edge film is nearly entirely atomized. In Fig. 11b the same film is seen from a different viewpoint: along the axis, centered at the impact zone. A certain perspective shortcut effect takes place in Fig.11, as droplets close to the viewpoint, i.e. on the plate edge, seem bigger than those far from it. Besides, the view contains an apparent accumulation of droplets from all plate depth (i.e. most droplets along depth of m of the plate are visible) which seems to contradict Fig. 11a if the said perspective effect is not taken into account.
The total CPU cost of the 3D, simulation presented here is approximated at CPUh, only twice the amount of 2D simulations presented in previous subsections. However, maximum refinement level is four times lower here, and simulated physical time is only about one tenth that of the 2D simulation.
4.5 The configuration results
In this section, we present results pertaining to the configuration. As mentioned above, this configuration enables us to focus on the airliquid impact study in more detail as, as long as the mean flow is considered, the configuration has an inherent symmetry. Placing the coated wall in the corner of a cubic domain, we use the SRR technique to coarsen the grid proportionally to the distance from the coated walls. This, as the simulations below confirm, has proven sufficient to dampen the turbulent flow far from the zone of interest and prevent backflows.
In Figure 12 we present a visualisation of the macroscopic shape of the interface for the simulation performed using a equivalent grid. This simulation corresponds to industrial parameters and is the analog of results mentioned above e.g. in Fig. 11. In three subfigures, instantaneous shapes are visible for (a) s (b) s and (c) s. Each of the pictures presents two separate view: an isometric one on the lefthandside, and a sideview (looking along axis) on the righthandside. The cutplane positioned at the domain wall is colored with Varying cell size in the vorticity cutplane is, of course, a consequence of employing the aforementioned SRR technique to limit adaptivity in regions above and below the nozzle. Full resolution is maintained inside the planar nozzle and within one mm of the coated wall. As visible in the r.h.s. images of Fig.12b and Fig.12c, the grid coarsening affects interfacial formations as well: the ejected droplets and ligaments are represented with a coarser grid the further from the coated wall.
Directly after the air contacts liquid in Fig.12a, we note a distinct imprint of the nozzle shape on the liquid. Three longitudinal bulges are formed: one below the nozzle, one directly opposing the air outlet and lastly, a small bulge is formed above the nozzle. A mere two milliseconds later, as shown in Fig.12b, the central bulge  whose liquid is ”trapped” by the airflow, has been completely atomized, turning it into a cloud of droplets and ligaments. (This is shown particularly well in the sideview.) This last result is consistent with that of (Yu et al., 2014), who have (in 2D) investigated a flow characterized by a higher of with lower density ratios. Their results show a distribution consistent with a cloud of droplets – with temporal averaging, it is displayed as a bulge.
The atomization process results in most of the liquid droplets being rejected out of the field of view. Some examples of fastmoving “glider” droplets are visible as traces just below the nozzle in Fig. 12b and c. In the meantime, the lower and upper bulges move away from the nozzle. In Fig. 12c, we note that the upper bulge has, by s advanced approx. mm upwards, and has been considerably smoothed. Compared to the lower bulge, there is almost no atomized material near the upper one. Meanwhile, as suggested by the righthandside view in Fig.12c, the material below the nozzle is partly stripped from the wall and immediately atomized. Fig. 12b suggests that most of the droplets in the impact zone originate from the atomized middle bulge material. Subsequently, the number of droplets below the nozzle in Fig. 12c is far smaller than visible in Fig.12b suggesting that atomization visible in Fig. 12b is a transient phenomenon.
Above the level of the nozzle and between the bulges, a thin film is formed, covered by a threedimensional wave structure as visible in Fig. 12c. At this resolution, we have which, compared with the prediction for (using thin film approximation (15)) or up to (using (18) means that thickness is about half of . Fortunately, within the Volume of Fluid framework, maintaining cells that hold a film thinner than is possible (Tryggvason et al., 2011), as well as advecting such structures. Thus, we conclude that in certain regions of the wall the thinner film is correctly represented.
Atomization of the film occurring at the first instance the airliquid contact might be investigated looking at the nondimensional numbers characterizing this interaction. While the film Reynolds number characterizes mostly film formation, we formulate the Weber number involving gas velocity, as follows:
(31) 
with standing for film thickness. Definition (31) is first applied to industrial parameters characterized by
Value  

Feeding the zeroflux thickness (26) into (31) one obtains If, instead, we settle upon using Groenveld’s thickness (31) yields Both these values seem consistent with a regime in which atomization might be expected ^{3}^{3}3For the “relaxed” configuration presented previously, using is more justified as the film Reynolds number is three times lower. Doing so, we obtain . Values of film thickness calculated using various definitions are given in Table 3, which contains also resulting values of the film Reynolds number as well as Weber number. As the atomization effect has not been reported previously (Myrillas et al., 2013; Ellen and Tu, 1984) we have decided to study it further. This is motivated by the fact that similar liquid breakup could be induced numerically, e.g. by inconsistent momentum transfers (Vaudor et al., 2017), curvature calculation errors, or not accounting for interactions between liquidgas interface and vortical structures in the latter phase (Aniszewski, 2016). Thus, we have included simulation configured as by decreasing air injection velocity. This configuration is characterized by and which, again, places the system just below the ”edge of criticality” as in a context of in our film formation study. We follow up with an examination of the flow at a decreased injection velocity.
Figure 13 presents the simulation roughly milliseconds after airliquid impact. In the Figure, white shaded isosurface represents the liquid interface, while several blue areas depict the (uncoated) moving wall^{4}^{4}4In Fig.13a, the bath level is overexposed (i.e. rendered as white) due to specific light positions in the visualization.. The far view presented in Fig. 13a confirms again the turbulent film character below the impact zone (denoted ”1” in the Figure). Interface geometry in the entrainment region 1 is comparable to Figure 6, and should not be associated with the airliquid interaction. The impact zone is visible above as an area with horizontal wrinkles (Fig. 13a:2 and Fig. 13b). Looking closer at the impact zone we note a small number of gaps (denoted 3) in the film, mainly close to the edges of the coated band. Defects may results from the expected film thickness being not fully resolved. There are however visible edge coating defects (denoted 4 in Fig. 13b) not likely associated with airflow. This is consistent with Fig. 6 and seems to suggest that not only increased resolutions are required in the neighbourhood of the coated edge, but possibly specific formulation of boundary conditions at the sharp solid edge (singularity).
Regarding the atomization phenomenon at instance of airliquid contact, comparing Fig. 13 with righthandside images in Fig. 12 we note the nearly complete lack of atomized structures for m/s and
Another look at the flow characterized by lower Weber number is provided in Figure 14. Three subfigures present the same flow as pictured in Figure 13 at instances of time with . Figure 14a presents interface geometry in the impact zone almost directly following the first airliquid contact^{5}^{5}5Note that Figure 13 corresponds to the instance of time situated between Fig.14b and c.. Differences caused by the decrease in the Weber number are instantly recognizable: no distinct horizontal liquid bulges are formed along the direction; instead, smallerwavelength disturbances are showing within a gradually broadening region, as seen in Fig. 14c. The sideviews included in the Figure show a significant decrease in the number of droplets, which we quantify below in Figure 16. Overall image of the flow is different than that for . Juxtaposing Fig. 12b with Fig. 14c we note the complete absence of the droplet cloud below the impact area. We suspect that in the low regime the film is merely disturbed by the airflow, while areas above the nozzle are continuously fed liquid; hence no permanent film thinning should be expected in such flows.
Indeed, looking at the film thickness profiles presented in Figure 15, we note that the averaged film thickness in the vicinity of the impact zone (the nozzle level is marked with an arrow around m) has been altered but not significantly diminished, except the area some mm above the nozzle. In the Figure, profiles are shown for two instances: (continuous line) and (black dots). Groenveld’s thickness is drawn in dashed line for comparison purposes. Profiles visible in Fig. 15 show clearly a bulge for . This formation can not be simply associated with the jet influence, as it is present as well in the profile prior to impact; it is more likely that it results from uneven coating. At this height, the film is characterized by dimples (Snoeier et al., 2008) – visible in Fig. 14a, also visible in 2D in Fig. 3 above m – and the coat is of threedimensional character, being on average thicker closer to the edge. This is thus displayed in the profiles below the impact zone. As for the consequences of the impact itself, oscillates in the vicinity of which is fully consistent with instantaneous images in Fig. 14ac and indicates alternating areas of thinning and thickening of the film. (Note that the zerolevel shift in Figure 15 is a correction for the wall thickness of m.) By representing and using gas velocity as in (16), we note their magnitudes for are one order below that for which in the context of (16) and (15) amounts to higher . While at first sight it would seem consistent with results presented in Fig. 14 and Fig. 15, a far longer simulation would be required to establish the actual postimpact film thickness (by widening the area in which thinner film would be established) which was not the objective of this parameter study. As mentioned, thickness is diminished for , with the thinner area coasted by two slight bulges. These formations are visible in Fig.14 as horizontal sets of dimples above the impact zone. In our opinion, both the film thinning for
as well as the bulges are artifacts of the collision of a large horizontal vortical structure with the film. Similar effect should be observed in a longer timescale. Namely, individual, spatially distinct ”craters” are probably created on the film surface by individual vortical structures  separated by distances resulting from the jet flapping frequency.
Figure 16 presents droplet volume distribution for (pink bars) and (yellow bars). In the Figure, minimum cell volume ( at the finest grid level) is denoted with a black vertical line. Clearly, both simulations involve a significant number of ’subgrid’ VOF ”debris” – gridcells containing nonzero fraction function values that cannot be geometrically reconstructed. This is due to the fact that, firstly, turbulent airflow contributes to droplet breakup which continues until grid resolution becomes insufficient. Secondly, the SRR technique makes this mechanism act much more often which can be seen e.g. in the r.h.s. image of Figure 12c. Focusing our attention on the resolved droplets, we note in Fig. 16 that at low injection velocity there is about resolved droplets in total (yellow bars) which is qualitatively different than at higher air velocity (red bars).
An attempt to characterize the influence of impinging gas flow onto the internal velocities of the liquid film is presented in Figure 17. In the Figure, we are looking at the approximated vertical component of liquid velocity obtained using the VOF fraction function , which is equal to inside the liquid. In other words the product disappears in the gas phase, and Fig. 17 shows its profile in the direct neighbourhood of the impact zone ( i.e. two centimeters below and above the impact zone). Two simulations are included with and m/s. Since two flows have slightly different characteristic time scales due to higher in the later, we have compared instantaneous profiles at the instance corresponding to the impact zone width approximately equal m (as e.g. in Fig. 14c). The curves have been averaged over a sampling window mm thick. For both flows, preimpact velocity distribution in the analysing window is similar with , which is caused by a film thinning in the analysing window slightly below the impact zone (i.e. the film is not perfectly flat even before the impact). After the impact, in case of the high simulation (inverted triangles in Fig. 17 one immediately observes the downward flow caused by the gas in the impact zone. Strong upward movement is visible above it. In the case of lower the average values remain positive, suggesting the air knife wiping is far weaker for chosen injection parameters. Note that the overall character of the profile curves drawn using the averaging chosen in Fig. 17 remains comparable to “raw” velocity profiles as presented in Fig. 18.
This concludes our investigation of the influence of the Weber number on atomization process  we conclude that the atomization effect visible at higher is a correct result. The simulated airliquid system responds as expected to the decrease in dimensional number, while other simulation parameters (e.g. grid resolution) are kept constant.
We finish our analysis of the simulations with a brief remark on the computational efficiency. Thus, the approximated computational cost for the Basilisk simulations presented e.g. in Figure 16 was CPUhours (for simulation with m/s) and CPUhours (for m/s simulation).
4.6 The configuration results
The configuration includes a geometry corresponding to a subset of the configuration (like the latter, which itself is a subset of . This time only the area directly adjacent to the central section of the impact zone is included, as the domain size is ten percent that of In this way we are able to resolve e.g. the flow within the nozzle even at refinement levels lower than that used while solving discussed previously. For example at levels of refinement, we have
which is sufficient even for a rudimentary representation of the boundary layer within the nozzle. The boundary conditions however do not involve the liquid metal bath. Instead, the film is predefined at chosen thickness; for the results presented in this subsection thickness is chosen. Resulting film thickness will thus depend only on air knife pressure and shear (Hocking et al., 2011), as there is no flux through the film. Simulations presented in the subsection have been carried out for . Thus, some film wiping is expected to take place after the airliquid interaction, however the coating will likely not be fully removed. To put the simulations into perspective: a thicker film () is included in the configuration, with that is similar to ones used in Fig. 13 and 14.
Figure 18 presents example profiles for the type simulation with In this simulation, gravitational coating of the wall was replaced with a predefined film of thickness Plots are instantaneous, which explains profile roughness; ensemble averaging, due to high associated cost, is not a viable option for this and other configurations. Symmetry of the vertical velocity component in Fig. 18a indicates that air is evenly distributed upwards and downwards of the injection zone – meaning that no jet flapping has occurred yet at , which is ms. This, as shown below, is consistent with film profiles taken at that same instant in time. The lack of flapping phenomenon seems confirmed by the symmetric distribution visible in Fig. 18b.
Profiles presented in Fig. 18 have however been averaged spatially, namely they are sampled using a window which has a one millimeter span in and spans the whole domain in the direction. The same applies to pressure profiles presented in Figure 19. Logarithmic profiles visible in Fig. 19a exhibit a visible progression from the moment of airliquid contact roughly at Total pressure growth is suppressed after , thus we turn our attention to a normalized quantity in Fig. 19b, similarly to what we have done in the “aironly” study presented in Fig. 8. A decrease is visible converging close to the value of predicted by (Tuu and Wood, 1996). Note that density used to prepare Fig. 19b is fraction functionweighted (effectively densityweighted) , namely
(32) 
which is a necessary step to include the phase mixture presence in the sampling window. Improvement in the profile are likely the effect of a much higher relative resolution within the nozzle. Basing on the turbulence characteristics and visualizations obtained in simulations for configuration (not shown), we are concluding that the results are not yet fully converged and substantially longer simulations are needed in order to make sure turbulent channel flow (Stolz and Adams, 1999) is approximated correctly.
These suspicions are substantiated in Figure 20, showcasing the film thickness profiles for the type simulation. Two groups of profiles are shown: instantaneous profiles sampled along axis using the entire depth (span, or ) in Fig. 20a, and profiles sampled using only a narrow strip in the middle of the film or in Fig. 20b. In both subplots of Fig. 20 the zeroflux thickness has been marked using a dotdashed line.
In is important to note that in the configurations film thickness is predefined, as and the coating process is not included. Thus, since film thickness before impact is artificially kept higher than in and configurations, as the latter use gravitational coating. This decision is motivated by the decreased computational cost associated with representing a film of predefined thickness; however simulations using are planned for future research.
Returning to Figure 20, one observes that narrow profiles in Fig. 20b exhibit thickness close to the initial thickness along entire span except the centre of the impact zone. At the same time, profiles taken along entire span (Fig. 20a) describe a visibly thinner film. This, along with the analysis of macroscopic film geometry, confirms that the film is not thinned evenly along the entire span of the domain. In other words, a stratification appears in the nozzle flow suggesting higher mean velocities close to the domain edges. Similar behavior might be observed in a channel flow which is not fully developed, and small scale turbulent structures are mostly present close to domain edges. However, it seems that in this particular case film thinning occurred mostly close to the edges and did not take place in a manner similar to and configurations. This difference, in the context of this paper, should be explained by much lower spatial resolution of the and simulations compared to ( points in for versus for ). For the continued research, inlet conditions for the air in configurations will be generated in such way that fully turbulent channel flow is obtained throughout the nozzle.
Inspecting Figure 20b we find the film thickness in the impact zone at approximately As mentioned above regarding the configuration, values up to order of magnitude higher than those expected at are expected at this injection velocity. However, a more involved analysis of the configuration, which is ”synthetic” in that it doesn’t involve gravitational film formation (or ), makes it clear than the analysis presented in Sect. 2.3 doesn’t apply to this particular case. For example, as and only could be accepted as a ”steady state” solution. Thus, the configuration is useful mostly to examine the impact phenomenon, appearance of atomization or the width of the impact zone which, as seen from Fig. 20b could be approximated, for at which is consistent with the empirical formulation (11) of (Elsaadawy et al., 2007).
To conclude the analysis of , it is possible a new configuration could be proposed, similar to but involving periodic boundary conditions at the extends of the domain, gravity and a more varied initial thickness distribution possibly oscillating around Groenveld’s thickness Such a configuration would in turn allow one to mimic the physics described by (9) or (13) while keeping the computational cost low. However, obtaining the proper definition of boundary conditions for liquid entering and leaving the domain could be nontrivial; leaving it as an interesting option for a future research on this subject.
5 Conclusions
In this paper, we have presented a novel set of simulations of a very demanding, twophase fluid flow whose characteristics closely correspond to that of air and liquid metal (Zinc). The boundary conditions correspond to the air knife jetwiping process in hotdip coating. In many aspects this is a pioneering work: to our knowledge, the only similar calculations published have described a twodimensional case with RANS/LES performed for the airflow (Myrillas et al., 2013) or investigated the film formation only (Snoeier et al., 2008) or, possibly, included a predefined film (i.e. not formed gravitationally). Similarly, for reasons of numerical stability (Vaudor et al., 2017), virtually all preceding attempts included a much decreased density ratio between the phases. Obviously, multitude of practical applications of similar results exist e.g. in metallurgical and automotive industries, however they are strictly proprietary and can not be consulted by general public.
None of these simplifications apply here: calculation accuracy for the methodology presented here is limited only by available computational resources dictating the grid resolution. Full resolution of turbulent flow (i.e. below the Kolmogorov scale) is still too expensive. However, thanks to grid adaptivity, we are able to achieve DNS in limited areas: this claim can be further justified considering the Hinze scale , defined as the ratio of turbulent kinetic energy and surface tension, and estimated for the industrial parameters at m. At the grid resolution used in most cases presented here ( levels of refinement) we obtain proving that is resolved. We can thus claim energy transfer from gas to liquid is at largely resolved, although of course we do not hold such claims for the turbulent flow in the gas itself.
Our results show that – as expected in metal foundry practice – the airflow inflicts a pressure gradient at the liquid layer, and “punctuates” it to a degree controlled by the in the air, and a properly defined Weber number. This gradient restricts the liquid feeding from reservoir, thinning the deposit. Our calculations fall short of resolving the upper film thickness in the full and geometries, case being resolved the best, unfortunately the latter includes no gravitational flow. However, the case clearly displays the thinning effect, it is also observable in the case performed with decreased
We have observed levels of atomization of the liquid metal that were not previously reported in the literature. This phenomenon, to our knowledge, has also not been observed experimentally, which leads us to believe it is a purely transient effect, taking place only as a consequence of the initial gasliquid impact event. It is predicted for , while for values closer to one, liquid wrinkling is observed, visible mostly in the configuration. The appearance of atomization seems thus predicted correctly, instead of being induced numerically. An additional observation is that the liquid material is “milled” (atomized) by the airknife before falling back to the reservoir, and that liquidliquid collisions are aplenty. This is already visible at the level (Figures 9b and 11), and suggest that the coatthinning mechanism is far more turbulent in its nature than known previously.
The three geometries introduced in the paper focus on the twoand threedimensional coat formation and nozzle interactions (), coat thinning and edge effects () and the character of gasliquid impact (). For future research – apart from the aforementioned improvements to the configuration – we would envision working preferable with the configurations in two and three dimensions, with increasing resolutions (and simulated timespans), preferably until full resolution of the thickness is feasible.
Acknowledgements
All the computations for this work have been performed using the French TGCC “Irene JoliotCurie” and CINES ”Occigen” supercomputers. Graphs and visualizations have been performed using the Basilisk View^{6}^{6}6See www.basilisk.fr/src/view.h Blender (Blender Online Community, 2019) and Gnuplot (Williams et al., 2010).
The authors would like to thank Gretar Tryggvason for helpful discussions concerning the final state of this paper. We also thank Stefan Kramel, Stanley Ling, Sagar Pal, Maurice Rossi and Daniel Fuster for interesting discussions on the subject. We thank the participants of the “DIPSI 2018” workshop – held in Bergamo and organized by G.E. Cossali and S. Tonini – for their comments pertaining the early stages of this work. WA wishes to thank A.A.S. for proofreading the manuscript.
Wojciech Aniszewski would like to dedicate all his work on this paper to the dear memory of Antonina Gilska (nee Hessler, 19362019).
References

Afkhami et al. (2018)
Afkhami, S., Buongiorno, J., Guion, A., Popinet, S., Saade, Y., Scardovelli,
R., Zaleski, S., 2018. Transition in a numerical model of contact line
dynamics and forced dewetting. Journal of Computational Physics 374, 1061 –
1093.
URL http://www.sciencedirect.com/science/article/pii/S0021999118305357  Afkhami and Bussmann (2008) Afkhami, S., Bussmann, M., 2008. Height functions for applying contact angles to 3D VOF simulations. International Journal for Numerical Methods in Fluids.
 Aniszewski (2011) Aniszewski, W., 2011. Large eddy simulation of turbulent twophase flow. Ph.D. thesis, Czestochowa University of Technology, Czestochowa, Poland.
 Aniszewski (2016) Aniszewski, W., 2016. Improvements, testing and development of the ADM subgrid surface tension model for twophase LES. Journal of Computational Physics 327, 389–415.
 Aniszewski et al. (2014) Aniszewski, W., Ménard, T., Marek, M., 2014. Volume of fluid (VOF) type advection methods in twophase flow: A comparative study. Computers & Fluids 97 (0), 52 – 73.

Bajpai (2018)
Bajpai, P., 2018. Chapter 7  coating. In: Bajpai, P. (Ed.), Biermann’s
Handbook of Pulp and Paper (Third Edition), third edition Edition. Elsevier,
pp. 159 – 176.
URL http://www.sciencedirect.com/science/article/pii/B9780128142387000076  Bell et al. (1989) Bell, J. B., Colella, P., Glaz, H. M., 1989. A secondorder projection method for the incompressible navierstokes equations. J. Comput. Phys. 85, 257–283.

Blender Online Community (2019)
Blender Online Community, 2019. Blender  a 3D modelling and rendering
package. Blender Foundation, Blender Institute, Amsterdam.
URL http://www.blender.org  Brandt (1984) Brandt, A., 1984. Multigrid Techniques: 1984 guide with applications to fluid dynamics. GMDStudien 85, Bonn, Germany.
 Cheng (1994) Cheng, H.C., 1994. Wave evolution on a falling film. Annu. Rev. Fluid Mech. 26, 103–136.
 Delhaye (1973) Delhaye, J., 1973. Jump conditions and etropy sources in twophase systems. local instant formulation. International Journal of Multiphase Flow 1, 395–409.
 Ellen and Tu (1984) Ellen, C., Tu, C., 1984. An analysis of jet stripping of liquid coatings. Journal of Fluids Engineering 106, 399–404.

Elsaadawy et al. (2007)
Elsaadawy, E., Hanumanth, G., Balthazaar, A., McDermid, J., Hrymak, A., Forbes,
J., Jun 2007. Coating weight model for the continuous hotdip galvanizing
process. Metallurgical and Materials Transactions B 38 (3), 413–424.
URL https://doi.org/10.1007/s1166300790372 
Fuster and Popinet (2018)
Fuster, D., Popinet, S., 2018. An allmach method for the simulation of bubble
dynamics problems in the presence of surface tension. Journal of
Computational Physics 374, 752 – 768.
URL http://www.sciencedirect.com/science/article/pii/S0021999118305187  Gauntner et al. (1970) Gauntner, J., Livingood, J., Hrycak, P., 1970. Survey of literature on flow characteristics of a single turbulent jet impinging on a flat plate. Tech. rep., NASA.
 Gibou and Min (2012) Gibou, F., Min, C., 2012. Efficient symmetric positive definite secondorder accurate monolithic solver for fluid/solid interactions. Journal of Computational Physics 231 (231), 3246–3263.
 Groenveld (1970) Groenveld, P., 1970. Laminar withdrawal with appreciable inertial forces. Chemical Engineering Science 25, 1267–1273.
 Hirth and Nichols (1979) Hirth, C., Nichols, B., 1979. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics 39, 201–225.
 Hocking et al. (2011) Hocking, G., Sweatman, W., Fitt, A., Breward, C., 2011. Deformations during jetstripping in the galvanizing process. Journal of Engineering Mathematics 70, 297–306.
 Lacanette et al. (2006) Lacanette, D., Gosset, A., Vincent, S., Buchlin, J.M., Arquis, E., 2006. Macroscopic analysis of gasjet wiping: Numerical simulation and experimental approach. Physics of Fluids 18.

Lakehal (2010)
Lakehal, D., 2010. LEIS for the prediction of turbulent multifluid flows
applied to thermalhydraulics applications. Nuclear Engineering and Design
240 (9), 2096 – 2106, experiments and CFD Code Applications to Nuclear
Reactor Safety (XCFD4NRS).
URL http://www.sciencedirect.com/science/article/pii/S0029549309005895  Landau and Levich (1942) Landau, L., Levich, B., 1942. LandauLevich film. Acta Physicochim. URSS 17.
 LinLin et al. (2016) LinLin, Z., Hui, G., ChuiJie, W., 2016. Threedimensional numerical simulation of a bird model in unsteady flight. Comput Mech 58, 1–11.

Menon (2015)
Menon, E. S., 2015. Chapter five  fluid flow in pipes. In: Menon, E. S. (Ed.),
Transmission Pipeline Calculations and Simulations Manual. Gulf Professional
Publishing, Boston, pp. 149 – 234.
URL http://www.sciencedirect.com/science/article/pii/B9781856178303000055  Moody (1944) Moody, L. F., 1944. Friction factors for pipe flow. Transactions of the ASME 66 (8), 671–684.
 Myrillas et al. (2013) Myrillas, K., Rambaud, P., Mataigne, J.M., Gardin, P., Vincent, S., Buchlin, J.M., 2013. Numerical modelling of gasjet wiping process. Chemical Engineering and Processing: Process Intensification 68, 26 – 31.
 Ozdemir and Whitelaw (1992) Ozdemir, B., Whitelaw, J., 1992. Impingement of an axisymmetric jet on unheated and heated flat plates. Journal of Fluid Mechanics 240, 503–532.
 Peskin (2002) Peskin, C. S., 2002. The immersed boundary method. Acta Numerica 11, 479–517.
 Pilliod (1996) Pilliod, J., 1996. A secondorder unsplit method for modeling flames in twodimensional compressible flow. Ph.D. thesis, University of California, Davis.
 Popinet (2003) Popinet, S., 2003. Gerris: a treebased adaptive solver for the incompressible euler equations in complex geometries. Journal of Computational Physics 190, 572–600.
 Popinet (2009) Popinet, S., 2009. An accurate adaptive solver for surfacetension driven interfacial flows. Journal of Computational Physics 228, 5838–5866.

Popinet (2015)
Popinet, S., Dec. 2015. A quadtreeadaptive multigrid solver for the
Serre–Green–Naghdi equations. Journal of Computational Physics 302,
336–358.
URL https://hal.archivesouvertes.fr/hal01163101  Puckett and Saltzman (1992) Puckett, E., Saltzman, J., 1992. A 3D adaptive mesh refinement algorithm for multimaterial gas dynamics. Physica D 60, 84–104.
 Rudman (1998) Rudman, M., 1998. A volumetracking method for incompressible multifluid flows with large density variations. International Journal for Numerical Methods in Fluids 28, 357–378.
 Snoeier et al. (2008) Snoeier, J., Ziegler, J., Andreotti, B., Fermiger, M., Eggers, J., 2008. Thick films of viscous fluid coating a plate withdrawn from liquid reservoir. Physical Review Letters 100, 24502–24504.
 Spiers et al. (1973) Spiers, R., Subbraman, C., Wilkinson, W., 1973. Free coating of a newtonian liquid onto a vertical surface. Chemical Engineering Science 29, 389–396.

Stolz and Adams (1999)
Stolz, S., Adams, N. A., 1999. An approximate deconvolution procedure for
largeeddy simulation. Physics of Fluids 11 (7), 1699–1701.
URL https://doi.org/10.1063/1.869867  Takeishi et al. (1995) Takeishi, Y., Yamauchi, A., Miyauchi, S., 1995. Gas wiping in hotdip coating process. Testu to Hagane (in Japanese) 81, 643–648.
 Tharmalingnam and Wilkinson (1978) Tharmalingnam, S., Wilkinson, W., 1978. The coating of newtonian liquids onto a rotating roll. Chemical Engineering Science 33, 1481–1487.
 Tryggvason et al. (2011) Tryggvason, G., Scardovelli, R., Zaleski, S., 2011. Direct Numerical Simulations of GasLiquid Multiphase Flows. Cambridge University Press.
 Tuck (1983) Tuck, E. O., 1983. Continuous coating with gravity and jet stripping. The Physics of Fluids 26 (9), 2352–2358.
 Tuu and Wood (1996) Tuu, C., Wood, D., 1996. Wall pressure and shear stress measurements beneath an impinging jet. Experimental Thermal and Fluid Science 13, 354–373.
 Tyliszczak (2014) Tyliszczak, A., 2014. A highorder compact difference algorithm for halfstaggered grids for laminar and turbulent incompressible flows. Journal of Computational Physics 66, 438–467.
 Tyliszczak et al. (2008) Tyliszczak, A., Boguslawski, A., Drobniak, S., 2008. Quality of LES predictions of isothermal and hot round jet. Quality and Reliability of Large Eddy Simulation, ERCOFTAC Series 12, 259–270.

van Hooft et al. (2018)
van Hooft, J. A., Popinet, S., van de Wiel, B. J. H., 2018. Adaptive cartesian
meshes for atmospheric singlecolumn models: a study using basilisk 180216.
Geoscientific Model Development 11 (12), 4727–4738.
URL https://www.geoscimodeldev.net/11/4727/2018/  Vaudor et al. (2017) Vaudor, G., Ménard, T., Aniszewski, W., Berlemont, A., 2017. A consistent mass and momentum flux computation method for two phase flows. application to atomization process. Computers and Fluids 152, 204–216.
 Williams et al. (2010) Williams, T., Kelley, C., many others, March 2010. Gnuplot 4.4: an interactive plotting program. http://gnuplot.sourceforge.net/.
 Yang et al. (2018) Yang, S., An, Y., Liu, Q., 2018. A twodimensional layeraveraged numerical model for turbidity currents. Geological Society, London, Special Publications 477.
 Yu et al. (2014) Yu, T., Park, J., Moon, H., Shim, J., You, D., Kim, D., Ovsyannikov, A., August 2014. Numerical simulation of liquidlayer breakup on a moving wall due to an impinging jet. In: Proceedings of the Summer Program 2014. Center of Turbulence Research.
Comments
There are no comments yet.