Planar Jet Stripping of Liquid Coatings: Numerical Studies

07/17/2019 ∙ by Wojciech Aniszewski, et al. ∙ ProtonMail UPMC 0

In this paper, we present a detailed example of numerical study of flm 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 (tens 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.



There are no comments yet.


page 3

page 4

page 9

page 11

page 12

page 13

page 14

page 15

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 “air-knives”. 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 zero-flux and LL-type 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 film-air 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 one-fluid 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 two-dimensional 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 three-dimensional code with very high spatio-temporal resolutions and grid adaptivity.

1.2 Problem Specification

The investigated configuration is visible in Figure 1. As we can see in the side-view, the nozzle-band distance is measured at 10mm in industrial configuration. Nozzle diameter is mm. The proportions in the two-dimensional 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.

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.

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 left-hand-side is recognized as in 2D, while the r-h-s of Fig. 1 depicts

Conf. nozzles
Table 1: Distinguishing features of the , and initial conditions (in all dimensions in meters)).

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 ”air-knife” 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.

Figure 2: Schematic renderings of the and simulation configurations. Outer domain boundaries are not visible, nozzles are visible in black.

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 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 airknife-liquid 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 air-liquid interactions.

The third introduced configuration, is shown at the bottom of Figure 2. It is a ”synthetic” sub-problem, designed for a fully academic investigation of the air-metal impact phenomenon, and the initial stages of the two-phase flow post-impact. 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 nozzle-film gap and the coat. The coated wall is not present except as a boundary condition on the direction. In this configuration, film is pre-defined (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 in-nozzle flow as well as the gas-liquid impact.

2 Description of the Flow

2.1 Governing Equations

In all cases presented here, full Navier-Stokes equations:


are solved, assuming also incompressibility of the flow:


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 as

Body force (gravity) is taken into account and represented by We will occasionally refer to the directions “up” and “down” which in both two- and three-dimensional 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 one-fluid 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


Integrating (3) one obtains a parabolic profile inside the vertically moving film


with arbitrary Using this profile, one can define a dimensionless flux


and dimensionless film thickness


subsequently establishing a following dependency between the two:


Note that knowing the upward-moving wall velocity and liquid properties, finding is possible from (6) given that has been pre-computed 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 Liquid-air 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 air-liquid 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 condition111At this stage, can be seen as a order-correct prediction of film thickness, Hocking uses a value of one micrometer., and

for the pressure. This is supplemented by a transport condition for thickness


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 Navier-Stokes equations to:


In the above, is the Stokes number – in the cited work,

and is of minus-fourth 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:


and by (Elsaadawy et al., 2007) for the shear stress:


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 thin-film 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


We integrate (12) to get the flux equation:


Note that in (13) the zero-flux 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


where We rephrase the above result by introducing a pressure to shear strain ratio defined as

with this, (14) becomes


Further on, by approximating pressure and strain, for example by


we can represent as

with and being pressure- and shear strain-related 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 thin-layer-approximated thickness to shear-stress; the actual thickness approximations for using that formulation will be given below. In the general, non-laminar 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


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




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, post-impact 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 in-house, GPL-licensed 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/oct-tree 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 MPI-type (multi-node) operation. Most recently, Basilisk has been applied e.g. to model compressible flows connected to bubble dynamics (Fuster and Popinet, 2018), propose single-column 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 multi-equation solver, i.e. a broad framework that allows choosing equations to be solved, making it de-facto a multi-physics code. Secondly, it contains a built-in parser/lexer providing “targeted”, minimal re-compilations for the configuration currently used.

Navier-Stokes 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 face-centered values defined for

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 Advective fluxes are obtained using the Bell-Collela-Glaz advection scheme (Bell et al., 1989). Discretizations are generally finite-differencing up to second order unless noted otherwise. The Runge-Kutta 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 Height-Functions 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 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 and fields’ error as the refinement criteria with and error thresholds, respectively.

3.1 Ensuring Momentum Conservation in Two-Phase Flow

The momentum-conserving methods (Vaudor et al., 2017) derive from a variant of VOF (Hirth and Nichols, 1979) method originally proposed in (Rudman, 1998) to treat two-phase flows with considerable density ratios. General idea is that instead of a simple incompressibility assumption


we instead write the mass transport equation in full, as is done in compressible formulation (Pilliod, 1996):


using also the conservative form of the Navier-Stokes equations (not shown) (Vaudor et al., 2017), which contain the momentum term


Subsequently, in implementation, we calculate density from the fraction function definition:


instead of the other way around as it is done traditionally (Hirth and Nichols, 1979). The way in which the momentum-conserving methods stand out from traditional two-phase Navier-Stokes 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 non-trivial if staggered grid (Tyliszczak, 2014) discretizations are used, and can be solved either by grid-cell subdivisions (Rudman, 1998) or using sub-fluxes 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 large-scale systems involving solids (Lin-Lin 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:


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 (Lin-Lin et al., 2016). A moving wall can be prescribed by using a non-zero () right-hand 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 no-slip 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/oct-tree data structure in Basilisk results in subdivisions of cells into four/eight sub-cells 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. 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 sub-areas while outside of them lower maximum level is forced. This tactic of refinement situates the presented simulation between the block-based (Lakehal, 2010) and point-based (Popinet, 2009) mesh refinement.

4 Results

4.1 Planning of the Simulations

The full air-knife configuration poses numerous challenges for computational simulation. Firstly, it comprises of sub-processes – 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 CPU-cost (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 sub-problems. For reasons of brevity we will only present here a subset of the obtained results, namely:

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

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

  3. Studies of film formation and airknife-liquid interaction with ”relaxed” and industrial parameters.

  4. 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 zero-flux thickness of the film in relaxed parameters being fourteen times that of its value in industrial parameters.

Unit (kg/m) (kg/m) (Pas) (Pas) (m/s) (m) (m) (m/s)
Table 2: Parameters for the discussed simulations in industrial and ”relaxed” variants).

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 grid-points 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 air-knives (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


where is liquid density, is the upward-moving band (wall) velocity, and is a zero-flux 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


Using (26) and calculating from (25) we arrive at values of and for industrial parameters. Indeed, one could say simulations prove that industrial-class 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:


While this value is three times smaller than , one could expect the film to be at least in the intermittent regime.

Figure 3: Configuration (no air injection). Interface geometry at chosen values (with in (a) dec and (b) log). The dashed line is m.

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 upward-moving 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: Configuration (two-dimensional, no air injection). (a) profiles through the film at varying values taken from Fig. 3 (lines), Groenveld’s prediction using (4) (points).

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 right-hand end-point 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.

Figure 5: Configuration film formation study; the flow at s. (a) Actual VOF-reconstructed liquid-gas interface geometry, colored by the velocity component. Inset (b): liquid-air interface shown in gray with the isosurface drawn in turquoise to approximately delimit the stagnation height.

Using the -equivalent grid, we have performed a three-dimensional 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).

Figure 6: Coating in configuration (no air injection) at s. Interface colored by the vertical velocity component.

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 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 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).

Figure 7: Coating in configuration: averaged velocity profiles for at varying values. (a) (vertical) velocity component; (b) (transverse) velocity component.

We continue our look at the physics of the three-dimensional 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 three-dimensional 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 two-dimensional 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 three-dimensional 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, 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 which is locally equivalent to a grid of points 222Due to CPU time and memory restrictions we have not continued this simulation into the injection stages..

4.3 Single-phase Impinging Jet Study

A brief study has been performed on the velocity and pressure profiles obtained in boundary layers of a (two-dimensional), single-phase 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 air-wall 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 ensemble-averaging in order to ensure smoothness of profiles (in most situations 15 simulations were ensemble-averaged). Time-spans 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 8: Study of an impinging jet (single-phase flow using two-dimensional ): (a) Velocity profile near the wall, simulation (ensemble averaged, blue) and (Ozdemir and Whitelaw, 1992) ((28) brown); (b) Ensemble-averaged mean pressure in the same simulation.

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:


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


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 above-nozzle film thickness is thus between and for this configuration with thin layer assumptions. By using an iterative computation with the Colebrook-White 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 thin-layer 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 Three-dimensional Wiping Simulations. Configuration

In this subsection, we present the first results of the three-dimensional simulations performed. This includes geometry specification given in Table 1 and Figure 1a. This first presented simulation ran at refinement level; with , translating to grid-cell size of


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.

Figure 9: (a) with nozzle locations sketched as a shaded area. Color: (b) The simulation in its final stages.

To further examine the influence of the air-knives, we will proceed to three-dimensional, 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 post-processing and only for orientation purposes. Additionally, a cut-plane is positioned in the back-drop (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 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 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).

Figure 10: The simulation parameter study: (a) , (b) with industrial parameters.

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 zero-flux 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 ).

Figure 11: The simulation at s. Views: (a) isometric and (b) side view (along axis). The back-drop cut-plane colored by vorticity.

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 back-flow 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 (A-B 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 m-wide area above and below the air-knives.

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 air-liquid 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.

Figure 12: The simulation at (a) s (b) s and (c) s.

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 sub-figures, 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 left-hand-side, and a side-view (looking along axis) on the right-hand-side. The cut-plane positioned at the domain wall is colored with Varying cell size in the vorticity cut-plane 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 side-view.) 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 fast-moving “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 right-hand-side 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 three-dimensional 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 air-liquid contact might be investigated looking at the non-dimensional numbers characterizing this interaction. While the film Reynolds number characterizes mostly film formation, we formulate the Weber number involving gas velocity, as follows:


with standing for film thickness. Definition (31) is first applied to industrial parameters characterized by

Table 3: Values of film thickness and the resulting dimensionless numbers for the industrial air-knife configuration. Values for and are calculated using the average values of and .

Feeding the zero-flux 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 333For 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 liquid-gas 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: The simulation at 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). bulges created in the high withdrawal, air impact area, coating defect in the impact zone , coating defect above the impact zone.

Figure 13 presents the simulation roughly milliseconds after air-liquid impact. In the Figure, white shaded isosurface represents the liquid interface, while several blue areas depict the (un-coated) moving wall444In Fig.13a, the bath level is over-exposed (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 air-liquid 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 air-liquid contact, comparing Fig. 13 with right-hand-side images in Fig. 12 we note the nearly complete lack of atomized structures for m/s and

Figure 14: The simulation at: (a) s, (b) s,(c) s.

Another look at the flow characterized by lower Weber number is provided in Figure 14. Three sub-figures 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 air-liquid contact555Note 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, smaller-wavelength disturbances are showing within a gradually broadening region, as seen in Fig. 14c. The side-views 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.

Figure 15: Film thickness profiles for the simulation.

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 three-dimensional 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. 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 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 post-impact 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: Droplet size distribution in the configuration for varying air injection density.

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 ’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 resolved droplets in total (yellow bars) which is qualitatively different than at higher air velocity (red bars).

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

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, pre-impact 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 air-liquid 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 CPU-hours (for simulation with m/s) and CPU-hours (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 air-liquid 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: : (a) velocity component (b) component of shear at varying moments in time.

Figure 18 presents example profiles for the type simulation with In this simulation, gravitational coating of the wall was replaced with a pre-defined 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.

Figure 19: : Pressure profiles.

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 air-liquid 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 “air-only” 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 function-weighted (effectively density-weighted) , namely


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.

Figure 20: : Film thickness profiles for an example simulation: (a) sampled using entire span; (b) using only a narrow strip around

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 zero-flux thickness has been marked using a dot-dashed line.

In is important to note that in the configurations film thickness is pre-defined, 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 non-trivial; 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, two-phase fluid flow whose characteristics closely correspond to that of air and liquid metal (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 (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 gas-liquid 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 liquid-liquid collisions are aplenty. This is already visible at the level (Figures 9b and 11), and suggest that the coat-thinning mechanism is far more turbulent in its nature than known previously.

The three geometries introduced in the paper focus on the two-and three-dimensional coat formation and nozzle interactions (), coat thinning and edge effects () and the character of gas-liquid 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 time-spans), preferably until full resolution of the thickness is feasible.


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 View666See 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, 1936-2019).


  • 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.
  • 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 two-phase flow. Ph.D. thesis, Czestochowa University of Technology, Czestochowa, Poland.
  • Aniszewski (2016) Aniszewski, W., 2016. Improvements, testing and development of the ADM- sub-grid surface tension model for two-phase 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 two-phase 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.
  • Bell et al. (1989) Bell, J. B., Colella, P., Glaz, H. M., 1989. A second-order projection method for the incompressible navier-stokes 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.
  • Brandt (1984) Brandt, A., 1984. Multigrid Techniques: 1984 guide with applications to fluid dynamics. GMD-Studien 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 two-phase 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 hot-dip galvanizing process. Metallurgical and Materials Transactions B 38 (3), 413–424.
  • Fuster and Popinet (2018) Fuster, D., Popinet, S., 2018. An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension. Journal of Computational Physics 374, 752 – 768.
  • 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 second-order 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 jet-stripping 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 gas-jet 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 thermal-hydraulics applications. Nuclear Engineering and Design 240 (9), 2096 – 2106, experiments and CFD Code Applications to Nuclear Reactor Safety (XCFD4NRS).
  • Landau and Levich (1942) Landau, L., Levich, B., 1942. Landau-Levich film. Acta Physicochim. URSS 17.
  • Lin-Lin et al. (2016) Lin-Lin, Z., Hui, G., Chui-Jie, W., 2016. Three-dimensional 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.
  • 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 gas-jet 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 second-order unsplit method for modeling flames in two-dimensional compressible flow. Ph.D. thesis, University of California, Davis.
  • Popinet (2003) Popinet, S., 2003. Gerris: a tree-based 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 surface-tension driven interfacial flows. Journal of Computational Physics 228, 5838–5866.
  • Popinet (2015) Popinet, S., Dec. 2015. A quadtree-adaptive multigrid solver for the Serre–Green–Naghdi equations. Journal of Computational Physics 302, 336–358.
  • 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 volume-tracking 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 large-eddy simulation. Physics of Fluids 11 (7), 1699–1701.
  • Takeishi et al. (1995) Takeishi, Y., Yamauchi, A., Miyauchi, S., 1995. Gas wiping in hot-dip 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 Gas-Liquid 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 high-order compact difference algorithm for half-staggered 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 single-column models: a study using basilisk 18-02-16. Geoscientific Model Development 11 (12), 4727–4738.
  • 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.
  • Yang et al. (2018) Yang, S., An, Y., Liu, Q., 2018. A two-dimensional layer-averaged 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 liquid-layer breakup on a moving wall due to an impinging jet. In: Proceedings of the Summer Program 2014. Center of Turbulence Research.