On the drying kinetics of non-spherical particle-filled polymer films: A numerical study

During the coating and drying of thin polymer-particle composites, the particle geometry has a big impact on the prediction of concentration profiles in the dry film. In this work, a plate-like geometry is used to evaluate the mass transport of the particles with the aspect ratio as a variable. The experimental determination of the viscosity and sedimentation rates allows to simulate concentration profiles in the wet film while drying. A previous simulation model was automated to describe the drying of the plate-like particles — polyvinyl alcohol-water material system using COMSOL with the initial concentration, aspect ratio, Péclet number, and sedimentation number as input parameters. The results are summarized in drying regime maps, which show an increase in the evaporation regime when the aspect ratio decreases due to lower particle mobility. This shows the importance of the geometry while predicting the particle distribution in the dry film and designing coating and drying processes.


| INTRODUCTION
The use of polymer-particle composites has gained the attention of the process engineering community in the last decade, because of the improvements that embedded particles can offer to different types of functionalized polymer films. Some examples include OLEDs, 1 fuel cell membranes, 2,3 lithium-ion batteries, 4 and biosensors. 5 During the production of said functionalized films, the location of the components along the film height is important to ensure the quality of the final product. Cardinal et al. 6 summarize the different possibilities for spherical particles to accumulate unidimensionally in charts for easy component distribution predictions based on the Péclet (Pe C ) and sedimentation (N S ) numbers, defined in Equation (1): Here, E 0 is the initial solvent evaporation rate, h 0 the initial film height, D C,0 and U C,0 are the diffusion coefficient and the sedimentation rate of a single particle.
During the drying process, the particles can either accumulate at the top (evaporation regime), sink to the bottom next to the substrate (sedimentation regime), or remain evenly distributed along the film height (diffusion regime). Corresponding drying regime maps are calculated using the conservation equations of the components and the borders are defined by the criterion that the particle concentration reaches 90% of the maximum packing concentration.
If the geometry of the particles is changed, three main aspects must be taken into consideration in order to properly model the mass transport during drying: the change of the maximum packing concentration, the orientation of the particles while drying, and the changing influence of the concentration dependency on the viscous forces acting on the particles during diffusion and/or sedimentation.
In a previous work, the influence of the geometry was studied to experimentally determine the location of the particles during drying and extrapolate the results to the dry film. It was found that flat particles align themselves with the substrate after the coating process due to shear stress, 7,8 and remain horizontal during drying and sedimentation. 9 This orientation increases the drag force, making the vertical diffusion and sedimentation of plate-like particles considerably slower in comparison to spherical particles with the same volume. The change of the particle geometry also affects the way the diffusion coefficients of single particles and the sedimentation rates are calculated. Both parameters are obtained by balancing the forces which make the particles move, including gravity, buoyancy, and viscous resistance or drag. 10 The drag coefficient is calculated by integrating the stress exerted from the fluid onto the particle surface, which it is dependent on the particle orientation as shown by Perrin. 11 The concentration dependency of sedimentation and diffusion can be modeled using expressions for the relative viscosity of polymer-solvent-systems analogous to Batchelor 12 via the sedimentation coefficient In the literature, different models can be found for the concentration dependency. In this manuscript, a general model from Russel 13 K ϕ C ð Þ¼ 1 À ϕ C ð Þ ÀK2 is utilized. This model reduces to the Batchelor 12 expression at small particle concentrations. The exponent K 2 is a main subject of the present work. Depending on the particle shape, this parameter changes its value and thus the mass transport properties of the particles.
An analytical expression for K 2 for spherical particles can be found in the literature. 14 It is derived from the creeping flow and the Oseen 15 equations as a function of the ratio between particle radius and particle center-to-center distance. The equations consider the hydrodynamic interactions and the potential between the particles, which translate to the hindrance caused by the increase of the particle concentration. A more detailed explanation regarding the modeling of K 2 can be found elsewhere. 16,14 In Table 1, there is a range of experimentally obtained K 2 for different material systems using different measurement techniques. He et al. 17 use this approach to obtain a similar function for ellipsoids dependent on the aspect ratio r p (see Figure 1) as shown in Unfortunately, this expression does not concord with experimental values from literature using different measuring techniques (see Table 1 and Figure 1). Values for non-spherical particles with different aspect ratios show a trend of lower mobility at lower aspect ratio for plate-like particles. For this reason, the literature data were fitted using Equation (3) to easily predict the sedimentation rate of plate-like particles. Expression (3) is based on the potential between two particles by modifying the expression from Weeks et al. 18 to better describe the data in Table 1, where A is the minimum value given by the spheres, B the rate of decrease and C defines a threshold value at which K 2 decreases rapidly. For a better understanding of the particle aspect ratio influence on the sedimentation rate, some exemplary particle geometries with the same volume are added as sketches to Figure 1.  Table 1. It can be seen that at r p ¼ 0:1 and lower, the particle mobility decreases dramatically. A proposed model fitted from the experimental values is plotted in green with the equation ÀK 2 r p ð Þ¼A þ B Á r p þ C ð Þ À4 À r p þ C ð Þ À2 h i and the parameters A ¼ 6:18,B ¼ 1:78 Á 10 À3 , C ¼ 5:22 Á 10 À2 . Some particles with r p at 1,0:5,0:1, and 0:03 are inserted to illustrate how the particles look like at different aspect ratios using the same volume be small and therefore it would be easier to reach the desired concentration at lower Péclet numbers. This also affects the drying regime map borders, because the discriminant of the regime zones is dependent on ϕ C,max (see Figure 4 and Equation (26)).
The aim of this work is to build on the previous simulation with quasi-binary models and expand it for different geometries, thus creating new regime maps which allow the formulation of coating suspension, the sizing of industrial dryers, and the selection of the drying conditions. The automated plotting of the regime maps is done utilizing a simple trial and error numerical method with small convergence times.
The code utilized is provided in the supporting information. The dependence of the drying regime borders on the aspect ratio of the particles is investigated, and it is discussed how the particle mobility affects the particle accumulation and the critical values of the dimensionless numbers. The aspect-ratio-dependent K 2 exponent is calculated using the approach of Equation (3) to facilitate the automatization. Finally, the influence of the aspect ratio on the drying regime maps is compared with the influence of the initial concentration of the components. consider just diffusion. The interactions between polymer and particle and solvent and particle are treated jointly as interactions between polymer-solvent-solution and particle. The equations are then written with the volume fractions of the three components ϕ C , ϕ P , and ϕ S as the variables to be solved. The subscript C is used for the particles (colloid), whereas P and S for polymer and solvent, respectively. The concentration profiles are therefore a function of space x and time t. The selected material system to simulate is silica flake particles dispersed in polyvinyl alcohol (PVA) and water. Particle-size distribution and stratification are not considered in this work.
The volume flux of the particles Γ C is composed of the diffusion D C and sedimentation velocity U C , as shown in Equation (4): The concentration-dependent particle diffusion coefficient (4) is modeled as the product of the sedimentation coefficient K ϕ C ð Þ, the derivative of the compressibility factor 26 Z ϕ C ð Þ with respect to the concentration times the concentration with a maximum packing concentration of ϕ C,max , and the temperaturedependent diffusion coefficient of a single particle D C,0 with a drag The sedimentation velocity of the particles is modeled similarly, using the sedimentation coefficient and the settling velocity of a single particle U C,0 with the density difference between the particles and the fluid, the gravity acceleration g and the particles displaced volume (6).

as shown in Equation
The change of the geometry is acknowledged by modifying the coefficients for diffusion/sedimentation and drag. For plate-like particles, ξ i is defined by Perrin's correlations, 11 which are dependent on the fluid viscosity, the aspect ratio of the particles and their orientation (either parallel or orthogonal to the transport direction). As previously mentioned, the particles are aligned with the substrate after the coating process and barely change their orientation during drying.
The different values used for the exponent K 2 were obtained using the fit of Equation (3). An example of the experimental determination of K 2 is displayed in Figure 3, where the concentration dependence of the sedimentation rates is fitted using Equation (6).
F I G U R E 2 Graphical representation of the drying of polymerparticle composites. At the top of the film the solvent evaporates; the particles in the bulk sediment and diffuse, the polymer also diffuses along the film and at the bottom of the film the substrate is considered impermeable The polymer and solvent volume fluxes are calculated using the quasi-binary approach relative to the particle movement, using the following set of equations: The diffusion coefficient of the solvent in the polymer is modeled using expression (9) according to Jeck et al. 27 The parameters A, B, and C are fitted experimentally, the values for the PVAwater system are provided below the expression. The expression has been proven to be useful in the past for describing the concentration profiles of polymer composites at other polymer molecular weight. 25 The conservation equation system is presented in Equation (10), with the Robin boundary conditions (11) to (13) x x The evaporation of the solvent is modeled with a viscous boundary layer at the phase boundary using the integrated expression (14) for _ n evap S 28 : where β S,air is the mass-transfer coefficient from the phase boundary to the top of the gas phase mass transport boundary layer.
Depending on the drying process, this coefficient can be calculated from Nusselt/Sherwood number correlations found in the literature. 29 e ρ gas is the molar density of the gas phase, e y The concentration at the phase boundary is modeled according to with U C,0 ¼ 317 μm=s at an acceleration of 12g and K 2 ¼ À44:90: Right diagram: Zero shear viscosity of polyvinyl alcohol at different concentrations and exponential fit function binary interaction parameter χ P,S is obtained from also from Jeck et al. 27 ln a S ¼ 1 À ϕ P ð Þþϕ P þ χ P,S Á ϕ 2 P : The film thickness is obtained from the mass balance at the phase boundary. It is supposed that the drying occurs only at the top of the film and the film shrinkage is one dimensional.
The non-dimensionalization of the system (10) is done with the characteristic length h 0 and the characteristic time The initial evaporation rate is calculated using expression (17).
This leads to the inclusion of the Péclet and sedimentation numbers into the equation system and thus the dimensionless equation systems (18) and (19). The transient change of the drag coefficient due to the increasing viscosity during drying is acknowledged by the binary relative viscosity η rel ϕ b P À Á based on the initial polymer concentration. for The concentration-dependence of the zero-shear viscosity of the aqueous PVA was experimentally determined and the data were fitted using the following model. The measured viscosity is plotted in Figure 3.
It is important to notice that the dimensionless evaporation rate E in Equation (21) does not require the mass-transfer coefficient to be calculated, making it independent of the air velocity.
The dimensionless boundary conditions are as it follows: The simulation is carried out with the initial concentration, temperature, aspect ratio, Pe C and N S as input parameters. The Péclet number for the polymer Pe P is obtained from the Pe C and the ratio of the diffusion coefficient of solvent in polymer and the diffusion coefficient of particles in the polymer-solvent solution.

| Numerical simulation
Similarly to a previous work, 24 Figure 4. The code utilized is given in the supporting information.

| Materials
Aqueous polymer solutions were prepared by dissolving PVA (99% hydrolyzed, Sigma-Aldrich) at a molecular weight of e M PVA ¼ 89,000 À 98,000 g=mol in distilled water (Sigma-Aldrich). The particles used were glass flakes (GF001, Glassflake Ltd.), composed mostly of SiO 2 , with an average particle size of 30 μm and a thickness between 1.0 and 1.3 μm, resulting in an aspect ratio of r p ¼ 0:03. The densities of the polymer and particles given by the manufacturer were ρ P ¼ 1190 kg=m 3 and ρ C ¼ 2600 kg=m 3 , respectively.

| Viscosity measurements
The viscosity was measured using a 20-mm cone-plate geometry at

| Sedimentation measurements
Several samples were prepared using the method described above It must be noted that the simulation cannot run properly using some combination of Pe C and N S . If it were to obtain an accumulation of polymer and particles at the top of the film, high Pe C and low N S numbers should be given to acquire such regime, due to the high diffusion velocity of the polymer in contrast to the particle. Unfortunately, to that extreme the program does not converge.
In Figure 6, a drying regime map at r p ¼ 0:10 is displayed. It is drawn utilizing the method described in the previous section, with the condition that the 90% of ϕ C,max must be reached without exceeding the overlap concentration of the polymer ϕ b P ¼ 0:13. The initial volume fractions given for the resulting maps are ϕ S,0 ¼ 0:95 and , initial concentration of each component ϕ i,0 , the particle shape r p , maximum packing concentration ϕ C,max , K 2 , and temperature. Using the input parameters, the program solves the partial differential equation system and verifies if the particle concentration at the top corresponds to 90% of ϕ C,max by calculating the error ε. If the given Pe n C difference is less than the tolerance, then the program starts again with the next N i S value, if not, a new value for the Pe n C is calculated using the secant method until it reaches convergence. The numerical approach does not change for the sedimentation regime border, only the error is calculated with the concentration at the bottom of the film ϕ C τ max ,0 ð Þ dimensionless numbers as shown in Equation (27). These numbers set the change when the evaporation dominates the diffusion Pe C,crit , the sedimentation dominates the evaporation N S,crit , and the sedimentation overcomes the diffusion Pe S,crit . Graphically the Pe S,crit can be described as the ordinate intercept of the line from the border between diffusion and sedimentation regime at high N S using Equation (28). log The influence of the aspect ratio can be seen in Figure 7. It is immediately noticed that there is a trend for a decreasing aspect ratio, which lowers the evaporation-diffusion border, thus expanding the evaporation regime for a lower aspect ratio. The map at r p ¼ 0:03 could not be fully calculated for the sedimentationdiffusion border with the same conditions as given for the other borders, that it is why the tolerance was lowered to achieve convergence.
The region in which the tolerance is lowered is marked with a dashed line.
The explanation of why such change of the borders occurs lies on the values of the exponent K 2 , which reduces with a lower aspect ratio. A lower aspect ratio at identical particle volume results in the particles getting flatter and flatter, which would mean for the fluid around the particles a greater obstacle to overcome, and thus a greater viscous resistance to flow. And when the evaporation rate increases, the particles cannot diffuse fast enough to compensate the gradient, thus expanding the evaporation regime at a low r p . As mentioned in the previous section, certain numerical values are not possible to simulate for the two main dimensionless numbers. The same principle applies for the sedimentation coefficient K ϕ C ð Þ, which at r p ¼ 0:03 is a function raised to the power of $45.
As depicted in Figure 7 the borders follow a constant value at lower N S and when they approach values in the neighborhood of the N S,crit , it starts to diverge, due to the increase in the sedimentation and diffusion driving forces. From higher N S until the N S,crit a similar behavior is observed, the border follows the line given by Equation (28) and diverge form that path when the N S lies under 10, caused F I G U R E 5 Volume fraction profiles as a function of the normalized film height at different dimensionless times, the left, middle, and right columns correspond to the concentration profiles of particles, solvent, and polymer, respectively. Calculations were done with an aspect ratio of r p ¼ 0:1. The first row corresponds to a Pe C ¼ 10 and N S ¼ 1 and an accumulation of the particles at the film surface can be observed. The row in the middle was calculated with Pe C ¼ 1 and N S ¼ 1. Even at lower Pe numbers a particle accumulation at the surface is obtained due to an increase of the viscosity. The row at the bottom corresponds to Pe C ¼ 10 and N S ¼ 10, at this high sedimentation numbers an accumulation at the bottom of the coating happens almost immediately by increment on the diffusion. The maximum in each border is the turning point for the Evaporation and Sedimentation regimes, in which 90%ϕ C,max is reached at both top and bottom of the film.
The quantification of the border shifts is shown in Figure 8, where the critical Pe C,crit and N S,crit are plotted as a function of the aspect ratio. The border between evaporation and diffusion regime corresponds to the Pe C,crit , the implicit border between evaporation and sedimentation regime is quantified with the N S,crit , and Pe S,crit defines the separation between the diffusion and sedimentation regime. In the studies from Cardinal et al. 6 it was stated that Pe S,crit is only dependent on the initial concentration of the particles. This is also shown by our results, since, all the borders at different r p overlap. The region 8log N S ð Þ 0, 1 ½ shows a different trend, because sedimentation driving forces have not yet fully overcome the evaporation and the drag provoked by the flat geometry. That it is why N S at r p > 0:03 needs to be one order of magnitude greater in order to follow the same constant trend. However, for r p ¼ 0:03 we cannot be certain, due to the limitation of the simulation routine mentioned above.
As it can be noticed in Figure 8, when the aspect ratio increases, the Péclet number follows an increasing trend, and on the other hand the sedimentation number decreases when particles get rounder. The trend of Pe C,crit can be explained by the values of K 2 , allowing a higher mobility as the aspect ratio of the plate-like particles increases (see Table 1). A higher mobility will imply higher diffusion coefficients, not allowing the particles to accumulate at the top of the film. Experimentally at r p ¼ 0:1 and r p ¼ 0:08 (see Table 1) it seems that the trend is not being followed, but those values are very close to each other and due to the different measuring techniques used to determine K 2 , it should be interpreted as a small deviation from different literature sources rather than a break of the trend.
At N S,crit , is assumed that the evaporation rate is approximately the sedimentation of the particles at x ¼ h t 0 ð Þ. In other words, F I G U R E 6 Drying regime map with initial solvent volume fraction ϕ S,0 ¼ 0:95, a particle volume fraction in the dry film ϕ C,dry ¼ 0:33, and an aspect ratio r p ¼ 0:10. The three different regimes are illustrated and marked with the corresponding letter; evaporation regime (E): the particles accumulate at the top, sedimentation regime (S): the particles accumulate at the bottom, diffusion regime (D): the particle concentration does not reach 90% of ϕ C,max and therefore the particles remain evenly distributed. The borders between regimes are marked in color using the critical dimensionless numbers as follows: D-E: regime with Pe C,crit in red, E-S: regime with N S,crit in blue, and S-D: regime with the line whose ordinate intercept corresponds to the Pe S,crit in green F I G U R E 7 Dependence of the drying regime borders on particle aspect ratio. The lower the aspect ratio r p , the lower the boundaries are located. The borders were simulated using the data available in the literature 23 and experimentally obtained. The dashed area for an r p ¼ 0:03 is due to the lower tolerance in the simulation to obtain convergence F I G U R E 8 Influence of the r p on the critical values of the Pe C and N S for plate-like particles. Showing an increase with the Pe C,crit and a decrease with the N S,crit (see Figure 6) The equation predicts that if the sedimentation coefficient K decreases (aspect ratio decreases), N S,crit will increase at a constant initial particle concentration ϕ C,0 , leading to an expansion of the evaporation regime. Similar as for Pe C,crit , this phenomenon is due to the decrease of mobility at a lower aspect ratio, making it harder for the particles to match the speed of the evaporation front and therefore the region to obtain accumulation at the top of the film is extended to higher sedimentation numbers.
Once the influence of the aspect ratio was determined, it was compared with the influence of the other input parameters, which are not dependent on the physical properties of the material system. In Figure 9, the influence of the initial volume fraction of the solvent ϕ S,0 and the influence of the volume fraction of the particles in the dry film ϕ C,dry are displayed. In the right diagram, ϕ C,max was kept constant, while varying ϕ S,0 . The volume fraction of the polymer is easily determined by adding the volume fractions of all components and equating it to one, knowing the volume fraction of the particles in the absence of a solvent. The reduction of the solvent leads to the obvious conclusion of the expansion of the evaporation regime, due to the greater initial particle concentration, thus reducing the evaporation needed to reach ϕ C,max . It must be pointed out that the change of the initial concentration changes the borders almost in the same order of magnitude as changing the aspect ratio. On the contrary, the variation of ϕ C,dry hardly changes the critical dimensionless numbers.
The main reason for the evaporation regime expansion and the low influence of the particle concentration lies in the low mobility of the plate-like particles, which intensifies the effect of the drying driving forces by raising the evaporation rate E 0 , allowing the phase boundary to move faster than the particles can diffuse.
Examples with units on the usage of the drying charts at different aspect ratios can be found in the supporting information. There, it is shown how the change in the drying conditions and particle properties, by means of the Pe C and N S , can be used to obtain a desired particle distribution in dry film.

| SUMMARY AND CONCLUSIONS
In this article, a simulation model based on past works 25,33 was developed and automated to obtain new drying regime maps for plate-like particles in water-PVA, utilizing the initial volume fractions of all components and the particle aspect ratio. The input parameters were systematically varied to quantify their influence on the drying regime borders. The component concentration profiles while drying could be calculated with the values for the particle sedimentation, which were obtained experimentally 24 and from the literature. 23,13 It was determined that these values could not be reproduced utilizing the analytical expression derived by He et al. 17 The decrease of the exponent K 2 is caused by the change of the particle aspect ratio and therefore the increase of the viscous resistance due to drag. One of the major challenges while simulating the drying regime maps is the scarce experimental data of sedimentation rates of non-spherical particles with narrow particle-size distribution. In future works, polymer particles F I G U R E 9 Influence of initial solvent and particle volume fraction in the dry film on the drying regime maps at r p ¼ 0:1. In the left diagram, the borders were drawn using a constant particle volume fraction in the dry film of ϕ C,dry ¼ 0:333 and the initial solvent volume fraction is modified, showing a reduction of the diffusion regime. This is due to the increase of the particle concentration, making it easier to fulfill the condition of the evaporation regime at lower Pe C . On the right, the ϕ C,dry is varied with a constant ϕ S,0 ¼ 0:95, which hardly produces any shift of the regime borders due to the low initial concentration ϕ C,0 with different aspect ratios will be obtained using synthesis methods with high r p purity found in the literature, 34,35 which will allow to expand the data base on sedimentation coefficients.
The regime maps were drawn as a function of the aspect ratio and the shifts of the borders were quantified using the critical values of Pe C , N S , and Pe S . The decrease of particle mobility resulted in the expansion of the evaporation regime (decrease of Pe C,crit ) and a small reduction of the sedimentation regime (increase of N S,crit ). Unfortunately, the physics of the material system limits the range of aspect ratio that could be completely simulated utilizing the method exposed in the present work. This is due to the extremely low values of K 2 .
It could appear that experimentally for the aspect ratios r p ¼ 0:1 and r p ¼ 0:08 (see Table 1) the decreasing trend inverses itself, but it must be noted that the values are close to each other and their respective K 2 s are also similar. The overall difference lies under 10% and a significant change can only be seen at logarithmic scale. Moreover, the reported values with uncertainty intervals overlap themselves, meaning that the reason why K ϕ C , r p ¼ 0:1 ð Þ < K ϕ C , r p ¼ 0:08 ð Þ lies in the different methods and experimental deviations while determining K 2 . It is important to mention, that the prediction on the particle concentration relies heavily on how the mobility of the particle is obtained. Furthermore, the influence of the aspect ratio was compared with the influence of the initial solvent volume fraction and the particle volume fraction in the dry film. It was determined that the decrease of solvent volume fraction also expands the evaporation regime, due to the lower mobility of the particles. The increase of Pe C,crit by the increasing ϕ S,0 from 0:925 to 0:95 is in the same order of magnitude or equivalent to increasing the aspect ratio r p from 0:10 to 0:15. This could be later considered while developing a drying process; instead of increasing the initial solvent composition, the aspect ratio of the particles could be changed. Variations of the particle concentration in the dried film did not produce significant changes in the same order of magnitude on the regime map borders.
The simulative results provided in the present manuscript generate new knowledge about the influence of particle geometry while drying particle-polymer composites, highlighting the connection between aspect ratio and particle mobility. These results generate more alternatives while designing drying process strategies, sizing industrial dryers, or formulating coating dispersions.