Modeling of flow of two-phase mixture in curved channel pipeline
V. E. Lyalin^{1} , A. N. Krasnov^{2}
^{1}Kalashnikov Izhevsk State Technical University, Izhevsk, Russia
^{2}Ufa State Petroleum Technological University, Ufa, Russia
^{1}Corresponding author
Vibroengineering PROCEDIA, Vol. 7, 2016, p. 188-194.
Received 5 August 2016; accepted 8 August 2016; published 31 August 2016
This article describes the three-dimensional modeling of steady motion of gas and condensed phase particles of varying severity in a curvilinear orthogonal coordinate system gas channels. A numerical method for solving the equations that describe the motion of the particles, based on the difference scheme and the application of the two-level iterative process. Built multivariate mathematical model of viscous steady flow of gas and the equilibrium of hydrate particles in axisymmetric taking into account transport and energy dissipation. Obtained pressure field in the pipe by solving the model equations using SIMPLE method. Numerical calculations of the velocity field of the gas and the dispersed phase, the trajectory of the mass flow deposited on particles of different sizes with the wall during the wet gas in the curved pipe. The calculation of particle trajectories in different kinds of process equipment show at the option of either abrasive wear parts located in a stream, or an intense build-up and accumulation of hydrates. The results can be applied in the analysis of impact of solid particles of various sizes on the elements of gas pipelines.
Keywords: model of gas flow, curved pipe channels, dispersed gas particles.
1. Introduction
According to studies hydrocarbon hydrates are white crystalline solids, like snow, and at seal similar to ice [1]. The condensed phase in the flow of natural gas hydrates in addition may contain other solid impurities (slag, sand, etc.) Mechanical impact at high velocity a mixture of gas and particles can be exposed structural elements and stop measuring gas equipment. Assessment of the impact can be made on the basis of solving the equations of motion of a two-phase mixture in the pipeline elements.
2. The equations of two-phase flow of gas in the curved pipes
The system of equations describing the steady flow of viscous heat-conducting gas in an arbitrary coordinate system is as follows:
where $\rho $ – density of the gas; $\mathbf{V}$ – The velocity vector; $e$ – Specific internal energy; $\mathbf{P}$ – Pressure tensor related to pressure $p$ and viscous stress tensor $\mathbf{S}$ relation $\mathbf{P}=p\mathbf{I}+\mathbf{S}$. The viscous stress tensor is defined as:
where $\mathbf{I}$ – the unit tensor, $\mu $ – dynamic viscosity coefficient. Tensor components $\mathbf{S}$ are shown in many literatures, for example in [2]. Heat flow $\mathbf{q}$ is determined by the Fourier law: $\mathbf{q}=-\lambda \nabla T$, where $\lambda $ – coefficient of thermal conductivity, $T$ – temperature. For the system of Eq. (2), you must add the equation of state of the gas in the form $p=\rho RT$, where $R$ – the gas constant. In the calculation of flows with low speeds (up to Mach numbers $M<\text{0.3}$), it is advisable to consider for incompressible. In this case, in the Eq. (1) $\rho =const$ and the equation for the energy can not be seen due to the small change in the flow temperature.
Together with the Eq. (1) it is necessary to solve the equation for the turbulent viscosity. It was used model of turbulent viscosity [3].
The equation for the turbulent viscosity is as follows:
where ${\nu}_{T}$ – turbulent viscosity; $\mu =\rho ({\nu}_{T}+{\nu}_{0})$, ${\mathrm{\Gamma}}_{T}$ – A term describing the generation and dissipation of turbulent viscosity [3]; ${\nu}_{0}$ – Molecular viscosity.
To calculate the flow in curved pipes need to solve the problem of flow in a three-dimensional setting. Solving such problems is difficult enough. To evaluate the mechanical interaction of a condensed phase with the walls will be considered two-dimensional setting. The system of Eqs. (1), (2) can be written in a flat or cylindrical coordinates with axial symmetry ($x$ – longitudinal, $r$ – transverse coordinates). For the curved regions of the flow plane $x$, $r$ applied coordinate transformation taking curvilinear computational domain into a rectangle. The difference mesh is constructed using a complex boundary element method [4].
To write equations used curvilinear orthogonal coordinate $\xi $, $\eta .$ System variables $\xi $, $\eta $ associated with Cartesian coordinates $x$, $r$ as follows: $\xi =\xi (x,r)$, $\eta =\eta (x,r)$. The components of the metric tensor are equal ${g}_{11}={g}_{22}=D$. In the case of orthogonal coordinates ${x}_{\xi}={r}_{\eta}$, ${x}_{\eta}=-{r}_{\xi}$, $D={x}_{\xi}^{2}+{r}_{\xi}^{2}$.
The equations of continuity, momentum and energy in the coordinate system $\xi $, $\eta $ are of the form [5]:
$(J{r}^{\gamma}\rho uu{)}_{\xi}+(J{r}^{\gamma}\rho vu{)}_{\eta}+{r}^{\gamma}\rho vu{J}_{\eta}-{r}^{\gamma}\rho vv{J}_{\xi}=-J{r}^{\gamma}{p}_{\xi}+Jr{A}_{1},$
$(J{r}^{\gamma}\rho uv{)}_{\xi}+(J{r}^{\gamma}\rho vv{)}_{\eta}+{r}^{\gamma}\rho vu{J}_{\xi}-{r}^{\gamma}\rho vv{J}_{\eta}=-J{r}^{\gamma}{p}_{\eta}+Jr{A}_{2},$
$\left(J{r}^{\gamma}\rho u\right(p+E){)}_{\xi}+(J{r}^{\gamma}\rho v(p+E){)}_{\eta}=\mathrm{\Phi},$
where $J=\sqrt{D}$; $u$, $v$ – Projections of the velocity vector $\mathbf{v}$ on the axis $\xi $, $\eta $; $E=\rho \left(e+\mathrm{0,5}\cdot {V}^{2}\right)$ – the total energy; $T$ – temperature. Axisymmetric flow corresponds to $\gamma =\text{1}$, a flat – $\gamma =\text{0}$. The coefficients in the right sides of the equations in the coordinate system $\xi $, $\eta $ is written as follows:
Odds ${Q}_{1}$, ${Q}_{2}$, ${Q}_{3}$ describe the exchange of momentum and energy between phases:
${Q}_{2}=n{B}_{u}\left(u-{u}_{s}\right),\mathrm{}\mathrm{}\mathrm{}\mathrm{}{Q}_{3}=n{B}_{u}\left(v-{v}_{s}\right),$
where, the index $s$ belongs to the dispersed phase; $G$ – The mass evaporation rate drops; ${c}_{s}$, ${\rho}_{l}$ – heat capacity and density of droplets of substance; $n$ – Condensed number of particles per unit volume; ${B}_{u}$, ${B}_{T}$ – Heat transfer coefficients and the resistance between the phases.
Together with the Eq. (3) it is necessary to solve the equation for the turbulent viscosity. In the plane $x$, $r$ orthogonal difference grid built complex boundary element method.
The basis of the solution of a system of Eq. (3) laid SIMPLE method [6], which is kind of focused on the use of curvilinear coordinates, implemented in [5, 7].
The equations for the dispersed phase is also recorded in a curvilinear coordinate system:
$(J{r}^{\gamma}{\rho}_{s}{u}_{s}{u}_{s}{)}_{\xi}+(J{r}^{\gamma}{\rho}_{s}{v}_{s}{u}_{s}{)}_{\eta}+{r}^{\gamma}{\rho}_{s}{v}_{s}{u}_{s}{J}_{\eta}-{r}^{\gamma}{\rho}_{s}{v}_{s}{v}_{s}{J}_{\xi}={J}^{2}{r}^{\gamma}{Q}_{2},$
$(J{r}^{\gamma}{\rho}_{s}{u}_{s}{v}_{s}{)}_{\xi}+(J{r}^{\gamma}{\rho}_{s}{v}_{s}{v}_{s}{)}_{\eta}+{r}^{\gamma}{\rho}_{s}{v}_{s}{u}_{s}{J}_{\xi}-{r}^{\gamma}{\rho}_{s}{v}_{s}{v}_{s}{J}_{\eta}={J}^{2}{r}^{\gamma}{Q}_{3},$
$(J{r}^{\gamma}{\rho}_{s}{u}_{s}{E}_{s}{)}_{\xi}+(J{r}^{\gamma}{\rho}_{s}{v}_{s}{E}_{s}{)}_{\eta}={J}^{2}{r}^{\gamma}{Q}_{1},\mathrm{}\mathrm{}(J{r}^{\gamma}{n}_{s}{u}_{s}{)}_{\xi}+(J{r}^{\gamma}{n}_{s}{v}_{s}{)}_{\eta}=0.$
Resistance and heat transfer coefficients are determined by the formulas [6]:
${B}_{T}=\frac{{B}_{u}Nu{c}_{p}}{3\mathrm{P}\mathrm{r}{c}_{s}f},Nu=2+\mathrm{0,459}\mathrm{R}{\mathrm{e}}_{}^{\mathrm{0,55}}\mathrm{P}{\mathrm{r}}^{\mathrm{0,33}}.$
3. The boundary conditions
On the permeable boundaries of the fixed costs of the gas and condensed phase, the enthalpy of the phase and the initial particle size. On impermeable boundaries puts a condition of sticking. At the outlet of the pipe for subsonic and where the pressure was set. On the axis of symmetry – the symmetry condition.
Changing the number n of particles per unit volume is determined by the speed lag of gas (including the interaction of the dispersed phase and gas). We write the equations of momentum and energy transfer:
where ${\varsigma}_{s}$ – coordinate in the direction from the $s$th path; ${W}_{s}=\sqrt{({U}_{s}{)}^{2}+({V}_{s}{)}^{2}}$; ${U}_{s}={u}_{s}{\xi}_{x}+{v}_{s}{\xi}_{r}$; ${V}_{s}={u}_{s}{\eta}_{x}+{v}_{s}{\eta}_{r}$; ${\mathbf{f}}_{s}=({u}_{s},{v}_{s},{T}_{s})$; ${\mathbf{C}}_{s}=\left[\begin{array}{l}{A}_{1s}\\ {A}_{1s}\\ {A}_{1T}\end{array}\right]$; ${\mathbf{B}}_{s}=\left[\begin{array}{c}{A}_{1s}u+{a}_{x}g\\ {A}_{1s}v+{a}_{y}g\\ {A}_{1Ts}T\end{array}\right]\mathrm{}\mathrm{}.$
Driving integration of the system Eq. (5) is shown in Fig. 1.
The integration is performed over the interval of the particle trajectory. The segment begins with the diagonal cell comprising an anchor point $\left(i,j\right)$, at 0 particles with known parameters, and ends at the node $(i,j)$, where the parameters to be calculated. To integrate the use of implicit Runge-Kutt second order of approximation:
An index of 0 corresponds to the beginning of the segment of the trajectory. Required parameters ${\phi}_{s0}$ at the point 0 determined by the formula:
using known values ${\phi}_{s}$ at the ends of the diagonal. The values ${\alpha}_{s}$, $\mathrm{\Delta}{\varsigma}_{s}$, $iu$, $jv$ are given by:
$\mathrm{\Delta}{\varsigma}_{s}=\frac{2}{{\alpha}_{s,i+iu,j}+{\alpha}_{s,i,j+jv}},\mathrm{}\mathrm{}\mathrm{}\mathrm{\Delta}{\xi}_{i}={\xi}_{i}-{\xi}_{i+iu},\mathrm{}\mathrm{}\mathrm{}\mathrm{\Delta}{\eta}_{j}={\eta}_{j}-{\eta}_{j+jv},$
$iu=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(-{U}_{s,i,j}),\mathrm{}\mathrm{}\mathrm{}\mathrm{}jv=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(-{V}_{s,i,j}).$
The equation for ${n}_{s}(\xi ,\eta )$ Eq. (4) is solved by countercurrent scheme [6]. The difference analog of the Eq. (4) has the form:
where ${a}_{{W}_{i,j}}=\mathrm{m}\mathrm{a}\mathrm{x}[\delta {\eta}_{j}{y}_{i-1/2,j}{D}_{i-1/2,j}{U}_{i-1/2,j},0]$, ${a}_{{E}_{i,j}}=\mathrm{m}\mathrm{a}\mathrm{x}[-\delta {\eta}_{j}{y}_{i+1/2,j}{D}_{i+1/2,j}{U}_{i+1/2,j},0]$, ${a}_{{S}_{i,j}}=\mathrm{m}\mathrm{a}\mathrm{x}[\delta {\xi}_{i}{y}_{i,j-1/2}{D}_{i,j-1/2}{V}_{i,j-1/2},0]$, ${a}_{{N}_{i,j}}=\mathrm{m}\mathrm{a}\mathrm{x}[-\delta {\xi}_{i}{y}_{i,j+1/2}{D}_{i,j+1/2}{V}_{i,j+1/2},0]$,
${a}_{{p}_{i,j}}={a}_{{W}_{i,j}}+{a}_{{E}_{i,j}}+{a}_{{N}_{i,j}}+{a}_{{S}_{i,j}}$.
The boundary conditions for the system Eq. (6).
Since the Eq. (6) recorded along the paths, which are characteristics that only the boundary conditions are placed on parts of the border, where the particles start their movement. Here you can set the initial parameters of the particles: ${u}_{0s}$, ${v}_{0s}$ – the initial velocity of the particles; ${T}_{0s}$ – The initial temperature of the particles.
For the Eq. (4) applies, counterflow circuit, so the boundary conditions are set similarly.
Thus, to solve the equations that describe the motion of particles, obtained difference scheme that monitors the direction of flow. To solve differential equations has the following two-level iterative process.
As an initial approximation for the particle velocities are taken gas velocity. Then, depending on the course of nature, there are two variants of the algorithm for solving this problem.
1st option. Within this field, there is no areas with a return (circulation) over. In this case, we realized marching along the longitudinal coordinate algorithm. Iterations on points held on the transverse coordinate. As the first value ${B}_{si,j}$, ${C}_{si,j}$ is taken ${B}_{s0}$, ${C}_{s0}$. Then iterate to refine the position of the point 0 on the diagonal of each cell and solved the Eq. (2-6). After that solved the Eq. (5) under the scheme Eq. (6) at each point $i$, $j$ and specified values ${B}_{si,j}$, ${C}_{si,j}$. In this embodiment, the parameters in the computer memory is stored on the particle layer on the layer $i$ and $i-\mathrm{}$1.
2nd option. The computational domain contains zones with reverse flow. In this case, point by point iteration are conducted throughout the computational domain. computer memory requirements in this embodiment is considerably higher than in the first embodiment. In order to achieve the residual from the equation of continuity ~10-4 takes only about 15 global iterations.
The presented numerical method allows the calculation of the velocity field of the gas and the dispersed phase, the trajectories of particle motion (by solving the equations $d\eta /d\xi ={V}_{s}/{U}_{s}$ when found ${U}_{s}$, ${V}_{s}$), the mass flow of particles deposited on the walls $G$.
Fig. 1. A scheme for integrating the equations of motion of particles
4. The current results of the calculations in complex configuration channels
For the numerical solution of the problem of motion of the gas and particles at the site of bending the pipe [4], a difference mesh, as shown in Fig. 2. Lines $\xi =const$ are longitudinal, a line $\eta =const$ orthogonal thereto.
The movement of the gas phase was calculated at a Reynolds number 50000. The vector velocity field is shown in Fig. 3.
As follows from the results of calculations, the flow field is not symmetric. First, during urged against the wall of a small radius, and a loan in the opposite direction. As a result, the flow deflection from the center of the eddy currents are formed, as shown in Fig. 3. The velocity profile is characteristic of a turbulent profile.
Particle motion is calculated for an equivalent particle diameter from 10 microns to 500 microns. The trajectories of particles of different sizes are shown in Figs. 4-7.
Particles of small size (10 microns) monitor the gaseous phase line current and the walls of virtually interacting. Larger particles 100 microns in diameter collide with the pipe wall after bending.
Large particles (300 and 500mkm) have a more direct path to the bending portion of the tube and collide with the wall with a high intensity. Low speed 300 micron particles can be trapped vortex gas flow and to be in it for a while. Large particles with a diameter of 500 mm, almost all come through the vortex. Only particles from the wall region with a very low speed to make a whirlwind return movement.
Fig. 2. Curved orthogonal net curved pipe
Fig. 3. The velocity field in the curved section of the pipe
Fig. 4. The particle trajectories diameter 10 mm
Fig. 5. The particle trajectories a diameter of 100 microns
Fig. 6. The particle trajectories diameter of 300 microns
Fig. 7. The particle trajectories diameter of 500 microns
As a result of solving the problem of the motion of a two-phase mixture of particles and gas through the curved portion of the pipe calculated mass flows falling on the wall of particles of different sizes. The results are shown in Fig. 8.
Fig. 8 shows the distribution of massive flows of precipitating particles $G$ in several sizes. Coordinate $S$ measured at the tube wall. Mass flow rate $G$ is related to the mass flow rate of particles on a straight section of pipe (${\rho}_{s}{W}_{s}$). These data can be used to assess mechanical impact on the particles pipeline wall at different degrees of drying of natural gas.
Narrowing the flow cross section can be various kinds of control and measurement devices. The computational domain and the difference curvilinear grid shown in Fig. 9.
Fig. 8. Distribution of mass flows deposited on particles of the wall of the pipe
Fig. 9. The difference in net contraction
Fig. 10. The flow field in the constriction
Fig. 11. The particle trajectories diameter 10 mm
Fig. 12. The particle trajectories diameter of 500 microns
Fig. 13. Mass flow deposited particulate
Calculation of the gas phase of motion parameters. The results are illustrated in Fig. 10.
Before narrowing of the pipe cross-section is observed stagnation zone with poor circulation flow in the wall region. After opening for jet nature is the presence of extensive areas (up to several pipe diameters) of the reverse flow. The movement of the condensed particles is shown in Fig. 11 and 12.
Small particles in this device follow the gas and do not fall into an obstacle.
Coarse particles have an intense effect on the structural member. precipitation intensity distribution for frequent diameter of 100-500 microns is shown in Fig. 13.
Here, the coordinate of $S$ corresponds to a distance from the pipe wall. Particles larger than 100 microns fall almost entire surface with obstacles fairly high intensity.
5. Conclusions
Figs. 6-8 demonstrate the possibility of increased erosion ash pipelines in areas of bending and rotation. The calculation of particle trajectories in the technological equipment of various kinds of shows at the option of either abrasive wear parts located in a stream, or an intense build-up and accumulation of hydrates.
The calculated results can be applied in the analysis of impact of solid particles of various sizes on the elements of gas pipelines.
References
- Strizhov I. N., Hodanovich I. E. Gas Production. Institute of Computer Science, Moscow-Izhevsk, 2003, p. 376, (in Russian). [Search CrossRef]
- Oran Je., Boris Dzh. Numerical Simulation of Reacting Flows. Second Edition Cambridge University Press, 2001, p. 640. [Search CrossRef]
- Sekundov A. N. The use of differential equations for the turbulent viscosity and analysis of non-self plane currents. Bulletin of the Academy of Sciences USSR, Vol. 5, 1971, p. 114-127, (in Russian). [Search CrossRef]
- Gromadka II T., Lej Ch The Complex Boundary Element Method in Applied Sciences. Mir, Moscow, 1990, p. 303, (in Russian). [Search CrossRef]
- Benderskij B. Ja, Tenenev V. A. Experimental and numerical study of flows in axisymmetric channels of complex shape with blowing. Bulletin of the Russian Academy of Sciences, Vol. 2, 2001, p. 24-28, (in Russian). [Search CrossRef]
- Patankar S. Numerical Methods for Solving Problems of Heat Transfer and Fluid Dynamics. Jenergoizdat, Moscow, 1984, p. 150, (in Russian). [Search CrossRef]
- Benderskij B. Ja, Tenenev V. A. Spatial subsonic flows in areas with complex geometry. Math Modeling, Vol. 13, Issue 8, 2001, p. 35-39, (in Russian). [Search CrossRef]
- Amimul Ahsan Two Phase Flow, Phase Change and Numerical Modeling. InTech, 2011, p. 596. [Search CrossRef]