- Email: [email protected]
Numerical simulation of stirred tanks using a hybrid immersed-boundary method
Numerical simulation of stirred tanks using a hybrid immersed-boundary method Shengbin Di, Ji Xu, Qi Chang, Wei Ge PII: DOI: Reference:
S1004-9541(16)30507-9 doi: 10.1016/j.cjche.2016.05.031 CJCHE 584
To appear in: Received date: Revised date: Accepted date:
29 August 2015 9 March 2016 18 March 2016
Please cite this article as: Shengbin Di, Ji Xu, Qi Chang, Wei Ge, Numerical simulation of stirred tanks using a hybrid immersed-boundary method, (2016), doi: 10.1016/j.cjche.2016.05.031
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
IP
T
Numerical simulation of stirred tanks using a hybrid immersed-boundary method*
1
SC R
Shengbin Di(狄升斌)1,2, Ji Xu(徐骥)1, Qi Chang(常麒) 1,3，Wei Ge(葛蔚)1,** State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese
Academy of Sciences, Beijing 100190, China 2
NU
University of Chinese Academy of Sciences, Beijing 100049, China
3
MA
School of Chemical Engineering and Technology, Tianjin University, Tianjin 300072, China
Article history: Received 29 August 2015
CE P AC
TE
Accepted 18 March 2016
D
Received in revised form 9 March 2016
ABSTRACT Conventionally, multiple reference frame (MRF) method and sliding mesh (SM) method are used in the simulation of stirred tanks, however, both methods have limitations. In this study, a hybrid immersed-boundary (IB) technique is developed in a finite difference context for the numerical simulation of stirred tanks. IBs based on Lagrangian markers and solid volume fractions are used for moving and stationary boundaries, respectively, to achieve optimal efficiency and accuracy. To cope with the high computational cost in the simulation of stirred tanks, the technique is implemented on computers with hybrid architecture where central processing units (CPUs) and graphics processing units (GPUs) are used together. The accuracy and efficiency of the present technique are first demonstrated in a relatively simple case, and then the technique is applied to the simulation of turbulent flow in a Rushton stirred tank with large eddy simulation (LES). Finally the
ACCEPTED MANUSCRIPT proposed methodology is coupled with discrete element method (DEM) to accomplish particle-resolved simulation of solid suspensions in small stirred tanks. It demonstrates that the
T
proposed methodology is a promising tool in simulating turbulent flow in stirred tanks with complex
IP
geometries.
SC R
Keywords: immersed-boundary method; CPU-GPU hybrid computing; stirred tank; large eddy
AC
CE P
TE
D
MA
NU
simulation
ACCEPTED MANUSCRIPT 1. Introduction Stirred tanks are widely used in chemical and process industries, and hence deserve experimental
and
computational
investigation.
Experimental
T
comprehensive
IP
measurement techniques such as laser-doppler velocimetry (LDV) and particle image
SC R
velocimetry (PIV), have significantly contributed to our understanding of the flow structure in stirred tanks. However, most experimental studies have limitations, such as time and cost, scale-up difficulties, and the difficulty of directly measuring some
NU
quantities like local velocity especially for high-viscosity fluids[1]. Computational fluid dynamics provides an alternative tool for unveiling the details of fluid flow and
MA
heat transfer in stirred tanks, but faces several challenges. The main difficulties lie in the relative motion between the impeller and vessel and the complex geometry of the
D
impeller [2]. To address these difficulties, several numerical methods have been
TE
developed.
CE P
Conventionally, boundary-conforming nonorthogonal grids are used to handle the complex geometries [3]. In the so-called body-fitted approaches [4], the grid lines follow the boundaries exactly and the boundary conditions can therefore be
AC
implemented directly. However, the generation of such grids is both difficult and time consuming, and it has to be repeated or modified during the simulation for moving boundaries. Furthermore, the discretized governing equations for body-fitted grids are much more complex than those for Cartesian grids, which increase the computational cost and decrease the numerical stability considerably. As a result, such a method is seldom directly applied to the simulation of stirred tanks. A non-inertial reference frame fixed to the moving body can be used to avoid mesh regeneration for moving boundaries, although the governing equations have to take a more complicated form. Multiple reference frames (MRF) [5] and sliding mesh (SM) models [6] are typical examples of this approach. They are among the most commonly used methods for the simulation of stirred tanks and available in mainstream commercial codes [7-10]. Both
ACCEPTED MANUSCRIPT methods solve the governing equations in two domains: an inner domain containing the rotating impeller and an outer domain containing the stationary tank wall and baffles.
T
In SM, the meshes are allowed to shear and slide to accommodate the relative motion
IP
at the domain interface, while in MRF, only a stationary solution of the flow field is provided for a given position of impeller with respect to the baffles [11]. The main
SC R
disadvantages of SM are the explicit acceleration terms added to the momentum equation and the related numerical problem at the domain interface. For complex
computationally very demanding [12].
NU
impeller designs and intermeshing extruders, both MRF and SM might be infeasible or
MA
Immersed boundary (IB) methods, first proposed by Peskin [13], are another category of methods for complex geometries that can be based on simple and regular Cartesian grids, and have received much attention recently [14-16]. IB methods have
D
higher computational efficiency when compared with body-fitted methods [16-19]. In
TE
IB methods, the Cartesian grids do not necessarily conform to the boundaries. The
CE P
effects of the boundaries on the fluid are represented by additional body forces in the momentum equations. Although the number of grid points increases more quickly than in body-fitted methods as the Reynolds number increases, the per-grid-point-operation
AC
count is notably reduced owing to the absence of grid transformation [20]. In recent years, IB methods have been applied in the simulation of complex fluid flow in stirred tanks because of their efficiency and capacity to deal with moving boundaries directly [2, 21, 22]. Unlike SM model, there is no constraint on the prescribed motion of the rigid bodies in IB methods [23]. Among various IB methods [15, 16, 24-26], two typical approaches are the bases of this study. In the first approach, the boundaries are represented by a set of Lagrangian points which are handled with direct forcing methods such as that proposed by Uhlmann [16]. No-slip boundary conditions are assigned to these Lagrangian markers. The positions of the Lagrangian markers are updated in each time step if the boundaries are moving. As only the surfaces of the immersed bodies should be
ACCEPTED MANUSCRIPT represented by the markers, the number of markers is relatively small compared with the total number of Cartesian grid points. Each Lagrangian point is connected with a
T
number of Cartesian grids around it. Mapping of the variables between Lagrangian and
IP
Eulerian locations is performed employing a regularized delta function. This technique effectively reduces non-physical force oscillations existing in moving-boundary
SC R
simulations[27], however, it is not necessary for stationary-boundary problems. Furthermore, when the IBs lie close to the borders of the computational domain, where
NU
the separation is smaller than the support width of the delta function, it is difficult to find enough Cartesian grid points for velocity interpolation and force distribution when
MA
using this technique.
In the second approach, the body force is evaluated directly at the IB using the desired velocity of the solid part. However, the fluid and solid velocities are smoothly
D
interpolated using a solid volume fraction at the boundary grid points. When the solid
TE
part is stationary, the solid fraction at the computational grid points does not change
CE P
during the computation. As a result, the evaluation of the solid fraction needs to be carried out only once. The method of Kajishima et al. [28] belongs to this category. It is noted that the two approaches described above are more suitable for moving and
AC
stationary boundaries, respectively (For simplicity, they are referred as the moving and stationary approaches hereafter). Both moving and stationary boundaries, such as the boundaries of rotating impellers and stationary tank walls, are present in a stirred tank simultaneously. However, most simulations use only one approach for all boundaries [22, 29, 30]. Tyagi et al. [30] proposed a method to combine the advantages of the IB method to represent moving rigid geometries and the multi-block structured curvilinear mesh to represent complex domains. However, it is still a challenge to generate body-fitted mesh for complex geometries such as the helical coils inside a tank. To this regard, we will propose in this work a hybrid IB technique in which the complex geometries are dealt with using different IB methods on a Cartesian grid.
ACCEPTED MANUSCRIPT On the other hand, one practical difficulty in simulating complex flows in stirred tanks is the high computation cost. Consequently, we implement the IB method on a
T
CPU-GPU hybrid architecture for massively parallel computing. Algorithms having
IP
good parallelism are therefore desirable. Since the fluid and IBs are discretized on two separate grids, their parallelization should be handled differently. The parallelization
SC R
for Eulerian grids is easy because the grid points are fixed in time and are assigned to the same subdomain for the whole computation period [31]. However, the Lagrangian
NU
markers move throughout the computational domain, so a given IB point may move across the subdomain border during one time step. Furthermore, the interactions
MA
between the fluid and IBs are another source of difficulty in the design of a distributed algorithm [32]. To cope with these problems, we developed a partial velocity interpolation and force distribution scheme, which is described in section 2..
D
To validate the proposed method, we first test its accuracy and robustness in two
TE
numerical experiments, incorporating large eddy simulation (LES) which is found to be
CE P
effective for turbulent flow in stirred tanks [33, 34]. Then, as an application, the discrete element method (DEM) is coupled with the present method to carry out direct numerical simulation (DNS) of fluid-solid systems in stirred tanks. Satisfactory results
AC
reveal that complex flows in stirred tanks can be accomplished using hybrid IB methods on supercomputers. It will become a powerful tool for optimal design of industrial stirred tanks if the method can be further developed. Finally, main conclusions are drawn with perspectives on wider applications of the present method.
2. Numerical method 2.1. Governing equations In the current work, we focus on incompressible fluid with constant properties. The governing equations of this method describe the conservation of momentum and mass of the flow. They are established for the fluid phase and the solid objects (or boundaries) separately, and coupled through the force terms describing the interactions between them.
ACCEPTED MANUSCRIPT 2.1.1.
Fluid-phase equations
The governing equations for unsteady incompressible flow in the IB methods take
u f u u p f 2 u + f t
(1) (2)
SC R
u 0
IP
f
T
the general form[16, 17]
where f and f are the fluid density and dynamic viscosity, respectively, while u, and p are the velocity and pressure, respectively. They are basically the Navier–Stokes
NU
equations with additional body force terms f determined by the no-slip boundary conditions at the fluid–solid interface. f is non-zero only in the vicinity of fluid–solid
2.1.2.
Solid-phase equations
MA
interfaces.
D
Only rigid bodies are considered in the present study. Given a prescribed motion, the
TE
desired velocity field within the solid region is [16]
ud ut x x0
(3)
CE P
where ut , and x0 are the translational and rotational velocity and center coordinates of a rigid body. Otherwise, if the motion of the rigid solid domain is due to
AC
the fluid force, the motion is controlled by Newton’s equations for linear and angular momenta of a rigid body, d ms vs dt d I s s dt
Fs Gs
(4)
Ts N s
(5)
where v s and s are the linear velocity at the mass center and the angular velocity of the rigid body, respectively. Fs and Ts represent the hydrodynamic force and torque acting on the rigid body due to the fluid. Gs and Ns are the external force and torque acting on the rigid body. ms is the mass of the rigid body and Is is the moment of inertia of the rigid body.
2.2. Large eddy simulation
ACCEPTED MANUSCRIPT In LES, the large-scale unsteady turbulent motions are resolved, while the effects of the smaller scale motions are modelled [35]. In the previous studies [22, 33], the well-known Smagorinsky model [36] was usually used:
T
t CS2 2 2Sij Sij
where xyz
, Sij 1 2 ui x j u j xi and CS is the Smagorinsky
SC R
13
(6)
IP
1/ 2
constant. However, near the walls or the quiescent flow regions, the subgrid-scale stresses should vanish, which is not automatically guaranteed in the standard
NU
Smagorinsky model [33]. In this study, we adopted a Lagrangian dynamic subgrid-scale model [37] in which CS is defined as:
MA
CS
FLM FMM
(7)
where FLM and FMM can be calculated by [37]
TE
D
n 1 FLM x H Lij M ij
n 1 FMM x M ij M ij
n 1
n 1
x 1 FLMn x u n t
n x 1 FMM x u n t
(8) (9)
CE P
̂ ̂ represents filtering at ̅|𝑆𝑖𝑗 ̅ − 4|𝑆̅|𝑆𝑖𝑗 ̅ ] , (∙) where 𝐿𝑖𝑗 = 𝑢̅̂ ̅𝑗 − 𝑢̅𝑖 𝑢̅𝑗 , 𝑀𝑖𝑗 = 2∆2 [|𝑆 𝑖𝑢
AC
scale 2 and
t T n , 1 t T n
n n T n 1.5 FLM FMM
1/8
where H x is the ramp function ( H x x if x 0, and zero otherwise). The turbulent kinetic energy k is for both resolved and unresolved scales, but Hartmann et al. [38] pointed out that the subgrid-scale kinetic energy in a midway baffle plane is about 0.1% of that at the grid-scale near the impeller and less in the rest domains. As a result, only the resolved kinetic energy is considered in the present study, i.e. k 1 2 u . 2
2.3. Determination of the source term f According to the direct forcing method, the boundary body force f in Eq. (1) can be simply calculated as [26]
(10)
ACCEPTED MANUSCRIPT ud un u u + p 2 u, on f f f f t 0, elsewhere
(11)
where t is the time step, ud is the desired velocity at an arbitrary point in the solid
IP
T
object, un is the velocity computed at step n, and is the domain in the vicinity of the IB. In general, the surface of the solid region where un+1= ud does not coincide with the
SC R
Cartesian grid points, especially when a staggered grid is applied, an interpolation procedure is thus needed.
NU
As mentioned in the introduction, a hybrid technique for determining the source term f is proposed in this paper. That is, the scheme of Uhlmann [16] is used for moving
MA
boundaries and that of Kajishima et al. [28] is used for stationary boundaries. The details are described below.
2.4. IB method for moving boundaries (moving approach)
D
Uhlmann [16] proposed a direct forcing method where the IBs are represented by a
TE
set of Lagrangian markers. This method allows a smooth transfer between the Eulerian
CE P
and Lagrangian representations while at the same time avoiding strong restrictions of the time step [16]. The force term is evaluated at the Lagrangian markers to suppress the oscillations of the hydrodynamic force due to insufficient smoothing when dealing
AC
with moving boundaries; i.e., F f
Ud U t
(12)
where the uppercase letters denote the variables at the Lagrangian points and U d is the desired velocity at the location on the fluid–solid interface and can be obtained from Eq. (3) The intermediate velocity U can be evaluated by interpolation from its Eulerian counterpart u , U u h x X h3 x
where h is the spacing of the Eulerian grid points and u can be computed from Eq. (1) without body force f. h is the discrete delta function usually constructed in the form [39]
(13)
ACCEPTED MANUSCRIPT h x X
1 x X h3 h
y Y z Z h h
(14)
where x, y, z and X, Y, Z are the components of x and X, respectively, and is the
T
one-dimensional discrete delta function:
if 0 r 1 if 1 r 2
SC R
IP
1 2 8 3 2 r 1 4 r 4r 1 φ r 5 2 r 7 12 r 4r 2 8 0
(15)
if 2 r
NU
To obtain the forcing term in Eq. (1), the body force acting on Lagrangian markers should be distributed to the Cartesian grid points: Nm
MA
f F h x X Vl
(16)
l 1
where Nm is the total number of markers and Vl is the volume assigned to the lth
D
Lagrangian marker; the definition refers to Uhlmann [16].
TE
2.5. IB method for stationary boundaries (stationary approach) Kajishima et al. [28] proposed an extremely simple scheme for modeling fluid–solid
CE P
interactions. At the interface, the fluid and solid velocities are smoothly connected using the solid volume fraction at each Cartesian grid point as a weight factor: (17)
AC
u 1 u ud
where represents the solid volume fraction at the computational grid point, and u and ud are the velocities within the fluid and solid domains, respectively. For the stationary tank wall, ud 0 , and the evaluation of the desired velocity at a Cartesian grid point can then be defined as
u 1 u
(18)
Similar to the case for the direct forcing method based on the no-slip boundary condition, the body force term at the boundary grid point should be
f f ud u t f u t
(19)
ACCEPTED MANUSCRIPT As the treated boundary is stationary, the volume fraction at the Cartesian grid points will not change during the computation. This means the volume fraction needs to be evaluated only once, usually before the simulation by a separate code.
IP
T
Obviously, accurate determination of the solid volume fraction is key to this method. The discrete approximation approach [40] is employed in this study. The grid
SC R
points lying sufficiently close to the interior of the IB have = 1 directly. Then, for grid points intersecting with the IBs, a series of marker points are generated uniformly
NU
in each of the three directions, with an interval much smaller than the interval of Cartesian grid points h (e.g., h/M, where M is an integer), and hence, the total number
MA
of marker points within a grid cell is M3. A ray tracing technique [41] is then employed to count the number of markers within the solid region. If this number is n, we immediately have = n/M3. Larger M gives higher accuracy but requires more
D
computational time. However, as the computation is only carried out once, large M is
TE
affordable in obtaining a fairly accurate estimation of .
CE P
Finally, in the context of the hybrid IB technique in this work, the body force at a Cartesian grid point can be represented as a sum of the two components for moving
AC
and stationary boundaries:
Nm
f F h x X Vl f u t l 1
Naturally, the two terms are nonzero only at different times and positions except for some short and localized contact.
2.5 Parallel implementation The hybrid IB technique described above is implemented and parallelized on a supercomputer with both multi-core CPUs and many-core GPUs. The Navier–Strokes equations for the fluid flow are solved on a fixed Cartesian grid, which can be parallelized straightforwardly. In the parallel implementation, the regular rectangle computational domain is divided into a number of sub-domains, and each process (running on a CPU core and the corresponding GPU) is assigned to a sub-domain. The
(20)
ACCEPTED MANUSCRIPT data exchange between sub-domains is communicated through inter-process data transfer in the industry-standard MPI protocol [42].
T
However, the IBs crossing subdomain borders may pose challenges, which will be
IP
tackled in this section. As the solid volume fraction of each grid point is stored locally, no additional difficulty is associated with the parallelization of the stationary approach.
SC R
Meanwhile, for the moving approach, the moving Lagrangian markers at the borders between subdomains cause a variety of complications. The interpolation of the velocity
NU
on markers (Eq. (13)) requires the exchange of extra flow information at the subdomain borders. In general, to minimize the communication, tagging and searching
MA
of Lagrangian points and construction of the delta functions should be done efficiently
AC
CE P
TE
D
[23] .
Figure 1. Schematic representation of the interpolation and distribution operations at the subdomain borders: (a) two layers of ghost Cartesian grid points, (b) one layer of ghost grid points.
The two steps for processing the IB markers, force distribution and velocity interpolation, both loop over the same indexes determined by the support of the regularized delta function, which has a width of four grid points in this study. Taking
ACCEPTED MANUSCRIPT the two-dimensional case in Figure 1(a) for example, we see that the IB maker m interacts with the surrounding 4×4 Eulerian grid points (point 1 to point 16).
T
Considering the staggered grid used, in the extreme case shown in Figure 1(a) when the
IP
Lagrangian markers are very close to the subdomain border, Eulerian grid points 1–8 are ghost grid points whose values should be transferred from the neighboring process.
SC R
Therefore, in general, two layers of ghost grid points are needed to perform the interpolation and distribution steps while only one layer is needed if a second order
during the IB stage in this extreme case.
NU
discretization scheme is applied. Therefore, two additional communications are needed
MA
■ Interpolation step: the data of the ghost Cartesian grid points (grid points 1–8) should be transferred from the adjacent process to compute the velocities of the Lagrangian
the case of the fluid solver.
D
markers (Eq. (13) ). One extra layer of ghost grid points is needed when compared with
TE
■ Distribution step: the body force source of Lagrangian markers should be spread to
CE P
the surrounding connecting Cartesian grid points (grid points 1–16). The velocity of the Cartesian grid points should be updated by adding the source terms to them. As a result, the calculated body force source term of Cartesian grid points 1–8 should be
AC
further transferred back to the neighboring process and added to the corresponding grid points.
Obviously,
these
two
additional
communications
greatly complicate
the
parallelization and increase the communication tremendously. To simplify the above procedure, a partial velocity interpolation and force distribution method is developed in this study, where the Lagrangian markers interact with only some of the 4×4 Eulerian grid points at the subdomain border. In this strategy, only one layer of ghost grid points is required and the interpolation and distribution steps are modified accordingly. It should be noted that only the Lagrangian markers assigned to the present process are computed. In Figure 1(b), for example, only Eulerian grid points 5–16 are used for the velocity interpolation and the body
ACCEPTED MANUSCRIPT force of m is only distributed to grid points 9–16. As a result, overlapping of one ghost grid point is sufficient and the body force does not need to be transferred to the
T
neighboring process any more. Two factors should be incorporated into the
IP
interpolation and distribution steps to ensure that the velocities are interpolated exactly
Cartesian grid points [43, 44]: 1 fv
u x X h
xv
Nm
f l 1
where f v
and Fl are defined by
x X h
xv
h
3
l
x X h
x f
h
3
l
D
Fl
3
l
F h x X l Vl Fl
MA
fv
h
NU
Ul
SC R
and the body forces exerted by the Lagrangian markers are fully transferred to the
TE
Here v is the Cartesian domain for the interpolation of the intermediate velocity at the Lagrangian markers (in the extreme case of Figure 1(b), it includes grid points 5–
CE P
16), and f is the domain for the body force spreading (in the extreme case of Figure 1(b), it includes grid points 9–16). Naturally, when a Lagrangian marker lies more than two grid point widths from the subdomain border, v and f will consist of all 16
AC
surrounding Eulerian grid points, and f v 1 and Fl 1 .Applying this strategy, the communication becomes more efficient. This strategy is derived from and an improvement of the half-distributed forcing strategy for IB methods, which was demonstrated to be feasible in our previous publication [44]. Furthermore, the forces and torques added to the fluid through the distribution procedures (Eq. (16)) are not changed [44]. From the discussions above, we conclude that the proposed partial interpolation and distribution strategy will not lead to additional communications. As a result, the parallelization is much simpler and more efficient.
3. Validation
(21) (22)
(23) (24)
ACCEPTED MANUSCRIPT To evaluate the efficiency of the hybrid IB technique and the accuracy of the proposed parallelization strategy, we simulated the flow field in two stirred tanks
T
with this method. The main results are reported in the follows.
IP
3.1.A simple stirred tank
The method is first applied to a cylindrical unbaffled stirred tank with eight straight
SC R
blades at the mid-height as an impeller. This tank has been studied experimentally [45] and numerically [2]. A schematic illustration of the stirred tank is presented in Figure 2, where the height of the tank H=10 cm. The impeller rotates at n=100 rpm, which
NU
leads to an impeller Reynolds number Re ND2 1024 , in terms of the turbine
MA
diameter D, turbine frequency N and fluid kinematic viscosity . Direct numerical simulation is possible for the present simple stirred tank in which geometrical dimensions and rotation speed are both small. The rectangular computational domain is
D
discretized using a uniform Cartesian grid of 256×256×256 grid points, which is finer
TE
than that used in the DNS by Verzicco et al. [2]. The computational domain is divided
CE P
into 2×2×2 subdomains in the parallel computing. No-slip boundary condition is imposed at the bottom of the tank, while free-slip condition is assigned at the top. Within the computational domain, no-slip boundary conditions are enforced for the
AC
cylindrical tank wall and the impeller using the hybrid IB technique described above.
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 2. Schematic illustration of the stirred tank (adapted from Dong et al. [45] and
D
Verzicco et al. [2] ).
TE
We compare the computational efficiencies of the two IB methods described above in handling the stationary and moving IBs, as listed in Table 1. Note that, if a moving
CE P
boundary is handled by applying the stationary approach (Case 1), the grid volume fraction should be evaluated in each time step because the position of the boundary changes with time. In a stirred tank, the surface area of the rotating impeller is
AC
generally small compared with the stationary wall, and the total number of Lagrangian markers is limited. Conversely, if the tank wall is also simulated with the moving approach (Case 2), the number of Lagrangian markers will increase by approximately an order of magnitude (from 15,611 to 132,183 for the present stirred tank), and the computational time will increase sharply. From Table 1, we see that the time required for the hybrid IB technique is much less than the traditional approach to a single IB for both stationary and moving boundaries. Table 1. Comparison of the efficiency of the two IB methods in the case of a stirred tank. SA and MA denote the stationary approach and moving approach, respectively. Data are the time taken for 100 time steps. In Cases 1 and 3, a small value of M (= 3) is used when evaluating the grid solid fraction in the stationary approach.
Case
Case 1
Case 2
Case 3
ACCEPTED MANUSCRIPT Impeller
Wall
Impeller
Wall
Impeller
Wall
SA
SA
MA
MA
MA
SA
Time (s)
889.646
193.708
92.144
T
Furthermore, the averaged velocity components are compared with experimental
IP
data [45], as shown in Figure 3. Data are evaluated at two different axial locations z/H
SC R
= 0.5 and 0.7. A satisfactory agreement is obtained for all three velocity components in these three cases, except for small deviations in the axial velocity at z/H = 0.7 and in the radial velocity at z/H = 0.5 in Case 1. This may be ascribed to the approximation of
AC
CE P
TE
D
MA
NU
the solid volume fraction due to the small M used.
Figure 3. Comparison of the computed (lines) average velocity components with the experimental data from Dong et al. [45] (symbols ▲) along the radial direction at different axial locations z/H = 0.5 (a, b, c) and 0.7 (d, e, f). ua, ur and ut represent axial, radial and tangential velocity components, respectively.
ACCEPTED MANUSCRIPT
3.2. Turbulent flow in a Rushton stirred tank
T
To verify our method in more practical applications, we then simulated a tank stirred
IP
by a Rushton turbine, and LES is incorporated to describe the turbulent flow in it [21,
SC R
22, 33, 46]. The problem has been studied by Hartmann et al. [38] who presented the comparison of LES and RANS simulation with LDA measurements. The stirred tank is encapsulated by a standard cylindrical vessel 150 mm in diameter
NU
(T = 150 mm) with four equispaced baffles. The geometry is sketched in Figure 4. The density of the working fluid (𝜌) is 1039 kg ∙ 𝑚−3 and the dynamic viscosity (𝜇) is
MA
15.9mPa ∙ s. The rotating speed is 2672 rpm, resulting in an impeller Reynolds number Re = 7300. This Reynolds number is not within the fully turbulent regime, however, the flow is fully turbulent near the turbine [29]. A uniform Cartesian grid of 5123 is
D
used in the computation and divided into 43 subdomains. Starting from the quiescent
TE
condition, the stirred tank flow requires approximately 100 turbine revolutions to reach
CE P
a statically steady state [29]. The start-up was first computed on the 2563 grid for about 100 revolutions and then interpolated onto the 5123 grid for another 15 revolutions. Then, the flow variables were collected from the last 10 revolutions. In the
AC
computation, the timestep was selected according to the tip speed U tip and the grid width h , i.e., tUtip h 0.2 . The computation requires approximately 1122 timesteps for a full revolution on the grid of 2563 and 2244 timesteps on the grid of 5123.
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 4. Cross-section of the tank driven by a Rushton turbine [33].
D
To compare the present LES results with experimental and numerical data from
TE
other publications, the phase averaged variables are evaluated. Figure 5 shows the
CE P
experimental, RANS and LES results of the phase averaged flow field and their associated kinetic energy. Figure 5(d) is the present LES results, Figure 5(c) is the LES results using a standard Smagorinsky subgrid-scale model with constant CS 0.1 [33].
AC
The main features of the flow field are the two large circulation loops on both the top and bottom of the impeller. The working fluid is pumped radially towards the tank wall, partly goes upward above the impeller and partly goes downward below the impeller. Finally, the fluid flows into the impeller region along the axis. In this figure, the RANS and LES for both subgrid-scale models show good agreement to the experiment. The magnitudes of the kinetic energy of the simulation results are comparable with that of the experiment. As can be seen, the kinetic energy in the flow jet area is much higher than the rest of the tank. A quantitative comparison of phase-averaged radial and tangential velocity components at three radial locations is shown in Figure 6. A satisfactory agreement is obtained between the simulation and the experiment except a small deviation for radial
ACCEPTED MANUSCRIPT velocity at 𝑟⁄𝑇 = 0.317 and the tangential velocity at location 𝑟⁄𝑇 = 0.183. All the simulation results overestimate the tangential velocity at location 𝑟⁄𝑇 = 0.183, but
T
the present LES result is better than that of RANS and LES(s). From the comparison,
IP
we can see that the results of RANS are obviously worse than those from LES, especially for the tangential velocity at 𝑟⁄𝑇 = 0.25 and 𝑟⁄𝑇 = 0.317. LES with
SC R
different subgrid-scale models get similar results. Although the advantage is not obvious, the results from the present LES is a bit better in general. For the present LES,
NU
the radial velocity at 𝑟⁄𝑇 = 0.183 and the tangential velocity at 𝑟⁄𝑇 = 0.317 are almost identical to the experimental data, while bigger deviation is seen for the LES
MA
with standard subgrid-scale model (LES(s)).
Figure 7 shows the axial profile of the turbulent kinetic energy at three different radial locations. The deviation of the kinetic energy between the simulation and the
D
experiment is much larger than that of velocity. Underestimation is observed for all the
TE
simulations at the impeller tip (i.e. 𝑟⁄𝑇 = 0.183). RANS predictions significantly
CE P
underestimate the kinetic energy at all the three locations, mesh refinements are not expected to improve substantially the prediction [10]. LES with the standard subgrid-scale model overestimate at locations away from the impeller region (i.e.
AC
𝑟⁄𝑇 = 0.25 and 𝑟⁄𝑇 = 0.317). While the peak values of the kinetic energy predicted by the present LES are comparable with that of the experiment at locations away from the impeller region (i.e. 𝑟⁄𝑇 = 0.25 and 𝑟⁄𝑇 = 0.317), but the width of the profile is smaller. From this study, we can find that the predictions of the present LES-IB method compare satisfactorily in general against the experimental data, and it is at least as accurate as other LES-IB methods. It suggests that this method can be promising for the simulation of stirred tanks at large scale with complex geometries and at high Reynolds numbers.
AC
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 5. Phase averaged plot of the velocity vector field and turbulent kinetic energy. (a)
Experiment, phase averaged turbulent kinetic energy are only measured in the impeller region, (b) RANS, (c) LES with Smagorinsky model, (d) Present LES model. (a)-(c) are taken from Hartmann et al. [38].
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 6. Phase averaged axial profile of the radial (left) and tangential (right) velocity
AC
components at three different radial locations r/T = 0.183 (top), r/T = 0.25 (middle) and r/T = 0.317 (right).
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 7. Phase averaged axial profile of the kinetic energy of the random velocity fluctuations
TE
D
at three different radial locations r/T = 0.183, r/T = 0.25 and r/T = 0.317.
4. Application
of
the
method
to
particle-resolved
CE P
simulation of stirred tanks Solid suspensions in liquid are widely encountered in technical systems. There is
AC
no universal approach for modelling and simulation of such problems, the appropriate method depends on the flow regime and on the required details and accuracy [47], but limited by the computational cost, Eulerian-Eulerian (two fluid model, TFM) and Eulerian-Lagrangian (discrete particle model, DPM) methods are used most often [48, 49]. However, empirical correlations are needed to account for unresolved parts of the suspensions physics in these two models and we cannot obtain the complex flow field around the solid particles. Discrete element method (DEM) , originally developed by Cundall and Strack [50], is a soft particle model which allows overlap between the contacting particles [50, 51]. A linear spring-dashpot force model is used in the implementation. In this study, IBM is coupled with the particle simulator to handle the solids suspension inside stirred tanks. The motion of particles and the particle-particle collisions are treated by DEM
ACCEPTED MANUSCRIPT [51], while no-slip conditions at the surfaces of the particles are assigned by using IBM. This method is usually called particle-resolved simulation [49], where the
T
Eulerian grid width is approximately smaller by an order of magnitude. In this method,
IP
no correlations are required, both the particle-particle and particle-fluid interactions are modelled in a realistic way. There is obviously a computational penalty for resolving
SC R
down to finer scales and the physical size of the domain that can be simulated gets limited. As a result, smaller stirred tanks compared to laboratory scale are studied here.
NU
It should be mentioned that the lubrication force is also considered in the collisions to account for condition when particles get in very close proximity [52].
MA
During the collision, a thin lubrication layer is formed between the particles and the fluid is squeezed out of this gap when the particles approach and is pushed back into the gap during rebound. Viscous forces hence become important in the approaching
D
and rebounding process and can lead to sizeable dissipation. This feature might be
TE
resolved by an adaptive local grid refinement [53], but this usually results in
CE P
substantially increased computation time. An appropriate kind of “subgrid-scale model” accounting for the film between the surfaces should be implemented when it is too thin to be resolved by the Eulerian grid. This contribution can be quantitatively
AC
described by [49, 54]
Flub
Ri2 R 2j 1 1 n uij s s1 6 f f 2 Ri R j s1 s0 Ri2 R 2j 1 1 6 f f n uij s1 s s0 2 Ri R j s s0 s s0 0
where n is the unit vector from particle i to particle j, uij ui u j , s0 0.2R0 and s1 2 104 R0 .
The physical parameters in the computation are identical to the experiments [55] and the tank size is smaller than that in experiment in consideration of the computation cost. A schematic diagram of the stirred tank is shown in Figure 8, and the numerical and physical parameters are shown in Table 2.
(25)
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 8. Schematic illustration of the lab-scale stirred tank. left: side view, right: top view.
MA
𝑫 = 𝟐𝟒 𝐦𝐦，𝐛 = 𝟗. 𝟏𝟕 𝐦𝐦.
Table 2．The numerical and physical parameters in the simulation of solid suspension in a
Value
Tank size
483
TE
Parameter
D
lab-scale stirred tank.
mm
Parameter
Value
Reference length
1
mm -6
0.75
mm
Reference mass
1.15×10
Particle radius (impeller)
0.5
mm
Reference velocity
0.1
ms-1
Particle density
2485
kgm-3
Reynolds number (Reimp)
15850
--
1150
-3
Eulerian grids (∆h)
0.09375
mm
CE P
Particle radius (freely)
Fluid density Gravity Viscosity
4.389e
AC
Particle mass
kgm -6
kg -2
9.81
0.87×10
ms -6
m s
2 -1 -1
-5
kg
∆t
1×10
Markers (per particle)
1008
--
Markers (impeller)
69141
--
s
Rotating velocity
50.14
rads
GPUs used
64
--
Particle numbers (freely)
4731
--
Particle numbers (impeller)
2790
--
Three cases with different geometric configurations are carried out in this study, the setup are listed in Table 3. The average solid volume fractions in these three cases are the same, 𝜙̅ = 10.4%. The no-slip boundary condition is enforced on all surfaces including the impeller, the tank wall, the baffles and the particles using the hybrid method of this study. It is difficult to judge the relative location between the particles and the impeller, for the convenience to deal with particle-impeller collisions, the impeller is represented by the combination of particles with diameter of 𝑑 = 1 mm.
ACCEPTED MANUSCRIPT As a result, there are two sets of representations of impeller: combination of particles during particle-impeller collision; and Lagrangian markers during fluid-solid
T
interactions. The particles are initially located at both the upside and the underneath of
IP
the impeller, as shown in Figure 9.
Table 3. The setup of the three cases in the simulation of solid suspension in a stirred tank.
Baffle
PBTD
1
Yes
Case 2
PBTU
2
Yes
Case 3
PBTU
Case 1
No
1
PBTU: up-pumping pitched blade turbine
15850
NU
PBTD: down-pumping pitched blade turbine
2
Reimp
SC R
Impeller
MA
The computational grid in this simulation is 512×512×512, the resolution of the present grid is not fine enough for DNS, as a result LES with wall-layer model is
D
applied in the simulations. A resolution of 16 grid spacing over a particle diameter is
TE
used to resolve the flow around a solid sphere. CFL is around 0.2 and it takes about 4177 timesteps for a full revolution.
CE P
We investigated the startup of the particle suspension in the three cases. As shown in Figure 9, we can conclude from the qualitative analysis that PBTU is more conducive
AC
to particle suspension than PBTD, and more uniform particle distribution in the radial direction is obtained when baffles are used. This difference is shown more clearly in Figure 10 with the distribution of solids concentration, and more quantitatively in Figure 11 using a global uniformity index (𝜉) [56]
1 2 1 1 NC
1 2
i 1 i 1 NC
2 i 2
where 𝜙𝑖 is the solid fraction in ith grid. The index is conveniently defined so that as 𝜉 → 0 the uniformity of the suspension is at its minimum; and when 𝜉 → 1 the solids are uniformly distributed within the vessel. Figure 12 shows the phase averaged velocity of the liquid phase, solid phase and the solid velocity between the two. The averaged solid velocity is calculated by
(26)
ACCEPTED MANUSCRIPT 𝑢̅𝑠𝑜𝑙𝑖𝑑 = (∑ 𝜙 ∙ 𝑢𝑠𝑜𝑙𝑖𝑑 )⁄∑ 𝜙 . Using cylindrical coordinates, the radial-axial 2D velocity maps were obtained by azimuthally averaging the 3D velocity data and
T
projecting them onto the 2D radial-axial plane. Whilst being two-dimensional, these
IP
flow patterns for both the fluid and the solid are qualitatively consistent with the experimental results of Guida et al. [56], which can be considered as a verification of
AC
CE P
TE
D
MA
NU
SC R
the present results.
Figure 9. Evolution of the solids suspension process: instantaneous realizations of particles locations at four different moments as indicated for the three cases.
AC
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 10. Time-averaged solids volume fraction in the vertical middle plane (left) and in the horizontal plane with z/H = 0.1 (right).
SC R
IP
T
ACCEPTED MANUSCRIPT
AC
CE P
TE
D
MA
NU
Figure 11. Time variation of global uniformity index in the three cases.
Figure 12. Phase averaged velocity of the liquid phase, solid phase and the difference between the two.
ACCEPTED MANUSCRIPT 5. Conclusions An efficient method for the DNS/LES of complex flow in stirred tanks on fixed
T
Cartesian grids is developed in the framework of the finite difference method. The
IP
moving and stationary boundaries are handled respectively by the direct forcing method using Lagrangian markers and a method based on solid volume fractions. This
SC R
hybrid technique is more efficient than using a single method applied to both stationary and moving boundaries.
NU
Meanwhile, the method is further implemented on a supercomputer with CPU–GPU hybrid architecture. The existence of the IB complicates the parallelization. To
MA
minimize the amount of information that needs to be communicated through the inter-subdomains, an easy partial velocity interpolation and force distribution strategy is developed for the parallelization of the IB, which greatly reduces the communication
D
and simplifies the implementation procedure.
TE
The accuracy and robustness of the present method were validated in two relatively
CE P
simple simulations: fluid flow in a simple stirred tank and turbulent flow in a Rushton tank. Both validation cases demonstrated, respectively, the efficiency of the hybrid IB technique and the capability of the combined LES and IBM implementation.
AC
Furthermore, hybrid IB technique is coupled with DEM to study the particle suspension inside the stirred tanks. In this method, no correlations are required; both the particle-particle and particle-fluid interactions are modelled in a realistic way. As a result, this method can be applied to study the fundamental mechanism of solid suspensions in complex flows. In summary, the present method can be used as a promising tool for simulating turbulent flows in stirred tanks with complex geometries.
ACCEPTED MANUSCRIPT NOMENCLATURE Smagorinsky constant
D
turbine diameter, m
f
additional force term, Nm-3
fv
factor
Fl
factor
F
Lagrangian point force, Nm-3
Fs
hydrodynamic force, N
g
gravity force, ms-2
Gs
external force, N
h
Eulerian grid width, m
H
tank height, m
Is
inertia momentum, kg·m2
k
turbulent kinetic energy , m2s-2
ms
mass, kg
L
length, m
M,n
integer
n
unit vector between two particles, (m,m,m)
N
turbine frequency , s-1
Nm
total number of markers
Ns
torque, N·m
p
pressure, Nm-2
r
radial location, m
Re ̅ 𝑆𝑖𝑗
Reynolds number
t
time, s
u ̃ 𝒖
IP SC R NU
MA
D
TE
CE P
Rate of strain tensor, s-1
AC
Ts
T
Cs
torque, N·m velocity, ms-1 intermediate Eulerian velocity, ms-1
ud
desired velocity , ms-1
un
velocity at nth time step, ms-1
ut
translational velocity, ms-1
U ̃ 𝑼
velocity at Lagrangian marker, ms-1
Ud
desired velocity at the fluid-solid interface, ms-1
Utup
impeller tip speed, ms-1
vs
linear velocity , ms-1
x
Eulerian position vector, m
x0
center coordinate of a rigid body, m
X
Lagrangian position vector, m
z
axial location, m
intermediate Lagrangian velocity, ms-1
ACCEPTED MANUSCRIPT
Greek letters solid volume fraction density of fluid, kgm-3
f
dynamic viscosity, Nsm-2
t
kinematic viscosity, m2s-1
a time ratio in LES model
global uniformity index
rotational velocity, rads-1
s
rotational velocity of the rigid body, rads-1
domain
t
time step, s
Vl
volume of the lth Lagrangian marker, m3
h
Dirac delta function
Operators
partial time derivative, s-1
gradient operator, m-1
divergence operator, m-1
2
Laplace operator, m-1
filter width, m
CE P
TE
D
t
AC
IP SC R
NU
MA
one-dimensional delta function
T
f
ACCEPTED MANUSCRIPT
AC
CE P
TE
D
MA
NU
SC R
IP
T
REFERENCES
AC
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
AC
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
AC
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
ACCEPTED MANUSCRIPT 2015-0417
AC
CE P
TE
D
MA
NU
SC R
IP
T
Graphic abstract
S1004-9541(16)30507-9 doi: 10.1016/j.cjche.2016.05.031 CJCHE 584
To appear in: Received date: Revised date: Accepted date:
29 August 2015 9 March 2016 18 March 2016
Please cite this article as: Shengbin Di, Ji Xu, Qi Chang, Wei Ge, Numerical simulation of stirred tanks using a hybrid immersed-boundary method, (2016), doi: 10.1016/j.cjche.2016.05.031
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
IP
T
Numerical simulation of stirred tanks using a hybrid immersed-boundary method*
1
SC R
Shengbin Di(狄升斌)1,2, Ji Xu(徐骥)1, Qi Chang(常麒) 1,3，Wei Ge(葛蔚)1,** State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese
Academy of Sciences, Beijing 100190, China 2
NU
University of Chinese Academy of Sciences, Beijing 100049, China
3
MA
School of Chemical Engineering and Technology, Tianjin University, Tianjin 300072, China
Article history: Received 29 August 2015
CE P AC
TE
Accepted 18 March 2016
D
Received in revised form 9 March 2016
ABSTRACT Conventionally, multiple reference frame (MRF) method and sliding mesh (SM) method are used in the simulation of stirred tanks, however, both methods have limitations. In this study, a hybrid immersed-boundary (IB) technique is developed in a finite difference context for the numerical simulation of stirred tanks. IBs based on Lagrangian markers and solid volume fractions are used for moving and stationary boundaries, respectively, to achieve optimal efficiency and accuracy. To cope with the high computational cost in the simulation of stirred tanks, the technique is implemented on computers with hybrid architecture where central processing units (CPUs) and graphics processing units (GPUs) are used together. The accuracy and efficiency of the present technique are first demonstrated in a relatively simple case, and then the technique is applied to the simulation of turbulent flow in a Rushton stirred tank with large eddy simulation (LES). Finally the
ACCEPTED MANUSCRIPT proposed methodology is coupled with discrete element method (DEM) to accomplish particle-resolved simulation of solid suspensions in small stirred tanks. It demonstrates that the
T
proposed methodology is a promising tool in simulating turbulent flow in stirred tanks with complex
IP
geometries.
SC R
Keywords: immersed-boundary method; CPU-GPU hybrid computing; stirred tank; large eddy
AC
CE P
TE
D
MA
NU
simulation
ACCEPTED MANUSCRIPT 1. Introduction Stirred tanks are widely used in chemical and process industries, and hence deserve experimental
and
computational
investigation.
Experimental
T
comprehensive
IP
measurement techniques such as laser-doppler velocimetry (LDV) and particle image
SC R
velocimetry (PIV), have significantly contributed to our understanding of the flow structure in stirred tanks. However, most experimental studies have limitations, such as time and cost, scale-up difficulties, and the difficulty of directly measuring some
NU
quantities like local velocity especially for high-viscosity fluids[1]. Computational fluid dynamics provides an alternative tool for unveiling the details of fluid flow and
MA
heat transfer in stirred tanks, but faces several challenges. The main difficulties lie in the relative motion between the impeller and vessel and the complex geometry of the
D
impeller [2]. To address these difficulties, several numerical methods have been
TE
developed.
CE P
Conventionally, boundary-conforming nonorthogonal grids are used to handle the complex geometries [3]. In the so-called body-fitted approaches [4], the grid lines follow the boundaries exactly and the boundary conditions can therefore be
AC
implemented directly. However, the generation of such grids is both difficult and time consuming, and it has to be repeated or modified during the simulation for moving boundaries. Furthermore, the discretized governing equations for body-fitted grids are much more complex than those for Cartesian grids, which increase the computational cost and decrease the numerical stability considerably. As a result, such a method is seldom directly applied to the simulation of stirred tanks. A non-inertial reference frame fixed to the moving body can be used to avoid mesh regeneration for moving boundaries, although the governing equations have to take a more complicated form. Multiple reference frames (MRF) [5] and sliding mesh (SM) models [6] are typical examples of this approach. They are among the most commonly used methods for the simulation of stirred tanks and available in mainstream commercial codes [7-10]. Both
ACCEPTED MANUSCRIPT methods solve the governing equations in two domains: an inner domain containing the rotating impeller and an outer domain containing the stationary tank wall and baffles.
T
In SM, the meshes are allowed to shear and slide to accommodate the relative motion
IP
at the domain interface, while in MRF, only a stationary solution of the flow field is provided for a given position of impeller with respect to the baffles [11]. The main
SC R
disadvantages of SM are the explicit acceleration terms added to the momentum equation and the related numerical problem at the domain interface. For complex
computationally very demanding [12].
NU
impeller designs and intermeshing extruders, both MRF and SM might be infeasible or
MA
Immersed boundary (IB) methods, first proposed by Peskin [13], are another category of methods for complex geometries that can be based on simple and regular Cartesian grids, and have received much attention recently [14-16]. IB methods have
D
higher computational efficiency when compared with body-fitted methods [16-19]. In
TE
IB methods, the Cartesian grids do not necessarily conform to the boundaries. The
CE P
effects of the boundaries on the fluid are represented by additional body forces in the momentum equations. Although the number of grid points increases more quickly than in body-fitted methods as the Reynolds number increases, the per-grid-point-operation
AC
count is notably reduced owing to the absence of grid transformation [20]. In recent years, IB methods have been applied in the simulation of complex fluid flow in stirred tanks because of their efficiency and capacity to deal with moving boundaries directly [2, 21, 22]. Unlike SM model, there is no constraint on the prescribed motion of the rigid bodies in IB methods [23]. Among various IB methods [15, 16, 24-26], two typical approaches are the bases of this study. In the first approach, the boundaries are represented by a set of Lagrangian points which are handled with direct forcing methods such as that proposed by Uhlmann [16]. No-slip boundary conditions are assigned to these Lagrangian markers. The positions of the Lagrangian markers are updated in each time step if the boundaries are moving. As only the surfaces of the immersed bodies should be
ACCEPTED MANUSCRIPT represented by the markers, the number of markers is relatively small compared with the total number of Cartesian grid points. Each Lagrangian point is connected with a
T
number of Cartesian grids around it. Mapping of the variables between Lagrangian and
IP
Eulerian locations is performed employing a regularized delta function. This technique effectively reduces non-physical force oscillations existing in moving-boundary
SC R
simulations[27], however, it is not necessary for stationary-boundary problems. Furthermore, when the IBs lie close to the borders of the computational domain, where
NU
the separation is smaller than the support width of the delta function, it is difficult to find enough Cartesian grid points for velocity interpolation and force distribution when
MA
using this technique.
In the second approach, the body force is evaluated directly at the IB using the desired velocity of the solid part. However, the fluid and solid velocities are smoothly
D
interpolated using a solid volume fraction at the boundary grid points. When the solid
TE
part is stationary, the solid fraction at the computational grid points does not change
CE P
during the computation. As a result, the evaluation of the solid fraction needs to be carried out only once. The method of Kajishima et al. [28] belongs to this category. It is noted that the two approaches described above are more suitable for moving and
AC
stationary boundaries, respectively (For simplicity, they are referred as the moving and stationary approaches hereafter). Both moving and stationary boundaries, such as the boundaries of rotating impellers and stationary tank walls, are present in a stirred tank simultaneously. However, most simulations use only one approach for all boundaries [22, 29, 30]. Tyagi et al. [30] proposed a method to combine the advantages of the IB method to represent moving rigid geometries and the multi-block structured curvilinear mesh to represent complex domains. However, it is still a challenge to generate body-fitted mesh for complex geometries such as the helical coils inside a tank. To this regard, we will propose in this work a hybrid IB technique in which the complex geometries are dealt with using different IB methods on a Cartesian grid.
ACCEPTED MANUSCRIPT On the other hand, one practical difficulty in simulating complex flows in stirred tanks is the high computation cost. Consequently, we implement the IB method on a
T
CPU-GPU hybrid architecture for massively parallel computing. Algorithms having
IP
good parallelism are therefore desirable. Since the fluid and IBs are discretized on two separate grids, their parallelization should be handled differently. The parallelization
SC R
for Eulerian grids is easy because the grid points are fixed in time and are assigned to the same subdomain for the whole computation period [31]. However, the Lagrangian
NU
markers move throughout the computational domain, so a given IB point may move across the subdomain border during one time step. Furthermore, the interactions
MA
between the fluid and IBs are another source of difficulty in the design of a distributed algorithm [32]. To cope with these problems, we developed a partial velocity interpolation and force distribution scheme, which is described in section 2..
D
To validate the proposed method, we first test its accuracy and robustness in two
TE
numerical experiments, incorporating large eddy simulation (LES) which is found to be
CE P
effective for turbulent flow in stirred tanks [33, 34]. Then, as an application, the discrete element method (DEM) is coupled with the present method to carry out direct numerical simulation (DNS) of fluid-solid systems in stirred tanks. Satisfactory results
AC
reveal that complex flows in stirred tanks can be accomplished using hybrid IB methods on supercomputers. It will become a powerful tool for optimal design of industrial stirred tanks if the method can be further developed. Finally, main conclusions are drawn with perspectives on wider applications of the present method.
2. Numerical method 2.1. Governing equations In the current work, we focus on incompressible fluid with constant properties. The governing equations of this method describe the conservation of momentum and mass of the flow. They are established for the fluid phase and the solid objects (or boundaries) separately, and coupled through the force terms describing the interactions between them.
ACCEPTED MANUSCRIPT 2.1.1.
Fluid-phase equations
The governing equations for unsteady incompressible flow in the IB methods take
u f u u p f 2 u + f t
(1) (2)
SC R
u 0
IP
f
T
the general form[16, 17]
where f and f are the fluid density and dynamic viscosity, respectively, while u, and p are the velocity and pressure, respectively. They are basically the Navier–Stokes
NU
equations with additional body force terms f determined by the no-slip boundary conditions at the fluid–solid interface. f is non-zero only in the vicinity of fluid–solid
2.1.2.
Solid-phase equations
MA
interfaces.
D
Only rigid bodies are considered in the present study. Given a prescribed motion, the
TE
desired velocity field within the solid region is [16]
ud ut x x0
(3)
CE P
where ut , and x0 are the translational and rotational velocity and center coordinates of a rigid body. Otherwise, if the motion of the rigid solid domain is due to
AC
the fluid force, the motion is controlled by Newton’s equations for linear and angular momenta of a rigid body, d ms vs dt d I s s dt
Fs Gs
(4)
Ts N s
(5)
where v s and s are the linear velocity at the mass center and the angular velocity of the rigid body, respectively. Fs and Ts represent the hydrodynamic force and torque acting on the rigid body due to the fluid. Gs and Ns are the external force and torque acting on the rigid body. ms is the mass of the rigid body and Is is the moment of inertia of the rigid body.
2.2. Large eddy simulation
ACCEPTED MANUSCRIPT In LES, the large-scale unsteady turbulent motions are resolved, while the effects of the smaller scale motions are modelled [35]. In the previous studies [22, 33], the well-known Smagorinsky model [36] was usually used:
T
t CS2 2 2Sij Sij
where xyz
, Sij 1 2 ui x j u j xi and CS is the Smagorinsky
SC R
13
(6)
IP
1/ 2
constant. However, near the walls or the quiescent flow regions, the subgrid-scale stresses should vanish, which is not automatically guaranteed in the standard
NU
Smagorinsky model [33]. In this study, we adopted a Lagrangian dynamic subgrid-scale model [37] in which CS is defined as:
MA
CS
FLM FMM
(7)
where FLM and FMM can be calculated by [37]
TE
D
n 1 FLM x H Lij M ij
n 1 FMM x M ij M ij
n 1
n 1
x 1 FLMn x u n t
n x 1 FMM x u n t
(8) (9)
CE P
̂ ̂ represents filtering at ̅|𝑆𝑖𝑗 ̅ − 4|𝑆̅|𝑆𝑖𝑗 ̅ ] , (∙) where 𝐿𝑖𝑗 = 𝑢̅̂ ̅𝑗 − 𝑢̅𝑖 𝑢̅𝑗 , 𝑀𝑖𝑗 = 2∆2 [|𝑆 𝑖𝑢
AC
scale 2 and
t T n , 1 t T n
n n T n 1.5 FLM FMM
1/8
where H x is the ramp function ( H x x if x 0, and zero otherwise). The turbulent kinetic energy k is for both resolved and unresolved scales, but Hartmann et al. [38] pointed out that the subgrid-scale kinetic energy in a midway baffle plane is about 0.1% of that at the grid-scale near the impeller and less in the rest domains. As a result, only the resolved kinetic energy is considered in the present study, i.e. k 1 2 u . 2
2.3. Determination of the source term f According to the direct forcing method, the boundary body force f in Eq. (1) can be simply calculated as [26]
(10)
ACCEPTED MANUSCRIPT ud un u u + p 2 u, on f f f f t 0, elsewhere
(11)
where t is the time step, ud is the desired velocity at an arbitrary point in the solid
IP
T
object, un is the velocity computed at step n, and is the domain in the vicinity of the IB. In general, the surface of the solid region where un+1= ud does not coincide with the
SC R
Cartesian grid points, especially when a staggered grid is applied, an interpolation procedure is thus needed.
NU
As mentioned in the introduction, a hybrid technique for determining the source term f is proposed in this paper. That is, the scheme of Uhlmann [16] is used for moving
MA
boundaries and that of Kajishima et al. [28] is used for stationary boundaries. The details are described below.
2.4. IB method for moving boundaries (moving approach)
D
Uhlmann [16] proposed a direct forcing method where the IBs are represented by a
TE
set of Lagrangian markers. This method allows a smooth transfer between the Eulerian
CE P
and Lagrangian representations while at the same time avoiding strong restrictions of the time step [16]. The force term is evaluated at the Lagrangian markers to suppress the oscillations of the hydrodynamic force due to insufficient smoothing when dealing
AC
with moving boundaries; i.e., F f
Ud U t
(12)
where the uppercase letters denote the variables at the Lagrangian points and U d is the desired velocity at the location on the fluid–solid interface and can be obtained from Eq. (3) The intermediate velocity U can be evaluated by interpolation from its Eulerian counterpart u , U u h x X h3 x
where h is the spacing of the Eulerian grid points and u can be computed from Eq. (1) without body force f. h is the discrete delta function usually constructed in the form [39]
(13)
ACCEPTED MANUSCRIPT h x X
1 x X h3 h
y Y z Z h h
(14)
where x, y, z and X, Y, Z are the components of x and X, respectively, and is the
T
one-dimensional discrete delta function:
if 0 r 1 if 1 r 2
SC R
IP
1 2 8 3 2 r 1 4 r 4r 1 φ r 5 2 r 7 12 r 4r 2 8 0
(15)
if 2 r
NU
To obtain the forcing term in Eq. (1), the body force acting on Lagrangian markers should be distributed to the Cartesian grid points: Nm
MA
f F h x X Vl
(16)
l 1
where Nm is the total number of markers and Vl is the volume assigned to the lth
D
Lagrangian marker; the definition refers to Uhlmann [16].
TE
2.5. IB method for stationary boundaries (stationary approach) Kajishima et al. [28] proposed an extremely simple scheme for modeling fluid–solid
CE P
interactions. At the interface, the fluid and solid velocities are smoothly connected using the solid volume fraction at each Cartesian grid point as a weight factor: (17)
AC
u 1 u ud
where represents the solid volume fraction at the computational grid point, and u and ud are the velocities within the fluid and solid domains, respectively. For the stationary tank wall, ud 0 , and the evaluation of the desired velocity at a Cartesian grid point can then be defined as
u 1 u
(18)
Similar to the case for the direct forcing method based on the no-slip boundary condition, the body force term at the boundary grid point should be
f f ud u t f u t
(19)
ACCEPTED MANUSCRIPT As the treated boundary is stationary, the volume fraction at the Cartesian grid points will not change during the computation. This means the volume fraction needs to be evaluated only once, usually before the simulation by a separate code.
IP
T
Obviously, accurate determination of the solid volume fraction is key to this method. The discrete approximation approach [40] is employed in this study. The grid
SC R
points lying sufficiently close to the interior of the IB have = 1 directly. Then, for grid points intersecting with the IBs, a series of marker points are generated uniformly
NU
in each of the three directions, with an interval much smaller than the interval of Cartesian grid points h (e.g., h/M, where M is an integer), and hence, the total number
MA
of marker points within a grid cell is M3. A ray tracing technique [41] is then employed to count the number of markers within the solid region. If this number is n, we immediately have = n/M3. Larger M gives higher accuracy but requires more
D
computational time. However, as the computation is only carried out once, large M is
TE
affordable in obtaining a fairly accurate estimation of .
CE P
Finally, in the context of the hybrid IB technique in this work, the body force at a Cartesian grid point can be represented as a sum of the two components for moving
AC
and stationary boundaries:
Nm
f F h x X Vl f u t l 1
Naturally, the two terms are nonzero only at different times and positions except for some short and localized contact.
2.5 Parallel implementation The hybrid IB technique described above is implemented and parallelized on a supercomputer with both multi-core CPUs and many-core GPUs. The Navier–Strokes equations for the fluid flow are solved on a fixed Cartesian grid, which can be parallelized straightforwardly. In the parallel implementation, the regular rectangle computational domain is divided into a number of sub-domains, and each process (running on a CPU core and the corresponding GPU) is assigned to a sub-domain. The
(20)
ACCEPTED MANUSCRIPT data exchange between sub-domains is communicated through inter-process data transfer in the industry-standard MPI protocol [42].
T
However, the IBs crossing subdomain borders may pose challenges, which will be
IP
tackled in this section. As the solid volume fraction of each grid point is stored locally, no additional difficulty is associated with the parallelization of the stationary approach.
SC R
Meanwhile, for the moving approach, the moving Lagrangian markers at the borders between subdomains cause a variety of complications. The interpolation of the velocity
NU
on markers (Eq. (13)) requires the exchange of extra flow information at the subdomain borders. In general, to minimize the communication, tagging and searching
MA
of Lagrangian points and construction of the delta functions should be done efficiently
AC
CE P
TE
D
[23] .
Figure 1. Schematic representation of the interpolation and distribution operations at the subdomain borders: (a) two layers of ghost Cartesian grid points, (b) one layer of ghost grid points.
The two steps for processing the IB markers, force distribution and velocity interpolation, both loop over the same indexes determined by the support of the regularized delta function, which has a width of four grid points in this study. Taking
ACCEPTED MANUSCRIPT the two-dimensional case in Figure 1(a) for example, we see that the IB maker m interacts with the surrounding 4×4 Eulerian grid points (point 1 to point 16).
T
Considering the staggered grid used, in the extreme case shown in Figure 1(a) when the
IP
Lagrangian markers are very close to the subdomain border, Eulerian grid points 1–8 are ghost grid points whose values should be transferred from the neighboring process.
SC R
Therefore, in general, two layers of ghost grid points are needed to perform the interpolation and distribution steps while only one layer is needed if a second order
during the IB stage in this extreme case.
NU
discretization scheme is applied. Therefore, two additional communications are needed
MA
■ Interpolation step: the data of the ghost Cartesian grid points (grid points 1–8) should be transferred from the adjacent process to compute the velocities of the Lagrangian
the case of the fluid solver.
D
markers (Eq. (13) ). One extra layer of ghost grid points is needed when compared with
TE
■ Distribution step: the body force source of Lagrangian markers should be spread to
CE P
the surrounding connecting Cartesian grid points (grid points 1–16). The velocity of the Cartesian grid points should be updated by adding the source terms to them. As a result, the calculated body force source term of Cartesian grid points 1–8 should be
AC
further transferred back to the neighboring process and added to the corresponding grid points.
Obviously,
these
two
additional
communications
greatly complicate
the
parallelization and increase the communication tremendously. To simplify the above procedure, a partial velocity interpolation and force distribution method is developed in this study, where the Lagrangian markers interact with only some of the 4×4 Eulerian grid points at the subdomain border. In this strategy, only one layer of ghost grid points is required and the interpolation and distribution steps are modified accordingly. It should be noted that only the Lagrangian markers assigned to the present process are computed. In Figure 1(b), for example, only Eulerian grid points 5–16 are used for the velocity interpolation and the body
ACCEPTED MANUSCRIPT force of m is only distributed to grid points 9–16. As a result, overlapping of one ghost grid point is sufficient and the body force does not need to be transferred to the
T
neighboring process any more. Two factors should be incorporated into the
IP
interpolation and distribution steps to ensure that the velocities are interpolated exactly
Cartesian grid points [43, 44]: 1 fv
u x X h
xv
Nm
f l 1
where f v
and Fl are defined by
x X h
xv
h
3
l
x X h
x f
h
3
l
D
Fl
3
l
F h x X l Vl Fl
MA
fv
h
NU
Ul
SC R
and the body forces exerted by the Lagrangian markers are fully transferred to the
TE
Here v is the Cartesian domain for the interpolation of the intermediate velocity at the Lagrangian markers (in the extreme case of Figure 1(b), it includes grid points 5–
CE P
16), and f is the domain for the body force spreading (in the extreme case of Figure 1(b), it includes grid points 9–16). Naturally, when a Lagrangian marker lies more than two grid point widths from the subdomain border, v and f will consist of all 16
AC
surrounding Eulerian grid points, and f v 1 and Fl 1 .Applying this strategy, the communication becomes more efficient. This strategy is derived from and an improvement of the half-distributed forcing strategy for IB methods, which was demonstrated to be feasible in our previous publication [44]. Furthermore, the forces and torques added to the fluid through the distribution procedures (Eq. (16)) are not changed [44]. From the discussions above, we conclude that the proposed partial interpolation and distribution strategy will not lead to additional communications. As a result, the parallelization is much simpler and more efficient.
3. Validation
(21) (22)
(23) (24)
ACCEPTED MANUSCRIPT To evaluate the efficiency of the hybrid IB technique and the accuracy of the proposed parallelization strategy, we simulated the flow field in two stirred tanks
T
with this method. The main results are reported in the follows.
IP
3.1.A simple stirred tank
The method is first applied to a cylindrical unbaffled stirred tank with eight straight
SC R
blades at the mid-height as an impeller. This tank has been studied experimentally [45] and numerically [2]. A schematic illustration of the stirred tank is presented in Figure 2, where the height of the tank H=10 cm. The impeller rotates at n=100 rpm, which
NU
leads to an impeller Reynolds number Re ND2 1024 , in terms of the turbine
MA
diameter D, turbine frequency N and fluid kinematic viscosity . Direct numerical simulation is possible for the present simple stirred tank in which geometrical dimensions and rotation speed are both small. The rectangular computational domain is
D
discretized using a uniform Cartesian grid of 256×256×256 grid points, which is finer
TE
than that used in the DNS by Verzicco et al. [2]. The computational domain is divided
CE P
into 2×2×2 subdomains in the parallel computing. No-slip boundary condition is imposed at the bottom of the tank, while free-slip condition is assigned at the top. Within the computational domain, no-slip boundary conditions are enforced for the
AC
cylindrical tank wall and the impeller using the hybrid IB technique described above.
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 2. Schematic illustration of the stirred tank (adapted from Dong et al. [45] and
D
Verzicco et al. [2] ).
TE
We compare the computational efficiencies of the two IB methods described above in handling the stationary and moving IBs, as listed in Table 1. Note that, if a moving
CE P
boundary is handled by applying the stationary approach (Case 1), the grid volume fraction should be evaluated in each time step because the position of the boundary changes with time. In a stirred tank, the surface area of the rotating impeller is
AC
generally small compared with the stationary wall, and the total number of Lagrangian markers is limited. Conversely, if the tank wall is also simulated with the moving approach (Case 2), the number of Lagrangian markers will increase by approximately an order of magnitude (from 15,611 to 132,183 for the present stirred tank), and the computational time will increase sharply. From Table 1, we see that the time required for the hybrid IB technique is much less than the traditional approach to a single IB for both stationary and moving boundaries. Table 1. Comparison of the efficiency of the two IB methods in the case of a stirred tank. SA and MA denote the stationary approach and moving approach, respectively. Data are the time taken for 100 time steps. In Cases 1 and 3, a small value of M (= 3) is used when evaluating the grid solid fraction in the stationary approach.
Case
Case 1
Case 2
Case 3
ACCEPTED MANUSCRIPT Impeller
Wall
Impeller
Wall
Impeller
Wall
SA
SA
MA
MA
MA
SA
Time (s)
889.646
193.708
92.144
T
Furthermore, the averaged velocity components are compared with experimental
IP
data [45], as shown in Figure 3. Data are evaluated at two different axial locations z/H
SC R
= 0.5 and 0.7. A satisfactory agreement is obtained for all three velocity components in these three cases, except for small deviations in the axial velocity at z/H = 0.7 and in the radial velocity at z/H = 0.5 in Case 1. This may be ascribed to the approximation of
AC
CE P
TE
D
MA
NU
the solid volume fraction due to the small M used.
Figure 3. Comparison of the computed (lines) average velocity components with the experimental data from Dong et al. [45] (symbols ▲) along the radial direction at different axial locations z/H = 0.5 (a, b, c) and 0.7 (d, e, f). ua, ur and ut represent axial, radial and tangential velocity components, respectively.
ACCEPTED MANUSCRIPT
3.2. Turbulent flow in a Rushton stirred tank
T
To verify our method in more practical applications, we then simulated a tank stirred
IP
by a Rushton turbine, and LES is incorporated to describe the turbulent flow in it [21,
SC R
22, 33, 46]. The problem has been studied by Hartmann et al. [38] who presented the comparison of LES and RANS simulation with LDA measurements. The stirred tank is encapsulated by a standard cylindrical vessel 150 mm in diameter
NU
(T = 150 mm) with four equispaced baffles. The geometry is sketched in Figure 4. The density of the working fluid (𝜌) is 1039 kg ∙ 𝑚−3 and the dynamic viscosity (𝜇) is
MA
15.9mPa ∙ s. The rotating speed is 2672 rpm, resulting in an impeller Reynolds number Re = 7300. This Reynolds number is not within the fully turbulent regime, however, the flow is fully turbulent near the turbine [29]. A uniform Cartesian grid of 5123 is
D
used in the computation and divided into 43 subdomains. Starting from the quiescent
TE
condition, the stirred tank flow requires approximately 100 turbine revolutions to reach
CE P
a statically steady state [29]. The start-up was first computed on the 2563 grid for about 100 revolutions and then interpolated onto the 5123 grid for another 15 revolutions. Then, the flow variables were collected from the last 10 revolutions. In the
AC
computation, the timestep was selected according to the tip speed U tip and the grid width h , i.e., tUtip h 0.2 . The computation requires approximately 1122 timesteps for a full revolution on the grid of 2563 and 2244 timesteps on the grid of 5123.
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 4. Cross-section of the tank driven by a Rushton turbine [33].
D
To compare the present LES results with experimental and numerical data from
TE
other publications, the phase averaged variables are evaluated. Figure 5 shows the
CE P
experimental, RANS and LES results of the phase averaged flow field and their associated kinetic energy. Figure 5(d) is the present LES results, Figure 5(c) is the LES results using a standard Smagorinsky subgrid-scale model with constant CS 0.1 [33].
AC
The main features of the flow field are the two large circulation loops on both the top and bottom of the impeller. The working fluid is pumped radially towards the tank wall, partly goes upward above the impeller and partly goes downward below the impeller. Finally, the fluid flows into the impeller region along the axis. In this figure, the RANS and LES for both subgrid-scale models show good agreement to the experiment. The magnitudes of the kinetic energy of the simulation results are comparable with that of the experiment. As can be seen, the kinetic energy in the flow jet area is much higher than the rest of the tank. A quantitative comparison of phase-averaged radial and tangential velocity components at three radial locations is shown in Figure 6. A satisfactory agreement is obtained between the simulation and the experiment except a small deviation for radial
ACCEPTED MANUSCRIPT velocity at 𝑟⁄𝑇 = 0.317 and the tangential velocity at location 𝑟⁄𝑇 = 0.183. All the simulation results overestimate the tangential velocity at location 𝑟⁄𝑇 = 0.183, but
T
the present LES result is better than that of RANS and LES(s). From the comparison,
IP
we can see that the results of RANS are obviously worse than those from LES, especially for the tangential velocity at 𝑟⁄𝑇 = 0.25 and 𝑟⁄𝑇 = 0.317. LES with
SC R
different subgrid-scale models get similar results. Although the advantage is not obvious, the results from the present LES is a bit better in general. For the present LES,
NU
the radial velocity at 𝑟⁄𝑇 = 0.183 and the tangential velocity at 𝑟⁄𝑇 = 0.317 are almost identical to the experimental data, while bigger deviation is seen for the LES
MA
with standard subgrid-scale model (LES(s)).
Figure 7 shows the axial profile of the turbulent kinetic energy at three different radial locations. The deviation of the kinetic energy between the simulation and the
D
experiment is much larger than that of velocity. Underestimation is observed for all the
TE
simulations at the impeller tip (i.e. 𝑟⁄𝑇 = 0.183). RANS predictions significantly
CE P
underestimate the kinetic energy at all the three locations, mesh refinements are not expected to improve substantially the prediction [10]. LES with the standard subgrid-scale model overestimate at locations away from the impeller region (i.e.
AC
𝑟⁄𝑇 = 0.25 and 𝑟⁄𝑇 = 0.317). While the peak values of the kinetic energy predicted by the present LES are comparable with that of the experiment at locations away from the impeller region (i.e. 𝑟⁄𝑇 = 0.25 and 𝑟⁄𝑇 = 0.317), but the width of the profile is smaller. From this study, we can find that the predictions of the present LES-IB method compare satisfactorily in general against the experimental data, and it is at least as accurate as other LES-IB methods. It suggests that this method can be promising for the simulation of stirred tanks at large scale with complex geometries and at high Reynolds numbers.
AC
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 5. Phase averaged plot of the velocity vector field and turbulent kinetic energy. (a)
Experiment, phase averaged turbulent kinetic energy are only measured in the impeller region, (b) RANS, (c) LES with Smagorinsky model, (d) Present LES model. (a)-(c) are taken from Hartmann et al. [38].
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 6. Phase averaged axial profile of the radial (left) and tangential (right) velocity
AC
components at three different radial locations r/T = 0.183 (top), r/T = 0.25 (middle) and r/T = 0.317 (right).
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 7. Phase averaged axial profile of the kinetic energy of the random velocity fluctuations
TE
D
at three different radial locations r/T = 0.183, r/T = 0.25 and r/T = 0.317.
4. Application
of
the
method
to
particle-resolved
CE P
simulation of stirred tanks Solid suspensions in liquid are widely encountered in technical systems. There is
AC
no universal approach for modelling and simulation of such problems, the appropriate method depends on the flow regime and on the required details and accuracy [47], but limited by the computational cost, Eulerian-Eulerian (two fluid model, TFM) and Eulerian-Lagrangian (discrete particle model, DPM) methods are used most often [48, 49]. However, empirical correlations are needed to account for unresolved parts of the suspensions physics in these two models and we cannot obtain the complex flow field around the solid particles. Discrete element method (DEM) , originally developed by Cundall and Strack [50], is a soft particle model which allows overlap between the contacting particles [50, 51]. A linear spring-dashpot force model is used in the implementation. In this study, IBM is coupled with the particle simulator to handle the solids suspension inside stirred tanks. The motion of particles and the particle-particle collisions are treated by DEM
ACCEPTED MANUSCRIPT [51], while no-slip conditions at the surfaces of the particles are assigned by using IBM. This method is usually called particle-resolved simulation [49], where the
T
Eulerian grid width is approximately smaller by an order of magnitude. In this method,
IP
no correlations are required, both the particle-particle and particle-fluid interactions are modelled in a realistic way. There is obviously a computational penalty for resolving
SC R
down to finer scales and the physical size of the domain that can be simulated gets limited. As a result, smaller stirred tanks compared to laboratory scale are studied here.
NU
It should be mentioned that the lubrication force is also considered in the collisions to account for condition when particles get in very close proximity [52].
MA
During the collision, a thin lubrication layer is formed between the particles and the fluid is squeezed out of this gap when the particles approach and is pushed back into the gap during rebound. Viscous forces hence become important in the approaching
D
and rebounding process and can lead to sizeable dissipation. This feature might be
TE
resolved by an adaptive local grid refinement [53], but this usually results in
CE P
substantially increased computation time. An appropriate kind of “subgrid-scale model” accounting for the film between the surfaces should be implemented when it is too thin to be resolved by the Eulerian grid. This contribution can be quantitatively
AC
described by [49, 54]
Flub
Ri2 R 2j 1 1 n uij s s1 6 f f 2 Ri R j s1 s0 Ri2 R 2j 1 1 6 f f n uij s1 s s0 2 Ri R j s s0 s s0 0
where n is the unit vector from particle i to particle j, uij ui u j , s0 0.2R0 and s1 2 104 R0 .
The physical parameters in the computation are identical to the experiments [55] and the tank size is smaller than that in experiment in consideration of the computation cost. A schematic diagram of the stirred tank is shown in Figure 8, and the numerical and physical parameters are shown in Table 2.
(25)
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 8. Schematic illustration of the lab-scale stirred tank. left: side view, right: top view.
MA
𝑫 = 𝟐𝟒 𝐦𝐦，𝐛 = 𝟗. 𝟏𝟕 𝐦𝐦.
Table 2．The numerical and physical parameters in the simulation of solid suspension in a
Value
Tank size
483
TE
Parameter
D
lab-scale stirred tank.
mm
Parameter
Value
Reference length
1
mm -6
0.75
mm
Reference mass
1.15×10
Particle radius (impeller)
0.5
mm
Reference velocity
0.1
ms-1
Particle density
2485
kgm-3
Reynolds number (Reimp)
15850
--
1150
-3
Eulerian grids (∆h)
0.09375
mm
CE P
Particle radius (freely)
Fluid density Gravity Viscosity
4.389e
AC
Particle mass
kgm -6
kg -2
9.81
0.87×10
ms -6
m s
2 -1 -1
-5
kg
∆t
1×10
Markers (per particle)
1008
--
Markers (impeller)
69141
--
s
Rotating velocity
50.14
rads
GPUs used
64
--
Particle numbers (freely)
4731
--
Particle numbers (impeller)
2790
--
Three cases with different geometric configurations are carried out in this study, the setup are listed in Table 3. The average solid volume fractions in these three cases are the same, 𝜙̅ = 10.4%. The no-slip boundary condition is enforced on all surfaces including the impeller, the tank wall, the baffles and the particles using the hybrid method of this study. It is difficult to judge the relative location between the particles and the impeller, for the convenience to deal with particle-impeller collisions, the impeller is represented by the combination of particles with diameter of 𝑑 = 1 mm.
ACCEPTED MANUSCRIPT As a result, there are two sets of representations of impeller: combination of particles during particle-impeller collision; and Lagrangian markers during fluid-solid
T
interactions. The particles are initially located at both the upside and the underneath of
IP
the impeller, as shown in Figure 9.
Table 3. The setup of the three cases in the simulation of solid suspension in a stirred tank.
Baffle
PBTD
1
Yes
Case 2
PBTU
2
Yes
Case 3
PBTU
Case 1
No
1
PBTU: up-pumping pitched blade turbine
15850
NU
PBTD: down-pumping pitched blade turbine
2
Reimp
SC R
Impeller
MA
The computational grid in this simulation is 512×512×512, the resolution of the present grid is not fine enough for DNS, as a result LES with wall-layer model is
D
applied in the simulations. A resolution of 16 grid spacing over a particle diameter is
TE
used to resolve the flow around a solid sphere. CFL is around 0.2 and it takes about 4177 timesteps for a full revolution.
CE P
We investigated the startup of the particle suspension in the three cases. As shown in Figure 9, we can conclude from the qualitative analysis that PBTU is more conducive
AC
to particle suspension than PBTD, and more uniform particle distribution in the radial direction is obtained when baffles are used. This difference is shown more clearly in Figure 10 with the distribution of solids concentration, and more quantitatively in Figure 11 using a global uniformity index (𝜉) [56]
1 2 1 1 NC
1 2
i 1 i 1 NC
2 i 2
where 𝜙𝑖 is the solid fraction in ith grid. The index is conveniently defined so that as 𝜉 → 0 the uniformity of the suspension is at its minimum; and when 𝜉 → 1 the solids are uniformly distributed within the vessel. Figure 12 shows the phase averaged velocity of the liquid phase, solid phase and the solid velocity between the two. The averaged solid velocity is calculated by
(26)
ACCEPTED MANUSCRIPT 𝑢̅𝑠𝑜𝑙𝑖𝑑 = (∑ 𝜙 ∙ 𝑢𝑠𝑜𝑙𝑖𝑑 )⁄∑ 𝜙 . Using cylindrical coordinates, the radial-axial 2D velocity maps were obtained by azimuthally averaging the 3D velocity data and
T
projecting them onto the 2D radial-axial plane. Whilst being two-dimensional, these
IP
flow patterns for both the fluid and the solid are qualitatively consistent with the experimental results of Guida et al. [56], which can be considered as a verification of
AC
CE P
TE
D
MA
NU
SC R
the present results.
Figure 9. Evolution of the solids suspension process: instantaneous realizations of particles locations at four different moments as indicated for the three cases.
AC
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
Figure 10. Time-averaged solids volume fraction in the vertical middle plane (left) and in the horizontal plane with z/H = 0.1 (right).
SC R
IP
T
ACCEPTED MANUSCRIPT
AC
CE P
TE
D
MA
NU
Figure 11. Time variation of global uniformity index in the three cases.
Figure 12. Phase averaged velocity of the liquid phase, solid phase and the difference between the two.
ACCEPTED MANUSCRIPT 5. Conclusions An efficient method for the DNS/LES of complex flow in stirred tanks on fixed
T
Cartesian grids is developed in the framework of the finite difference method. The
IP
moving and stationary boundaries are handled respectively by the direct forcing method using Lagrangian markers and a method based on solid volume fractions. This
SC R
hybrid technique is more efficient than using a single method applied to both stationary and moving boundaries.
NU
Meanwhile, the method is further implemented on a supercomputer with CPU–GPU hybrid architecture. The existence of the IB complicates the parallelization. To
MA
minimize the amount of information that needs to be communicated through the inter-subdomains, an easy partial velocity interpolation and force distribution strategy is developed for the parallelization of the IB, which greatly reduces the communication
D
and simplifies the implementation procedure.
TE
The accuracy and robustness of the present method were validated in two relatively
CE P
simple simulations: fluid flow in a simple stirred tank and turbulent flow in a Rushton tank. Both validation cases demonstrated, respectively, the efficiency of the hybrid IB technique and the capability of the combined LES and IBM implementation.
AC
Furthermore, hybrid IB technique is coupled with DEM to study the particle suspension inside the stirred tanks. In this method, no correlations are required; both the particle-particle and particle-fluid interactions are modelled in a realistic way. As a result, this method can be applied to study the fundamental mechanism of solid suspensions in complex flows. In summary, the present method can be used as a promising tool for simulating turbulent flows in stirred tanks with complex geometries.
ACCEPTED MANUSCRIPT NOMENCLATURE Smagorinsky constant
D
turbine diameter, m
f
additional force term, Nm-3
fv
factor
Fl
factor
F
Lagrangian point force, Nm-3
Fs
hydrodynamic force, N
g
gravity force, ms-2
Gs
external force, N
h
Eulerian grid width, m
H
tank height, m
Is
inertia momentum, kg·m2
k
turbulent kinetic energy , m2s-2
ms
mass, kg
L
length, m
M,n
integer
n
unit vector between two particles, (m,m,m)
N
turbine frequency , s-1
Nm
total number of markers
Ns
torque, N·m
p
pressure, Nm-2
r
radial location, m
Re ̅ 𝑆𝑖𝑗
Reynolds number
t
time, s
u ̃ 𝒖
IP SC R NU
MA
D
TE
CE P
Rate of strain tensor, s-1
AC
Ts
T
Cs
torque, N·m velocity, ms-1 intermediate Eulerian velocity, ms-1
ud
desired velocity , ms-1
un
velocity at nth time step, ms-1
ut
translational velocity, ms-1
U ̃ 𝑼
velocity at Lagrangian marker, ms-1
Ud
desired velocity at the fluid-solid interface, ms-1
Utup
impeller tip speed, ms-1
vs
linear velocity , ms-1
x
Eulerian position vector, m
x0
center coordinate of a rigid body, m
X
Lagrangian position vector, m
z
axial location, m
intermediate Lagrangian velocity, ms-1
ACCEPTED MANUSCRIPT
Greek letters solid volume fraction density of fluid, kgm-3
f
dynamic viscosity, Nsm-2
t
kinematic viscosity, m2s-1
a time ratio in LES model
global uniformity index
rotational velocity, rads-1
s
rotational velocity of the rigid body, rads-1
domain
t
time step, s
Vl
volume of the lth Lagrangian marker, m3
h
Dirac delta function
Operators
partial time derivative, s-1
gradient operator, m-1
divergence operator, m-1
2
Laplace operator, m-1
filter width, m
CE P
TE
D
t
AC
IP SC R
NU
MA
one-dimensional delta function
T
f
ACCEPTED MANUSCRIPT
AC
CE P
TE
D
MA
NU
SC R
IP
T
REFERENCES
AC
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
AC
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
AC
CE P
TE
D
MA
NU
SC R
IP
T
ACCEPTED MANUSCRIPT
ACCEPTED MANUSCRIPT 2015-0417
AC
CE P
TE
D
MA
NU
SC R
IP
T
Graphic abstract