CFD Analysis of Turbulent Fibre Suspension Flow

: 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


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.

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: where is the viscosity of the suspending fluid, the local strain rate, the additional viscosity due to the presence of fibres and is the unit tensor. The fibre orientation tensors and are the second and fourth order statistical moments of fibre orientation, respectively, and are defined as where Ψ is the orientation distribution function (ODF) and is the fibre orientation vector. The additional viscosity due to fibre presence depends on fibre properties and fibre concentration and may be computed for a semi-dilute suspension as Shaqfeh and Fredrickson (1990) [18] = 8 3 ln 1 + ln ln 1 + 1.4389 (4) where is the number density of fibres, is the half-length of fibres and the volume concentration of fibres. Semi-dilute suspensions are here referred to as those satisfying ≫ 1 and < 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] where is the velocity, the density of the suspending fluid and is pressure. It may be remarked that the body force term 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]: where is the fibre orientation vector, the fibre position and and are the gradient operators in orientation and translation space, respectively, and and 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: where μ is the viscosity of the suspending fluid, S the local strain rate, μ the additional viscosity due to the presence of fibres and δ is the unit tensor. The fibre orientation tensors a and a are the second and fourth order statistical moments of fibre orientation, respectively, and are defined as a = 〈p p 〉 = p p Ψ(p)dp , (8) a = 〈p p p p 〉 = p p p p Ψ(p)dp (9) where Ψ is the orientation distribution function (ODF) and p is the fibre orientation vector.

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 and , 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.

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; 3. Fibres may be treated as rigid so that the Equation (4) and Equation (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 and , 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]: 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 as ∼ N( ̅ , ) 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]: 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 as where is the slope, is a bulk velocity of the flow and 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.

Explicit a priori Orientation Tensors
For a given ODF variance and flow field direction ( , ̅ ) the modelled ODF is known and all orientation tensor components can be computed, as an example for ,: where the orientation vector components consist of a product of functions of and , respectively, so that the tensor can be written on the form ( ) ( ). The double integral can therefore be separated into two integrals and are then functions of the flow field angles and ̅ and are therefore denoted ( ) and ( ̅ ), respectively. The values for ( ) and ( ̅ ) 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 − , , which may be interpreted from An expression for the tensor component is then constructed as 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 can then be computed by linear interpolation between the constructed expressions. The full procedure is applied to all orientation tensor components and . As a remark, only a number of terms up to , = 3 is needed to obtain a good fit for the Fourier series. An example can be seen for the component 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 and do not need to be computed by numerical integration in each cell and iteration, as they are determined explicitly from the fluid flow. −

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 userdefined 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 and from the explicit expressions constructed, described above. The six unique components of the fibre stress term are then computed as and stored in the six first UDS variable slots. The UDF then computes the three stress divergences terms / with the gradients of 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 , Equation (4), is rewritten to a function of fibre volume fraction and fibre aspect ratio . Using that, the volume concentration is the number of fibres per unit volume times the volume of a fibre: an expression for as a function of and reads: so that Equation (4) may be written as Since information on either β and or and may be missing in some cases, the model may contain uncertainties. It was therefore decided to introduce a constant, called , to account for this. The constant > 0 was introduced in the expression for additional viscosity due to fibres and Equation (24) was thus modifies as µ = µ .

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]: where is the yield stress, ̇ the local shear rate, ̇ the critical shear rate and is called the consistency index. The local shear rate is defined as Fock (2010) [23] ̇= (27) and the critical shear rate is defined as where 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] where and 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.

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 − 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 ≤ 0, plotted with the experimental curves = 0 from Claesson (2012) [24] where 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.005mPa and density = 998.2kg/m , corresponding to a standard case with a temperature of 20 °C (Mörtstedt and Hellsten, 2010). Furthermore, the mass concentration = 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 was needed to compute . An approximation was therefore employed. A relation between volume concentration and mass concentration can be found in Derakhshandeh (2011) [6], reading: where is the fibre density, the amount of water absorbed in the fibre, the volume of the hollow channel in the fibres per unit length and 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: ≈ 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 and 285 kg/m , 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 650kg/m . The fibre density was therefore approximated as the average of the dry density and the density of water as = 1 2 + ≈ 1 2 (1000 + 290) = 650 (32) Different parameter settings were used for simulations with the ODF model, both with the ODF standard deviation constant and varying linearly with turbulent intensity . The constant was also introduced, as shown in Table 1, and was introduced in the expression for additional fibre viscosity (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 − 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: 1 . Vary a and b; 2. Vary µ0.

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 = 222,000Pa, = 1.95, = 100Pa ⋅ 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 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 gave a notable impact on the results. The computed recirculation zones are visualised by the data points where the streamwise velocity fulfils ≤ 0 and the experimental zone is visualised by the level curve = 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 normalised by the bulk velocity , hereby defined as the mean streamwise velocity at the step and at a height of 20mm.
It can be seen from Figures 5 and 6 that for both the models, the experimental velocity profiles adhere reasonably well for the first part downstream up to /ℎ ≈ 5 . Further downstream, the obtained computational results do not synchronize with that of experimental results. The difference from the experimental data downstream /ℎ > 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 Figures 5 and 6. Looking at the recirculation zones, in Figures 5 and 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.   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 Figures 8 and 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.

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.