Planar jet stripping of liquid coatings: Numerical studies☆
Introduction
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 the 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 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 (Landau and Levich, 1942). Analytic solutions were found e.g. by (Groenveld, 1970) who focused on withdrawal with “appreciable” inertial forces (relatively high Reynolds number (Re ) flows) or (Spiers et al., 1973) who have modified the withdrawal theory of Landau and Levich, obtaining improved predictions for the film thickness that were also confirmed experimentally. Later, (Snoeijer 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.
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 (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 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 (Hocking et al., 2011).
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 time averaging (Myrillas et al., 2013), certain instantaneous effects, such as the formation 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 (Ellen and Tu, 1984) who have shown analytically that not only is the pressure gradient acting on the film, but also the surface shear stress terms play 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 that only the pressure gradients play a role in the film deformation – and adopting the lubrication approximation for the film flow. The work (Takeishi et al., 1995) provided certain numerical solutions for velocity and shear profiles at the film-air interface during wiping (using a glycerine solution as the coating liquid).
The authors of (Hocking et al., 2011) 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 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 (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 adopted e.g. by (Lacanette et al., 2006), 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 (Myrillas et al., 2013) performed a study very similar to (Lacanette et al., 2006) – 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 (Yu et al., 2014), 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.
The investigated configuration is visible in Fig. 1. Dimensions visible in the leftmost illustration pertain to 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 Fig. 1 (left) with grayed lines.
As we can see in the side-view (Fig. 1 left), the nozzle-band distance dnf is measured at mm in the industrial configuration. The 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 the upward band velocity is in most cases taken at m/s. Except where noted, we have decided to choose liquid zinc as the coating liquid. Properties of 30Zn are assumed, that is surface tension density and viscosity Properties of the surrounding gas - which in all cases is air - are density ρa ≈ 1.22[kg/m3] 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. Two geometries considered will be termed Gi 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, G1,14 stands for the first configuration at 214-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 uinj) are varied between configurations, it will be designated in parenthesis (e.g. ) stands for the G2 configuration on a 211-equivalent grid with the air injection velocity equal to 42 m/s). Using the above terminology, we can now revisit Fig. 1: the configuration presented on the left-hand-side is recognized as G2 in 2D, while the r-h-s of Fig. 1 depicts the three-dimensional G1.
Our departure point is the full “industrial” configuration G1, visible in Fig. 1 on the right. As sketched in Fig. 1, we orient the geometry so that y is the vertical direction, and air injection takes place along the 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 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 Table 1, the thickness hw of the coated band A is kept at 1mm. The position of the coated wall along the x axis is given by xwall in the Table; it is centered in the G1 configuration, moved leftmost in G2. In all cases, we impose the upward wall velocity (m/s). Due to the fact that the extent (depth) of the coated wall is smaller than the nozzle depth, the G1 configuration allows the air issuing from both nozzles to collide. This ends the description of the G1 configuration.
An additional configuration is rendered in Fig. 2. As with Fig. 1, note that the rendering is not fully up-to-scale: dimensions used in actual simulations are given in Table 1. This G2 configuration has been created from G1 by including only half of it and a symmetry boundary condition at the direction. In other words, the G2 configuration is a three-dimensional realization of the sketch presented on the left-hand-side of Fig. 1. Going into G2, the depth (extent) of the coated wall has also been slightly decreased (from 15 to 5 cm) to limit computational cost of the simulation. Still, in the G2 configuration the film is formed gravitationally and the airknife-liquid interaction is maintained. Since the coated wall is now centered at only half of its thickness (x-span) is included in the G2 configuration, which makes G2 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. Of course, the G2 includes only a single nozzle.
Section snippets
Governing equations
In all the cases presented henceforth, the full Navier-Stokes equations: are solved, assuming the flow to be incompressible:
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 asGravity is taken into account and represented by the body force fg.
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 allows for the local adaptive mesh refinement (AMR) (Puckett and Saltzman, 1992) using the quad/oct-tree type mesh – regular, structured cubic meshes without refinement are also possible. The code is optimised for speed and
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, the results of which will be presented in Section 4.2. this lets us compare the obtained thickness with the analytical prediction, as well as
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 (Myrillas et al., 2013) or investigated the film
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 Karnakov7. 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 well-resolved areas, using the HF method with full sized
CRediT authorship contribution statement
Wojciech Aniszewski: Conceptualization, Methodology, Software, Formal analysis, Investigation, Data curation, Writing - original draft, Writing - review & editing, Visualization. Youssef Saade: Formal analysis, Investigation, Writing - review & editing, Software. Stéphane Zaleski: Supervision, Project administration, Funding acquisition, Conceptualization, Investigation, Methodology. Stéphane Popinet: Software, Supervision, Conceptualization, Methodology.
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowlgedgments
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 View8, 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
References (41)
- et al.
Transition in a numerical model of contact line dynamics and forced dewetting
Journal of Computational Physics
(2018) Improvements, testing and development of the ADM-τ sub-grid surface tension model for two-phase LES
J. Comput. Phys.
(2016)- et al.
Volume of fluid (VOF) type advection methods in two-phase flow: a comparative study
Computers & Fluids
(2014) Chapter 7 - coating
- et al.
A second-order projection method for the incompressible navier-stokes equations
J. Comput. Phys.
(1989) - et al.
Modeling of a falling-film photocatalytic reactor: Fluid dynamics for turbulent regime
Applied Mathematical Modelling
(2016) - et al.
An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension
Journal of Computational Physics
(2018) - et al.
Efficient symmetric positive definite second-order accurate monolithic solver for fluid/solid interactions
J. Comput. Phys.
(2012) Laminar withdrawal with appreciable inertial forces
Chem. Eng. Sci.
(1970)LEIS for the prediction of turbulent multifluid flows applied to thermal-hydraulics applications
Nuclear Engineering and Design
(2010)
Impingement of an axisymmetric jet on unheated and heated flat plates
J. Fluid Mech.
An accurate adaptive solver for surface-tension driven interfacial flows
J. Comput. Phys.
A quadtree-adaptive multigrid solver for the Serre–Green–Naghdi equations
Journal of Computational Physics
A 3D adaptive mesh refinement algorithm for multimaterial gas dynamics
Physica D
Height functions for applying contact angles to 3D VOF simulations
Int. J. Numer. Methods Fluids
Multigrid techniques: 1984 guide with applications to fluid dynamics
Wave evolution on a falling film
Annu. Rev. Fluid Mech.
Jump conditions and etropy sources in two-phase systems. local instant formulation.
Int. J. Multiphase Flow
Cited by (11)
Numerical simulation of the accumulation process in modulated liquid jets of defined volume using the volume-of-fluid methodology
2023, International Journal of Multiphase FlowA fourth-order accurate adaptive solver for incompressible flow problems
2022, Journal of Computational PhysicsEvolution of waves in liquid films on moving substrates
2023, Physics of Fluids
- ☆
Wojciech Aniszewski would like to dedicate all his work on this paper to the dear memory of Antonina Gilska (nee Hessler, 1936–2019).