Next Article in Journal
The Fast Discrete Interaction Approximation Concept
Next Article in Special Issue
Electroviscous Effects in Stationary Solid Phase Suspensions
Previous Article in Journal
Heat-Dissipation Performance of Nanocomposite Phase-Change Materials in a Twin-Heat-Source System
Previous Article in Special Issue
New Generalized Viscosity Model for Non-Colloidal Suspensions and Emulsions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

CFD Analysis of Turbulent Fibre Suspension Flow

1
Department of Civil, Environmental and Natural Resources of Engineering, Luleå University of Technology, 971 87 Luleå, Sweden
2
VOLVO AB, Gropegårdsgatan 2, 417 15 Göteborg, Sweden
3
ÅF Industry AB, Grafiska vägen 2, 401 51 Göteborg, Sweden
*
Author to whom correspondence should be addressed.
Fluids 2020, 5(4), 175; https://doi.org/10.3390/fluids5040175
Submission received: 16 August 2020 / Revised: 25 September 2020 / Accepted: 1 October 2020 / Published: 8 October 2020
(This article belongs to the Special Issue Fluid Mechanics of Suspensions and Emulsions)

Abstract

:
A new model for turbulent fibre suspension flow is proposed by introducing a model for the fibre orientation distribution function (ODF). The coupling between suspended fibres and the fluid momentum is then introduced through the second and fourth order fibre orientation tensors, respectively. From the modelled ODF, a method to construct explicit expressions for the components of the orientation tensors as functions of the flow field is derived. The implementation of the method provides a fibre model that includes the anisotropic detail of the stresses introduced due to presence of the fibres, while being significantly cheaper than solving the transport of the ODF and computing the orientation tensors from numerical integration in each iteration. The model was validated and trimmed using experimental data from flow over a backwards facing step. The model was then further validated with experimental data from a turbulent fibre suspension channel flow. Simulations were also carried out using a Bingham viscoplastic fluid model for comparison. The ODF model and the Bingham model performed reasonably well for the turbulent flow areas, and the latter model showed to be slightly better given the parameter settings tested in the present study. The ODF model may have good potential, but more rigorous study is needed to fully evaluate the model.

Graphical Abstract

1. Introduction

The Intergovernmental Panel on Climate Change (IPCC) has concluded that the emission of greenhouse gases does have an effect on the climate and that the levels are now the highest in history [1]. In Sweden, the pulp and paper industry consumes 52% of the power in the whole industrial sector [2]. Achieving increased energy efficiency in the industry sector is therefore highly relevant in contributing to the development towards a climate neutral society, as the use and supply of energy causes approximately 60% of the emission of greenhouse gases [3].
A crucial step to design and optimise sustainable and energy-efficient methods for the different processes used within the pulp and paper industry is understanding the flow of pulp suspensions.
Fibre modelling is a complex physical problem and fibre suspension rheology depends on various parameters such as fibre concentration and fibre type, but also the flow regime itself [4]. For instance, in turbulent flows, suspensions are referred to as “fluidised” since the turbulent fluctuations affects the flow in such a way that the suspension behaves much like a Newtonian fluid [5]. On the other hand, in flows of decaying turbulence, for instance, fibres form local concentrations of fibres sticking together, called flocs, which have an impact on the suspension rheology [6]. In papermaking, the decaying turbulence flow is the most common type of flocculating flow [4].
Because of the complex phenomena governing fibre suspension rheology, existing modelling approaches also span a wide range of complexity and level of detail in description. Fibres may be tracked individually, modelled as linked rigid segments in a flow field [7,8]. The suspension may be as a non-Newtonian fluid exhibiting a certain yield stress [9]. The suspension may also be described by the coupling between fluid momentum and the orientation distribution of fibres [10,11,12,13,14,15,16,17,18] The two latter approaches are macroscopic models that describe the fibre suspension flow as the collective result of fibres effecting the rheology, without a detailed description of individual fibres. In the present work, a macroscopic model for turbulent fibre suspension flow suitable for engineering applications was desired. Especially, softwood fibres used in papermaking processes were of interest. The orientation distribution approach was used to propose a new model and derive a computationally efficient implementation of the same.

2. Theory

A constitutive equation for the deviatoric stress tensor in the fibre suspension of rigid fibres was used by Shaqfeh and Fredrickson (1990) [18] as below:
τ i j = 2 μ S i j + μ f ( a i j k l 1 3 δ i j a k l ) S k l
where μ is the viscosity of the suspending fluid, S i j the local strain rate, μ f the additional viscosity due to the presence of fibres and δ i j is the unit tensor. The fibre orientation tensors a i j and a i j k l are the second and fourth order statistical moments of fibre orientation, respectively, and are defined as
a i j = p i p j Ψ = p i p j Ψ ( p ) d p  
a i j k l = p i p j p k p l Ψ = p i p j p k p k Ψ ( p ) d p
where Ψ is the orientation distribution function (ODF) and p is the fibre orientation vector. The additional viscosity due to fibre presence μ f depends on fibre properties and fibre concentration and may be computed for a semi-dilute suspension as Shaqfeh and Fredrickson (1990) [18]
μ f = μ 8 π n l 3 3 [ ln ( 1 C V ) + ln ( ln ( 1 C V ) ) + 1.4389 ]
where n is the number density of fibres, l is the half-length of fibres and C V the volume concentration of fibres. Semi-dilute suspensions are here referred to as those satisfying n l 3 1 and n l 2 d < 1 . By applying the constitutive relation in Equation (4) to the fluid momentum equations, a modified version of the incompressible Navier–Stokes equation is obtained as Yang (2013) [10]
u i t + u j u i x j = 1 ρ p x i + μ ρ 2 u i x j x j + μ f ρ x j [ a i j k l S k l 1 3 ( δ i j a k l ) S k l ]
where u i is the velocity, ρ the density of the suspending fluid and p is pressure. It may be remarked that the body force term ρ f i was left out by Yang (2013) [10] but may be included by the addition on the right-hand side of Equation (5). To close the system of equations for the fibre suspension flow, the probability density for the fibre orientation distribution, the ODF Ψ is needed. The evolution of Ψ in time and space can be described by a probability density transport equation, called the Fokker–Planck equation, reading Zhang, (2014) [11]:
t Ψ ( p , x , t ) + p F p + x F x
where p is the fibre orientation vector, x the fibre position and p and x are the gradient operators in orientation and translation space, respectively, and F p and F x are the flux densities in the respective spaces. The flux densities require complex modelling in order to close Equation (6).
Also, a constitutive equation for the deviatory stress tensor in the fibre suspension of rigid fibres was used by (Shaqfeh and Fredrickson, 1990) [18] as below:
τ ij = 2 μ S ij + μ f ( a ijkl 1 3 δ ij a kl ) S kl
where μ is the viscosity of the suspending fluid, S ij the local strain rate, μ f the additional viscosity due to the presence of fibres and δ ij is the unit tensor. The fibre orientation tensors a ij and a ijkl are the second and fourth order statistical moments of fibre orientation, respectively, and are defined as
a ij = p i p j Ψ = p i p j Ψ ( p ) dp   ,
a ijkl = p i p j p k p l Ψ = p i p j p k p k Ψ ( p ) dp
where Ψ is the orientation distribution function (ODF) and p is the fibre orientation vector.

3. Method

Solving the fluid flow from the modified Navier–Stokes Equation (5) and the Fokker–Planck Equation (6) requires discretization and iteration on the two equations. It also requires the computation of the orientation tensors a i j and a i j k l , which have 21 unique components in total, in order to obtain the coupling of the equations. This needs to be executed through numerical integration in each computational cell and in each iteration. Shortly, solving the suspension flow comes with high computational cost. Solving the Fokker–Planck is also a complex problem itself. In the present work, a new model for the ODF is instead proposed. In addition, an implementation where explicit expressions for the orientation tensors as functions of the fluid flow are constructed prior to simulations is derived from the model with the use of Fourier series. The method substantially reduces the computational effort needed to solve the suspension flow, whilst keeping a lot of the detail in the anisotropic description of fibre stresses provided by the orientation distribution approach.

3.1. Modelling the ODF

In several numerical studies of the orientation distributions of fibres in turbulent suspension flows, slender fibres have been observed as tending to align with the fluid flow streamlines [12,14,15,16,17]. Some of these results have also been confirmed with experimental data [16,17]. Lin [16] also found that turbulent fluctuations have a randomizing effect on the alignment and that the alignment is faster the higher the aspect ratio of the fibres is. Therefore, the model developed in the present work was based on three main assumptions, namely;
  • Fibres are of high aspect ratio causing them to align according to the fluid flow direction with negligible delay;
  • The deviation of the alignment from the fluid flow direction increases with a magnitude of turbulent fluctuations; and
  • Fibres may be treated as rigid so that the Equations (4) and (5) may be used to described the stresses acted on the fluid by the suspended fibres.
The validity of the last assumption for softwood fibres may be discussed, as it leaves out possible degrees of freedom such as fibre bending, twisting and hooking, which causes contact forces affecting fibre suspension rheology [4]. It was, though, assumed that the turbulence would cause these effects to be off small importance compared to the interaction between fibres and turbulent fluctuations, and that the orientation distribution approach itself was sufficient to describe the suspension rheology.
A statistical model is proposed were the ODF is explicitly determined by the local velocity direction and the turbulent intensity in the fluid flow. Using the model, the orientation tensors could in turn be computed for all possible cases and explicit expressions can be constructed for all orientation tensor components a i j and a i j k l , as functions of the flow field and prior to simulations. The explicit expressions may then be used to compute the fibre stress term in Equation (5) directly from the flow field. In this way, no additional differential equations need to be solved, and the numerical integration for the orientation tensors in each computational cell and iteration is also avoided. A large amount of computational effort is thus saved.
Fibre orientation is described using the angular components of a spherical coordinate system, which can be seen in Figure 1, and the fibre orientation vector then reads [12]:
p = [ cos ϕ sin θ sin ϕ sin θ cos θ ]
The fibre orientation angles ϕ and θ are assumed to be independently normally distributed with means corresponding to the local velocity direction and with a common variance σ 2 as
ϕ N ( ϕ ¯ , σ 2 )
θ N ( θ ¯ , σ 2 )
where ϕ ¯ and θ ¯ are the angles corresponding to the local direction of fluid velocity. The probability density functions (pdf) of the distributions of ϕ and θ , respectively, then read (Rice, 2007) [19]:
Ψ ϕ ( ϕ ) = 1 σ 2 π e ( ϕ ϕ ¯ ) 2 2 σ 2
Ψ θ ( θ ) = 1 σ 2 π e ( θ θ ¯ ) 2 2 σ 2
The modelled ODF, being the joint distribution of ϕ and θ , is then constructed as
Ψ ( ϕ , θ ) = Ψ ϕ ( ϕ ) Ψ θ ( θ )
ODF variance is connected to turbulence through the turbulent intensity of the flow. The randomizing effect on the fibre distribution caused by turbulent velocity fluctuations is modelled by letting the ODF standard deviation σ be a linear function of the turbulent intensity I as
σ = m I = m U B ( 2 3 k )
where m is the slope, U B is a bulk velocity of the flow and k is the turbulent kinetic energy. A minimum value of 0.05 rad is defined considering perfect alignment being non-realistic. Since the expressions for the orientation tensors are to be constructed prior to simulations for a range of values of σ , a limited interval is needed to be defined and thus a maximum for σ needs to be defined. A maximum of 0.5 rad is defined being considered to cover the interval of realistic values observed in the literature [12,14,15,16,17]. The construction of explicit expressions for the orientation tensors is described below.

3.2. Explicit a Priori Orientation Tensors

For a given ODF variance σ 2 and flow field direction ( ϕ ¯ ,   θ ¯ ) the modelled ODF is known and all orientation tensor components can be computed, as an example for a i j , :
a i j = p i p j Ψ ( ϕ , θ ) d ϕ d θ = ( f ( ϕ ) Ψ ϕ ( ϕ ) d ϕ ) ( g ( θ ) Ψ θ ( θ ) d θ ) F ( ϕ ¯ ) G ( θ ¯ ) ,
where the orientation vector components p i consist of a product of functions of ϕ and θ , respectively, so that the tensor p i p j can be written on the form f ( ϕ ) g ( θ ) . The double integral can therefore be separated into two integrals and are then functions of the flow field angles ϕ ¯ and θ ¯ and are therefore denoted F ( ϕ ¯ ) and G ( θ ¯ ) , respectively. The values for F ( ϕ ¯ ) and G ( θ ¯ ) are computed through numerical integration for a series of discrete values for ϕ ¯ and θ ¯ , respectively, covering all possible fibre directions. The angles then cover the interval [ π 2 , π 2 ] , which may be interpreted from Figure 1. Then, Fourier series expansions are applied and fitted to F ( ϕ ¯ ) and G ( θ ¯ ) , respectively, on the form Nordling and Österman, (2006) [20]:
F ( ϕ ¯ ) F ¯ ( ϕ ¯ ) a 0 2 + n = 1 N ( a n cos n π ϕ ¯ L + b n sin n π ϕ ¯ L )  
G ( θ ¯ ) G ¯ ( θ ¯ ) c 0 2 + n = 1 M ( c n cos n π θ ¯ L + d n sin n π θ ¯ L )  
An expression for the tensor component is then constructed as
a i j ( ϕ ¯ , θ ¯ ) = F ¯ ( ϕ ¯ ) G ¯ ( θ ¯ )
To allow for the ODF standard deviation to vary throughout the flow field, the procedure is repeated for a discrete series of values for σ and the local value of a i j can then be computed by linear interpolation between the constructed expressions. The full procedure is applied to all orientation tensor components a i j and a i j k l . As a remark, only a number of terms up to M , N = 3 is needed to obtain a good fit for the Fourier series. An example can be seen for the component a 1123 in Figure 2. A similarly good fit was obtained for all tensor components and for the range of values for the ODF standard deviation σ used.
Using the model derived with the Fourier series implementation, the fibre stress term in Equation (5) can be computed directly from the flow field. The ODF Ψ does not need to be solved and the 21 unique components a i j and a i j k l do not need to be computed by numerical integration in each cell and iteration, as they are determined explicitly from the fluid flow. −

3.3. Model Implementation

In the present work, the model was intended for simulation with the commercial Computational Fluid Dynamics CFD software Ansys Fluent [21]. The model described above was therefore implemented with the use of so-called user-defined functions (UDF), which are user-created subroutines to the built-in solver in Fluent written in the C programming language. In total, three UDFs were created and used to implement the ODF model for the Fluent framework. So-called user-defined scalar (UDS) variables, which are additional transport equations defined by the user, were used for the storage of variables and since Fluent automatically computes the gradient of all UDS variables. The solution of the UDS transport equation, however, was turned off. The first UDF computes the local values of orientation tensor components a i j   and a i j k l from the explicit expressions constructed, described above. The six unique components of the fibre stress term τ i j f are then computed as
τ i j f = ( a i j k l 1 3 δ i j a k l ) S k l
and stored in the six first UDS variable slots. The UDF then computes the three stress divergences terms τ i j f / x j with the gradients of τ i j f obtained from the Fluent solver and stores them in the remaining three UDS variable slots. Three additional UDFs were then used to add the three terms as momentum sources in the respective momentum equations.
The expression for μ f , Equation (4), is rewritten to a function of fibre volume fraction C V and fibre aspect ratio β . Using that, the volume concentration is the number of fibres per unit volume times the volume of a fibre:
C V = n π d 2 4 2 l = n π 2 ( 2 l β ) 2 2 l = 2 π n l 3 β 2
an expression for n l 3 as a function of C V and β reads:
n l 3 = C V β 2 2 π
so that Equation (4) may be written as
μ f = μ 4 β 2 C V 3 [ l n ( 1 C V ) + l n ( l n ( 1 C V ) ) + 1.4389 ]
Since information on either β and C V or n and l may be missing in some cases, the model may contain uncertainties. It was therefore decided to introduce a constant, called D f , to account for this. The constant D f > 0 was introduced in the expression for additional viscosity due to fibres and Equation (24) was thus modifies as
µ f = µ   .   D f 4 β 2 C V 3 [ l n ( 1 C V ) + l n ( l n ( 1 C V ) ) + 1.4389 ]

3.4. Comparison with Bingham Viscoplastic Fluid Model

The performance of the ODF model with an existing fibre suspension model was compared to the Bingham viscoplastic fluid model.
In the non-Newtonian rheology described by the model, the effective viscosity is infinite for low stress levels [22]. A more numerically favourable implementation was used by Fock (2010) [23] to simulate the laminar fibre suspension flow using the Herschel–Bulkley model in Fluent with the power law index set to 1. The effective viscosity was then computed as (ANSYS Inc. (b), 2013) [21]:
μ = { τ y γ ˙ c ( 2 γ ˙ γ ˙ c ) + k ,    γ ˙ γ ˙ c   τ y γ ˙ + k ,    γ ˙ > γ ˙ c
where τ y is the yield stress, γ ˙ the local shear rate, γ ˙ c the critical shear rate and k is called the consistency index. The local shear rate is defined as Fock (2010) [23]
γ ˙ = S i j S i j
and the critical shear rate is defined as
γ ˙ c = τ y μ 0
where μ 0 is a model parameter that should be large compared to the viscosity of the suspending water. The yield stress was computed by the empirical relation [4,6,23]:
τ y = a C m b
where a and b are model constants that vary from suspension to suspension. The Herschel–Bulkley in the Fluent model is not available for use alongside turbulence models, which was desired. In the present work, a UDF was therefore created to compute the local viscosity as a function of the shear rate, which allows for the use of the Bingham model with turbulence models.

3.5. Simulations

Simulations with the new ODF model was performed on the geometry (see Figure 3), by applying a backward facing step method in order to validate the computed results through comparison with experimental data obtained by Claesson (2012) [24] for softwood fibres.
The impact of changes in consistency or flow rate on flow structure after a backward facing step was investigated by comparing the reattachment length, at different fluid and flow properties. When the reattachment length was investigated, the point at the wall where the velocity changed from negative to positive was compared as Claesson (2012) [24].
To minimise the introduction of modelling errors in addition to the ones from the fibre modelling, an accurate turbulence model was desired. The computational resources were, however, limited to a single laptop workstation with a four core Central Processing Unit CPU and direct numerical simulation (DNS) was therefore considered unfeasible. Wall-resolved large eddy simulation (LES) was also considered unaffordable as fine spatial and temporal resolution was needed in the near-wall region (Davidson, 2014) [25]. Detached eddy simulations (DESs) were instead used, where the turbulent boundary layer was resolved by unsteady Reynolds-averaged Navier–Stokes equations RANS with the k ω Shear Stress Transport SST model and the outer region was resolved by LES (Davidson, 2014) [25], as the model was considered sufficiently accurate.
With the geometry resembling the experiments of Claesson (2012) [24], a numerical mesh was generated accordingly. These simulations were conducted for pure water. The computed recirculation zone results were compared to the experiments. The level of detail of the mesh was then gradually increased and the procedure was repeated. A total number of three meshes were constructed and the final one was identified to be sufficient enough. The three meshes and the corresponding computed recirculation zones u 0 , plotted with the experimental curves u = 0 from Claesson (2012) [24] where u is the streamwise velocity, can be seen in Figure 3. It can be observed that the changes in geometry have a clear impact on the computed results for the flow case in hand.
All meshes were generated with prism boundary layers and polyhedral cells of sizes suitable for DES. A symmetry boundary was used in the middle of the channel along the streamwise direction, which effectively reduced the number of elements by half. The final mesh consisted of approximately 180,000 cells. A grid independence study was conducted on the mesh using simulations of pure water flow. A fine mesh with approximately 434,000 cells and a coarse mesh with approximately 34,500 cells was created. The study showed that the 180,000 cell mesh had sufficient resolution to resolve the flow with the DES model.
Simulations with the ODF model were designed to match the conditions in Claesson (2012) [24]. As data on experimental temperature conditions were not available, the simulations were carried out for water with viscosity μ = 1.005 mPa and density ρ = 998.2 kg / m 3 , corresponding to a standard case with a temperature of 20 °C (Mörtstedt and Hellsten, 2010). Furthermore, the mass concentration C m = 0.015 was used, and the aspect ratio β = 50 was used as a suitable value since softwood fibres have length of the millimetre scale and a diameter of a few microns [6].
The mass concentrations for the different flow cases were given in Claesson (2012) [24]. However, the unknown volume concentration C V was needed to compute μ f . An approximation was therefore employed. A relation between volume concentration and mass concentration can be found in Derakhshandeh (2011) [6], reading:
C V = C m ( 1 ρ f + X w ρ + V L ) ρ b
where ρ f is the fibre density, X w the amount of water absorbed in the fibre, V L the volume of the hollow channel in the fibres per unit length and ρ b the bulk density of the fibre suspension, which is approximately equal to the water density because of the low concentration. Now, since the main part of these quantities are unknown, it was assumed that they could be included in an overall approximation of the fibre density, such that:
C V ρ ρ f C m
Claesson (2012) [24] used never-dried softwood pulp from a Swedish paper mill consisting of a mixture of fibres from pine and spruce. The values of dry densities of wood from pine and spruce are 290   kg / m 3 and 285   kg / m 3 , respectively [26]. However, since the fibres were soaked by a significant amount of water, it was approximated that their water content was around 50%. The fibre density was calculated as the average of the wood’s dry density and the density of water, which is approximately equal to 650 kg / m 3 . The fibre density was therefore approximated as the average of the dry density and the density of water as
ρ f =   1 2 ( ρ f + ρ w ) 1 2 ( 1000 + 290 ) k g m 3 = 650   k g m 3
Different parameter settings were used for simulations with the ODF model, both with the ODF standard deviation constant and varying linearly with turbulent intensity I . The constant D f was also introduced, as shown in Table 1, and was introduced in the expression for additional fibre viscosity μ f   (Equation (24)). The aim was to investigate the sensitivity of changing the fibre stress term or, in extension, the sensitivity of modifying the fibre data, such as aspect ratio and volume concentration. The full set of parameter settings used for the ODF model is shown in Table 1. A similar parameter study was performed for the Bingham model and the ideal setting was then used for comparison with the ODF model.
To avoid overfitting, due to only validating the fibre models with one experimental flow, the same were further validated by simulating a turbulent channel flow. Computed results were then compared to the experimental data from Xu and Aidun (2005) [27] for softwood fibres. The turbulent flow was modelled with DES in the same way as for the backwards facing step case. The flow conditions and fibre data used were synchronized with the experimental conditions in Xu and Aidun (2005) [27], where the fibre used in the experiments was natural flexible cellulous wood fibre with an average length of 2.3 mm and an average diameter of about 35 μm and the channel dimensions were 150 cm long, 5.08 cm wide and 1.65 cm high.
The simulations with the Bingham model over the backwards facing step were initialized by starting from stabilised transient simulations of pure water using DES with k ω SST turbulence model. The validation of the model was conducted in two separate steps. In the first step, three combinations of the parameters a and b were tested and chosen based on the results in Derakhshandeh et al. (2010) [28], to identify the best combination compared to the experimental data. In the first step the constant µ0 was kept with a same value as Ford et al. (2006) [9]. In the second step of the validation, the effect on varying µ0, by changing the same by a factor of approximately 2, were studied, combined along with the best combination of the parameters a and b and in the first step. In total, five parameter settings were used in the two steps, which are listed in Table 2. The two validation steps are summarized in short below:
  • Vary a and b;
  • Vary µ0.

4. Results

For the Bingham model, different parameter settings were tested based on the literature [9,28,29,30,31]. The final and best performing was the one used to compare the ODF model with the parameters a = 222,000   Pa , b = 1.95 , μ 0 = 100   Pa s . The best performing parameter setting for the ODF model was the one where the standard deviation was computed with the slope 1.25 and D f was set to 1 , i.e., setting R6 according to Table 1. It could be seen though, that changing the settings both for the ODF standard deviation σ and the constant D f gave a notable impact on the results.
The computed recirculation zones are visualised by the data points where the streamwise velocity fulfils u 0 and the experimental zone is visualised by the level curve u = 0 . In order for A, the sketch of the lines and the recirculation zone area are shown in Figure 4.
The experimental data in Claesson (2012) [24] were normalized with and were defined as the mean streamwise velocity at the step and at height. As a remark, the profiles were only plotted up to the height in the same way as the experimental profiles.
As per Figure 4, the computed results from the backwards facing step flow were then visualised by plotting the velocity profiles along a set of wall-normal lines near the centre of the channel downstream of the step, that matched the measurement positions of Claesson (2012) [24]. The recirculation zone, i.e., the zone where the streamwise velocity was negative, was also plotted. All computed results and experimental data were obtained for fibre suspension with mass consistency 1.5 % .
The best parameter settings for the respective models were identified by comparing the computed profiles and recirculation zones with the experimental ones. In Figure 5, the computed velocity profiles and recirculation zone is shown along with the corresponding experimental data with the best parameter setting for the ODF model and in Figure 6 for the Bingham model. The velocity profiles show the streamwise velocity U normalised by the bulk velocity U B , hereby defined as the mean streamwise velocity at the step and at a height of 20 mm .
It can be seen from Figure 5 and Figure 6 that for both the models, the experimental velocity profiles adhere reasonably well for the first part downstream up to x / h 5 . Further downstream, the obtained computational results do not synchronize with that of experimental results. The difference from the experimental data downstream x / h > 5   may be due to: a. geometry and b. fibre modelling errors, where the flow downstream is likely to be more turbulent in nature. It can also be noted that, further downstream, the turbulence decay mechanisms such as for example, flocculation, can have significant effects on the real flow pattern. Such mechanisms are poorly captured by the fibre models. Therefore, it can be expected that the fibre models deviate from the experimental results in an intensely turbulent region. This can be validated by judging from the velocity profiles shown in Figure 5 and Figure 6. Looking at the recirculation zones, in Figure 5 and Figure 6, it can be noted that the results obtained by the Bingham model simulations adheres better than the results from the ODF model compared to the experiments. This is due to the fact that the ODF model results slightly underestimate the length of the recirculation zone.
In Figure 7 are shown the computed recirculation zones for pure water, the Bingham model and for the ODF model. The results are compared to the experimental recirculation zone curve for a fibre mass consistency of 1.5%. It can be noted from the same that the fibre models (non-Newtonian fluids) have a better numerical prediction compared to the flow of pure water (Newtonian fluid).
Simulations were also performed for a turbulent channel flow to further validate models with their respective best settings from the backwards facing step simulations. Velocity profiles and turbulent intensity profiles for the ODF model, the Bingham model and the experimental data in Xu and Aidun (2005) [27] can be seen for Reynolds number 37,000 in Figure 8 and for Reynolds number 92,000 in Figure 9. The Reynolds numbers are based on the properties of pure water. It can be observed that the experimental data are clearly better resembled for the larger Reynolds number, i.e., for the more turbulent flow. The result indicates that the observation from the channel flow, that the fibre models work for turbulent fluidised flow regimes but not as well for low turbulence flows, may be correct. For the larger Reynolds number, the ODF model and the Bingham model show practically identical velocity profiles and the turbulence intensity is even slightly closer for the Bingham model.
Furthermore, Figure 7, where computed recirculation zones for simulations with pure water (left), the Bingham model (middle) and the ODF model (right) are plotted along with the experimental recirculation zone curve for 1.5 % fibres from Claesson (2012) [24], Figure 8, the normalised velocity profiles (left) and turbulent intensity profiles (right) for fibre suspension consistency 1 % and Reynolds number 37,000 for the Bingham model ( ), the ODF model (-- × ) and experiments from Xu and Aidun (2005) [27] were shown, whereas Figure 9 shows normalised velocity profiles (left) and turbulent intensity profiles (right) for fibre suspension consistency 1 % and Reynolds number 92,000 for the Bingham model ( ), the ODF model (-- × ) and experiments from Xu and Aidun (2005) [27]. It is obvious that it is hard to identify explicitly why the turbulent flow modulation effect is due to the addition of fibers, and it is not clear that the phenomenon can be physically reproduced by comparing only with the experimental values of multiphase flow, as shown in Figure 8 and Figure 9, as the Bingham model (middle) and the ODF model (right) showed a more synchronized result compared to the single-phase flow as indicated in Figure 7 where computed recirculation zones for simulations with pure water (left) are shown. Hence, only multiphase flow was taken into consideration.
It can be discussed if the assumption of rigidity that is employed in the ODF model is at all correct for softwood fibres. One motivation for using the approach even if incorrect is the assumption of turbulent flow. It may be correct to use the ODF model if the flow is sufficiently turbulent, so that the fibre effects on the flow pattern mainly consist of the interaction between the fibres and the turbulent flow field, rather than complex fibre–fibre interaction mechanisms. Such effects are assumed to be of less importance because of the fluidisation of the suspension flow, as turbulent fluctuations rip any major flocs apart. In addition, by trimming in the model parameters with experimental data, a collective representation of all complex mechanisms that are not explicitly described by the model are also included in the applied model.
As shown in the Figure 9, the results obtained by applying both Bingham and ODF models deviated from the experiments in some areas. However, it can be concluded that the Bingham model predicted the flow slightly better than the ODF model. The Bingham model also provides both a simpler and a cheaper implementation. It may thus be questionable if there is a reason to use the more complex ODF model to simulate turbulent fibre suspension flow. It should be stressed, though, that only a very limited set of parameter settings for the ODF model could be studies due to the limited computational resources. A more rigorous parameter study where more parameter settings are simulated and a number of experimental flow cases are used for validation may be used to optimise the model even more, and the anisotropic stress description provided through the ODF and the orientation tensors could serve to provide an accurate description of turbulent fibre suspension flow. At this stage, however, an interesting result is that for a sufficiently turbulent flow, a simple model such as the Bingham model may suffice to describe a fibre suspension in certain applications.

5. Conclusions

A new model for turbulent fibre suspension flow was proposed. In the model, the fibres are coupled to the fluid momentum through the orientation distribution function (ODF). A model for the ODF was introduced based on empirical observations in the literature from numerical and experimental results. From the modelled ODF, explicit expressions for all components of the fibre orientation tensors were derived, which were then used to compute additional stress term that arises in the fluid momentum equations due to the presence of fibres. Compared to solving the transport equation for the ODF, a very cheap implementation was obtained by the methodology. The ODF model was validated with experimental data by running simulations for the flow over a backwards facing step, also identifying the parameter settings resembling the experimental data the best out of the ones used. Simulations were also performed using a Bingham plastic model for comparison. The results showed that the ODF model performed well for the turbulent areas of the flow and not as well for less turbulent areas. The same was observed for the Bingham model. Both models’ performance were similar but the Bingham model performed slightly better, at least given the limited set of parameter settings tried in the present work. The result was further validated by the simulation of a turbulent channel flow with both models. It can be concluded that the ODF model may have potential to provide both a relatively cheap and detailed description of the turbulent flow of a fibre suspension. However, more rigorous studies are needed to optimise the model. At the current stage, a simple model such as the Bingham model gives sufficient description of fibre suspension flow for certain applications if the suspension is turbulent and fluidised.

Author Contributions

Conceptualization, V.S., Ö.J. and L.-O.L.; methodology, V.S. and A.L.; software, V.S., A.L. and T.P.; validation, V.S.; formal analysis, V.S. and T.P.; investigation, V.S., A.L. and Ö.J.; resources, V.S., T.P. and Ö.J.; data curation, V.S.; writing—original draft preparation, V.S and A.L.; writing—review and editing, V.S., T.P. and Ö.J.; visualization, V.S.; supervision, V.S.; project administration, V.S.; funding acquisition, V.S. and Ö.J.; All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Åforsk and Energimyndigheten grant number 2013, 2014 and 42606-1.

Acknowledgments

This project was carried out at Luleå University of Technology and ÅF-Industry AB in collaboration with ÅForsk, Energimyndigheten (Swedish Energy Agency), Holmen, SCA and Stora Enso.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. IPCC. Climate Change. 2014. Available online: http://www.ipcc.ch/pdf/assessment-report/ar5/syr/SYR_AR5_LONGERREPORT_Corr2.pdf (accessed on 23 January 2015).
  2. Swedish Energy Agency. En Samlad Bild Över Energiläget i Sverige. 2014. Available online: http://www.energimyndigheten.se/Press/Nyheter/En-samlad-bild-over-energilaget-i-Sverige/ (accessed on 23 January 2015).
  3. European Commission. Progress towards Achieving the Kyoto and 2020 Objectives; Progress Report; European Comission: Brussels, Belgium, 2014. [Google Scholar]
  4. Kerekes, R.J. Rheology of fibre suspensions in papermaking: An overview of recent research. Nord. Pulp. Pap. Res. J. 2006, 21, 598–612. [Google Scholar] [CrossRef]
  5. Bennington, C.P.J.; Kerekes, R.J. Power requirements for pulp suspension fluidization. TAPPI J. 1996, 79, 253–258. [Google Scholar]
  6. Derakhshandeh, B.; Kerekes, R.J.; Hatzikiriakos, S.G.; Bennington, C.P.J. Rheology of pulp fibre suspensions: A critical review. Chem. Eng. Sci. 2011, 66, 3460–3470. [Google Scholar] [CrossRef]
  7. Lindström, S.B.; Uesaka, T. Simulation of the motion of flexible fibers in viscous fluid flow. Phys. Fluids 2007, 19, 113307. [Google Scholar] [CrossRef]
  8. Andrić, J.; Fredriksson, S.T.; Lindström, S.B.; Sasic, S.; Nilsson, H. A study of a flexible fiber model and its behavior in DNS of turbulent channel flow. Acta Mech. 2013, 224, 2359–2374. [Google Scholar] [CrossRef] [Green Version]
  9. Ford, C.; Ein-Mozaffari, F.; Bennington, C.P.J.; Taghipour, F. Simulation of mixing dynamics in agitated pulp stock chests using CFD. AlChE J. 2006, 52, 3562–3569. [Google Scholar] [CrossRef]
  10. Yang, W.; Shen, S.; Ku, X. A new model of turbulent fibre suspension and its application in the pipe flow. Can. J. Chem. Eng. 2013, 91, 992–999. [Google Scholar] [CrossRef]
  11. Zhang, F. Instability in Settling Fibres: A Numerical Study. Ph.D. Thesis, Engineering Mechanics, Royal Institute of Technology KTH Department of Mechanics, Stockholm, Sweden, 2014. [Google Scholar]
  12. Krochak, P.J.; Olson, J.A.; Martinez, M. Fiber suspension flow in a tapered channel: The effect of flow/fiber coupling. Int. J. Multiph. Flow 2009, 35, 676–688. [Google Scholar] [CrossRef]
  13. Moosaie, A.; Manhart, M. Direct Monte Carlo simulation of turbulent drag reduction by rigid fibers in a channel flow. Acta Mech. 2013, 224, 2385–2413. [Google Scholar] [CrossRef]
  14. Lin, J.Z.; Liang, X.Y.; Zhang, S.L. Numerical simulation of fiber orientation distribution in round turbulent jet of fiber suspension. Chem. Eng. Res. Des. 2012, 90, 766–775. [Google Scholar] [CrossRef]
  15. Lin, J.-Z.; Sun, K.; Lin, J. Distribution of orientations in fibre suspension flowing in a turbulent boundary layer. Chin. Phys. Lett. 2005, 22, 3111. [Google Scholar]
  16. Lin, J.; Liang, X.; Zhang, S. Fibre orientation distribution in turbulent fibre suspensions flowing through an axisymmetric contraction. Can. J. Chem. Eng. 2011, 89, 1416–1425. [Google Scholar] [CrossRef]
  17. Olson, J.A.; Frigaard, I.; Chan, C.; Hämäläinen, J.P. Modeling a turbulent fibre suspension flowing in a planar contraction: The one-dimensional headbox. Int. J. Multiph. Flow 2004, 30, 51–66. [Google Scholar] [CrossRef]
  18. Shaqfeh, E.S.G.; Fredrickson, G.H. The hydrodynamic stress in a suspension of rods. Phys. Fluids 1990, 2, 7–24. [Google Scholar] [CrossRef]
  19. Rice, J.A. Mathematical Statistics and Data Analysis, 3rd ed.; Thomson Brooks/Cole: Belmont, CA, USA, 2007. [Google Scholar]
  20. Nordling, C.; Österman, J. Physics Handbook for Science and Engineering, 84th ed.; Studentlitterature AB: Lund, Sweden, 2006. [Google Scholar]
  21. ANSYS Inc. ANSYS Fluent User’s Guide; Canonsburg; ANSYS Inc.: Canonsburg, PA, USA, 2013. [Google Scholar]
  22. Irgens, F. Rheology and Non-Newtonian Fluids; Springer: Berlin/Heidelberg, Germany, 2014. [Google Scholar]
  23. Fock, H.; Rasmuson, A.; Wikström, T. CFD modeling of non-newtonian fluid mixing accounting for transient changes in local solids concentration—Application to an agitated pulp stock chest. Nord. Pulp Pap. Res. J. 2010, 25, 56–64. [Google Scholar] [CrossRef]
  24. Claesson, J.; Wikström, T.; Rasmuson, A. Flow of concentrated fiber suspensions over a backward facing step studied using LDA. Nord. Pulp. Pap. Res. J. 2012, 27, 653–661. [Google Scholar] [CrossRef]
  25. Davidson, L. Fluid Mechanics, Turbulent Flow and Turbulence Lecture Notes; Department of Applied Mechanics, Chalmers University of Technology: Gothenburg, Sweden, 2014. [Google Scholar]
  26. Mörtstedt, S.-E.; Hellsten, G. Data och Diagram, Energi-och Kemitekniska Tabeller, 7th ed.; Liber AB: Stockholm, Sweden, 2010. [Google Scholar]
  27. Xu, H.; Aidun, C.K. Characteristics of fiber suspension flow in a rectangular channel. Int. J. Multiph. Flow 2005, 31, 318–336. [Google Scholar] [CrossRef]
  28. Derakhshandeh, B.; Hatzikiriakos, S.G.; Bennington, C.P.J. The apparent yield stress of pulp fiber suspensions. J. Rheol. 2010, 54, 1137–1154. [Google Scholar] [CrossRef]
  29. Cotas, C.; Garcia, F.; Rasteiro, M.G.; Asendrych, D. Modelling of concentrated fibre suspension pipe flow with low-Reynolds-number k-ε turbulence models: New damping function. Nord. Pulp. Pap. Res. J. 2017, 32, 132–147. [Google Scholar] [CrossRef]
  30. Pamidi, T.R.K.; Johansson, Ö.; Löfqvist, T.; Shankar, V. Comparison of two different ultrasound reactors for the treatment of cellulose fibers. Ultrason. Sonochem. 2020, 62, 104841. [Google Scholar] [CrossRef]
  31. Zheng, E.; Rudman, M.; Kuang, S.; Chryss, A. Turbulent coarse-particle suspension flow: Measurement and modelling. Powder Technol. 2020, 373, 647–659. [Google Scholar] [CrossRef]
Figure 1. Angular components of the spherical coordinates used to describe fibre orientation.
Figure 1. Angular components of the spherical coordinates used to describe fibre orientation.
Fluids 05 00175 g001
Figure 2. Numerically integrated values ( ) and Fourier series values (--) of F ( ϕ ¯ ) (left) and G ( θ ¯ ) (right) of the fourth order orientation tensor component a 1123 . Orientation distribution function (ODF) standard deviation σ = 0.1   r a .
Figure 2. Numerically integrated values ( ) and Fourier series values (--) of F ( ϕ ¯ ) (left) and G ( θ ¯ ) (right) of the fourth order orientation tensor component a 1123 . Orientation distribution function (ODF) standard deviation σ = 0.1   r a .
Fluids 05 00175 g002
Figure 3. Backwards facing step mesh with increasing geometry detail (left) and the corresponding computed and experimental recirculation zones for pure water.
Figure 3. Backwards facing step mesh with increasing geometry detail (left) and the corresponding computed and experimental recirculation zones for pure water.
Fluids 05 00175 g003
Figure 4. Two-dimensional model of the backwards step with the lines used to compare the computed profiles with experimental profiles (−) and the area where the recirculation zone was plotted (--).
Figure 4. Two-dimensional model of the backwards step with the lines used to compare the computed profiles with experimental profiles (−) and the area where the recirculation zone was plotted (--).
Fluids 05 00175 g004
Figure 5. Velocity profiles (left) and recirculation zone (right) for the best ODF model setting. Experimental results from Claesson (2012) [24] are marked (--).
Figure 5. Velocity profiles (left) and recirculation zone (right) for the best ODF model setting. Experimental results from Claesson (2012) [24] are marked (--).
Fluids 05 00175 g005
Figure 6. Velocity profiles (left) and recirculation zone (right) for the best Bingham model setting. Experimental results from Claesson (2012) [24] are marked (--).
Figure 6. Velocity profiles (left) and recirculation zone (right) for the best Bingham model setting. Experimental results from Claesson (2012) [24] are marked (--).
Fluids 05 00175 g006
Figure 7. Computed recirculation zones for simulations with pure water (left), the Bingham model (middle) and the ODF model (right), plotted with the experimental recirculation zone curve for 1.5 % fibres from Claesson (2012) [24].
Figure 7. Computed recirculation zones for simulations with pure water (left), the Bingham model (middle) and the ODF model (right), plotted with the experimental recirculation zone curve for 1.5 % fibres from Claesson (2012) [24].
Fluids 05 00175 g007
Figure 8. Normalised velocity profiles (left) and turbulent intensity profiles (right) for fibre suspension consistency 1 % and Reynolds number 37,000 for the Bingham model ( ), the ODF model (-- × ) and the experiments (o) from Xu and Aidun (2005) [27].
Figure 8. Normalised velocity profiles (left) and turbulent intensity profiles (right) for fibre suspension consistency 1 % and Reynolds number 37,000 for the Bingham model ( ), the ODF model (-- × ) and the experiments (o) from Xu and Aidun (2005) [27].
Fluids 05 00175 g008
Figure 9. Normalised velocity profiles (left) and turbulent intensity profiles (right) for fibre suspension consistency 1% and Reynolds number 92,000 for the Bingham model ( ), the ODF model (-- × ) and the experiments (o) from Xu and Aidun (2005) [27].
Figure 9. Normalised velocity profiles (left) and turbulent intensity profiles (right) for fibre suspension consistency 1% and Reynolds number 92,000 for the Bingham model ( ), the ODF model (-- × ) and the experiments (o) from Xu and Aidun (2005) [27].
Fluids 05 00175 g009
Table 1. Parameter settings used for the ODF model in the simulations in the simulations of backwards facing step flow.
Table 1. Parameter settings used for the ODF model in the simulations in the simulations of backwards facing step flow.
SettingR1R2R3R4R5R6R7R8
σ 0.10 0.25 0.05 2.5 I 5.0 I 1.25 I 1.25 I 1.25 I
D f 1.0 1.0 1.0 1.0 1.0 1.0 2.0 0.5
Table 2. Parameter settings used for the Bingham model in the simulations of the backwards facing step flow.
Table 2. Parameter settings used for the Bingham model in the simulations of the backwards facing step flow.
SettingB1B2B3B4B5
A (Pa)810,000495,000222,000222,000222,000
B2.502.331.951.951.95
µ0 (Pa.s)10010010050200
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Shankar, V.; Lundberg, A.; Pamidi, T.; Landström, L.-O.; Johansson, Ö. CFD Analysis of Turbulent Fibre Suspension Flow. Fluids 2020, 5, 175. https://doi.org/10.3390/fluids5040175

AMA Style

Shankar V, Lundberg A, Pamidi T, Landström L-O, Johansson Ö. CFD Analysis of Turbulent Fibre Suspension Flow. Fluids. 2020; 5(4):175. https://doi.org/10.3390/fluids5040175

Chicago/Turabian Style

Shankar, Vijay, Anton Lundberg, Taraka Pamidi, Lars-Olof Landström, and Örjan Johansson. 2020. "CFD Analysis of Turbulent Fibre Suspension Flow" Fluids 5, no. 4: 175. https://doi.org/10.3390/fluids5040175

APA Style

Shankar, V., Lundberg, A., Pamidi, T., Landström, L. -O., & Johansson, Ö. (2020). CFD Analysis of Turbulent Fibre Suspension Flow. Fluids, 5(4), 175. https://doi.org/10.3390/fluids5040175

Article Metrics

Back to TopTop