Analysis of the Shear Stresses in a Filling Line of Parenteral Products: The Role of Tubing

: Parenteral products appear to be sensitive to process conditions in bioprocessing steps, such as interfacial stress and shear stress. The combination of these elements is widely believed and proven to inﬂuence product stability, but the deﬁned roles of these players in the product damage process have not yet been identiﬁed. The present work addresses a current industrial problem, by focusing on the analysis of shear stress on protein-based therapeutics ﬂowing in tubing by means of Computational Fluid Dynamics simulations. The purpose of this article is not to pinpoint the mechanism triggering the damage of the product, but it represents the ﬁrst step towards wider experimental investigations and introduces a new strategy to quantify the average shear stress. The ﬁeld of scale-down approaches, used to scale the commercial process down to the laboratory level, is also explored. Since quality control is critical in the pharmaceutical realm, it is essential that the scale-down approach preserves the same stress exposure as the commercial scale, which in the present work is considered to be that resulting from shear effects. Therefore, a new approach for scaling down the commercial process is proposed, which has been compared with traditional approaches and shown to provide greater representativeness between the two scales.


Introduction
Fill-finish is widely considered the most critical step in the drug manufacturing process [1]. It involves a series of operations that can turn a purified drug product into the finished product, i.e., the final dosage form [2]. Fill-finish of protein-based therapeutics poses unique challenges, such as maintaining the stability and efficacy of the product throughout the process [3,4].
Proteins are polymer chains composed of amino acids [5], which are connected by free rotating bonds, and therefore, present many possible conformations, among which only a few stable ones grant the specificity of their biological functions [6]. Protein structural damage holds a great deal of biological importance; it can result in loss of the compact, native conformation, and, consequently, of its biological activity [2,7]. Several stressors have been identified and verified as protein denaturants, i.e., high temperatures, pressures [8], and mechanical forces [9,10]. In addition, structural destabilization appears to result from solid-liquid interfaces [11][12][13][14], air-liquid interfaces [15][16][17], and viscous stress, e.g., shear stress [18,19], which commonly occur in fill-finish manufacturing [20]. The precise role of these factors in the product damage process has not yet been universally defined and accepted, although extensively investigated [21][22][23]. According to some authors, shear stress can cause product damage regardless of the presence of interfaces [24,25]. Another theory holds that the damage is primarily due to the interface, possibly amplified by high shear stress that promotes product turnover in areas in contact with the interface [22,[26][27][28]. Furthermore, although the precise effect of exposure duration has yet to be determined, special attention should be paid to shear history. Charm and Wong [29], for example, experimentally demonstrated that catalase activity diminishes when exposure to shear increases.
Given these assumptions, it becomes paramount to precisely quantify the distribution of shear stress at a commercial scale, and control its effect during the scale-down of bioprocessing steps.
Shear stress quantifies the hydrodynamic drag force acting on a fluid [30], and is generated by its velocity gradients du dx . Shear stress has been broadly studied in the literature [31][32][33]. In the case of a Newtonian fluid under laminar conditions, the Newton's law of viscosity holds [34]: where µ is the dynamic fluid viscosity. It results in shear stress having force units per unit area (Pa). Shear stress cannot easily be measured. Some velocimetry techniques can indirectly assess the wall shear stress starting from the velocity profile measurements, relying on the Doppler shift of light scattered by particles [35]. Another popular indirect device is represented by the thermal sensor, which works on the transduction of temperature to voltage [36]. Ultimately, some probes have been developed for real-time direct measurement [37], but they show high sensitivity with respect to external noise [38], which limits their use on commercial scales. An alternative to indirect and direct measurement of the shear stress is represented by Computational Fluid Dynamics (CFD) modeling [39], which uses mathematical modeling and numerical methods to analyze and solve problems in fluid dynamics [40]. CFD is experiencing significant growth in pharmaceutical applications, as it can accelerate product and process development, facilitate optimization of existing processes, energy savings, and efficient product and process design [41].
In the biopharmaceutical field, process validation is not only a regulatory requirement, but provides an economic value; through knowledge and control of the process, batch failures are minimized [42]. In this scenario, scale-down models hold appeal in the industry; they are indispensable for process development, characterization, optimization, and validation [43]. More specifically, they allow accurate characterization of the commercial process by performing process evaluations at the laboratory level, resulting in a design space of the process that enables the control of the pharmaceutical product's Critical Quality Attributes (CQA). Such analyses are not feasible or cost-prohibitive at the commercial unit. In general, two processes are similar if they possess a geometric, kinematic, and dynamic similarity. Geometric similarity requires that linear dimensions of the original and the scaled-down system are scaled by the same ratios [34]. In turn, kinematic similarity requires geometric similarity, and the same ratio scales the velocities. Finally, dynamic similarity requires both geometric and kinematic similarity, and adds the condition that characteristic forces scale by the same ratio: specific dimensionless groups in the differential equations and boundary conditions must be the same. To provide an example, Reynolds number (Re) can be used to scale industrial processes; by definition, it is derived to describe the physical system in a scale-independent manner [44].
A common practice in pharmaceutical operations prescribes the replication in a laboratory unit of the worst-case scenario of the commercial line, in order to identify and employ a reasonable safety margin [45]. However, this method is time-and product-consuming, and it is not possible to make predictions in case of changes from the standard conditions of the line. A typical critical operating parameter in pipe flows is the average fluid velocity (u) as a function of the pump speed.
There is still a lack of in-depth understanding of the shear stresses generated when the product goes through the line's components. This could represent the first step towards further experimental investigations. In the present work, numerical modeling of the flow through a pipe was investigated under laminar and turbulent conditions. First, traditional scale-down and worst-case scenario approaches were tested numerically, aiming to study the result of their application to shear stress. Afterwards, a new scale-down approach was proposed in order to scale the commercial unit to the laboratory one while using the shear stress distribution as a scaling parameter. This last approach proved successful in ensuring representativeness between commercial and laboratory scales in terms of shear stress exposure.

Governing Equations and Theoretical Background
The continuity equation and the Navier-Stokes equation express the conservation equations of mass and momentum, and are used to estimate the pressure and velocity fields; in the case of incompressible fluids, they read as follows: where u i is the ith component of the fluid velocity, p is the fluid pressure, and ρ and ν are its density and kinematic viscosity, respectively. Laminar and turbulent flow regimes, which both occur in filling lines, were examined under steady-state conditions. The geometry of the flow is displayed in Figure 1, where a typical pipe flow is shown. There is still a lack of in-depth understanding of the shear stresses generated when the product goes through the line's components. This could represent the first step towards further experimental investigations.
In the present work, numerical modeling of the flow through a pipe was investigated under laminar and turbulent conditions. First, traditional scale-down and worst-case scenario approaches were tested numerically, aiming to study the result of their application to shear stress. Afterwards, a new scale-down approach was proposed in order to scale the commercial unit to the laboratory one while using the shear stress distribution as a scaling parameter. This last approach proved successful in ensuring representativeness between commercial and laboratory scales in terms of shear stress exposure.

Governing Equations and Theoretical Background
The continuity equation and the Navier-Stokes equation express the conservation equations of mass and momentum, and are used to estimate the pressure and velocity fields; in the case of incompressible fluids, they read as follows: = 0 (2) where is the ith component of the fluid velocity, is the fluid pressure, and and are its density and kinematic viscosity, respectively.
Laminar and turbulent flow regimes, which both occur in filling lines, were examined under steady-state conditions. The geometry of the flow is displayed in Figure 1, where a typical pipe flow is shown. For such case study, laminar flow occurs when Re is lower than 2100 while turbulent when it is greater than 4000. When dealing with turbulence, various modeling approaches can be used, named Reynolds Averaged Navier-Stokes (RANS), Large Eddy Simulation (LES), and Direct Numerical Simulation (DNS). The first approach resolves the mean flow, whereas the fluctuating components resulting from turbulence are modeled [46]. LES and DNS methods, on the other hand, are intrinsically more accurate, but entail a significant increase in the computational demand [47]. The RANS approach was therefore adopted under the given conditions. Within this framework, the − Shear Stress Transport (SST) model was used to characterize the turbulence, being a two-equation model that successfully applies in a large variety of flows [48]. It utilizes the traditional − model in the inner region of the boundary layer and switches to the − model in free shear flows, thus resulting in improved capabilities [40]. The model involves the introduction of new transport equations, namely for turbulent kinetic energy ( ) and specific turbulent dissipation rate ( ). The latter is expressed as: where is the turbulent dissipation rate. The transport equations for the turbulent properties are expressed as [48,49]: For such case study, laminar flow occurs when Re is lower than 2100 while turbulent when it is greater than 4000. When dealing with turbulence, various modeling approaches can be used, named Reynolds Averaged Navier-Stokes (RANS), Large Eddy Simulation (LES), and Direct Numerical Simulation (DNS). The first approach resolves the mean flow, whereas the fluctuating components resulting from turbulence are modeled [46]. LES and DNS methods, on the other hand, are intrinsically more accurate, but entail a significant increase in the computational demand [47]. The RANS approach was therefore adopted under the given conditions. Within this framework, the κ − ω Shear Stress Transport (SST) model was used to characterize the turbulence, being a two-equation model that successfully applies in a large variety of flows [48]. It utilizes the traditional κ − ω model in the inner region of the boundary layer and switches to the κ − ε model in free shear flows, thus resulting in improved capabilities [40]. The model involves the introduction of new transport equations, namely for turbulent kinetic energy (κ) and specific turbulent dissipation rate (ω). The latter is expressed as: where ω is the turbulent dissipation rate. The transport equations for the turbulent properties are expressed as [48,49]: ∂ω ∂t where F 1 is a blending function, ν T is the turbulent kinematic viscosity, S is the invariant measure of the strain rate, and P κ is a turbulent production term. Each regime's relevant equations are given in two separate Sections 2.1 and 2.2.

Laminar Regime
The radial velocity profile in a circular pipe under laminar conditions is well known in the literature, and depends on the distance from the center (r): where u is the average fluid velocity, and R the tubing radius.
Combining Equations (1) and (7), the shear stress profile for a Newtonian fluid also varies linearly with respect to the distance from the center, and reaches its maximum value at the tubing wall: where Q is the volumetric flowrate (m 3 s −1 ).

Turbulent Regime
Under turbulent conditions, a chaotic and random state of motion develops in which pressure and velocity undergo significant fluctuations with time [50,51]. The velocity is hence decomposed into a steady mean value with a fluctuating component [40]. The mean radial velocity can be calculated as follows [52]: where n is a function of Re [53] and u max is expressed as: Furthermore, the shear stress in such conditions is the sum of two contributions, named viscous, and turbulent shear stresses. The new turbulent term results from the Reynolds stresses and it is not straightforward to determine [54]. However, the total shear stress is still a linear function of the distance from the center and can be described as: where C f refers to the skin friction factor. This parameter can be computed starting from the expression of the mean velocity profile for the logarithmic-law correlation, known as Prandtl's friction law [50]. Nonetheless, the application of such an equation is not straightforward and, therefore, several approximations are available for defined Re ranges. Among those, the Blasius equation was chosen [55], being valid in the Re range of interest (2.1·10 3 -10 5 ) [34].

Numerical Set Up
In this paper, an analysis of the shear distribution during flowing operation was conducted thanks to CFD simulations. The finite volume method-based open-source code OpenFOAM 9 (https://github.com/OpenFOAM/OpenFOAM-9 (accessed on 1 January 2022)) was used for both laminar and turbulent regimes. The simulations were carried out An incompressible and Newtonian fluid was considered [56], and gravity was neglected. The peristaltic pump is commonly used in bioprocessing lines to drive the fluid across the process. It works by alternating compression and relaxation of the tubing, drawing content in and moving product away from the pump [57], thus creating a pulsating flow [58]. Since the time-averaged properties in pulsating flow are those of the steady flow under laminar and turbulent conditions [59], steady-state simulations were implemented rather than transient ones because of the more straightforward set up and lower computational cost needed.
The OpenFOAM utility blockMesh was utilized to build the tubing grid, whose length was chosen in order to ensure a fully developed flow under laminar and turbulent conditions, respectively [46]. A structured O-grid topology was chosen for the layout and smooth tubing was considered. Moreover, a traditional scale-down approach (same Re among the scales) and worst-case scenario approach (same u) were tested numerically, named, respectively Approach 1 and Approach 2, aiming to study their effect on shear stress. The commercial internal diameter was chosen equal to 9.53 mm while the laboratory tubing was smaller, i.e., 4.76 mm, in such a way as to decrease the product consumption. These dimensions were chosen because they are those generally used industrially for pharmaceutical processes. Regarding the velocity, 0.208 m s −1 and 1.145 m s −1 were chosen for laminar and turbulent simulations, respectively.
The domains were discretized by structured grids comprised of hexahedral cells. Geometric grid expansions were used to achieve local refinement of the mesh near the boundaries, as shown in the cross-sectional tubing area of Figure 2.

Numerical Set Up
In this paper, an analysis of the shear distribution during flowing operation was conducted thanks to CFD simulations. The finite volume method-based open-source code OpenFOAM 9 (https://github.com/OpenFOAM/OpenFOAM-9 (accessed on 1 January 2022)) was used for both laminar and turbulent regimes. The simulations were carried out on an Intel(R) Xeon(R) Gold 5118 CPU at 2.30 GHz and 128 GB of RAM. The simulations' post-processing was conducted via Python 3.8 scripts and the Paraview application.
An incompressible and Newtonian fluid was considered [56], and gravity was neglected. The peristaltic pump is commonly used in bioprocessing lines to drive the fluid across the process. It works by alternating compression and relaxation of the tubing, drawing content in and moving product away from the pump [57], thus creating a pulsating flow [58]. Since the time-averaged properties in pulsating flow are those of the steady flow under laminar and turbulent conditions [59], steady-state simulations were implemented rather than transient ones because of the more straightforward set up and lower computational cost needed.
The OpenFOAM utility blockMesh was utilized to build the tubing grid, whose length was chosen in order to ensure a fully developed flow under laminar and turbulent conditions, respectively [46]. A structured O-grid topology was chosen for the layout and smooth tubing was considered. Moreover, a traditional scale-down approach (same Re among the scales) and worst-case scenario approach (same ) were tested numerically, named, respectively Approach 1 and Approach 2, aiming to study their effect on shear stress. The commercial internal diameter was chosen equal to 9.53 mm while the laboratory tubing was smaller, i.e., 4.76 mm, in such a way as to decrease the product consumption. These dimensions were chosen because they are those generally used industrially for pharmaceutical processes. Regarding the velocity, 0.208 m s −1 and 1.145 m s −1 were chosen for laminar and turbulent simulations, respectively.
The domains were discretized by structured grids comprised of hexahedral cells. Geometric grid expansions were used to achieve local refinement of the mesh near the boundaries, as shown in the cross-sectional tubing area of Figure 2. Within the framework of the − SST model, insensitive wall functions were adopted.
is a non-dimensional scalar quantity used to define the distance from the wall as follows: [60]  Within the framework of the κ − ω SST model, y + insensitive wall functions were adopted. y + is a non-dimensional scalar quantity used to define the distance from the wall as follows: [60] y + = y u t ν (13) where y is the absolute distance from the wall, and u t the friction velocity which, in turn, is derived from the wall shear stress. Under such conditions, the boundary conditions for the turbulent properties at the wall were chosen as reported in Table 1 [49], where β is a turbulent model constant whose default value is 0.075 [49], and y w is the distance to the first cell center normal to the wall. Defined boundary conditions were specified for the other variables. At the wall, a Dirichlet boundary condition, i.e., no-slip condition, was imposed for the velocity, while at the outlet a Neumann boundary condition, i.e., zero-gradient condition, was set. A uniform velocity was then imposed at the inlet of the tubing, since the dynamics would develop once in the fully developed region. Regarding pressure, a zero-gradient condition was chosen for the wall and the inlet of the geometry, while a uniform value (equal to 0) was used for the tubing outlet; it needs to be remembered that having employed an incompressible solver, these are values of relative pressure.
A Semi-Implicit Method for Pressure Linked Equations algorithm was used through the implementation of the simpleFoam solver, applying under-relaxation factors equal to 0.9. A Geometric Agglomerated Algebraic Multigrid (GAMG) solver was implemented for pressure, while a smooth solver coupled with the Gauss-Seidel smoother was chosen for the other variables. Regarding the convergence criterion, the threshold values for the scaled residuals were set to 10 −6 for pressure and velocity. Along with the convergence criterion, monitoring of values of interest was performed to ensure they reached an accurate and steady solution. To this end, the analytical equations reported in Sections 2.1 and 2.2 were utilized for the comparison.
A grid independence analysis was conducted by increasing the cell number in order to find the optimal grid, resulting from a trade-off between computational cost and solution's accuracy. Such grid independence analysis was done by varying the number of cells in the blocks forming the O-grid topology and the number of cells in the z-direction. Results of the grid independence study for flow through a 9.53 mm diameter tubing in laminar and turbulent flows are shown in Figure 3, where the maximum shear stress was chosen as monitored parameter. For the sake of clarity, these values were validated by comparison to the analytical equations reported in Sections 2.1 and 2.2. Under laminar conditions, an observation can be drawn by comparing the meshes corresponding to approximately 3.3·10 4 and 1.1·10 5 cells (Figure 3a). Increasing the number of cells by an order of magnitude does not lead to a significant improvement in the accuracy of the maximum shear stress (i.e., 0.1%) and, therefore, the computationally less demanding mesh can be used for the analysis. Similar considerations can be drawn under turbulent conditions (Figure 3b). Around 3.3 • 10 cells were used for the laminar simulations, and 1.8 • 10 for the turbulent simulations. Given the same fluid dynamic regime, the mesh resolution between the commercial and laboratory scales was preserved. Around 3.3·10 4 cells were used for the laminar simulations, and 1.8·10 5 for the turbulent simulations. Given the same fluid dynamic regime, the mesh resolution between the commercial and laboratory scales was preserved.

Shear Stress Distribution
As described in Sections 2.1 and 2.2, shear stress under both laminar and turbulent conditions varies linearly with the distance from the tubing center. This is shown in Figure 4, where the colored lines are representative of particle's trajectories uniformly released across the diameter in the tubing's fully developed region. This provides a neat visualization of the shear stress profile through the cross-sectional area. Around 3.3 • 10 cells were used for the laminar simulations, and 1.8 • 10 for the turbulent simulations. Given the same fluid dynamic regime, the mesh resolution between the commercial and laboratory scales was preserved.

Shear Stress Distribution
As described in Sections 2.1 and 2.2, shear stress under both laminar and turbulent conditions varies linearly with the distance from the tubing center. This is shown in Figure  4, where the colored lines are representative of particle's trajectories uniformly released across the diameter in the tubing's fully developed region. This provides a neat visualization of the shear stress profile through the cross-sectional area. In bioprocessing steps, it is a common practice to refer to the maximum shear stress to characterize the hydrodynamic force acting on the product [26]. However, this phenomenon is experienced only by the fluid flowing close to the wall, as previously explained in Section 2. In addition, the fluid residence time varies through the sectional area of the tubing and reaches the lowest values at the wall, where the shear stress is the maximum. These points reveal how the effective average shear stress on the product might be In bioprocessing steps, it is a common practice to refer to the maximum shear stress to characterize the hydrodynamic force acting on the product [26]. However, this phenomenon is experienced only by the fluid flowing close to the wall, as previously explained in Section 2. In addition, the fluid residence time varies through the sectional area of the tubing and reaches the lowest values at the wall, where the shear stress is the maximum. These points reveal how the effective average shear stress on the product might be a function of the volumetric flowrate that it experiences. Not surprisingly, the volumetric flowrate varies along the radial coordinate of the tubing, reaching its maximum value at the tubing center.
Therefore, Equation (14) was developed, in which the average shear stress is obtained by summing the contributions of the local shear stress weighted on the relative volumetric flowrate. For the sake of clarity, the local properties are collected by considering each cell of a tubing slice in the fully developed region. To provide an example, Figure 5 shows a typical plot for the shear stress distribution under laminar conditions; considering an arbitrary shear threshold that is 80% of the maximum value, i.e., 1.44 ·10 −1 Pa, Figure 5 indicates that only a small percentage of the overall product, i.e., around 12%, can experience it. Sturges' rule was applied to calculate the optimum number of bins for all the histograms hereafter presented [61]. typical plot for the shear stress distribution under laminar conditions; considering an arbitrary shear threshold that is 80% of the maximum value, i.e., 1.44 • 10 Pa, Figure 5 indicates that only a small percentage of the overall product, i.e., around 12%, can experience it. Sturges' rule was applied to calculate the optimum number of bins for all the histograms hereafter presented [61].
Similar considerations can be drawn for Figure 6, representative of turbulent conditions.   Similar considerations can be drawn for Figure 6, representative of turbulent conditions. flowrate. For the sake of clarity, the local properties are collected by considering each cell of a tubing slice in the fully developed region. To provide an example, Figure 5 shows a typical plot for the shear stress distribution under laminar conditions; considering an arbitrary shear threshold that is 80% of the maximum value, i.e., 1.44 • 10 Pa, Figure 5 indicates that only a small percentage of the overall product, i.e., around 12%, can experience it. Sturges' rule was applied to calculate the optimum number of bins for all the histograms hereafter presented [61].
Similar considerations can be drawn for Figure 6, representative of turbulent conditions.

Scale-Down Approaches
Approaches 1 and 2 were numerically tested to transition from commercial to laboratory size. In addition, a new method was proposed to scale down pipe tubing using the shear stress distribution as a scaling parameter and will be further referred as Approach 3. Both laminar and turbulent regimes were investigated and are presented in two separate Sections 3.3 and 3.3.2.

Laminar Regime
As the shear stress distribution is the desired scale-down parameter, through manipulation of Equation (8), the inlet velocity at the laboratory level results from that used on the commercial scale as follows: where the suffixes c and l stand for commercial and laboratory, respectively. The tubing length is calculated as to ensure the same shear stress history, i.e., same residence time: Processes 2023, 11, 833 9 of 15 Equation (15) presents a major benefit: while scaling to the laboratory unit, i.e., decreasing the tubing diameter, the fluid velocity is also diminished. This strategy preserves the fluid dynamics regime to be preserved among the plant scales.

Turbulent Regime
The wall shear stress introduced in Equation (11) was set as a scaling parameter.
By making the skin friction factor explicit, Equation (18) shows the relation which enables the wall shear stress to be constant between commercial and laboratory units.
It is worth pointing out that the skin friction factor was calculated using Blasius approximation, as presented in Section 2.2; therefore, Equation (18) is applicable only within its validity range (2.1·10 3 -10 5 ). A more widely applicable equation could be obtained by exploiting Prandtl's equation [34], where the skin friction factor is a non-linear function of the Re; however, its application is not trivial. Furthermore, when using Equation (18), care must be taken with the characteristic Re; when decreasing the tubing size at the laboratory level, the fluid velocity decreases, and so does Re. Hence, it is essential to check the Re number of the laboratory scale to avoid working under different fluid dynamics regimes; to avoid this concern, for example, a higher tubing diameter can be selected in the laboratory.
In a similar fashion to laminar flow, the tubing length under turbulent conditions is chosen to preserve the shear history between the scale, and thus Equation (16) is followed.

Results
Approaches 1 and 2 explained in Section 1 were numerically tested under laminar and turbulent conditions. Through CFD simulations, it was found that their application does not maintain the same distribution of shear stress between the scales. Figure 7 presents the results of the application of Approach 1 under both laminar and turbulent conditions. At the laboratory level, where the diameter, and hence the radius, is half that of the commercial scale, the maximum shear stress is four times higher. Therefore, the laboratory unit would not be representative of the commercial line from the perspective of the shear stress exerted on the product. Similarly, the application of Approach 2 is shown in Figure 8. It results that under laminar conditions (Figure 8a), the maximum shear stress at the laboratory level would be overestimated by a factor of 2, being the ratio between the commercial and laboratory diameter. On the other hand, under turbulent conditions (Figure 8b), the difference between the shear stress distributions does not appear to be as pronounced; however, referring to Equations (11) and (12), this disagreement increases as the diameter chosen at the laboratory scale decreases. Similarly, the application of Approach 2 is shown in Figure 8. It results that under laminar conditions (Figure 8a), the maximum shear stress at the laboratory level would be overestimated by a factor of 2, being the ratio between the commercial and laboratory diameter. On the other hand, under turbulent conditions (Figure 8b), the difference between the shear stress distributions does not appear to be as pronounced; however, referring to Equations (11) and (12), this disagreement increases as the diameter chosen at the laboratory scale decreases. Similarly, the application of Approach 2 is shown in Figure 8. It results that under laminar conditions (Figure 8a), the maximum shear stress at the laboratory level would be overestimated by a factor of 2, being the ratio between the commercial and laboratory diameter. On the other hand, under turbulent conditions (Figure 8b), the difference between the shear stress distributions does not appear to be as pronounced; however, referring to Equations (11) and (12), this disagreement increases as the diameter chosen at the laboratory scale decreases.  Figure 9a,b; a good consistency is observed between the shear stress distributions in the two scales. Furthermore, under turbulent conditions (Figure 9c,d), the maximum shear stress is aligned between the laboratory unit and the commercial process, and so is the shear stress distribution. The slight discrepancy between the two units could be due to the fact that the average velocity at the laboratory unit is assessed using Blasius approximation and not the original equation. This choice was made to derive an equation easy to be implemented. Furthermore, under turbulent conditions (Figure 9c,d), the maximum shear stress is aligned between the laboratory unit and the commercial process, and so is the shear stress distribution. The slight discrepancy between the two units could be due to the fact that the average velocity at the laboratory unit is assessed using Blasius approximation and not the original equation. This choice was made to derive an equation easy to be implemented. Table 2 summarizes the results of the tested approaches between the commercial and laboratory scales in terms of the maximum shear stresses (σ max , Pa) and the effective average shear stress (σ, Pa) calculated according to Equation (14) using the operating conditions provided in Section 3.1. Again, it can be observed that Approach 3 is the one that best ensures representativeness between commercial and laboratory units in terms of shear stress exposure. Table 2. Comparison between the tested approaches under laminar and turbulent flow in terms of maximum shear stress (σ max , Pa) and effective average shear stress (σ, Pa) between commercial and laboratory scales.

Conclusions
In the present work, CFD modeling was used to analyze the shear stress distribution experienced by a model product when flowing through a pipe under laminar and turbulent conditions. Numerical validation of the results was performed using the well-known analytical equations for laminar and turbulent flows. It was found that the average shear stress could be effectively calculated as a function of the volumetric flowrate that it experiences. Furthermore, a new scale-down strategy was proposed and compared with the most conventional approaches. While the traditional (Approach 1) and worst-case scenario (Approach 2) approaches were not able to ensure representativeness between commercial and laboratory scales in terms of shear stress distribution, the new approach (Approach 3) was successful. From the perspective of scale-down, maintaining the fluid dynamic regime between the commercial and laboratory scales is an essential prerequisite during the application of any scale-down approach. By applying the new scale-down method, this is ensured under laminar conditions; in turbulent flow, however, a decrease in tubing diameter can lead to a decrease in Re. Therefore, it is essential to carry out a careful control of the Re at the laboratory scale and, if necessary, choose a higher tubing diameter to remain in the turbulent regime.
This work introduces a novelty in the pharmaceutical field, where scale-down is never completed considering the process from the product's perspective. However, unlike other industrial processes, product quality control is critical; therefore, it is essential that the scale-down strategy preserves the same shear stress exposure as the commercial scale, which in the present work was considered to be that resulting from shear effects.
The following steps of this study will/could involve (i) performing experimental tests to study the effect of such shear stress distribution on product stability, (ii) wider investigating the actual shear stress distribution in pulsating flow, typical of peristaltic pumps, and (iii) extending the analysis reported in this paper to all the other components commonly present in a filling line.