Particle-resolved simulations of liquid fluidization of rigid and flexible fibers

Abstract

Particle-resolved, three-dimensional, time-dependent simulations of rigid and flexible cylinders fluidized by a liquid flow in fully periodic domains have been performed by means of the lattice-Boltzmann method supplemented with immersed boundaries. The solids volume fraction ranges from 0.10 to 0.48 and the length-over-diameter aspect ratio of the cylinders from 4 to 12. The bending stiffness of the cylinders is the third major input parameter. The resulting Reynolds numbers based on the average slip velocity of the cylinders and their equivalent diameter range from 6 to 70. It is shown that increasing the flexibility—that is, decreasing the bending stiffness—reduces the Reynolds number, an effect that is most pronounced for low solids volume fractions and long cylinders. As for rigid cylinders, the distribution of the orientation relative to the direction of gravity of the flexible cylinders is a pronounced function of the solids volume fraction and the aspect ratio. Flexibility tends to somewhat randomize the orientation distribution, which could explain the effect of flexibility on the slip velocity and thus the Reynolds number.

1 Introduction

Fiber suspensions are encountered in a number of process engineering applications. Two examples are pulp processing for paper production [1] and processing of bagasse—which is a fibrous residue from the production of sugar from sugarcanes—for, e.g., making biofuels [2]. Fiber suspensions are complex systems. In order to characterize them, one needs an extensive number of parameters and properties such as the fiber volume fraction, the size distribution of the fibers as well as their mechanical and surface properties. In principle, all these have an impact on the processability of the suspensions. Being able to anticipate the flow behavior of specific fiber suspensions has relevance for process design and optimization. In this paper, we make an attempt to quantify fiber suspension flow by means of numerical simulation. In the simulations, the objective is to account for the suspension properties in as much detail as possible. The main fiber properties we consider in the simulations are their shape that we quantify as the length over diameter aspect ratio and their bending stiffness. The continuous-phase fluid the fibers are suspended in is assumed to be a Newtonian liquid.

The applications we have in mind concern non-Brownian fibers with diameters typically in the range from 0.1 to 1 mm. On these length scales, solids and fluid inertia is appreciable with fiber-based Reynolds and Stokes numbers significantly larger than 1. There is significant experimental and computational literature on fibers with (sub)micrometer diameters suspended in liquids, mostly inspired by biological applications, for which inertia is not significant—see the excellent review by du Roure et al. [3] and references therein. Such systems allow for computations based on creeping flow assumptions, i.e., solving Stokes equations [4].

Collective behavior of suspensions involving inertial fibers has been studied experimentally in great detail with a view to fiber orientation and (anisotropic) velocity fluctuations [5]. Orientation distributions, hindered settling speeds and the formation of large-scale structures are pronounced functions of solids volume fraction and aspect ratio [6].

In recent papers [7,8,9], we have reported on particle-resolved simulations of inertial fiber suspensions where we solve for the flow dynamics through the lattice-Boltzmann (LB) method and impose no-slip conditions at the fiber surfaces by applying the immersed boundary method. An important feature of this approach is its ability to study the behavior of fibers in complex flow. As an example, in [8] suspensions of cylindrical particles were placed in a small container that was agitated by a revolving impeller. Flexibility of fibers has been introduced in [9]. There it was shown that fiber deformation is strongest near the impeller and leads to a reduced chance of fiber–impeller collisions with increase in fiber flexibility with consequences for the power required for stirring.

A limiting factor of particle-resolved simulation is the relatively large computational requirement per particle. Our current computational approach only allows for simulations with hundreds of rigid or flexible fibers. Even laboratory-scale systems usually have many orders of magnitude more fibers. A more generic way to employ particle-resolved simulation is the use of fully periodic boundary conditions, thereby mimicking the behavior of a suspension in an unbounded domain. In this paper, we consider fully periodic, three-dimensional domains containing Newtonian liquid through which fibers settle. The same systems can be interpreted as fluidized if we observe the fibers in a reference frame that moves with the average vertical fiber velocity. Such simulations provide information about the microstructure the fibers organize themselves in, the average relative velocity between liquid and solids, and the velocity fluctuation levels. Simulations of this kind with spherical particles have generated insights in hydrodynamic forces on spheres as a function of solids volume fraction and other dimensionless quantities such as Reynolds and Stokes numbers [10], constitutive equations for granular matter [11] and scalar dispersion in the interstitial liquid of dense suspensions [12]. In the previous work, we also have been studying hindered settling in periodic domains of particles that have the shape of red blood cells [13] and cylinders [7]. The latter reference is the starting point of the current paper.

The principal aim of this paper is to investigate, through simulation, to what extent flexibility of cylinders has an impact on the way they organize themselves while settling/being fluidized and how this influences their settling speed/slip velocity. We anticipate that flexibility effects are a strong function of the aspect ratio. Given that our previous work on cylinders [7] only went up to an aspect ratio (cylinder length over diameter) of 4, we first set the stage in the present paper by presenting results for rigid cylinders with an aspect ratio up to 12 that we subsequently compare to results of flexible cylinders. The three main degrees of freedom in our parameter space are the solids volume fraction, the aspect ratio and the flexibility of the cylinders.

In the next section, we define—in terms of dimensionless numbers—the flow systems and flow conditions. The computational methods have been described in detail previously [7, 9] and will be described briefly in this paper with an emphasis on the underlying assumptions. Section 4 begins with an account of tall rigid fibers. The orientation distributions and hindered settling speeds of rigid fibers form the baseline for discussing the results with flexible fibers. Main conclusions are reiterated in the last section of this paper.

2 Flow systems

Simulations have been performed in rectangular, fully periodic, three-dimensional flow domains of volume $V=nx\cdot ny\cdot nz$ with $nx\text{, } ny \text{ and } nz$ the domain size in x, y and z-direction, respectively. The domains contain Newtonian liquid (density $\rho$, kinematic viscosity $\nu )$ and a number n of solid cylinders with density $\rho _{p} >\rho$, diameter d and length $\ell$. The bending stiffness of the cylinders is $EI_{\mathrm{cs}}$ with E Young’s modulus and $I_{\mathrm{cs}}$ the moment of inertia of the cross-sectional area of the cylinder ($I_{cs} ={\pi d^{4}} / {64}$). Gravity points in the negative z-direction in the Cartesian coordinate system that we use: $\mathrm{\mathbf{g}}=-g\mathrm{\mathbf{e}}_{\mathrm{\mathbf{z}}}$.

The above input parameters are combined in a set of dimensionless groups: (1) aspect ratio $\ell / d$; (2) overall solids volume fraction $\left\langle \phi \right\rangle =\frac{n\pi \ell \mathrm{d}^{2}}{4V}$; (3) density ratio $\gamma ={\rho _{p} } / \rho$; (4) Archimedes number $\text{ Ar }=\left( {\gamma -1} \right) {gd^{3}} / {\nu ^{2}}$; and (5) a flexibility parameter $\chi =\frac{\left( {\rho _{p} -\rho } \right) g\mathrm{d}\ell ^{4}}{EI_{\mathrm{cs}} }$. The rationale for the definition of the dimensionless group $\chi$ comes from the endpoint deflection of a clamped beam under its own (net) weight $\delta =\frac{\left( {\rho _{p} -\rho } \right) g\tfrac{\pi }{4}\mathrm{d}^{2}\ell ^{4}}{8EI_{\mathrm{cs}} }$ so that $\chi \approx 10\delta / d$. As an example, a fiber with flexibility parameter $\chi =10$ deflects under its weight by $\delta \approx d$ on the one end if the other end is clamped.

As global output parameters we consider, in the first place, the average value of the Reynolds number based on the slip velocity between solid and liquid, and the equivalent fiber diameter $d_{e} =\root 3 \of {\tfrac{3}{2}\mathrm{d}^{2}\ell }$ is defined as $\text{ Re }=\frac{\overline{\left| {\left\langle {u_{z} } \right\rangle -\left\langle {u_{pz} } \right\rangle } \right| } \,\mathrm{d}_{e} }{\nu }$ with $\left\langle {u_{z} } \right\rangle$ the volume-averaged superficial fluid velocity in the z-direction and $\left\langle {u_{pz} } \right\rangle$ the average velocity of the fibers in z-direction. The overbar indicates time averaging. In the second place, we consider a measure for the root-mean-square fiber deflection scaled by the diameter: ${\Delta w} / d$ with $\left( {\Delta w} \right) ^{2}=\overline{\left\langle {\frac{1}{\ell }\int \nolimits _0^\ell {\left| {\mathrm{\mathbf{w}}} \right| ^{2}\mathrm{d}x_{1} } } \right\rangle }$; as before, $\langle \rangle$ denotes volume averaging. In the latter expression, $x_{1}$ is the coordinate along the cylinder centerline and $\mathrm{\mathbf{w}}\left( {x_{1} } \right)$ the deflection vector. In the third place, we study the distribution of orientation angles $\varphi$ between the fibers’ average centerline and the direction of gravity.

In the simulations discussed in this paper, the liquid properties, the fiber diameter and the fiber density are fixed so that the density ratio and Archimedes number have fixed values: $\gamma =2$ and $\hbox {Ar}=864$, respectively. As an example, a physical system with these characteristics has fibers of $d=1~\hbox {mm}$ and density $\rho _{p} =2000~\hbox {kg}/\hbox {m}^{\mathrm {3}}$ in a liquid with kinematic viscosity $\nu =3.4\cdot 10^{\mathrm {-6}}\,\hbox {m}^{\mathrm {2}}$/s and $\rho =1,000~\hbox {kg}/\hbox {m}^{\mathrm {3}}$ and experiences earth’s gravity ($g=9.8~\hbox {m}/\hbox {s}^{\mathrm {2}}$). The aspect ratio $\ell / d$ is varied in the range from 4 to 12, and the overall solids volume fraction $\left\langle \phi \right\rangle$ is in the range from 0.10 to 0.48. For rigid fibers $\chi =0$, while the flexible fibers have flexibility parameters up to $\chi =60$.

3 Simulation method and numerical settings

The numerical procedures have been built around the lattice-Boltzmann (LB) method that we use as a solver for the Navier–Stokes equations [14, 15]. The method provides us with a solution of the three-dimensional, time-dependent liquid flow on a uniform cubic grid with spacing $\Delta$ evolving in time with a time step $\Delta t$. The no-slip condition at solid surfaces (i.e., the surfaces of the deformed cylinders) is imposed by means of an immersed boundary method [16]. In this method, the surface of a deformed cylinder is represented by a set of closely spaced marker points. The typical spacing between marker points is $0.5\Delta$. We then locally apply forces on the fluid such that its velocity at each marker point matches the solid surface velocity at that point. Given that action equals minus reaction, the local forces on the fluid that maintain no slip on a fiber surface can be interpreted as a distribution of hydrodynamic forces acting on the fiber over its surface. This distribution of forces—in turn—determines the translation, rotation and deformation of the fiber.

In addition to hydrodynamic forces, the dynamics and deformation of fibers are determined by collisions with other fibers. For this, a contact model has been applied akin to discrete function representations (DFR) in the discrete element method [17]. The contact model first detects the proximity of marker points on different cylinders. Once the distance between two marker points lying on two different cylinders falls below a certain threshold value, a contact force is activated. The force has a linear elastic, repulsive contribution and a damping (or lubrication) contribution. The former depends on the distance between the marker points, and the latter depends on the relative velocity of the marker points. For determining the direction of the contact forces, each marker point is equipped with the outward normal vector at the (deformed) cylinder surface. This collision procedure has been detailed in previous papers [7,8,9]. Exactly, the same contact model coefficients as given in [7] have been used in the present paper.

Once we have the force distribution over each fiber as a result of fluid–solid and solid–solid interactions, we are ready to update the fiber’s linear and rotational velocity as well as its deformation. This is done in a two-step process. In the first step, we integrate the force distribution over the fiber surface to determine the total force and torque on the fiber which then enables updating the linear and rotational velocity according to Newton’s second law and the Euler equation, respectively. The Euler equation is applied in a reference frame attached to each cylinder. The rotational transformations involved in transferring information between the overall, inertial Cartesian coordinate system and the cylinder reference frames are facilitated by quaternions [18]. For this, the orientation of each cylinder is monitored by a quaternion which is updated each time step according to its angular velocity and the exact solution due to [19].

In the second step, the distribution of forces over the cylinder is first corrected for inertial effects [9]. The corrected force distribution is then used to determine the bending moment along the centerline of each cylinder and subsequently to determine the deflection along the centerline. This all is based on linear elastic theory under the assumption of small deflections and quasi-static deformation [20]. Determining the bending moment distribution and deflection along the centerline involves solving boundary value problems of second-order ordinary differential equations which we do by means of a finite difference method. The resulting systems of linear algebraic equations are solved directly through Gauss elimination (given that the number of nodes along a cylinder centerline—and thus the size of the linear systems—is 60 at most). Finally, the deflection along the centerline is used to bend the cylinder which implies a displacement of the marker points used for the immersed boundaries and collision detection. It therefore is the deformed (bended) fiber shape that interacts with the interstitial liquid and collides with other cylinders.

The small-deflections approximation, as used when determining the deformation, justifies two other approximations. In the first place, for simplifying the algorithm, when solving the Euler equation the moment of inertia tensor of the undeformed cylinder is used. In the second place, when applying the immersed boundary method, we do not include the velocity associated with the change in deformation over time when locally forcing the fluid velocity to the cylinder surface velocity at the marker points.

In previous papers [7, 9], a number of verification tests on the simulation procedure as described above have been performed. We refer to these papers to justify the numerical settings of the simulations described in the present paper. The spatial resolution has been set such that $d=16\Delta$. Comparing with results obtained with $d=24\Delta$ showed—for rigid cylinders—good agreement in terms of the average settling velocity Reynolds numbers as well as Reynolds numbers based on velocity fluctuations [7]. The simulation domains also need to be sufficiently large in order to avoid unphysical self-interaction of fibers over the periodic boundaries. As before [7], the domain size in the direction of gravity is two times the domain size in the lateral directions: $nz=2nx=2ny$. In [7], it was shown that as long as $nx\ge 12d$, the domain size has a minor effect on average slip velocity Reynolds numbers. However, it has a more pronounced effect on solids velocity fluctuation levels. In the present paper, the cylinders are taller than previously. We make sure that the height of the domain is at least three times the cylinder length: $nz\ge 3\ell$ For the longest fibers (with $\ell / d=12$), we have domains with $nz=576\Delta$ and $nx=ny=288\Delta$. The above verifications dealt with rigid fibers. Verification tests for flexible settling cylinders [9] also show that $d=16\Delta$ is an appropriate resolution.

The number of nodes used to represent the bending moments and deflection distributions over the length of a fiber has been set such that there are five nodes over a distance equal to one cylinder diameter. This implies that the fibers with aspect ratio $\ell / d=12$ are divided into 60 segments. Sensitivity of the deflection levels with regard to the number of segments has been checked in [9] for a single, settling cylinder. It showed negligible differences between results obtained with 5 and 7.5 segments per distance d.

The time step in the simulations is such that ${d^{2}} / \nu =6400\Delta t$. This is a small time step, instigated by maintaining a low Mach number so that the LB method simulates effectively an incompressible flow. Over one time step, a cylinder moves over not more than 0.005 times its diameter. A typical simulation runs for $2\cdot 10^{\mathrm {5}}$ to $4\cdot 10^{\mathrm {5}}$ time steps.

The simulations run in a sequential manner on a standard linux compute server. Given the grid size (up to $288\times 288\times 576$) and the number of time steps, this results in long run times (up to 4.5 months). Parallelization based on domain composition clearly is an option but would require quite intricate additional coding (think of dealing with one fiber that is simultaneously present in a number of subdomains) that we leave for the future.

4 Results

4.1 Rigid fibers

The structure of the fiber suspension is a pronounced function of the solids volume fraction, as can be seen in Fig. 1. In this figure, we see single realizations of rigid fibers with $\ell / d=12$ in the periodic domains. The snapshots are taken after the system has reached a dynamic steady state. In the densest case with $\left\langle \phi \right\rangle =0.48$, the fibers have a very strong preference for an orientation aligned with the direction of gravity, i.e., along the z-axis, a phenomenon also known from experimental work [5]. The leaner the suspension, the more random the orientations get. This effect is quantified in Fig. 2 where we show distribution functions of the angle $\varphi$ between the fiber centerline and the direction of gravity (the z-axis). The data for these distributions have been collected from all fibers in the simulation over at least 1000 instantaneous realizations separated in time by $0.02{d^{2}} / \nu$. The case with $\left\langle \phi \right\rangle =0.48$ and $\ell / d=12$ is an exceptional case with a high peak in the angle distribution function at $\varphi \approx 0.042\pi$. For this solids volume fraction, also the shorter cylinders (with $\ell / d=8$ and 4) settle mostly vertically, although with a wider angle distribution than the $\ell / d=12$ cylinders. The orientations are more random for the lower solids volume fractions where, for $\left\langle \phi \right\rangle =0.10$ and $\ell / d=4$, the angle distribution function approaches the $\sin \varphi$ curve which indicates that the cylinders do not have a preferential orientation. For the longer cylinders at $\left\langle \phi \right\rangle =0.10$, (near) horizontal settling, i.e., $\varphi \approx \pi / 2$, is still underrepresented as compared to a random isotropic distribution.

Fig. 1

Instantaneous realizations of liquid-fluidized rigid cylinders with $\ell / d=12$ and overall solids volume fraction $\left\langle \phi \right\rangle =0.48,~0.29,~0.10$ (left to right) after dynamic steady state was reached. Fully periodic conditions with a cylinder displayed if its center falls inside the $nx\times ny\times nz$ flow domain (indicated by the rectangular box). Far left: definition of coordinate system

Fig. 2

Distribution of the angles $\varphi$ between the centerline of the cylinders and the vertical. Effect of overall solids volume fraction $\left\langle \phi \right\rangle$ (left to right) and aspect ratio $\ell / d$. The solid curve in each panel is $\sin \varphi$ which represents an isotropic random distribution of orientations

Fig. 3

Time series of the volume-averaged Reynolds number $\text{ Re }={\left| {\left\langle {u_{z} } \right\rangle -\left\langle {u_{pz} } \right\rangle } \right| d_{e} } / \nu$. Top: $\left\langle \phi \right\rangle =0.48$; bottom $\left\langle \phi \right\rangle =0.20$. Time averaging starts at ${t\nu } / {d_{e}^{2} }=1.5$. Rigid cylinders with aspect ratios as indicated

In a previous paper [7], we have reported on settling speeds of cylinders with aspect ratios in the range of ${0.5\le \ell } / d\le 4$. An important observation was that hindered settling follows the dependence with the solids volume fraction according to Richardson and Zaki [21]. In terms of the Reynolds number (defined above): $\mathrm{Re}=\mathrm{Re}_{\infty } \left( {1-\left\langle \phi \right\rangle } \right) ^{N}$. Since the Reynolds number is based on the equivalent particle diameter $d_{e}$, it is meaningful to compare with results for spherical particles: $N=4.45\mathrm{Re}_{\infty }^{-0.1} \quad \text{ for } 1<\mathrm{Re}_{\infty } <500$ [21]. The exponent N for cylinders appeared to be consistently lower than for spheres [7]. In the current paper, we report on cylinders with ${4\le \ell } / d\le 12$.

Before presenting results on settling speed Reynolds numbers as a function of $\ell / d$ and $\left\langle \phi \right\rangle$, we first discuss how this Reynolds number has actually been determined. Figure 3 shows time series of the volume-averaged Reynolds number for a number of conditions. Usually, at time equal to zero, a simulation starts from a zero velocity state with the cylinders placed horizontally in a regular array such that they are not overlapping and not touching and with some randomness built in (a few randomly picked holes in the grid and/or the fibers having small random orientation angles relative to their neighbors). We then switch on net gravity on the cylinders and an opposing body force on the fluid so that the domain is overall force balanced [7, 11]. The systems then go through a clearly transient stage, the length and form of which depend on the specific conditions as well as the initial configuration of particles. We carefully monitor this stage, mainly in terms of the volume-averaged settling Reynolds number, but also by looking at snapshots of fiber configurations (as in Fig. 1) and by assessing whether the distribution of angles (as in Fig. 2) becomes steady.

After approximately one viscous timescale ${d_{e}^{2} } / \nu$, a dynamic steady state sets in. We then start time-averaging the volume-averaged Reynolds number, the values of which are reported in Fig. 4. It is clear from Fig. 3 that $\left\langle \phi \right\rangle =0.20$ and $\left\langle \phi \right\rangle =0.48$ show significantly different fluctuation levels. There are two reasons for this: in the first place because volume averaging involves averaging over more particles when $\left\langle \phi \right\rangle$ increases, and in the second place because there is not much room for the particles to fluctuate in a very dense suspension. Solids velocity fluctuations have been discussed in the previous paper [7] and will not be discussed here any further. The levels of fluctuation as observed in Fig. 3, however, have an impact on the accuracy of the time-averaged value of the Reynolds number [22]. A coarse estimate for this can be expressed as $\varepsilon _{\mathrm{Re}} =\sigma _{\mathrm{Re}} \sqrt{{2T_{f} } / T}$ with $\sigma _{\mathrm{Re}}$ the standard deviation of the fluctuating Re signal, $\varepsilon _{Re}$ the uncertainty in the Reynolds number time-averaged over a time window T and $T_{f}$ an integral flow timescale. For the latter, we take ${T_{f} =d_{e} } / {u_{s} }$, i.e., the equivalent fiber diameter over the slip velocity. If we nondimensionalize $T_{f}$ with the viscous time $T_{\nu } ={d_{e}^{2} } / \nu$, we get ${{T_{f} } / {T_{\nu } }=1} / {\mathrm{Re}}$. The averaging times are of the order of $T / {T_{\nu } }\approx 5$ (see Fig. 3). Then, the uncertainty in the time-averaged Reynolds number becomes $\varepsilon _{\mathrm{Re}} =\sigma _{Re} \sqrt{2 / {\left( {5\mathrm{Re}} \right) }}$. With reference to the cases displayed in Fig. 2, this implies a relative uncertainty in the average Reynolds numbers (${\varepsilon _{\mathrm{Re}} } / {\mathrm{Re}}$) of 1 to 2%.

Fig. 4

Hindered settling Reynolds number Re (volume as well as time averaged) as a function of porosity $1-\left\langle \phi \right\rangle$ in a double-logarithmic (“Richardson and Zaki”) representation. Rigid cylinders with $\ell / d=4,6,8,10,12$. The dashed lines are the best fits through the three data points on the right for each aspect ratio

Fig. 5

Impressions of liquid-fluidized flexible cylinders with $\ell / d=12$ and overall solids volume fraction $\left\langle \phi \right\rangle =0.10,~0.29,~0.48$ (left to right) and flexibility parameter $\chi =2.3$ (top) and 7.8 (bottom). Fully periodic conditions with a cylinder displayed if its center falls inside the $nx\times ny\times nz$ flow domain (indicated by the rectangular box). Far left: definition of coordinate system

Figure 4 shows average (over volume & time) Reynolds numbers as a function of solids volume fraction for five different aspect ratios. The results for $\ell / d=4$ have been presented in a previous paper [7]. The way the data are presented in Fig. 4, i.e., $1-\left\langle \phi \right\rangle$ versus Re on a double-logarithmic scale, allows for quickly assessing the validity of Richardson and Zaki [21] scaling in an analogy with spherical particles. We are not aware of correlations describing hindered settling of cylinders. The overall increase in Re with $\ell / d$ is due to an increase in particle weight and thus slip velocity as well as an increase in $d_{e}$.

According to Fig. 4, the Richardson and Zaki correlation $\mathrm{Re}=\mathrm{Re}_{\infty } \left( {1-\left\langle \phi \right\rangle } \right) ^{N}$ ceases to fit the data points for $\ell / d>6$: clearly, the data points in Fig. 4 for $\ell / d\ge 8$ do not follow a straight line with the deviations from the straight trend lines increasingly strong for small $1-\left\langle \phi \right\rangle$, i.e., for the dense suspensions. We attribute this to the preferential alignment of the fibers with the direction of gravity, which is particularly significant for the longer cylinders and high solids volume fractions.

4.2 Flexible fibers

Impressions of flexible fiber fluidization are shown in Fig. 5. As for Fig. 1, these are instantaneous realizations at a moment in time after a dynamically steady state has set in. The bending deformations of the fibers are clearly visible, with deformations increasing with increasing flexibility parameter $\chi$. A striking difference when comparing Fig. 1 (rigid fibers) with Fig. 5 (flexible fibers) is the higher levels of randomness in the orientations of flexible fibers at $\left\langle \phi \right\rangle =0.48$.

Flexibility has an impact on the average Reynolds number, see Fig. 6 where we show results of all simulations conducted in this study. Mostly, the Reynolds numbers for flexible fibers are lower than their rigid counterparts. The extent to which the Reynolds number changes with flexibility depends on the solids volume fraction $\left\langle \phi \right\rangle$, the flexibility parameter $\chi$, as well as the aspect ratio $\ell / d$. This means that although the length of the fibers is contained in $\chi$, its effect cannot be captured solely by $\chi$. In order to interpret the results for Re as a function of $\left( {\left\langle \phi \right\rangle ,\chi ,\ell / d} \right)$, we show in Figs. 7 and 8 data on the levels of deformation of the fibers and orientation of the fibers, respectively.

Fig. 6

Hindered settling Reynolds numbers Re of flexible cylinders as a function of their flexibility $\chi$ for solids volume fractions $\left\langle \phi \right\rangle =0.48,~0.40,~0.29,~0.20$ and 0.10 as well as aspect ratios $\ell / d=12,~10,~8,~6$ and 4

Fig. 7

Average cylinder deformation ${\Delta w} / d$ (defined in the text) as a function of flexibility $\chi$ for solids volume fractions $\left\langle \phi \right\rangle =0.48,~0.40,~0.29,~0.20$ and 0.10 as well as aspect ratios $\ell / d=12,~10,~8,~6$ and 4

Fig. 8

Distribution of the angles $\varphi$ between the centerline of the flexible cylinders and the vertical. Effect of overall solids volume fraction $\left\langle \phi \right\rangle$ (left to right), the aspect ratio $\ell / d$ (top to bottom) as well as cylinder flexibility $\chi$. The solid curve in each panel is $\sin \varphi$ which represents an isotropic random distribution of orientations

For the shorter fibers with $\ell / d=4$ and 6, there is virtually no effect of flexibility on Re. This is despite the fact that flexibility does change the orientation distribution of short fibers, as shown (for $\ell / d=4)$ in Fig. 8. In the dense suspension with $\left\langle \phi \right\rangle =0.48$, flexibility weakens the preferential vertical orientation, whereas in the lean suspension with $\left\langle \phi \right\rangle =0.10$, flexibility favors the fibers to get horizontal.

Flexible fibers with $\ell / d\ge 8$ generally slow down as a results of flexibility with this effect being most pronounced for the two ends of the solids volume fraction spectrum (i.e. $\left\langle \phi \right\rangle =0.10$ and 0.48). For $\left\langle \phi \right\rangle =0.48$, Reynolds numbers reduce by as much as 20% for $\ell / d=12$ and 10 with the reduction leveling off with increasing $\chi$ quickly. We note that also fiber deformation tends to level off with $\chi$, see Fig. 7. Long fibers in dense suspension ($\left\langle \phi \right\rangle =0.48)$ tend to randomize their orientation as compared to rigid fibers. This was already observed in the instantaneous realizations (Figs. 1 and 5) and has been quantified in the angle distributions: compare the lower left panel of Fig. 8 to the left panel of Fig. 2 (for rigid fibers). As we hypothesized when discussing hindered settling of rigid fibers, preferential vertical orientation increases the Reynolds number. Therefore, one might expect that a wider distribution of orientation angles as a result of flexibility reduces the Reynolds number.

For lean suspensions: The reduction in Re with flexibility at $\left\langle \phi \right\rangle =0.10$ for the longer fibers is of the order of 15%. Also here we believe most of this is due to a randomization of the orientation, see the lower right panel of Fig. 8. There the most flexible fibers closely follow the $\sin \varphi$ function which indicates an isotropic distribution. The rigid fibers with $\ell / d=12$ at $\left\langle \phi \right\rangle =0.10$ have a clear preference for smaller angles $\varphi$ and thus vertically oriented settling.

For the intermediate solids volume fractions, $0.20\le \left\langle \phi \right\rangle \le 0.40$, the influence of flexibility on Re is modest (variations within 10%) and on some occasions (e.g., $\left\langle \phi \right\rangle =0.29)$ show a slight increase in Re with increasing flexibility parameter $\chi$. It is noted that the orientation distributions for suspensions with $\left\langle \phi \right\rangle =0.29$ are not much affected by flexibility (see Fig. 8).

In all cases, average fiber deformation, ${\Delta w} / d$ as defined above, increases monotonically with $\chi$, as shown in Fig. 7. However, deformation tends to level off at the high end of $\chi$. This seems logical for the dense suspensions where the fibers have limited space to maneuver and them bending strongly is hindered by fibers in their direct proximity. It is less obvious why and to what extent the leveling off occurs for lower $\left\langle \phi \right\rangle$. Capturing this would need more data points at higher $\chi$ values, specifically for $\ell / d=12$ at $\left\langle \phi \right\rangle =0.10$ and 0.20. This we leave for future work. As an overall trend, fiber deformation increases with a decreasing solids volume fraction. This is due to a reduction in hindrance by other fibers when the solids volume fraction decreases.

5 Summary and conclusions

In this paper, we have applied a numerical methodology developed in earlier papers [7, 9] to fluidization/sedimentation of rigid and flexible fibers with length-over-diameter aspect ratios $\ell / d$ up to 12 in a Newtonian liquid. The simulation explicitly accounts for the (deformed) shape of the fibers by imposing no-slip at their surface. Hydrodynamic and contact forces move (translate and rotate) and deform the fibers.

Fiber deformation gets stronger with decreasing solids volume fraction when the fibers have more free space around them. In dense suspensions, long fibers very strongly align with the direction of gravity, which makes them settle relatively fast, departing significantly from a Richardson- and Zaki-type scaling of hindered settling speeds. Flexibility of the fibers to some extent opposes the preferential orientation and—as a result—reduces the settling speed. In more dilute suspensions, settling speeds also get reduced as a result of flexibility. This likely is related to the fact that flexibility promotes horizontal fiber orientations.

For fibers with $\ell / d$ smaller than 8, flexibility only has minor impact on the Reynolds numbers based on settling velocity.

In the future, we very much would like to compare the results obtained in this paper with experimental data. Very relevant experiments of a relatively simple form would be measuring the sedimentation velocity of fibers in suspension as a function of the degrees of freedom explored in this paper: solids volume fraction, aspect ratio and flexibility. A precise characterization of the fibrous material—specifically determining its bending stiffness—would be essential for such experiments.