Simulation of Multi-Phase Flow in Autoclaves Using a Coupled CFD-DPM Approach

Abstract

In this work, a numerical simulation study on the mixing characteristics of multiphase flow in an autoclave was carried out using CFD technology. The Eulerian–Eulerian model and discrete phase model (DPM) were employed to investigate the solid holdup, critical suspension speed, nonuniformity of solid suspension, gas holdup distribution, bubble tracks, and residence time during stirring leaching in the autoclave. Experiments validate the accuracy of the numerical model, and the experimental values correspond well with the simulation results. The numerical simulation results show that the solid–liquid mixing is mainly affected by the axial flow, the best agitation speed is 400 rpm, and increasing the speed further cannot make the mixture more homogenous and buildup occurred above the autoclave. The calculated critical suspension speed is 406 rpm, which is slightly lower than that obtained from the empirical formula. The gas phase is mainly concentrated in the vortex area above the blade. When the gas phase is in a completely dispersed state (N = 300 rpm), the average residence time of the bubbles is 5.66 s.

1. Introduction

The autoclave is a pressurized stirred tank reactor that is widely used in the metal extraction processes of aluminum, zinc, lead, nickel, and complex refractory minerals in pressurized hydrometallurgy. In the pressure leaching process of metal, the gas used is generally oxygen. Taking the oxygen pressure leaching of zinc concentrate as an example, the way that oxygen participates in the leaching reaction is to dissolve in the leaching solution, so the leaching reaction is still a solid–liquid reaction process in nature [1]. Thus, the solid–liquid mixing effect in the autoclave has a great influence on the leaching reaction. Additionally, the solubility of oxygen in the leaching system is another crucial factor affecting the leaching reaction. According to the two-film theory of Lewis and Whitman [2], the volumetric mass transfer coefficient k_La can be used to quantify the gas–liquid mass transfer process. Among them, “k_L” is the liquid mass transfer coefficient, which is mainly affected by the density and viscosity of the liquid medium and the gas–liquid affinity. Therefore, most of the current research on enhancing oxygen mass transfer focuses on increasing the gas–liquid interface area “a” which is easier to achieve [3,4,5,6]. However, limited by the aggressive and airtight conditions in the autoclave, it is difficult to quantitatively describe the multiphase mixing effect of solid–liquid or gas–liquid in the autoclave [7,8]. With the development of computer technology, the numerical simulation technology based on computational fluid dynamics (CFD) has been quite mature so far [9,10,11,12]. Computer numerical simulation is the primary method for expanding the reactor design and studying the characteristics of the flow field. It can obtain the characteristic parameters and information of the flow field that are difficult to obtain by traditional experimental methods.

Simulation research of the multiphase flow field in stirred tank reactors conducted by specialists and academics focuses primarily on the investigation of flow field characteristics and multiphase mixing characteristics. Ramírez-Torres et al. [13] performed an analysis of the computational hydrodynamic continuous stirred-tank behavior. The findings indicate that the optimal impeller configuration can reduce energy use by up to 9%. Gu et al. [14] used a CFD-PBM coupled model to study the hydrodynamic characteristics of the stirred tank. It was discovered that a self-similarity impeller can improve the uniformity of gas holdup distribution and facilitate the gas dispersion process. Cheng et al. [15] conducted a simulation study on the macroscopic mixing characteristics of gas–liquid–liquid in a turbo-propeller stirred tank by using the Euler model based on CFD software. The reliability of the model is verified by the experimental data of immiscible liquid–liquid mixing. Using the tracer, the mixing times and flow traces of the experiment are compared with the simulation results, and the findings demonstrate that the Reynolds stress model can accurately predict the flow field and mixing time. In order to calculate the local velocities, local concentrations, and concentration gradients inside the adsorber, García-Hernández et al. [16] carried out a 3D numerical simulation of the adsorption process. Several kinetic and transport characteristics were recorded while varying the stirring rate (30–200 rpm) and the initial pyridine concentration (99–487 mg/L). Fan et al. [17] studied the effect of adding baffles on the structure of the circulating flow in the stirring tank. The results show that the baffled agitated tank shows a larger overall circulation capacity for flow rates of different liquid levels, and the root mean square of the axial and radial pulsation velocities is stronger than that of the unbaffled tank. The addition of baffles enhances the overall fluctuation. Overall, the majority of simulation studies of multiphase flow concentrate on the mixing characteristics and optimization of the impeller or tank construction parameters. However, what is neglected by the majority of research is that the residence time of the bubbles is also critical as it relates to the mass transfer process.

On the basis of CFD theory, the finite volume approach is utilized to create a multiphase flow model in the autoclave, and the flow field characteristics and multiphase mixing characteristics are investigated. The distribution of solid holdup, critical suspension speed, nonuniformity of solid suspension, gas holdup distribution, bubble tracks, and residence time of bubbles were studied to describe the main characteristic parameters of multiphase flow in the autoclave. This work can provide a theoretical basis for enhanced leaching in the pressure leaching process and theoretical support for the design and enlargement of the autoclave.

2. Mathematical Model

Whether simulating gas–liquid or solid–liquid multiphase flow, describing the turbulent flow characteristics of the system is crucial. At present, there are three main types of methods for solving turbulent flow: direct numerical simulation (DNS) [18,19,20], large eddy simulation (LES) [21,22,23], and Reynolds averaged Navier–Stokes equations (RANS) [24,25,26]. Direct numerical simulation and large eddy simulation demand greater computational resources and time compared with the Reynolds-averaged Navier–Stokes equations. As a result, the RANS method is the primary method for solving the turbulent characteristics of multiphase flow [27,28,29,30,31]. People are primarily interested in the time-averaged velocity and turbulent flow characteristics in practical engineering applications [32,33,34,35]. Hence, the purpose of the Reynolds average is to obtain the average of the macroscopic changes in physical quantity.

For multiphase systems, the continuity equation and momentum equation can be expressed as:

Continuity equation

$$\frac{\partial \left({\rho }_n{\alpha }_n\right)}{\partial t}+\nabla \cdot \left({\rho }_n{\alpha }_n{u}_n\right)=S_n$$

(1)

Momentum conservation equation

$$\frac{\partial }{\partial t}\left({\rho }_n{\alpha }_n{u}_n\right)+\nabla \cdot \left({\rho }_n{\alpha }_n{u}_n{u}_n\right)=-{\alpha }_p\nabla P+\nabla \cdot \left({\alpha }_p{\tau }_p\right)+F_p+{\rho }_n{\alpha }_{n}g$$

(2)

where n represents the phase, ${\tau }_p$ is the viscous stress, $F_p$ is the interphase momentum transfer, u is the instantaneous flow velocity, and $S_n$ is the interphase mass transfer since this work does not consider an interphase mass transfer, so $S_n=0$.

The RANS is obtained by time-averaging the above basic conservation Equations (1) and (2):

Continuity equation

$$\frac{\partial \left({\tilde{\rho }}_n{\tilde{\alpha }}_n\right)}{\partial t}+\nabla \cdot \left({\tilde{\rho }}_n{\tilde{\alpha }}_n{\tilde{u}}_n\right)=0$$

(3)

Momentum conservation equation

$$\frac{\partial }{\partial t}\left({\tilde{\rho }}_n{\tilde{\alpha }}_n{\tilde{u}}_n\right)+\nabla \cdot \left({\tilde{\rho }}_n{\tilde{\alpha }}_n{\tilde{u}}_n{\tilde{u}}_n\right)=-{\tilde{\alpha }}_n\nabla \tilde{P}+\nabla \cdot \left({\tilde{\alpha }}_n{\tau }_n\right)+{\tilde{F}}_n+{\tilde{\rho }}_n{\tilde{\alpha }}_{n}g$$

(4)

where $\tilde{\alpha }$ is the instantaneous phase volume fraction, all other variables are instantaneous values after time averaging, and $\tilde{P}$ is assumed to be the pressure field common to all phases.

2.1. Eulerian Model

In order to solve the governing equations of the multiphase flow model, the conservation equation of the phase volume fraction must be introduced. The volume V_n of the n phase in the Euler model is defined as Equation (5):

$$V_n={\displaystyle {\int }_V{\alpha }_{n}dV}$$

(5)

where ${\alpha }_n$ is the volume fraction of the n phase, ${\displaystyle {\sum }_{i=1}^n{\alpha }_n}=1$.

The continuous phase equation of the n phase in the multiphase flow calculation is Equation (6):

$$\frac{\partial }{\partial t}\left({\alpha }_n{\rho }_n\right)+\nabla ·\left({\alpha }_n{\rho }_n{\overrightarrow{v}}_n\right)={\displaystyle {\sum }_{i=1}^n{\dot{m}}_{in}}$$

(6)

where ${\overrightarrow{v}}_n$ is the velocity of the n phase and ${\dot{m}}_{in}$ is the mass transfer from the i phase to the n phase.

For the n phase, the momentum conservation equation is as Equation (7):

$$\begin{array}{l}\frac{\partial }{\partial t}\left({\alpha }_n{\rho }_n{\overrightarrow{v}}_n\right)+\nabla ·\left({\alpha }_n{\rho }_n{\overrightarrow{v}}_n{\overrightarrow{v}}_n\right)=-{\alpha }_n\nabla p+\nabla ·{\stackrel{=}{\tau }}_n\\ +{\displaystyle {\sum }_{i=1}^n({\overrightarrow{R}}_{in}+{\dot{m}}_{in}·{\overrightarrow{v}}_{in}})+{\alpha }_n{\rho }_n\left({\overrightarrow{F}}_n+{\overrightarrow{F}}_{life,n}+{\overrightarrow{F}}_{v_m,n}\right)\end{array}$$

(7)

where ${\stackrel{=}{\tau }}_n$ is the pressure strain tensor of the n phase and can be represented by Equation (8):

$${\stackrel{=}{\tau }}_n={\alpha }_n{\mu }_n\left(\nabla {\overrightarrow{v}}_n+{\overrightarrow{v}}_n^T\right)+{\alpha }_n({\phi }_n-\frac{2}{3}{\mu }_n)\nabla ·{\overrightarrow{v}}_n\stackrel{=}{I}$$

(8)

where ${\phi }_n$ and ${\mu }_n$ are the shear viscosity and bulk viscosity of the n phase, respectively. ${\overrightarrow{R}}_{in}$ is the interphase force, ${\overrightarrow{F}}_n$ is the external body force, ${\overrightarrow{F}}_{life,n}$ is the lift, and ${\overrightarrow{F}}_{v_m,n}$ is the virtual mass force.

2.2. Discrete Phase Model

Fluent solves the discrete phase force balance by solving the interphase force equation in Lagrangian coordinates and judging the discrete phase particle trajectory [36]. The force balance of particles in the fluid phase can be expressed by Equation (9):

$$m_d\frac{d{\overrightarrow{u}}_d}{dt}=m_d{F}_D\left(\overrightarrow{u}-{\overrightarrow{u}}_d\right)+m_d\frac{\overrightarrow{g}({\rho }_d-\rho )}{{\rho }_d}+\overrightarrow{F}$$

(9)

where $m_d$ is the particle mass, $\overrightarrow{u}$ and ${\overrightarrow{u}}_d$ are the fluid phase and particle phase velocities, respectively, $\overrightarrow{F}$ is the additional force, and $m_d{F}_D\left(\overrightarrow{u}-{\overrightarrow{u}}_d\right)$ is the drag force on the particle, where $F_D$ can be defined as:

$$F_D=\frac{18\mu }{{\rho }_d{d}_d^2}\frac{C_{d}Re}{24}$$

(10)

where $\mu$ is the fluid dynamic viscosity, $Re$ is the relative Reynolds number, and $C_d$ is the drag coefficient related to the Reynolds number and particle shape factor.

For some additional forces on the particle, it can be ignored relative to other forces on the particle. However, its virtual mass force and pressure gradient force cannot be simplified for calculation when the fluid phase density is much greater than the particle mass ($\rho \gg {\rho }_d$).

The virtual mass force is the force acting on the particle when it is accelerated by the fluid around the particle, its expression is Equation (11):

$$F_v=m_d{C}_{vm}\frac{\rho }{{\rho }_d}\frac{d}{dt}(u-u_d)$$

(11)

where $C_{vm}$ is the virtual mass force coefficient and the default value is 0.5.

The pressure gradient force refers to the additional force due to the fluid pressure gradient near the particle and its expression is Equation (12):

$$F_p=m_d\frac{\rho }{{\rho }_d}{u}_d\frac{\partial u}{\partial x}$$

(12)

According to the definition of stirring Reynolds number, its expression is Equation (13):

$$Re=I^{2}N\rho /\mu$$

(13)

where $\rho$ is the liquid phase density and $\mu$ is the liquid phase dynamic viscosity. Since the stirring speed considered in this paper is 100–500 rpm, the range of Reynolds numbers considered is 1.6 × 10⁴–8.3 × 10⁴.

3. Solution Strategy

3.1. Computational Domain and Mesh Distribution

In this work, SpaceClaim is used for modeling, and the computational domain is shown in Figure 1a. During modeling, the autoclave body is split into the main fluid domain and the stirring fluid domain, and the BOI geometry is constructed at the stirring fluid domain to encrypt the stirring fluid region locally. The structural parameters are shown in

Table 1. Structural parameters of autoclave.

Structure Parameters	Value	Structure Parameters	Value
Tank height	26 cm	Diameter of ring sparger	19.2 cm
Diameter of autoclave	25 cm	Width of outlet	0.23 cm
Bottom ellipse height	5 cm	Blade width	1.7 cm
Diameter of impeller	10 cm	Blade thickness	0.17 cm
Axial location of impeller	6.8 cm	Chord angel	30°

. The type of impeller is a pitched blade turbine, which is the most commonly used axial flow impeller and is widely used in multiphase circulating mixing systems such as those in biology, metallurgy, and the chemical industry. The meshing is performed by the fluent meshing model, and the grid scheme is poly-hexcore [37]. This grid type employs a hybrid grid configuration and it has the advantage of drastically reducing the number of cells and conserving computational resources. The final generated mesh is shown in Figure 1b. The geometric construction of the pitched blade turbine and the generated mesh are shown in Figure 1c and Figure 1d, respectively.

3.2. Boundary Conditions

In this work, the MRF (multiple reference frame) method is used [38,39] for the rotation modeling, and the dynamic and static regions are calculated using different reference frames. This method considers the impeller and its surrounding area to be a moving area, whereas the rest of the tank is considered to be a static area. The transport equations in different regions of the MRF method are separately solved, and information transfer at the region boundaries is taken into account. The tank walls and impeller surfaces were treated as no-slip boundaries with standard wall functions. The gas enters the domain of fluids via the gas sparger and the speed is controlled to be 1 × 10⁻⁶ kg/s. The diameter of solid particles is controlled to be 0.002 m, and the density is 2000 kg/m³. Both solid–liquid and gas–liquid simulations have pressure conditions of 1 MPa.

In solid–liquid multiphase flow simulations, only the Schiller–Naumann drag is considered. The SIMPLE algorithm was used for pressure–velocity coupling. In addition, the first-order upwind discretization scheme was applied for convection terms. Due to the low density of the gas phase and the significant difference between it and the density of water, the bubble–particle phase of the gas–liquid simulation must account for the Saffman lift force and the pressure gradient force. The phase-coupled SIMPLE method was used to solve the coupling between pressure and velocity. The calculation is transient, time step is 0.01 s, and convergence residual variable is set to 10⁻⁴. All calculations were run on the same workstation, which included a 3.0-GHz CPU (AMD Ryzen Threadripper 2990WX 32-Core Processor) and 64 GB of RAM.

4. Results and Discussion

4.1. Grid Independence

While a smaller grid discretization error results in more precise simulations, it also requires more processing power and slows down convergence. In general, it is required to strike a balance between computation precision and calculation amount, and grid independence verification is used to determine the ideal number of grids [40]. The fundamental concept is to create a compromise between calculation accuracy and grid number by evaluating the deviation of a variable under different grid numbers of models under different grid number variables. In this work, the averaged flow field velocities near the impeller tip are compared using different grid numbers while all other variables remain unchanged. All simulations were carried out at N = 150 rpm. The flow field velocity changes for five grid numbers are shown in

Table 2. Grid independence analysis.

Grid Solutions	Mesh Number	Velocity (m/s)
Solution 1	264,853	0.56
Solution 2	373,245	0.61
Solution 3	422,658	0.74
Solution 4	456,258	0.75
Solution 5	468,523	0.75

. As the number of grids increases, the change in the velocity of the flow field decreases gradually. Considering the computational cost issue, the grid finally adopts solution 3.

4.2. Experimental Verification

This work uses PIV [41] technology to experimentally verify the computational model. The equipment used is the PIV system produced by TSI Company in the United States, which consists of four parts: a tracer particle adding device, a double-pulse laser+ light sheet system, a CCD camera, and an information processing system. Through preliminary experiments on various impeller speeds and different planes, the stirring speed was selected as 150 rpm, and the exposure time delay was determined to be 1000 μs. The tracer uses a PSP tracer with a diameter of 50 μm, which is close to the density of water and has good followability and stability. The velocity component is calculated by DANTEC software, the captured image is set as a query area of 32 × 32 pixels, and the data of 200 images is time-averaged to obtain a reasonable velocity flow field. Comparing the PIV experimental data with the CFD simulation results, the streamline is shown in Figure 2.

The simulation-predicted value and the experimental value of the liquid velocity are compared in Figure 3. It can be observed that the simulation results of the CFD model are reliable and in good agreement with the experimental data.

4.3. Solid Holdup Distribution

Figure 4 shows the volume fraction distribution of the solid phase in the autoclave at different times when the stirring speed is 200 rpm. It is clear that the axial flow in the autoclave has the greatest influence on solid-phase dispersion, and the solid phase diffuses with the large circulating flow above the blades. Because of the impact of the blade jet, the solid phase diffuses to the upper section of the autoclave body in the first 4 s. Within 4–8 s, the solid-phase particles are affected by the large circulating flow above the autoclave and spread throughout the autoclave body. After about 10 s, the solid-phase distribution in the autoclave reached the stable state. Due to the low stirring speed, there was some solid-phase accumulation and insufficient solid-phase dispersion at the bottom of the autoclave after the steady state was reached. The volume fraction of most of the solid phase in the autoclave was about 2.2%.

Figure 5 shows the volume fraction distribution of the solid phase in the autoclave when the system reaches a steady state under different agitation speed conditions. When the speed is low (200 rpm to 300 rpm), solid-phase deposition occurs at the bottom of the autoclave and the solid-phase distribution above the autoclave is low, resulting in poor mixing. When the rotational speed was increased to 400 rpm, the solid-phase distribution in the autoclave became uniform, and the solid-phase distribution in the upper portion improved. When the stirring speed was increased to 500 rpm and 600 rpm, the accumulation of solid-phase suspensions appeared above the autoclave. This is due to the high turbulent kinetic energy in the autoclave caused by the high stirring speed. The lift force generated by the mainstream axial flow in the autoclave is large, and the solid phase is back-accumulated in the lower flow rate area above the autoclave.

In order to describe the solid-phase distribution inside the autoclave more visually, the solid-phase volume distribution at different locations inside the autoclave was monitored. Figure 6a is the solid-phase volume fraction curve in the horizontal direction at the Y = 0 section at the bottom of the autoclave. When the stirring speed is low, the volume fraction of the solid phase at the bottom of the autoclave will peak, which indicates that there is an accumulation of solid phase at the bottom of the autoclave. When the agitation speed is high, however, the solid holdup at the autoclave’s bottom tends to be uniform, and the phase volume fraction decreases. Due to the high speed, the majority of solids in the autoclave gather above the impeller.

Figure 6b is the solid-phase volume fraction curve in the axial direction at r = 0.3 m on the Y = 0 section. When the stirring speed is low (200 rpm to 300 rpm), the volume fraction of the solid phase increases and the part of the solid holdup in the autoclave decreases. The volume percentage of the solid phase above the autoclave grows as the rotation speed increases, and a peak appears at Z = 0.8 m, suggesting that the solid-phase aggregates at this position and the overall solid holdup uniformity in the autoclave declines, which is not conducive to mixed flow.

4.4. Critical Suspension Speed

The critical suspension speed is an important basis for judging the solid–liquid dispersion performance of the autoclave. It was first proposed by Zwietering [42], and it is also the most ideal operating condition for the autoclave. In general, the critical suspension speed is determined by the speed at which the solid phase remains at the bottom for less than 12 s. The empirical formula is as shown in Equation (14). We determined the average solid volume fraction at the Z = 0.05 m level at the bottom of the autoclave at 200 rpm, 300 rpm, 400 rpm, 500 rpm, and 600 rpm, respectively, by examining the distribution of the solid holdup at the bottom of the autoclave at various speeds. The critical suspension speed corresponds to the junction of the two tangents, which is drawn at the maximum slope and the minimum slope of the curve. Figure 7 is a schematic diagram for judging the critical suspension agitation speed under different agitation speed conditions.

$$N_{js}=K{\left(\frac{g({\rho }_p-\rho )}{\rho }\right)}^{0.45}\frac{{(100C_V{\rho }_p/\rho )}^{0.13}{d}_p^{0.2}{v}^{0.1}}{D^{0.85}}$$

(14)

With the increase in stirring speed, the average solid-phase volume fraction in the bottom area of the autoclave decreased. Then, make the maximum slope and the minimum slope tangent, and from the figure, it can be found that the coordinates of the intersection point O is (406, 0.0322). Therefore, the critical suspension speed of the autoclave is N_sp = 406 rpm. Take the shape factor of the autoclave as K = 1.4, the solid-phase particle density is 2000 kg/m³, the diameter is 0.1 mm, the solid holdup is C_v = 0.042, the blade diameter is 0.1 m, and the dynamic viscosity of water is 1.005 × 10⁻³ N∙s/m². After calculation, the result obtained is N_js = 432 rpm, which is slightly higher than the simulated data.

4.5. Nonuniformity of Solid Suspension

Under actual working conditions, there are velocity and pressure differences in the flow field inside the autoclave, and the solid–liquid distribution cannot be completely uniform. Calculating the homogeneity of solid-phase dispersion is one of the criteria used to judge the quality of multi-phase mixing in the autoclave. The more uniform the solid-phase dispersion, the better the autoclave’s stirring effect and the more conducive mass transfer between phases. Take 5 positions of Z = 0.1 m, 0.3 m, 0.5 m, 0.7 m, and 0.9 m in the Z direction of the autoclave section, and calculate the average solid-phase volume distribution at different positions. Then, use Equation (15) to calculate the inhomogeneity of the solid-phase distribution at this speed. Take a stirring speed of 500 rpm as an example. Figure 8a is the solid-phase volume fraction curve at different horizontal positions of the Y = 0 section of the autoclave when the stirring speed of the autoclave is 500 rpm. With the increased position of the autoclave, the volume distribution of the solid phase shows an increasing trend, and the volume fraction of the solid phase is the highest at Z = 0.7 m, while the solid-phase content decreases near the stirring axis.

$$\sigma =\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(\frac{{\alpha }_i}{\alpha }-1\right)}^2}}$$

(15)

where α_i is local solid holdup and α is the average solid holdup of the autoclave.

The inhomogeneity of the solid-phase suspension at the stirring speed is determined by calculating the average solid-phase volume fraction corresponding to each curve. Figure 8b shows the calculated suspension unevenness at different agitation speeds. In the autoclave, as the stirring speed increased, the non-uniformity of the solid-phase suspension first decreased and subsequently increased. At N = 400 rpm, the inhomogeneity of the solid-phase suspension is the lowest (σ = 8.65 × 10⁻²), the solid–liquid mixing in the autoclave is the most uniform, and the stirring effect is the best.

4.6. Gas Holdup Distribution

The impeller we use is a typical axial-flow impeller. During the ascent of the gas phase, it will interact with the blade jet in the blade region. When the stirring speed is high, the shear force generated by the blade jet is large, and the gas phase will diffuse to the bottom of the autoclave under the influence of drag force before ascending along the wall. Figure 9 shows the distribution of gas holdup in the autoclave at different stirring speeds.

When the stirring speed is N = 100 rpm, the air bubbles are less affected by the blade jet. The gas phase distribution is a transition state between the gas flooding state and the remaining carrier gas state and is only partially affected by the flow field. The effect of phase separation and dispersion is poor, and there is a phenomenon of gas aggregation in the vortex area of the blade, which is not conducive to gas phase diffusion and mass transfer. When the stirring speed was increased to 300 rpm, however, the autoclave’s gas phase was entirely dispersed. Under the action of the blade jet, the gas phase is spread throughout the entire autoclave body, the distribution is more uniform, and the range of the dead zone is reduced. Figure 10 shows the gas phase distribution curve in the axial direction of the autoclave at different stirring speeds. With the increase in agitation speed, the overall gas holdup in the autoclave increases. This is due to the existence of circulation in the gas phase. The residence time thus becomes longer, and the gas–liquid mixing effect is good. In addition, the gas phase is mostly concentrated in the flow field’s vortex core region. As the blade speed increases, the gas phase concentration in the vortex core region becomes more apparent.

4.7. Bubble Tracks and Residence Time

The DPM obtains the particle tracks by calculating the force balance and motion equation of each particle phase. Figure 11 shows the trajectory of the bubble movement in the autoclave when the stirring speed is N = 400 rpm. Via monitoring of the overall gas holdup, it is determined that the flow field has achieved a steady state when the bubble dispersion time is 5.6 s. However, once the flow field had stabilized, we proceeded to observe bubble dispersion to determine the maximum bubble residence time. When the agitation speed is high, the bubble movement is mainly affected by the drag force. The bubbles are introduced by the ring sparger and propelled by the blade’s jet during the ascent process. Under the influence of flow field drag, they rotate around the axis. Under the influence of the axial flow blade, gas phase diffusion is comparatively adequate. The gas phase will diffuse to the bottom of the autoclave and all parts of the autoclave due to the action of axial flow, and the overall gas–liquid dispersion effect of the autoclave is improved.

The bubble residence time evaluates the autoclave’s gas phase circulation performance, which has a direct impact on the gas–liquid mass transfer time. This is carried out by capturing the bubble particles that reach the top of the autoclave under different stirring speeds and counting their particle information to obtain the bubble residence time at different agitation speeds, as shown in

Table 3. Bubble residence time at different agitation speeds.

Agitation Speed (rpm)	Minimum Residence Time (s)	Maximum Residence Time (s)	Average Residence Time (s)	Variance
100	1.56	10.16	5.03	1.89
200	1.45	10.30	5.24	2.07
300	1.37	11.13	5.66	2.18
400	1.34	11.28	5.75	2.22

. With the increase in stirring speed, the minimum residence time of bubbles decreases and the maximum residence time increases. This is due to the growing impact of turbulence on bubble movement and the increase in residence time variation. It demonstrates that the regularity of bubble residence time distribution is diminished and that the residence time is affected by high dispersion of turbulent shear. The average residence time of bubbles increases with the increase in agitation speed, which is conducive to gas phase dispersion and mass transfer. At 300 and 400 rpm, the bubble residence time is closer. This is due to the fact that the gas phase has reached a condition of complete dispersion, and increasing the agitation speed has little influence on the dispersion of the gas phase.

Figure 12 shows the distribution of bubble particles in the autoclave at different times. When t = 2 s, the bubble particles had just arrived at the blade area and were starting to disperse under the influence of the blade jet; when t = 4 s, the bubbles had diffused throughout the autoclave, and a small amount of bubble particles began to overflow. When t = 6 s, the bubble particles further increased, the overall concentration of bubble particles above the autoclave increased, and the phenomenon of bubble particle circulation began to appear. When t = 8 s, the phase of bubble particles increases, and the particles with a higher residence time mainly appear in the area near the stirring shaft. This may be due to the fact that this portion of the region is the recirculation region of some major circulations on the flow field, and some bubbles have circulating diffusion in this region. In addition to slowing the ascent of bubble particles, the recirculation zone increases the bubbles’ residence time in this region.

The particles are captured on the pressure outlet boundary of the autoclave, and the residence time distribution of the bubble particles can be quantitatively calculated using statistical particle information (Figure 13). We counted the residence time of some particles. The residence time distribution of the bubble particles shows a tailing phenomenon, and the peak value of the residence time appears at t = 5 s.

5. Conclusions

In this work, a multiphase flow model was used to simulate the multiphase mixing characteristics in the autoclave. The simulation studies on the solid–liquid mixing characteristics and gas–liquid dispersion in the autoclave were carried out, respectively. The conclusions are as follows:

(1)

The solid-phase volume distribution clouds were analyzed at different times. The solid–liquid mixing was greatly influenced by the axial flow of the autoclave. The solid-phase deposition at the bottom of the autoclave decreased as the stirring speed increased. The solid-phase distribution in the autoclave became more uniform when the rotation speed was increased to 400 rpm. If the agitation speed is increased again, the phenomenon of solid-phase accumulation above the autoclave will arise, which is not conducive to solid-phase dispersion and mass transfer in the autoclave, and the autoclave’s power consumption will increase.

(2)

At various speeds, the average solid-phase volume fraction at the bottom of the autoclave is used to determine the critical suspension speed of solid–liquid mixing. After processing and calculating the simulated data, N_sp = 406 rpm was found to be the critical suspension speed for mixing solids and liquids in the autoclave. Its value is marginally lower than the critical suspension speed calculated using the empirical equation. By examining the suspension inhomogeneity at various stirring speeds, the solid-phase suspension inhomogeneity exhibited a decreasing and subsequently increasing pattern. The nonuniformity of the solid suspension is lowest when the stirring speed is 400 rpm (σ = 8.65 × 10⁻²).

(3)

The contours of the volume distribution of the gas phase in the autoclave under various agitation speed conditions were examined. The gas phase dispersion in the autoclave is totally diffused at 300 rpm stirring speed. With the increase in agitation speed, the gas content of the autoclave increases, and the gas phase is mainly concentrated in the vortex area above the blades.

(4)

Analysis was conducted on the bubble tracks and residence time. The drag force of the flow field mostly influences the movement of the bubbles. In the vicinity of the agitator shaft, air bubbles remain longer in the fluid recirculation zone. At the exit, we measured the time of bubble residence. The average bubble residence time increases as stirring speed increases. When the gas phase is in a completely dispersed state (N = 300 rpm), the average bubble residence time is 5.66 s, and increasing the stirring speed cannot effectively improve the gas phase residence time.