EULER-LAGRANGIAN MODEL OF PARTICLE MOTION AND DEPOSITION EFFECTS IN ELECTROSTATIC FIELDS BASED ON OPENFOAM

In order to study the powder coating process of metal substrates, a comprehensive, numerical 3D Eulerian-LaGrangian model, featuring two particle sub-models, has been developed. The model considers the effects of electro-static, fluid-dynamic and gravity forces. The code has been implemented in C++ within the open source CFD platform OpenFoam®, is transient in nature with respect to the applied LaGrangian particle implementation and the electro-static field calculation and is stationary regarding fluiddynamic phenomena. Qualitative validation of the developed solver has already been achieved by comparison to simple coating experiments and will hereby be presented alongside a thorough description of the model itself. Upon combining knowledge of the relevant dimensionless groups and the numerical model, a dimensionless chart, representing all possible states of coating [1], was populated with comprehensive, exemplary cases, which are shown here as well.


INTRODUCTION
The powder coating process is widely spread and commonly used in industry.However, knowledge about detailed phenomena and concerning parameter-effect relations remains of predominantly empirical nature.Numerical modelling and recently introduced characterization methods [1] can help to deepen the understanding of all involved material-and process parameters and can thus lay the foundation for future, knowledge based improvement of the technology.
In addition to previous modelling efforts (see e.g.: [2] to [5]), this work presents a new Eulerian-LaGrangian 3D model, which has the ability to grasp detailed particle deposition effects on the substrate surface, as well as phenomena regarding large numbers of moving particles, within the coating chamber.

Over View
While chapter 1.2 of this paper provides a very brief introduction into the powder coating technology by discussing an exemplary coating set-up, chapter 2 presents the OpenFoam® [9] based, 3D modelling approach of the electro-static (see chapter 2.1), fluid-(see chapter 2.2) and particle-dynamic (see chapter 2.3) conditions in the vicinity of a coating pistol and a coated substrate.In order to validate the numerical model, coating experiments have been conducted under defined conditions on a simple metal plate.A comparison between experimentally-and numerically obtained results is presented in chapter 3. To draw a bottom line, chapter 4 shows the application of the validated simulation tool in a qualitative parameter study.A dimensionless charting method, lain out in [1], is thus combined with more concrete, comprehensive visualizations of the distinct coating states.

Exemplary Powder Coating Set-Up and Basics of Field Formation
A sketch of one possible, exemplary powder coating set-up is depicted in Figure 1.Thereby an airflow laden with e.g.thermoplastic or thermoset polymer coating particles with diameters D p between 5µm and 300µm, is led past a high voltage electrode and a baffle nozzle [6].In combination with a grounded, metallic substrate of zero potential, the negatively charged electrode causes the formation of an electric potential field Ψ(x,y,z) (V).The potential field in turn is the cause of an electric force field (x,y,z) (N/C) whose field lines originate from the metallic substrate and enter the electrode.
Figure 1: Sketch of exemplary powder coating set-up with coating pistol (left), grounded substrate (right) and a sketch of acting forces on the negatively charged coating particles, [7].
Due to the presence of negatively charged oxygen ions O - 2 [1], the coating particles are charged up to a maximum negative surface specific charge q p (C/m 2 ), as they are carried past the electrode.After charging, the particles are impacted by the spatial electric field.Upon entering the space between the charging zone and the substrate, which is referred to as the coating chamber, three main sorts of acting forces determine particle motion and particle-substrate deposition behaviour: gravity forces  !, electric forces  !" and fluid-drag forces  !(see Figure 1, right).

THE 3D OPENFOAM MODEL OF COATING PARTICLE MOTION AND DEPOSITION
A numerical 3D model of the coating process within a coating chamber has been developed.It considers LaGrangian particle motion, the effects of electro-static, fluiddynamic and gravity forces.Two separate particle sub models have been created: A detailed large LaGrangian particle model and a simpler, but more efficient small LaGrangian particle model.The code is implemented in C++ within the open source CFD platform OpenFoam® [9], is transient in nature with respect to the applied LaGrangian particle implementation and stationary regarding fluid-dynamic phenomena.Within the small LaGrangian solver, particle charges do not impact the electro-static field, which can then considered to be stationary.Within the large LaGrangian solver however, the electric field changes transiently with respect to the particle borne charge distribution.

Geometry, Mesh and Flow Simulation
In a slightly simplified approach, the stationary, turbulent flow analysis within a coating chamber can be decoupled from the electro-static field-and particle motion calculation.Thus the OpenFoam® standard solver pimpleFoam [11], which neither considers particles, nor electric field effects, was applied to conduct the CFD simulation.A SST-k-ω turbulence model [10] was chosen.
A CAD model of an accompanying experimental set-up (chapter 4), involving the coating of a small metal plate within a testing chamber, was created (Figure 2, left) and meshed, using the OpenFoam® based meshing utility SnappyHexMesh [9].The resulting, structured hexahedral mesh features 5e5 cells and is shown in Figure 2, right.
Figure 2: CAD model of the coating chamber (Figure 8) including the substrate (left) and structured, hexahedral mesh (right).
While Figure 2, left includes definitions regarding the boundary patches, Table 1 sums up the chosen boundary conditions in terms of velocity vector field U and scalar pressure field p.As seen in Table 1, simulations have been conducted for various pistol inlet velocities between 1.00m/s and 20.0m/s.Figure 3 shows an exemplary snap shot of a vector field plot of the flow field between pistol inlet (left) and metallic substrate (right).Note the slight downdraft at the substrate.It is caused by a ventilation system beyond the outlet (Figure 2) of the model, which is represented by a fixed outlet velocity of 0.88m/s.

Electro-Static Field Model
An electric field is formed between the negatively charged electrode and the grounded substrate.Mobile, spatial charges, expressed as charge densities ρ el,tot (C/m 3 ), carried by either coating particles or ions, might have a certain impact on the electric field formation (Eqn. 2 and Eqn. 4).However, the governing electric effects remain to be mainly static in nature such that electro-dynamics do not play a role in this context.Thus the static version of Faraday's Law of Induction (Eqn. 1) and Gauss' Law concerning Electric Fields (Eqn.2), lead to the insignificance of effects related to the magnetic field  (Ns/Cm), the introduction of an electric potential field Ψ (Nm/C), the electric force field  (N/C) (Eqn.3) and Poisson's Equation (Eqn.4).

∇ ×
In Eqn. 2 and Eqn. 4, ε 0 (C 2 /Nm 2 ) is the electric field constant.The field-impact of charges, carried by coating particles, can be established by combining calculated particle charge densities ρ el,p with the ionic space charge density field ρ el,I .Thus the total space charge density field is the local superposition of the two fields (Eqn.5).An integration of charge impact, based on ion-and/or electron-distribution, requires the introduction of a full conservation equation for ρ el,I (Eqn.8).The latter can be derived from continuity, charge conservation (Eqn.6), Gauss' Law (Eqn.2) and a formulation for charge flux density  (C/sm 2 ) (Eqn.7), [12].
Eqn. 5 In Eqn. 7 and Eqn. 8, b (Cs/kg) is an expression for charge mobility under the influence of the electric field, D el (m 2 /s) is a charge diffusion coefficient, which includes the effect of repulsion of equally signed charges and  a is the airflow velocity Effects related to mobile spatial charges, described by Eqn. 6 to Eqn. 8, are mainly relevant in regions of very high spatial charge densities close to the electrode.Since the focus of our considerations lies mainly on the coating chamber, the inclusion of these effects has as of now been omitted within the model.
Eqn. 1 to Eqn. 4 have been implemented within the framework of the pimpleFoam solver [11] and lead to a simplified calculation of the Ψ-and consequentially the  field.
Figure 4 shows a snapshot of calculated electric field lines between the electrode and a metallic plate.Once the electric field is calculated and the position  !,! , of a particle i with surface area A p,i (m 2 ) and surface specific charge q p,i (C/m 2 ) , is known, the electric force effect on that particle can be obtained by Eqn. 9.

Particle Dynamic Models
Two LaGrangian particle dynamic models have been developed and applied within the given task of modelling electro-static powder coating processes: A large and a small LaGrangian model.In this context the terms large and small refer to particle sizes in relation to the applied finite-volume cell sizes.

Large vs. Small Lagrangian Particle Solver
The small particle solver has been created to deal with larger amounts of small particles, meaning multiple particles per fluid mesh cell.Its main area of application is the open space of the coating chamber, where the coating particle cloud is to be studied.The large particle solver, on the other hand, has been tailored to model smaller amounts of highly resolved, large particles, meaning individual particles spanning across multiple fluid mesh cells.Its main application is the modelling of detailed deposition-and interaction phenomena, close to the substrate.The basics of the large particle solver have been published in [13] to [15].It has the capability to handle nonspherical shape effects, as well as full bi-directional particle-flow and particle-particle interactions.Within this work, the large particle solver has been extended to be able to handle particle borne charges and their interactions with the electric field.The small particle solver, on the other hand, neglects particle-particle interactions as well as particle effects on the flow and on the electric field.Figure 5 presents a head-on comparison of the two solvers.
Figure 5: Head-on comparison of the large particle solver (see also [13] to [15]) and the small particle solver.Note that, in addition to the PISO [9] algorithm, the PIMPLE [11] algorithm has been used for flow solution as well.

Implementation of LaGrangian Particle Motion
While the large particle solver is based on the principles of the icoLagrangianFoam solver [16], the small particle solver is a combination of OpenFoam's solidParticle class with the pimpleFoam solver [11].In both cases, the original ideas behind drag forceand the numeric implementation of LaGrangian particle motion are similar [16].

Implementation of Fluid Drag
Three modes of particle drag force calculation, accounting for varying states of turbulence, are implemented.In any case, Eqn. 10 is used to calculate the fluid drag force F f at time step t i .However, the calculation of particle relaxation time τ p , depending on the drag coefficient c d , varies from case to case.Table 2 sums up the underlying equations, where µ a (Pas), ρ a (kg/m 3 ) and Re p (-) are the dynamic viscosity of air, the density of air and the particle Reynolds number respectively.

Numeric Implementation of LaGrangian Particle Motion
The numeric implementation of LaGrangian particle motion is simply about resolving the particle momentum equation (Eqn.11), based on Newton's second law, by discretizing the particle acceleration term.In Eqn.11 m p (kg) is the particle mass.
The procedure to move the numerical particle forward is this: First, in order to obtain a more stable implementation, the new particle velocity at time t i+1 ,  p (t i+1 ) instead of  p (t i ) is used to discretize the acting fluid drag force, according Eqn.12.Then the new particle velocity at time t i+1 can be expressed and calculated according Eqn. 13.Finally the new particle position at time t i+1 ,  p (t i+1 ) is computed according Eqn.14.Note that Δt p , used in Eqn. 13 and Eqn.14, is a particle sub time step, which can be considerably smaller than the prevailing fluid time step Δt f , [8].

Qualitative Model Results
In order to qualitatively validate the solver's functionality, simplified test geometries have been chosen.An exemplary snapshot of a case, modelled by the large particle solver is seen in Figure 6.Thereby several dozen, highly resolved, positively charged, large particles are dragged by a fluid flow from the left towards the right, while being attracted by cylindrical, negatively charged fibres.Since the coupling between particle borne, spatial charges ρ el,P and the electric field (Eqn. 2 and Eqn.4) is activated, the impact of particles on the electric potential field can be visualized in terms of isopotential planes of the Ψ-field.In effect the simulation shows that loaded particles tend to be deposited on free spaces along the fibres, thus form a repulsive electric shield around the oppositely charged fibre and consequentially repel oncoming particles, which then penetrate deeper into the fibre structure towards the right.Figure 7 presents a snapshot of a small particle solver demonstration case.Hereby the real-life set-up of our testing facility (Figure 8 and for mesh see Figure 2) has been re-created.Coating particles, including particle velocity vectors are visualized.They are moving from the electrode towards the metallic plate and deposit there.Electric field lines are visualized as well.It can be seen that particle-substrate deposition occurs on both sides of the plate.

EXPERIMENTS AND VALIDATION
Efforts to validate the small particle solver (chapter 2), have been under taken and are discussed in the following.

Experimental Set-Up
In order to achieve some degree of validation of the small LaGrangian particle model, we began by choosing the simplest possible substrate geometry: a quadratic metal plate.As shown in Figure 8, the plate was placed in the centre of a small, experimental coating chamber, air draft was initiated by a ventilation system beyond the outlet filter and a cloud of coating particles was injected by applying a flat-jet-nozzle coating pistol.

Evaluating the Experiments
After concluding the coating procedure, the coated substrate was placed in front of a Coat-Master measurement device [17] (Figure 9, left).Then a measurement pattern, yielding 23 measurement locations on the plate's front-and 23 locations on the plate's back-side, was defined (Figure 9, right).By thermally scanning these selected locations, the measurement device could retrieve a profile of relative coating thicknesses on the front-side (Figure 12, left) and on the back-side (Figure 10 and Figure 11) of the plate.

Solver Validation
Alongside the experiments, coating solver simulations have been carried out within a re-created experimental set-up (Figure 2).The small LaGrangian particle solver cannot match real-life particle numbers, but is capable of grasping several thousands of individual, representative particles.The numerically obtained deposition behaviour of these particles was monitored, evaluated and compared to measured data from the Coat-Master measurement.Figure 10 shows a direct comparison of measured, relative coating thickness values within the 23 measurement locations and a qualitative 1:1 comparison to the computed, visualized deposition behaviour on the back-side of the plate.Both results show the same, characteristic window frame effect.This known phenomenon in the coating industry is about the occurrence of noticeably thicker coatings towards the edges of a substrate, [6].Figures 11 and Figure 12 present a full quantitative comparison of experimentally obtained, relative thicknesses (left) and simulated relative thicknesses, derived from representative particle numbers (right).Figure 11 shows results from the back-side of the plate, while Figure 12 depicts conditions at the front-side.Note that in Figure 11 and Figure 12 the respective measurement locations A3 + A4 and A2 + A3 are left blank.Here the holding device was attached during the coating and no particles were deposited.Comparisons of experimental and simulated coating behaviour, shown in Figure 10 to Figure 12, yield good qualitative correspondence.Even though quantitative matches cannot be achieved, the simulation grasps the main tendencies and outstanding phenomena, such as: the increasingly strong window frame towards the bottom of the back-side (Figure 11), or the locations of maximum, relative coating thickness on the front-side (Figure 12).Considering the variety of variable parameters involved in the coating process, as well as the measurement procedure, the over-all degree of similarity is quite promising.

RESULTS
So far the most valuable results have been achieved by combining numerical output achieved by the OpenFoam® solver (chapter 2), with new methods to characterize the given coating state.Thus we were capable of visualizing average force effects on different particle size classes along their respective flow paths (chapter 4.1) in triangle notation [1].In addition and most importantly the new solver has been used to interpret and understand the dimensionless chart of coating states [1] in terms of its meaning for particle cloud shapes and actual coating results (chapter 4.2).

Particle Flow Paths and Acting Forces in Triangle Chart Notation
The coating solver was used to study the relation of acting forces on particles of different size classes on their way from the electrode to the metallic plate.For each investigated size class, five particles were injected at selected locations, close to the electrode.These locations were chosen such that they represent extreme cases of possible flow paths: at the centre, the top, the bottom, to the very left and to the very right of the injection cone.Results of this numerical investigation were averaged for all five injection-positions per size class, plotted in triangle chart notation and are shown in Figure 13.The triangle chart notation has been introduced in [1].Note that the chart dimensions π el , π F and π g are dimensionless groups representing the impact of electric-, fluid drag-and gravity forces respectively.As seen in Figure 13, the force relations at injection positions are always located rather towards the corner of fluid drag dominance (blue).As particles penetrate further into the coating chamber, they shift towards gravitational-(yellow) and from there towards electro-static dominance (red).Larger particles rather tend towards gravitational dominance, than smaller ones.Close to the substrate, being the final particle destination, electro-static forces tend to dominate.

Dimensionless Chart of Coating States, Populated by Model Results
The chart of coating states, published in [1], is a method to characterize all possible parameter states of the powder coating procedure for a given geometry.The newly developed coating solver has been used to re-create exemplary situations of various process parameter combinations for each zone within the chart of coating states.Upon qualitative inspection of the results, it became clear that the shape of the coating particle cloud is an important indicator for the over-all state of the coating process.The results of these inspections are compiled within in Figure 14.While states of electrostatic dominance (red) exhibit spherical coating cloud shapes, gravitationally dominated states (yellow, green) show more concentrated particle jets, which do not even hit the target but merely fall towards the bottom.Flow dominated states (blue, orange) on the other hand, cause clouds that are oriented along the streamlines of the air velocity field.The tendency of particles hitting the bottom edge of the plate within those flowdominated zones, can easily be explained by the state of the calculated flow field, which exhibits a downdraft towards the outlet filter (Figure 3).The combination of the dimensionless charting method and the numerical solver shows very plausible results throughout the entire spectrum of governing effects.In addition to merely evaluating particle cloud shapes, the solver can also produce predictions for representative particle distributions on the front-side and back-side of the exemplary metallic plate.The respective results are shown in Figures 15 and Figures    The numerically obtained coating results help to understand the nature of certain phenomena, experienced in the actual application.Thus it becomes obvious, that higher electrode potentials might lead to stronger window-frames (Figure 16) on the back-side of the substrate, while a more evenly distributed coating thickness is to be expected on the front-side (Figure 15).In general the results shown in Figure 15 and Figure 16, point out that the only chart zone, providing adequate coating results at all, is the region of electro-static dominance (red).In real life and for further investigations this zone is certainly the focal area, which needs to be resolved in more detail.

CONCLUSION
It has been shown that the combination of numerical modelling and new characterization methods constitutes a powerful tool to thoroughly understand the coating process.Where the chart of coating states (published in [1]) focuses more on the big picture, the developed OpenFoam® solver is capable of re-creating detailed, individual scenarios.The solver's main aspects, including the two LaGrangian particle sub models have been laid out in this paper.The small LaGrangian solver has been experimentally validated and can now help to get a plausible idea of what varying coating parameter states actually mean in terms of real life application.Qualitative investigations of coating particle cloud shapes and coating particle deposition patterns have served as examples.
The presented results can now form the starting point for conducting knowledge based parameter improvements for any conceivable powder coating set-up.

Figure 3 :
Figure 3: Vector field plot of the velocity field for a set pistol inlet velocity of 1.0m/s.Coloration corresponds to local velocity vector magnitude.

Figure 4 :
Figure 4: Snapshot of visualized electric field lines sourcing from metallic substrate and sinking at the negatively charged electrode.

Figure 6 :
Figure 6: Snapshot of qualitative demonstration of large particle solver functionality.Large, positively charged particles coming in from the left, dragged by fluid flow and attracted to negatively charged fibres.Particle borne charges are coupled to the spatial potential field.Visualization shows iso-potential planes within the Ψ-field.

Figure 7 :
Figure 7: Snapshot of re-created experimental set-up case.Small particles with velocity vectors are moving from electrode towards metallic plate.Situation is depicted from the side (left) and from behind the substrate (right).

Figure 8 :
Figure 8: Flat-jet-nozzle coating pistol (top, left), experimental coating chamber with outlet filter (bottom, left) and snapshot of on-going coating procedure of a metal plate (right).

Figure 10 :
Figure 10: Qualitative comparison of measured, relative thicknesses (left) and computed visualized particle distributions (right) on the back-side of the metallic plate.Measurements and simulations show the same, characteristic window frame effect.

Figure 11 :
Figure 11: Direct comparison of 23 relative thickness measurements (left) and corresponding simulations (right) at the back-side of the plate.

Figure 12 :
Figure 12: Direct comparison of 23 relative thickness measurements (left) and corresponding simulations (right) at the front-side of the plate.

Figure 13 :
Figure 13: Triangle chart notation of calculated averaged forces acting on moving particles of size classes: 20µm (red), 40µm (blue), 60µm (green) and 100µm (purple).Each dot stands for another particle location along the flow path.

Figure 14 :
Figure 14: Dimensionless chart of coating states according to [1], with exemplary numerical model results, showing the basic shape of coating particle clouds within the various state-zones. 16.

Figure
Figure 15: Dimensionless chart of coating states according to [1] with exemplary numerical model results, showing coating particle distributions on the front-side of a metallic plate.

15 :
Figure 15: Dimensionless chart of coating states according to [1] with exemplary numerical model results, showing coating particle distributions on the front-side of a metallic plate.

Figure 16 :
Figure 16: Dimensionless chart of coating states according to [1] with exemplary numerical model results, showing coating particle distributions on the back-side of the plate.

Table 1 :
Boundary patches with boundary conditions for velocity-and pressure fields.