text:
Wojciech Aniszewski
email: aniszewski@dalembert.upmc.fr

Youssef Saade

Stéphane Zaleski

Stéphane Popinet

Planar Jet Stripping of Liquid Coatings: Numerical Studies

Keywords: Coating, Film formation, turbulence-interface interaction, simulation, Volume-of-Fluid

In this paper, we present a detailed example of numerical study of film formation in the context of metal coating. Subsequently we simulate wiping of the film by a planar jet. The simulations have been performed using Basilisk, a grid-adapting, strongly optimized code. Mesh adaptation allows for arbitrary precision in relevant regions such as the contact line or the liquid-air impact zone, while coarse grid is applied elsewhere. This, as the results indicate, is the only realistic approach for a numerical method to cover the wide range of necessary scales from the predicted film thickness (hundreds of microns) to the domain size (meters). The results suggest assumptions of laminar flow inside the film are not justified for heavy coats (liquid zinc). As for the wiping, our simulations supply a great amount of instantaneous results concerning initial film atomization as well as film thickness.

Introduction

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 [START_REF] Takeishi | Gas wiping in hot-dip coating process[END_REF], photography, painting and manufacturing of materials [START_REF] Bajpai | Chapter 7 -coating[END_REF], where the need arises to control the thickness of the deposit. One of the means to establish this control is by the use of an airflow, for example with flat planar jets known as "air-knives". These, located horizontally above the coat reservoir, will act by wiping the film in a controlled manner. However, the effect of the jets is not fully predictable when the airflow issuing from them becomes turbulent, especially around the product edges. The significant kinetic energy of the incoming turbulent airflow may cause unwanted coat atomization or defective coating around the product edges, forcing the operators to lower injected air velocity below certain thresholds -these are in practice found empirically. There is a sustained need for studies of such a configuration for the purposes such as process optimization.

Forming of the liquid film -the basis of the coat formation procedure -has been studied both experimentally and analytically by many authors, starting with the -now classical -results of [START_REF] Landau | Landau-Levich film[END_REF]. Analytic solutions were found e.g. by [START_REF] Groenveld | Laminar withdrawal with appreciable inertial forces[END_REF] who focused on withdrawal with "appreciable" inertial forces (relatively high Reynolds number (Re ) flows) or [START_REF] Spiers | Free coating of a newtonian liquid onto a vertical surface[END_REF] who have modified the withdrawal theory of Landau and Levich, obtaining improved predictions for film thickness that were also confirmed experimentally. Later, [START_REF] Snoeijer | Thick films of viscous fluid coating a plate withdrawn from liquid reservoir[END_REF] investigated extensively the film formation regimes in which bulges are formed, focusing on the transition between zero-flux and LL-type films.

As mentioned, in the process of coating, liquid is drawn from a reservoir onto a retracting sheet, forming a coat. The latter is characterized by phenomena such as longitudinal thickness variation (in 3D) or waves akin to that predicted by Kapitza & Kapitza [START_REF] Cheng | Wave evolution on a falling film[END_REF] (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 Re f > 2000, 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, rapid stages of film formation [START_REF] Hocking | Deformations during jet-stripping in the galvanizing process[END_REF].

As mentioned, significant Reynolds numbers in the air are expected in the wiping stage. Although the airflow effects on the coat can be studied using the time averaging [START_REF] Myrillas | Numerical modelling of gas-jet wiping process[END_REF], certain instantaneous effects, such as forming of bulges, edge effects or film defects 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 [START_REF] Ellen | An analysis of jet stripping of liquid coatings[END_REF] 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. [START_REF] Tuck | Continuous coating with gravity and jet stripping[END_REF] derived analytical expressions for a dependency between jet airflow velocity and resulting film thickness -assuming only the pressure gradients plays role in film deformation -and adopting the lubrication approximation for the film flow. The work [START_REF] Takeishi | Gas wiping in hot-dip coating process[END_REF] provided certain numerical solutions for velocity and shear profiles at the film-air interface during wiping (using a glycerine solution as the coating liquid).

In 2017 the authors of [START_REF] Hocking | Deformations during jet-stripping in the galvanizing process[END_REF] have analysed the problem numerically using a simplified model -including empirically determined shape functions -and a method of lines to study the modified equations of [START_REF] Tuck | Continuous coating with gravity and jet stripping[END_REF]. They concluded e.g. that disturbances of the coating (as bulges/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 [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF]) approach, i.e. we solve a complete set of Navier Stokes equations describing the flow in both phases (in the one-fluid formulation [START_REF] Delhaye | Jump conditions and etropy sources in two-phase systems. local instant formulation[END_REF]) with proper boundary conditions, if permitted by the computational code used. A similar approach has previously been adapted e.g. by [START_REF] Lacanette | Macroscopic analysis of gas-jet wiping: Numerical simulation and experimental approach[END_REF], however their 2006 paper was limited to the two-dimensional Large Eddy Simulation (LES) approach. Still, they were able to recover the pressure profiles of an impinging jet, or predict splashing will take place below the impingement area. The authors of [START_REF] Myrillas | Numerical modelling of gas-jet wiping process[END_REF] performed a study very similar to [START_REF] Lacanette | Macroscopic analysis of gas-jet wiping: Numerical simulation and experimental approach[END_REF] -but substituting dipropylene glycol for the 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 [START_REF] Yu | Numerical simulation of liquid-layer breakup on a moving wall due to an impinging jet[END_REF], yielding information e.g. about certain droplet trajectories after impact. In this paper, we continue such a numerical approach, this time applying a three-dimensional code with very high spatio-temporal resolutions and adaptive mesh refinement.

This paper is structured as follows. In the further parts of the Introduction, we outline the geometrical specification of the setup as well as its physical parameters. Section 2 deals with the mathematical description of the flow at hand. In Section 3, we briefly describe the computational methods chosen for the study. Subsequently, Section 4 presents all the results obtained from simulations, and the conclusions are presented in Section 5.

Problem Specification

The investigated configuration is visible in Figure 1. Dimensions visible in the leftmost illustration pertain our target (or "industrial") configuration. The coating liquid is drawn from the reservoir C at the bottom, and deposits on the vertical band A as the latter moves upwards. Subsequently, air injected from the nozzles B collides with the coated band A, interacts with the film deposit, and leaves the flow domain Ω below and above the nozzle(s); outlets are drawn in Figure 1 (left) with grayed lines.

As we can see in the side-view (Fig. 1 left), the nozzleband distance d n f is measured at d n f = 10mm in the industrial configuration. Nozzle diameter d is 1mm. The proportions in the two-dimensional schematic are forgone for presentation purposes, hence the vertically elongated domain shape is slightly more visible in the 3D rendering (Fig. 1 right). Gravity is taken into account, and upward band velocity is in most cases taken at u wall = 2m/s. Except where noted, we have decided to choose liquid zinc as the coating liquid. Properties of 30 Zn are assumed, that is surface tension σ = 0.7[N/m], density ρ l = 6500[kg/m 3 ] and viscosity µ l = 3.17

• 10 -3 [Pa • s].
Properties of the surrounding gas -which in all cases is air -are density ρ a ≈ 1.22[kg/m 3 ] and viscosity µ a = 2.1

• 10 -5 [Pa • s].
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. Two geometries considered will be termed G i with i = 1, 2. 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, G 1,14 stands for the first configuration at 2 14 -equivalent refinement level. Most of the distinguishing features of the two geometries have been delineated in Table 1. In case other quantities (such as injection velocity u in j ) are varied between configurations, it will be designated in parenthesis (e.g. G 2,11 (u in j = 42) stands for the G 2 configuration on a 2 1 1equivalent grid with the air injection velocity equal to 42 m/s). Using above terminology, we can now revisit Figure 1: the configuration presented on the left-hand-side is recognized as G 2 in 2D, while the r-h-s of Fig. 1 depicts the three-dimensional G 1 .

Our departure point is the full "industrial" configuration G 1 , visible in Fig. 1 on the right. As sketched in Figure 1, we orient the geometry so that y is the vertical direction, and air injection takes place along x axis with nozzles extended in the z directions. As visible in Table 1 this configuration involves both "air-knife" nozzles; additionally there are outlet areas at the z+, z-and y+ 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 1, the thickness h w of the coated band A is kept at 1mm. The position of the coated wall along the x axis is given by x wall in the Table ; it is x-centered in the G 1 configuration, moved leftmost in the G 2 . In all cases, we impose the upward wall velocity u wall = 2 (m/s). Due to the fact that the z-extent (depth) of the coated wall is smaller than the nozzle depth, the G 1 configuration allows the air issuing from both nozzles to collide. This ends the description of the G 1 configuration.

Two additional configurations are rendered in Figure 2. As with Figure 1, note that rendering is not fully up-to-scale: dimensions used in actual simulations are given in Table 1. The G 2 configuration has been created from G 1 by including only half of it and a symmetry boundary condition at the xdirection. In other words, the G 2 configuration is a threedimensional realization of the sketch presented on the lefthand-side of Fig. 1. Going into G 2 , the depth (z-extent) of the coated wall has also been slightly decreased (from 15 to 5 centimeters) to limit computational cost of the simulation. Still, in the G 2 configuration the film is formed gravitationally and the airknife-liquid interaction is maintained. Since the coated wall is now centered at x = 0, only half of its thickness (x-span) is included in the G 2 configuration, which makes G 2 less suited for studies e.g. of the edge effects of the coated band. Instead, more computational resources can be directed at studying the Conf.

L x × L y × L z h w x wall u wall f g # nozzles d G 1 0.25 × 0.65 × 0.25 1 • 10 -3 0.125 2 9.81 2 1 • 10 -3 G 2 0.512 3 0.5 • 10 -3 0 2 9.81 1 1 • 10 -3
Table 1: Distinguishing features of the G 1 and G 2 initial conditions (all units are meters)). air-liquid interactions. Of course, the G 2 includes only a single nozzle.

Description of the Flow

Governing Equations

In all the cases presented henceforth, the full Navier-Stokes equations:

∂u ∂t + ∇ • (u ⊗ u) = 1 ρ (∇ • (µD -pI) + σ nκδ S ) + f g , (1)
are solved, assuming the flow to be incompressible:

∇ • u = 0. (2)
In (1), u stands for the velocity vector and p signifies pressure. The liquid properties are designated by µ and ρ for viscosity and density, respectively. Symbols I and D represent unitary and rate of strain tensors, respectively, with D defined as

D = ∇u + ∇ T u.
Gravity is taken into account and represented by the body force f g . Capillary forces are represented in (1) by σnκδ s where σ is the surface tension coefficient, κ is the curvature of the interface S , and δ S is a Dirac function defined only in the cells containing an interface. We adopt the one-fluid approach [START_REF] Delhaye | Jump conditions and etropy sources in two-phase systems. local instant formulation[END_REF], in which density and viscosity can change at S , and a pressure jump is possible there in case of non-zero surface tension. We will occasionally refer to the directions "up" and "down" which in both 2D and 3D simulations are to be associated with the yaxis. Moreover, we will occasionally denote the fluid properties with suffixes l and g (liquid/gas).

Film Formation and the Air-Knife Theory

In this subsection, we shall present a simplified theoretical description of the air-knife process, that allows us to obtain an estimation of the final coating thickness. We employ the thin film approximation, which, followed by an order of magnitudes analysis and the associated simplifications, brings us to the following form of (1):

µ l ∂ 2 u y ∂x 2 = ρ l g + ∂p ∂ y . ( 3
)
Comparing magnitudes of pressure and shear stress imposed by the air flow with that of surface tension often results in dropping the latter from the model. The boundary conditions associated with (3) are the "no-slip" restriction at x = 0, where u y = u wall , and an imposed shear stress on the surface (x = h) by the air flow:

τ xy = µ l ∂u y ∂x x=h . (4)
Upon integrating (3) and applying the boundary conditions, one obtains a parabolic profile inside the vertically moving film:

u y (x) = - ρ l g + ∂ y p 2µ l x(2h -x) + τ xy µ l x + u wall . (5)
The volume flux through the film, by unit length is therefore:

q = h 0 u y (x)dx = - ρ l g + ∂ y p 3µ l h 3 + τ xy 2µ l h 2 + u wall h. (6)
There is a single length scale in the problem h 0 = [µ l u wall /(ρ l g)] 1/2 . Using this length, one can define a dimensionless flux,

Q : = q u wall h 0 , (7)
a dimensionless film thickness,

T : = h h 0 , (8)
and a dimensionless shear stress,

S : = τ xy h 0 µ l u wall . (9)
The dimensionless effective gravitational acceleration writes as:

G : = 1 + ∂ y p ρ l g . ( 10
)
Subsequently, the dependency amongst these dimensionless groups arises from (6) and writes as:

Q = T 1 - GT 2 3 + S T 2 2 . ( 11
)
Note that (11) can also be used to study the gravitational film formation, that is the formation of a film by withdrawal, before the action of the air-knives (S = 0 and G = 1). We will employ this to estimate the Groenveld's thickness (denoted T G here) of the film when studying its formation in Section 4.2 (see also [START_REF] Groenveld | Laminar withdrawal with appreciable inertial forces[END_REF]). In addition, one can readily derive the zeroflux thickness by additionally imposing Q = 0, thus obtaining T 00 = √ 3. Our main interest however, is to find an estimate of the final coating thickness T + due to the thinning effect of the air-knives. Examining the right-hand side of (11), we see that for a fixed value of the flux Q, the cubic function for the thickness T admits two positive solutions T -and T + . This is valid for every value of y since G and S are functions of that coordinate. There exists a transitional critical point at a certain y = y c , where the cubic function admits a double positive root T c . Thus, for y < y c , the physical thickness would be T -whereas for y > y c , the selected thickness should be T + . In the spirit of what has been done in [START_REF] Hocking | Deformations during jet-stripping in the galvanizing process[END_REF], the condition for this critical transition is:

∂Q ∂T = 1 -G c T 2 c + S c T c = 0. ( 12
)
Knowing the pressure and the shear stress profiles, we expect this critical thickness to occur at y < y c , where S c < 0 and G c > 0 . Therefore, a positive solution to ( 12) is:

T c = S c 2G c         1 -1 + 4G c S 2 c         . (13)
We then approximate the pressure imposed by the air flow and the shear stress as

p ∼ c p ρ g u 2 in j , τ xy ∼ c s ρ g u 2 in j , (14)
where c p and c s are coefficients to be determined from numerical simulations. In addition, we define a length scale l 0 designating the area of the film over which the jet has the most effect. It is typically larger than the nozzle diameter due to the downstream spread of the jet, and can be readily determined from the numerical simulations. Hence, we approximate the gradient of pressure as ∂ y p ∼ c p ρ g u 2 in j /l 0 . In an attempt to find the dominant term in the effective gravitational acceleration G, we find that ∂ y p/(ρ l g) ∼ c p (ρ g /ρ l )Fr where the Fr = u 2 in j /(gl 0 ) is the Froude number, typically very large for the standard operational parameters. Therefore, ∂ y p/(ρ l g) 1, and the effective gravitational acceleration can be written as G c ∂ y p/(ρ l g). We now seek the simplification of solution (13). For that purpose, we approximate the term

G c /S 2 c ∼ (c p /c 2 s )(u wall /u in j )/(mRe l 0 ) (15)
where, m = µ g /µ l is the gas to liquid viscosity ratio, and Re l 0 = ρ g u in j l 0 /µ g is the Reynolds number based on the previously defined length scale l 0 . At the standard operating conditions, Re l 0 is large and u in j u wall . Therefore, typically, G c /S 2 c < 1, leading to a rough estimation of the critical thickness:

T c - 1 S c = 1 | S c | . ( 16
)
Mass is conserved; hence, for a steady solution, the flux is constant throughout the film, and evaluated from (11) to be:

Q c = G c 3S 3 c - 1 2S c - 1 2S c , (17)
employing the previous order of magnitudes analysis. Downstream in the film, and far enough from the jet impingement area, the thickness, denoted as T + has reduced due to the air-knife effect. The film is usually thin enough so that a constant velocity, that of the rising wall, may be assumed throughout. In (11), this yields G = S = 0. And the final coating thickness would therefore be:

T + = Q c = 1 2 | S c | . ( 18
)
This falls in line with the findings of [START_REF] Hocking | Deformations during jet-stripping in the galvanizing process[END_REF], where the end result of a detailed calculation yields a ratio T c /T + 2. Our estimates of the coating thickness will be given below (see Section 4.3), as well as summarized for the industrial parameters in Table 3.

Computational methods

In the research presented here we have applied the "Basilisk" computational code [START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations[END_REF], which is an in-house, GPL-licensed code whose main developer is one of the present authors (SP). It is a descendant of the "Gerris" code [START_REF] Popinet | An accurate adaptive solver for surface-tension driven interfacial flows[END_REF] and as the latter, it allows for the local adaptive mesh refinement (AMR) [START_REF] Puckett | A 3D adaptive mesh refinement algorithm for multimaterial gas dynamics[END_REF] using the quad/oct-tree type mesh -regular, structured cubic meshes without refinement are also possible. The code is optimised for speed and capable of both OpenMP (single node) and MPI (multi-node) parallelism. Most recently, Basilisk has been applied e.g. to model compressible bubble dynamics [START_REF] Fuster | An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF], propose certain single-column models in meteorological simulations [START_REF] Van Hooft | Adaptive cartesian meshes for atmospheric single-column models: a study using basilisk 18-02-16[END_REF], or simulate turbidity currents [START_REF] Yang | A two-dimensional layer-averaged numerical model for turbidity currents[END_REF]. It is a multi-equation solver, making it de-facto a multi-physics code. Basilisk supplies a macro language (the Basilisk C) built atop the C99 standard of the C programming language. Upon compilation of the simulation case files, they are first parsed into clean C by the builtin parser/lexer (qcc). Then, the compilation is finalized with libc as the only dependence. This approach greatly enhances code portability.

Navier-Stokes equations are solved using the projection method [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF] with a procedure similar to that applied in Gerris [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension driven interfacial flows[END_REF]. Centered discretization is used for all scalar and vector fields, with additional face-centered values defined for u which are used e.g. to ensure divergence-free condition during mesh refinement. For consistency reasons, advection term of ( 1) is defined and calculated on cell faces as is ∇p. Advective fluxes are obtained using the Bell-Collela-Glaz scheme [START_REF] Bell | A second-order projection method for the incompressible navier-stokes equations[END_REF]. All discretizations are finite-differencing up to second order, unless noted otherwise. The Runge-Kutta scheme is used for time advancement. Optimisations of Poisson equation's solution are achieved by implementing the Multigrid (MG) method [START_REF] Brandt | Multigrid Techniques: 1984 guide with applications to fluid dynamics[END_REF].

The Volume of Fluid (VOF) method [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF] is used to track the interface using geometric interface reconstruction [START_REF] Aniszewski | Volume of fluid (VOF) type advection methods in two-phase flow: A comparative study[END_REF]. In this method, fraction function C (equal to one or zero in either phases ) is passively advected with the flow. Grid cells with fractional C 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 C, obtained by averaging in the interfacial cells. Interface curvature is also computed from C, using the Height-Functions method [START_REF] Afkhami | Height functions for applying contact angles to 3D VOF simulations[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension driven interfacial flows[END_REF] with a proper treatment close to solution boundaries. Depending on the local resolution, κ may be obtained by the HF method using full stencils (optimal) or one of two fallback options: parabola fitting of height functions and parabola fitting of the VOF centroids.

Basilisk's procedure for local mesh adaptation employs the wavelet transform of a given scalar field to assess the latter's discretization error. If the error is above the user-specified threshold, the grid is locally refined by subdividing it onto four (quad-tree) and eight (octree) sub-cells and performing a prolongation of the courser-mesh scalar onto children cells to obtain their initial values (the inverse process is termed restriction). For the simulations presented herein, we use u and C fields' error as the refinement criteria with 1 • 10 -3 and 1 • 10 -2 error thresholds, respectively.

Ensuring Momentum Conservation in Two-Phase Flow

The momentum-conserving methods [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. application to atomization process[END_REF] derive from a variant of the VOF [START_REF] Hirth | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] method originally proposed in [START_REF] Rudman | A volume-tracking method for incompressible multifluid flows with large density variations[END_REF] to treat the twophase flows with considerable density ratios. The idea is that instead of the simple incompressibility assumption

∇ • u = 0, (19)
one writes the mass transport equation in full, as is done in compressible formulation [START_REF] Pilliod | A second-order unsplit method for modeling flames in two-dimensional compressible flow[END_REF]:

∂ρ ∂t + ∇ • (ρu) = 0, (20)
using also the conservative form of the Navier-Stokes equations (not shown) [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. application to atomization process[END_REF], which contain the momentum term

∇ • (ρu ⊗ u). (21)
Subsequently, in implementation, we calculate density from the fraction function definition:

ρ = ρ l C + (1 -C)ρ g , (22)
which is an reversal of the traditional approach [START_REF] Hirth | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF]. Expression ( 22) is nontrivial only in the interface cells. The way in which the momentum-conserving methods differ from the traditional two-phase Navier-Stokes equation models using VOF is that subsequently, the ρ(C)U products found in both ( 20) and ( 21) are calculated consistently in the same control volumes. This can be non-trivial if staggered grid discretizations are used, and can be solved either by gridcell subdivisions [START_REF] Rudman | A volume-tracking method for incompressible multifluid flows with large density variations[END_REF] or using sub-fluxes of fraction function [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. application to atomization process[END_REF]. Thus consistency between transports of mass and momentum are ensured numerically, resulting in a far more robust computation.

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 large-scale systems involving solids [START_REF] Lin-Lin | Three-dimensional numerical simulation of a bird model in unsteady flight[END_REF]. It is a strongly simplified variant of the Immersed Boundary Method (IBM) of Peskin [START_REF] Peskin | The immersed boundary method[END_REF], 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:

∀x ∈ Ω ∪ Γ : u(x) = 0, (23)
that is, all velocity components are set to zero within the solid. As long as no provisions are needed for x ∈ Γ, the crude approximation provided by ( 23) yields satisfactory results [START_REF] Lin-Lin | Three-dimensional numerical simulation of a bird model in unsteady flight[END_REF]. A moving wall can be prescribed by using a non-zero (u wall ) right-hand side in (23). Note however, that pressure p 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 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 [START_REF] Gibou | Efficient symmetric positive definite second-order accurate monolithic solver for fluid/solid interactions[END_REF]. 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 G 1 setup -the domain interior is merely a sum of cuboids: it contains no inclined nor curved surfaces. The no-slip condition at the surface of the substrate wall moving with velocity u wall can be reasonably approximated using (23).

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/oct-tree 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 0-level (parent/root) cell with four/eight 1-level sub-cells and so on. If the refinement criterion is locally fulfilled, Basilisk will keep refining the grid until it reaches 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. The 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 sub-areas while outside of them lower maximum level is forced. This tactic of refinement situates the presented simulation between the block-based [START_REF] Lakehal | LEIS for the prediction of turbulent multifluid flows applied to thermal-hydraulics applications[END_REF] and point-based [START_REF] Popinet | An accurate adaptive solver for surface-tension driven interfacial flows[END_REF] mesh refinement. Its drawback is the increase in the globally calculated numerical dissipation.

Results

Simulation Planning

The full, three-dimensional airknife configuration poses numerous challenges for reasons of code stability, CPU cost or the wide range of simulated physical scales. Due to this challenging character, we have tackled the case progressively, including the following steps.

We have commenced with the film formation studies, both in two and three dimensions, results of which will be presented in Section 4.2. this lets us compare the obtained thickness with the analytical prediction, as well as shed some light on the turbulent, or transitory, character of the film formation itself.

Subsequently, we include a study of the gas dynamics, namely of the impinging air jet without including the liquid.

Here, the goal is to validate the single-phase Basilisk result against the known approximation of the u profile. Results can be found in Section 4.3.

Once the two above phenomena are investigated, we present, in Section 4.4 the full, three-dimensional configuration results. This Section contains additionally a relaxed/industrial parameters comparison for G 1 . Subsequently, the G 2 results are presented in Section 4.5, featuring the discussion on varying injection velocity u in j in this configuration.

Film formation studies

As mentioned in Section 2.2, the dimensional zero-flux film thickness can be readily obtained by √ 3h 0 leading to

h 00 = 3µ l u wall ρ l g . ( 24
)
Using ( 24) we arrive at h 00 = 5.46 • 10 -4 µm for the industrial parameters. With this value, we can make a rough attempt at assessing the Reynolds number characterizing the liquid film formation in the industrial process, that is

Re f (h 00 ) = ρ l h 00 u w µ l = 2240. ( 25
)
This value is close to a critical Re f delineating the laminar and turbulent film formation regimes (see also Table 2). Another estimate is possible if the h G thickness is used. Note that its dimensionless counterpart T G was mentioned above in the context of (11). For the high-Re regime, [START_REF] Groenveld | Laminar withdrawal with appreciable inertial forces[END_REF] proposes T = 0.52 and Q = 0.47 which results in h G = 1.63 • 10 -4 µm.

(26)

To get (26), one uses ( 7) and ( 8) and returns to dimensional quantities through a multiplication by h 0 . Once h G is used to calculate the Reynolds number, a value of Re f = 672 is found, so three times smaller than Re f (h 00 ). To complete the picture we denote that the definition of Re f introduced by (Colina-Marquez et al., 2016) results in Re f = 2600. To summarize, these estimates of the Reynolds number suggest a possibility of a turbulent character of the film formation process, as the numbers are between second and third order. One could conclude that the coating process used in the industry is a man-made system at the edge of criticality. This is confirmed in the results presented below. We start our investigation or the film formation regime with the two-dimensional case. Figure 3 presents the interface geometry at six instances of time where t ∈ [0, 0.14]. For this 2D simulation, the u wall velocity is implemented as a boundary condition; the simulation used the 2 14 -equivalent grid. Fig. 3b presents the same interfaces, but using log x for the horizontal axis. First off, it is visible from Fig. 3 that the uptake of the film by the moving wall is somewhat faster than expected from the value of u wall . We attribute this to imposed contact angle ofthe VOF interface [START_REF] Afkhami | Transition in a numerical model of contact line dynamics and forced dewetting[END_REF]. Wavy character of the film becomes visible towards the t = 0.14s mark, especially visible in Fig. 3b. Groenveld's thickness prediction h G is marked with a dashed line in both subfigures, a good overall fit is observed at least in the upper parts of the film. Apparent bulges on the film for t > 7.8 • 10 -2 s are consistent with a transitory or turbulent regime, we will examine this effect more closely in three dimensions below.

At this grid resolution, the Basilisk grid cells are of the size ∆x = L x /2 14 ≈ 39µm for the 2D simulation depicted in Fig- Figure 4 presents the creation of a boundary layer in time (in the t > 0.1457s range) of the same flow. The profiles have been sampled at h = 0.14m or 0.04m above the reservoir. The velocity profile remains parabolic, however it clearly becomes steeper for t > 0.1s with an apparent plateau extending for x > 5 • 10 -4 suggesting a detachment of the layer adjacent to the plate [START_REF] Snoeijer | Thick films of viscous fluid coating a plate withdrawn from liquid reservoir[END_REF]. In Fig. 4 we additionally compare the profile for t = 0.1457s with analytical expression (5) (dots). Consistency is visible especially very close to the wall, suggesting that the final profiles lend themselves well to those assumed in [START_REF] Groenveld | Laminar withdrawal with appreciable inertial forces[END_REF], as hinted previously by Figure 3. Additionally, it is observable in Fig. 4 (inset) that the transition through the profile described by (5) occurs for t ∈ [7.85 • 10 -2 s, 1.45 • 10 -1 s] which coincides with the moment the boundary layer starts forming. This serves as an argument that the film evolution is reasonably well described by the high-Re theory. Moreover, in Fig. 4 profiles are sampled only for C > 0 (i.e. inside the liquid film). Thus, for each of the lines, the abscissa of its right-hand end-point corresponds to the film thickness h(y) at y = 0.14m. As one can observe e.g. for t ≈ 7.8 • 10 -2 s we have h(0.14, t) ≈ 6 • 10 -3 m whereas for t = 1.457 • 10 -1 s the thickness drops, suggesting a bulge has passed over the point and retracted.

We now shift our attention to three-dimensional cases. Using the 2 12 3 -equivalent grid, we have performed a threedimensional simulation G 2,12 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, t = 1.45 • 10 -1 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 com- ponent u y . As we can observe -from the u y values the film is colored with -distinct liquid boundary layer develops directly adjacent to the wall, traveling with velocity u wall . This is fully consistent with liquid velocity profiles presented in Fig. 4 for t = 1.45 • 10 -1 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 (Fig. 5b) an isosurface for the zero vertical velocity u y = 0 , rendered in turquoise against the gray interfacial surface. (Note that u y = 0 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 t = 1.48 • 10 -1 s as 0.13m e.g. 0.03m 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 0.01m (or eighty times the grid size at 12 levels of grid refinement).

Another simulation is presented, for the G 1 configuration, in Figure 6. Even using a slightly less refined grid (11 levels of refinement, or 2048 3 -equivalent), we still observe a wrinkling of the interface as well, mostly on the coated band edge. In this stage of the flow, 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 u wall = 2. The surface of the film undergoes progressive distortion starting from the side of coated band. This applies especially to the coated x+ and x-walls, in which wrinkling ap- pears progressively further from the band edges. The turbulent nature of the film is suggested by surface disturbances, along with fully three-dimensional character of the wrinkles/waves. To our knowledge, this is the first published result of a 3D coating simulation including the edge, and Re f is far higher than the previously published 2D results [START_REF] Lacanette | Macroscopic analysis of gas-jet wiping: Numerical simulation and experimental approach[END_REF][START_REF] Myrillas | Numerical modelling of gas-jet wiping process[END_REF].

We continue our examination of the physics of threedimensional film formation with Figure 7, which contains velocity profiles for the vertical component (u y ) and the transverse component (u z ) along the wall height -only the height range of y ∈ [0, 0.2] is included, as all t values included are smaller than t = 0.15s. At that time, the liquid reaches roughly to y = 0.3m, consistent with Fig. 3. Four instantaneous profiles are presented with separate point types. Each of the profiles in Figure 7 has been z-averaged so they represent information from the entire width of the coated plate. Also, profiles include data only for x < 0.001m (across the film), in other words the measurement window includes only the direct proximity of the coated wall. In Fig. 7a, we observe a transition from a rather smooth u y profile at t = 1.5 • 10 -2 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 u y < 0 velocities are present in this region further from the wall -this however has not been captured with the profile measurement window. Average u y z 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 u y z at t = 0.15s (red color in Fig. 7a). This curve, although calculated using a three-dimensional simulation, is comparable with Fig. 4 (curve for t = 0.1457s). If, using the latter of the mentioned curves, one calculates a mean value (for x ∈ [0, 1]) of u y , it is equal to 1.22 m/s. This value should be at least comparable with Fig. 7a taken for t = 0.15s and y = 0.14m; in fact, we find u y z (0.14) ≈ 1.1 which is within ten percent of the two-dimensional simulation. The slight discrepancy might be attributed to the z-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 three-dimensional character of the film is found in Figure 7b, showcasing this time the profiles of the transverse velocity component u z . While close to the beginning of the flow at t = 0.015s (blue squares) this component is nearly zero (approximately two-dimensional flow), u z 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 z = 0). 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 possible that at this Re f values, three-dimensional 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 2 13 which is (locally) equivalent to a grid with 8192 3 points 1 .

Single-phase Impinging Jet Study

This subsection presents results of a 2D numerical study of the G 2 configuration (Figure 2). However, for the test presented here, the liquid phase (along with the reservoir) is removed, and we focus on the gas phase. Calculations presented here have been motivated by the need to find coefficients c s and c p in ( 14), which such a simplified (and cheaper) configuration allows. Additionally, it provides validation for the Basilisk code as shown below.

The simulation used the "industrial" parameters and a 2 11equivalent grid. Results are presented in Figure 8. As we can see from Fig. 8b, the airflow develops symmetrically with the nozzle in its axis. The turbulent character of the flow is not immediately visible in the Figure due to the application of temporal averaging -the technique was used only in this simulation, as discussed below.

Figure 8a presents the profile of the velocity component u y normalized by the mean value u m , taken at 3 • d n f /4 distance from the stagnation point. To plot this profile, we have applied a combination of time-and ensemble-averaging in order to ensure smoothness (15 simulations were ensemble-averaged). Time-spans used for temporal averages were equal to the time needed for largest vortices to leave the flow domain. The curve is accompanied by a fit to the analytic prediction presented by [START_REF] Ozdemir | Impingement of an axisymmetric jet on unheated and heated flat plates[END_REF], namely:

u u max = γ β x/x 0.5 β γ-1 • exp - x/x 0.5 β γ , ( 27
)
where we assume u max = 200m/s as per problem specification, while coefficients γ and β are taken 1.32 and 0.73 respectively. Value of x 0.5 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 0.75mm. A good fit between the simulation and the analytic curve is observed.

Figure 8b 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 ( 27) and the existence of boundary layer. According to [START_REF] Tuu | Wall pressure and shear stress measurements beneath an impinging jet[END_REF], peak pressure evolves with the distance from the nozzle as p s

1 2 ρ g |u| 2 ≈ 7d d n f -1 . ( 28
)
Values presented in Fig. 8b are about half of predicted by ( 28), meaning that the potential core is not resolved well enough for G 2,11 . This warrants an increase in refinement, and is one of the reasons for using at least 12 levels of refinement in majority of the simulations. By fitting the simplified p air and τ yy curves, constructed using ( 14) to the results obtained in this subsection, we were able to establish the values of the c s and c p coefficients mentioned in Section 2.2. For the industrial parameters, they amount to c p = 0.09, c s = 0.00325. These values are plugged into (18) through ( 15), and the resulting T + value is then made dimensional as h + = T + h 0 . We thus arrive at the dimensional coating thickness h + value roughly between 20 and 50 microns. In this section, results for the full G 1 configuration are discussed. In most cases, the simulations ran at the 2 12 -equivalent resolution, i.e. ∆x ≈ 100µm. Thus the grid size is below both h 00 and h G , but above the expected h + values. However, due to specifics of the geometric VOF method [START_REF] Aniszewski | Volume of fluid (VOF) type advection methods in two-phase flow: A comparative study[END_REF], it is possible to represent a film with thickness below ∆x, including advection of such a film with the uniform velocity defined in that cell2 .

The G 1 results within this Section will include variations in boundary conditions, which we refer to as industrial and relaxed parametrers. For reference, Table 2 contains parameters for both industrial and relaxed parametrisations of the considered problem. Most important differences between them include an order of magnitude lower liquid density and higher u w in the relaxed set: 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 u in j is decreased twofold in relaxed parametrisation.) This results for example in the zero-flux h 00 thickness of the film in relaxed parameters being fourteen times that of its value in industrial parameters. Additional difference between the relaxed and industrial configurations is the coated plate thickness, it is held at 5mm for the relaxed variant and 1mm 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 (23) formulation as more grid-points end up contained in the Ω region.

We begin with a discussion of the G 1,12 simulation, which is first presented in Figure 9 displaying the interface geometry at t = 0.162s. This is the stage of the flow at which the film has already been formed on the coated band, and right after the airflow, issuing from the planar nozzles, impacts it. For comparative purposes, the flow for the relaxed configuration is shown in Fig. 9b. For this image, the value of t is obviously different, and the images have been matched via the width of the impact zone. In Figure 9, nozzle locations are drawn using shading and black outlines in this view -this is done in post-processing and only for orientation purposes. Additionally, a cut-plane is positioned in the back-drop (parallel to z = 0 coordinate) colored by vorticity.

For the industrial configuration, at t ≈ 0.16, 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 three-dimensional: 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 emanating from the 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 [START_REF] Myrillas | Numerical modelling of gas-jet wiping process[END_REF].

The obvious difference between the relaxed and industrial parametrisations for the wiping process is the degree of visible atomization. Comparing the adimensional numbers as seen in Table 2, the airflow Re is three times higher in the industrial case (based on u in j and d). The effects are demonstrated in Figure 10, which displays the relaxed case in the earlier (than in Figure 9) stages of the macroscopic flow evolution3 . Significant atomization level is visible compared to Fig. 9a. Continuing the flow analysis for the industrial case, we turn our attention to Figure 11 which displays film geometry at approximately 0.17s. By this time, the lower bulge (A) starts forming (below the nozzle level) leading to the onset of backflow into the reservoir. The distance between the points A and B in Fig. 11 can be referred to as the impact zone, at t ≈ .017s we estimate its width at 0.04m. The wavy edge film structure visible in Fig. 11 is an early effect film formation. However, inspecting Figure 11b we may conclude that the edge film is nearly entirely atomized in the impact zone. In Fig. 11b the same film is seen, looking parallel to the z 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.

The total CPU cost of the 3D, 2 12 simulation of the G 1 configuration is approximated at 122000 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. The G 1 configuration includes the both airknives and the edge of the coated band -which, potentially, can produce an abundance of information about the physics of the flow. However, its CPU cost is nearly prohibitive if longer physical times were to be simulated. Thus, we include a more detailed analysis using the G 2 configuration, results of which will be presented below.

The G 2 configuration results

In this section, we present results pertaining to the G 2 configuration. It enables us to focus on the air-liquid impact study in more detail as, as long as the mean flow is considered, the G 1 configuration has an inherent symmetry. In G 2 , 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. Thus, we simulate only a single nozzle and one side of the coated band, decreasing the CPU cost. All other assumptions are carried over from the G 1 configuration.

In Figure 12, we present a visualisation of the macroscopic shape of the interface for the G 2 simulation performed using a 2 12 -equivalent grid. This simulation corresponds to industrial parameters (and is the G 2 -analog of the G 1 results mentioned above e.g. in Fig. 11). In three sub-figures, instantaneous shapes are visible for (a) t = 1.656 • 10 -1 s (b) t = 1.677 • 10 -1 s and (c) t = 1.747 • 10 -1 s. Each of the pictures presents two separate view: an isometric one on the left-hand-side, and a side-view (looking along z axis) on the right-hand-side. The cut-plane positioned at the z-domain wall is colored with ω. Varying cell size in the vorticity cut-plane is, of course, a consequence of employing the SRR technique to limit adaptivity in regions above and below the nozzle. Full resolution is maintained everywhere inside the planar nozzle and within 1 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 the liquid in Fig. 12a, we note a distinct imprint of the nozzle shape on it. 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 [START_REF] Yu | Numerical simulation of liquid-layer breakup on a moving wall due to an impinging jet[END_REF], who have investigated (in 2D) a flow characterized by a higher We of 13.5 with lower density ratios. Their results show a C 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 andc. In the meantime, the lower and upper bulges move away from the nozzle. In Fig. 12c, we note that the upper bulge has, by t = 1.747 • 10 -1 s advanced approx. 10mm 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 right-handside 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 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 three-dimensional wave structure as visible in Fig. 12c.

Atomization of the film occurring at the first instance the air-liquid contact might be investigated looking at the nondimensional numbers characterizing this interaction. While the film Reynolds number Re f characterizes mostly film formation, we formulate the Weber number We involving gas velocity, as follows:

We = ρ g u 2 g h σ , (29)
with h standing for film thickness. Definition ( 29) is first applied to industrial parameters characterized by u in j = 200. Feeding the zero-flux thickness ( 24) into (29) one obtains We (h 00 ) = 38.1. If, instead, we settle upon using Groenveld's thickness h G , (29) yields We = 11.4. Both these values seem consistent with a regime in which atomization might be expected4 . 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 [START_REF] Ellen | An analysis of jet stripping of liquid coatings[END_REF][START_REF] Myrillas | Numerical modelling of gas-jet wiping process[END_REF] 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 [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. application to atomization process[END_REF], curvature calculation errors, or not accounting for interactions between liquid-gas interface and vortical structures in the latter phase [START_REF] Aniszewski | Improvements, testing and development of the ADM-τ sub-grid surface tension model for two-phase LES[END_REF]. Thus, we have included simulation configured as G 2,12 (u in j = 42) by decreasing air injection velocity. This configuration is characterized by We (h 00 ) = 1.68 and We (h G ) = 0.5 which, again, places the system just below the "edge of criticality" as in a context of Re f in our film formation study. We follow up with an examination of the flow at a decreased injection velocity.

Figure 13 presents the G 2,12 (u in j = 42) simulation roughly 2 milliseconds after air-liquid impact. In the Figure, white shaded isosurface represents the liquid interface, while several blue areas depict the (un-coated) moving wall5 . 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 air-liquid interaction. The impact zone is visible above as an area with horizontal wrinkles (marked "2" Fig. 13a). 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 right-hand-side images in Fig. 12 we note the nearly complete lack of atomized structures for u in j = 42m/s and We ≈ 1.

Another look at the low Weber number flow is provided by Figure 14. Three sub-figures present the same flow as pictured in Figure 13 Overall image of the flow is different than that for We > 10. 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-We 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 z-averaged film thickness in the vicinity of the impact zone (the nozzle level is marked with an arrow around y = 0.3m) has been altered but not significantly diminished, except the area some 15mm above the nozzle. In the Figure, h(y) profiles are shown for two instances: t = 1.687 • 10 -1 (continuous line) and t = 1.783 • 10 -1 (black dots). Groenveld's thickness h G is drawn in dashed line for comparison. Profiles visible in Fig. 15 show clearly a bulge for y ∈ [0.265, 0.295]. 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 [START_REF] Snoeijer | Thick films of viscous fluid coating a plate withdrawn from liquid reservoir[END_REF] -visible in Fig. 14a, also visible in 2D in Fig. 3 above y = 0.3m -and the coat is of three-dimensional character, being on average thicker closer to the z+ edge. This is reflected in the profiles below the impact zone. As for the consequences of the impact itself, h(y) oscillates in the vicinity of y = 0.3 which is fully consistent with instantaneous images in Fig. 14a-c and indicates alternating areas of thinning and thickening of the film. (Note that the zero-level shift in Figure 15 is a correction for the wall thickness of 5 • 10 -4 m.)

By representing p air (y) and τ yy using gas velocity as in ( 14), we note their magnitudes for u in j = 42 are one order below that for u in j = 200, which in the context of ( 14) amounts to a higher h + . 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 post-impact film thickness h + (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 y ∈ [0.31, 0.32], 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 y ≈ 0.314 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 the droplet volume distribution for G 2,12 (u in j = 200) (pink bars) and G 2,12 (u in j = 42) (yellow bars). In the Figure, minimum cell volume ((∆x) 3 at the finest grid level) is denoted with a black vertical line. Clearly, both simulations involve a significant number of 'sub-grid' VOF "debris"grid-cells containing non-zero 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 15 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 the liquid velocity obtained using the VOF fraction function C, which is equal to 1 inside the liquid. In other words the product C(x)u y disappears in the gas phase, and Fig. 17 shows its profile in the direct neighbourhood of the impact zone (y ∈ [0.28, 0.32]) i.e. two centimeters below and above the impact zone). Two simulations are included with u in j = 42 and 200 m/s. Since two flows have slightly different characteristic time scales due to higher We in the later, we have compared instantaneous profiles at the instance corresponding to the impact zone width approximately equal 0.01m (as e.g. in Fig. 14c). The curves have been averaged over a sampling window 1 mm thick. For both flows, pre-impact velocity distribution in the analysing window is similar with u y ≈ 1, 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-We 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 u in j the average u y values remain positive, suggesting the air knife wiping is far weaker for chosen injection parameters.

This concludes our investigation of the influence of the Weber number on atomization process -we conclude that the atomization effect visible at higher We is a correct result. The simulated air-liquid system responds as expected to the decrease in dimensional number, while other simulation parameters (e.g. grid resolution) are kept constant.

x We finish our analysis of the G 2 simulations with a brief remark on the computational efficiency. Thus, the approximated computational cost for the G 2,12 Basilisk simulations presented e.g. in Figure 16 was 4.67 • 10 5 CPU-hours (for simulation with u in j = 42m/s) and 4.8 • 10 5 CPU-hours (for u in j = 200m/s simulation).

Conclusions

In this paper, we have presented a novel set of simulations of a very demanding, two-phase fluid flow whose characteristics closely correspond to that of air and liquid Zinc. The boundary conditions correspond to the air knife jet-wiping process in hot-dip coating. In many aspects this is a pioneering work: to our knowledge, the only similar calculations published have described a two-dimensional case with RANS/LES performed for the airflow [START_REF] Myrillas | Numerical modelling of gas-jet wiping process[END_REF] or investigated the film formation only [START_REF] Snoeijer | Thick films of viscous fluid coating a plate withdrawn from liquid reservoir[END_REF] or, possibly, included a predefined film (i.e. not formed gravitationally). Similarly, for reasons of numerical stability [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. application to atomization process[END_REF], 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 substantiated e.g. by considering the Hinze scale l H , defined as the ratio of turbulent kinetic energy and surface tension, and estimated for the industrial parameters at l H ≈ 1.76 • 10 -3 m. At the resolution used in most cases presented here (12 levels of refinement) we obtain l H /∆x ≈ 14, proving that l H 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 Re in the air, and a properly defined Weber number. This gradient restricts the liquid feeding from reservoir, thinning the deposit. Our calculations possibly fall short of resolving the upper film thickness h + in the full G 1 and G 2 geometries due to lack of grid resolution. However, the G 2 case clearly displays the thinning effect, it is also observable in the G 1 case performed with decreased ρ l .

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 We ≈ 38, while for We values closer to one, liquid wrinkling is observed. 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 liquid-liquid collisions are aplenty. This is already visible at the 2 12 level (Figures 9b and 11), and suggest that the coat-thinning mechanism is far more turbulent in its nature than known previously.

The two geometries introduced in the paper focus on the two-and three-dimensional coat formation and nozzle interactions (G 1 ), coat thinning and edge effects (G 2 ) and the character of gas-liquid impact. For future research we would envision working preferable with the G 2 configurations in two and three dimensions, with increasing resolutions (and simulated timespans), preferably until full resolution of the h + thickness is feasible.

Disclaimer

On the 9th of August 2019, that is while this paper was already in its review stages, a potentially serious bug was reported in the Basilisk solver by the user Petr Karnakov 7 . The bug caused a redundant multiplication of the surface tension force by the σ coefficient (i.e. σ 2 instead of σ), which was activated only in the areas of strongly under-resolved interface. The wellresolved areas, using the HF method with full sized stencils, were not impacted, nor the secondary backup solution using parabola fitting of height functions. Only the tertiary backup method, i.e. parabola fitting of VOF centroids introduced the additional multiplication. The bug was present in the Basilisk code repository between December 13, 2017, and August 19, 2019.

7 https://groups.google.com/d/msg/basiliskfr/8Ub0LKQWyDc/NtBsk7p2EwAJ

In this paper, the simulations of 2D and 3D film formation presented in Section 4.2 are not affected by the bug, as they were performed before Dec 13, 2017. However, the buggy surface tension calculation might have influenced results of the Section 4.4 in which we study the air wiping process.

To verify that, we have substantially elongated the redaction process by rerunning part of the simulations, and assessing the result obtained with the bug-free code. In particular, the G 2,12 (u in j = 200) simulation presented in Fig. 12 and others, has been repeated using a bug-free version of Basilisk. The results of the repeated simulations have proven to be statistically and qualitatively identical to those already presented. We thus stand by the validity of the results presented in the paper, including the wiping process simulations discussed above.

Figure 1 :

 1 Figure 1: The coating configuration in two (left, half of the geometry visible) and three (right) dimensions. A -upward moving band; B -the air-knives or flat jet nozzles; C -liquid zinc containers. Note that outer domain walls are invisible in 3D rendering.

Figure 2 :

 2 Figure 2: Schematic renderings of the G 2 configuration. Outer domain boundaries are not visible, nozzles are visible in black.

Figure 3 :

 3 Figure 3: Configuration G 2,14 (no air injection). Interface geometry at chosen t values (with x in (a) dec and (b) log). The dashed line is h G = 163µ m.

ure 3 .

 3 This translates to four grid points inside the film with thickness h G and allows us to get a very rough estimation of the velocity profile inside the forming film, as discussed in context of the next Figure.

Figure 4 :

 4 Figure 4: Configuration G 2,14 (two-dimensional, no air injection); u y (x) profiles through the film at varying t values taken from Fig. 3 (lines), Groenveld's prediction (5) (points).

Figure 5 :

 5 Figure 5: Configuration G 2,12 film formation study; the flow at t = 0.148s. (a) Actual VOF-reconstructed liquid-gas interface geometry, colored by the u y velocity component. Inset (b): liquid-air interface shown in gray with the u y = 0 isosurface drawn in turquoise to approximately delimit the stagnation height. Scale bar corresponds to 0.05m in both sub-figures.

Figure 6 :

 6 Figure 6: Coating in G 1,11 configuration (no air injection) at t ≈ 0.32s. Interface colored by the vertical velocity component.

Figure 7 :

 7 Figure 7: Coating in G 2,13 configuration: z-averaged velocity profiles for x ∈ [0, 0.001] at varying t values. (a) u y (vertical) velocity component; (b) u z (transverse) velocity component. The shaded plot areas correspond to those parts of the sampling window completely below the bath level.

  displays the spatial distribution of mean pressure p normalized by the dynamic flow pressure, or p = p 1 2 ρ g |u| 2 1 Due to CPU time and memory restrictions we have not continued this simulation into the air injection stages.

Figure 8 :

 8 Figure 8: Study of an impinging jet flow using two-dimensional G 2,11 ): (a) Velocity profile near the wall, simulation (ensemble averaged, blue) and (Ozdemir and Whitelaw, 1992) ((27) brown); (b) Ensemble-averaged mean pressure in the same simulation; flow is from left to right, the nozzle is visible as a two rectangles close to the coordinate origin.

  4.4. Three-dimensional Wiping Simulations. Configuration G 1 .

  Re (u in j , d) Re f (h 00 , u wall ) Unit (kg/m 3 ) (kg/m 3 ) (Pa•sParameters for the discussed simulations in industrial and "relaxed" variants).

Figure 9 :

 9 Figure 9: (a) G 1,12 with nozzle locations sketched as a shaded area. Color: ω 0. (b) The G 1,12 (u in j = 75, u wall = 4) simulation in its final stages.

Figure 10 :

 10 Figure 10: The G 1,12 u in j = 75, u wall = 4 simulation ("relaxed" parameters).

Figure 11 :

 11 Figure 11: The G 1,12 simulation (industrial parameters) at t = 0.164s. Views: (a) isometric and (b) side view (along z axis). The back-drop cut-plane colored by vorticity.

Figure 12 :

 12 Figure 12: The G 2,12 u in j = 200 simulation at (a) t = 1.656 • 10 -1 s (b) t = 1.677 • 10 -1 s and (c) t = 1.747 • 10 -1 s.

Figure 13 :

 13 Figure 13: The G 2,12 u in j = 42 simulation at t = 1.756 • 10 -1 s. (a) Isometric view; (b) zoom into the impact region at the same time instant. Navy blue color indicates the wall (where coating is absent). 1 : bulges created in the high Re f withdrawal, 2 : air impact area, 3 : coating defect in the impact zone 2, 4 : coating defect above the impact zone.

  at instances of time with t ∈ [0.1736, 0.178]. Figure 14a presents interface geometry in the impact zone almost directly following the first air-liquid contact 6 . Effects of the decrease in the Weber number are instantly recognizable: no distinct horizontal liquid bulges are formed along the z direction; instead, smaller-wavelength 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.

Figure 14 :

 14 Figure 14: The G 2,12 u in j 42 simulation at: (a) t = 1.736 • 10 -1 s, (b) t = 1.752 • 10 -1 s,(c) t = 1.78 • 10 -1 s.

Figure 15 :

 15 Figure 15: Film thickness profiles for the G 2,12 (u in j 42) simulation.

Figure 16 :

 16 Figure 16: Droplet size distribution in the G 2 configuration for varying air injection density.

Figure 17 :

 17 Figure 17: Average liquid velocities (u y component) in the impact zone.

Table 3 :

 3 Values of film thickness and the resulting dimensionless numbers for the industrial air-knife configuration. In the rightmost column, the values for Re and We are calculated using the average value of h + . Thickness given in m.

		h 00	h G	h +
	Value	5.46 • 10 -4 1.63 • 10 -4 ≈3 • 10 -5
	Re f (•)	2.24 • 10 3	6.72 • 10 2	1.23 • 10 2
	We (•)	3.8 • 10 1	1.14 • 10 1	2.09

Naturally, higher resolutions are required to represent velocity variation within the film.

Note that in Figure10, the shape of the nozzles is distinguishable by looking at the cut-plane in the back.

For the "relaxed" configuration presented previously, using h 00 is more justified as the film Reynolds number is three times lower. Doing so, we obtain We (h 00 ) = 7.57.

In Fig.13a, the bath level liquid is over-exposed (i.e. rendered as white) due to specific light positions in the visualization.

Note that Figure13corresponds to the instance of time situated between Fig.14b and c.

See www.basilisk.fr/src/view.h

Acknowledgements

All the computations for this work have been performed using the French TGCC "Irene Joliot-Curie" and CINES "Occigen" supercomputers. Graphs and visualizations have been performed using the Basilisk View 8 , 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. Last but not least, we thank the anonymous reviewers of the paper for their valuable suggestions.

Wojciech Aniszewski would like to dedicate all his work on this paper to the dear memory of Antonina Gilska (nee Hessler, 1936(nee Hessler, -2019)).