Simulation of slot-coating of nanocellulosic material subject to a wall-stress dependent slip-velocity at die-walls

Bio-based nanocellulosic materials are non-toxic, renewable, exhibit excellent barrier properties, and are suitable candidates for sustainable food packaging applications. Sizing and designing coating parameters for slot-coating process using nanocellulose suspensions is challenging due to complex shear-thinning rheology and the presence of a water-rich boundary layer, effecting significant apparent slip at the wall. Previous studies have shown that the flow inside the coating bead can be complex, with occasional stagnation regions and a rheological model incorporating yield stress which should be considered while analyzing slot coating of nanocellulosic flows. This work extends earlier investigations by including the effects of the particle depleted water-rich boundary layer. The suspension is modeled as a Casson fluid with a shear-thinning viscosity, and the particle depletion at the wall is represented by an infinitely thin layer modeled as a local shear-dependent nonlinear slip law. The resulting two-phase flow equations are solved using a Finite Volume Method (FVM) coupled with the Volume of Fluid (VoF) method for tracking the free surface interface. It is observed that slip alters the flow’s dynamics in the coating bead, and the effect of slip cannot be ignored, especially at high shear rates. For thin films, the presence of slip enhances the flow, leading to more material coated on the substrate. In contrast, for thicker coatings, apparent slip leads to an augmentation in stagnant, non-yielded regions, potentially generating uneven surfaces.


Introduction
Widespread industrial and academic interest in nanocellulosic materials 1 as a bio-based alternative for food packaging 2 is driven due to renewability, the abundance of cellulose, and proven barrier properties against oxygen, 3 mineral oils, 4 and grease. 5 From an operative perspective, a thorough understanding of the flow and rheological properties of nanocellulosic materials is essential to control and optimize high throughput manufacturing and efficient automation of coating processes. 6 In coating applications, where the desired product is a thin film (either free-standing or on a substrate), control of rheological properties is directly relevant to the final product properties. For example, in slot coating, the suspension rheology is critical for favorable coating processing conditions and will control the drying energy requirements, whereas the resulting spatial and orientation distribution of nanocellulosic particles will influence the final product characteristics. The distribution of nanocellulose fibrils in an aqueous suspension is extremely complex and leads to distinct rheological characteristics. Recent studies from our group 7,8 show the processing viability of slot-die coating and the applicability in multi-layer packaging barriers for gases and grease.
Rheological data of nanocellulose depends on the morphology of the fibrils, including fibrillation, aspectratio, length, shape, 9 and flexural properties. 10 Apart from that, the fibrils' surface charge can influence the rheology to a certain extent. 11 These factors are dictated by the process life cycle that is utilized to prepare the material. 12 In general, nanocellulosic suspensions are non-Newtonian yield-stress materials and are usually analyzed in the purview of standard models such as Herschel-Bulkley 13 or a power law model [14][15][16] . However, at high shear rates generated in slot-coating processes, the viscosity field predicted by these models may not be necessarily physically accurate. 17 One of the most prominent causes of inaccuracy in predicting the coating conditions is the presence of apparent slip at the walls, due to the presence of a thin wall-bounded particle-depleted layer on the order of a few microns where the shear is highly nonuniform. 16,18 One study 19 suggests that apparent slip may arise due to particle migration as a result of the random Brownian motion and two competing forces: hydrodynamic interactions between the fibrils and the wall, and fibril-fibril interactions, while other studies argue that fibrils migrate from the high shear rate, near-wall regions toward lower-shear bulk flow until the shearrate gradients are balanced by the concentration gradients. 20 Although the effect of apparent wall slip may not be directly relevant and can be neglected in lower concentrations, the contribution of nonlinear wall slip becomes more critical as the solid content is increased and must be addressed along with other essential modeling challenges stemming from the presence of yield stress and highly shear-thinning behavior. 21,22 Various approaches for modeling the effective deficit in momentum transfer resulting from the apparent slip layer have been suggested in the past; including molecular dynamics (MD) 23 method and Lattice Boltzmann method (LBM). 24 However, the most effective mode is to employ the traditional flow solving techniques, which discretize and solve the macroscopic Navier-Stokes (NS) equations on a continuum mechanics framework. As a first approximation, the particle-depleted layer can be modeled as an infinitely thin layer and can be effectively applied by replacing the no-slip boundary condition in the wallparallel direction by a wall-shear stress-dependent nonlinear Navier slip condition. This approach is motivated by the knowledge that, provided there is an empirical model of the slip behavior from standard experiments, the model can provide accurate stress and velocities at the wall without calculating the small-scale particle interactions, which can be prohibitively cumbersome for a slot-coating setup at industrial production rates. The current work aims to provide a study of the mechanisms on how slip can influence the flow of nanocellulosic materials in realistic industrially applicable conditions and continues our recent works in the same area. 17,25 Precisely, herein, the simulation results of slot-die flow for a 3% micro fibrillated cellulose (MFC) is presented with the knowledge of previously reported experimental data. 16 The article is organized as follows. A physical description of the slot-die coating setup, the computa-tional method, the rheology model along with the boundary conditions, are described in ''Definition of the physical problem'' section. The results of the numerical simulations are presented in ''Results and discussion'' section, followed with a summary of the work and concluding remarks.

Definition of the physical problem
In the pilot coating setup developed by Kumar et al. 26 , micro-fibrillated cellulose is pumped into the slot-die through an elongated distribution chamber, which distributes the suspension along the narrow slot gap channel to generate a uniform flow. As shown in Fig. 1, the coating head is installed at a 3 o'clock position about the backup roller and is offset (L o ) 7 mm downward from its centre. The radius (R) of the backup roller is 56.7 mm, the slot length (L s ) is 34 mm, and slot-gap (h) is kept fixed at 500 lm. Upon exiting the slot-gap, the suspension spans the gap between the adjacent die lips and the coating substrate, which is translating past with a velocity V s , attached to the backup roller. To simplify the computational setup, the region comprising the coating bead and the slot gap excluding the distribution chamber (highlighted region with dashed lines in Fig. 1) is adopted as domain of interest, and a uniform inlet velocity is prescribed at the feed-slot inlet. Both the upstream lip length (L u )  The upstream length (L uu ) is sufficiently extended to 15 mm in order to reduce the influence of artificial outlet boundary conditions into the domain, while the downstream length, L dd , is prolonged to 5 mm to assure the flow field is fully developed along the coating direction.
The minimum slot web gap (SWG) between the substrate and the downstream lip is varied, as shown in Table 1. A relative offset position of the roller creates a converging geometry in the coating gap along the downstream direction, allowing the coat weight to be controlled independently of the slot die feed flow rate by adjusting the SWG. The idea for operating the slotcoater in an unconventional metering mode is to decouple the slot inflow rate from the coat weight. However, the process in this slot-die setup is based on the forced rejection of the coating liquid, similar to knife or bar coating, and is usually classified as a postmetered method. 27,28 In other words, the unconventional metering described earlier can also be characterized as a slot-fed knife coating system. Such a setup generates high enough shear rates inside the slot to fluidize the material, thereby facilitating low process viscosities. The excess coating material falls down due to gravity along the upstream opening and is collected for reuse. If the coating layer is very thin or the drying capacity of the equipment is adequate, the entire material fed to the slot die can be converted to coating. In this case, the coat weight is controlled by the ratio of the substrate speed to the inlet velocity (V s =V in ).

Mathematical formulation
The fluid flow is considered as two-phase, laminar and isothermal. The governing equations are the conservation of mass and the unsteady Navier-Stokes equations, and is assumed to be without mass-transfer across the gas-liquid interface.
The incompressible mass conservation equation is represented as: where u is the velocity vector. The conservation of linear momentum is expressed as: Above, t, p, q, and T are the time, pressure, density, and viscous stress-tensor, respectively. Acceleration due to gravity is denoted by g. The force due to surface tension can be incorporated as a source term per unit volume and is denoted here as F b . The viscous stress tensor follows T ¼ l ru þ ru T À Á , for dynamic viscosity, expressed as l. To capture the free-surface interface, the Volume of Fluid (VoF) 29 method introduces the indicator field (a), describing volume fraction of either of the fluids, bounded by 0 a 1. The physical properties are calculated using the weighted average function. For example, if a l indicates the liquid volume fraction, and the corresponding liquid/gas density and viscosity are denoted as q l /q g and l l /l g , respectively, then the average density and viscosity are defined as The free surface of the suspension is tracked by solving the transport equation for the volume fraction as, The geometric VoF isoAdvector scheme by Roenby et al. 30 was employed for interface reconstruction due to its superior results in interface sharpness, boundedness, and shape preservation, compared to the conventional artificial surface compression method.
In equation (2), the surface tension term is considered as a volumetrically distributed force active only in the cells containing the interface as, Here r is the surface tension and j is the local curvature of the interface determined using the phasefraction field distribution as,

Rheology model
In this study, the rheology of MFC suspension is numerically described using a Casson model 31 , where shear stress (s) can be expressed as a function of shear rate (_ c) as, A Casson model is preferred here over a power law or a Herschel-Bulkley model to numerically describe the constitutive equation because the former includes the yield stress (s 0 ) and the asymptotic viscosity (g p ) at high shear rates. Furthermore, a reasonable assumption is made so that the rheology of the suspension is free of time-dependent fluctuations. The bi-viscosity regularization model 32 was employed to deal with the discontinuous material behavior of unyielding regions in the flow domain. Finite viscosity (l 0 ) is assumed when the fluid functions as a solid (s\s 0 ). As the flow yields, the local viscosity is estimated using the Casson model using equation (8).
To mimic the particle-depleted slip layer, the wall slip velocity is modeled as a function of wall shear stress, and the relationship between the two quantities is prescribed as a nonlinear power law: The parameters for the Casson model, viz. the yield stress (s 0 ¼ 255:4 Pa) and the asymptotic viscosity (g p ¼ 9:55 mPas), and for the slip model (S ¼ 3:31 and m ¼ 4:63) were computed from the rotational and slot viscometry measurements reported for 3% concentration of mechanically produced MFC. 16 These parameters were obtained by performing an optimized set of CFD studies of the slot-viscometer setup by systematically varying the slot gap, so that the differences in the experimental and predicted numerical pressure drops were minimal. From Fig. 2, we can see that the effect of nonlinear slip becomes prominent at high shear rates and the experimental data deviates from the Casson fit. Similar deviations of the experimental fitted model at high shear-region are also observed in earlier studies 16,33 Equations (2) to (7) are solved subject to the boundary conditions as given below. Here n and t denote the local normal and tangent to the boundary surface, respectively.
-BC-1: Along the velocity inlet boundary plane, a fixed uniform velocity profile (V in ) is prescribed. -BC-2: Along the artificial outflow boundary plane in the downstream, the flow is assumed to be fully developed (n Á ru ¼ 0), with fixed ambient pressure (p ¼ p 0 ). -BC-3: Along the upstream boundary plane, a direction-specific inlet-outlet velocity condition is prescribed, if the direction of flow is pointed out of the domain, a zero-gradient boundary condition is implemented, while a prescribed velocity based on the flux in the boundary normal direction is assigned in case of reverse flow. In addition, the pressure along the boundary is prescribed to be at ambient pressure (p ¼ p 0 ). -BC-4: At solid walls, for the reference no-slip case, the no-slip and no-penetration conditions are applied. (u ¼ 0), for the cases where slip is accounted for, the no-slip boundary conditions are replaced by the nonlinear slip law given in equation (9). -BC-5: For the moving substrate, the velocity component parallel to the moving wall is set to the prescribed substrate velocity (u Á t ¼ V s ), while the surface normal velocity component is set to zero (u Á n ¼ 0).

Numerical solution
The governing system of equations is solved in a segregated method employing the pressure implicit method for pressure linked equations (PIMPLE) algorithm for pressure-velocity coupling using Open-FOAM. 34 The accuracy of the implementation of the nonlinear slip boundary condition was validated against analytical solutions known in the literature. 25 For the computations, a Crank-Nicolson temporal scheme with a blending coefficient of 0.9 was adopted. For computing the gradients, a Gauss linear scheme was adopted, while the Laplacian operators were discretized using a Gauss linear corrected scheme. The convective terms in the momentum equation were estimated using a limited Gauss linear scheme, with the scheme switching to an upwind scheme in the regions of high velocity gradients. Additional details on the implementation of the numerical schemes can be found in the OpenFOAM user's guide. 35 The system of linear equations appearing in the discretization of the pressure terms was solved adopting the generalized geometric algebraic multigrid (GAMG) solver, along with the combined  DIC Gauss-Seidel smoother. A self-adapting time step based on the maximum Courant number (Co max ¼ 0:5) has been used for the time step computations. A relatively small maximum Courant number criteria ensure small time steps and guarantee that the discretization errors due to the time scheme are minimal. Additionally, to ensure the proper velocitypressure coupling, a momentum predictor was enabled with the isoAdvector scheme as it was essential for maintaining the accuracy of the simulations.
All simulations were performed at a fixed substrate velocity to inlet-velocity ratio of ðV s =V in ¼ 0:83Þ, so that the only parameters that were varied in this work are substrate velocity V s and the slot web gap (SWG), as depicted in Table 1. To clarify and compare the role of the depleted layer, all cases were run with the nonlinear slip along with the reference no-slip condition. For resolving the dynamic behavior of the moving contact line present on the backup roller, a dynamic contact angle boundary condition present in the existing framework of OpenFOAM is applied. A grid independence test for the case with slot web gap SWG = 150 lm and substrate speed V s ¼ 0:5 m/s was conducted to ensure that the numerical simulations were not dependent on the grid resolution. Three different meshes with a total element count of 50234, 21840, and 15698 were analyzed. The generated numerical results varied within 3% and were considered independent of the grid resolution. In all cases investigated, the finest mesh with 50234 elements was selected.

Results and discussion
The converging geometry ensures a more favorable pressure-gradient downstream as the flow resistance in a pressure-driven flow scales with the third power of the gap distance. The impact of slip on the pressure drop can be illustrated using Fig. 3, where pressure contours of the coating bead in the downstream gap at SWG = 600 lm are compared for the prescribed substrate speeds of V s ¼ 0:167 m/s and V s ¼ 0:5 m/s, respectively. As expected, the pressure drop is attenuated appreciably in the presence of slip, compared to the no-slip case. A similar behavior is noticed for the other slot web gaps as well.
The average pressure predicted by the numerical model at the inlet boundary indicates a measure of pressure drop inside the slot as the outlet boundaries are kept at ambient pressure. Fig. 4 depicts the relative difference in pressure drops between the corresponding slip and no-slip cases, and the effect of slip is more prominent at lower substrate speeds. As the substrate velocity is increased, the contribution from the pressure-driven flow is much less than that from the sheardriven flow. From Fig. 4, for the case V s ¼ 0:5, it can be seen that among different SWGs, the difference in inlet pressure remains small, as the suspension is already fluidized due to high shear rate. In addition to that, higher substrate speeds at narrower SWGs are found to lead to large fluctuations in the pressure fields, especially near to the moving substrate, as shown in Fig. 3b, leading to abrupt disruptions at the moving contactline in the upstream meniscus, letting air entrainment into the coating bead. Such disturbances are undesir- able as they may eventually lead to variations in the coating thickness causing coating defects. In Fig. 5, we examine the slip velocity generated at the downstream die-wall, scaled by the substrate velocity (the lip boundary from y ¼ 16 mm to y ¼ 21 mm in Fig. 3). A representative case for V s ¼ 0:5 m/s is only shown here as the behavior of the profiles are closely matched for the lower substrate velocities. The difference in slip velocities at the wall for the cases with SWG = 300 lm and SWG = 450 lm is marginal among them except near the slot exit (y [5). However, for the SWG = 150 lm we observe that the slip velocities at the gap entrance are much lower, as the shear from the substrate drives the flow entirely.
Additionally, as the pressure-driven flow dominates, the slip layer at the die-wall develops very quickly and at the die-lip exit, the velocities reach 90% of the substrate velocity, yielding a full plug-like profile. It is to be noted that for the largest SWG, the slipinfluenced region remains confined to a narrow layer as the resulting slip velocities are less than 60% of the substrate speed. This may not be desirable as some of the regions inside the gap may remain unyielded, and the flocculated structures inside the MFC suspensions continue intact, possibly leading to inhomogeneous coatings.
Understanding the development of the upward flow in the coating bead inside the converging gap enables us to further elucidate the role of slip. Figure 6 represents the non-dimensional vertical component of velocity scaled with the substrate velocity (U(y)/U(s)) as a function of the local non-dimensional horizontal position (x=x max ) sampled at three different sections in the converging geometry, viz. 1 mm and 3 mm above the slot gap exit, and at the die lip exit (at y ¼ 17 mm, y ¼ 19 mm and at y ¼ 21 mm, respectively, in Fig. 3). In the case of SWG = 150 lm, the low-viscosity zone generated alongside the moving substrate due to the shear-thinning nature of the suspension reduces its ability to transfer momentum away from the substrate. Consequently, a high-viscosity, stagnant region is formed at the downstream gap entrance along the lip wall. However, as the fluid moves along the gap, the pressure-driven flow dominates, leading to a fully developed profile at the flow exit. This, in turn, generates velocities higher than the moving substrate in the bulk of the flow.
As the SWG is increased, low shear rates mean higher momentum from the substrate is transferred across the gap by the fluid and the velocities at the lip exit are comparable to the substrate velocity. Notice that the slip velocity generated at the upstream die wall lowers the local shear stresses and leads to a more plug-like velocity profile at the exit, and the influence of the slip layer becomes limited, with just 20% of the gap for the case with SWG = 600 lm. For lower substrate velocities, the pattern appears similar and is not shown here for brevity, but the slip layer's influence is much more diminished. For instance, for the case with SWG = 600 lm and V s ¼ 0:333 m/s, the slip layer extends to 15% of the slot gap, while for the case with V s ¼ 0:167 m/s, the slip layer stretches just 8% of the slot gap. Figure 7 compares the model predictions for the final wet coating layer thickness at the lip exit for the investigated slot web gaps at different substrate speeds. The wet coating layer thicknesses for the no-slip at SWG = 150 lm is controlled by the SWG at the lip exit, except for the lowest substrate velocity, at 156 lm. However, as we employ slip at the boundary, the wet layer thickness increases by ca. 12 À 15%. For instance, at V s ¼ 0:167 m/s and 3% consistency of MFC, slip at the wall results in an increase in dry coating layer thickness from 4:7 lm to 5:3 lm. This excess material might cause runnability issues due to bleeding, i.e., accumulation of coating material outside the top lip. However, with an increase in the substrate speed, the wet-layer thickness approaches that of the SWG, for no-slip cases. Furthermore, the percentage difference between the final wet coating thickness for the slip and no-slip cases, by and large, remains constant.
Increasing the SWG to 450 lm shows a similar trend but the difference in the final wet layer thicknesses between the slip and no-slip cases is 6 À 8%, whereas for SWG = 600 lm, the influence of slip in the final coating thickness is hardly perceivable (less than 1.5%). In contrast to the larger ones, for small SWGs, the slip at the wall could considerably change the final amount of suspension coated on the substrate. Also, the marginal decrease in the coat weight with increasing substrate speed could be resulting from the weakened ability to pull the material downstream at higher shear rates. Due to the shear-thinning behavior of suspension, the downstream shear-flow induced by the moving substrate cannot overcome the pressure-driven flow.
At high consistencies, processing MFC suspension through narrow slots becomes challenging due to the presence of yield stress. We expect the flocculated fluid to be fully fluidized as it flows through the narrow slot gap. However, earlier investigations 16,36,37 suggest that there is a presence of non-deforming plugs in pipe flows and narrow slots. It will be useful to investigate non-yielded, stagnant regions (s\s 0 ) for the flows inside a coating bead. The challenge of mitigating the non-yielded regions (recalling that yield-stress strongly depends on the fiber concentration, s 0 $ OðConc 3 Þ 22 ) is even more acute when trying to increase the nanocellulose concentration to minimize the amount of water to be dried. To get a quantitative picture of the stagnant regions in the flow, the region in the downstream gap from the slot gap exit to the die lip exit is chosen, and the estimated percentage unyielded area is shown in Fig. 8. At lower SWGs, the percentage unyielded area is reduced considerably as the overall shear stress level rises in the coating bead gap bridged between the substrate and the coating head. However, at a fixed SWG, increasing the substrate speeds leads to

Conclusions
Production of nanocellulose films using slot coating is one of the favored coating methods due to its potential scalability for large-scale commercial production. However, coating very thin films at relatively large substrate speeds is challenging as the suspension exhibits complex boundary phenomena, such as particle depleted boundary layer, leading to apparent slip at the wall. In this study, a model describing the rheological behavior of nanocellulose suspension, incorporating the particle depletion at the wall is proposed. The model is implemented in an opensource CFD program and numerically modeled for an existing mini-pilot slot coating setup operating with an unconventional metering mode. The model analyzed here is for a 3% consistency suspension of microfibrillated cellulose (MFC), described as shear thinning, yield stress, Casson fluid. The resulting set of governing equations were solved using a two-phase Volume of Fluid method using OpenFOAM. The results reveal that the slip at the wall leads to an increase in the wet coating thickness, especially at smaller slot gaps. The flow dynamics in the coating bead and the evolution of the velocity profile through the coating gap strongly depend on the slot web gap. For wider slot web gaps, the slip velocities generated in the downstream gap are relatively lower and the region influenced by slip is limited. On the other hand, as the slot web gap decreases, an increased shear rate favorably improves the viscosity field downstream, mitigating unyielded regions. The numerical results suggest that the presence of slip at the wall marginally enhances the coatability of the nanocellulose suspension using the slot coater for thinner films, but may lead to an increase in unyielded regions for thicker films due to the reduced shear rates, augmenting local viscosities.
Acknowledgments KSLA, The Royal Swedish Agricultural Academy within the Tandem Forest Values supported this work via grant TFV I 2018 for the project, 'SimNano-Numerical simulation of coating processes for nanocellulosic materials.' AS thanks Pro2BE, the research environment on Processes and Products for a circular forest-based bioeconomy, hosted at Karlstad University.
Funding Open access funding provided by Abo Akademi University (ABO).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/. x (mm) 4 6 8 x (mm)