An investigation of forces on a representative surface in a pulp flow through rotating and non-rotating grooves

Softwood pulp flow in rotating and non-rotating grooves is numerically simulated in the present study to investigate the fluid flow and the forces acting on a representative surface mounted in the groove. The viscosity of softwood pulp with various consistencies is available from the measurements reported in the literature providing the opportunity to examine the effects of fiber consistency on the velocity and pressure distribution within the groove. The simulations are carried out in OpenFOAM for different values of gap thickness, angular velocity and radial positions from which the pressure coefficient and shear forces values are obtained. It is found that the shear forces within the gap increase linearly with the angular velocity for all fiber consistencies investigated and in both grooves. Also, this behavior can be successfully predicted by modeling the gap flow as a Couette flow in a two-dimensional channel. Meanwhile, a more detailed analysis of the flow kinetic energy close to the stagnation point using Bernoulli’s principle is carried out to provide a better understanding of the pressure coefficient variation with angular velocity in the non-rotating groove. A comparison of pressure coefficients obtained numerically with those calculated by considering the compression effects revealed that the comparison effects are dominating in the pulp flow within the groove.


Introduction
The forces acting on the pulp flow in the refining process are of major interest to the pulp and paper industry due to their effect on the quality of the final product. A detailed investigation of these forces provides the required information for the adjustment of the operating parameters in order to reach the desired quality. Also, such a study provides information that can be used to reduce energy consumption and increase the production rate in the process. During the refining process or more specifically in a low consistency (LC) refiner, the concentrated fiber suspension or pulp flows within the rotating and non-rotating grooves. These two types of grooves are located on the rotor and stator plates that are separated by a very small gap. By passing through these grooved plates, the fibers are mechanically treated in order to improve their bonding characteristics. Due to the high importance of this process in the pulp and paper industry, it has been the subject of different studies in the literature [1][2][3][4][5]. However, the number of studies is limited when it comes to the investigation of forces within the grooves. Martinez and Kerekes [6] derived the shear forces acting on a rotating groove indirectly from the measurements of the torque as a function of time in a single bar refiner. One result is that the shear forces increase with a decrease in the gap thickness. Prairie [7] performed an experimental investigation of the normal and shear forces on the bars in a refiner using a piezo-ceramic force sensor. Normal force distributions as a function of time showed a distinct peak at the beginning of each bar passing event followed by a sudden decrease. The peak in the shear force distribution is observed to be significantly smaller but still with a sharp decrease in a similar manner as for the normal forces. Using a piezo-ceramic force sensor similar to that in Prairie [7], Prairie et al. [8] measured the values of normal and shear forces over a range of specific edge loads (SEL) from 0.3 to 0.5 J/m. SEL is an empirical concept that is defined as the ratio of the power supplied to the refiner to the number of the bars on the rotor and stator and its rotational velocity. A wide range of variables was studied including peak normal and shear forces, coefficient of friction and shear work in each bar passing event.
It is shown that peak shear and normal forces increase with an increase in SEL. Recently, Aigner et al. [9] extended the earlier experimental works by using a high-resolution rotary encoder along with a piezo-ceramic force sensor to investigate the effect of the gap on the bar-force profile. Using the encoder, they were able to obtain the force distribution for different rotor-stator positions during each bar passing event. In addition to these experimental works, a few theoretical studies on the forces are also reported in the literature with the purpose to quantify the refining forces acting on pulp flow. Kerekes and Senger [10] theoretically characterized the shear forces on the bars as a function of SEL. Correlations for the forces acting on the pulp mass and the fiber were developed based on the geometrical parameters of fibers and bars. In his later paper, Kerekes [11] compared the results of the proposed correlations with the forces measured in different refiners and found that the predicted and measured forces are in agreement. Although being successful in different applications such as characterizing the refining action (Kerekes [12]) and quantifying the fiber shortening (Berg et al. [1]), these correlations do not provide much information on the pulp flow interaction with the groove and the gap region and its effect on the pressure distribution and the forces within the groove. The recent theoretical work of Berna et al. [13] is an attempt to categorize the parameters affecting the forces in a system of grooves. The normal forces are derived as the product of three functions, h 1 (r ) , h 2 (gap) and h 3 (pulp) . The h 2 (gap) and h 3 (pulp) originate from the fiber-fiber interaction and fibers compression in the gap region while h 1 (r ) originates from the carrying phase. Here r is radial position and ω is angular velocity. While the compressive stress associated with the h 2 (gap) and h 3 (pulp) has been extensively studied in different contexts [14][15][16], the behavior of the h 1 (r ) term is still unknown. Based on the theoretical work of Eriksen et al. [17,18] and experimental measurements of Olender et al. [19], Berna et al. [13] proposed that h 1 (r ) = water (r ) 2 . Here, it should be noted that the normal forces reported in Eriksen et al. [17,18] and Olender et al. [19] are mainly obtained in the presence of fibers and the correlation presented in Eriksen et al. [17,18] for the normal force on the bars considers specifically the effect of fibers compression. Therefore, further studies are needed to isolate the effect of the carrying phase in the absence of fibers on the h 1 (r )-term to draw an accurate conclusion for the parameters affecting the forces in this term. Motivated by this, the present study aims to investigate the effect of the carrying phase on the pressure distribution and shear forces in rotating and non-rotating grooves. Therefore, a representative surface is considered over which the pressure coefficient and shear forces are measured in different conditions. The softwood pulp is considered as working fluid and the system is modeled in open-source software, OpenFOAM. Such a study provides the opportunity to further examine the assumptions made earlier to model forces within an LC refiner and their dependency on the angular velocity and radial position. Using the numerical simulations, it is also possible to investigate the effect of gap thickness on the normal and shear forces caused by the carrying phase within the refiner. At the same time, it is attempted to present simple analytical formulations for predicting flow features specifically in the gap region and evaluate their applicability by comparing them with the results of numerical simulations. The mass and momentum equations governing the incompressible flow of softwood pulp in the stator or non-rotating groove are as: where is the apparent dynamic viscosity of pulp flow which is given in Table 1 for different consistencies, C, of fibers in the pulp [20]. Consistency is defined as the ratio of dry fibers mass to the mass of water in the pulp flow and, as can be seen in Table 1, viscosity increases with an increase in consistency. The larger fiber consistency entails a larger fiber-fiber and fiber-fluid interaction within the pulp resulting in higher viscosity, as shown in Table 1. According to the dry fiber densities ( ≈ 1500 ) reported in the literature, the average values of pulp density in the fiber consistencies studied here (1-4%) are very close to that of water, and therefore, the pulp density is considered constant in the present study. Meanwhile, a rotating frame of reference is attached to the rotor in the rotating groove simulation, as a result of which, the Coriolis inertial and centrifugal forces appear in the momentum equation as follows: where and are the rotational and position vectors, respectively. The R in Eq. (3) represents the relative velocity vector which is related to the absolute velocity as: The normal and shear forces over the representative surface are computed using the following correlations: where is density, f is the surface area vector, p is the pressure and dev is the deviatoric stress tensor. As can be seen from Eq. (5), the values of normal forces are a function of the absolute value of the pressure at the measuring point which can be different from case to case making the forces obtained from Eq. (5) case dependent. In order to generalize the results, a non-dimensional pressure coefficient is defined as: where p ref is the reference pressure within the refiner and A RS denotes the surface area of the representative surface.
It should be noted that the absolute values of normal forces are used in the present study to calculate the pressure coefficient from Eq. (7). Figure 1 illustrates the geometry of the grooves used as a numerical model for stator and rotor cavities in an LC refiner. The cross section of the groove along with the corresponding paraments is brought in Fig. 1b. The upper plate rotates with an angular velocity of in the case of non-rotating groove while the groove is stationary. In the rotating groove simulations, the groove itself rotates with an angular velocity of around the y-axis while the upper plate is fixed. Figure 1c and d illustrates, respectively, the moving parts in rotating and non-rotating cases with red outlines while the fixed parts are colored with a dark color in these figures. The dimensions of the groove cross section are similar to those in the refiner experimentally studied in Wittberg et al. [21] with a groove width of b, gap thickness of h and a groove depth H whose values in the base case simulations are presented in Table 2. The groove's length along the z-direction is set to, L, in the present study (see Fig. 1d), with an inlet and outlet located at z = 0.26 m and z = 0.3 m, respectively. In order to derive the normal and shear forces in the geometries presented in Fig. 1, a representative surface (RS) is defined on the stationary parts regardless of the case studied. The representative surface's location along the x-and y-axis is selected given the importance of the region close to the edge of the groove in the gap region in the refining process. It is in this region where the fibers get stapled and then fibrillated in an LC refiner [20]. Also, in order to investigate the effect of radial position on the forces exerted on the representative surface, simulations are carried out for different radial positions of the representative surface (see Fig. 1e) in Sect. 5.2. The Table 1 The dynamic viscosity of pulp for different fiber consistencies [20] C (%) 1. The different values of µ (Table 1) are fed into the model by changing the transport properties. Kondora and Asendrych [20] showed that using a non-Newtonian correlation for viscosity does not influence the results for the softwood pulp considered here, and therefore, they have used a Newtonian model with the viscosities presented in Table 1. The accuracy of such modeling has been confirmed also by Wittberg et al. [21] by comparing their numerical results with the experimental measurements. A second order method with explicit non-orthogonal correction is used for the discretization of diffusion terms. The convection terms are discretized using a central differencing scheme while the second order backward differencing scheme is used to discretize the unsteady terms. The Courant-Friedrich-Lewy (CFL) number in the simulations is between 0 and 1 based on the time steps. The data are sampled after the disappearance of the initial transient effects and after the flow reached steadystate conditions. The details of the numerical simulation are provided in Table 3. The total pressure boundary conditions are used at the inlet and outlet of the groove to establish a pressure difference across the grooves. The two gap faces along in the x-direction are considered periodic through which the pulp flows into the neighboring groove and no-slip boundary condition is considered over the other boundaries shown in Fig. 1.

Physical model
The simulations are carried out for four different mesh sizes in = 1000 rpm and C = 1% . It can be seen in Fig. 2a that the z-velocity component distribution along the channel height remains unchanged at x = 0.0047 m, z = 0.280 m for an increase in grid resolution after 344,700. This case has been selected since it has the highest Reynolds number ( Re b = 5490 ) among all cases studied here. The Reynolds number for the groove is defined as Re b = U av ⋅ ∕ in which the U av is the average velocity magnitude along the groove length, = 4 ⋅ A∕P wett is the characteristics length and is the kinematic viscosity of the pulp flow. Based on the values given in Tables 1 and 2, the Reynolds numbers for other cases are calculated and presented in Table 4. Due to the significance of the averaged variables such as the pressure coefficient and the shear forces in the present study, the grid independency of these variables is also studied in Fig. 2b. It can be seen that these variables only change slightly for a mesh resolution finer than 686,000 suggesting the adequacy of this grid size for further simulations.
The validation of the present solver is also carried out in earlier studies such as Jouybari et al. [22] and Jouybari et al. [23] for non-Newtonian pulp flow in channel and groove in both turbulent and laminar flows. The verification of the present results is performed by comparing the x-velocity distribution obtained using the present solver in a L × L cavity with the lid velocity of unity with the results of Erturk et al. [24] for Reynolds numbers of 1000 and 2500. As can be seen in Fig. 3a, there is a good agreement between the x-velocity along the height of the cavity calculated with the present solver and those presented in Erturk et al. [24]. Since the pulp flow within the groove can be considered as a combination of cavity and channel flows, a comparison of the friction factor calculated by the present solver with those obtained via Direct Numerical Simulation (DNS) of flow in a channel in the literature is performed in Fig. 3b. It is found that the present results are in a good agreement with the modified Blasius correlation in a channel and successfully follows the DNS results of Abe et al. [25].

The effects of fiber consistency and angular velocity
In order to disclose the effect of pulp flow rheological properties on the normal and shear forces acting on the representative surface, theoretical analysis as well as extensive sets of numerical simulations are carried out. The simulations are performed for fiber consistencies in the range of 1-4 percent and angular velocities in the range of 700-1000 rpm.   Figure 4 illustrates the pressure built up at the stagnation point at the upper right corner of the cavity. The presence of such stagnation points is of high importance in the paper industry applications where the pulp flow is passing through the grooves of LC refiner. It is at this point where the fibers get, as it is called in the literature, stapled and refined by passing the rotor grooves. The magnitude of the pressure in the stagnation point strongly affects the force that is exerted on the pulp fibers in this area. It is therefore practically important that this pressure does not exceed the maximum value that fibers or flocs can tolerate. Based on the Bernoulli principle, it can be assumed that the stagnation pressure at the corner of the cavity is the sum of dynamic pressure which is originated mainly from the rotary motion of the upper plate and the static pressure which is assumed to vary linearly between the inlet and outlet. As a result, and by introducing several assumptions, a simple relation for the stagnation pressure would be: Assumptions include that the point for which the static and dynamic pressures are measured on the right-hand side of Eq. (8) is located on the same streamline which ends up to the stagnation point and the small effect of work caused by irreversibility between this point and stagnation point. In order to fulfill the latter, it is assumed that the point is very close to the stagnation point while it is considered to be located adjacent to the upper wall so that its velocity is practically the velocity of the upper wall.

Shear forces
Since the gap thickness is so small, it can be assumed that the velocity varies linearly along the gap thickness. Considering the flow within the gap as a Couette flow in a two-dimensional parallel plate channel, Eq. (6) can be simplified as: which suggests that the shear forces vary linearly with the rotational velocity. In order to compare the accuracy of this assumption, the slope of the lines fitted to the shear forces data will be compared to that from Eq. (9) in the following. To this end, Eq. (9) can be rewritten as  in which has the unit of rpm compared to rad/s in Eq. (9). The values of m from Eq. (10) are calculated in different angular velocities and presented in Table 5.

Pulp flow in non-rotating groove
For the pulp flow simulation within the non-rotating groove, the groove in Fig. 1 is considered fixed while the upper plate is rotating with a constant angular velocity around the y-axis. The variation in pressure coefficient on the representative surface located at x = 1.75 mm, y = 6.9 mm, z = 263.5 mm is shown in Fig. 5. As can be observed, the pressure coefficient in the non-rotating groove increases with the angular velocity of the upper plate in a linear fashion. This can be attributed to the increase in the pressure at the upper right corner of the groove cross section with an increase in angular velocity of the upper plate, as illustrated in Fig. 4. The values of the pressure calculated from Eq. (8) along with those obtained at the stagnation point from the numerical simulation for a consistency of 4% are plotted in Fig. 6. The results for C = 4% are selected for comparison since, as will be shown in Fig. 7 in the following, momentum diffusion depth is the largest in this case. Therefore, the last assumption in the deriving of Eq. (8), on the proximity of fluid velocity to that of the upper wall, is better fulfilled in this case enabling a quantitative comparison in Fig. 6 in addition to the qualitative comparison of trends. In addition to the good agreement between the analytical correlation of Eq. (8)   and the numerical results, it can be observed from Eq. (8) that the stagnation pressure has a quadratic behavior in a broader range of angular velocities than that considered in the present simulations. Therefore, the linear behavior of the pressure coefficient in Fig. 5 can be due to the range of angular velocities considered here. These results are supported by the analytical and experimental observation of Berna et al. [13] for the normal forces in the refiner. An interesting observation can also be made by investigation of the pressure coefficient variation with fiber consistency in Fig. 5. As presented in Table 1, the larger the fiber consistency, the higher the fluid's viscosity which leads to higher kinematic viscosity. Therefore, the momentum diffusion from the upper plate into the groove increases as the fiber consistency increases resulting in the higher velocity of the fluid at the gap inlet and larger pressure at the upper right corner of the cavity. This increase in momentum diffusion is reflected in the thickness of diffusion depth shown with the dotted lines in Fig. 7 for different consistencies. Figure 8 illustrates the variation in shear force acting on the representative surface with the angular velocity of the upper plate and fiber consistency in the non-rotating cavity. As can be seen, the shear forces vary linearly with the upper plate angular velocity and increase with an increase in the fluid viscosity. This observation is in agreement with the theoretical correlation presented in the previous section and suggests that the flow features in the gap region are close to that in a Couette flow. Moreover, the slopes of the lines fitted to the numerical simulation in Fig. 8 are close to those calculated analytically from Eq. (10) (see Table 5) which further supports the accuracy of the theoretical analysis presented above.

Pulp flow in rotating groove
Pulp flow in the rotor is simulated assuming that the upper plate is fixed, and the groove is rotating with a constant angular velocity around the y-axis. This rotation causes two extra terms in the momentum equation, as shown in Eq. (3). In Fig. 9, the variations of C p ′ and v with fiber consistency and angular velocity in the base case simulations are plotted. In this case, the values of C p ′ in Fig. 9a are negative and they decrease with an increase in the angular velocity of the groove. The negative sign of C p ′ indicates that the pressure over the representative surface is smaller than p ref or the reference pressure/ working pressure. Additionally, it can be observed that the values of C p ′ in Fig. 9a decreases with an increase in the angular velocity. As the groove sweeps the pulp flow over the fixed upper wall, the fluid velocity increases in the gap region. Owing to this increase, the static pressure decreases causing the negative values of C p ′ . As the angular velocity increases, this decrease in static pressure is augmented, and therefore, a decreasing trend is observed for C p ′ versus in Fig. 9a. As shown in Fig. 9a the values of C p ′ are smaller for lower fiber consistencies. This can be attributed to a smaller effect of the fixed upper wall on the momentum diffusion within the groove in lower consistencies, as shown in Fig. 10, resulting in higher flow velocity adjacent to the upper wall. This smaller momentum diffusion depth is illustrated with the black dotted line in Fig. 10. Therefore, an increase in pressure drop is observed in the gap region as the groove sweeps this flow in lower consistencies. Figure 9b illustrates the variation in shear forces on the representative surface with angular velocity for different fiber consistencies. It is found that the slopes of the lines fitted to the data in Fig. 9b are close to those calculated from the theoretical correlation presented in Table 5, suggesting that the dependency of shear forces can be prescribed with the Couette flow assumption and the effect of centrifugal and Coriolis inertial forces on the flow in the gap region can be negligible specifically for the lower consistencies and small gap thickness, h = 0.2 mm.

Variation in normal and shear forces as a function of radial position
The effect of radial position, z-axis in Fig. 1, along the groove length on the normal and shear forces acting on the representative surface has been under question due to its importance on the work done on the pulp flow and the quality of final products in the pulp and paper industry. Therefore, the location of the representative surface is varied along the groove's length in the z-direction as z = 0.2745 m, z = 0.2855 m and z = 0.2965 m to examine the effect of representative surface location on the pressure coefficient and shear forces. As shown in Fig. 11a for C = 4% in the case of the non-rotating groove, the pressure coefficient increases along the groove length in the + z-direction which can be attributed to the larger tangential velocity of the plate at larger radius, and therefore, larger stagnation pressure at the upper right corner of the cavity. The values of pressure at the stagnation point from the numerical simulation are compared in Fig. 11b with those calculated using Eq. (8) for C = 4%. It can be observed that the pressure increases linearly with angular velocity in both cases despite the quadratic behavior in Eq. (8) suggesting that the linear behavior of C p ′ in Fig. 11b can be due to the scale of the problem in the radial direction. Meanwhile, the shear forces in Fig. 12 exhibit a weak deviation from the linear behavior predicted by the theoretical analysis. This suggests that the pressure variation along the gap in the x-direction slightly affects the velocity of the pulp flow in the gap, and therefore, affects the pure linear behavior of shear force. Following Eq. (9) and in the same way that Eq. (10) is derived, a linear correlation between shear forces and radial position can be obtained as: The values of m r from Eq. (11) are compared with the slopes of the lines fitted to the data in Fig. 12 and it can be seen that despite the deviation mentioned above there is an excellent agreement between these values in Table 6. Before proceeding to the variation in pressure coefficient along the grooves in the rotating cases, it is interesting to investigate its variation with angular velocity in different radial positions. As can be seen in Fig. 13, the radial position is crucial for variation in the pressure coefficient with angular velocity. The pressure coefficient increases with an increase in angular velocity at larger radii such as at z = 0.2855 m and z = 0.2965 m while an inverse trend is observed at smaller radii. This difference in behavior can be attributed to the centrifugal forces whose effect is strengthened by moving in the + z-direction. Fig. 8 The shear forces variation in the representative surface as a function of fiber consistency and angular velocity for the non-rotating groove The pressure coefficient and shear force as a function of the radial position of the representative surface in the rotating groove are inspected in Fig. 14. The results indicate a quadratic behavior for pressure coefficient with respect to the radial position which is in agreement with the previous theoretical studies of Berna et al. [13]. This behavior may be originated from the increase in the pressure due to centrifugal force that increases along the radius in a quadratic behavior. The behavior of pressure forces with respect to the radial position (z-direction in the present simulations) can also be explored through an order of magnitude analysis of the terms in the momentum equation of rotating groove in Eq. (3). Starting from the Taylor number which characterizes the importance of centrifugal forces relative to viscous forces as: which in the case of a channel or pipe is rewritten as Ta = 4 2 R 2 ∕ 2 . The R in the Taylor number shows the characteristic linear dimension and is the kinematic viscosity. The square root of Taylor number is the rotational Reynolds number which is at least of the order of 10 3 using the geometrical properties of the groove given in Table 2. Therefore, it can be concluded that the centrifugal forces are dominant on the LHS of Eq. (3) and the momentum equation reduces to:  After integration, the following correlation is obtained confirming the quadratic behavior of pressure with respect to the radial position shown by z in the present study: The quadratic behavior of the pressure distribution along the length of the groove supports the earlier numerical results presented in Fig. 14a for the pressure coefficient. Unlike the pressure coefficient, the shear forces show a linear behavior in Fig. 14b for all ω reflecting the dominancy of shear stress induced by the walls in the narrow (14) p(z) = 2 z 2 2 + constant Fig. 11 a The pressure coefficient variation in the representative surface with radius and angular velocity in the non-rotating groove and b dynamic pressure variation from Eq. (8) for C = 4% Fig. 12 The shear force variation in the representative surface with radius and angular velocity in the non-rotating groove Fig. 13 The pressure coefficient variation in the representative surface with radius and angular velocity in the rotating groove gap. However, as can be seen in Table 6, the slope of the lines fitted to the shear forces data in Fig. 14b deviates from those obtained theoretically from Eq. (11). The deviations are 29.7, 33, 33.7 and 33.7% for = 700 rpm, = 800 rpm, = 900 rpm and = 1000 rpm, respectively. Given the fact that the differences decrease with a decrease in the angular velocity, the deviation of numerical data from the theoretical calculations can be attributed to the larger effects of centrifugal forces in the rotating groove.

The effects of gap thickness
Another important parameter that affects C p ′ and v in the grooves is the thickness of the gap region, h . The variations of C p ′ and v with h in a non-rotating groove are shown in Fig. 15 where both C p ′ and v decrease with an increase in the gap thickness. A larger gap allows the flow to retain its kinetic energy as it enters the gap by deviating from the stagnation point at the corner of the groove cavity. At the same time, the velocity gradient decreases with an increase in the gap thickness resulting in lower amount of shear forces. It can also be seen in Fig. 15a that the C p ′ increases with an increase in \omega at a constant gap thickness. This can be attributed to the increase in the stagnation pressure as the upper plate's angular velocity is increased when the gap thickness is constant. Figure 16 illustrates the values of C p ′ and v calculated for the representative surface located at z = 0.2855m for four different values of h in the rotating groove. As h increases in Fig. 16a, C p ′ gradually decreases due to a decrease in the pressure at the entrance of the gap, as discussed following Fig. 15a. It can be seen that the range of variation of C p ′ increases with an increase in \omega. This can be attributed to the larger pressure drop in the gap region at higher velocities as the groove sweeps the flow over the upper wall in larger gaps. However, in the case of narrow gaps, the increase in the pressure at the stagnation point overwhelms this pressure drop resulting in a larger pressure coefficient at larger angular velocities. The variation in shear forces with the gap thickness in different angular velocities in Fig. 16b reveals that the effect of gap thickness decreases as it becomes larger so that after h = 0.2 mm the variation of v with h is significantly smaller than its variation between h = 0.1 mm and h = 0.2 mm. The variation in shear forces with the gap thickness, h, observed in Fig. 16b is in agreement with that predicted by Eq. (9). For a constant \omega, it can be rewritten as: The values of v calculated by Eq. (15) are also close to the v values obtained via numerical simulation confirming the previous conclusion on the appropriateness of the Couette flow assumption between two parallel plates to model the carrying phase within the gap.

Conclusion
Numerical simulation of softwood pulp flow in rotating and non-rotating grooves is carried out in the present study with the aim to investigate the values of pressure coefficients and shear forces on a representative surface in various operating conditions. A large number of simulations are implemented in OpenFOAM to scrutinize the effects of angular velocity, fiber consistency, radial position and (15) v = A RS 2 r (rpm) 60 ⏟⏞⏞⏞⏞⏟⏞⏞⏞⏞⏟ constant 1 h Fig. 16 a The pressure coefficient b shear force variation in the representative surface with gap thickness and angular velocity in the rotating groove gap thickness on the pressure distribution and shear forces in the gap region. With the help of theoretical analysis, it is shown that the pressure coefficient can have a quadratic behavior with respect to angular velocity in the nonrotating groove despite its linear behavior in the range of angular velocities studied here. Meanwhile, the shear forces vary linearly with the angular velocity which is in agreement with the theoretical analysis presented based on the Couette flow assumption for the flow within the gap region. Unlike the pressure coefficient in the non-rotating gap, this quantity acquires negative values in the rotating groove indicating that the pressure over the representative surface in the gap region is smaller than working pressure considered in the pressure coefficient correlation. However, the shear forces follow the same trend as that observed in the non-rotating groove. When it comes to the variation in pressure coefficients and shear forces with the radial position in the non-rotating groove, a comparison with theoretical analysis showed that the pressure coefficient can have a quadratic behavior though a linear trend is observed in the range of radial positions considered here. The shear forces revealed a weak quadratic behavior but with a slope very close to that calculated from theoretical analysis. Meanwhile, a clear quadratic behavior is observed for the variation in pressure coefficient with the radial position in the rotating groove due to the centrifugal forces that increase with an increase in radius. Investigation of the gap thickness effect showed that the pressure coefficients decrease with an increase in the gap thickness in both rotating and non-rotating grooves with different trends. This difference can be attributed to the presence of centrifugal forces acting on the pulp flow in the rotating groove. A similar trend is observed for the variation in shear forces with the gap thickness in the rotating and nonrotating cases. It is also shown that a Couette flow model in a two-dimensional channel successfully predicts the effect of gap thickness on the shear forces within the gap. Finally, a comparison between the pressure forces obtained here and those arising from the compression effects calculated using the model presented in Erikson et al. [18] shows that pressure buildup due to compressions effects is significantly larger than the pressure values calculated over the representative surface in the present simulations. The presented results also reveal that the pulp flow within the gap region behaves differently compared to the flow within the groove, and therefore, a fiber level modeling is needed to investigate the flow parameters within the gap which is the subject of future studies.