Multiphysics Simulation of Particle-Surface Interaction and its Effect on Powder Patterns and Process Optimization

In industrial processes, powder coating is widely utilized to attain functional or aesthetic surface properties on manufactured parts. A Eulerian-Lagrangian Multiphysics solver has been developed within the OpenFOAM framework in order to simulate and optimize such processes with regards to coating efficiency and homogeneity. In the scope of this study, the powder particle-substrate and particle-particle interactions that occur on the surface of a substrate during the coating process are investigated. This is instigated by the observation that some particles glide over the substrate, rather than sticking to the substrate upon first contact. The phenomenon is governed by the balance of pressure, fluid shear stress traction, electrostatic particle-particle repulsion and gravity forces on the substrate. On the basis of experimental data previously gathered, it is demonstrated that the surface interactions are essential to predicting the coating outcome accurately enough, such as to serve as basis for later process-optimization steps. Furthermore, a dimensional analysis illustrates the weight of the individual force contributions on the overall force balance.


INTRODUCTION
Powder coating is a widely used industrial process to paint manufactured goods, falling under the general category of surface coating methods to attain functional properties (e.g., corrosion resistance) or aesthetics qualities.The process involves micron range powder particles being electrically charged in a coating pistol, where they are transported along with a carrier airflow.After emerging from the pistol, the particles are driven by electrostatic forces as well as fluid drag across the coating chamber, towards the electrically grounded substrate, where they either get blown off or stick to the surface to be cured, consequently coating it [12].The process is governed by the Multiphysics interactions of gravity, electrostatics and air-flow forces with the individual powder particles.Due to the complexity of the interactions, often ad-hoc heuristic, trial-and-error testing methodology is utilized in the industry in order to determine acceptable process parameters.Particularly during process design-and (re-) adjustment phases, this leads to the intense consumption of resources in terms of energy, labour and powder material.More significantly, this approach does not necessarily lead to an optimized process, translating to loss of efficiency throughout the whole lifecycle of the process.
Hence, there is a great improvement potential through the use of a predictive simulation tool to tackle the issues outlined.The most essential requirement from such a predictive tool would be accuracy, to be validated on the basis of experimental data generated by actual coating processes.Multiphysics simulation methods can be employed to develop such predictive tools.One particular class of such tools focuses on Eulerian-Lagrangian simulation approaches.In such efforts the Eulerian carrier airflow field within the coating chamber as well as the Eulerian electric field emanating from the coating pistol electrode are numerically calculated.Then the resultant forces within the field are imposed onto individual Lagrangian powder particles to predict their trajectory and subsequently the coating pattern on the substrate.In previous contributions the authors introduced, a Eulerian-Lagrangian Multiphysics solver, developed within the OpenFOAM framework specifically for the simulation of surface coating processes [1]- [5] adapted from one developed for the filtration industry [6]- [8].Similar approaches have been described in the literature [13]- [16].In the work of Shah et.al. [13], a constant surface-charge-to-mass-ratio is imposed based on the particle diameter and later on expanded into a parameter study by Li et.al. [14].In contrast Böttner et.al. [15] and Ye et.al. [16] have relied on experimental data for the surface charging.The aforementioned works have not given particular attention to interactions that occur on the substrate surface once the particles have reached it; hence the possibility of particles gliding alongside its surface is not considered.This simplification conflicts with observations in real coating applications and renders the prediction of the final coating pattern impossible.Considering that the assessment of the homogeneity of the coating is crucial in industrial applications, this issue needs to be addressed.In the current study, the interactions on the substrate surface are taken into account through a force balance depending on pressure, fluid shear stress, the electric field state on the substrate along with gravity and electric particleparticle repulsion force resulting from the agglomeration of particles.This allows the particles to glide along the surface allowing their redistribution to predict the final coating pattern with increased accuracy.

METHODS
This study is based upon the Multiphysics, OpenFOAM based Eulerian-Lagrangian powder coating solver, described within [1]- [5].The solver is capable of considering two-way-(fieldparticle; particle-field) and three-way coupling effects (field-particle; particle-particle; particle-field) between particles and their surroundings.However, this study focuses on a oneway coupled scenario between fluid and particles as well as separately treated particle-particle interaction effects.Thus, in order to simulate the interaction of prevailing phenomena, the solver computes flow-and electrostatic field separately in a Eulerian frame of reference.The gravity field is then superimposed.Only in the next step are the trajectories of the injected particles calculated in a Lagrangian frame of reference and based on their location and the force balance prevailing thereupon.

Governing Equations
The physics simulated by the solver is governed by: i) the fundamental momentum balance equations for fluid flow, ii) a steady state version of Maxwell equations and iii) a force balance for each particle.Given the typical application range, the carrier air-flow can be safely assumed to be incompressible for which the continuity equation and the 3D-momentum balance can be represented as in Equations ( 1) and ( 2) respectively [17].

𝛻𝛻 ⋅ 𝑈𝑈 � �⃗ = 0
(1) where the variable  � �⃗ represents the velocity field of the carrier air, the   its effective kinematic viscosity and   its density.
For the electrostatic field, the governing equations are the static Maxwell equations without magnetic effects as represented by the Poisson-type Equation (3) [18].
In Equation (3),  represents the electric potential which is related to the electric field  �⃗ through Equation ( 4),   the volumetric charge density and  0 the permittivity of free space.The corona formed by the electrode within the coating pistol features potentials in the range of tens of kilovolts and is the dominant factor in defining the properties of the electric field.As a consequence, the impact of space charges can be neglected and Equation (3), simplifies into a Laplace-type equation as represented in Equation (5).
The particle motion can be calculated based on fluid drag, the electrostatic forces, gravity and particle inertia.The fluid-drag force is estimated from the empirical relations for a spherical particle in isolation [24].
Re  �1 + 0.15Re  0.687 � (7) In the case of the isolated cylinder, the drag force depends on the particle Reynolds number Re  given in Equation (6), where   ����⃗ represents relative fluid-particle velocity,   its diameter and   the kinematic viscosity of air.In Equation (7), the empirical drag force coefficient   is represented based on Re  , which in turn allows the calculation of the magnitude of the drag force   ����⃗ in accordance with Equation (8), where   is the projected area of the spherical particle.The formulation of the fluid drag force along with the electrostatic-and gravitational force, give the final expression for computing the acceleration of individual particles   ����⃗.
In Equation ( 9), the electric field  �⃗ is multiplied by the surface electric charge density of the particle ′′ to attain the electrostatic force, and the final right-hand term represents the effect of gravity  ⃗, where   and   are fluid-and particle density respectively.

Substrate Surface Interactions
Equations ( 5)-( 9) describe the motion of the particles within the coating chamber as they emerge from the coating pistol and make their way to the substrate.Once they impinge on the substrate, it can be observed experimentally that they either stick immediately or glide along the surface until they eventually stick or get blown off.This phenomenon can be attributed to the ratio of the sum of forces on the particle acting normal to the substrate  ⃗  to that acting tangential to it  ⃗  as depicted in Figure 1.If this ratio is above a critical threshold, the particle would stick.If the ratio is below that threshold the particle would either glide along the surface or get blown off.
Figure 1: A particle on the substrate surface is exposed to a total force which can be decomposed into its normal-and its tangential components.
In order to quantify the particle-substrate force components ratio, the following effects are taken into account: i) pressure force; ii) fluid traction; iii) electrostatic force; iv) gravity and v) electrostatic particle-particle repulsion due to particle accumulation on the substrate.The pressure force  ⃗  is always normal to the substrate surface area vector   ���⃗ as expressed in Equation ( 10) and is obtained from the flow field.In Equation (10),  is the value of the fluid pressure field at the cell-centre of the finite volume cell engulfing the particle at its current location on the substrate surface.
The fluid traction force,  ⃗  , on the substrate surface can also be attained from the flow field as formulated in Equation (11).Here  is the deviatoric part of the shear rate tensor.
In the OpenFOAM framework  is calculated in accordance with Equation ( 12) [21].
It should be noted here that the effective viscosity allows the turbulent flow effects to be factored into the shear stress rate calculation.In addition, the last term on the right-hand side of Equation ( 12) would vanish for an incompressible fluid as the trace of the velocity gradient would amount to zero in accordance with Equation (1).The fluid traction force represents the force acting on the particle due to the velocity gradients near the substrate surface and is predominantly a force in the tangential direction.The computation of the electrostatic force acting on any particle at the substrate surface  ⃗  first requires a result for the local electrostatic field  �⃗ .This can be derived from the solution of the electric potential Equation ( 5) and subsequently its relation to the electric field according to Equation (4).The particle force itself is then related to the surface area of the particle as well as to its surface charge density ′′ through Equation (13).
The electrostatic force contribution also acts in the normal direction to the surface.The electrostatic particle-particle repulsion force,  ⃗  , resulting from the dense accumulation of particles on the surface of the substrate can be modelled in a similar fashion.Hence it can be expressed analogously as in Equation (14).
The difference to Equation (13) lies in the calculation of the electric repulsion field,  �⃗  , which depends on the amount of charge accumulated at a certain location on the substrate.To calculate the amount of charge accumulated on the substrate, its geometric discretization is utilized.As the simulation progresses, the number of particles in each discretized surface element is monitored, and local particle charges are summed up according to Equation (15).
where,   is a factor indicating whether the  ℎ particle is within the bounds of the  ℎ surface element.If the  ℎ particle is in fact located within the  ℎ surface element, its value is counted as one, and its charge is added to the total charge of the  ℎ surface element, and otherwise its value is zero.This operation can be interpreted as a transfer of the particle charges from a Lagrangian field to the substrate surface as a Eulerian field as illustrated in Figure 2. Once the total surface charges on the substrate surface elements are attained, they can be interpolated onto the vertices of each element as illustrated in Figure 3.This interpolation is based on the inverse distance of the surrounding elements of the vertex, as formulated in Equation (16).
where the charge at the vertex ,   , depicted in Figure 3 is obtained from the summation of all the neighbouring charge values,   , multiplied by a scaling factor  , .This scaling factor is based on the inverse of the distance of the neighbour cell centre at location   ����⃗, to the coordinates of the vertex   ����⃗ as expressed in Equation (17).
After the charges are attained at the surface element vertices, they are treated as point charges and utilized to calculate the repulsion electric field,  �⃗  , on the surface elements in accordance with Equation (18).
In summary, (i) the charges on the discretized substrate surface elements are summed up, (ii) they are interpolated onto the surface element vertices, which are (iii) treated as point charges to calculate the repulsion electric field back at the surface elements.Once a particle impinges on a specific surface element, the repulsion force can be calculated as formulated in Equation ( 14).This method has been devised to utilize the built-in interpolation functions of OpenFOAM so that the calculation of the repulsion, to be updated at each iteration, comes at minimal computational performance cost.Finally, the total force acting on a particle can be calculated by summing up all mentioned force contributions, with the addition of the gravitational force as in Equation (19).
Figure 3: The interpolation from the surface elements to the vertices.A zoom into the geometrically discretized substrate with the calculated charge distribution demonstrates the interpolated vertex charges based on the neighbouring surface element values.

Numerical Methods
Equations ( 1), ( 2) and ( 5) form the basis of the numerical simulations to resolve the electric field and the flow field within the coating chamber and subsequently lead to the forces an individual particle will experience.Since the focus of this study is on the behaviour of the particles on the substrate, a static configuration where the coating pistol is at a fixed distance of 15 cm with respect to a 10 x 10 cm aluminium plate was used.The dimensions of the coating chamber utilized were a 0. The flow field was simulated using the Finite Volume Method (FVM) [19] within the OpenFOAM framework [20]- [21].The turbulence of the carrier airflow is modelled with the RAS k-ω shear stress transport model [22] using the SIMPLE algorithm [23].The underlying solver also incorporates electrodynamic effects and forms the basis of the software coatSim [9]- [10], which has the capability to simulate pistols performing periodic translational and rotational motions directly on the cloud with the Massive Simultaneous Cloud Computing (MSCC) technology [11].In the scope of this study the software has been utilized to set up the simulation case without any motion.
Applying said process parameters, the simulations were run with either 25'000 or 250'000 particles with the sticking force ratio threshold (ratio of normal forces to tangential ones) set to either: off (stick upon impact), 5 or 10.
Figure 4: The coating configuration for the simulations.The coating pistol was positioned 15 cm from the centre of the plate within the coating chamber with the dimensions 0.4m x 0.3m x 0.7m.The coating chamber was vented on the backside, while besides the pistol flow, there was flow inlet on the front surface.

Experimental Data Collection
An experimental campaign was previously carried out to validate simulation efforts [2].In the scope of that campaign, 216 separate measurements with the Coatmaster 3D [25] measurement device were performed on 69 different deposition experiments.The Coatmaster 3D technology relies on thermal optics, where the surface to be measured is heated in a pulsed manner to assess the coating thickness.The key advantage of the method is that it is contactless and provides the detailed thickness distribution on the substrate without requiring it to be cured in the oven.Hence it is particularly suited for collecting large amounts of experimental data.
The experimental data was assessed based on the methods described in [3].In the scope of this study, data collected at 35 kV with the substrate coating pistol distance of 15 cm on the 10 x 10 cm aluminium plate was used for comparisons.

RESULTS
The flow and the electrostatic fields obtained from the simulation along with the resultant pressure, fluid shear traction and electrostatic fields on the substrate surface are illustrated in Figure 5.It can be observed that the pressure field is more pronounced at the centre of the plate, while the fluid shear traction shows the opposite trend.The electric lines run perpendicular to the plate, and they are especially pronounced at the edges of the plate indicative of the window frame effect.The window frame effect describes regions of thicker coating prevalent at the edges of the plate, usually observed on the back surface (opposite to pistol facing surface).The pressure force and the electrostatic forces are the major normal components of the total particle force on the plate surface, while the fluid shear traction is the main tangential contribution.Particle-particle repulsion forces add additional tangential components, while gravity contributes both normal as well as tangential components in relation to substrate surface orientation.In the case at hand, the plate-like substrate was positioned upright, such that gravity acted tangentially to its surface.Figure 6 illustrates the simulation results regarding the total force distribution for a 30 µm particle on the surface of the plate-shaped substrate, for the cases of 25'000 and 250'000 computed particles.
The normal force patterns look quite similar for the cases of 25'000 and 250'000 computed particles, with the differences appearing on the tangential force patterns due to the increased effect of the repulsion force for the case of 250'000 computed particles.The resultant coating thickness field obtained on the pistol-facing surface (front) of the plate substrate can be illustrated as in Figure 7.It can be observed that the force ratio threshold starts to have a significant effect only at around the value of 10.The coating patterns for the case of sticking upon impact and sticking above the threshold of 5 are quite similar with largest discrepancies occurring around the edges of the plate.In order to relate the final coating pattern to the surface forces acting on the coating particles, it proved beneficial to non-dimensionalise the latter.Thus, each force-component was normalized by the force required to accelerate a particle of 30 µm by 1 m/s 2 .
Figure 7: The coating thickness profiles obtained on the plate surface; the top row with 25'000, the bottom row with 250'000 computed particles, the leftmost column for particles sticking upon impact, the middle column for particles sticking at force ratio above 5 and the rightmost above 10.
The non-dimensionalised force distribution on the substrate plate is illustrated in Figure 8.It can be observed that the magnitude of the normal component of the total particle force is larger than that of the tangential component on the whole front surface of the plate, where the ratio gets its lowest values at the edges of the plate (see Figure 8 a).Of the normal forces, the pressure force is mostly effective at the centre of the plate.Due to the horizontal slit configuration of the nozzle the flow splits into two jets impinging on the plate producing two distinct stagnation points with high pressure (see Figure 5 a), dropping to low values towards top-and bottom edges (see Figure 8 c).The second significant normal-force contribution besides pressure comes from electrostatic forces.Electrostatic force contributions dominate especially at the edges of the plate but are weakened in magnitude towards the centre (see Figure 8 e).In this case, gravity acts purely tangentially to the substrate surface and of course is perfectly homogeneous (see Figure 8 b).The shear traction force is lowest in magnitude at the centre of the plate and increases in magnitude towards the edges (see Figure 8 d) The repulsion force is largest at the sites where the particle concentration is the largest (see Figure 8 f and Figure 9) Finally, the coating pattern obtained experimentally can be compared with that of the simulations as illustrated in Figure 9.It can be observed that an ellipsoid region highlighted in purple in Figure 9, can be attained only with the consideration of a high enough particlesticking-threshold criterion regarding the ratio of normal-to tangential forces acting on the surface particles.In other words: the simulation yields higher correspondence to experiments, if the particles are given the freedom to move across the substrate surface.However, it can also be observed that at a higher threshold, two strongly coated regions emerge (see Figure 9, highlighted in red), which are not reflected in the experimental data.The latter fact yields room for further investigation.

CONCLUSION
The main conclusion of this current work is that the particle force effects on the surface of the substrate must be taken into consideration in order to provide better predictions of coating patterns as a basis for further process optimization.This is especially important for processes where optimization is synonymous with maximized homogeneity of coating thickness.More specifically this means that in order to retrieve high quality simulation-based predictions in terms of coating pattern distribution, coating particles must be granted the freedom to move across the substrate surface well after initial impact.Furthermore, it has been shown that the introduction of a particle-sticking-threshold criterion, based on the ratio of normal to vertical particle force components, if chosen high enough, proves beneficial in this context.The particle force balance on the substrate surface includes flow dependent forces namely i) the pressure force and ii) the fluid shear traction force, iii) the field-based electrostatic force, iv) the electric particle-particle repulsion forces as well as v) the gravitational force.
Assuming steady-state flow and a static coating setup, where the pistol and the substrate are fixed in relative position, flow dependent forces, electrostatic field-based forces and gravity forces remain stationary.However, particle-particle repulsion forces will increase with the number of particles sticking to the surface of the substrate as their integral charge value on the surface directly affects the repulsion effect.In the study carried out, it can be observed from Figure 8 that the gravitational force is indeed the most dominant tangential force, even with the higher number of particles being injected, e.g., 250'000.Within a real-life application, this relation would be reversed as normally several grams of powder particles would be injected, corresponding to a particle count ranging in the 100s million.It would not be feasible to simulate this many particles numerically, but the amount of charge calculated on the substrate surface in accordance with Equation ( 16) should be scaled to reflect the effect of a specified mass of particles being injected.Looking at the locations of the highest repulsion force in Figure 8, such a scaling could disperse the unrealistic accumulation of particles depicted in Figure 9.In this study, the number of particles being injected does not produce a significant change on the final coating pattern of the substrate as illustrated in Figure 7, since even at 250'000 computed particles, the repulsion force is lower than that of gravity in terms of magnitude, except at localized regions.Without the scaling of the charge accumulated on the surface to reflect a given mass being injected, increasing the number of particles should have a larger impact on the repulsion force and subsequently on the final coating pattern.As future work, the ratio of normal to tangential forces describing the particle sticking criterion should be carried out to best fit the entire experimental data available, where the scaling to actual injected mass is taken into consideration.

Figure 2 :
Figure2: The transfer of the particle distribution on the substrate, to its charge distribution.On the left, the particles accumulated on the substrate are depicted.The substrate surface is discretized into mainly square elements.On the right, is the particle charge distribution obtained in each element integrated from the individual particles and their surface charge density values.
4 x 0.3 m inlet and a length of 0.7 m as illustrated in Figure 4.The carrier air-flow rate in the coating pistol was set to 3 m 3 /h, the total ventilation flow of the chamber was 216 m 3 /h.Hence a volumetric flow rate of 3 m 3 /h was imposed on the pistol air-flow inlet boundary while the remaining 213 m 3 /h of flow was imposed on the chamber inlet boundary.The outlet flow was set as inlet-outlet boundary, where the flow adjusts itself based on the mass conservation.The voltage imposed on the electrode in the pistol was set as 35 kV and the pistol nozzle and the substrate were set as electrically grounded.

Figure 5 :
Figure 5: a) shows the flow field contours obtained from the simulation, while b) shows the resultant pressure and fluid shear traction magnitudes normalized by the air density; c) shows the electric field lines directed normal to the substrate surface and d) shows the resultant electric field on the substrate surface.

Figure 6 :
Figure 6: The top row shows force ratios on the plate substrate front surface in a) based on 25'000 computed particles and in b) based on 250'000 computed particles.The bottom row shows, the total normal (left) and tangential (right) force component magnitudes contributing to the respective force ratios.

Figure 8 :
Figure 8: The non-dimensionlised force components magnitude in a) the ratio of the normal to the tangential forces on the substrate, in b) the gravitational force, in c) the pressure force, in d) the fluid shear traction force, in e) the electrostatic force and in f) the repulsion force.

Figure 9 :
Figure 9: The comparison of simulation and experimental coating profiles.On the left, the simulation with the particles sticking upon impact, in the middle the experimental result, on the right the simulation with particles sticking above a force ratio of 10.