Jia Cheng
| Chunhui Zhang* | Huizhe Sun |Shaoxiong Zhang | Jiazhuang Zhao
© 2024 The authors. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).
OPEN ACCESS
With the increasing demands for groundwater management in industrial and environmental engineering, the research and application of air sparging remediation technology have become crucial for enhancing multiphase flow control. Particularly under complex geological conditions with temperature gradients, understanding and predicting the behavior of multiphase flows are essential for improving the efficiency and safety of remediation technology. Currently, most studies focus on the analysis of multiphase flows under static conditions, with fewer assessments on the impact of dynamically changing temperatures and pressures, thus limiting the accuracy and reliability of models in practical applications. Addressing this research gap, this study proposes a fully transient temperature-pressure field coupled model, along with a corresponding iterative solution algorithm, to accurately simulate the dynamic characteristics of underground multiphase flows under varying temperature gradients. The progress of this research extends our understanding of the mechanisms of multiphase flow in air sparging, supporting technological innovation and environmental protection efforts in related fields. It also offers new solutions for resource development and pollutant management under complex geological conditions.
air sparging remediation technology, multiphase flow, temperature gradients, temperature-pressure field coupled model, iterative solution algorithm, groundwater management
In many industrial and environmental engineering fields, air sparging remediation technology has become crucial due to its effectiveness in managing and controlling groundwater flow and its associated pollutants. Especially in areas such as mine drainage, soil remediation, and geothermal energy development, this technology has a broad application prospect [1, 2]. However, the complexity of underground multiphase flow, especially under the influence of temperature gradients, presents significant challenges for accurately predicting and controlling fluid behavior [3-7]. Since temperature changes can significantly affect the physical properties of fluids, such as density and viscosity, as well as the interactions between fluids, the impact of temperature gradients on multiphase fluid systems has become a key issue that must be thoroughly researched [8, 9].
The application of air sparging remediation technology in environmental engineering and resource development is becoming more frequent, thus, understanding how temperature gradients affect the behavior of multiphase flows is crucial [10-12]. This not only helps to improve the efficiency and safety of existing technologies but may also provide a theoretical basis for the development of new technologies. Moreover, with the increasing requirements for environmental protection, designing efficient and environmentally friendly sparging systems has become particularly urgent [13-16]. Therefore, researching the impact of temperature gradients on multiphase flow, especially in the context of air sparging technology, has significant practical significance for advancing technology in this field [17, 18].
Although there are several studies focused on the simulation and prediction of multiphase flows, most research has concentrated on steady-state analysis and has not fully considered the dynamic changes of temperature and pressure over time and their impact on the characteristics of multiphase flows [19-22]. In addition, the models used in existing research often overlook the complexity of changes in groundwater physical parameters with temperature and pressure, limiting the applicability and accuracy of models under actual complex geological conditions [23, 24].
Given the background and limitations of existing research, the goal of this paper is to deepen the understanding of the mechanisms of multiphase flow under the influence of temperature gradients in air sparging remediation. Firstly, the paper constructs a fully transient temperature-pressure field coupled model that fully considers the changes in the physical parameters of groundwater in the vadose zone over time, temperature, and pressure, providing a more accurate method for simulating and analyzing multiphase flows. Secondly, to efficiently handle the computational challenges caused by complex boundary conditions and nonlinear factors, this study proposes a novel model iterative solution algorithm. Through these two core research contents, this paper not only supplements the shortcomings of existing research but also provides reliable theoretical and technical support for solving problems in practical engineering, having significant academic and practical value.
2.1 Assumption conditions
In constructing the air sparging transient temperature-pressure mathematical model, to simplify the model and make it solvable, a series of assumptions are usually made. The following are the assumptions that may apply to this study: 1) It is assumed that the permeability and other related physical properties of the underground medium are uniform in all directions, that is, the medium is isotropic. 2) The permeability and porosity of the underground medium are uniformly distributed in space, ignoring the heterogeneity of geological structures. 3) Fluid flow is time-varying, that is, considering the impact of time factors on fluid flow. 4) The medium's fluid may be partially saturated, that is, there is two-phase flow of gas and liquid. 5) It is assumed that there is a consistent temperature gradient throughout the study area, and this gradient remains constant over time or its change can be predicted. 6) It is assumed that fluid and solid media can quickly reach thermodynamic equilibrium in the transient process, ignoring non-equilibrium heat conduction and convection processes. Figure 1 provides a schematic diagram of the slope model.
Figure 1. Schematic diagram of slope model
2.2 Mass and quantity conservation equations
The groundwater system typically contains gas, liquid, and solid phases, and in multiphase flow, the behavior of one phase affects the others. To accurately simulate the movement and interaction of these phases in the underground environment, mass conservation equations must be established for each phase.
In the air sparging transient temperature-pressure mathematical model, the substances in the gas, liquid, and solid phases specifically refer to the gas components flowing in the underground medium, including oxygen, nitrogen, carbon dioxide, and other gases, the water in the underground medium including groundwater and possible other liquid pollutants, and the underground medium itself, including soil, sand, clay, rock, etc. For the gas phase, the mass conservation equation describes the change in gas mass per unit volume over time, including density changes caused by pressure and temperature changes, and the mass balance of gas inflow and outflow. The liquid phase mass conservation equation describes the movement of liquids in pore spaces, including percolation, adsorption, desorption, and interaction with the gas phase. Although the solid phase is usually considered immobile, in some cases, it may also participate in mass exchange. The solid phase mass conservation equation considers the conservation of substances in the solid medium, as well as the interaction between the solid phase and the liquid and gas phases. Assuming the flow channel area is represented by S, the densities of the gas, liquid, and solid phases are represented by ϑh, ϑ1, and ϑa, the volume fractions of the gas, liquid, and solid phases are represented by βh, β1, and βa, the actual flow rates of the gas, liquid, and solid phases are represented by ch, c1, and ca, and the unit solid phase production rate is represented by wa, the following equations provide the expressions for the mass conservation equations of the gas, liquid, and solid phases:
$\frac{\partial }{\partial y}\left( {{\vartheta }_{h}}{{\beta }_{h}}S \right)+\frac{\partial }{\partial x}\left( {{\vartheta }_{h}}{{\beta }_{h}}{{c}_{h}}S \right)=0$ (1)
$\frac{\partial }{\partial y}\left( {{\vartheta }_{1}}{{\beta }_{1}}S \right)+\frac{\partial }{\partial x}\left( {{\vartheta }_{1}}{{\beta }_{1}}{{c}_{1}}S \right)=0$ (2)
$\frac{\partial }{\partial y}\left( {{\vartheta }_{a}}{{\beta }_{a}}S \right)+\frac{\partial }{\partial x}\left( {{\vartheta }_{a}}{{\beta }_{a}}{{c}_{a}}S \right)={{w}_{a}}$ (3)
In the air sparging transient temperature-pressure mathematical model, the flow patterns of air-liquid two-phase flow mainly depend on factors such as the relative flow rates of gas and liquid, flow direction, the porous structure of the medium, and fluid properties. In underground multiphase flow systems, common flow patterns include bubble flow, droplet flow, slug flow, annular flow, string flow, and unstable flow.
In the study of air sparging remediation technology, in addition to the mass conservation equations, momentum conservation equations are also an essential part of simulating multiphase flow. The momentum conservation equation, sometimes referred to as the Navier-Stokes equation, describes the change in fluid particle velocity over time and the interaction between various forces in multiphase flow. The flow patterns of air-liquid two-phase flow in the air sparging transient temperature-pressure mathematical model primarily depend on factors such as the relative flow rates of gas and liquid, flow direction, the porous structure of the medium, and fluid properties. Common flow patterns in underground multiphase flow systems include bubble flow, droplet flow, slug flow, annular flow, string flow, and unstable flow. The momentum conservation equation can also predict the flow patterns formed under different temperature and pressure conditions and the impact of flow pattern changes on permeability properties. Assuming gravity acceleration is represented by h, multiphase flow frictional pressure drop is represented by ods, the following equation provides the expression for the momentum conservation equation:
$\begin{align} & \frac{\partial }{\partial y}\left( {{\vartheta }_{h}}{{\beta }_{h}}{{c}_{h}}S+{{\vartheta }_{1}}{{\beta }_{1}}{{c}_{1}}S+{{\vartheta }_{a}}{{\beta }_{a}}{{v}_{a}}S \right)+\frac{\partial }{\partial x}\left( {{\vartheta }_{h}}{{\beta }_{h}}{{c}_{h}}S+{{\vartheta }_{1}}{{\beta }_{1}}{{c}_{a}}S+{{\vartheta }_{a}}{{\beta }_{a}}{{v}_{a}}S \right) \\ & +\left( {{\vartheta }_{h}}{{\beta }_{h}}+{{\vartheta }_{1}}{{\beta }_{1}}+{{\vartheta }_{a}}{{\beta }_{a}} \right)+hS+\frac{\partial \left( oS \right)}{\partial x}+S\frac{\partial {{o}_{ds}}}{\partial x}=0 \\\end{align}$ (4)
2.3 Air sparging temperature mathematical model
Figure 2. Situation of multiphase flow after model stabilization
This paper considers the air sparging process as a heat exchanger with certain boundary conditions, where the entire circulation process can be divided into several regions radially, each described by different temperature models. Figure 2 shows the situation of multiphase flow after model stabilization. Specifically, it includes the injection zone model, near-well zone model, far-well zone model, transition zone model, and drainage zone model. Assuming density is represented by ϑ, fluid flow rate is represented by w, air sparging frictional power is represented by W, temperature is represented by Y, convective heat transfer coefficient is represented by g, thermal conductivity is represented by j, specific heat capacity is represented by v, axial coordinate is represented by x, radial coordinate is represented by e, in the subscripts, liquid phase is represented by m, gas phase is represented by h, injection zone is represented by o, near-well zone is represented by q, far-well zone is represented by s, transition zone is represented by v, drainage zone is represented by d, well inner wall is represented by ou, well outer wall is represented by op, solid phase inner wall is represented by vu, solid phase outer wall is represented by vp.
The injection zone is where gas is injected into the well, undergoing heat exchange with the surrounding formation. Key parameters include the temperature of the injected gas, flow rate, initial ground temperature, solid thermal conductivity, porosity, and temperature. The temperature equation of the fluid in this area is given by the following equation:
$\begin{align} & {{W}_{o}}-{{\vartheta }_{1}}{{w}_{1}}{{v}_{1}}\frac{\partial {{Y}_{o}}}{\partial x}-2\tau {{e}_{ou}}{{g}_{ou}}\left( {{Y}_{o}}-{{Y}_{q}} \right) ={{\vartheta }_{1}}{{v}_{1}}\tau e_{ou}^{2}\frac{\partial {{Y}_{o}}}{\partial y} \\\end{align}$ (5)
Adjacent to the injection zone is the near-well zone, where intense heat exchange occurs between the fluid and the solid. Key parameters may involve well wall temperature, formation thermal diffusivity, thermal capacity of the gas and formation, and the well's geometric parameters. The temperature equation of the fluid in this area is given by the following equation:
$\begin{align} & {{j}_{q}}\frac{{{\partial }^{2}}{{Y}_{q}}}{\partial {{x}^{2}}}+\frac{2{{e}_{op}}{{g}_{uop}}}{e_{op}^{2}-e_{ou}^{2}}\left( {{Y}_{o}}-{{Y}_{q}} \right) +\frac{2{{e}_{ou}}{{g}_{ou}}}{e_{op}^{2}-e_{ou}^{2}}\left( {{Y}_{o}}-{{Y}_{q}} \right)={{\vartheta }_{q}}{{v}_{q}}\frac{\partial {{Y}_{q}}}{\partial y} \\\end{align}$ (6)
In the far-well zone, the temperature gradually stabilizes as the distance from the well increases, and heat exchange progressively weakens. Key parameters may include the formation's thermal diffusivity, thermal conductivity, the flow rate of the pore fluid, and its temperature. The temperature equation of the fluid in this region is given by the following equation:
$\begin{align} & {{W}_{s}}+\left( {{\vartheta }_{h}}{{v}_{h}}{{w}_{h}}{{\beta }_{h}}+{{\vartheta }_{1}}{{v}_{1}}{{w}_{1}}{{\beta }_{1}}+{{\vartheta }_{a}}{{v}_{a}}{{w}_{a}}{{\beta }_{a}} \right)\frac{\partial {{Y}_{s}}}{\partial x} +2\tau {{e}_{vu}}{{g}_{vu}}\left( {{Y}_{v}}-{{Y}_{s}} \right)+2\tau {{e}_{op}}{{g}_{op}}\left( {{Y}_{q}}-{{Y}_{s}} \right) \\ & =\left( {{\vartheta }_{h}}{{v}_{h}}{{\beta }_{h}}+{{\vartheta }_{1}}{{v}_{1}}{{\beta }_{1}}+{{\vartheta }_{1}}{{v}_{1}}{{\beta }_{a}} \right)\tau \left( e_{vu}^{2}-e_{op}^{2} \right)\frac{\partial {{Y}_{s}}}{\partial y} \\\end{align}$ (7)
The transition zone, located between the near-well and far-well zones, may experience a process where the gas brought from the injection zone gradually equilibrates with the formation water. Key parameters may include the flow rates of gas and liquid, heat exchange coefficients, pressure, and temperature gradients. The temperature equation of the fluid in this region is given by the following equation:
$\begin{align} & {{j}_{v}}\frac{{{\partial }^{2}}{{Y}_{v}}}{\partial {{x}^{2}}}+\frac{2{{e}_{vp}}{{g}_{vp}}}{e_{vp}^{2}-e_{vu}^{2}}\left( {{Y}_{s}}-{{Y}_{v}} \right) +\frac{2{{e}_{vu}}{{g}_{vu}}}{e_{vp}^{2}-e_{vu}^{2}}\left( {{Y}_{s}}-{{Y}_{v}} \right)={{\vartheta }_{v}}{{v}_{v}}\frac{\partial {{Y}_{v}}}{\partial y} \\\end{align}$ (8)
The drainage area is where groundwater is extracted, possibly carrying energy from deep geothermal sources. Key parameters might include the drainage temperature, drainage rate, formation temperature, and heat exchanged with the gas. The temperature equation of the fluid in this region is given by the following equation:
$\frac{{{\partial }^{2}}{{Y}_{d}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{Y}_{d}}}{\partial {{e}^{2}}}+\frac{1}{e}\frac{\partial {{Y}_{d}}}{\partial e}=\frac{{{\vartheta }_{d}}{{V}_{d}}}{{{j}_{d}}}\frac{\partial {{Y}_{d}}}{\partial y}$ (9)
During the air sparging process, calculating the convective heat transfer coefficient of gas-liquid two-phase flow is a complex issue, as it involves interactions between gas and liquid under different flow patterns and their impact on thermal transfer properties. Different flow patterns have distinct flow characteristics and heat exchange mechanisms. To improve the accuracy of on-site calculations, it's usually necessary to choose an empirical formula for gas-liquid two-phase convective heat transfer that suits the specific flow pattern. Initially, the flow pattern of gas-liquid two-phase flow needs to be determined, usually through parameters such as Reynolds number, flow rate, and gas-liquid volume fraction. The flow characteristics under different flow patterns have significantly different impacts on the heat transfer coefficient. Empirical formulas need to be compared with experimental data to verify their accuracy. For specific application scenarios like air sparging, empirical formulas that have been tested under similar conditions and show low error should be chosen. The Dittus-Boelter equation is suitable for turbulent single-phase flow, while the Chisholm equation is used for calculating the pressure drop in two-phase flow, which can be combined with the calculation of heat transfer coefficients. The Lockhart-Martinelli relationship provides a relationship between two-phase flow pressure drop and single-phase flow pressure drop, which can be combined with the calculation of heat transfer coefficients. The Martinelli-Nelson method can be used to estimate the two-phase convective heat transfer coefficient under specific conditions. The final choice should be based on the match between experimental data and site conditions.
3.1 Boundary conditions
When constructing the air sparging transient temperature-pressure mathematical model, temperature and pressure boundary conditions are key factors for the correct representation of the model. These boundary conditions must accurately reflect the actual physical processes to enable the model to reliably predict temperature and pressure changes during transient processes. For temperature boundary conditions, it's necessary to assume a constant value for the temperature boundary condition at the ground surface, set the temperature boundary condition of the injection well to the temperature of the injected gas, assume the model's far-field temperature boundary condition as the initial formation temperature, and consider the thermal exchange between the well wall, surrounding rock, and fluid for the well wall temperature boundary condition. For pressure boundary conditions, the wellhead's pressure boundary condition is usually set to atmospheric pressure, the injection well's pressure boundary condition is determined based on the injected gas pressure, the far-field's pressure boundary condition is set to the formation's original hydrostatic pressure, and the surrounding pressure boundary condition of the drainage well is based on the operational pressure of the drainage well.
3.2 Model discretization
In numerical simulation, different solving methods are used for different equations and problems, mainly based on their advantages and disadvantages in terms of stability, accuracy, and computational efficiency. Since the air sparging transient heat conduction problem often requires handling thermal diffusion processes over long-time scales and simulating physical processes over extended periods, this paper selects the Crank-Nicolson fully implicit scheme for solving the air sparging transient temperature model. The Crank-Nicolson method is unconditionally stable, which is particularly important for transient heat conduction problems, and it is a second-order accurate method in time discretization, meaning it offers higher accuracy in time advancement. For the air sparging transient pressure model, this paper chooses a four-point difference explicit scheme for solving. Explicit methods are generally simpler and faster to compute at each time step, especially when the problem can tolerate smaller time steps. This method is also easier to capture the transient characteristics of pressure and is stable under conditions with smaller spatial and temporal scale variability.
Regarding grid division, for the simulation of air sparging multiphase flow transient temperature-pressure fields, this paper selects fixed axial, fixed time, and variable radial grids. That is, a fixed-spacing grid is allowed in the axial direction where temperature and pressure changes are relatively gentle. To ensure consistency with the typical time scales related to physical processes, a fixed time step is used for simulation. In the radial direction, where there can be larger gradients of temperature and pressure, using a variable radial grid provides higher resolution near the wellbore.
In the fully implicit treatment of the air sparging transient temperature model, it's important to consider the dependence on unknown future temperature values, meaning that at each time step, all dependencies in time and space focus on the value to be calculated at that moment. For spatial variation, a first-order upwind scheme is used to handle changes in the direction of fluid flow. This method helps prevent unrealistic oscillations during calculation and is more stable and accurate when dealing with directional flow issues. For changes over time, a two-point backward difference method is used, allowing consideration of both current and previous temperature values at each time step. This method allows for larger time steps without sacrificing stability, which is particularly important for simulating transient behavior that may occur over longer time scales. For spatial second derivatives, a three-point centered difference method is used, a standard difference format that balances accuracy and stability. It approximates calculations based on the temperature values of nodes and their neighbors, effectively simulating diffusion effects.
In the fully implicit finite difference format, all time terms and spatial derivatives are constructed based on the estimated values at future moments, meaning that each time step's computation requires solving a linear or nonlinear system of equations involving future temperatures of all nodes. The injection zone simulation focuses on temperature changes near the injection well because the injected fluid may have different temperatures, so special attention is needed for the thermal changes brought by fluid injection. The fully implicit format predicts the temperature for the next time step, including temperature changes caused by the injected fluid and thermal conduction from the surrounding rock, with its finite difference format as follows:
$\begin{align} & \left( {{W}_{w}} \right)_{k}^{b}-\left( {{\vartheta }_{1}}{{w}_{1}}{{v}_{1}} \right)_{k}^{b}\frac{Y_{1,k}^{b}-Y_{1,k-1}^{b}}{\Delta x} -2\tau {{e}_{ou}}\left( {{g}_{ou}} \right)_{k}^{b}\left( Y_{1,k}^{b}-Y_{2,k}^{b} \right) =\left( {{\vartheta }_{1}}{{v}_{1}} \right)_{k}^{b}\tau e_{ou}^{2}\frac{Y_{1,k}^{b}-Y_{1,k}^{b-1}}{\Delta y} \\\end{align}$ (10)
The temperature changes in the near-well zone are influenced by the dynamics of the fluid inside the well and the thermal transfer of the surrounding rock. The fully implicit method calculating the temperature distribution in this area must consider well wall conditions, the thermal properties of the fluid, and potential thermal convection. Its finite difference format is:
$\begin{align} & {{j}_{q}}\frac{Y_{2,k+1}^{b}-2Y_{2,k}^{b}+Y_{2,k-1}^{b}}{\Delta {{x}^{2}}} +\frac{2{{e}_{op}}\left( {{g}_{op}} \right)_{b}^{k}}{e_{op}^{2}-e_{ou}^{2}}\left( Y_{3,k}^{b}-Y_{2,k}^{b} \right) \\ & +\frac{2{{e}_{ou}}\left( {{g}_{ou}} \right)_{b}^{k}}{e_{op}^{2}-e_{ou}^{2}}\left( Y_{1,k}^{b}-Y_{2,k}^{b} \right)={{\vartheta }_{q}}{{v}_{q}}\frac{Y_{2,k}^{b}-Y_{2,k}^{b-1}}{\Delta y} \\\end{align}$ (11)
In areas far from the well, changes in the temperature field are primarily controlled by thermal conduction of the underground rocks. In the fully implicit format, temperature changes in the far-well zone can consider thermal transfer effects from a broader area, and time-related changes are implicitly predicted. Its finite difference format is:
$\begin{align} & \left( {{W}_{s}} \right)_{k}^{b}+\left( Z \right)_{k}^{b}\frac{Y_{3,k+1}^{b}-Y_{3,k}^{b}}{\Delta x} +2\tau {{e}_{vu}}\left( {{g}_{vu}} \right)_{k}^{b}\left( Y_{4,k}^{b}-Y_{3,k}^{b} \right) +2\tau {{e}_{op}}\left( {{g}_{op}} \right)_{k}^{b}\left( Y_{2,k}^{b}-Y_{3,k}^{b} \right) =\left( {{Z}_{2}} \right)_{k}^{b}\tau \left( e_{vu}^{2}-e_{op}^{2} \right)\frac{Y_{3,k}^{b}-Y_{3,k}^{b-1}}{\Delta y} \\\end{align}$ (12)
The transition area can be seen as an intermediate zone of thermal conduction and convection. In this area, the fully implicit method needs to handle the interaction of these two mechanisms, ensuring the continuous transition of temperature and accurate description of physical phenomena. Its finite difference format is:
$\begin{align} & {{j}_{v}}\frac{Y_{4,k+1}^{b}-2Y_{4,k}^{b}+Y_{4,k-1}^{b}}{\Delta {{x}^{2}}} +\frac{2{{e}_{vp}}\left( {{g}_{vp}} \right)_{k}^{b}}{e_{vo}^{2}-e_{vu}^{2}}\left( Y_{5,k}^{b}-Y_{4,k}^{b} \right) +\frac{2{{e}_{vu}}\left( {{g}_{vu}} \right)_{k}^{b}}{e_{vo}^{2}-e_{vu}^{2}}\left( Y_{3,k}^{b}-Y_{4,k}^{b} \right) ={{\vartheta }_{v}}{{v}_{v}}\frac{Y_{4,k}^{b}-Y_{4,k}^{b-1}}{\Delta y} \\\end{align}$ (13)
Temperature changes in the drainage area may be affected by extraction operations and fluid replacement. The fully implicit format here aims to predict thermal convection effects brought by fluid movement and thermal exchange between rock and fluid. Its finite difference format is:
$\begin{align} & \frac{Y_{u,k+1}^{b}-2Y_{u,k}^{b}+Y_{u,k-1}^{b}}{\Delta {{x}^{2}}}+\frac{Y_{u,k+1}^{b}-2Y_{u,k}^{b}+Y_{u-1,k}^{b}}{\Delta {{e}^{2}}} +\frac{1}{{{e}_{u}}}\frac{Y_{u,k+1}^{b}-Y_{u,k}^{b}}{\Delta e}=\frac{{{\vartheta }_{d}}{{v}_{d}}}{{{j}_{d}}}\frac{Y_{u,k+1}^{b}-Y_{u,k}^{b-1}}{\Delta y} \\\end{align}$ (14)
where,
${{Z}_{1}}={{\vartheta }_{h}}{{v}_{h}}{{w}_{h}}{{\beta }_{h}}+{{\vartheta }_{1}}{{v}_{1}}{{w}_{1}}{{\beta }_{1}}+{{\vartheta }_{a}}{{v}_{a}}{{w}_{a}}{{\beta }_{a}}$ (15)
${{Z}_{2}}={{\vartheta }_{h}}{{v}_{h}}{{\beta }_{h}}+{{\vartheta }_{1}}{{v}_{1}}{{\beta }_{1}}+{{\vartheta }_{a}}{{v}_{a}}{{\beta }_{a}}$ (16)
The explicit treatment of the air sparging transient pressure model means directly using currently known information to predict the state at the next time step, without the need to solve a complex system of equations. For the first-order spatial derivatives, a first-order upwind scheme is chosen, which selects difference points based on the direction of material flow, ensuring material conservation and the physical reasonableness of the numerical solution. In the air sparging model, this means that the pressure gradient in the direction of fluid flow will be calculated based on the upwind node of the flow direction. For the first-order time derivatives, a four-point centered difference is used, which calculates using information from two points in time past, one point in time past, the current point, and a future time point.
In explicit methods, time derivatives are approximated using known values from several previous time steps, while spatial derivatives are approximated using spatial node values at the current time step. Due to the conditional restrictions of explicit methods, a sufficiently-small time step must be chosen to ensure the stability of the numerical solution. In this case, the finite difference format of various parts of the pressure model predicts the state at the next time step based on current and past information. When dealing with the pressure models of liquid, gas, and solid phases, the construction principle of the explicit finite difference format is based on the laws of mass conservation and momentum conservation. The flow of each phase can be described by fluid dynamics equations, and in multiphase flow, these equations often become more complex due to the need to consider interactions and energy exchanges between different phases. In the explicit finite difference method, these continuous equations are converted into discrete equations to estimate future state values at grid points.
For the liquid phase, the model typically includes Darcy's law, which describes the flow of liquid through porous media. The explicit finite difference format combines Darcy's law with the continuity equation of the fluid. The following equation provides the finite difference format expression for the liquid phase pressure model:
$\left( {{\vartheta }_{1}}{{\beta }_{1}}{{c}_{1}} \right)_{k}^{b}-J\left( {{\vartheta }_{1}}{{\beta }_{1}}{{c}_{1}} \right)_{k-1}^{b}=\frac{\Delta x}{2\Delta y}\left[ J\left( {{\vartheta }_{1}}{{\beta }_{1}} \right)_{k-1}^{b-1}+ \right.\left. \left( {{\vartheta }_{1}}{{\beta }_{1}} \right)_{k}^{b-1}-J\left( {{\vartheta }_{1}}{{\beta }_{1}} \right)_{k-1}^{b}-\left( {{\vartheta }_{1}}{{\beta }_{1}} \right)_{k}^{b} \right]$ (17)
Gas flow can be simulated using a similar method, but the compressibility of gas means that its density can significantly vary, especially with large changes in pressure. Therefore, the gas phase pressure model must include the state equation, and it may also need to consider the impact of temperature. The following equation provides the finite difference format expression for the gas phase pressure model:
$\left( {{\vartheta }_{h}}{{\beta }_{h}}{{c}_{h}} \right)_{k}^{b}-J\left( {{\vartheta }_{h}}{{\beta }_{h}}{{c}_{h}} \right)_{k-1}^{b}=\frac{\Delta x}{2\Delta y}\left[ J\left( {{\vartheta }_{h}}{{\beta }_{h}} \right)_{k-1}^{b-1}+ \right.\left. \left( {{\vartheta }_{h}}{{\beta }_{h}} \right)_{k}^{b-1}-J\left( {{\vartheta }_{h}}{{\beta }_{h}} \right)_{k-1}^{b}-\left( {{\vartheta }_{h}}{{\beta }_{h}} \right)_{k}^{b} \right]$ (18)
Solid phase flow is usually described by geomechanical equations, which lean more towards elasticity theory and plasticity theory than fluid dynamics. The pressure model for the solid phase involves rock deformation and the stress-strain relationship. The following equation provides the finite difference format expression for the solid phase pressure model:
$\begin{align} & \left( {{\vartheta }_{a}}{{\beta }_{a}}{{c}_{a}} \right)_{k}^{b}-J\left( {{\vartheta }_{a}}{{\beta }_{a}}{{c}_{a}} \right)_{k-1}^{b}=\frac{\Delta x}{2\Delta y}\left[ J\left( {{\vartheta }_{a}}{{\beta }_{a}} \right)_{k-1}^{b-1}+ \right.\left. \left( {{\vartheta }_{a}}{{\beta }_{a}} \right)_{k}^{b-1}-J\left( {{\vartheta }_{a}}{{\beta }_{a}} \right)_{k-1}^{b}-\left( {{\vartheta }_{a}}{{\beta }_{a}} \right)_{k}^{b} \right] \\ & +\frac{\Delta z}{2{{S}_{k}}}\left[ \left( {{w}_{a}} \right)_{k-1}^{b}+\left( {{w}_{a}} \right)_{k}^{b} \right] \\\end{align}$ (19)
To construct the finite difference format of the momentum conservation equation, the continuous momentum conservation equation needs to be discretized in space and time. The computational domain is divided into a finite number of grid cells, and the velocity of the fluid and other relevant physical quantities are defined at each grid point or at the center of each grid cell. Spatial derivatives are calculated based on the values at neighboring grid points. Time derivatives can be approximated through the changes in velocity at different time steps. Inserting the above spatial and temporal discretized derivatives into the continuous momentum conservation equation yields the difference equation at each grid point. The following provides the finite difference format expression for the momentum conservation equation:
$\begin{align} & \frac{\Delta x}{2\Delta y}\left( JT_{1k-1}^{b-1}+T_{1k}^{b}-JT_{1k-1}^{b}-T_{1k}^{b} \right)+JT_{2k-1}^{b}-T_{2k}^{b}-\frac{h\Delta x}{2}\left( JT_{3k-1}^{b}+T_{3k}^{b} \right) \\ & -\frac{\Delta x}{2}\left[ \left( J\frac{\partial {{o}_{ds}}}{\partial x} \right)_{k-1}^{b}+\left( \frac{\partial {{o}_{ds}}}{\partial x} \right)_{k}^{b} \right] \\\end{align}$ (20)
where,
$J={{{S}_{k-1}}}/{{{S}_{k}}}\;$ (21)
${{T}_{1}}={{\vartheta }_{h}}{{\beta }_{h}}{{c}_{h}}+{{\vartheta }_{1}}{{\beta }_{1}}{{c}_{1}}+{{\vartheta }_{a}}{{\beta }_{a}}{{c}_{a}}$ (22)
${{T}_{2}}={{\vartheta }_{h}}{{\beta }_{h}}c_{h}^{2}+{{\vartheta }_{1}}{{\beta }_{1}}c_{1}^{2}+{{\vartheta }_{a}}{{\beta }_{a}}c_{a}^{2}$ (23)
${{T}_{3}}={{\vartheta }_{h}}{{\beta }_{h}}+{{\vartheta }_{1}}{{\beta }_{1}}+{{\vartheta }_{a}}{{\beta }_{a}}$ (24)
3.3 Solution algorithm
Figure 3. Steps for solving the coupled model
The solution to the air sparging multiphase flow transient temperature-pressure field coupled model is achieved through an iterative process, as detailed in Figure 3 and described below:
(1) In the process of establishing the numerical model, the entire physical domain of air sparging is first divided into small control volumes according to certain rules, with each grid cell representing a certain physical space and characteristics. Time is also discretized into multiple time steps, which can be fixed or dynamically adjusted as needed. In terms of spatial division, the grid must be detailed enough to capture gradients in pressure and temperature while also considering computational efficiency. In terms of time, the selection of time steps needs to balance computational accuracy and stability.
(2) The thermal properties of different fluids typically depend on pressure and temperature. This step requires calculating the thermal properties of each component under the current state based on state equations and empirical formulas. The calculation of volume fractions needs to be determined based on the phase equilibrium relationship of the fluids, which often involves complex multiphase flow and phase transition physical processes.
(3) The transient temperature equation typically includes the energy conservation equation, which accounts for conduction, convection, thermal radiation, and heat exchange between fluids. By substituting the current temperature Y0k, pressure o0k, thermal properties, and volume fractions of each grid cell into the transient temperature equation and using an appropriate numerical method, the temperature distribution Yk for the next time step y=y0+∆y is calculated.
(4) The transient pressure model describes the flow of fluids through porous media, involving Darcy's law and the fluid continuity equation. The pressure field for the next time step is calculated using the temperature distribution and the thermal properties of the fluids at the time step y=y0+∆y.
(5) For each grid point in the model, local fluid parameters such as flow velocity, pressure, temperature, and water content are calculated. This step is key to achieving temperature and pressure coupling and needs to consider the interactions between fluid motion and heat transfer.
(6) A predefined convergence criterion is used to determine if the pressure field obtained from the current iteration meets the accuracy requirements. If not, adjustments may be needed, including grid division, time step size, boundary conditions, etc., and the process returns to Step 4 to recalculate the pressure field until it satisfies the accuracy requirements for air sparging pressure ok.
(7) Similar to the pressure field, the temperature field also needs to meet certain convergence criteria, i.e., ||ok-o0k||<γo and ||Yk-Y0k||<γY. If these criteria are not met, it may be necessary to return to Step 3, possibly adjusting related parameters in the temperature model, such as the methods for calculating thermal properties and thermal diffusivity, and recalculate the temperature field.
(8) Once the pressure and temperature fields at y=y0+∆y meet the accuracy requirements, these values are used as the initial conditions for the next time step. Then, the iterative process begins again from Step 1 until the computations for the entire simulation period are completed.
Figure 4 shows the layout of groundwater monitoring points during the experiment.
Figure 4. Layout of groundwater monitoring points
Table 1. Percentage drop in groundwater level downstream of the air injection point
|
Distance |
20kPa |
25kPa |
||||||
|
0.4m |
0.8m |
1.2m |
1.6m |
0.4m |
0.8m |
1.2m |
1.6m |
|
|
4.5 |
10.2% |
16.9% |
20.8% |
23.8% |
16.2% |
24.8% |
26.8% |
31.2% |
|
5 |
7.5% |
14.8% |
21.6% |
21.5% |
14.5% |
23.1% |
24.7% |
28.9% |
|
5.5 |
6.7% |
13.5% |
18.9% |
22.3% |
13.6% |
24.9% |
24.8% |
28.2% |
|
6 |
5.4% |
13.4% |
19.2% |
22.7% |
12.8% |
24.5% |
25.1% |
27.3% |
|
6.5 |
5.2% |
12.1% |
18.5% |
22.3% |
12.4% |
23.8% |
25.3% |
26.3% |
|
7 |
5.1% |
11.4% |
16.8% |
18.9% |
8.8% |
23.8% |
23.9% |
25.2% |
|
7.5 |
4.3% |
9.5% |
16.2% |
17.8% |
8.8% |
23.4% |
22.8% |
23.4% |
|
Average Drop |
6.3% |
13.1% |
18.9% |
21.3% |
12.4% |
24.0% |
24.8% |
27.2% |
Table 1 shows the percentage drop in groundwater level downstream of the air injection point under different pressures (20kPa and 25kPa) at various distances (from 4.5m to 7.5m) and downstream distances (from 0.4m to 1.6m). The table reveals that at both pressure levels, the drop in groundwater level gradually decreases as the distance from the injection point increases, indicating that the pressure effect diminishes with distance. For a given pressure, the drop in groundwater level is greater closer to the injection point and decreases as the distance increases, indicating that the pressure transmission effect attenuates with distance. At a pressure of 25kPa, the drop in groundwater level is generally higher than at 20kPa, demonstrating that the drop in groundwater level is positively correlated with the applied pressure—the higher the pressure, the greater the drop. For 20kPa pressure, the maximum drop occurs at a distance of 1.6m; while at 25kPa pressure, the maximum drop appears at different distances but mostly at 1.2m or 1.6m. This may suggest that there is an optimal distance where the pressure-induced water level drop is maximized. At 25kPa pressure, the average drop (13.1% to 26.5%) is significantly higher than at 20kPa pressure (6.3% to 22.6%), further confirming the direct impact of pressure magnitude on groundwater level drop. It can be concluded that by quantifying the water level drop under different pressures and distances, the model effectively describes the propagation of pressure in the groundwater system and its impact on water levels. The proposed method can provide accurate predictions of groundwater level drops under different pressures and distances, aiding in better decision-making in water resource management and engineering practices.
Figure 5. Change in average drop of downstream water level (air sparging water cutting effect) with air injection pressure and time
Data provided in Figure 5 shows how the average drop in downstream water level changes with time under different air injection pressures. This data can be used to analyze the air sparging water cutting effect under various pressure and time conditions. It is observed that in the initial period (especially at lower pressures), the water level drop remains at 0% throughout the measurement time, indicating that the water cutting effect is not significant at this pressure. As pressure increases (from 15kPa to 40kPa), the average drop in water level gradually increases over time, showing a more significant water cutting effect. At this point, the water level drop also becomes less variable, which may be due to the groundwater system needing time to respond to the increase in air injection pressure. As time continues to increase, the water level drop becomes more pronounced, especially at higher pressures. At higher pressures (e.g., 35kPa and 40kPa), the water level drop tends to saturate after a period, meaning that increasing time no longer significantly increases the water level drop. At 40kPa pressure, the water level drop slightly decreases after 64 hours, possibly indicating that the groundwater system has reached a dynamic equilibrium state or due to other factors such as soil pore pressure redistribution or changes in permeability characteristics. It can be concluded that the downstream water level drop increases with the increase in air injection pressure, indicating that pressure is a key factor affecting the water cutting effect. The model constructed can capture the impact of pressure and time on the water cutting effect, showing effective prediction for changes in groundwater level.
Table 2. Groundwater level drop under different air injection pressures
|
Distance from the Top Boundary of the Slope |
15kPa |
20kPa |
25kPa |
30kPa |
|
4.5 |
7.1% |
11.6% |
12.2% |
- |
|
5 |
6.7% |
11.7% |
15.3% |
18.6% |
|
5.5 |
5.2% |
8.1% |
10.1% |
15.4% |
|
6 |
4.8% |
7.1% |
9.7% |
13.8% |
|
6.5 |
3.6% |
5.3% |
9.1% |
12.7% |
|
7 |
3.7% |
5.8% |
8.9% |
11.3% |
|
7.5 |
3.2% |
5.8% |
8.7% |
11.3% |
|
Average Drop |
4.9% |
7.9% |
10.6% |
11.9% |
1) At larger Marangoni number
2) At smaller Marangoni number
Figure 6. Temperature distribution under different Marangoni numbers and pressures
Table 2 displays the relationship between the drop in groundwater level under different air injection pressures and the distance from the top boundary of the slope. By examining the groundwater level drop under various distances and pressure combinations, the effectiveness of the proposed air sparging transient temperature-pressure field coupled model can be assessed. The table indicates that within a certain range, as the air injection pressure increases (from 15kPa to 30kPa), the average drop in groundwater level also increases. This trend remains consistent across the entire measurement range, showing that pressure is an important factor affecting the drop in water level. For a given pressure, the water level drop decreases with increasing distance from the top boundary of the slope. This may mean that areas closer to the boundary are more sensitive to air injection pressure. At all pressure levels, as the distance increases, the trend of decreasing water level drop becomes more gradual, indicating that the distance effect may tend to saturate after a certain distance. At lower pressures (15kPa and 20kPa), the water level drop decreases more rapidly with increasing distance. At higher pressures (25kPa and 30kPa), the trend of decreasing water level drop is relatively gentle, indicating that the water cutting effect is less sensitive to distance under high pressure. It can be concluded that the water level drop decreases with increasing distance from the top boundary of the slope, but this decrease has a saturating distance effect. The proposed method can effectively predict water level drop under different pressures and distances, which has obvious practical value for groundwater interception and drainage management in actual engineering.
The Marangoni effect, caused by the temperature gradient of liquid surface tension, is especially important in multiphase flow as it can significantly influence the behavior of liquid interfaces. From Figure 6, it can be observed that under smaller and larger Marangoni numbers, the temperature distribution changes accordingly with increasing pressure. At smaller Marangoni numbers, with increasing radial coordinate, the temperature distribution shows a trend of first increasing and then decreasing under all pressure conditions, reaching a peak at radial coordinate 0.7. Increasing pressure causes a slight increase in temperature at various radial coordinate positions, but this increase is very minimal. At radial coordinate 1, the temperature distribution under different pressures converges to nearly the same value. Under larger Marangoni numbers, the peak of the temperature distribution is higher, and it increases more significantly with the radial coordinate. The temperature difference under different pressures is greater than in cases with smaller Marangoni numbers, especially at intermediate radial coordinate positions (0.5 to 0.8). Despite differences in temperature distribution under different pressures, at radial coordinate 1, the temperature distributions under different pressures also tend to the same value. It can be concluded that for both Marangoni number scenarios, the increase in pressure has some impact on the temperature distribution, but this impact is relatively small, especially when the Marangoni number is smaller. Since the model can capture the changes between temperature and pressure, this indicates that the constructed transient temperature-pressure field coupled model can effectively simulate the dynamic changes in temperature in multiphase flow, particularly under different Marangoni numbers and pressure conditions.
Further, based on Figure 7, a comparative analysis of two types of critical heat flux density under different Marangoni numbers and ambient temperatures is conducted. The figure indicates that at lower pressures, the critical heat flux density decreases with an increase in ambient temperature, likely due to the saturation temperature of the liquid approaching as the ambient temperature rises. For different Marangoni numbers (1200, 2400, 3600), the critical heat flux density generally shows a downward trend as the number increases, suggesting that higher Marangoni numbers might reduce the thermal stability of the system. At medium pressures, the trend of change in critical heat flux density is similar to that at lower pressures, gradually decreasing with increasing ambient temperature. Under medium pressure, for the same ambient temperature and Marangoni number, different critical values (Critical 1 and Critical 2) show varying heat flux densities, due to other system conditions or parameters (such as fluid velocity, liquid properties, etc.) changing under different critical states. As pressure increases, for Marangoni number 1200's two critical points, heat flux density slightly increases at lower ambient temperatures but is lower under higher pressure than lower pressure at higher ambient temperatures. Similar to medium pressure, high pressure also shows different heat flux densities for critical 1 and 2, further indicating that critical heat flux density is affected by multiple factors. It can be concluded that changes in pressure have a complex impact on heat flux density, depending on specific operational conditions and ambient temperature. The model proposed in this paper can capture trends of heat flux density with changes in ambient temperature and pressure, which is crucial for ensuring the safety and efficiency of heat exchange systems.
1) At low pressure
2) At medium pressure
3) At high pressure
Figure 7. Two types of critical heat flux density under different Marangoni numbers and ambient temperatures
This research work represents a significant advancement in the field of air sparging remediation technology, focusing on the impact of temperature gradients on multiphase flow and proposing and validating a fully transient temperature-pressure field coupled model. The paper successfully constructed a fully transient temperature-pressure field coupled model, innovatively considering the complexity of changes in the physical parameters of groundwater in the cyclic slope body with time, temperature, and pressure, which is not commonly seen in previous research. Using this model, the paper simulated and analyzed multiphase flow problems, providing an accurate theoretical basis for exploring the impact of temperature gradients on multiphase flow in air sparging remediation technology. A novel model iterative solution algorithm was also proposed to address computational challenges caused by complex boundary conditions and nonlinear factors, enhancing the efficiency and accuracy of the solution process.
Experimental data show that as air injection pressure increases and the duration of air injection extends, the average drop in downstream water level gradually increases, indicating that the model can effectively capture the impact of pressure and time variables on the water cutting effect. Data also reveals that the water level drop decreases with increasing distance from the top boundary of the slope, further validating that the model can accurately reflect complex groundwater flow phenomena.
The fully transient temperature-pressure field coupled model proposed in this paper effectively integrated the theory of multiphase flow with the impact of temperature gradients, offering new perspectives and methods for understanding and predicting multiphase flow phenomena in air sparging remediation technology. At the same time, the novel iterative solution algorithm provides a powerful computational tool for handling complex groundwater flow problems. The research results prove the effectiveness of the model and method, not only having significant theoretical value for scientific research but also aiding in guiding groundwater management and control measures in practical engineering applications, especially in multiphase flow scenarios where temperature effects need to be considered.
This paper was supported by the Central Leading Local Science and Technology Development Fund Program (Grant No.: 236Z5405G) and the Natural Science Foundation of Hebei Province (Grant No.: E2021210092).
[1] Xie, W., Shang, Y., Lü, Q., Jiang, H., Wei, Z. (2018). Experimental study of groundwater level variation in soil slope using air-injection method. Geotechnique Letters, 8(2): 144-148. https://doi.org/10.1680/jgele.18.00005
[2] Xie, W., Shang, Y., Wu, G., Wei, Z. (2018). Investigation of the formation process of a low-permeability unsaturated zone by air injection method in a slope. Engineering Geology, 245: 10-19. https://doi.org/10.1016/j.enggeo.2018.08.005
[3] Qiao, M., Xing, Z., Fu, J., Liu, J. (2023). Multiphase flow physics of room temperature liquid metals and its applications. Science China Technological Sciences, 66(6): 1483-1510. https://doi.org/10.1007/s11431-022-2295-4
[4] Melikov, H.K., Ismayilov, S.Z., Suleymanov, A.A., Mammadli, N.F. (2022). Diagnosing multiphase flow regime in multilayered reservoir by distributed temperature sensor measurements. SOCAR Proceedings, 1: 47-55. https://doi.org/10.2118/189034-MS
[5] Liu, Z., Li, B., Xiao, L., Gan, Y. (2022). Modeling progress of high-temperature melt multiphase flow in continuous casting mold. Acta Metallurgica Sinica, 58(10): 1236-1252. https://doi.org/10.11900/0412.1961.2022.00175
[6] Li, X.B., Hao, X.Y., Zhang, H.N., Zhang, W.H., Li, F.C. (2023). Review on multi-parameter simultaneous measurement techniques for multiphase flow–Part A: Velocity and temperature/pressure. Measurement, 223: 113710. https://doi.org/10.1016/j.measurement.2023.113710
[7] Wei, N., Pei, J., Xue, J., Jiang, L., Li, H., Sun, W., Zhang, Y. (2022). Modeling of multiphase flow during underbalanced drilling considering velocity, temperature and pressure. Energy Reports, 8: 2574-2587. https://doi.org/10.1016/j.egyr.2022.01.160
[8] Park, H.D., Paul, M., Hammond, G.E., Valocchi, A.J. (2022). Newton trust-region methods with primary variable switching for simulating high temperature multiphase porous media flow. Advances in Water Resources, 168: 104285. https://doi.org/10.1016/j.advwatres.2022.104285
[9] Wei, C., Tabjula, J.L., Sharma, J., Chen, Y. (2023). The modeling of two-way coupled transient multiphase flow and heat transfer during gas influx management using fiber optic distributed temperature sensing measurements. International Journal of Heat and Mass Transfer, 214: 124447. https://doi.org/10.1016/j.ijheatmasstransfer.2023.124447
[10] George, B., Pierson, L., Barrows, R. (2020). Research in using polyurethane foam to mitigate rock bolt installation difficulties. In 54th U.S. Rock Mechanics/Geomechanics Symposium, Golden, Colorado, USA, pp. 1-5
[11] George, B., Pierson, L., Barrows, R. (2018). Research in using Polyurethane Foam to Mitigate Rock Bolt Installation Difficulties. Transportation Research Record, 2672(52): 307-315. https://doi.org/10.1177/03611981187903
[12] Turkaya, S., Toussaint, R., Eriksen, F. K., Zecevic, M., Daniel, G., Flekkøy, E.G., Måløy, K.J. (2015). Bridging aero-fracture evolution with the characteristics of the acoustic emissions in a porous medium. Frontiers in Physics, 3: 70. https://doi.org/10.3389/fphy.2015.00070
[13] Liu, X., Luo, O., Huang, X., Zhu, H. (2022). Optimization algorithm of building drainage pipe system based on isotope tracking method. In International Conference on Advances in Materials and Manufacturing, Bhopal, India, pp. 861-871. https://doi.org/10.1007/978-981-99-2921-4_77
[14] Zhou, Y., Li, Y., Yan, Z., Wang, H., Chen, H., Zhao, S., Acharya, K. (2023). Microplastics discharged from urban drainage system: Prominent contribution of sewer overflow pollution. Water Research, 236: 119976. https://doi.org/10.1016/j.watres.2023.119976
[15] Zhu, D.F., Cao, H., Luo, G.Y., Pan, H., Mei, J.L. (2018). Application of interception and drainage anti-floating system in treatment of uplift accidents. Yantu Gongcheng Xuebao/Chinese Journal of Geotechnical Engineering, 40(9): 1746-1752. https://doi.org/10.11779/CJGE201809023
[16] Zhao, C., Wang, M., Cao, H., Qu, G., Wang, Y., Guo, Y. (2022). Water control principle of lateral water cutoff curtain and numerical analysis of its water interception effect in Yuanbaoshan open-pit coal mine. Meitiandizhi Yu Kantan/Coal Geology and Exploration, 50(7): 10-17. https://doi.org/10.12363/issn.1001-1986.21.12.0767
[17] Yu, J., Yu, S., Guo, X. (2015). Flow and heat transfer in high-pressure bypass valve. Journal of Central South University (Science and Technology), 46(12): 4693-4699. https://doi.org/10.11817/j.issn.1672-7207.2015.12.041
[18] Yang, H., Li, J., Zhang, G., Zhang, H., Guo, B., Chen, W. (2022). Wellbore multiphase flow behaviors during gas invasion in deepwater downhole dual-gradient drilling based on oil-based drilling fluid. Energy Reports, 8: 2843-2858. https://doi.org/10.1016/j.egyr.2022.01.244
[19] Liu, Z., Li, B., Xiao, L., Gan, Y. (2022). Modeling progress of high-temperature melt multiphase flow in continuous casting mold. Acta Metallurgica Sinica, 58(10): 1236-1252. https://doi.org/10.11900/0412.1961.2022.00175
[20] Lu, W., Yang, Z., Duan, Y. (2023). Numerical investigation on bubble growth characteristics of flow boiling in a microchannel. Journal of Engineering Thermophysics, 44(6): 1612-1616.
[21] Pan, Z., Trusler, J.M. (2022). Interfacial tensions of systems comprising N2, 7 mass% KI (aq): Decane and iododecane at elevated pressures and temperatures. Fluid Phase Equilibria, 556: 113364. https://doi.org/10.1016/j.fluid.2021.113364
[22] Si, Y., Li, S., Chen, Q., Meng, B., Wang, C., Tian, B. (2023). Heat transfer effects on multiphase Richtmyer-Meshkov instability of dense gas-particle flow. Physics of Fluids, 35(5): 053339. https://doi.org/10.1063/5.0149563
[23] Ji, D., Xu, J., Lyu, X., Li, Z., Zhan, J. (2022). Numerical modeling of the steam chamber ramp-up phase in steam-assisted gravity drainage. Energies, 15(8): 2933. https://doi.org/10.3390/en15082933
[24] Jiang, Z., Yang, G., Shen, Q., Li, S., Liao, J., Wang, H., Zhang, G., Li, Z., Liu, Z. (2023). Pore scale study of ice melting and multiphase flow in gas diffusion layer with porosity gradient in PEMFC. Applied Thermal Engineering, 231: 120951. https://doi.org/10.1016/j.applthermaleng.2023.120951
