LES of Particle-Laden Flow in Sharp Pipe Bends with Data-Driven Predictions of Agglomerate Breakage by Wall Impacts

Abstract

The breakage of agglomerates due to wall impact within a turbulent two-phase flow is studied based on a recently developed model which relies on two artificial neural networks (ANNs). The breakup model is intended for the application within an Euler-Lagrange approach using the point-particle assumption. The ANNs were trained based on comprehensive DEM simulations. In the present study the entire simulation methodology is applied to the flow through two sharp pipe bends considering two different Reynolds numbers. In a first step, the flow structures of the continuous flow arising in both bend configurations are analyzed in detail. In a second step, the breakage behavior of agglomerates consisting of spherical, dry and cohesive silica particles is predicted based on the newly established simulation methodology taking agglomeration, fluid-induced breakage and breakage due to wall impact into account. The latter is found to be the dominant mechanism determining the resulting size distribution at the bend outlet. Since the setups are generic geometries found in dry powder inhalers, important knowledge concerning the effect of the Reynolds number as well as the design type (one-step vs. two-step deflection) can be gained.

1. Introduction

Efficient and reliable methods for predicting particle-laden flows in pipe bends are necessary for advancing different technologies. A prominent application from the pharmaceutical industry involves the development of dry powder inhalers (DPIs). These are devices used to deliver drug particles to the human respiratory tract. Such particles are typically cohesive and tend to naturally exist in form of agglomerates. However, large agglomerates are mainly intercepted in the upper airways, i.e., in the nose or the throat, and thus can not deliver drugs to the lung tissue. In order to effectively treat the central and small airways, the diameter of the inhaled particles should be within a certain range smaller than about 6 $\mathsf{\mu }$m [1]. Hence, a successful DPI design aims at delivering a fine particle size distribution by promoting intense turbulent fluid-agglomerate interactions and strong agglomerate-wall impactions to ensure adequate deagglomeration performance [2].

Experimental research has been the primary method used for evaluating the breakage behavior of agglomerates in wall-bounded turbulent flows in general and specifically in configurations relevant for DPIs (see, e.g., [2,3,4,5]). Therein, the discharged particles are typically assessed in terms of a characteristic diameter and the fraction of particles possessing sizes in the desired size range. Complementary numerical simulations are carried out to elucidate the breakage results obtained in the experiments [2,4,5,6,7]. However, the role of these simulations has been traditionally limited to predicting the Lagrangian trajectories of point particles through the turbulent flow field without taking collisions, agglomerations or breakage events into account [8]. Such a restricted approach offers basic knowledge on the nature of the fluid-particle interactions and on the character of the particle-wall impactions. However, it can not fully replace experiments, since it does not directly capture the evolution of the particle size distribution.

More advanced simulation methodologies propose the coupling between fluid flow solvers and the discrete element method (DEM) for tracking particles whilst accounting for the inter-particle and particle-wall interactions [9,10,11]. These microscale interactions are spatially and temporally resolved in DEM as functions of the overlap between the interacting surfaces mimicking their deformation. Thus, the agglomeration of particles and the breakage of agglomerates directly follows from the equation of motion of the particles allowing the prediction of the particle size development. However, this method is well-known to be computationally challenging. This is especially true for dense flow systems involving many agglomerates comprising a large number of primary particles. In addition, the computational burden immensely increases if the flow around the agglomerates structure is fully resolved by an adequately fine grid in order to obtain the most accurate description of the fluid-particle interactions.

For these reasons, considerable effort has focused on introducing and improving cost-efficient modeling techniques as viable alternatives to the most detailed methods discussed above. Reasonable CPU-times are achieved by giving up the continuous tracking of the structures of agglomerates. Thus, agglomerates are assumed to be spheres possessing mass-equivalent diameters, which allows to track agglomerates as single particles. In addition, the inter-particle and particle-wall collisions are commonly treated based on the computationally cheap hard-sphere model [12,13]. Based on first principles and on the knowledge offered by the accurate approaches, models describing phenomena like agglomeration [14,15,16], particle-wall adhesion [17,18], and agglomerate breakage due to fluid stresses [19,20] and due to wall impactions [21,22,23,24] can be utilized.

Such a recently developed data-driven wall-impact breakage model relies on artificial neural networks [23,24]. These networks were trained based on DEM results for the wall impact of dry agglomerates consisting of monodisperse silica particles. The aim of the present study is to apply the complete and efficient methodology based on the large-eddy simulation technique to investigate the deagglomeration of agglomerates in pipe bend configurations. Two sharp pipe bend geometries are considered taking two Reynolds number for each bend into account leading to four independent cases. These bend geometries and flow conditions are inspired by the experimental work of Adi et al. [5] and the corresponding numerical study of Tong et al. [10]. The overall setup is relevant for dry powder inhalers employing carrier-free powder formulation (i.e., drug only). The emphasis is put on characterizing the breakage behavior of agglomerates under different inhaler designs and flow conditions. In addition, the ability of the wall-impact breakage model to deliver important insights at manageable costs is highlighted.

The remainder of this paper is structured as follows: First, a brief theoretical background on the adopted Euler-Lagrange framework is given in Section 2. Afterwards, the wall-impact breakage is discussed in Section 3. In Section 4 the setups of the pipe bend flows are described. The results are analyzed and discussed in Section 5. Finally, the most important conclusions are summarized in Section 6.

2. Euler-Lagrange Simulation Methodology

The Eulerian-Lagrangian framework comprises the prediction of the continuous phase in the Eulerian frame of reference and the disperse phase in the Lagrangian frame of reference. In the following, the most important features of both parts are briefly described. For more details the reader is referred to previous publications (see, e.g., [13,15,16,18,19,20,25]).

2.1. Description of the Continuous Phase

Since the flow is assumed to be turbulent, the continuous phase is predicted based on the large-eddy simulation (LES) technique. For this purpose, the 3-D filtered Navier-Stokes equations for an incompressible fluid [26,27,28] are solved. The applied finite-volume scheme is formulated for curvilinear, block-structured grids applying a cell-centered, non-staggered variable arrangement. The last circumstance makes the use of the momentum interpolation technique by Rhie and Chow [29] necessary in order to avoid pressure-velocity decoupling. By approximating the surface and volume integrals by the midpoint rule and interpolating the cell face values linearly, a second-order accurate spatial discretization is achieved. The time-marching scheme relies on a classic predictor-corrector scheme, where the predictor step advances the momentum equations in time using a low-storage Runge-Kutta scheme. In the corrector step the required pressure correction is determined by the solution of the underlying Poisson equation based on an incomplete LU decomposition method [30]. The resulting time-stepping scheme is second-order accurate. For mimicking the non-resolved scales the Smagorinsky model [31] extended by the van Driest damping function for the near-wall region is chosen. The classical Smagorinsky constant $C_S=0.1$ is set. As will be shown below, the present LES predictions are wall-resolved. That means that the near-wall grid resolution is sufficiently fine to apply Stokes’ no-slip condition at solid walls which makes the use of other wall treatments superfluous. The feedback effect of the disperse phase on the fluid flow can be taken into account by the well-known particle-source-in-cell method by Crowe et al. [32]. However, in the present case the volume fraction of the particles is so low (see Section 4.4) that this effect can be neglected.

2.2. Description of the Disperse Phase

The single particles and nearly spherical agglomerates are tracked through the continuous flow field based on Newton’s second law for translational momentum. The structure of the particles and agglomerates are not resolved. Instead the point-particle approach assuming a hard-sphere model is applied. The following forces are taken into account: The drag force is determined assuming Stokes flow with the correction factor provided by Schiller and Naumann [33]. For the lift force due to velocity gradients (Saffman) the model by Crowe et al. [12] and Mei [34] is used. If the particle is rotating, its rotation is described by Newton’s second law for the angular momentum. The viscous torque acting on the particles due to rotation relative to the fluid [12,35] is taken into account. Under this circumstance the Magnus lift force has to be considered, which is achieved by the model of Rubinow and Keller [36], Crowe et al. [12] and Oesterlé and Bui Dinh [37]. Although the added-mass force typically plays a negligible role in cases in which the density of the particles is much larger than the density of the continuous phase, it is presently considered by the model of Brennen [38] and Kuerten [39]. Finally, the pressure gradient force resulting from a pressure gradient in the fluid surrounding the particle is modeled based on the model by Maxey and Riley [40]. Both the added-mass and the pressure gradient forces are not ignored in the present case because the general-purpose code includes these forces by default and the execution of the corresponding calculations only marginally increases the total computational effort. Note that gravity and buoyancy are not accounted for in the present application. Please refer to Appendix A of Breuer and Khalifa [20] for a detailed description of all forces. In addition, the effect of the unresolved scales in LES on the motion of the particles is described using the stochastic Langevin model by Breuer and Hoppe [41].

The equations of motion are integrated based on a fourth-order Runge-Kutta scheme. The first integration is carried out in the physical space providing the new velocities of the particles, whereas the second integration determining the positions of the particles on the block-structured curvilinear grid is done in the computational space. This measure avoids CPU-time intensive search algorithms, since an explicit relation between the position of the particle and the index of the cell containing the particle exists [42]. The procedure provides a highly efficient tracking scheme capable of predicting the paths of millions of particles [13,15,25,43]. Particle-wall collisions are described by an inelastic hard-sphere model including Coulomb’s law of friction [12,32,44,45]. A deterministic collision detection algorithm based on virtual cells [13,25,46] is applied. Similar to the particle-wall collisions the inter-particle collisions are described by an inelastic hard-sphere model again including Coulomb’s law of friction. Since particle-particle collisions are the prerequisite for agglomeration, these processes can be taken into account based on a momentum-based agglomeration model by Breuer and Almohammed [15]. The model takes the cohesive van-der-Waals forces acting between the collision partners described by dry, electrostatically neutral particles into account. Thus, the total collision impulse includes both impulsive and cohesive contributions. For the latter a relation based on the van-der-Waals force and the time difference of the collision intervals between the restitution and the compression phase was derived. If the cohesive force is sufficiently large, the agglomeration condition is satisfied which means that the collision partners stick together after the impact and move with a common velocity which results from the analysis of the collision event. The arising agglomerates are modeled as porous spheres defining an equivalent diameter and a packing fraction which determines the density of the agglomerate. The value of the packing fraction is obtained from the ratio of the volume of the agglomerated particles to the volume of a convex hull enwrapping the agglomerate assuming a sphere-like structure [19,20].

Replacing agglomerates by single effective spheres is a simplifying assumption, which is frequently used in the literature. For instance, in the Euler-Euler methodologies the breakage is typically described by rate functions relying on equivalent diameters. In the Euler-Lagrange context, Harshe et al. [47] investigated the motion of agglomerates with sizes varying between 10 and 1000 primary particles in a uniform shear flow taking the detailed structure of agglomerates into account. Open (dendritic) and sphere-like agglomerate structures were considered. It was demonstrated that the simplifying effective sphere approach is useful to reproduce the behavior of realistic agglomerates. Specifically, sphere-like agglomerates can be well represented by equivalent spheres, whereas ellipsoids are more suited for open structures. Another important study on this issue was presented by Dietzel and Sommerfeld [48] who performed particle-resolved lattice-Boltzmann simulations at low Reynolds numbers (i.e., Stokes regime) to investigate the drag force coefficient of agglomerates possessing different open and sphere-like structures. It was shown that the resolved drag force acting on the detailed structure can be well reproduced based on a sphere with an area equal to the projected surface area of the convex hull wrapping the agglomerated particles. The area equivalent diameter based on the convex hull resulted in drag coefficients which are nearly independent of the morphology and the orientation of the agglomerates [48,49].

From a practical point of view, the advantages of the simplifying assumption of spherical agglomerates are manifold [22,23]. For instance, computationally demanding procedures are avoided by giving up the tracking and updating of the interconnections between the particles forming an agglomerate. In addition, modeling agglomerates as single spherical entities allows to apply the same closure relationships for the fluid forces (e.g., drag and lift and their corresponding drag and lift coefficients) developed for spherical primary particles to predict the motion of agglomerates. In contrast, tracking individual particles of agglomerates requires sophisticated modifications on the relations for the fluid forces used in the context of the point-particle approximation in order to account for the shadow effect caused by closely-packed particles (see, e.g., [50] and the references therein). However, replacing agglomerates by single spheres means that the interactions between the particles forming an agglomerates (i.e., internal forces) are not taken into account. Hence, the fragmentation of agglomerates, e.g., by wall impacts, needs to be described by additional models, since it can not be directly evaluated.

Agglomerates might breakup while transported through the flow field. Besides wall-impact or collision breakage, three other mechanisms can be responsible: The drag stress, the rotary stress, and the turbulent stress. For these fluid-induced stresses the models recently proposed [19,20] are applied. To evaluate the possibility of breakup, these stresses are individually compared with the strength of the agglomerate computed based on the relation by Rumpf [51]. If breakup is detected, agglomerates are replaced by a number of fragments determined based on the responsible breakage mechanism. In addition, post-breakage velocities are assigned to the fragments to describe their motion after breakage. Furthermore, a time-lag preventing that the breakup rate is dependent on the chosen numerical time-step size is defined. All issues, i.e., the breakage condition for each individual mechanism, the resulting size distribution of the fragments, the corresponding post-breakage velocities and the time-lag, are motivated and derived based on physically meaningful arguments supported by the theory and experimental observations. Note that the adopted methodology to replace the realistic structure of an agglomerate by a sphere is also true for fragments generated either by wall-impact breakage or any of the other breakup mechanisms mentioned above. For more details on these issues, the reader is referred to Breuer and Khalifa [19,20].

Note that no particle-wall adhesion model is applied. The motivation behind this choice is to keep the analysis on the overall deagglomeration behavior focused on the interplay between breakage and agglomeration without additional complexities.

3. Wall-Impact Breakage Model

The data-driven wall-impact breakage model adopted in the present study relies on the results of extensive DEM simulations of wall-impact events of single agglomerates under controlled conditions [23,52]. To approximate the underlying relations between the impact conditions and the analyzed breakage measures included in the DEM database, two ANNs are utilized [24]. Three important tasks are tackled by the model: (i) The determination of the onset of breakage, (ii) the prediction of the resulting fragment size distribution, and (iii) the determination of the velocity vectors of the fragments arising due to breakage. For the sake of clarity, a brief description of the model is given in the following. However, the full details can be found in the original studies [23,24,52].

3.1. DEM Wall-Impact Database

The wall-impact DEM simulations are carried out using the open-source code LIGGGHTS [53]. The most important features of a wall-impact breakage event are schematically depicted in Figure 1. Agglomerates of seven size classes $N_p^{\mathrm{tot}}$ comprising 5 to 10${}^3$ dry silica particles are prepared from a large spherical DEM packing of particles [52]. To study the influence of the strength of the agglomerates, three different primary particle diameters in the range between 1 and 5 $\mathsf{\mu }$m are taken into account. The cohesion between particles is accounted for based on the van-der-Waals force model suggested by Hamaker [54] assuming dry particles. Note that the physical properties and impact parameters used in the DEM simulations are identical to those adopted in the Euler-Lagrange simulation within the hard-sphere context [23,24,52]. The computational domain used in the wall-impact DEM simulation defines a 3-D cuboid of vacuum possessing a wall-type boundary condition on two opposite sides and periodic boundary conditions on the other four sides.

To cover all possible impact scenarios, i.e., shear, oblique and normal impact, the impact angle ${\Theta }_{\mathrm{imp}}$ defined between the trajectory of the agglomerate and the impacted wall (see Figure 1) is varied between a flat 0.2${}^{°}$ and a normal 90${}^{°}$ angle. The range of the impact velocities is set to achieve the full breakage scenario from zero to full fragmentation. To ensure meaningful statistics, each impact event is repeated multiple times while changing the relative orientation of the agglomerate with respect to the wall, i.e., the position of the impact point on the surface of the agglomerate. In total about $19·{10}^3$ independent impact events are included in the database.

In order to quantify the outcomes of the DEM simulations, suitable characteristic measures are introduced. Concerning the modeling tasks (i) and (ii), two measures are introduced to describe whether breakage takes place and to determine the corresponding fragment size distribution. The first one is the fragmentation ratio [52]:

$$\mathrm{FR}={\displaystyle \frac{N_{\mathrm{frag}}^{\max }-1}{N_p^{\mathrm{tot}}-1}},$$

(1)

relating the maximum number of fragments detected during the simulation $N_{\mathrm{frag}}^{\max }$ to the maximum number of fragments an agglomerate can produce in case of full disintegration, i.e., the number of particles comprising the agglomerate $N_p^{\mathrm{tot}}$. The subtraction of unity appearing in the numerator and the denominator means that the donating agglomerate is not taken into account. Hence, FR varies in the range [0,1] expressing cases of no breakage to full fragmentation, respectively. The second measure is the fragment size parameter [52]:

$${\zeta }_i={\displaystyle \frac{N_p^{\mathrm{tot}}-{\displaystyle \sum _{j=1}^i}{N}_p^j}{N_p^{\mathrm{tot}}-i}},$$

(2)

which relates the number of primary particles which are not part of the i largest fragments to the maximum number of primary particles that can exist beside the i largest fragments in the case of a full fragmentation. However, the analysis is restricted to the three largest fragments. This means that the derived model can be used to predict the size of the three largest fragments. If more than three fragments appear and the third largest one is not a primary particle, the size of the remaining fragments is randomly set while conserving the mass of the original agglomerate and the total number of fragments predicted by FR.

Concerning the modeling task (iii) dealing with the velocity vectors of the fragments after breakage, the analysis is based on four main quantities. The first is the velocity ratio $v_{\mathrm{ratio}}=|{\mathit{v}}_{\mathrm{cm}}^{\mathrm{fr}}|/|{\mathit{v}}_{\mathrm{imp}}|$ which relates the magnitude of the velocity of the center of mass of the fragment $|{\mathit{v}}_{\mathrm{cm}}^{\mathrm{fr}}|$ to the magnitude of the impact velocity $|{\mathit{v}}_{\mathrm{imp}}|$. The inclination of the velocity vector of a fragment with respect to the impacted wall is described by the reflection angle $\alpha$, which can vary in the range [0${}^{°}$,90${}^{°}$]. The third quantity is the spreading angle $\beta$ expressing the angle of revolution of the velocity vector of the fragment ${\mathit{v}}_{\mathrm{cm}}^{\mathrm{fr}}$ with respect to the projection $\mathit{P}$ of the impact velocity vector onto the plane of the wall (see Figure 1). The fourth quantity investigated is the kinetic energy ratio:

$${\mathrm{ER}}_{\mathrm{trans}}=\frac{{\displaystyle \sum _{i=1}^{N_{\mathrm{frag}}^{\max }}}\frac{1}{2}{m}_i{v}_i^2}{\frac{1}{2}{m}_{\mathrm{ag}}{v}_{\mathrm{imp}}^2},$$

(3)

expressing the sum of the translational kinetic energy of the fragments normalized by the incident kinetic energy of the agglomerate.

Since the quantities $v_{\mathrm{ratio}}$, $\alpha$, and $\beta$ are analyzed for each fragment separately, they can be collectively described based on probability density functions (PDF). Such an analysis requires a large number of fragments under equal impact conditions. For this purpose, each impact event is repeated up to 64 times while changing the relative orientation of the agglomerate with respect to the wall. However, this large number of necessary repetitions restricted the analysis for task (iii) to a narrower set of impact conditions in comparison to tasks (i) and (ii) described above. Thus, the selected set of impact conditions takes a single primary particle diameter of $d_p=5.08$ $\mathsf{\mu }$m and three agglomerate size classes $N_p^{\mathrm{tot}}=50$, 100 and 200 into account. Nevertheless, nine impact angles are still considered and the impact velocity is varied to achieve cases with fragmentation ratios ranging between 0.05 and 1.

Figure 2 exemplarily depicts the PDFs corresponding to fragments obtained by 64 breakage events of agglomerates comprising $N_p^{\mathrm{tot}}=100$ particles and impacting the wall with an angle of ${\Theta }_{\mathrm{imp}}=0.2^{°}$. The impact velocity leads to values of FR of about 0.47 in each of the 64 events. The PDFs are approximated by the two-parameter Weibull PDF:

$$f\left(x\right)=\left\{\begin{array}{cc}\frac{k}{\lambda }{\left(\frac{x}{\lambda }\right)}^{k-1}{\mathrm{e}}^{-{(x/\lambda )}^k}\hfill & x\ge 0,\hfill \\ 0\hfill & x<0,\hfill \end{array}\right.$$

(4)

implying that for each set of impact conditions, the velocity of the fragments is modeled by means of three Weibull PDFs defining three pairs of shape k and scale $\lambda$ parameters. Note that the kinetic energy ratio ${\mathrm{ER}}_{\mathrm{trans}}$ is additionally introduced in the analysis in order to ensure the conservation of energy in the Euler-Lagrange application [23,24].

3.2. Regression by Neural Networks

The mapping between the impact conditions and the corresponding values of the breakage measures provided by the DEM database is modeled by means of two feed-forward multilayer percepton ANNs [24]. The training of both networks is performed using the software MATLAB ^® R2020a relying on the backpropagation algorithm adopting the Levenberg-Marquardt optimization [55] with Bayesian regularization (BR) [56]. The cost function is chosen as the (regularized) mean squared error (MSE).

The first network is devoted to the modeling tasks (i) and (ii), i.e., the prediction of the fragmentation ratio FR and the fragment size parameters ${\zeta }_i$. The topology of this network is 4–10–6–4 indicating two hidden layers. The number of hidden layers and the number of nodes in these layers was found after a variety of tests. The four input parameters are the impact conditions $d_p$, $N_p^{\mathrm{tot}}$, ${\Theta }_{\mathrm{imp}}$, and $v_{\mathrm{imp}}$. To facilitate the training, the values of the input parameters are scaled into the range [0,1] by linear normalization. The output parameters are FR, ${\zeta }_1$, ${\zeta }_2$, and ${\zeta }_3$. The logistic sigmoid activation function is applied in all three layers. Such a bounded function is not critical for the output layer since the output parameters of this network are naturally bounded in the interval [0,1]. As typical for training based on BR, the database is divided into training ($90\%$) and testing ($10\%$) datasets. It is noteworthy that a further division into a validation dataset is not necessary, since the generalization of the network is achieved by the Bayesian regularization instead of validation [57]. As depicted in Figure 3 the training of this network is concluded with a training MSE of $5.49·{10}^{-4}$ and a testing MSE of $5.63·{10}^{-4}$. The corresponding coefficient of determination $R^2$ for both datasets is equal to 0.997. Such results indicate a well-trained network with a high generalization ability.

The second network is assigned to the modeling task (iii), i.e., the prediction of the post-breakage velocity of the fragments. It is worth mentioning that a second network is inevitable since the analysis of the DEM results concerning the velocity of the fragments is conducted based on a restricted set of impact conditions as explained in Section 3.1. The network has a topology of 3–6–7, where the three input parameters are ${\Theta }_{\mathrm{imp}}$, $N_p^{\mathrm{tot}}$, and FR. The output parameters are the shape k and the scale $\lambda$ parameters of the three Weibull PDFs describing $v_{\mathrm{ratio}}$, $\alpha$, and $\beta$ and the fragment energy ratio ${\mathrm{ER}}_{\mathrm{trans}}$. The input parameters as well as the output parameters are normalized into the range [0,1]. The activation function in the hidden layer is the logistic sigmoid function, whereas a rectified linear activation unit (ReLU) is applied at the output layer. The database is divided into training ($87\%$) and testing ($13\%$) datasets. As depicted in Figure 4, the training of this network leads to a MSE of $8.21·{10}^{-4}$ and $8.80·{10}^{-4}$ on the training and the testing dataset, respectively. The corresponding $R^2$ values are 0.91 and 0.92. Overall, the achieved training performance is reasonable.

4. Definition of the Setup: Particle-Laden Flow through Pipe Bends

In order to study the role of impact angles on the deagglomeration performance of powders applied in dry powder inhalers (DPI), Adi et al. [5] experimentally investigated the transport of agglomerates in five different pipe bend designs for three different flow rates. These different configurations should mimic potential designs of DPIs. The objective is to provide an aerosol with a finely dispersed powder so that the active pharmaceutical ingredients can sufficiently penetrate the lower airways [58]. Agglomerate breakage due to wall impaction might help to assure the desired properties of the agglomerates. Two of these geometrical configurations are borrowed in the present study and investigated for two different flow rates. In the following, the flow configurations, the computational setups and the properties of the agglomerates are explained in detail.

4.1. Flow Configuration

The deflection of the flow in both pipe bends is 90 degrees and both pipes have a diameter of $D=1.1·{10}^{-2}\mathrm{m}$. In the first configuration the deflection is achieved by two straight intersecting pipes forming a sharp 90${}^{°}$ bend at the outer and the inner corners (see Figure 5a), whereas the second setup is composed of three straight intersecting pipes with an angle of 45${}^{°}$ at the two junctures as depicted in Figure 5b. Similar to the previously described 90${}^{°}$ pipe bend, the junctures are not rounded but sharp. The lengths of the pipe segments before the first and after the second bend are about $9D$, while the middle part has a length of approximately $2D$. For the 90${}^{°}$ bend the pipe length upstream and downstream is $10D$. Both upstream and downstream lengths are not critical since on the one hand a fully developed flow field is provided at the inlet (see Section 4.2) and on the other hand the recirculation regions in the vicinity of the bend are far away from the outlet.

The considered fluid is air at atmospheric pressure (101,325 Pa) with a constant operating temperature of 295.15 K [5]. Since the Mach number in the bends is very low, an incompressible fluid is assumed. Accordingly, the fluid properties are assumed to be constant with the following values: Density ${\rho }_f=1.196\mathrm{kg}/\mathrm{m}{}^3$ and kinematic viscosity ${\nu }_f=1.511·{10}^{-5}\mathrm{m}{}^2/\mathrm{s}$. In order to investigate the effect of the Reynolds number on the agglomerate breakage behavior in the bends, two different volumetric flow rates of 60 and 120 L/min are considered leading to Re${}_L$ = 7660 for the low Reynolds number case and Re${}_H$ = 15,320 for the high Reynolds number case, respectively. So in total four different cases (two geometries at two Re numbers based on the bulk velocity and the pipe diameter) are considered.

4.2. Computational Setup for the Fluid Flow

The computational domains are already defined in Figure 5. Since the filtered Navier-Stokes equations for an incompressible fluid are solved, boundary conditions at all boundaries (wall, outlet, inlet) are required. At the walls, which are assumed to be ideally smooth, Stokes’ no-slip condition and the impermeability condition are applied. That is possible since a wall-resolved simulation methodology is used as will be explained below. At the outlet a convective boundary condition is applied, where the convection velocity is set to the bulk velocity. At the inlet of the pipe the turbulent flow is assumed to be statistically fully developed mimicking the situation that the turbulent flow passes through a long entrance pipe before reaching the bends. In order to generate appropriate instantaneous inflow data, two supplementary simulations of pipe flows are carried out for the two Re numbers considered.

The straight pipe has a diameter D and a length of $2\pi D$. Assuming a statistically fully developed flow, periodic boundary conditions can be applied in these auxiliary simulations, which sole task is to guarantee reliable inflow data for the main simulations. One arbitrary cross-section of the pipe is chosen and the time-resolved flow velocities are stored for a time series of 900,000 time steps which comprises a time interval of 900 and 630 dimensionless time units for the low and high Re number, respectively. The prerequisite for the direct usage of the stored data of the pipe flow as inflow conditions for the bends are that identical time-step sizes and identical grid point distributions in the cross-sections of both grids are chosen. For the former a dimensionless time step of Δ$t^{*}=1.0·{10}^{-3}$ is set for the Re${}_L$ case, whereas Δ$t^{*}=7.0·{10}^{-4}$ is set at Re${}_H$.

The block-structured grids applied for the pipes and the bends have an O-type topology in the cross-section consisting of four equal outer blocks and one central block (see Figure 6c). The maximal cell spacings of the grids in the streamwise direction (Δ$x_{\max }^{+}$), the maximum circumferential spacing at the outer pipe radius (Δ${\left(r\phi \right)}_{\max }^{+}$) and the distance of the cell center of the first cell to the wall (Δ$r_{cc}^{+}$) are chosen based on the ranges suggested by Piomelli and Chasnov [59]. The selected values in wall units are 30, 15 and 0.5 for the low-Re case, whereas for the high-Re case the streamwise and wall-normal resolutions are identical and solely the circumferential spacing is slightly larger, i.e., Δ${\left(r\phi \right)}_{\max }^{+}$ = 21 instead of 15. These values are valid for both geometries investigated. Note that the grid is stretched in the radial direction towards the center of the pipe with a mild geometric stretching factor of 1.05. For the conversion of the physical units into wall units, the friction velocity $u_{\tau }={({\tau }_w/{\rho }_f)}^{0.5}$ was determined by means of an estimation of the wall-shear stress ${\tau }_w$ with the Darcy friction factor according to Moody [60].

The number of grid points in the streamwise direction, in one cross-section and the finally resulting total number of points are listed in

Table 1. Number of grid points of the block-structured grids used for LES of the bend flows.

Geometry	Re	Streamwise	Cross-Section	Total
90${}^{°}$	7660	898	4725	4.2 million
	15,320	1498	7680	11.5 million
45${}^{°}$	7660	1077	4725	5.0 million
	15,320	1637	7680	12.6 million

. For the generation of the grids, the commercial meshing software Pointwise V18.2 is used, which allows to improve the grid quality by smoothening the grid lines based on an elliptic partial differential equation. Especially for the grids of the 90${}^{°}$ pipe bends the improvement is significant since the grid lines of the initial grids had a sharp 90${}^{°}$ kink at the bend (see Figure 6). Moreover, during the smoothening process the distance of the cell centers of the first cell to the wall is kept nearly equal to the initial values. Thus, an appropriate grid resolution in the wall-normal direction is ensured also after the smoothening process. The quality of the generated grids can be demonstrated based on the equiangle skewness and the cell non-orthogonality index, which are provided directly by the meshing software Pointwise V18.2. The equiangle skewness is defined as the maximum ratio of the angle included in the cell to the angle of an equilateral element. The non-orthogonality index refers to the angle between a line connecting the centroids of the two cells and the normal of the shared face. The maximum value of this measure for all faces of a given cell is represents the cell. The distributions of both measures for all four grids are depicted in Figure 7. As expected the quality of the grids for the 45${}^{°}$ pipe bends is higher since the deflection of the flow takes place in two steps.

For parallelization purposes based on domain decomposition and MPI, the geometric block structure is mapped to a parallel block structure consisting of 126 and 280 blocks for the low-Re and high-Re case, respectively.

4.3. Properties of the Agglomerates

Spherical, dry and cohesive silica particles possessing the properties listed in

Table 2. Properties and other characteristic parameters of the considered silica particles (SiO₂).

Parameter	Unit	Value
Primary particle diameter $d_p$	m	0.97 $·{10}^{-6}$
Primary particle density ${\rho }_p$	kg·m${}^{-3}$	2000
Poisson’s ratio $\nu$	-	0.17 ${}^d$
Modulus of elasticity E	N/m${}^2$	$7.2·{10}^{10}$${}^d$
Hamaker constant H	J	$2.148·{10}^{-20}{}^a$
Min. inter-particle distance ${}_0$	m	$4.0·{10}^{-10}{}^b$
Normal restitution coefficient $e_n$	-	0.97 ${}^c$
Tangential restitution coefficient $e_t$	-	0.44 ${}^c$
Static friction coefficient ${\mu }_s$	-	0.94 ${}^c$
Kinetic friction coefficient ${\mu }_{kin}$	-	0.092 ${}^c$

${}^a$ Schubert [63], ${}^b$ Krupp [64], ${}^c$ Foerster et al. [65] for soda-lime-silica glass, ${}^d$ AZO-Materials [66].

are studied [23,52,61]. The primary particles have a diameter of 0.97 $\mathsf{\mu }$m, which corresponds to powder A in the cited references and is thus the powder with the smallest diameter leading to agglomerates of the largest strength [51,62]. This is due to the fact that the ratio of the cohesive van-der-Waals force to the inertia force increases with decreasing particle diameters. The restitution and friction coefficients in Table 2 are necessary to describe the particle-particle and particle-wall collisions as explained in Section 2.2.

Note that in the original experimental study for the investigation of the agglomerate breakage behavior [5] mannitol agglomerates were studied. Thus, the properties used here deviate from those in [5]. The reason is that the wall-impact breakage model depends on ANNs which are trained based on results obtained by DEM simulations for silica agglomerates. Therefore, a completely new training of the networks including the very CPU-time intensive data generation process based on thousands of DEM simulations would be required otherwise.

The initial agglomerates which are released into the flow field are idealized as spherical particles. This classic assumption is advantageous since the standard expressions for the fluid-particle forces and torque can be applied. In the present study, the agglomerates are composed of $N_p^{\mathrm{tot}}=1200$ primary silica particles. The equivalent diameter of the agglomerates is determined based on a mass balance taking the reduced agglomerate density by the packing fraction into account. The resulting values are $d_{\mathrm{ag}}\approx 19.17\mathsf{\mu }\mathrm{m}$ and ${\rho }_{\mathrm{ag}}\approx 311.01\mathrm{kg}/\mathrm{m}{}^3$. The latter relies on the packing fraction which depends on the assumed arrangement and properties of the primary particles. To determine the packing fraction a series of nearly spherical model agglomerates consisting of 2 to 10${}^3$ primary particles and possessing fractal dimensions close to 3 was generated by means of Monte-Carlo simulations [20]. Furthermore, the coordination number describing the average number of contacts per primary particle was evaluated in order to predict the strength of the agglomerates. The packing fraction and coordination number were rescaled according to DEM-based scaling models [67] leading to $f_{\mathrm{pack}}=0.156$ and $k_c=2$ for agglomerates consisting of 1000 primary particles. Since for large agglomerates both properties converge to (nearly) constant values, the same values are applied for the slightly larger agglomerates. More details about this structural model can be found in [19,20].

The Stokes number $St$ of the agglomerates at the moment of release is about 0.67 or 1.34 for the low and the high Reynolds number, respectively. Note that this Stokes number is built on the convective characteristic time scale ${\tau }_f=D/\left(2U_b\right)$. If the viscous time scale ${\tau }_f^{+}={\nu }_f/u_{\tau }^2$ is used instead, the corresponding values are 9.84 and and 31.53. Note that in this case, $u_{\tau }$ refers to the friction velocity computed in the periodic pipe simulation at the corresponding Re number.

4.4. Simulation Procedure

The simulations are started without releasing agglomerates into the flow to first establish a fully developed flow field. After this stage is reached, the temporal averaging of the flow is started and conducted over a sufficiently long time interval to obtain the mean flow field and the turbulent statistics. Subsequently, the agglomerates are released into the flow at the inlet of the pipes assuming the same velocity as the flow. For the low-Re cases one agglomerate is released after every 1000 time steps, while for the high-Re cases due to the smaller time step size the interval is 1500 time steps. These time increments are chosen in order to keep the dimensionless time between the released agglomerates similar for both investigated Reynolds numbers. Thereby the interactions between agglomerates released at different time steps should be kept similar. The release position is varied over the whole cross-section of the inlet, while a certain distance to the walls of the pipe is always ensured to avoid early breakage after the release. All simulations are run for the same dimensionless time interval which means that also nearly the same number of agglomerates were released in all cases, i.e., $N_{\mathrm{agg}}^{\mathrm{rel}}\approx 340$.

5. Results and Discussion

In this section the results of the Euler-Lagrange simulations are presented as follows. First, the inflow data used in the simulations of the pipe bend flows are validated in Section 5.1. Subsequently, the time-averaged flow fields of the 90${}^{°}$ and the 45${}^{°}$ pipe bends are presented in Section 5.2 for both Reynolds numbers. Finally, the breakage behavior of the agglomerates in all cases is analyzed in Section 5.3.

5.1. Validation of the Inflow Data

The aim of the simulations of the pipe flow with periodic boundary conditions is to provide fully developed inflow data for the simulations of the flow through the bends. In order to ensure suitable inflow conditions, the flow in the straight pipe is evaluated based on a comparison with experimental measurements. Durst et al. [68] conducted experiments to measure the mean velocity and the velocity fluctuations of a turbulent fully developed pipe flow. The flow in a straight pipe was investigated using laser-Doppler anemometry. Different Reynolds numbers based on the pipe diameter and the bulk velocity were taken into account. Two cases (Re = 7442 and 13,500) are comparable to those considered in the present study (Re = 7660 and 15,320).

The periodic pipe simulations are carried out for a sufficiently long time interval while averaging the flow field in time and in the streamwise direction. To facilitate the comparison with the experimental measurements, a transformation of the velocity field into a cylindrical coordinate system instead of the Cartesian coordinate system used in the predictions is done. In addition, circumferential averaging is performed.

Figure 8 compares the LES results and the experimental measurements. Figure 8a,c depict the profiles of the averaged streamwise velocity normalized by the bulk velocity $\langle u\rangle /U_b$ as a function of the normalized radius of the pipe $r/D$ for both Re${}_L$ and Re${}_H$, respectively. Although it is obvious that the profiles are symmetric, both halves are shown for the LES predictions since this is proof that the statistics are converged. As expected, the maximum streamwise velocity is found at the center of the pipe. The velocity profiles get flatter with increasing Re number leading to larger velocity gradients at the wall. Overall, a reasonable agreement between predictions and measurements is observed at both Re numbers.

Furthermore, the root-mean-squared (RMS) velocity fluctuations in streamwise direction normalized by the friction velocity ${\langle u^{\prime }\rangle }^{\mathrm{RMS}}/u_{\tau }$ are shown in Figure 8b,d for Re${}_L$ and Re${}_H$, respectively. The maximum fluctuations are found with a sharp peak in the region close to the wall. Towards the center of the pipe, the fluctuations decrease and reach a minimum at the center line. Some discrepancies between the measured and predicted peaks are recognizable especially for the Re${}_H$ case. The peaks of the streamwise velocity fluctuations are higher in the simulations than in the experiment. Moreover, it is worth mentioning that the fluctuations in the other two directions, i.e., the radial and tangential directions are additionally compared and found to be in reasonable agreement. For the sake of brevity these results are not provided here. In summary, a satisfactory agreement between the measurements and the simulation results exists ensuring that the inflow data applied in the bend flow simulations are appropriate.

5.2. Fluid Flow in Pipe Bends

The turbulent flow in sharp pipe bends is well known to exhibit separation zones at the inner and outer corners of the kink as well as secondary flow structures of the first kind [69] in form of two counter-rotating helical vortices downstream of the bend [70]. These structures, which are known as Dean vortices [71], arise due to the counteracting centrifugal force and the radial pressure gradient inside the bend. The high-velocity fluid in the core region overcomes the effect of the adverse pressure gradient owing to the high centrifugal force and moves towards the outer wall of the bend. That does not apply for the slower moving fluid, which is driven towards the inner wall leading to the development of the secondary flow. It is worth mentioning that the secondary flow was found to be unsteady, since a single Dean vortex alternately dominates the other one at a low frequency [70]. The switching between the two stable states is well known as the swirl-switching phenomenon, which was found to occur in smooth as well as sharp bend flows [72].

In the following, the predicted flow fields of all cases considered are discussed focusing on the mean velocity field and the second-order statistics of the turbulent flow.

5.2.1. Flow in the 90${}^{°}$ Pipe Bend

The time-averaged velocity magnitude normalized by the bulk velocity and the corresponding streamlines are depicted in Figure 9a,b for the sharp 90${}^{°}$ pipe bend at Re${}_L$ and Re${}_H$, respectively. In both cases the fully-developed velocity profile upstream of the bend is consistent with the flow in a straight pipe. Further downstream an adverse pressure gradient in the streamwise direction develops towards the bends (not depicted here) leading to the growth of the boundary layer thickness at the outer side of the bend before the flow finally detaches. The separation of the flow gives rise to recirculation zones at the outer corners consisting of a series of asymmetric, counter-rotating vortices. In the present case two counter-rotating vortices could be resolved as depicted in the upper zoom in Figure 9a. Consequently, the flow cross-section is tapered, which causes the flow to be accelerated at the inner kink to globally conserve mass.

Inside the bend the deflected flow is driven to the outer side of the bend by the centrifugal forces. Therefore, the flow is accelerated in the axial direction leading to the reattachment of the boundary layer at the outer wall. At the inner side of the bend the flow detaches since it can not follow the sharp bend. Consequently, a large separation zone is established on the inner wall directly after the bend. This zone is occupied by slowly moving fluid and possesses a reverse-flow region in its core. Obviously, a strong shear layer arises between the fast central region and the recirculation zone.

In Figure 10 the contours of the downstream streamwise velocity $\langle v\rangle /U_b$ in a cross-section one pipe diameter downstream of the bend are shown. Additionally, the streamlines of the secondary flow are superimposed. As visible, the lowest streamwise velocity is located at the center of the separation zone. The Dean vortices, which are rotating around the two pressure minima, are nearly symmetric for the time-averaged flow field. That is a hint that the time-averaging period was sufficiently long.

The resolved turbulent kinetic energy $\langle k\rangle =1/2(\langle u^{\prime }{u}^{\prime }\rangle +\langle v^{\prime }{v}^{\prime }\rangle +\langle w^{\prime }{w}^{\prime }\rangle )$ normalized by $U_b^2$ is shown in Figure 11a,b for the Re${}_L$ and the Re${}_H$ cases, respectively. As expected, the strongest velocity fluctuations are found in the shear layer between the recirculation region and the core flow, and inside the recirculation region. Moreover, the resolved Reynolds shear stress $\langle u^{\prime }{v}^{\prime }\rangle /U_b^2$ is depicted in Figure 12a,b. The distributions can be attributed to the strong velocity gradients at the edge of the shear layer, where the fluid velocity increases strongly towards the outer wall. As visible in Figure 9 showing the velocity magnitude in the symmetry plane as well as Figure 10 depicting the streamwise velocity in a cross-section downstream on the bend, the velocity gradients also remain high a few diameters downstream of the bends. That explains why the Reynolds shear stress remains high far behind the bends.

In general, the results obtained at the two Reynolds numbers considered are qualitatively very similar. Moreover, the described flow patterns are in close agreement with the experimental observations [70,73] and numerical results [74,75] reported in the literature concerning sharp 90${}^{°}$ bend flows. However, a direct comparison can not be carried out, since on the one hand a different range of Re (up to $\mathrm{Re}\approx 1.2·{10}^5$) was investigated in these studies. On the other hand, a sharp duct bend with a square cross-section was taken into account in the studies employing comparable Re [73,75].

5.2.2. Flow in the 45${}^{°}$ Pipe Bend

In Figure 13 the time-averaged velocity magnitude for the 45${}^{°}$ pipe bends is shown for both Reynolds numbers. Due to the relatively sharp change of the flow direction an adverse pressure gradient develops in the upstream part. This leads to a recirculation region at the outer corner of the first bend which is found to be smaller than its corresponding 90${}^{°}$ counterpart. The flow attaches again at the outer wall in the middle segment between the two bends owing to the centrifugal force. At the inner corner of the first bend the flow detaches giving rise to a separation zone over the inner side of the middle segment. Hence, the observed acceleration of the fluid at the first bend is attributed to the reduction of the cross-section area of the flow.

A similar situation occurs at the second bend. An axial adverse pressure gradient leads to a vortical structure in the outer corner and the flow detaches again at the inner corner. Owing to the centrifugally induced pressure gradient, the high-velocity fluid is pushed to the outer side, whereas the low-velocity fluid is driven to the inner side downstream of the second bend. Hence, Dean vortices are formed as depicted by the streamlines of the secondary flow overlaying the contours of the downstream streamwise velocity $\langle v\rangle /U_b$ in Figure 14. The centers of these Dean vortices are closer to the symmetry plane and closer to the inner bend radius compared with the flow in the 90${}^{°}$ bend.

In Figure 15 the time-averaged resolved turbulent kinetic energy $\langle k\rangle /U_b^2$ is depicted for both Reynolds numbers. The largest magnitudes are again found in the shear layers but in contrary to the 90${}^{°}$ bend the turbulent kinetic energy is not very high in the recirculation regions itself. Furthermore, it has to be noted that the range of $\langle k\rangle /U_b^2$ depicted in Figure 15 comprises only half of the range shown in Figure 11. Obviously, redirecting the flow in two steps as done in the 45${}^{°}$ bend leads to lower velocity magnitudes and thus a weaker shear layer. As a direct consequence the production of turbulent kinetic energy is reduced.

The same applies to the time-averaged Reynolds shear stress $\langle u^{\prime }{v}^{\prime }\rangle /U_b^2$ shown in Figure 16. Again, the highest values are detected in the shear layers. Comparing the 45${}^{°}$ bend with the 90${}^{°}$ bend reveals that the deflection of the flow in two steps decreases the arising maxima of the Reynolds shear stress. In addition, in both geometries the size of the regions of high values of $\langle k\rangle /U_b^2$ and $\langle u^{\prime }{v}^{\prime }\rangle /U_b^2$ slightly decreases with increasing Reynolds number.

5.3. Breakage of Agglomerates

The breakage behavior of the agglomerates is evaluated based on an analysis of the data obtained within a dimensionless time interval of Δ$T^{*}=350$, which starts at the release of the first agglomerate.

Table 3. Summary of the main breakage results obtained in the pipe bends. The evaluation time interval is Δ$T^{*}=350$ starting at the release of the first agglomerate.

Bend	Case	Wall Impact	Rotation	${\mathit{N}}_{\mathbf{break}}^{\mathbf{tot}}/{\mathit{N}}_{\mathbf{agg}}^{\mathbf{rel}}$	${\mathit{N}}_{\mathbf{agg}}^{\mathbf{event}}/{\mathit{N}}_{\mathbf{agg}}^{\mathbf{rel}}$	PPF
		[%]	[%]			[%]
90${}^{°}$	${\mathrm{Re}}_L$	97.0	3.0	52.8	289.1	84.7
90${}^{°}$	${\mathrm{Re}}_H$	98.9	1.1	54.1	357.2	93.5
45${}^{°}$	${\mathrm{Re}}_L$	98.6	1.4	69.2	320.1	75.0
45${}^{°}$	${\mathrm{Re}}_H$	98.3	1.7	74.5	516.0	88.9

summarizes some of the most important results obtained for all considered cases. To begin with, the percentage contributions of the breakage mechanisms listed in the third and fourth column of Table 3 show that the breakage scenario is overwhelmingly predominated by wall impactions followed by a weak contribution of the rotary stress mechanism. Note that the contributions of the other two fluid-induced breakage mechanisms (i.e., drag and turbulence) are not included since they are found to be zero or negligibly small. The most obvious explanation for this observation is that the agglomerates collide with the wall before other physical effects such as strong shear or high turbulence intensities can have a significant effect on them. In other words, due to the geometrical setup most agglomerates simply do not reach the flow region which might lead to breakage by drag or turbulence. An exception is the rotary stress mechanism since agglomerates not breaking immediately at the first wall impact experience high rotation rates due to the collision with the wall which might induce this breakage mechanism.

The total number of breakage events per released agglomerate $N_{\mathrm{break}}^{\mathrm{tot}}/N_{\mathrm{agg}}^{\mathrm{rel}}$ listed in the fifth column of Table 3 suggests that breakage is enhanced with rising Re for both bend geometries. This can be clearly attributed to the higher impact velocities experienced by the agglomerates at Re${}_H$ compared to the case at Re${}_L$. Furthermore, more breakage events take place in the 45${}^{°}$ bend than in the corresponding 90${}^{°}$ bend at the same Re. A closer analysis reveals that this trend is related to the higher number of agglomeration events $N_{\mathrm{agg}}^{\mathrm{event}}/N_{\mathrm{agg}}^{\mathrm{rel}}$ (see sixth column in Table 3) occurring in the 45${}^{°}$ bends. In more detail, due to the introduction of a second bend the fragments of the agglomerates breaking in the middle segment tend to accumulate for a limited time in the regions near the corners of the second bend, where the fluid velocities are low. In these zones agglomeration is promoted due to the low-velocity interactions between the particles. The relatively smoother deflection of the flow in the 45${}^{°}$ bends leads to impact events characterized by small impact angles ${\Theta }_{\mathrm{imp}}$ as can be deduced from the PDF distributions of ${\Theta }_{\mathrm{imp}}$ depicted in Figure 17. As a consequence, the resulting fragments bounce off from the wall with small angles and advance nearly parallel to the axis of the middle segment to reach the regions near the corners of the second bend. The newly formed agglomerates are very likely to break again since they are already in close vicinity of the walls leading to higher breakage statistics in the 45${}^{°}$ bend compared to the 90${}^{°}$ case.

To provide further insight into the effect of the shape of the bend on the breakage behavior, Figure 18 depicts the positions of the breakage events along the axis of the pipes. Figure 18a confirms that in the ${90}^{°}$ pipe bend the majority of breakage events takes place near the kink at both Re. However, the distributions of the breakage events in the regions upstream and downstream of the bend are not similar at the two Re numbers. Obviously, more wall breakage events take place in the upstream region at Re${}_H$, which is attributed to the higher flow velocity and hence the more forceful wall-impact events encountered by the agglomerates even before reaching the bend. In contrast, the weaker flow conditions in the upstream region at Re${}_L$ shifts the distribution of the breakage events towards the kink and the downstream part of the bend.

The trends are overall similar in the ${45}^{°}$ pipe bend as illustrated in Figure 18b. In this bend configuration, the distributions show peaks near the kinks of the two bends. Again, slightly higher numbers of breakage events are observed in the upstream region at Re${}_H$, whereas in the Re${}_L$ case more agglomerates survive the critical regions near the kinks to break further downstream on their way towards the outlet.

An important aspect in the context of dry powder inhalers is the capability of the investigated pipe bends to deagglomerate powders. A common way to assess the deagglomeration performance is to analyze the fraction of particles arriving at the outlet with a certain size. For this purpose, the primary particle fraction (PPF) is taken into account, which is defined as the number of single primary particles reaching the outlet of the pipe divided by the total number of primary particles (i.e., including those forming agglomerates) reaching the outlet. The last column in Table 3 lists the PPF values obtained for the different configurations. As anticipated, the fraction of single primary particles reaching the outlet increases with the Reynolds number in both bends. In addition, for equal Re more single primary particles are detected at the outlet of the 90${}^{°}$ than for the 45${}^{°}$ bend. This circumstance is referred to the enhanced agglomeration in the 45${}^{°}$ bend as discussed above.

In the original numerical and experimental studies [5,10] in which the present geometrical configurations were introduced, the deagglomeration performance was characterized based on a quantity comparable to PPF. This quantity was denoted the fine particle fraction (FPF) and measures the mass fraction of the particles leaving the pipe with a diameter smaller than 5 $\mathsf{\mu }$m. It was concluded that the 45${}^{°}$ pipe bend offers a better deagglomeration performance (i.e., higher FPF) than the 90${}^{°}$ pipe bend. This finding does not coincide with the present results. The discrepancy is obviously related to the difference in the material properties (mannitol vs. silica particle agglomerates) and the size of the initial agglomerates considered. Specifically, the cohesion between the particles in the mentioned studies is higher (i.e., higher Hamaker constant) and the initial agglomerates consist of a remarkably larger number of primary particles (i.e., higher inertia). Hence, different motion and breakage patterns are expected making a direct comparison difficult. As mentioned in Section 4.3 reproducing the same agglomerates as applied in [5,10] would have required considerable efforts to generate completely new DEM training data (i.e., thousands of new DEM simulations) and to repeat the entire training procedure for the ANNs used in the breakage model. Thus, the present results specifically describes the breakage of relatively small silica particle agglomerates.

The primary particle fraction, however, focuses on the primary particles and does not offer knowledge about the size distribution of the surviving agglomerates. In order to analyze this quantity, the agglomerate size achieved at the outlet of the bends is evaluated by means of cumulative mass-weighted distributions $Q_3$. The results are shown in Figure 19. Similar trends are obtained for the 90${}^{°}$ and the 45${}^{°}$ pipe bends as visible in Figure 19a,b, respectively. Namely, small agglomerates comprising $N_p^{\mathrm{tot}}=5$ particles or less constitute about 50% of the mass population at the outlet. In general, higher Re numbers lead to smaller agglomerates which is consistent with the higher number of breakage events recorded in Table 3. In addition, large agglomerates consisting of more than $N_p^{\mathrm{tot}}=100$ particles are more present in the Re${}_L$ cases as can be inferred from the steeper gradients of the blue curves in that size interval. The small upsurge at the end of the cumulative distributions implies that in all cases some initial agglomerates (i.e., $N_p^{\mathrm{tot}}=1200$ particles) manage to arrive intact at the outlet of the pipe in all four cases.

6. Conclusions

An Euler-Lagrange method to study the turbulent particle-laden flow in pipe bend configurations relevant for dry powder inhalers is demonstrated. The method employs wall-resolved LES to predict the turbulent fluid flow. In addition, the disperse phase is tracked through the flow field while adopting the point-particle approximation to determine the fluid-particle interactions. Practical, event-based models are applied to predict collisions and agglomerations of particles as well as the breakage of agglomerates due to various physical mechanisms. Special emphasis is put on an artificial neural-network based approach for modeling the wall-impact breakage of nearly spherical agglomerates.

The scope of the study is to analyze the deagglomeration of silica agglomerates initially consisting of 1200 primary particles in four configurations defining two pipe bend geometries and two Reynolds numbers. In the adopted efficient methodology the structures of agglomerates are replaced by effective spheres, which allows to consider a large number of agglomerates in a single simulation in contrast to other detailed Euler-Lagrange approaches accounting for the full structure of agglomerates. The main outcomes are as follows:

• The flow in sharp pipe bends is mainly characterized by separation zones at the kinks of the bends and secondary flow structures known as Dean vortices. The presence of the separation regions accelerates the flow in the core of the bend. The abrupt deflection of the flow direction in the 90${}^{°}$ bend leads to large separation regions, high velocities, and strong turbulent fluctuations. The two-step deflection of the flow in the 45${}^{°}$ bend results in weaker flow conditions.
• In all four investigated configurations (two geometries and two Reynolds numbers) deagglomeration is mainly attributed to the wall-impact breakage followed by the breakage due to rotation. Other fluid-induced stress mechanisms such as drag and turbulence are found to be irrelevant under the considered, rather low Reynolds numbers. The same finding was reported by Tong et al. [10] who studied the breakage of agglomerates employing the same pipe geometries and flow Reynolds numbers, albeit different particle properties. The geometry of the bend pipes also contribute to this circumstance, since breakage by wall impactions is promoted by the sharp kink before agglomerates reach the regions of high flow shear.
• In any of the two investigated bend geometries, increasing Re enhances breakage and leads to finer particle size distributions at the outlet of the pipe.
• Considering the same Re number, more breakage takes place in the 45${}^{°}$ than in the 90${}^{°}$ pipe bend. This is attributed to the favorable agglomeration conditions in the 45${}^{°}$ case leading to more re-agglomerations and thus a higher breakage possibility. However, concentrating on the outlet of the pipe, a better deagglomeration performance is attained by the 90${}^{°}$ bend. This is perceptible from the higher primary particle fractions in Table 3 and the finer agglomerate size distributions at the outlet of the 90${}^{°}$ bends (see Figure 19).

In the future, the investigations might be extended to collisions between particles/ agglomerates, since this effect can also play a significant role for breakage.