A higher-order numerical analysis to study the flow physics and to optimize the design of a short-dwell blade coaters for higher efficiency

Blade coaters are most commonly used for coating of paper and paperboard with higher efficiency. The efficiency of short-dwell blades coaters depends on many factors such as the properties of the coating material, design of the coating reservoir, the types of flow behaviour taking place inside the reservoir, etc. In this work, we have proposed an optimal design of the reservoir to improve the efficiency of short-dwell coaters. The reservoir has been modeled as flow inside a two-dimensional rectangular cavity. Incompressible Navier-Stokes equations in primitive variable formulation have been solved to obtain the flow fields inside the cavity. Spatial derivatives present in the momentum, and continuity equations are evaluated using a sixth-order accurate compact scheme whereas the temporal derivatives are calculated using the fourth-order Runge-Kutta method. The actual rate of convergence of the numerical scheme has been discussed in detail. In addition, the accuracy and stability of the used numerical method are also analysed in the spectral plane with the help of amplification factor and group velocity contour plot. The obtained numerical solutions have been validated with the existing literature. Four different aspect ratio cases (L/H = 3/4,4/3,4/5 and 5/4) have been considered for the simulations including the case of square cavity. It has been observed that L/H = 5/4 case provides best results among all others.


Introduction
Blade coaters are most commonly used in industry for paper and paperboard coating. Shortdwell blade coaters are the specifically designed blade coater used for high-speed and uniform coating. However, the trend to increase machine speed and reduction in coat weight has increased the difficulty of attaining uniform coat weight profiles on the blades. Experiments conducted with short-dwell coaters [1] indicate that the coated surface quality is highly depends on the flow inside the coating chamber. Therefor it is important to understand the flow behavior inside the chamber. In general the Reynolds numbers inside the cavity is varying from 300 to 1000. Further, Eklund, et al., in there numerical simulations [2] and studies from other researcher [3] demonstrated that two-dimensional recirculating vortices dominates inside the coating chamber. Therefor, considering the previous studies, the flow inside coater chamber has been modeled as a two-dimensional flow inside the lid driven cavity problem.

Problem Description
Flow inside the coater cavity has been modeled as a flow inside a rectangular cavity. The length and height of the cavity are varies from 0 → L and 0 → H, respectively as shown in Fig. 1. We have considered the cavity of different aspect ratios (L/H) for this problem. The grid point distributions for each aspect ration cases have been chosen in such a way that the grid spacing remain constant for every cases.

Formulation of the Problem
The unsteady incompressible Navier-Stokes equations have been solved using primitive variable formulations for the present problem. The unsteady non-dimensional incompressible momentum equations in the cartesian coordinate system can be written as, and the continuity equation, Variables appear in the Navier-Stokes equations (x, y, u, v, ρ, p) have been nondimensionalized with respect to corresponding reference quantities x L , y L , u u∞ , v v∞ , ρ ρ∞ , P ρ∞u 2 ∞ and Reynolds number Re = ρLu µ , where L is the characteristic length. The x-momentum and ymomentum equations have been used to obtain time-varying u-and v-component of the velocity field, respectively. Pressure has been obtained by solving the pressure Poisson equation, here D = ∂u ∂x + ∂v ∂y is known as dilatation.

Grid Transformation
The generalized grid transformation has been implemented to transform the physical governing equations into the computational domain. Uniform grid distribution in the computational domain is necessary for the easier implementation of a higher-order compact scheme. The generalized transformation of physical (x, y) coordinate systems to computational (ξ, η) coordinate are follows, ∂ ∂η This transformation of the governing equations has helped to solve the issue raised while solving problems with non-uniform grids in case of higher-order compact scheme.
The spatial derivatives present in the governing equations are obtained using a sixth-order compact (C6) finite-difference scheme and temporal derivatives are calculated using the fourthorder Runge-Kutta (RK4) method. Unsteady velocity fields have been obtained by solving the momentum equations and the corresponding pressure field has been calculated from the velocity field by solving the pressure Poisson equation. The SIMPLE algorithm [4] has been implemented to address the pressure velocity coupling issue raised while solving incompressible Navier-Stokes equations. BiCGSTAB algorithm [5] has been implemented in the pressure Poisson equation for faster convergence.  Accurate computation of convection dominated unsteady fluid flow problems have been solved using sixth-order compact difference scheme. The numerical properties have been evaluated considering the 1-D wave equation which is given by,

Analysis of used Numerical Method
where, c is speed of sound. Effectiveness of the used scheme has been shown here in Fig. 2. The x-and y-axis represents the non-dimensional wavenumber and ratio of numerical over exact wavenumber of the numerical methods. So, it is desirable that ratio of numerical over exact waveneumber should maintain the unit value for wide range of non-dimensional wavenumbers. Here in Fig. 2 we have compared the effectiveness of the used method (C6) with the other traditional methods (second-order central and fourth-order compact) mostly use in CFD. We have seen that used C6 method is capable calculating waveneumber up to kh ∼ 1.75 with higher accuracy.
Another two properties of the numerical method have been obtained to further establish the accuracy of the used method. Amplification factor and group velocity properties have been obtained for the used C6 method. Amplification property shows that how accurately the amplitude of a particular frequency has been captured comparing with actual physical amplitude of the signal. Whereas, the group velocity plot shows the property of numerically obtained phase velocity. More details about the theoretical and mathematical concepts are beautifully explained in literature. [9,10]  The amplification factor (|G|) and group velocity (|V g |) contour have been shown here in Fig. 2 and 3. Here, N c (= c∆t/h) in x-axis represent CFL number and kh in y-axis represent non-dimensional wavenumber. In both the contour, the desired combination of kh and N c is unit value. The used C6 method have unit |G| value for N c ≤ 0.4 for wide range of kh and |V g | contour shows a unit group velocity for kh ≤ 1.1 for N c ≤ 0.4. So, in present calculation we have considered this theoretical understanding and chosen ∆t and h values in such a way that the combination of N c and kh should falls within the neutrally stable values of |G| and |V g |.

Results and Discussions
In the present computation we have considered a variety of aspect ratios (L/H = 1/1, 3/4, 4/3, 4/5 and 5/4) and compared the velocity profiles. Time step ∆t = 0.001 has been considered for the simulations. We have validated our solver with the existing benchmark problems from the literature. The boundary conditions are,

Validation of the numerical solver
The developed numerical solver has been validate with the well known benchmark problem of lid driven cavity as shown in Fig. 1    As a benchmark problem, we have solved flow inside a lid-driven cavity. This is a classical benchmark problem used to validate numerical methods and formulation for laminar flows. Figure 1 shows a schematic of the lid-driven cavity problem. The top surface of the square cavity is moving with a unit velocity and the rest of the faces are fixed. We have performed the simulations considering three different Reynolds numbers for the present problem. Figure 4 shows streamlines of a fully-developed flow inside the cavity at Reynolds number 1000. Because of the constant velocity of the top plate and the no-slip condition, fluids inside the cavity rotates in the clockwise direction. The no-slip boundary condition has also been imposed on the cavity walls and hence the flow separation takes place near the bottom two corners of the cavity due to the viscous action. As a result, two anti-clockwise rotating vortices are observed at the bottom two corners.
We have compared the position of the centre of the three vortices (P C1 , P C2 , and P C3 ) with the literature and observed that the locations of all three vortices are analogues to the literature. The comparison with the literature has been provided in Table 1. Figure 5 shows the validation of the velocity profiles with the literature. Figure 5(a) shows the comparison of the u-component of velocity variation along x = 0 line at various Reynolds numbers for the present simulations and the velocity variation reported in the various literature. Similarly, v-velocity comparisons along y = 0 have been shown in Fig. 5(b).  cases is to understand the flow behavior for different aspect ratios. We are trying to find out the best cavity shape suitable for blade coater application. In blade coaters generally coating material from the coater releases from the top right corner of the cavity. The extrusion velocity of the coating material varies from application to application. However, higher momentum of the coating material near the top right corner helps in achieving high quality coating. The analysis of the momentum profile has been discussed later in this manuscript. Figure 6 shows the fully developed streamline pattern for various aspect ratio (L/H) cases. Although, multiple vortices are formed for all L/H cases, locations and numbers of the vortices are varies from case to case. Cavities having higher L/H(≥ 1.0) shows formation of one large primary vortex along with two secondary counter-rotating vortices, whereas lower L/H(< 1.0) cases shows only two large vortices. Table 2 gives the x and y coordinate location of the different vortices for all L/H cases. Figure 7 shows the fully developed velocity profiles along the center lines of the lid driven cavity. Figure 7 Figure 8 shows the resultant velocity profiles for all different L/H cases. Velocity profiles have been shown only for 25% of the cavity height (from top surface to 75% of total distance).  Figure 6. Streamlines for fully developed flow for different aspect ratios.

Actual order of accuracy of the numerical methods
In the absence of an exact solution, it is quite necessary to check the order of convergence of the numerical results. This is performed by taking divided differences of the average velocity values for step size h i.e. d(Avg. velocity)/dh. Figure 9 shows the decay of d(Avg. velocity)/dh as a function of h and is presented on log-log scale for Re = 1000. Here, the value of h on the x-axis is taken as the average of step sizes of the grids corresponding to the divided differences. It is very important to note that since the two graphs are almost parallel, so we can say that the slope of d(Avg. velocity)/dh tends to the slope of O(h 3 ) as h → 0. This confirms that our dividend difference results are almost third order accurate. So, the order the results provided by our proposed scheme are almost fourth-order accurate.