Open Access
Numéro
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 75, 2020
Numéro d'article 71
Nombre de pages 23
DOI https://doi.org/10.2516/ogst/2020065
Publié en ligne 19 octobre 2020

© L. Mao et al., published by IFP Energies nouvelles, 2020

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Nomenclature

Qg: Gas production

Q g - H 2 S $ {Q}_{\mathrm{g}-{\mathrm{H}}_2\mathrm{S}}$ : Released H2S

ρg : Densities of gas

ρl : Density of drilling fluid

vg : Velocity of gas

vl : Velocity of drilling fluid

Eg : Gas void fraction

El : Liquid holdup

A : Cross-sectional area

g : Local acceleration of gravity

( P z ) fr $ {\left(\frac{\mathrm{\partial }P}{\mathrm{\partial }z}\right)}_{\mathrm{fr}}$ : Friction pressure drop

c1 : Specific heat capacity of the drilling fluid in the drilling string

T1 : Drilling fluid temperature in the drilling string

up : Velocity in the x direction of the drilling fluid in the drilling string

vp : Velocity in the y direction of the drilling fluid in the drilling string

Γx: Heat transfer coefficient of the drilling fluid in the x direction

Γy: Heat transfer coefficient of the drilling fluid in the y direction

ρ2 : Density of the drilling string

T2: Drilling string temperature

h : Well depth

Ql : Mud displacement

Vpg : Pit gain at shut-in time

P0 : Wellhead pressure

Pp : Formation pressure

vsg : Superficial velocities of gas

vsl : Superficial velocities of liquid

voo : Gas drift velocity

σ : Surface tension

k1, k2: Correction coefficients

yi : Mole fractions of H2S in the gas

xi : Mole fractions of H2S in the liquid

φ i V $ {\phi }_i^V$ : Fugacity coefficients of H2S in the gas

φ i l $ {\phi }_i^{\mathrm{l}}$ : Fugacity coefficients of H2S in the liquid

kij : Interaction coefficient between H2S and hydrocarbon components in natural gas

Tci : Critical temperature of H2S

Pci : Critical pressure of H2S

Tri : Correspondent temperature of H2S

w : Acentric factor of H2S

D : Borehole size

d : Drill pipe diameter

vgr : Drift velocity

C0 : Distribution coefficient

vm : Mixed velocity of gas and liquid

Hr: Viscous friction head

D e : Equivalent diameter

n : Flow index of the drilling fluid

R e : Reynolds number of the mixed fluids

εe : Equivalent absolute roughness

Pb : Bottom hole pressure

rw : Hole radius

K : Reservoir permeability

c : System compressibility

μg : Gas viscosity

t : Overflow time

Wg : Molar mass of natural gas

R : Molar gas constant

P a : Casing pressure

Pla : Circulating pressure loss

Pma : Hydrostatic pressure created by the fluid column

R e : Mean Reynolds number of the mixed fluids

φ i v $ {\phi }_i^v$ : Fugacity coefficients of H2S in the gas

φ i l $ {\phi }_i^{\mathrm{l}}$ : Fugacity coefficients of H2S in the liquid

Z : Natural gas compression factor

1 Introduction

Gas kick and well-killing are the main concerns of well control during the drilling of oil and gas fields. Various multiphase flow models have been proposed to predict the flow pattern distribution in the wellbore after gas kick (Amaya-Gómez et al., 2019; Ansari et al., 1994; Bilicki and Kestin, 1987; Ebrahimi and Khamehchi, 2015; Feng et al., 2015; Guo et al., 2018; Hasan and Kabir, 1988; Kelessidis and Dukler, 1989; Li and Mouline, 1997; Zhang et al., 2019); among them, the most widely used in drilling engineering are the patterns of bubble, slug, churn, and annular flows (Osman, 2004; Raghavan, 1989; Sun et al., 2017). On the basis of these studies, many researchers have focused on the multiphase flow in the wellbore under different conditions. Nickens (1987) proposed a dynamic computer model by a solution of the mass and momentum balance equations for predicting the dynamical two-phase flow behavior and pressure response of the wellbore. Sun and Wang proposed a mathematical model to predict the multiphase flow in a wellbore for deep water wells (Wang et al., 2016). Their model has also been used to study the phase transition of hydrate during the exploitation of natural gas hydrate (Wang and Sun, 2009, 2014). Xu et al. (2018, 2019) developed a non-isothermal two-phase flow model to investigate the effect of major parameters on the two-phase flow behavior in the wellbore. They found that temperature, pressure, and solubility fields are mutually influential. The gas solubility effect and heat transfer effect influence gas kick characteristics significantly. Meng et al. (2015) proposed a transient gas–liquid drift flux model based on Shi slip relation for simulating gas kick/injection and solved the model numerically using finite volume method and advection upstream splitting method V scheme. Xie et al. (2014) established a transient model to study the changing rule of wellhead parameters during blowout and found that the gas outlet velocity increases as well depth increases. Hou et al. (2019) conducted an experiment to study the heat transfer for gas-liquid two-phase flow in wellbore and developed a new model of convective heat transfer coefficient. Fu et al. (2019) proposed the mathematic model of wellbore annulus transient water hammer with the consideration of transient multi-phase flow characteristics and pointed out that both the additional water hammer pressure and shut-in casing pressure generated by the closure of BOP are likely to cause formation at the shallow casing shoe damage.

With the increasing demand for natural gas, an increasing number of gas reservoirs that contain H2S will be drilled (Guo and Wang, 2016; Yin et al., 2015; Zhang et al., 2020; Zhu et al., 2011). The existing states of H2S are greatly affected by temperature and pressure. The state may transform from supercritical to gas in the wellbore as annular pressure and temperature decrease, and this condition influences annular pressure and flow pattern distribution. Therefore, the typical gas–liquid two-phase flow models may be unsuitable for predicting the multiphase flow behavior in H2S-containing natural gas wells. The “12·23” Gas Blowout Accident in Kaixian County, Chongqing has elicited scholarly attention on phase transition of H2S in well control. Some scholars have focused on multiphase flow behavior of H2S-containing natural gas wells. Sun and Wang et al. established a multiphase flow model that considers the phase change of the H2S in the drilling fluid; simulation results showed that, as the invasion gas moves up along the wellbore to the critical position of H2S, the released H2S may cause a rapid volume expansion, thereby increasing the blowout risk (Sun et al., 2013, 2018).

It can be seen that most of the existing studies are focused on the gas kick simulation of a well and, the well-killing simulation after gas kick is scarce, especially for an H2S-containing natural gas well. By applying casing pressure in the wellhead during well killing, the two-phase flow behavior is different from that in the overflow (gas kick). Thus, the phase change behavior of H2S may be different from that described in the abovementioned published literature. Thus, it is of significance to study the well-killing process of an H2S-containing natural gas well to offer a guidance for exploitation of H2S-containing gas reservoir.

Thus, the aim of this study is to establish a dynamical well-killing model considering an H2S solubility to simulate the overflow and well-killing process of a vertical H2S-containing natural gas well. The mass and momentum equations of the coupled model were solved using finite difference method, while the transient temperature prediction model were solved using finite volume method. The continuity and momentum equations are decoupled from temperature prediction model. The effect of H2S content, mud displacement, drilling fluid density, and initial overflow volume on the dynamical well-killing process of an H2S-containing natural gas well were obtained and analyzed in this work.

2 Mathematical model

2.1 Governing equations for mass and momentum

During drilling, gas enters the wellbore from the stratum when overflow occurs, and gas–liquid two-phase flow appears below the location the gas reaches. Given that two-phase flow exists in overflow and well-killing periods, a mathematical model is established in this work to obtain the two-phase flow behavior.

To simplify the calculation, the following assumptions are made in this work:

  1. The gas and drilling fluid flow in a vertical wellbore is considered one dimensional;

  2. The compressibility of the drilling fluid is ignored;

  3. Gas and liquid phases are continuous in the control unit;

  4. The influence of annulus eccentricity is disregarded;

  5. The geothermal gradient is regarded as constant.

Then, the continuity equations of gas and liquid phases are expressed as follows (Wei et al., 2018):(ρgvgEgA)z+(ρgEgA)t={Qg+Qg-H2S(overflow period)Qg-H2S(well-killing period),$$ \frac{\mathrm{\partial }({\rho }_{\mathrm{g}}{v}_{\mathrm{g}}{E}_{\mathrm{g}}A)}{\mathrm{\partial }z}+\frac{\mathrm{\partial }({\rho }_{\mathrm{g}}{E}_{\mathrm{g}}A)}{\mathrm{\partial }t}=\left\{\begin{array}{c}{Q}_{\mathrm{g}}+{Q}_{\mathrm{g}-{\mathrm{H}}_2\mathrm{S}}(\mathrm{overflow}\enspace \mathrm{period})\\ {Q}_{\mathrm{g}-{\mathrm{H}}_2\mathrm{S}}(\mathrm{well}-\mathrm{killing}\enspace \mathrm{period})\end{array}\right., $$(1)(ρlvlElA)z+(ρlElA)t=0.$$ \frac{\mathrm{\partial }({\rho }_{\mathrm{l}}{v}_{\mathrm{l}}{E}_{\mathrm{l}}A)}{\mathrm{\partial }z}+\frac{\mathrm{\partial }({\rho }_{\mathrm{l}}{E}_{\mathrm{l}}A)}{\mathrm{\partial }t}=0. $$(2)

On the basis of the law of conservation of momentum, the momentum equation can be obtained as follows:t(ρlvlEl+ρgvgEg)+z(ρlvl2El+ρgvg2Eg)+Pz+(ρlEl+ρgEg)g+(Pz)fr=0.$$ \frac{\mathrm{\partial }}{\mathrm{\partial }t}({\rho }_{\mathrm{l}}{v}_{\mathrm{l}}{E}_{\mathrm{l}}+{\rho }_{\mathrm{g}}{v}_{\mathrm{g}}{E}_{\mathrm{g}})+\frac{\mathrm{\partial }}{\mathrm{\partial }z}({\rho }_{\mathrm{l}}{v}_{\mathrm{l}}^2{E}_{\mathrm{l}}+{\rho }_{\mathrm{g}}{v}_{\mathrm{g}}^2{E}_{\mathrm{g}})+\frac{\mathrm{\partial }P}{\mathrm{\partial }z}+({\rho }_{\mathrm{l}}{E}_{\mathrm{l}}+{\rho }_{\mathrm{g}}{E}_{\mathrm{g}})g+{\left(\frac{\mathrm{\partial }P}{\mathrm{\partial }z}\right)}_{\mathrm{fr}}=0. $$(3)

2.2 Transient temperature prediction model

To accurately predict the phase change behavior of H2S in wellbore, a transient temperature prediction model is developed (Mao and Zhang, 2018). The unsteady two-dimensional convection–diffusion and unsteady two-dimensional diffusion equations are used to describe the heat transfer models as follows:

(1) Heat transfer model in the drilling string:(ρ1c1T1)t+x(ρ1c1upT1)+y(ρ1c1vpT1)=x(Γ1xc1T1x)+y(Γ1yc1T1y)+Sp.$$ \frac{\mathrm{\partial }({\rho }_1{c}_1{T}_1)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left({\rho }_1{c}_1{u}_{\mathrm{p}}{T}_1\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }y}\left({\rho }_1{c}_1{v}_{\mathrm{p}}{T}_1\right)=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left({\mathrm{\Gamma }}_{1x}{c}_1\frac{\mathrm{\partial }{T}_1}{\mathrm{\partial }x}\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }y}\left({\mathrm{\Gamma }}_{1y}{c}_1\frac{\mathrm{\partial }{T}_1}{\mathrm{\partial }y}\right)+{S}_{\mathrm{p}}. $$(4)

(2) Heat transfer model of drilling string:(ρ2c2T2)t=x(Γ2xc2T2x)+y(Γ2yc2T2y).$$ \frac{\mathrm{\partial }({\rho }_2{c}_2{T}_2)}{\mathrm{\partial }t}=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left({\mathrm{\Gamma }}_{2x}{c}_2\frac{\mathrm{\partial }{T}_2}{\mathrm{\partial }x}\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }y}\left({\mathrm{\Gamma }}_{2\mathrm{y}}{c}_2\frac{\mathrm{\partial }{T}_2}{\mathrm{\partial }y}\right). $$(5)

(3) Heat transfer model in the annular:(ρ3c3T3)t+x(ρ3c3uaT3)+y(ρ3c3vaT3)=x(Γ3xc3T3x)+y(Γ3yc3T3y)+Sa.$$ \frac{\mathrm{\partial }({\rho }_3{c}_3{T}_3)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left({\rho }_3{c}_3{u}_{\mathrm{a}}{T}_3\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }y}\left({\rho }_3{c}_3{v}_{\mathrm{a}}{T}_3\right)=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left({\mathrm{\Gamma }}_{3\mathrm{x}}{c}_3\frac{\mathrm{\partial }{T}_3}{\mathrm{\partial }x}\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }y}\left({\mathrm{\Gamma }}_{3\mathrm{y}}{c}_3\frac{\mathrm{\partial }{T}_3}{\mathrm{\partial }y}\right)+{S}_{\mathrm{a}}. $$(6)

(4) Heat transfer model of casing, cement sheath, and formation:(ρiciTi)t=x(ΓixciTix)+y(ΓiyciTiy).$$ \frac{\mathrm{\partial }({\rho }_i{c}_i{T}_i)}{\mathrm{\partial }t}=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left({\mathrm{\Gamma }}_{{ix}}{c}_i\frac{\mathrm{\partial }{T}_i}{\mathrm{\partial }x}\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }y}\left({\mathrm{\Gamma }}_{{iy}}{c}_i\frac{\mathrm{\partial }{T}_i}{\mathrm{\partial }y}\right). $$(7)

2.2 Initial and boundary conditions

(1) Overflow period

When overflow occurs, the wellhead is open and the parameters under normal drilling conditions can be used as the initial conditions of a well kick:{P(h,0)=PbQ(h,0)=QlEg(h,0)=0El(h,0)=1vl(h,0)=Ql/A.$$ \left\{\begin{array}{c}P\left(h,0\right)={P}_{\mathrm{b}}\\ Q\left(h,0\right)={Q}_{\mathrm{l}}\\ \begin{array}{c}{E}_{\mathrm{g}}(h,0)=0\\ {E}_{\mathrm{l}}(h,0)=1\\ {v}_{\mathrm{l}}(h,0)={Q}_{\mathrm{l}}/A\end{array}\end{array}.\right. $$(8)

The boundary conditions are shown as follows:{P(0,t)=P0Q(h,t)=Ql+Qg(h,t)Pg=Vpg.$$ \left\{\begin{array}{c}P\left(0,t\right)={P}_0\\ Q\left(h,t\right)={Q}_{\mathrm{l}}+{Q}_{\mathrm{g}}(h,t)\\ \mathrm{Pg}={V}_{\mathrm{pg}}\end{array}\right.. $$(9)

(2) Killing period

After the shut-in operation is performed, the flow pattern, gas void fraction, and annular pressure distributions at the shut-in time can be considered the initial conditions of the well-killing period. During the well-killing process, the bottom hole pressure is equal to the formation pressure, and the boundary conditions are shown as follows:{P(0,t)=PpQl(0,t)=QlQg(0,t)=0.$$ \left\{\begin{array}{c}P(0,t)={P}_{\mathrm{p}}\\ {Q}_{\mathrm{l}}(0,t)={Q}_{\mathrm{l}}\\ {Q}_{\mathrm{g}}(0,t)=0\end{array}\right.. $$(10)

2.3 Flow pattern discriminant

The main patterns of a gas–liquid two-phase flow are bubble, slug, churn, and annular flows (Hasan and Kabir, 1988). They can be discriminated by the following equations:

Bubble flow:vsgk1×(0.429×vsl+0.357×voo),$$ {v}_{\mathrm{sg}}\le {k}_1\times \left(0.429\times {v}_{\mathrm{sl}}+0.357\times {v}_{\mathrm{oo}}\right), $$(11)voo=1.53[g(ρl-ρg)σρl2]0.25.$$ {v}_{\mathrm{oo}}=1.53{\left[\frac{g({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}})\sigma }{{{\rho }_{\mathrm{l}}}^2}\right]}^{0.25}. $$(12)

Slug flow:vsg>k1×(0.429×vsl+0.357×voo).$$ {v}_{\mathrm{sg}}>{k}_1\times \left(0.429\times {v}_{\mathrm{sl}}+0.357\times {v}_{\mathrm{oo}}\right). $$(13)

Churn flow:{vsg<3.1[g(ρl-ρg)σρg2]0.25ρgvg2>25.41g(ρlvsl2)-39ifρlvsl2>74.4ρgvg2>0.0051(ρlvsl2)1.7ifρlvsl274.4.$$ \left\{\begin{array}{c}{v}_{\mathrm{sg}} < 3.1{\left[\frac{g({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}})\sigma }{{\rho }_{\mathrm{g}}^2}\right]}^{0.25}\\ {\rho }_{\mathrm{g}}{v}_{\mathrm{g}}^2>25.41g\left({\rho }_{\mathrm{l}}{v}_{\mathrm{sl}}^2\right)-39\to {if}{\rho }_{\mathrm{l}}{v}_{\mathrm{sl}}^2>74.4\\ {\rho }_{\mathrm{g}}{v}_{\mathrm{g}}^2>0.0051({\rho }_{\mathrm{l}}{v}_{\mathrm{sl}}^2{)}^{1.7}\to {if}{\rho }_{\mathrm{l}}{v}_{\mathrm{sl}}^2\le 74.4\end{array}\right.. $$(14)

Annular flow:vsg>3.1[g(ρl-ρg)σρg2]0.25.$$ {v}_{\mathrm{sg}}>3.1{\left[\frac{g({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}})\sigma }{{\rho }_{\mathrm{g}}^2}\right]}^{0.25}. $$(15)

2.4 Calculation of key parameters

2.4.1 H2S solubility

The solubility of H2S in the drilling fluid can be expressed as (Sun et al., 2018):xi=pyiφivpφil,$$ {x}_i=\frac{p{y}_i{\phi }_i^v}{p{\phi }_i^{\mathrm{l}}}, $$(16)P=RTV-b-aV(V+b)+b(V-b).$$ P=\frac{{RT}}{V-b}-\frac{a}{V(V+b)+b(V-b)}. $$(17)

With A=apR2T2$ A=\frac{{ap}}{{R}^2{T}^2}$, B=bpRT$ B=\frac{{bp}}{{RT}}$, and Z=PVRT$ Z=\frac{{PV}}{{RT}}$, equation (13) can be written asZ3-(1-B)Z2+(A-3B2-2B)Z-(AB-B2-B3)=0,$$ {Z}^3-(1-B){Z}^2+(A-3{B}^2-2B)Z-({AB}-{B}^2-{B}^3)=0, $$(18)a=ijxixjαiαj(1-kij),$$ a=\sum_i \sum_j {x}_i{x}_j\sqrt{{\alpha }_i{\alpha }_j}\left(1-{k}_{{ij}}\right), $$(19)b=ixi0.0778RTciPci,$$ b=\sum_i {x}_i\frac{0.0778R{T}_{\mathrm{ci}}}{{P}_{\mathrm{ci}}}, $$(20)αi=[1+m(1-Tri0.5)]2,$$ {\alpha }_i={\left[1+m(1-{T}_{\mathrm{ri}}^{0.5})\right]}^2, $$(21)αj=0.45724R2Tci2Pciαi,$$ {\alpha }_j=\frac{0.45724{R}^2{T}_{\mathrm{ci}}^2}{{P}_{\mathrm{ci}}}{\alpha }_i, $$(22)andm=0.37464+1.5226-0.26992w2.$$ m=0.37464+1.5226-0.26992{w}^2. $$(23)

Thus, the fugacity coefficient of a certain component is expressed as:lnφi=bib(Z-1)-ln(Z-b)-a22b(2jxiαija-bib)ln[Z+(1+2)bZ+122b].$$ \mathrm{ln}{\phi }_i=\frac{{b}_i}{b}\left(Z-1\right)-\mathrm{ln}\left(Z-b\right)-\frac{a}{2\sqrt{2}b}\left(\frac{2\sum_j {x}_i{{\alpha }_i}_j}{a}-\frac{{b}_i}{b}\right)\mathrm{ln}\left[\frac{Z+\left(1+\sqrt{2}\right)b}{Z+12\sqrt{2}b}\right]. $$(24)

The solubility of H2S in the liquid can be obtained by combining equation (16) with equation (24).

2.4.2 Drift flux model

In the drift flux model, the corresponding slip rate of the gas phases and distribution coefficient in case of different flow patterns are shown as follows (Wei et al., 2018):

Bubble flow:vgr=1.53(g(ρl-ρg)σρl2)0.25,$$ {v}_{\mathrm{gr}}=1.53{\left(\frac{g({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}})\sigma }{{{\rho }_{\mathrm{l}}}^2}\right)}^{0.25}, $$(25)C0=1.20+0.371×dD.$$ {C}_0=1.20+0.371\times \frac{d}{D}. $$(26)

Slug flow:vgr=0.35(gD(ρl-ρg)ρl)0.5,$$ {v}_{\mathrm{gr}}=0.35{\left(\frac{{gD}({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}})}{{\rho }_{\mathrm{l}}}\right)}^{0.5}, $$(27)C0=1.182+0.9×dD.$$ {C}_0=1.182+0.9\times \frac{d}{D}. $$(28)

Churn and annular flows:vgr=(0.35+0.22Dd)(g(D-d)(ρl-ρg)ρl)0.5,C0=1.$$ \begin{array}{c}{v}_{\mathrm{gr}}=\left(0.35+0.22\frac{D}{d}\right){\left(\frac{g\left(D-d\right)({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}})}{{\rho }_{\mathrm{l}}}\right)}^{0.5},\\ {C}_0=1.\end{array} $$(29)

Thus, the gas void fraction can be calculated by the following equations:Eg=vsgCovm+vgr,$$ {E}_{\mathrm{g}}=\frac{{v}_{\mathrm{sg}}}{{C}_o{v}_{\mathrm{m}}+{v}_{\mathrm{gr}}}, $$(30)vm=vsg+vsl,$$ {v}_{\mathrm{m}}={v}_{\mathrm{sg}}+{v}_{\mathrm{sl}}, $$(31)

2.4.3 Friction pressure drop

For single phase flow, the friction pressure drop can be obtained using the friction factor of the power-law fluid.Hr=2fvl2ρl/De,$$ H\mathrm{r}=2f{{v}_{\mathrm{l}}}^2{\rho }_{\mathrm{l}}/{D}_{\mathrm{e}}, $$(32){f=8vl2ρl[8vlDe3n+14n]n,if Re 20001f=2.1n0.75log[Re(f4)1-n/2]-0.2n1.2>2000.$$ \left\{\begin{array}{c}f=\frac{8}{{v}_{\mathrm{l}}^2{\rho }_{\mathrm{l}}}{\left[\frac{8{v}_{\mathrm{l}}}{{D}_{\mathrm{e}}}\frac{3n+1}{4n}\right]}^n,\hspace{1em}\mathrm{if\enspace Re\enspace }2000\\ \frac{1}{\sqrt{f}}=\frac{2.1}{{n}^{0.75}}\mathrm{log}\left[\mathrm{Re}{\left(\frac{f}{4}\right)}^{1-n/2}\right]-\frac{0.2}{{n}^{1.2}}>2000.\end{array}\right. $$(33)

For gas–liquid two-phase flow, the friction pressure drop can be calculated using the following equations established by Yin et al. (2017).

Bubble flow:Hr=2fvm2ρm/De,$$ H\mathrm{r}=2f{{v}_{\mathrm{m}}}^2{\rho }_{\mathrm{m}}/{D}_e, $$(34){1f=-4log(εe3.71De-5.05logARe),A=m×(εe2.56De)1.11+(7.149Re)0.898.$$ \left\{\begin{array}{c}\frac{1}{\sqrt{f}}=-4\mathrm{log}\left(\frac{{\epsilon }_{\mathrm{e}}}{3.71{D}_{\mathrm{e}}}-\frac{5.05\mathrm{log}A}{\mathrm{Re}}\right),\\ A=m\times {\left(\frac{{\epsilon }_e}{2.56{D}_{\mathrm{e}}}\right)}^{1.11}+{\left(\frac{7.149}{\mathrm{Re}}\right)}^{0.898}.\end{array}\right. $$(35)

Slug flow:Hr=2(1-Eg)fvm2ρm/De$$ H\mathrm{r}=2(1-{E}_{\mathrm{g}})f{{v}_{\mathrm{m}}}^2{\rho }_{\mathrm{m}}/{D}_e $$(36){1f=-4log(εe3.71De-5.05logARe)A=(εe2.56De)1.11+(7.149Re)0.898$$ \left\{\begin{array}{c}\frac{1}{\sqrt{f}}=-4\mathrm{log}\left(\frac{{\epsilon }_e}{3.71{D}_e}-\frac{5.05\mathrm{log}A}{{Re}}\right)\\ A={\left(\frac{{\epsilon }_e}{2.56{D}_e}\right)}^{1.11}+{\left(\frac{7.149}{{Re}}\right)}^{0.898}\end{array}\right. $$(37)

Churn and annular flows:Hr=2fvm2ρm/(DeEg2)$$ {Hr}=2f{v}_{\mathrm{m}}^2{\rho }_{\mathrm{m}}/({D}_{\mathrm{e}}{E}_{\mathrm{g}}^2) $$(38)f=0.079[1+75(1-Eg)]/(Reg0.25)$$ f=0.079\left[1+75(1-{E}_{\mathrm{g}})\right]/(\mathrm{Re}{g}^{0.25}) $$(39)

2.4.4 Gas production

Under the effect of pressure differential between bottom hole and formation pressures, the H2S-containing natural gas enters the wellbore. In this work, the following equation is adopted to calculate gas production (Elsharkawy, 2004):Qg=2.64×10-20Kh(Pp2-Pb2)(0.8+IntD)[(T-255)Zμg],$$ {Q}_{\mathrm{g}}=\frac{2.64\times 1{0}^{-20}\mathrm{Kh}({P}_{\mathrm{p}}^2-{P}_{\mathrm{b}}^2)}{(0.8+{\mathrm{Int}}_{\mathrm{D}})\left[(T-255)Z{\mu }_{\mathrm{g}}\right]}, $$(40)tD=max{10,1.47×10-9trw2(Kμg)}.$$ {t}_D=\mathrm{max}\left\{10,\frac{1.47\times 1{0}^{-9}t}{{{r}_w}^2}\left(\frac{K}{{c\phi }{\mu }_{\mathrm{g}}}\right)\right\}. $$(41)

2.4.5 Gas density

The gas density can be described by introducing a compression factor into the ideal gas equation, as defined in the equation below (Yin et al., 2017):ρg=PWgZRT.$$ {\rho }_{\mathrm{g}}=\frac{P{W}_{\mathrm{g}}}{{ZRT}}. $$(42)

2.4.6 Convective heat transfer coefficients

Convective heat transfer coefficients can be expressed as follows (Mao and Zhang, 2018):h=Nuλd.$$ h=\frac{\mathrm{Nu}\cdot \lambda }{d}. $$(43)

For laminar flow, Nu = 4.36; for turbulent flow, Nu = 0.023 Re0.8 Pr0.33.

3 Model solution

3.1 Solution for mass and momentum governing equations

The finite difference method is used in this work to solve the mass and momentum governing equations. Meanwhile, the influence of the boundary and initial conditions in each control region are considered in the solving model. The solution is completed by three steps: generating discrete grids, constructing discrete equations, and solving these equations.

(1) Generating discrete grids

In finite difference numerical calculation, the spatial domain is the wellbore. For the overflow simulation, the time domain is the entire time from the overflow to the shut-in. For the well-killing simulation, the time domain is from the shut-in to the time the lower boundary of two-phase reaches wellhead. The schematic of discrete grid nodes of the wellbore is shown in Figure 1a.

thumbnail Fig. 1

Schematic diagram of solution for the mass and momentum equations: (a) Schematic of discrete grid nodes of the wellbore; (b) Diagram of cell grid integration area K.

(2) Model discretization

Figure 1b shows the cell grid integration area K. The partial differential equation in the mathematical model can be written asXt+Yz=0.$$ \frac{\mathrm{\partial }X}{\mathrm{\partial }t}+\frac{\mathrm{\partial }Y}{\mathrm{\partial }z}=0. $$(44)

Then, by integrating this equation into area K, a curvilinear integral along the boundary of L can be obtained in accordance with Green’s theorem:(Xt+Yz)dtdz=LXdz-Ydt=L1+L2+L3+L4=-t(i)t(i+1)Y(t,z(j))dt+z(j)z(j+1)X(t(i+1),z)dz+t(i)t(i+1)Y(t,z(j+1))dt-z(j)z(j+1)X(t(i),z)dz=0.$$ \iint \left(\frac{\mathrm{\partial }X}{\mathrm{\partial }t}+\frac{\mathrm{\partial }Y}{\mathrm{\partial }z}\right)\mathrm{d}t\mathrm{d}z=\underset{L}{\int } X\mathrm{d}z-Y\mathrm{d}t=\underset{L1}{\int } +\underset{L2}{\int } +\underset{L3}{\int } +\underset{L4}{\int } =-\underset{t(i)}{\overset{t\left(i+1\right)}{\int }} Y(t,z(j))\mathrm{d}t+\underset{z(j)}{\overset{z\left(j+1\right)}{\int }} X(t(i+1),z)\mathrm{d}z+\underset{t(i)}{\overset{t\left(i+1\right)}{\int }} Y(t,z(j+1))\mathrm{d}t-\underset{z(j)}{\overset{z\left(j+1\right)}{\int }} X(t(i),z)\mathrm{d}z=0. $$(45)

The above-mentioned equation can be converted to the following equation by simplification:Yi+1j+1-Yi+1j=Δz2Δt(Xji+Xj+1i-Xji+1-Xj+1i+1),$$ {Y}_{i+1}^{j+1}-{Y}_{i+1}^j=\frac{\Delta z}{2\Delta t}\left({X}_j^i+{X}_{j+1}^i-{X}_j^{i+1}-{X}_{j+1}^{i+1}\right), $$(46)

(3) Numeralization of the continuity equation

For gas-phase continuity equation, we let{X=ρgEg,Y=ρgEgvg.$$ \left\{\begin{array}{c}X={\rho }_{\mathrm{g}}{E}_{\mathrm{g}},\\ Y={\rho }_{\mathrm{g}}{E}_{\mathrm{g}}{v}_{\mathrm{g}}.\end{array}\right. $$(47)

Thus, the gas-phase difference equation can be obtained by combining equations (46) and (47):[(ρgVgEg)i+1j+1-(ρgVgEg)ij+1]=0.5Δz/Δt[(ρgEg)ij+(ρgEg)ij+1-(ρgEg)ji+1-(ρgEg)i+1j+1]+0.5Δz(Qg-H2Si+1j+1+Qg-H2Si+1j)$$ \left[({\rho }_{\mathrm{g}}{V}_{\mathrm{g}}{E}_{\mathrm{g}}{)}_{i+1}^{j+1}-({\rho }_{\mathrm{g}}{V}_{\mathrm{g}}{E}_{\mathrm{g}}{)}_i^{j+1}\right]=0.5\Delta z/\Delta t\left[({\rho }_{\mathrm{g}}{E}_{\mathrm{g}}{)}_i^j+({\rho }_{\mathrm{g}}{E}_{\mathrm{g}}{)}_i^{j+1}-({\rho }_{\mathrm{g}}{E}_{\mathrm{g}}{)}_j^{i+1}-({\rho }_{\mathrm{g}}{E}_{\mathrm{g}}{)}_{i+1}^{j+1}\right]+0.5\Delta z\left({{Q}_{\mathrm{g}-{\mathrm{H}}_2\mathrm{S}}}_{i+1}^{j+1}+{{Q}_{\mathrm{g}-{\mathrm{H}}_2\mathrm{S}}}_{i+1}^j\right) $$(48)

Similarly, the liquid-phase difference equation can be obtained as follows:[(ρlVlEl)i+1j+1-(ρlVlEl)ij+1]=0.5Δz/Δt[(ρlEl)ij+(ρlEl)ij+1-(ρlEl)ji+1-(ρlEl)i+1j+1].$$ \left[({\rho }_{\mathrm{l}}{V}_{\mathrm{l}}{E}_{\mathrm{l}}{)}_{i+1}^{j+1}-({\rho }_{\mathrm{l}}{V}_{\mathrm{l}}{E}_{\mathrm{l}}{)}_i^{j+1}\right]=0.5\Delta z/\Delta t\left[({\rho }_{\mathrm{l}}{E}_{\mathrm{l}}{)}_i^j+({\rho }_{\mathrm{l}}{E}_{\mathrm{l}}{)}_i^{j+1}-({\rho }_{\mathrm{l}}{E}_{\mathrm{l}}{)}_j^{i+1}-{\left({\rho }_{\mathrm{l}}{E}_{\mathrm{l}}\right)}_{i+1}^{j+1}\right]. $$(49)

(4) Numeralization of the momentum equation

For mixed momentum equation, we let{X=ρgEgvg+ρlElvl,Y=ρgEgvg2+ρlElvl2+ρgEg+ρlElg+P+Pfr.$$ \left\{\begin{array}{c}X={\rho }_{\mathrm{g}}{E}_{\mathrm{g}}{v}_{\mathrm{g}}+{\rho }_{\mathrm{l}}{E}_{\mathrm{l}}{v}_{\mathrm{l}},\\ Y={\rho }_{\mathrm{g}}{E}_{\mathrm{g}}{{v}_{\mathrm{g}}}^2+{\rho }_{\mathrm{l}}{E}_{\mathrm{l}}{{v}_{\mathrm{l}}}^2+{\rho }_{\mathrm{g}}{E}_{\mathrm{g}}+{\rho }_{\mathrm{l}}{E}_{\mathrm{l}}\mathrm{g}+P+{P}_{\mathrm{fr}}.\end{array}\right. $$(50)

Thus, the mixed momentum difference equation can be obtained by combining equations (46) and (50):Pi+1j+1-Pi+1j=Δz2Δt{[(ρlvlEl+ρgvgEg)]ij+[(ρlvlEl+ρgvgEg)]i+1j+1-[(ρlvlEl+ρgvgEg)]i+1j-[(ρlvlEl+ρgvgEg)]i+1j+1}+(ρlvl2El+ρgvg2Eg)i+1j-0.5Δzg(ρlEl+ρgEg)i+1j+-(ρlvl2El+ρgvg2Eg)i+1j+1-0.5Δzg(ρlEl+ρgEg)i+1j+1-0.5Δz[(Pz)fr]i+1j-0.5Δz[(Pz)fr]i+1j+1.$$ {P}_{i+1}^{j+1}-{P}_{i+1}^j=\frac{\Delta z}{2\Delta t}\left\{[({\rho }_{\mathrm{l}}{v}_{\mathrm{l}}{E}_{\mathrm{l}}+{\rho }_{\mathrm{g}}{v}_{\mathrm{g}}{E}_{\mathrm{g}}){]}_i^j+[({\rho }_{\mathrm{l}}{v}_{\mathrm{l}}{E}_{\mathrm{l}}+{\rho }_{\mathrm{g}}{v}_{\mathrm{g}}{E}_{\mathrm{g}}){]}_{i+1}^{j+1}-[({\rho }_{\mathrm{l}}{v}_{\mathrm{l}}{E}_{\mathrm{l}}+{\rho }_{\mathrm{g}}{v}_{\mathrm{g}}{E}_{\mathrm{g}}){]}_{i+1}^j-[({\rho }_{\mathrm{l}}{v}_{\mathrm{l}}{E}_{\mathrm{l}}+{\rho }_{\mathrm{g}}{v}_{\mathrm{g}}{E}_{\mathrm{g}}){]}_{i+1}^{j+1}\right\}+({\rho }_{\mathrm{l}}{v}_{\mathrm{l}}^2{E}_{\mathrm{l}}+{\rho }_{\mathrm{g}}{v}_{\mathrm{g}}^2{E}_{\mathrm{g}}{)}_{i+1}^j-0.5\Delta {zg}({\rho }_{\mathrm{l}}{E}_{\mathrm{l}}+{\rho }_{\mathrm{g}}{E}_{\mathrm{g}}{)}_{i+1}^j+-({\rho }_{\mathrm{l}}{v}_{\mathrm{l}}^2{E}_{\mathrm{l}}+{\rho }_{\mathrm{g}}{v}_{\mathrm{g}}^2{E}_{\mathrm{g}}{)}_{i+1}^{j+1}-0.5\Delta {zg}({\rho }_{\mathrm{l}}{E}_{\mathrm{l}}+{\rho }_{\mathrm{g}}{E}_{\mathrm{g}}{)}_{i+1}^{j+1}-0.5\Delta z{\left[{\left(\frac{\mathrm{\partial }P}{\mathrm{\partial }z}\right)}_{{fr}}\right]}_{i+1}^j-0.5\Delta z{\left[{\left(\frac{\mathrm{\partial }P}{\mathrm{\partial }z}\right)}_{{fr}}\right]}_{i+1}^{j+1}. $$(51)

3.2 Solution for transient temperature prediction model

Figure 2 shows the grid diagram of the wellbore and formation. By using finite volume method, equations (4)(7) can be expressed as (Mao and Zhang, 2018):j=1aPTP=aNNTNN+aNTN+aWTW+aSTS+aETE+b,$$ j=1\hspace{1em}{a}_{\mathrm{P}}{T}_{\mathrm{P}}={a}_{\mathrm{NN}}{T}_{\mathrm{NN}}+{a}_{\mathrm{N}}{T}_{\mathrm{N}}+{a}_{\mathrm{W}}{T}_{\mathrm{W}}+{a}_{\mathrm{S}}{T}_{\mathrm{S}}+{a}_{\mathrm{E}}{T}_{\mathrm{E}}+b, $$(52)j=2aPTP=aNTN+aWTW+aSTS+aETE+b,$$ j=2\hspace{1em}{a}_{\mathrm{P}}{T}_{\mathrm{P}}={a}_{\mathrm{N}}{T}_{\mathrm{N}}+{a}_{\mathrm{W}}{T}_{\mathrm{W}}+{a}_{\mathrm{S}}{T}_{\mathrm{S}}+{a}_{\mathrm{E}}{T}_{\mathrm{E}}+b, $$(53)j=3aPTP=aSSTSS+aNTN+aWTW+aSTS+aETE+b,$$ j=3\hspace{1em}{a}_{\mathrm{P}}{T}_{\mathrm{P}}={a}_{\mathrm{SS}}{T}_{\mathrm{SS}}+{a}_{\mathrm{N}}{T}_{\mathrm{N}}+{a}_{\mathrm{W}}{T}_{\mathrm{W}}+{a}_{\mathrm{S}}{T}_{\mathrm{S}}+{a}_{\mathrm{E}}{T}_{\mathrm{E}}+b, $$(54)j4aPTP=aNTN+aWTW+aSTS+aETE+b.$$ j\le 4\hspace{1em}{a}_{\mathrm{P}}{T}_{\mathrm{P}}={a}_{\mathrm{N}}{T}_{\mathrm{N}}+{a}_{\mathrm{W}}{T}_{\mathrm{W}}+{a}_{\mathrm{S}}{T}_{\mathrm{S}}+{a}_{\mathrm{E}}{T}_{\mathrm{E}}+b. $$(55)

thumbnail Fig. 2

The grid diagram of the wellbore and formation.

By using under-relaxation iteration method, discretized scheme of the heat transfer control equations are obtained as follows:j=1Ti,1t+Δt=Ti,1t+ω[(aNNi1t+ΔtTi-2,1t+Δt+aNi1t+ΔtTi-1,1t+Δt+aSi1t+ΔtTi+1,1t+Δt+aEi1t+ΔtTi,2t+Δt+bi1t+Δt)-aPi1t+ΔtTi,1t]aPi1t+Δt,$$ j=1\hspace{1em}{T}_{i,1}^{t+\Delta t}={T}_{i,1}^t+\frac{\omega \left[\left({{a}_{\mathrm{NN}}}_{i1}^{t+\Delta t}{T}_{i-\mathrm{2,1}}^{t+\Delta t}+{{a}_{\mathrm{N}}}_{i1}^{t+\Delta t}{T}_{i-\mathrm{1,1}}^{t+\Delta t}+{{a}_{\mathrm{S}}}_{i1}^{t+\Delta t}{T}_{i+\mathrm{1,1}}^{t+\Delta t}+{{a}_{\mathrm{E}}}_{i1}^{t+\Delta t}{T}_{i,2}^{t+\Delta t}+{b}_{i1}^{t+\Delta t}\right)-{{a}_{\mathrm{P}}}_{i1}^{t+\Delta t}{T}_{i,1}^t\right]}{{{a}_{\mathrm{P}}}_{i1}^{t+\Delta t}}, $$(56)j=2Ti,2t+Δt=Ti,2t+ω[(aNi2t+ΔtTi-1,2t+Δt+aWi2t+ΔtTi,1t+Δt+aSi2t+ΔtTi+1,2t+Δt+aEi2t+ΔtTi,3t+Δt+bi2t+Δt)-aPi2t+ΔtTi,2t]aPi2t+Δt,$$ j=2\hspace{1em}{T}_{i,2}^{t+\Delta t}={T}_{i,2}^t+\frac{\omega \left[\left({{a}_{\mathrm{N}}}_{i2}^{t+\Delta t}{T}_{i-\mathrm{1,2}}^{t+\Delta t}+{{a}_{\mathrm{W}}}_{i2}^{t+\Delta t}{T}_{i,1}^{t+\Delta t}+{{a}_{\mathrm{S}}}_{i2}^{t+\Delta t}{T}_{i+\mathrm{1,2}}^{t+\Delta t}+{{a}_{\mathrm{E}}}_{i2}^{t+\Delta t}{T}_{i,3}^{t+\Delta t}+{b}_{i2}^{t+\Delta t}\right)-{{a}_{\mathrm{P}}}_{i2}^{t+\Delta t}{T}_{i,2}^t\right]}{{{a}_{\mathrm{P}}}_{i2}^{t+\Delta t}}, $$(57)j=3Ti,3t+Δt=Ti,3t+ω[(aSSi3t+ΔtTi+2,3t+Δt+aNi3t+ΔtTi-1,3t+Δt+aWi3t+ΔtTi,2t+Δt+aSi3t+ΔtTi+1,3t+Δt+aEi3t+ΔtTi,4t+Δt+bi3t+Δt)-aPi3t+ΔtTi,3t]aPi3t+Δt,$$ j=3\hspace{1em}{T}_{i,3}^{t+\Delta t}={T}_{i,3}^t+\frac{\omega \left[\left({{a}_{\mathrm{SS}}}_{i3}^{t+\Delta t}{T}_{i+\mathrm{2,3}}^{t+\Delta t}+{{a}_{\mathrm{N}}}_{i3}^{t+\Delta t}{T}_{i-\mathrm{1,3}}^{t+\Delta t}+{{a}_{\mathrm{W}}}_{i3}^{t+\Delta t}{T}_{i,2}^{t+\Delta t}+{{a}_{\mathrm{S}}}_{i3}^{t+\Delta t}{T}_{i+\mathrm{1,3}}^{t+\Delta t}+{{a}_{\mathrm{E}}}_{i3}^{t+\Delta t}{T}_{i,4}^{t+\Delta t}+{b}_{i3}^{t+\Delta t}\right)-{{a}_{\mathrm{P}}}_{i3}^{t+\Delta t}{T}_{i,3}^t\right]}{{{a}_{\mathrm{P}}}_{i3}^{t+\Delta t}}, $$(58)j4Ti,jt+Δt=Ti,jt+ω[(aNijt+ΔtTi-1,jt+Δt+aWijt+ΔtTi,j-1t+Δt+aSijt+ΔtTi+1,jt+Δt+aEijt+ΔtTi,j+1t+Δt+bijt+Δt)-aPijt+ΔtTi,jt]aPijt+Δt.$$ j\ge 4\hspace{1em}{T}_{i,j}^{t+\Delta t}={T}_{i,j}^t+\frac{\omega \left[\left({{a}_{\mathrm{N}}}_{{ij}}^{t+\Delta t}{T}_{i-1,j}^{t+\Delta t}+{{a}_{\mathrm{W}}}_{{ij}}^{t+\Delta t}{T}_{i,j-1}^{t+\Delta t}+{{a}_{\mathrm{S}}}_{{ij}}^{t+\Delta t}{T}_{i+1,j}^{t+\Delta t}+{{a}_{\mathrm{E}}}_{{ij}}^{t+\Delta t}{T}_{i,j+1}^{t+\Delta t}+{b}_{{ij}}^{t+\Delta t}\right)-{{a}_{\mathrm{P}}}_{{ij}}^{t+\Delta t}{T}_{i,j}^t\right]}{{{a}_{\mathrm{P}}}_{{ij}}^{t+\Delta t}}. $$(59)

Figure 3 shows the flow chart of model solution.

thumbnail Fig. 3

Flow chart of the model solution process.

3.3 Model coupling

Moreover, the continuity and momentum equations are decoupled from temperature prediction model. Consistent time step can be expressed as the following equation:Δt=Δzvg.$$ \Delta t=\frac{\Delta z}{{v}_{\mathrm{g}}}. $$(60)

When H2S is released from the drilling fluid, the state changes from supercritical to gas, which directly leads to a significant increase in gas void fraction. The gasification volume of H2S is influenced greatly by the temperature at a certain well depth. And the temperature also has a significant effect on the fluid properties (Xu et al., 2019). Therefore, before solving the mass and momentum equations, we must solve the transient temperature prediction model.

The well-killing simulation is based on the overflow condition at shut-in time; thus, the entire solution process can be divided into two steps: overflow and well-killing periods. For overflow period, the solution flowchart of overflow simulation is depicted in Figure 4. For well-killing period, the solution flowchart of well-killing simulation is depicted in Figure 5. As shown in the figure, depending on whether the gas reaches the wellhead, the solution process can be divided into two parts: gas rising and gas exhausting stages. The boundary positions of two-phase segment are key parameters in distinguishing the killing stage. The two-phase flow module used in the overflow period is also called to solve the two-phase flow behavior in the well-killing period.

thumbnail Fig. 4

Solution flowchart of overflow simulation for coupling model.

thumbnail Fig. 5

Solution flowchart of well-killing simulation.

4 Model validation

The basic data of an H2S-containing gas well Tiandong #5 were used to simulate the overflow process for verifying the validity of the dynamical well-killing model considering an H2S solubility. The overflow was found in Well Tiandong #5 at the drilling depth of 3570 m, and the pit gain was 1.65 m3 at 22:14. At 22:44, the wellhead emitted a loud noise, and a violent blowout accident occurred. The pit gain was recorded every 5 min, and the gas well was shut in when the pit gain reached 6 m3 at 22:44. Table 1 shows the basic parameters of Well Tiandong #5 (Blowout accident of TianDong #5 well, 1990). As shown in Figure 6a, the simulation results are in good agreement with the field data of Well Tiandong #5. When H2S is released, the pit gain shows a sudden increase in simulation results and field data. The slight difference between the two may be due to that the compressibility of the drilling fluid is disregarded.

thumbnail Fig. 6

Comparison results: (a) is the comparison between simulation and field data of Well Tiandong #5; (b) is the comparison between simulation and experimental results.

Table 1

Basic parameters of Well Tiandong #5.

Rader conducted a gas kick and well-killing experiment in an artificial well L.S.U. 7 and provided the experimental data (Rader et al., 1975). The basic data of experimental well L.S.U. 7 were used to simulate the killing process, as shown in Table 2, for validating the validity of the well-killing model. In the experiment, (1) nitrogen was injected into the well-bottom hole from the injection pipe until the pit gain was 1.59 m3, and gas was continued to inject during the shut-in operation; (2) injection of gas was stopped when the casing pressure reached the desired value; (3) well killing was started. Figure 6b shows the comparison between simulation and experimental results. As shown in Figure 6b, the simulation results of our model are in good agreement with the experimental results. The comparison between simulation results of our model and traditional “Gas column model” is also shown in Figure 6b. Large error exists in the simulation results of “Gas column model” because the friction and slip between gas and liquid are ignored.

Table 2

Basic data of experimental well L.S.U 7.

5 Case study

The basic parameters from Table 3 are used to simulate the killing process of an H2S-containing natural gas well in Sichuan. The number of spatial nodes is 60 and number of spatial nodes is 80 in the simulation. The wellbore configuration of this well is shown in Figure 7. When overflow behavior occurs, the invaded gas will push the drilling fluid with the same volume out of the wellbore, thereby leading to an increase in the pit gain. In drilling engineering, once the pit gain exceeds a threshold value, a shut-in operation must be performed. On the basis of field experience, the threshold value of pit gain is set as 3 m3 in this work.

thumbnail Fig. 7

Wellbore configuration of the well in Sichuan.

Table 3

Major computational parameters of the gas well in Sichuan.

Figure 8a shows the change in gas void fraction and solubility of H2S with the variation in well depth at shut-in time. As shown in the figure, gas rises to a location that is 1800 m below the wellhead, and the maximal gas void fraction is approximately 0.471. The flow pattern criterion proposed by Hasan and Kabir (1988) indicates that three flow patterns exist in the wellbore: bubble, slug, and single phase flows. Besides, the H2S solubility decreases as the well depth decreases. This phenomenon can be attributed to the decrease in annular temperature and pressure with the well depth, as shown in Figures 8b and 8c. When the solubility of H2S decreases to H2S concentration in the drilling fluid at about 500 m, the phase change of H2S occurs.

thumbnail Fig. 8

Change in gas void fraction, solubility of H2S, annular temperature and annular pressure with the variation in well depth at shut-in time: (a) Gas void fraction and solubility of H2S; (b) Annular temperature; (c) Annular pressure.

Figure 8b shows the drilling fluid temperature distribution in the drilling string and annulus at shut-in time. The annular temperature is always higher than that in the drilling string because the drilling fluid in the annulus is closer to the formation, thus, more heat can be gotten from the formation. Moreover, the annular temperature and geothermal temperature intersect at about 2250 m, which indicates that when the well depth is less than 2250 m, the drilling fluid in the annulus can be cooled by the formation. When the well depth is greater than 2250 m, the drilling fluid in the annulus can be heated by the formation. Figure 8c shows the change in annular pressure with the variation in well depth at shut-in time. As shown in Figure 8c, the annular pressure below 2000 m is lower than the formation pressure. This phenomenon can be attributed to an invasion of natural gas, which increases the differential pressure between bottom hole and formation. The increasing differential pressure will cause an increase in gas flow from formation to the wellbore and result in blowout accidents easily if not controlled.

Driller’s method is widely used in well killing in oil and gas fields (Feng et al., 2016a). Thus, the method is also used in the simulation of this work. Figure 9 shows the physical model of the killing process. During killing, the pressure balance relationship can be expressed as follows (Feng et al., 2016b):PpPb=Pa+Pla+Pma.$$ {P}_{\mathrm{p}}\le {P}_{\mathrm{b}}={P}_{\mathrm{a}}+{P}_{\mathrm{la}}+{P}_{\mathrm{ma}}. $$(61)

thumbnail Fig. 9

Physical model of the killing process.

Figure 10a shows the change in bottom hole pressure with the variation in time. As shown in the figure, the bottom hole pressure first decreases during the overflow period and then keeps constant (greater or equal to formation pressure) when well killing is conducted. Therefore, no gas flows into the wellbore from the formation during killing.

thumbnail Fig. 10

Change in bottom hole pressure and casing pressure the variation in time: (a) Bottom hole pressure; (b) Casing pressure.

With the change in hydrostatic pressure created by the fluid column, the casing pressure of the wellhead can be changed by adjusting the throttle opening to maintain the bottom hole pressure constant. Figure 10b shows the change in casing pressure with the variation in time. As shown in the figure, the entire killing process can be divided into six stages. Stage 1 is the initial stage of the killing process, which is determined by the flow pattern at shut-in time. Stage 2 indicates that the two-phase section in the wellbore is rising up toward the wellbore prior to reaching the wellhead. At this time, the flow pattern distribution in the wellbore is single phase flow, gas–liquid two-phase flow, and single phase flow in sequence. Moreover, the casing pressure increases with the increase in killing time. This phenomenon is attributed to the slippage and expansion of gas. The gas will expand as it rises up along the wellbore because the annular and formation temperatures and depth decrease (Sun et al., 2018). Thus, the hydrostatic pressure decreases accordingly and the casing pressure increases to keep the bottom hole pressure constant. Stage 3 indicates that the gas reaches the wellhead and single phase and gas–liquid two-phase flows exist in the wellbore. The total gas volume reaches a maximum at this point. Stage 4 indicates that all the gas has been pushed out of the wellbore and only single phase flow exists in the wellbore. This stage will last for some time until the weighted drilling fluid enters the annulus. Furthermore, the casing pressure keeps constant in Stage 4. Stage 5 indicates that the weighted drilling fluid enters the annulus and the hydrostatic pressure increases. Thus, the casing pressure decreases. Stage 6 indicates that the annulus is full of the weighted drilling fluid and the casing pressure decreases to 0. The well-killing operation ends.

Figure 10b also shows that, when all the gas is about to be pushed out of the well, H2S is released, and the curve drops slowly. This phenomenon is the result of casing pressure applied at the wellhead. Figure 8a shows that, when overflow behavior occurs, H2S will release at a location close to the wellhead if shut-in operation is not performed. However, H2S solubility decreases with the increase in pressure, and the casing pressure applied at the wellhead changes the pressure distribution in the wellbore (Fig. 11). As shown in Figure 11, the pressure of each location along the wellbore is higher than that of the overflow process during the killing process. Therefore, the H2S solubility of each location in the killing process increases compared with that of shut-in time. Although it reaches the wellhead, the H2S will not release until the H2S solubility is lower than its content in the drilling fluid. Therefore, when the casing pressure decreases, the H2S solubility decreases. As a result, the H2S transforms from a supercritical phase to a gas phase, gas volume expands suddenly, and gas void fraction increases abruptly. Thus, the casing pressure drops slowly.

thumbnail Fig. 11

Dynamical pressure distribution in the wellbore: (a) during gas kick; (b) during well killing.

5.1 Analysis under different H2S contents

Considering that the phase change of H2S occurs at low pressure and temperature and the phase change of H2S has a great effect on the pressure distribution (Sun et al., 2013), the diagrams of wellhead flow pattern when gas reaches wellhead during well killing under different working conditions are depicted to study the flows at wellhead.

Figure 12 shows the diagram of wellhead flow pattern when gas reaches wellhead during well killing with H2S contents of 0%, 2%, 4%, and 6%. As shown in the figure, the gas void fraction increases with the increase in H2S content. This phenomenon is caused by the increase in gasification volume of H2S as H2S content increases. When the well depth decreases, the annular pressure and annular temperature decrease, as shown in Figure 13. Besides, the higher the H2S content is, the lower the annular temperature will be as shown in Figure 13b. This is because, near the critical point of H2S, the viscosity of supercritical H2S decreases significantly (Vesovic et al., 1998). Therefore, in the upper well section, the convective heat transfer coefficient between the drilling fluid and the formation increases as viscosity of supercritical H2S decreases. Thus, more heat will be diffused from the drilling fluid in the annulus to the formation. As H2S content increases, this effect will be enhanced. Therefore, in the upper well section which is near the critical point of H2S, the annular temperature is lower at higher H2S content and, as a result, at the same well depth, the H2S solubility will be lower and more H2S will gasify.

thumbnail Fig. 12

Diagram of wellhead flow pattern when gas reaches wellhead during well killing with H2S contents of 0%, 2%, 4%, and 6%.

thumbnail Fig. 13

Change in annular pressure and temperature with the variation in well depth at 1300 s during well killing with H2S contents of 0%, 2%, 4%, and 6%: (a) Annular pressure; (b) Annular temperature.

Additional H2S will gasify at higher H2S content, and this condition changes the wellhead flow pattern. Churn flow appears at the wellhead with H2S content of 0%; however, when H2S is released, annular flow appears and the droplets in the bubble decrease as H2S content increases. Therefore, the hydrostatic pressure decreases with the increase in H2S content, and the casing pressure increases to keep the bottom hole pressure constant. Figure 14 shows that, when H2S is released, the casing pressure increases with H2S content. Thus, the annular pressure of the upper well section increases with the increase in H2S content, as shown in Figure 13a. However, H2S content has a small effect on entire wellbore pressure distribution, as shown in Figure 14.

thumbnail Fig. 14

Change in wellbore pressure and casing pressure with the variation in time and H2S contents of 0%, 2%, 4%, and 6%:

From comparison between Figures 8c and 13a, it can be concluded that well killing can rebuild the pressure distribution in the wellbore and stop the overflow behavior caused by the pressure differential between bottom hole and formation. However, for the stratum with low formation fracture pressure, the well-killing operation may cause annulus pressure to approach or exceed formation fracture pressure, thereby causing well leakage. This scenario should be highly regarded, especially for H2S-containing natural gas wells. Figure 13a shows that, when H2S is released, the annular pressure of the upper well section will exceed the formation fracture pressure. Moreover, a high H2S content indicates serious fracture of the stratum. Therefore, during the late period of well killing in H2S-containing natural gas wells, the casing pressure should be reduced appropriately to avoid stratum fracturing and prevent leakage accidents.

5.2 Analysis under different initial overflow volumes

Figure 15 shows the diagram of wellhead flow pattern when gas reaches wellhead during well killing with initial overflow volumes of 2, 4, and 8 m3. As shown in the figure, the gas void fraction increases with the increase in initial overflow volume because a large amount of gas enters the wellbore at high initial overflow volume. Moreover, the wellhead flow pattern with initial overflow volume of 2 m3 can be regarded as slug flow, whereas the annular flow appears at high initial overflow volumes of 4 and 8 m3. The gas phase becomes continuous at initial overflow volumes of 4 and 8 m3. The continuous phase volume increases as the initial overflow volume increases.

thumbnail Fig. 15

Diagram of wellhead flow pattern when gas reaches wellhead during well killing with initial overflow volumes of 2, 4, and 8 m3.

A high initial overflow volume indicates gas rises to a high location in the wellbore at shut-in time. When initial overflow volume reaches 8 m3, the gas enters the wellhead at shut-in time. This conclusion can also be drawn from Figure 16. Figure 16b shows the change in casing pressure with the variation in time and initial overflow volumes of 2, 4, and 8 m3. As shown in the figure, the surface casing pressure increases with the increase in initial overflow volume. This phenomenon can be attributed to the increase in hydrostatic pressure loss as gas volume increases. Thus, additional surface casing pressure should be applied to balance the bottom hole pressure at high initial overflow volume. Moreover, as shown in Figure 16b, the pressure change in the wellbore intensifies as initial overflow volume increases.

thumbnail Fig. 16

Change in wellbore pressure and casing pressure with the variation in time and initial overflow volumes of 2, 4, and 8 m3: (a) Wellbore pressure; (b) Casing pressure.

Besides, as the initial overflow volume increases, the gasification starting time of H2S increases and, the casing pressure is lower as shown in Figure 16a. This phenomenon can be attributed to the change of drilling fluid temperature. Figure 17a shows the change in annular temperature with the variation in well depth at 1300 s during well killing with different initial overflow volume. We can see that, the drilling fluid temperature is higher at higher initial overflow volume. This is because, as more gas enters the wellbore, more heat from the formation will be carried into the wellbore. Thus, the H2S solubility will be higher at the same well depth and less H2S will gasify. Therefore, when H2S gasifies, the casing pressure applied at wellhead is higher at lower initial overflow volume.

thumbnail Fig. 17

Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with initial overflow volumes of 2, 4, and 8 m3: (a) Annular temperature; (b) Annular pressure.

Figure 17b shows the change in annular pressure with the variation in well depth and initial overflow volumes of 2, 4, and 8 m3. As initial overflow volume increases, the gas expansion capacity is limited. Thus, pressure changed rapidly from 2 m3 to 4 m3, but slowly from 4 m3 to 8 m3. As shown in the figure, the annular pressure is high with high initial overflow volume. This phenomenon is caused by the increase in surface casing pressure as initial overflow volume increases. However, when initial overflow volume increases, the annular pressure may approach or even exceed the formation fracture pressure, which may lead to formation fracturing. Therefore, careful monitoring of the pit gain is necessary during the exploitation of H2S-containing natural gas wells. Once the pit gain exceeds a threshold value, the shut-in operation should be performed immediately to prevent further overflow. Furthermore, during the late period of well killing in H2S-containing natural gas wells with high initial overflow volume, the surface casing pressure should be decreased appropriately to avoid stratum fracturing and prevent leakage accidents. Early detection of gas kick should be more frequent to avoid severe overflow.

5.3 Analysis under different drilling fluid densities

Figure 18 shows the diagram of wellhead flow pattern when gas reaches wellhead during well killing with drilling fluid densities of 1.20, 1.25, and 1.30 g/cm3. As shown in the figure, the gas void fraction at wellhead increases with the increase in drilling fluid density. This phenomenon can be attributed to the decrease in casing pressure caused by the increase in drilling fluid density. The annular pressure is low at shut-in time with low drilling fluid density because the hydrostatic pressure is proportional to drilling fluid density. Therefore, to maintain the bottom hole pressure, the surface casing pressure should be higher at low drilling fluid density, as shown in Figure 19. During well killing, the annular pressure is high at low drilling fluid density, as shown in Figure 20a. The gas void fraction is high at high drilling fluid density because the gas expands violently at low pressure. Figure 18 shows that the drilling fluid density slightly affects the wellhead flow pattern, and annular flow is the common flow pattern under drilling fluid densities of 1.20, 1.25, and 1.30 g/m3.

thumbnail Fig. 18

Diagram of wellhead flow pattern when gas reaches wellhead during well killing with drilling fluid densities of 1.20, 1.25, and 1.30 g/cm3.

thumbnail Fig. 19

Change in wellbore pressure and casing pressure with the variation in time and drilling fluid densities of 1.20, 1.25, and 1.30 g/cm3: (a) Wellbore pressure; (b) Casing pressure.

thumbnail Fig. 20

Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with drilling fluid densities of 1.20, 1.25, and 1.30 g/cm3: (a) Annular pressure; (b) Annular temperature.

Figure 20b shows the change in annular temperature with the variation in well depth at 1300 s during well killing with different drilling fluid density. As shown in Figure 20b, for the deeper well section, the drilling fluid temperature increases with the decrease in drilling fluid density. This is because, when the same heat is obtained from the formation, the greater the drilling fluid density, the less the increase in temperature (Gao et al., 2017). For the upper well section, the drilling fluid temperature is closer to the formation temperature, the temperature differences decrease. Considering that H2S only gasifies in the upper well section, the effect of temperature difference caused by the change in drilling fluid density on gasification of H2S is small as shown in Figure 19b.

Furthermore, annular pressure of the upper well section will exceed the formation fracture pressure, and low drilling fluid density indicates serious fracture of the stratum. Thus, during the well killing in H2S-containing natural gas wells, higher drilling fluid density should be used to avoid stratum fracturing and prevent leakage accidents.

5.4 Analysis under different well-killing displacements

Figure 21 shows the diagram of wellhead flow pattern when gas reaches wellhead during well killing with displacements of 20, 30, and 40 L/s. As shown in the figure, the gas void fraction at wellhead increases slightly with the increase in displacement. The annular flow is the common wellhead flow pattern under displacements of 20, 30, and 40 L/s, and the number of droplets in the continuous gas decreases with the increase in displacement.

thumbnail Fig. 21

Diagram of wellhead flow pattern when gas reaches wellhead during well killing with displacements of 20, 30, and 40 L/s.

Figure 22a shows the change in annular temperature with the variation in well depth at 1300 s during well killing with different displacement. The lower the displacement is, the longer time for heat exchange with the formation will be. Thus, for the upper well section, the annular temperature is higher than formation, at higher displacement, the drilling fluid temperature is higher because lesser heat of drilling fluid will loss. For the deeper well section, the annular temperature is lower than formation, at higher displacement, the drilling fluid temperature is lower because lesser heat from formation will be exchanged to the drilling fluid. However, the effect of displacement on temperature distribution in annulus is small.

thumbnail Fig. 22

Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with displacements of 20, 30, and 40 L/s: (a) Annular temperature; (b) Annular pressure.

Figure 23 shows the change in wellbore pressure and casing pressure with the variation in time and displacements of 20, 30, and 40 L/s. As shown in the figure, the higher the displacement is, the faster the gas kick leaves wellbore and, the lower the casing pressure will be. Besides, H2S will gasify at an earlier time. This phenomenon can be attributed to the increased friction resistance due to the increase in displacement (Shi et al., 2014). Thus, the bottom hole pressure increases as displacement increases and, the casing pressure decreases as a result. Therefore, the annular pressure is lower at higher displacement in the upper well section as shown in Figure 22b. Besides, the annular pressure of the upper well section will exceed the formation fracture pressure, and low displacement indicates serious fracture of the stratum. Thus, during well killing in H2S-containing natural gas wells, higher well-killing displacement should be used to avoid stratum fracturing and prevent leakage accidents.

thumbnail Fig. 23

Change in wellbore pressure and casing pressure with the variation in time with displacements of 20, 30, and 40 L/s: (a) Wellbore pressure; (b) Casing pressure.

6 Conclusion

This work established a dynamical well-killing model considering an H2S solubility to simulate the well-killing process of a vertical H2S-containing natural gas well. The following important conclusions can be drawn:

  1. H2S will gasify near wellhead during well killing when casing pressure decreases. Near the critical point of H2S, the annular temperature decreases as H2S content increases and, the H2S solubility will be lower and more H2S will gasify. To balance the bottom hole pressure, when H2S releases, the casing pressure increases as H2S content increase.

  2. As initial overflow volume increases, the annular temperature, annular pressure and the casing pressure increase significantly. When H2S gasifies, the casing pressure applied at wellhead should be higher at lower initial overflow volume to balance bottom hole pressure.

  3. In the well-killing process, the annular pressure and temperature decrease as drilling fluid density increases and a lower casing pressure is needed for balancing bottom hole pressure. No significant effect of drilling fluid density on the gasification of H2S is found. As well-killing displacement increases, for the upper well section, the annular temperature increases while the annular pressure decreases. The casing pressure is lower at a higher displacement for higher friction resistance. Besides, H2S will gasify at an earlier time.

  4. When drilling for H2S-containing natural gas well, early detection of gas kick should be more frequent to avoid severe overflow. Besides, during well killing, higher displacement and density of drilling fluid should be considered to avoid stratum fracturing and prevent leakage accidents under the premise of meeting drilling requirements.

Acknowledgments

This work was supported by the National Key Research and Development Program of China (2019YFC0312303) and Sichuan Science and Technology Project (2019YFS0045).

References

  • Amaya-Gómez R., López J., Pineda H., Urbano-Caguasango D., Pinilla J., Ratkovich N., Muñoz F. (2019) Probabilistic approach of a flow pattern map for horizontal, vertical, and inclined pipes, Oil Gas Sci. Technol. - Rev. IFP Energies nouvelles 74, 67. [CrossRef] [Google Scholar]
  • Ansari A.M., Sylvester N.D., Sarica C., Shoham O., Brill J.P. (1994) A comprehensive mechanistic model for upward two phase flow in wellbores, SPE Prod. Facil. 5, 143–152. https://doi.org/10.2523/20630-MS. [CrossRef] [Google Scholar]
  • Bilicki Z., Kestin J. (1987) Transition criteria for two-phase flow patterns in vertical upward flow, Int. J. Multiphase Flow. 13, 3, 283–294. [CrossRef] [Google Scholar]
  • Blowout accident of TianDong well #5 (1990) Petrochina Southwest Oil and Gas field Company. [Google Scholar]
  • Ebrahimi A., Khamehchi E. (2015) A robust model for computing pressure drop in vertical multiphase flow, J. Nat. Gas. Sci. Eng. 26, 1306–1316. [Google Scholar]
  • Elsharkawy A.M. (2004) Efficient methods for calculations of compressibility, density and viscosity of natural gases, Fluid Phase Equilibr. 218, 1, 1–13. [CrossRef] [Google Scholar]
  • Feng J., Fu J., Chen P., Liu Z., Wei H. (2015) Predicting pressure behavior during dynamic kill drilling with a two-phase flow, J. Nat. Gas Sci. Eng. 22, 591–597. [Google Scholar]
  • Feng J., Fu J., Chen P., Luo J., Kou B. (2016a) Comparisons of the Driller’s method and the wait and weight method in deepwater well killing operation, Arab. J. Sci. Eng. 41, 7, 2699–2706. [Google Scholar]
  • Feng J., Fu J., Luo J., Liu Z., Wei H. (2016b) An advanced Driller’s Method simulator for deepwater well control, J. Loss Prevent. Proc. 39, 131–140. [CrossRef] [Google Scholar]
  • Fu J., Su Y., Jiang W., Li S., Chen Y. (2019) Wellbore annulus water hammer pressure prediction based on transient multi-phase flow characteristics, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 74, 84. [CrossRef] [Google Scholar]
  • Gao Y., Cui Y., Xu B., Sun B., Zhao X., Li H., Chen L. (2017) Two phase flow heat transfer analysis at different flow patterns in the wellbore, Appl. Therm. Eng. 117, 544–552. [Google Scholar]
  • Guo R., Chen Y., Waltrich P.J., Williams W.C. (2018) An experimental investigation on flow pattern map and drift-flux model for Co-Current upward liquid-gas two-phase flow in narrow annuli, J. Nat. Gas. Sci. Eng. 51, 65–72. [Google Scholar]
  • Guo X., Wang Q. (2016) A new prediction model of elemental sulfur solubility in sour gas mixtures, J. Nat. Gas Sci. Eng. 31, 98–107. [Google Scholar]
  • Hasan A.R., Kabir C.S. (1988) A study of multiphase flow behavior in vertical wells, Spede. 3, 2, 263–272. https://doi.org/10.2118/15138-PA. [Google Scholar]
  • Hou Z., Yan T., Li Z., Feng J., Sun S., Yuan Y. (2019) Temperature prediction of two phase flow in wellbore using modified heat transfer model: An experimental analysis, Appl Therm. Eng. 149, 54–61. [Google Scholar]
  • Kelessidis V.C., Dukler A.E. (1989) Modeling flow pattern transitions for upward gas-liquid flow in vertical concentric and eccentric annuli, Int. J. Multiphase Flow. 15, 2, 173–191. [CrossRef] [Google Scholar]
  • Li H.Z., Mouline Y. (1997) Chaotic bubble coalescence in non-newtonian fluids, Int. J. Multiphase Flow. 23, 713–723. [CrossRef] [Google Scholar]
  • Mao L., Zhang Z. (2018) Transient temperature prediction model of horizontal wells during drilling shale gas and geothermal energy, J. Petrol. Sci. Eng. 169, 610–622. [CrossRef] [Google Scholar]
  • Meng Y., Xu C., Na W., Gao L., Hongtao L., Mubai D. (2015) Numerical simulation and experiment of the annular pressure variation caused by gas kick/injection in wells, J Nat Gas Sci Eng. 22, 646–655. [Google Scholar]
  • Nickens H.V. (1987) A dynamic computer model of kick Well, SPE Drill. Eng. 2, 2, 158–173. https://doi.org/10.2118/14183-PA. [Google Scholar]
  • Osman E.S.A. (2004) Artificial neural network models for identifying flow regimes and predicting liquid holdup in horizontal multiphase flow, SPE Prod. Facil. 19, 1, 33–40. https://doi.org/10.2118/68219-MS. [CrossRef] [Google Scholar]
  • Rader D.W., Bourgoyne A.T., Ward R.H. (1975, May 1) Factors affecting bubble-rise velocity of gas kicks, Society of Petroleum Engineers. https://doi.org/10.2118/4647-PA. [Google Scholar]
  • Raghavan R. (1989) Well-test analysis for multiphase flow, SPE Form. Eval. 4, 4, 585–594. https://doi.org/10.2118/14098-PA. [CrossRef] [Google Scholar]
  • Shi H., Li G., Huang Z., Shi S. (2014) Properties and testing of a hydraulic pulse jet and its application in offshore drilling, Petrol. Sci. 11, 3, 401–407. [CrossRef] [Google Scholar]
  • Sun B., Gong P.B., Wang Z.Y. (2013) Simulation of gas kick with high H2S content in deep well, J. Hydrodyn. 25, 2, 264–273. [CrossRef] [Google Scholar]
  • Sun B.J., Sun X.H., Wang Z.Y., Chen Y.H. (2017) Effects of phase transition on gas kick migration in deepwater horizontal drilling, J. Nat. Gas Sci. Eng. 46, 710–729. [Google Scholar]
  • Sun B.J., Guo Y.H., Sun W.T., Gao Y., Hao L., Wang Z., Zhang H. (2018) Multiphase flow behavior for acid-gas mixture and drilling fluid flow in vertical wellbore, J. Petrol. Sci. Eng. 165, 388–396. [CrossRef] [Google Scholar]
  • Vesovic V., Assael M.J., Gallis Z.A. (1998) Prediction of the viscosity of supercritical fluid mixtures, Int. J. Thermophys. 19, 5, 1297–1313. [Google Scholar]
  • Wang N., Sun B., Wang Z., Wang J., Yang C. (2016) Numerical simulation of two phase flow in wellbores by means of drift flux model and pressure based method, J. Nat. Gas. Sci. Eng. 36, 811–823. [Google Scholar]
  • Wang Z.Y., Sun B.J. (2014) Deepwater gas kick simulation with consideration of the gas hydrate phase transition, J. Hydrodyn. 26, 1, 94–103. [CrossRef] [Google Scholar]
  • Wang Z., Sun B. (2009) Annular multiphase flow behavior during deep water drilling and the effect of hydrate phase transition, Petrol. Sci. 6, 1, 57–63. [CrossRef] [Google Scholar]
  • Wei N., Xu C.Y., Meng Y.F., Li G. (2018) Numerical simulation of gas-liquid two-phase flow in wellbore based on drift flux model, Appl. Math. Comput. 338, 175–191. [Google Scholar]
  • Xu Z.M., Song X.Z., Li G.S., Wu K., Pang Z., Zhu Z. (2018) Development of a transient non-isothermal two-phase flow model for gas kick simulation in HTHP deep well drilling, Appl. Therm. Eng. 141, 1055–1069. [Google Scholar]
  • Xu Z.M., Song X.Z., Li G.S., Zhu Z., Zhu B. (2019) Gas kick simulation in oil-based drilling fluids with the gas solubility effect during high-temperature and high-pressure well drilling, Appl. Therm. Eng. 149, 1080–1097. [Google Scholar]
  • Xie J., Zhang X., Tang Y., Wang Y., Shao Q., Yu B. (2014) Transient simulation of the blowing-out process of the air pockets in vertical wellbore, Appl. Therm. Eng. 72, 97–103. [Google Scholar]
  • Yin H., Liu P., Li Q., Wang Q., Gao D. (2015) A new approach to risk control of gas kick in high-pressure sour gas wells, J. Nat. Gas Sci. Eng. 26, 142–148. [Google Scholar]
  • Yin B., Gang L., Li X. (2017) Multiphase transient flow model in wellbore annuli during gas kick in deepwater drilling based on oil-based mud, Appl. Math. Model. 51, 159–198. [Google Scholar]
  • Zhang X., Li Q., Zheng L., Li X., Xu L. (2020) Numerical simulation and feasibility assessment of acid gas injection in a carbonate formation of the Tarim Basin, China. Oil Gas, Sci. Technol. 75, 1–17. [Google Scholar]
  • Zhang J., Li X., Tang X., Luo W. (2019) Establishment and analysis of temperature field of riserless mud recovery system, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 74, 19, 1–9. [CrossRef] [Google Scholar]
  • Zhu G., Zhang S., Huang H., Liang Y., Meng S., Li Y. (2011) Gas genetic type and origin of hydrogen sulfide in the Zhongba gas field of the western Sichuan Basin, China, Appl. Geochem. 26, 7, 1261–1273. [Google Scholar]

All Tables

Table 1

Basic parameters of Well Tiandong #5.

Table 2

Basic data of experimental well L.S.U 7.

Table 3

Major computational parameters of the gas well in Sichuan.

All Figures

thumbnail Fig. 1

Schematic diagram of solution for the mass and momentum equations: (a) Schematic of discrete grid nodes of the wellbore; (b) Diagram of cell grid integration area K.

In the text
thumbnail Fig. 2

The grid diagram of the wellbore and formation.

In the text
thumbnail Fig. 3

Flow chart of the model solution process.

In the text
thumbnail Fig. 4

Solution flowchart of overflow simulation for coupling model.

In the text
thumbnail Fig. 5

Solution flowchart of well-killing simulation.

In the text
thumbnail Fig. 6

Comparison results: (a) is the comparison between simulation and field data of Well Tiandong #5; (b) is the comparison between simulation and experimental results.

In the text
thumbnail Fig. 7

Wellbore configuration of the well in Sichuan.

In the text
thumbnail Fig. 8

Change in gas void fraction, solubility of H2S, annular temperature and annular pressure with the variation in well depth at shut-in time: (a) Gas void fraction and solubility of H2S; (b) Annular temperature; (c) Annular pressure.

In the text
thumbnail Fig. 9

Physical model of the killing process.

In the text
thumbnail Fig. 10

Change in bottom hole pressure and casing pressure the variation in time: (a) Bottom hole pressure; (b) Casing pressure.

In the text
thumbnail Fig. 11

Dynamical pressure distribution in the wellbore: (a) during gas kick; (b) during well killing.

In the text
thumbnail Fig. 12

Diagram of wellhead flow pattern when gas reaches wellhead during well killing with H2S contents of 0%, 2%, 4%, and 6%.

In the text
thumbnail Fig. 13

Change in annular pressure and temperature with the variation in well depth at 1300 s during well killing with H2S contents of 0%, 2%, 4%, and 6%: (a) Annular pressure; (b) Annular temperature.

In the text
thumbnail Fig. 14

Change in wellbore pressure and casing pressure with the variation in time and H2S contents of 0%, 2%, 4%, and 6%:

In the text
thumbnail Fig. 15

Diagram of wellhead flow pattern when gas reaches wellhead during well killing with initial overflow volumes of 2, 4, and 8 m3.

In the text
thumbnail Fig. 16

Change in wellbore pressure and casing pressure with the variation in time and initial overflow volumes of 2, 4, and 8 m3: (a) Wellbore pressure; (b) Casing pressure.

In the text
thumbnail Fig. 17

Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with initial overflow volumes of 2, 4, and 8 m3: (a) Annular temperature; (b) Annular pressure.

In the text
thumbnail Fig. 18

Diagram of wellhead flow pattern when gas reaches wellhead during well killing with drilling fluid densities of 1.20, 1.25, and 1.30 g/cm3.

In the text
thumbnail Fig. 19

Change in wellbore pressure and casing pressure with the variation in time and drilling fluid densities of 1.20, 1.25, and 1.30 g/cm3: (a) Wellbore pressure; (b) Casing pressure.

In the text
thumbnail Fig. 20

Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with drilling fluid densities of 1.20, 1.25, and 1.30 g/cm3: (a) Annular pressure; (b) Annular temperature.

In the text
thumbnail Fig. 21

Diagram of wellhead flow pattern when gas reaches wellhead during well killing with displacements of 20, 30, and 40 L/s.

In the text
thumbnail Fig. 22

Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with displacements of 20, 30, and 40 L/s: (a) Annular temperature; (b) Annular pressure.

In the text
thumbnail Fig. 23

Change in wellbore pressure and casing pressure with the variation in time with displacements of 20, 30, and 40 L/s: (a) Wellbore pressure; (b) Casing pressure.

In the text

Les statistiques affichées correspondent au cumul d'une part des vues des résumés de l'article et d'autre part des vues et téléchargements de l'article plein-texte (PDF, Full-HTML, ePub... selon les formats disponibles) sur la platefome Vision4Press.

Les statistiques sont disponibles avec un délai de 48 à 96 heures et sont mises à jour quotidiennement en semaine.

Le chargement des statistiques peut être long.