Simulation of non-dilute fibre suspensions using RBF-based macro–micro multiscale method

The multiscale stochastic simulation method based on the marriage of the Brownian Configuration Field (BCF) and the Radial Basis Function mesh-free approximation for dilute fibre suspensions by our group, is further developed to simulate non-dilute fibre suspensions. For the present approach, the macro and micro processes proceeded at each time step are linked to each other by a fibre contributed stress formula associated with the used kinetic model. Due to the feature of non-dilute fibre suspensions, the interaction between fibres is introduced into the evolution equation to determine fibre configurations using the BCF method. The fibre stresses are then determined based on the fibre configuration fields using the Phan–Thien–Graham model. The efficiency of the simulation method is demonstrated by the analysis of the two challenging problems, the axisymmetric contraction and expansion flows, for a range of the fibre concentration from semi-dilute to concentrated regimes. Results evidenced by numerical experiments show that the present method would be potential in analysing and simulating various suspensions in food and medical industries.


Introduction
The regime of a suspension flow is based on two parameters: the aspect ratio a r and the volume fraction . Concretely, a suspension can be either dilute, semi-dilute, or concentrated for a 2 r < 1 , 1 < a 2 r < a r , or a 2 r > a r , respectively. In dilute suspensions, the interaction of fibres can be neglected and the fibre evolution can be captured by the Jeffery equation [1]. Meanwhile, the physical description of the evolution of suspension configurations poses a challenge due to the necessity to take into account fibre interactions for the semi-dilute and concentrated suspensions. The fibre parameters for cases ranging from semi-dilute to concentrated suspensions are presented in Table 1.
Besides the bulk properties of a flow, the orientation of fibres in the flow is also considered. At a position in the flow, the orientation of fibres is illustrated by an ellipse's geometry with three cases: a circle, an ellipse, or a straight line (see Fig. 1). They are corresponding to a predominant direction of fibres in parallel with the ellipse's long major axis, or all fibres aligning with the line at a position, respectively.
For non-dilute concentrated suspensions, the fibre-fibre interaction is significant. Thus, it is necessary to take into consideration this interaction and one possible way is to introduce a diffusion term into Jeffery's equation Folgar and Tucker [2]. From the literature, the simulation of a fiber suspension is basically carried out through three following steps: (i) introduce a fiber stress component into the momentum conservation equation to include dynamic effects of fibers on the bulk properties of the flow, (ii) apply an appropriate motion equation to describe the evolution of fiber particles, as stated above, the Jeffery's equation is suitable for dilute suspension, whereas the Folgar-Tucker's equation is applicable for semi-dilute and concentrated ones; and (iii) determine the fiber contribution to stress (named fiber stress tensor) using a relevant constitutive equation as a function of the fibers' orientation.
Since the fiber stress tensor is essentially calculated from the fourth-order orientation tensor ⟨ ⟩ , the basic 1 3 difference of numerical methods to simulate fiber suspensions is the way to handle the fourth-order tensor. There are several approaches to process the fourth-order orientation tensor. One approach is to use a quadratic closure approximation to break the tensor ⟨ ⟩ into two second-order tensors ⟨ ⟩ Lipscomb et al. [3]. A fully alignment assumption is subsequently applied to calculate these second-order tensors. This approach was employed to successfully simulate fiber suspension flows through axisymmetric contraction and expansion geometries Lipscomb et al. [3], Chiba et al. [4] and Baloch and Webster [5]. Another one is to directly solve the evolution of the fourth-order tensor as presented in Advani and Tucker [6]. And last but not least, the Brownian Configuration Field (BCF) approach Hulsen, et al. [7], Tran-Canh and Tran-Cong [8,9] was successfully applied to the simulation of fiber suspensions by Fan et al. [10], Lu et al. [11], Dou et al. [12] and Eberle et al. [13]. Follow the approach, a high number of fiber configurational fields is initiated on each computational node and the fourth-order tensor will be averagely calculated. Over the last decade, the macro-micro multiscale methods based on the differentiated and integrated RBF approximations have been developed to simulate successfully a range of dilute polymer solutions and melts Tran et al., [14,15], Nguyen et al. [16,17]. In addition, a multiscale modelling based on the combination of IRBF, DAVSS and BCF idea has been also proposed for simulations of dilute fibre suspensions in Nguyen et al. [18]. Owing to the advantages of RBF-based high-order approximation schemes, the approach achieved high-order convergence and accuracy Nguyen et al. [16,17].
Aligning with this approach, high-order RBF-based Brownian configuration field (BCF) method for dilute fibre suspensions Nguyen et al. [17] by our group is further developed for non-dilute suspensions by introducing a diffusion term into the Jeffery equation to capture fibre-fibre interaction. For this scheme, the conservation equations using the vorticity-stream function form are discretised using the high-order integrated RBF scheme, whereas the fibre configurations governed by the Folgar-Tucker equation are advanced by the BCF approach. Such two macro-micro processes are coupled using the Phan-Thien and Graham model Phan-Thien and Graham [19] for the fibre stress.
Since the present method is the combination between the Stochastic Simulation Technique and a mesh-free numerical method, it is powerful for problems with moving boundaries, complex boundary or free surface and without the need of closed form constitutive equation. Results evidenced by numerical experiments show the present method would be potential in simulating and producing suspensions in food and medical industries.
This work is organised as follows: the governing equations in the dimensionless form are first presented in Sect. 2. An introduction of the discrete adaptive viscoelastic stress splitting (DAVSS) technique into the conservation equations is also summarised in this section. Section 3 presents a coupled macro-micro system of governing equations, followed by a short review of the present method together with its algorithm. The application of the present method in semi-dilute and concentrated suspensions is then demonstrated and discussed through several benchmark examples in Sect. 4.

Governing equations for non-dilute suspension flows
The dimensionless conservation equations for fibre suspension flows are given by [11] where u, p and t are the velocity, pressure and time, respectively; s = 2 the stress contribution of the Newtonian solvent; = 1 2 ∇ + (∇ ) T the rate of strain tensor; f the  In this work, the fibre stress f for both the semi-dilute and concentrated suspensions is determined using the modified Phan-Thien-Graham model as follows Phan-Thien and Graham [19] and Fan et al. [10]: where ⟨ ⟩ and ⟨ ⟩ are the second-and the fourth-order orientation tensors, respectively, in which is the unit direction vector of fibres. D r is the diffusion coefficient; ⟨( * )⟩ the average value of (*); and A, F and f ( ) are fluid parameters which are defined as functions of the volume fraction and the aspect ratio a r of the suspended fibres as follows [10]: where m is the maximum volume package and empirically determined as a linear function of the aspect ratio [20]: The second-and fourth-order orientation tensors, < > and < > in Eq. (3) are defined by where N f is the number of fibre configuration fields. More details can be found in Bird et al. [21,22]. The evolution of fibres' orientation in non-dilute suspensions is determined using the Folgar and Tucker equation Folgar and Tucker [2]: where D Dt (⋅) is the material time derivative of (.), the identity tensor and is the effective velocity gradient tensor and given by It is worth noting that the fibres' interaction is random collisions by the Brownian force (b) (t) with the properties ⟨ (b) (t)⟩ = 0 and ⟨ (b) (t+s) (b) (t)⟩ = 2D r (s) , where (s) is the Dirac delta function; D r = C i̇ the diffusion coefficient; ̇= √ 2( ∶ ) is the general strain rate and C i the interaction coefficient. In this work, C i is chosen as a constant as done in [10,11] for simplicity.
By introducing ( , t) = Q ( , t) , Eq. (7) is transformed into [23] where Q is the modulus of and √ 2C i̇d dt is the Brownian force (b) , a function of the Wiener process .

Vorticity-stream function formulation in the cylindrical coordinates
In this work, several 2-D problems are considered using the axisymmetric vorticity-stream function in the cylindrical coordinates (r, z: radial and axial directions). The relations between velocity ( u r , u z ), and vorticity and stream function Ψ are given by The stream function and vorticity transport equations can be then derived from Eqs. (1) and (2) using Eq. (10), respectively, as follows: where zz f , zr f , rz f and rr f are the components of the fibre/ suspension stress tensor f .
Using the DAVSS scheme, Eq. (12) is rewritten as follows Nguyen et al. [17]: Re where a is the adaptive viscosity and expressed by [11]: With = ∇ + (∇ ) T is twice the strain rate tensor.

Discrete governing equations of fibre configurations in 2-D axisymmetric coordinates
The effective velocity gradient is developed in 2-D axisymmetric coordinate (z, r) as follows: With Eq. (15), the evolution equation (9) is developed in z and r-directions as Using tensor products, the fibre stress tensor in Eq. (3) can be developed in the z-and r-coordinates as follows: rr

Numerical method of the present method for non-dilute fibre suspension flows
A macro-micro multiscale system included Eqs. (11), for the governing equations together with the stress formula in 2-D axisymmetric coordinates is discretised using the high-order radial basis function-based BCF method Nguyen et al. [17] as follows.
The vorticity equation (13) is temporally discretised using the semi-implicit scheme as follows: For the microscale, the evolution Eq. (16a & b) are explicitly discretised using the Euler-Maruyama method: where superscripts (n + 1) and n denote the two successive time steps at t n+1 = (n + 1)∆t and t n = n∆t, respectively; ∆t is the time step size.
At each time step, the vorticity and/or stream function equations, as well as the first and second derivatives of the field variables including the fibre stresses at the collocation points, are approximated using the 1-D integrated RBF (1D-IRBF) method which is presented in Sect. 3.1.

The 1D-IRBF-based spatial discretisation scheme
At a time t, the highest order derivative of dependent variable ω(x, t) (the second order in this work), the first-order derivatives and the function itself are decomposed as follows Mai-Duy et al. [24] and Thai-Quang et al. [25]: where w j (t) m j=1 is the RBF weights; g j (x) m j=1 the RBFs; m a c h o s e n n u m b e r ; G [1] j (x) = ∫ G [2] j (x)dx ; j (x) = ∫ G [1] j (x)dx and C 1 and C 2 are unknown integration constants at time t. In this work, the multi-quadric RBF (MQRBF) is used and given by where c j m j=1 and a j m j=1 are the RBF centres and widths, respectively. The centres are chosen to be the same as the data points x j in this work.
Equations (20a), (20b) and (20c) are evaluated at each and every collocation point and then rearranged to produce the following a set of algebraic equations: where Owing to the presence of integration constants, more additional constraints can be incorporated into the algebraic equation system through Eq. (22c) as follows: and ̂=̂̂ are additional constraints.
The conversion of the network-weight space into the physical space yields where ̂ −1 is the conversion matrix. Equation (24)

Algorithm of the present macro-micro multiscale method
The algorithm of the present multiscale method is detailed as follows.
• The fibre stresses are determined at a time step t i+1 from the fibre configuration at step t i (the initial one for the first step) using Eq. (17a, b & c), and the derivatives of stresses are then approximated using the high-order integral RBF technique; • Solve Eqs. (18) and (11) for the vorticity and stream function, respectively. Then, calculate the velocity of the current step using Eq. (10);

Numerical examples
The present method is used to simulate two challenging nondilute fibre suspension flows: flows through 4:1 and 4.5:1 axisymmetric contractions, and flows through 1:4 axisymmetric expansion for several fibre parameters and Reynolds numbers. The fibre parameters for cases ranging from semidilute to concentrated suspensions are presented in Table 1 of the introduction.

Flow through a circular tube
Although the flow through a circular tube of non-dilute fibre suspensions is not a focus of this work, this numerical example is the first investigation, because the velocity profile, the vorticity and stream function at the outlet by the simulation will be used as the boundary conditions at the inlet of the fibre suspension flow through an axisymmetric contraction and expansion presented in Sect. 4.2.
The flow through a circular tube is described in Fig. 2. The length and radius of the tube are L = 10 and R = 0.5 , respectively. This simulation is carried out with Re = 0 , Δt = 1E − 3 and the number of fibres N f = 1000.
The boundary conditions are given as follows: (i) non-slip boundary condition on the wall BC: u z = 0 and u r = 0 ; (ii) Experiences find that finer meshes near the centreline and the outlet are necessary for higher solution accuracy in these regions. Therefore, a designed grid is applied with To demonstrate the role of fibre parameters in the kinetic behaviour of the flow, the non-dilute fibre suspension flow is simulated for a range of a 2 r . The obtained results by the present method are discussed and compared with others as follows.
• Figure 3 presents the distribution of the axial velocity u z along the centreline of the tube for a range of a 2 r = {1, 5, 8, 10, 15, 20} . Their corresponding volume fractions ( ) and aspect ratios ( a r ) can be found in Table 1. Results show that in all cases of a 2 r , undershoot is observed near the entrance (z ∈ (0.41, 0.74)) . Moreover, the undershoot is more pronounced with increasing a 2 r . Undershoots reflect the effect of the isotropic configuration of fibres (see Fig. 1a) at the inlet. In fact, the isotropy of fibre configurations resists the development of the velocity (u z ) at the region near the inlet. The velocity then increases along the centreline towards the outlet with a gradual decrease of the isotropy of fibres' orientation.
Results in Fig. 3 also show that the velocity (u z ) on the centreline at the steady state reduces as and/or a r increases. • Figure 4 shows the impact of and a r on the profile of axial velocity (u z ) at several cross sections (z = {0, 0.5, 1, 2, 4, 10}) with a 2 r = 1 (semi-dilute) (Fig. 4a), a 2 r = 10 (concentrated) (Fig. 4b) and a 2 r = 20 (concentrated) (Fig. 4c). Numerical results show that the velocity profile becomes more plug-like with increasing and/or a r . The profiles of outlet velocity of flows for a range of a 2 r = {1, 5, 10, 15, 20} presented in Fig. 5 also depict that the velocity profiles are more plug-like for higher values of a r and/or . • Finally, the influence of a 2 r on the convergence measure CM for the velocity is reported in Fig. 6. Results show a significant influence of a 2 r on the convergence of the method, where the CM generally degenerates with increasing level of fibre concentration and/or aspect ratio.

The 4:1 axisymmetric contraction flows
The geometry of the 4:1 axisymmetric contraction flow is schematically described in Fig. 7 with the length and radius of the upstream tube L U and R U , respectively; the length and radius of the downstream L D and R D ; and the length of the salient vortex L v .
The contraction ratio ( ) and the dimensionless vortex length ( L * v ) are given by The boundary conditions for this problem are as follows: (i) At the inlet OA , the velocity, the vorticity and the stream function are the values at the outlet of the circular Poiseuille flows investigated in Sect. 4.1. Furthermore, Ψ z = 0 is also assigned; where ω w 1 , ω w 2 and ω w 3 are determined and updated using Eq. (11) with the known stream function at each time step; (iv) On the centreline OE , symmetric boundary condition for the velocity, i.e., The non-uniform grid designed in [17] is applied to the simulation of this example as follows:  Fig. 8).
Number of fibres N f = 1000 and the time step size Δt = 1E − 3 are applied in this simulation. However, a finer time step is used at the initial time for the numerical stability of the method in several cases of highly concentrated suspensions. In addition, unless otherwise stated, the simulation is carried out for the 4:1 contraction flow with L U = 6 ,R U = 1,L D = 5 , R D = 0.25 and Re = 0.
The effect of the volume fraction ( ) and the fibre aspect ratio ( a r ) on the flow pattern and fibres' orientation is first investigated. Results by the present method are discussed and compared with those by [11] using the BCF-finite-element method as follows.
- Figure 9 depicts the effect of the fibre parameters ( , a r ) on the vortex length ( L * v ) for two cases of fibre suspension: a r = 10 with a range of ∈ {0.01, 0.02, 0.05, 0.08, 0.10, 0.12, 0.15, 0.18, 0.20} a n d a r = 2 0 w i t h a r a n g e o f ∈ {0.01, 0.02, 0.05, 0.08, 0.10} . Numerical experiments show that the vortex length is more pronounced with the increase of or/and a r . This observation was reported in [11] but with a gap in the results by the present method as shown in Fig. 9.
-The impact of fibre parameters on the shape and length of the cortex at the salient corner for a range of a 2 r = {5, 10, 15, 20} is presented in Fig. 10 in which the salient vortex's size is more pronounced with increasing a 2 r . Furthermore, the present results also confirm that unlike the Newtonian flow, where the vortex boundary is concave, it is slightly convex for fibre suspension flows as reported in Lipscomb et al. [3], Chiba et al. [4] and Lu et al. [11]. - Figure 11 clearly shows the influence of a 2 r on the orientation of fibres around the abrupt contraction by an ellipse/line/circle (as defined in Fig. 1) at a position, whereas the black dashed lines represent the streamlines for four cases of a 2 r = {5, 10, 15, 20} . Several notable points are observed as follows.
• Fibre distribution in the main area of the flow tends to be less anisotropic with the increasing of a 2 r . Experimental results also depict that the region dominated by anisotropic fibres distribution is gradually shrinking from lower a 2 r = 5 for semidilute (Fig. 11a) to higher a 2 r = 20 for concentrated (Fig. 11d). This can be explained by the random interaction between fibres, which is stronger with increasing a 2 r , reduces the influence of the flow velocity on the fibres' orientation. • Meanwhile, the fibres distribution is more isotropic towards the centre of the vortex.
Finally, the present work also investigates the impact of Re on the vortex size for a 2 r = 12 . Numerical observation by Fig. 12 shows that the length of vortex lightly reduces with the increasing of Re for a range of Re ∈ {0, 1, 2, 5} . This is in good agreement with those in Abdul-Karem et al. [26].

The 1:4 axisymmetric expansion flow
The geometry of 1:4 axisymmetric expansion flow is schematically described in Fig. 13 in which L D = 4 and R D = 1 are the length and radius of the downstream tube, respectively; L U = 6 and R U = 0.25 the length and radius of the upstream tube, respectively; and L v the vortex length. The expansion ratio (β) and dimensionless vortex length ( L * v ) are defined as in Eq. (27) but for the expansion flow. The similar problem was studied by [11] using the finite-element method-based BCF. The boundary conditions of this problem are described as follows. where ω w 1 , ω w 2 and ω w 3 are determined using Eq. (11) with the known stream function at each time step; (iv) On the centreline OE , the symmetric boundary condition for the velocity is applied, i.e., A non-uniform grid is designed for the simulation. Finer grids in regions near the outlet and inlet and around the abrupt expansion and the centreline are necessary to produce accurate solutions. A detailed grid is described in Fig. 14  The time step size Δt = 1E − 3 is also used in this simulation as in the contraction flow problem. However, a smaller time step is used at the initial time for the simulation of several highly concentrated suspensions. Numerical results by the present method are discussed and compared with those by [11] and other publications as follows: The effect of fibre parameters ( , a r ) on the vortex length ( L * v ) of the 1:4 expansion flow is presented by Fig. 15 for two cases of fibre suspensions: (i) a r = 10 coupled with a range of ∈ {0.01, 0.02, 0.05, 0.08, 0.10, 0.12, 0.15, 0.18, 0.20} ; and τ rr f ) and the first normal stress difference ( τ rr f − τ zz f ) at the position z = 5.75 and r = 0.5 a n d ( i i ) a r = 20 c o u p l e d w i t h a r a n ge o f ∈ {0.01, 0.02, 0.05, 0.08, 0.10} . The numerical results by Fig. 15 are in very good agreement with those by [11]. Indeed, in the contraction flows the vortex length is more pronounced with increasing and/or a r (Fig. 10), while it changes insignificantly in the expansion flows ( L * v ≈ 0.18 … 0.26). Moreover, numerical results by the present method confirm that the vortex lengths in the expansion flows are smaller than those in the contraction flows with the same fibre parameters. However, while the vortex length is nearly unchanged, a small lip vortex is observed near the re-entrant corner since a sufficiently large value of a 2 r (15 & 20) as shown in Fig. 16 for a range of a 2 r = {5, 10, 15, 20} . This observation is notable, because lip vortex was reported only for the Newtonian and viscoelastic expansion flows with high expansion ratio and Weissenberg number by Baloch and Webster [5] and Baloch et al. [27] but not in any report on fibre suspension flows.
The orientation of fibres as well as streamlines (black dashed lines) around the expansion area for a range of a 2 r = {5, 10, 15, 20} is shown in Fig. 17. Results show that except fibres on the centreline, most fibres align with the flow direction (straight line or ellipse) in the upstream region, whereas the direction of fibres is more or less random along the streamline (a circle/circular ellipse) in the downstream region. In other words, near the vortex boundary fibres are not tangential to streamlines, which is in contrast to our numerical observation in the simulation of contraction flow (see Fig. 11). These findings are in good agreement with those of several numerical and experimental works Abdul-Karem et al. [26], Chiba and Nakamura [28] and Verweyst and Tucker [29].
The significant impact of inertia on the vortex size is also found by the present method. Indeed, results by Fig. 18 show a fast development of vortex with increasing Re = {0, 1, 2, 5} , which was reported in several experimental works by Abdul-Karem et al. [26] and Townsend and Walters [30]. Furthermore, the observation also shows, while the vortex boundary is convex referred to the centreline for the creeping flow ( Re = 0 , Fig. 18a), it is gradually concave for the flows at higher Re (Fig. 18b-d).
A simulation of the 1:4 expansion flow of non-dilute fibre suspensions with the fibre's aspect ratio a r = 10 and volume fraction = 0.15 using a range of the numbers of fibre configurations N f = {500, 1000, 2000} is also carried out to study the impact of the fibre configuration on the mechanical properties of suspension. Figure 19 presents the profiles of shear stress ( τ zr f ), normal stresses ( τ zz f and τ rr f ) and the first normal stress difference ( τ rr f − τ zz f ) at the position (z,r) = (5.75, 0.5) in the considered domain. The numerical results indicate that while τ rr f and τ zr f , about 12 (Fig. 18a) and 59 (Fig. 18c), respectively, are in good agreement with those of the similar problem by Verweyst and Tucker [11], τ zz f and τ rr f − τ zz f reported in Fig. 18b, d are smaller than those mentioned in the same publication. In addition, the smoothness of curve presenting fibre stresses, especially τ zz f , is improved with increasing the number of fibre configurations as described in Fig. 19.

Conclusions
The multiscale simulation method based on the combination of high-order integral RBF approximation and BCF idea is further developed for non-dilute fibre suspensions in this work. The research is to simulate semi-dilute and concentrated fibre suspension flows in which the evolution of fibre configurations governed by the Folgar-Tucker equation are determined by the BCF method and the fibre stress is approximated by the Phan-Thien-Graham model. The efficiency of the present method for the simulation of nondilute fibre suspension flows is based on both the accuracy of the method and the stability of the stochastic process. As an illustration of the present method, two challenging flows: the 1:4 axisymmetric expansion and 4:1 axisymmetric contraction flows are investigated for a large range of the fibre parameters from semi-dilute to concentrated regimes. Numerical results by the present method are in good agreement with those published by Lu et al. [11]. Moreover, several dynamics behaviours of suspension flows by the present work including the orientation of fibres are nearly similar to experimental observations published in Abdul-Karem et al. [26], Verweyst and Tucker [29] and Baloch et al. [5]. Finally, the present work finds the existence of a lip vortex in the 1:4 expansion flows of highly concentrated fibre suspensions.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.