Systematic Oil Flow Modeling in the Quasi-3D Approximation Yields Additional Terms that Allows for Variable Cross-Section Area Tubing
SSystematic Oil Flow Modeling in the Quasi-3D Approximation Yields AdditionalTerms that Allows for Variable Cross-Section Area Tubing
Edval J. P. Santos
Laboratory for Devices and Nanostructures – Nanoscale Engineering Group,Universidade Federal de Pernambuco, Recife-PE, BrazilE-mail: [email protected] or [email protected].
A systematic model development for oil flow in quasi-3D (1D + 2D) is presented. Our approachprovides a unified modeling scheme. Besides, additional terms are obtained, which allows for tubingarea variation along the flow direction. The area variation can be modeled as analytic function orrandom fluctuation, as it could be the result of deposits or tubing internal surface roughness. Theproposed approach can be used to obtain analytic solutions which provide physical insight into thephenomena under scrutiny, including the validation of software tools, sensor development and sensorplacement. One starts from conservation laws as given by kinetic theory and applies the transverseaveraging technique (TAT) to extract the one-dimensional approximation in formal grounds. Todemonstrate its application, the steady-state Ramey’s model, the Hasan’s transient model and asimple two-phase model are generated from the obtained equations.
Keywords: Conservation equation, Flow model, Two-phase flow, Ramey’s model, Hasan’s model.
I. INTRODUCTION
Many models have been developed to describe oil and gas flow in production wells or injection fluid in artificiallift wells in terms of flowrate, pressure profile, and temperature profile. A pioneer work in such development is theone dimensional phenomenological model for steady-state flow in injection wells reported by Henry J. Ramey Jr. [1]in 1962. Later this model was improved for transient flow. With the advancement of computers, 3D models havebeen developed and implemented as computer software. Despite the sophistication of computer tools, analytic modelsstills plays a role to offer insight into physical problems, to help in tool validation, sensor development and sensorplacement [2–15].In the oil and gas industry, balance of mass, momentum or energy are applied to a fluid parcel volume to build modelsfor the extracted fluid flow in production wells or the injection fluid in artificial lift wells. The model equations arecombined with equations of state and empirical equations. Depending on the model, many other ancillary parametersneed to be measured or modeled [16]. Full 3D coupled models are used during well drilling to monitor pressure surgeor swab to keep the drilling process within the safe window to avoid well fracture or blowout. During production, themodel is used for keeping the well in production, deciding about treatment, or deciding about shut-in for maintenance,safety precaution or economical reasons [1–10].Extracted fluid is a combination of many components, such as: oil, water, gas, and sediments. Depending on theviscosity of the oil, different viscosity models apply, such as: power law fluid, Bingham fluid, yield power law fluid,Casson fluid, among others. Range of crude oil density [17] is shown in Table I.
TABLE I: Crude oil density.API gravity Kg/m Light > < < > The simplest model of the extracted fluid is a mixture of ideal gas and incompressible Newtonian liquid. A two-phaseflow. Considering the two-phase mixture, the simplest approximation for the flow density, ρ F , is a linear combination.Thus, the density of the fluid can be estimated as follows. ρ F = ρ L H L + ρ G (1 − H L ) (1)in which, H L is the liquid fraction also known as liquid holdup.An equivalent expression can be written for the fluid velocity. ρ F (cid:126)u = ρ L (cid:126)u L H L + ρ G (cid:126)u G (1 − H L ) (2) a r X i v : . [ phy s i c s . g e o - ph ] N ov In the simplest model, liquid and gas velocities are equal. This is the basis of the Homogeneous Equilibrium Model(HEM). Unequal liquid and gas velocities results in phase sliding. As oil and gas move up to the wellhead as afunction of pressure profile, it exchanges thermal energy with the surroundings, resulting in a variable temperatureprofile in the wellbore. Variable temperature affects oil and gas density, and many other properties of the fluid. Inreturn, this causes a variation of the pressure profile, and the flow regime. One says, the flow equations describingsuch phenomena are coupled. Besides, along the tubing, the flow regime can vary, ranging from plug flow to laminarflow to turbulent flow. Plug flow is a single phase smooth flow. A characteristic of two-phase flow is the presence ofone or more interfaces separating the phases, a few examples of flow regime with approximate cylindrical symmetryare [18, 19]: • Stratified – at low gas and liquid flowrates, phase separated, multicomponent smooth flow. • Stratified/wavy – at higher gas and liquid flowrates, alternates between liquid and gas. • Slug – at medium gas and liquid flowrates, large bubbles. • Bubbly – at high liquid flowrates, bubbles of gas in the liquid phase. • Annular – at high gas flowrates, gas with liquid wetting the pipe wall.To solve the full 3D coupled problem of a multiphase flow is a difficult task even for multiphysics tools. Approximatemethods based on averaging techniques are useful to offer insight and to validate software tools. There are manyapproximate methods based on area averaging, such as: Homogeneous Equilibrium Model (HEM), Separated FlowModel (SFM), Drift Flux Model (DFM), and Two Fluid Model (TFM). In the Homogeneous Equilibrium Modelis assumed that fluid phases are in equilibrium with each other, so that any differences in velocity, pressure ortemperature rapidly disappear resulting in a single phase pseudo fluid. In the Separated Flow Model, the phases donot necessarily have the same velocity. Thus, there is a relative or drift velocity of the lighter phase over the heavierone, i.e, phase slip can occur with equal pressure and temperature. The Drift Flux Model is a particular case ofthe separated flow model in which the drift velocity is written as a function of external forces and fluid properties.The Two Fluid Model is the most complex model, each phase can achieve its own velocity, pressure and temperatureprofile, and is treated separately [18, 19].In kinetic theory, conservation laws play a key role in describing system behavior. Frequently, conservation ofparticles, conservation of momentum and conservation of energy equations are used to provide insight into manyphysical problems. Phenomena ranging from flow of charge carriers in a semiconductor device to flow of oil and gasin a well can be described within such theoretical framework [18–22].In a previous work, we have demonstrated how to break the problem into 1D + 2D dimensions by applying thetransverse averaging technique (TAT), and we have used this approximation to study the nonuniformly doped orshaped p - n junction diodes. In the occasion the TAT technique was applied to the conservation of particles equationfor holes and electrons [21, 22]. In TAT, the cross-section area is not necessarily held constant. In this paper theTAT technique is extended to tensor quantities and applied to the conservation of mass, conservation of momentumand conservation of energy equations from kinetic theory. Thus, the one-dimensional approximation of the threeequations is set on formal grounds without assuming constant cross-section area. As an application example, a setof equations is obtained to describe flow in an oil well. The area variation can be modeled as analytic function orrandom fluctuation, as it could be the result of deposits or tubing internal surface roughness.The goal of this approach is not to replace sophisticated models, but to provide a systematic model constructionapproach in 1D + 2D which can be used to provide insight into oil flow or to test software tools. In our approach, theconservation equations are simplified with the transverse averaging technique (TAT) to provide equations for oil flow.Under appropriate assumptions, the resulting equations are equivalent to equations obtained by studying a controlparcel volume [1, 2]. The paper is divided into six sections. This introduction is the first. Next, the conservationequations in fluids are discussed. In the third section, the TAT approximation is applied. In the fourth section, the setof equations for oil flow is obtained. In the fifth section, application examples are presented. Finally, the conclusions. II. CONSERVATION EQUATIONS IN FLUIDS
To build a mathematical model for the flow, the fluid is divided into small parts or parcels with identical mass, m , which are small enough to justify the continuum approximation and large enough as compared to the molecules.Thus, the exact number of parcels, N , is not well defined [18, 19, 23–25].The distribution of parcels or particles in the phase space is given by the probability density function, also nameddistribution function or particle density function, F B ( (cid:126)r, (cid:126)v, t ) = N f B ( (cid:126)r, (cid:126)v, t ). It is the probability density of finding aparcel at position between (cid:126)r and (cid:126)r + d(cid:126)r , with velocity between (cid:126)v and (cid:126)v + d(cid:126)v , at instant t . From the probability densityfunction one can extract properties of the system, such as fluid density, ρ F , momentum density, ρ F (cid:126)u , total pressure, P ,total flux of kinetic energy density, Φ q , and total kinetic energy density, e K . This is also named Boltzmann statisticalaverage [18–20]. ρ F ( (cid:126)r, t ) = (cid:90) m F B ( (cid:126)r, (cid:126)v, t ) d(cid:126)v (3) ρ F ( (cid:126)r, t ) (cid:126)u ( (cid:126)r, t ) = (cid:90) m(cid:126)vF B ( (cid:126)r, (cid:126)v, t ) d(cid:126)v (4) P ( (cid:126)r, t ) = (cid:90) m vv F B ( (cid:126)r, (cid:126)v, t ) d(cid:126)v (5) e K ( (cid:126)r, t ) = (cid:90) m ( (cid:126)v · (cid:126)v ) F B ( (cid:126)r, (cid:126)v, t ) d(cid:126)v (6)Φ q ( (cid:126)r, t ) = (cid:90) m ( (cid:126)v · (cid:126)v ) (cid:126)vF B ( (cid:126)r, (cid:126)v, t ) d(cid:126)v (7)in which, (cid:126)u is the parcel average velocity, P is the pressure tensor, and vv = (cid:126)v (cid:126)v T is the tensor dyadic product.The parcel velocity, (cid:126)v , can be written in terms of the average velocity, (cid:126)u , which is the fluid velocity. (cid:126)v = (cid:126)u + (cid:126)c (8)in which, (cid:126)c , is the velocity of the parcel in a reference frame moving with the fluid at the average velocity.In general, the dynamics of a set of classical particles is described by Boltzmann equation, from which one derivesconservation equations for particles, momentum, and energy [20]. • Conservation of mass density. ∂ρ F ∂t + ∇ · ( ρ F (cid:126)u ) = G N (9)in which, G N = G − R is the net generation, G is related to the source or creation of parcels, and R is relatedto the sink or destruction of parcels. • Conservation of momentum density. ∂ ( ρ F (cid:126)u ) ∂t + ∇ · P = (cid:126)f ext (10)in which, P is the total pressure tensor or total flux of momentum density. • Conservation of energy density. ∂e K ∂t + ∇ · Φ q = (cid:126)f ext · (cid:126)u (11)in which, e K is the total kinetic energy density, and Φ q is the total flux of kinetic energy density.Rewriting the velocity in terms of the average, as shown in Equation 8, and inserting in Equations 5, 6, and 7, onegets the following expressions. P = p + ρ F uu (12) e K = UV + 12 ρ F ( (cid:126)u · (cid:126)u ) (13)Φ q = Φ Q + e K (cid:126)u + p · (cid:126)u (14)in which, UV = nk B T . A. Conservation of particles
Defining the material or convective derivative,
DDt ≡ ∂∂t + (cid:126)u · ∇ , also known as the Eulerian operator. The materialderivative is the time variation of the quantity of interest as the fluid moves along a trajectory. Considering theequation of conservation of particles, and using the identity ∇ · ( ρ F (cid:126)u ) = ρ F ∇ · (cid:126)u + (cid:126)u · ( ∇ ρ F ), one can rewrite asfollows [23]. ∂ρ F ∂t + ∇ · ( ρ F (cid:126)u ) = Dρ F Dt + ρ F ∇ · (cid:126)u = G N (15)in which, G N is the net generation term. It can be used to model the reservoir. Above the reservoir, G N = 0.The material derivative is related to the creation or destruction of mass. For an incompressible fluid, ∇ · (cid:126)u = 0.Above the reservoir the flow can be assumed free of sink or sources, R = 0, G = 0, the material derivative is zero. Dρ F Dt = 0 (16) B. Conservation of momentum
Replacing Equation 12 in the equation for the conservation of momentum, Equation 10. ∂ ( ρ F (cid:126)u ) ∂t + ∇ · ρ F uu = −∇ · p + (cid:126)f ext (17)in which, uu = (cid:126)u (cid:126)u T is the tensor dyadic product, and ∇ · p (cid:12)(cid:12) i = (cid:80) j ∂p ij ∂x j .Replacing the identity ∇ · ( ρ F (cid:126)u(cid:126)u T ) = (cid:126)u ( ∇ · ρ F (cid:126)u ) + ( ρ F (cid:126)u · ∇ ) (cid:126)u , one gets a new equation for the conservation ofmomentum. ∂ ( ρ F (cid:126)u ) ∂t + (cid:126)u ( ∇ · ρ F (cid:126)u ) + ( ρ F (cid:126)u · ∇ ) (cid:126)u = −∇ · p + (cid:126)f ext (18)The pressure tensor can be replaced by the Cauchy stress tensor, σ , as follows [25]. − p = σ = σ τ τ τ σ τ τ τ σ (19)in which, σ ii are the normal stresses and τ ij are the shear or viscous stresses.If the external force per unit volume acting on the fluid is gravity, then (cid:126)f ext = ρ F (cid:126)g . Combining the conservation ofmomentum equation with the conservation of mass. ρ F ∂(cid:126)u∂t + ( ρ F (cid:126)u · ∇ ) (cid:126)u + (cid:126)u ( G − R ) = ∇ · σ + ρ F (cid:126)g (20)For a fluid free of sink or sources, R = 0, G = 0, one gets the Cauchy momentum equation. ρ F ∂(cid:126)u∂t + ( ρ F (cid:126)u · ∇ ) (cid:126)u = ∇ · σ + ρ F (cid:126)g (21)For a Newtonian fluid, the stress tensor is a combination of hydrostatic pressure and shear stress proportional tothe shear rate. σ ij = − pδ ij + β ( ∇ · (cid:126)u ) δ ij + η (cid:18) ∂u i ∂x j + ∂u j ∂x i (cid:19) (22)in which, β is the second viscosity coefficient and η is the shear viscosity.For an incompressible Newtonian fluid, ∇ · (cid:126)u = 0, the Cauchy momentum equation becomes the Navier-Stokesequations. ρ F ∂(cid:126)u∂t + ( ρ F (cid:126)u · ∇ ) (cid:126)u = −∇ p + η ( ∇ · ∇ ) (cid:126)u + ρ F (cid:126)g (23)For hydrostatic pressure, σ ij = − pδ ij . For an ideal fluid with constant density and no body forces, and replacingthe Eulerian operator, Equation 21 reduces to the Euler equation. ρ F D(cid:126)uDt = −∇ p (24) C. Conservation of energy
Similarly, with the expressions in Equations 12, 13, and 14, one can re-write the conservation of energy equation [20]. ∂U/V∂t + ∇ · ( (cid:126)uU/V ) + ∇ · Φ Q + ∂ ρ F ( (cid:126)u · (cid:126)u ) ∂t + ∇ · (cid:18) ρ F ( (cid:126)u · (cid:126)u ) (cid:19) (cid:126)u + ∇ · ( p · (cid:126)u ) = (cid:126)f ext · (cid:126)u (25)Replacing Equation 10. Using the identity ∇ · ( p · (cid:126)u ) = ( ∇ · p ) · (cid:126)u + p : ( ∇ ⊗ (cid:126)u ) T = ( ∇ · p ) · (cid:126)u + p : ∇ (cid:126)u and othervector identities. ∂U/V∂t + ∇ · ( (cid:126)uU/V ) + ∇ · Φ Q + p : ∇ (cid:126)u + ∂ ρ F ( (cid:126)u · (cid:126)u ) ∂t + ∇ · (cid:18) ( (cid:126)u · (cid:126)u ) 12 ρ F (cid:126)u (cid:19) = (cid:18) ∂ ( ρ F (cid:126)u ) ∂t + (cid:126)u ( ∇ · ρ F (cid:126)u ) + ( ρ F (cid:126)u · ∇ ) (cid:126)u (cid:19) · (cid:126)u (26)in which, p : ∇ (cid:126)u = (cid:80) ij p ij ∂u i /∂x j is a scalar. ∂U/V∂t + ∇ · ( (cid:126)uU/V ) + ∇ · Φ Q + p : ∇ (cid:126)u =12 (cid:126)u · (cid:126)u (cid:18) ∂ρ F ∂t + ∇ · ( ρ F (cid:126)u ) (cid:19) (27)Replacing Equation 9. ∂U/V∂t + ∇ · ( (cid:126)uU/V ) + ∇ · Φ Q + p : ∇ (cid:126)u =12 ( (cid:126)u · (cid:126)u ) G N (28)The oil well can be regarded as a thermodynamic system. The well plus surroundings can be considered a closedsystem, but the well alone is an open system, which exchanges thermal energy, mass and work with its surroundings,as shown in Figure 1. Unfortunately, the internal energy, U , is not a practical quantity. Good practical candidatesare the enthalpy, H , and the Gibbs free energy, G as they are functions of temperature variation, pressure variation,and mass exchange. Thus, one needs to re-write Equation 28. FIG. 1: The oil well as a thermodynamic closed system. dH = T dS + V dp + (cid:88) j µ j dN j (29) dG = − SdT + V dp + (cid:88) j µ j dN j (30)Experimentally, it is easier to determine the heat capacity at constant pressure, C p . Thus, it is more convenient touse the enthalpy. Writing the entropy, S , as a function of temperature and pressure and applying Maxwell identities. dH = ∂H∂T (cid:12)(cid:12)(cid:12)(cid:12) p,N dT + ∂H∂p (cid:12)(cid:12)(cid:12)(cid:12) T,N dp + (cid:88) j ∂H∂N j (cid:12)(cid:12)(cid:12)(cid:12) p,T dN j (31) ∂H∂T (cid:12)(cid:12)(cid:12)(cid:12) p,N = C p (32) ∂H∂p (cid:12)(cid:12)(cid:12)(cid:12) T,N = − µ JT C p (33)in which, µ JT = ∂T∂p (cid:12)(cid:12)(cid:12) H is the Joule-Thomson coefficient.It can be shown that the Joule-Thomson coefficient is related to the coefficient of thermal expansion. µ JT = 1 C p (cid:32) T ∂V∂T (cid:12)(cid:12)(cid:12)(cid:12) p,N − V (cid:33) = α p T − ρ F c p (34)in which, α p = V ∂V∂T (cid:12)(cid:12) p,N is the coefficient of thermal expansion at constant pressure, and c p = C p /m F is the specificheat capacity.One can define the enthalpy density to evaluate the exchange of thermal energy. H V = H/V = U/V + p + (cid:88) j µ j n j (35)in which, n j = N j /V . FIG. 2: Wellbore profile to apply the transverse averaging technique.
Dividing the enthalpy density by the fluid density, one gets the specific enthalpy, which is the enthalpy per unitmass.This is also known as the Local Instant Formulation combined with the Boltzmann Statistical Average [18]. Thethree conservation equations are coupled and must be solved simultaneously. Unfortunately, this is not an easy task.Thus, empirical models and assumptions are introduced to reduce the level of coupling.
III. QUASI 3D CASE - TRANSVERSE AVERAGING TECHNIQUE
Defining the average of a scalar, φ ( (cid:126)r, z ), and the average of a vector, (cid:126)V ( (cid:126)r, z ), over the cross-section area for a axiallysymmetric wellbore [21, 22]. Differently from the local instant formulation [18], there is no need to assume the areais constant along the longitudinal direction. φ ( z ) = 1 A ( z ) (cid:90) A ( z ) φ ( (cid:126)r, z ) dS (36) V ( z ) = 1 A ( z ) (cid:90) A ( z ) (cid:126)V ( (cid:126)r, z ) dS (37)It is more convenient to use the normal, ˆ n ( z ), to the profile curve, r ( z ), as shown in Figure 2.tan α = dr ( z ) dz = ⇒ cos α = 1 (cid:112) dr ( z ) /dz ) (38)Thus, one gets the following identities for the gradient and divergence average [21, 22]. (cid:90) A ( z ) ∇ φ dS = ∂∂z (cid:2) ˆ e z φ ( z ) A ( z ) (cid:3) + (cid:73) L ( z ) φ ˆ n cos α dl (39) (cid:90) A ( z ) ∇ · (cid:126)V dS = ∂∂z (cid:2) ˆ e z · V ( z ) A ( z ) (cid:3) + (cid:73) L ( z ) (cid:126)V · ˆ n cos α dl (40)The transverse averaging technique can be applied to the divergence of rank-2 tensors. To achieve this, one canstart with the divergence of a vector, shown in Equation 40, to obtain the following expression. (cid:90) A ( z ) ∇ · σ (cid:12)(cid:12) i dS = ∂∂z [ˆ e z · σ | i A ( z )]+ (cid:73) L ( z ) σ (cid:12)(cid:12) i · ˆ n cos α dl (41)Using such definitions, the conservation of particles equation, the conservation of momentum equation, and theconservation of energy equation can be re-written to calculate flow rate, pressure profile and temperature profile alongan oil well. Next, the equations are obtained for a single phase in the oil flow. A. Conservation of particles
For laminar flow, one can assume that ( ρ F (cid:126)u ) · ˆ n = 0, i.e., there is no flow in the perpendicular direction to thetubing. Integrating the equation of conservation of particles over the cross-section area, one can apply the transverseaveraging technique, as defined in Equation 36 and 37, and use the identity in Equation 40. ∂ ( ρ F A ( z )) ∂t + ∂ ( ρ F u z A ( z )) ∂z + (cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:58) (cid:73) L ( z ) ( ρ F (cid:126)u ) · ˆ n cos α dl = G N A ( z ) (42)The source, G , can be used to model the oil flow from the reservoir at the well bottom, and the sink, R , can beused to model the well top, as the oil is stored in a tank. Should be desired, this equation can be written in terms ofthe volumetric flow rate, Q V = u z A or mass flow rate, Q M = ρ F u z A .Expanding the average of the product in terms of the variation around the average, and noting that δ u z = δρ F = 0. ρ F u z = ρ F u z + δρ F δ u z (43)The simplest approximation is to neglect second order effects, thus ρ F u z ≈ ρ F u z . ∂ ( ρ F A ( z )) ∂t + ∂ ( ρ F u z A ( z )) ∂z = G N A ( z ) (44)Assuming the cross-section area is constant along the z -direction. ∂ρ F ∂t + u z ∂ρ F ∂z + ρ F ∂ u z ∂z = G N (45) Dρ F Dt + ρ F ∂ u z ∂z = G N (46)Defining compressibility at constant temperature, β T . β T = 1 ρ F Dρ F Dp = − V ∂V∂p (cid:12)(cid:12)(cid:12)(cid:12) T (47)The density variation can be replaced with the pressure variation. β T DpDt + ∂ u z ∂z = 1 ρ F G N (48)Using the bulk modulus definition, K b = 1 /β T , finally, one gets. DpDt + c a ρ F ∂ u z ∂z = c a G N (49)in which, c a = K b /ρ F is the square of the speed of sound. B. Conservation of momentum
Calculating the transverse averaging of Equation 18. ∂ρ F u z A ( z ) ∂t + ∇ · ( ρ F (cid:126)u(cid:126)u T ) (cid:12)(cid:12)(cid:12) z A ( z ) = ∇ · σ (cid:12)(cid:12)(cid:12) z A ( z ) + f ext,z A ( z ) (50)Assuming there is no momentum flow in the radial direction, and that the radial stress component is assumedconstant along the perimeter. (cid:73) L ( z ) σ (cid:12)(cid:12) z · ˆ n cos α dl = τ rz cos α πr ( z ) (51)Thus. ∂ρ F u z A ( z ) ∂t + ∂∂z [ˆ e z · ρ F uu | z A ( z )] + (cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:58) (cid:73) L ( z ) ρ F (cid:126)u(cid:126)u T (cid:12)(cid:12) z · ˆ n cos α dl = ∂∂z (cid:2) ˆ e z · σ (cid:12)(cid:12) z A ( z ) (cid:3) + τ rz πr ( z )cos α + f ext,z A ( z ) (52)Also, neglecting higher order effects, this equation can be further simplified. ∂ρ F u z A ( z ) ∂t + ∂∂z [ ρ F u z u z A ( z )] = ∂∂z [ σ zz A ( z )] + τ rz πr ( z )cos α + f ext,z A ( z ) (53) C. Conservation of energy
Finally, the transverse averaging technique is applied to the conservation of energy as described by Equation 28.There is exchange of energy, but no exchange of mass along the tubing, this results in the following equation. ∂ U A ( z ) /V∂t + ∂∂z (cid:2) u z U A ( z ) /V (cid:3) + p zz A ( z ) ∂ u z ∂z + ∂∂z (cid:2) ˆ e z · Φ Q A ( z ) (cid:3) + (cid:73) L ( z ) Φ Q · ˆ n cos α dl =12 u z G N A ( z ) (54)The flow of internal energy in the radial direction is zero. (cid:73) L ( z ) U A ( z ) /V (cid:12)(cid:12) z · ˆ n cos α dl = 0 (55)To extract the temperature profile, it is more convenient to replace the internal energy with the enthalpy density, H V = H /V . The enthalpy variation is a measure of thermal energy flow. H V = U /V + p + (cid:88) i µ i n i (56)0 FIG. 3: Wellbore structure for the one dimensional approximation. in which, n i = N i /V .Resulting in the following equation. ∂ H V A ( z ) ∂t + ∂ ( u z H V A ( z )) ∂z − ∂pA ( z ) ∂t − u z ∂pA ( z ) ∂z − (cid:88) i ∂ ( µ i n i A ( z )) ∂t − (cid:88) i ∂∂z [ µ i u z n i A ( z )] + ∂ (ˆ e z · Φ Q A ( z )) ∂z + (cid:73) L ( z ) Φ Q · ˆ n cos α dl =12 u z G N A ( z ) (57)Along the tubing, only thermal energy is exchanged, and it is assumed constant along the perimeter. Notice thatthe thermal energy transfer is positive flowing into the tubing. (cid:73) L ( z ) Φ Q · ˆ n cos α dl = − Φ Q,r cos α πr ( z ) (58) IV. DISCUSSION
One can apply the concepts introduced in the previous section to evaluate the one dimensional flow of oil. Thisapproximation can give an insight in model construction to extract the flow rate, pressure profile and temperatureprofile. Consider an isotropic homogeneous fluid with constant density in a production well. This is the simplifiedsingle phase flow.It is assumed that the oil well tubing displays a constant radius, r ( z ) = R . Thus, cos α = 1, and A ( z ) = πR . Ascan be seen in Figure 3, oil flows radially from the reservoir inside the drainage radius, r d , into the well. The externalforce per unit volume is partly due to the reservoir pressure, which forces the oil upwards, and suffers opposition fromgravity.The velocity profile, u z ( r ), in the tubing is not constant along the radius. Thus, applying the transverse averagingtechnique, one defines the average velocity over the cross-section area. u z = 1 πR (cid:90) π (cid:90) R u z ( r ) rdθdr = u (59)in which, u z ( r ) is the velocity profile along the radius.One approach is to combine the mass balance and momentum balance equations and solve the energy balanceequation separately with the help of thermodynamics to extract the temperature profile [4]. First, the problem issolved in steady-state. The state-state solution can be used as initial condition for the dynamic problem.1The estimated temperature profile is used to update temperature dependent properties, which are introduced inthe mass conservation and momentum conservation equations to calculate the pressure profile, production flow rate,or sound speed. Thus, combining the equations obtained in the previous section for the particular case of constanttubing radius, R , and constant chemical potential, µ i , one gets the following set. ∂ρ F ∂t + ∂ ( ρ F u ) ∂z = G N∂ρ F u ∂t + ∂ ( ρ F u u ) ∂z = ∂σ zz ∂z + τ rz R + ρ F ˆ e z · (cid:126)g sin θ ∂ H V ∂t + ∂ ( uH V ) ∂z = u G N + ∂p∂t + u ∂p∂z − ∂ (ˆ e z · Φ Q ) ∂z + Q,r R + (cid:80) i µ i (cid:16) ∂n i ∂t + ∂ ( u z n i ) ∂z (cid:17) (60)in which, θ is the tubing axis angle with respect to the horizontal and G N can be used to model the flow of oil fromthe reservoir into the well.According to Equation 22, the stress tensor components for a Newtonian fluid are as follows. σ zz = − p + 43 η ∂ u ∂z = − p (61) τ zr = η ∂ u ∂r = 0 (62)in which, u r = 0.In general, the multiphase flow in an oil well is not Newtonian, but to a first approximation the single phase flow oflight to medium crude oil can be considered Newtonian. Thus, the simplest approximation for oil/gas mixture flow isto model gas fraction as a compressible flow and oil fraction as a incompressible Newtonian fluid. To be more accurate,more sophisticated rheological models can be applied, such as: Ostwald-de Waele or power law fluid, Bingham fluidor plastic, Herschel-Bulkley or viscous-plastic or yield power law, among others.To account for pressure loss in the pipe due to friction, one can use the Darcy-Weisbach equation.2 τ zr R = − f D ρ F u R (63)in which, f D = 4 f F = 64 / Re for
Re < f F is the John T. Fanning frictionfactor.The third equation is used to calculate the temperature profile. The laminar flow assumption makes null anyconvective flow in the fluid layer. One should keep in mind the quantity for thermal energy flux, as presented inSection II is per unit area. One can assume that the thermal energy flows mainly by convection in the z-direction asthe fluid moves upwards. ∂ (ˆ e z · Φ Q ) ∂z = − h πR V dTdz (64)in which, h is the thermal convection coefficient ( J/ ( m .s.K )).In the radial direction, the flow of thermal energy is a combination of conduction and convection. The first layer isthe fluid inside the tubing. Next, one finds the tubing wall. Then, the annulus, the cement and the sand/rock layers,as shown in Figure 3. Thus, the radial direction can be modeled as a set of thermal layers, which can be modeled asa cascade of thermal RC-circuit. In steady-state, only the thermal resistance is considered [8].Φ Q A = ˙ Q r = U θ A ∆ T = R − θ ∆ T (65)in which, ˙ Q r is the rate of thermal energy transfer, U θ is the overall coefficient of thermal energy transfer, and R θ = UA is the thermal resistance. 1 U θ = (cid:88) i U θ,i (66)2in which, U θ,i is the coefficient of thermal energy transfer for the i -th well layer. • For linear conduction. R θ = ∆ zk A (67) • For radial conduction. R θ = ln r out /r in πk ∆ z (68) • For convection. R θ = 1 h A (69)For temperature transient, one needs to include the thermal capacity of each layer to model the energy transfer asa propagation through a cascade of thermal R θ C θ circuits. The obtained set of equations can be used in many flowconditions: no flow, steady-state flow, dynamic flow. V. APPLICATIONA. Steady-state: Ramey’s model
The most widely known model for wellbore thermal energy transmission for injection fluids is the Ramey’s model [1].This model assumes single phase steady-state vertical flow of a incompressible Newtonian fluid, materials propertiesare invariant with temperature, thermal energy transfer is radial only (ˆ e z · Φ Q = 0), and the energy flow is muchfaster than the fluid flow. The origin, z = 0, is set at the wellhead, with the z axis pointing downward. Under suchassumptions Equation 60 is simplified as follows. ∂ u ∂z = 00 = − dpdz + ρ F g d H V dz = dpdz + u (cid:104) Q,r R (cid:105) (70)The enthalpy per unit volume can be written as a function of temperature and pressure. d H V = ( C p /V ) dT + (1 − α p T ) dp (71)= ( C p /V ) dT − ( C p /V ) µ JT dp (72)For incompressible flow, α p = 0. Thus, µ JT = − ρ F c p = − VC P . d H V dz = ( C p /V ) dTdz + dpdz (73)Thus. d H V dz = ( C p /V ) dTdz = 1 u (cid:20) Q,r R (cid:21) (74)3Along the well, the equilibrium temperature at a distant position from the well is set by the geothermal profile.This is the surroundings temperature, T s . In the simplest case, the right side is constant, yielding a linear geothermalprofile. dT s dz = constant (75) T s ( z ) = az + b (76)in which, a is the geothermal gradient, and b is the surface temperature.As the tubing temperature, T t , is not known, it is best to write this temperature in terms of the surroundingstemperature, T s . To achieve this, one can subdivide the radial flow of thermal energy into two sections in series:tubing to casing and casing to surroundings.( R θ,tc + R θ,cs ) ˙ Q r = ( T t − T c ) + ( T c − T s ) (77)in which, ˙ Q = Φ Q × A . A R = u C p (cid:18) R θ,tc RA V + R θ,cs RA V (cid:19) (78) A R = m F u ∆ z C p m F (cid:18) πRU θ + R θ,cs A πR (cid:19) (79)in which, U θ is the overall coefficient of thermal energy transfer for all well layers from the liquid inside the tubing tothe casing.Defining the fluid density, ρ F = m F /V , the fluid injection rate, w = ( m F u / ∆ z ) = ρ F u A , and the specific heat atconstant pressure, c p = ( C p /m F ). A R = wc p (cid:18) πRU θ + R θ,cs A πR (cid:19) (80)According to Ramey’s model, a time-dependent dimensionless function f ( t ) is introduced to account for the thermalenergy flow from the casing into the surrounding formation. A R ( t ) = wc p (cid:18) πRU θ + f ( t )2 πk (cid:19) (81)in which, t is the injection time period.The thermal energy transfer is negative as it goes from the tubing to the surroundings. dT t dz = 1 u ( C p /V ) (cid:18) Q,r R (cid:19) = − wc P (cid:32) Q r R (cid:33) (82) dT t ( z, t ) dz = − T t ( z, t ) − T s ( z ) A R ( t ) = − T t A R ( t ) + ( az + b ) A R ( t ) (83)Solving this equation with the integrating factor method, one gets the Ramey’s model for the tubing temperatureprofile as a function of time for incompressible Newtonian liquid injection or a modified version for gas injection. • Newtonian liquid injection. T t ( z, t ) = az + b − A R ( t ) a + ( T − b + A R ( t ) a ) e − z/A R ( t ) (84)4 • Gas injection. T t ( z, t ) = az + b − A R ( t ) (cid:18) a + 1 c w c (cid:19) + (cid:20) T − b + A R ( t ) (cid:18) a + 1 c w c (cid:19)(cid:21) e − z/A R ( t ) (85)in which, T = T (0 , t ) is the constant injected fluid temperature at the surface, z = 0 and c w = 4 . J.g − o C − =778 . f t.lbf is the specific heat of water (mechanical equivalent of heat). B. Transient: simplified Hasan’s model
Firstly, consider a steady-state Newtonian oil flow [4, 5]. As in Ramey’s model the transfer of thermal energy isassumed radial only, ˆ e z · Φ Q = 0. The fluid density is not constant, thus Equation 60 is reduced as follows. ∂ ( ρ F u ) ∂z = 0( ρ F u ) ∂ u ∂z = − ∂p∂z + ρ F g sin θ ∂ ( uH V ) ∂z = u ∂p∂z + Q,r R (86)According to Equation 72, the enthalpy per unit volume can be written as a function of temperature and pressure.As the equation is presented in international units, SI, there is no need to introduce the g c conversion factor. Theflow of thermal energy towards the surroundings is accounted as negative. dTdz = (cid:18) µ JT + 1 ρ F c p (cid:19) ∂p∂z − wc p (cid:34) Q r R (cid:35) (87)To model the radial thermal energy flow, one can use Ramey’s model. dTdz = (cid:18) µ JT + 1 ρ F c p (cid:19) ∂p∂z − T t ( z, t ) − T s ( z ) A R ( t ) (88)Next, considering transient flow [4, 5] of a Newtonian fluid, Equation 60 is reduced as follows. ∂ρ F ∂t + ∂ ( ρ F u ) ∂z = 0 ∂ρ F u ∂t + ∂ ( ρ F u u ) ∂z = − ∂p∂z + ρ F g sin θ ∂ H V ∂t + ∂ ( uH V ) ∂z = ∂p∂t + u ∂p∂z + Q,r R (89)Combining the first two equations. ρ F ∂ u ∂t + ρ F u ∂ u ∂z = ρ F D u Dt = − ∂p∂z + ρ F g sin θ (90)Assuming constant mass flow rate, w = ρ F u A , and replacing the enthalpy per unit volume for the specific enthalpy, H m = H /m F . H V = m F V H m F = ρ F H m (91)5 uH V = wA H m F = wA H m (92) d H V = ρ F c p dT − ρ F c p µ JT dp (93) ∂ H V ∂t + wA ∂ H m ∂z = ρ F c p ∂T∂t − ρ F c p µ JT ∂p∂t + wc p A (cid:18) ∂T∂z − µ JT ∂p∂z (cid:19) (94)= ∂p∂t + u ∂p∂z + 2Φ Q,r R (95) ρ F c p ∂T∂t + wc p A ∂T∂z = µ JT ρ F c p (cid:18) ∂p∂t + wAρ F ∂p∂z (cid:19) + ∂p∂t + u ∂p∂z + 2Φ Q,r R (96) ∂T∂t + wAρ F ∂T∂z = µ JT (cid:18) ∂p∂t + wAρ F ∂p∂z (cid:19) +1 ρ F c P (cid:18) ∂p∂t + u ∂p∂z (cid:19) +1 ρ F c P (cid:18) Q,r R (cid:19) (97) DTDt = (cid:18) µ JT + 1 ρ F c P (cid:19) DpDt + 1 ρ F c P (cid:18) Q,r R (cid:19) (98)The final equation reduces to Equation 87 in steady-state. Similarly, one can use Ramey’s model for the radialthermal energy flow. DTDt = (cid:18) µ JT + 1 ρ F c P (cid:19) DpDt − T t ( z, t ) − T s ( z ) A R ( t ) (99) C. Two-phase flow
Two-phase flow is of great importance in many areas [18, 19]. As an example of two-phase flow, one can consider amixture of gas phase and liquid phase with liquid holdup, H L . The simplest model is to write the density as a linearcombination. ρ F = ρ L H L + ρ G (1 − H L ) (100)in which, H L is the liquid fraction also known as liquid holdup, and the density is a function of temperature.An equivalent expression can be written for the fluid velocity. ρ F (cid:126)u = ρ L (cid:126)u L H L + ρ G (cid:126)u G (1 − H L ) (101)6Should the gas and liquid velocity be the same, one has the pseudo-fluid approximation as found in the HomogeneousEquilibrium Model (HEM). Should the velocities be allowed to differ one has the Separated Flow Model (SFM) orDrift Flux Model (DFM), depending on the closure equations.To get the temperature profile, one can use the following approximation.1 A ∂ H V A∂t ≈ ρ F c p ∂T∂t − µ JT ρ F c p ∂p∂t (102)1 A ∂ ( uH V A ) ∂z ≈ ρ F c p u ∂T∂z − µ JT ρ F c p u ∂p∂z (103)Above the reservoir there is no flow from the surroundings into the tubing, thus G N = 0. From Equations 44, 53and 57, one can see that new terms are obtained should the area change with time or position. ∂ρ F ∂t = − ∂∂z ( ρ F u ) − ρ F D ln ADt ∂p∂z = − Dρ F u Dt − ρ F u ∂ u ∂z + τ rz r ( z ) cos α + ρ F ˆ e z · (cid:126)g sin θ − ρ F u D ln ADt − p ∂ ln A∂z
DTDt = (cid:16) µ JT + ρ F c P (cid:17) DpDt − Aρ F c P ∂ (ˆ e z · Φ Q A ) ∂z + ρ F c P Q,r r ( z ) cos α + pρ F c P D ln ADt (104)Examples of area models. • Constant area, A = A . • Exponential deposition, A = A e − tτ . • Power law, A = A (cid:16)(cid:80) Nn =0 (cid:0) zL (cid:1) n (cid:17) . • Combined, A = A (cid:16)(cid:80) Nn =0 (cid:0) zL (cid:1) n (cid:17) e − tτ .The model equations are combined with equations of state, empirical equations, and ancillary parameters [16]. Itcan also be applied under the assumptions of the Two Fluid Model (TFM) by writing one set of equations for eachphase and including equations for momentum and energy transfer between phases. VI. CONCLUSION
The Transverse Averaging Technique (TAT) is applied to the conservation equations for particles, momentum andenergy to provide the 1D+2D approximation in formal grounds. In this work, TAT is extended to tensor quantities.Within this approximation technique, a systematic model construction for oil and gas flow is presented, which canbe used to examine the various contributions and to calculate analytic expressions for a variety of problems. Theequations obtained for the flowrate and pressure profile are equivalent to the equations obtained by other methods,such as HEM, except that the generation term is maintained. The net generation term can be used to model theoil flow from the reservoir into the tubing. The energy equation is presented in terms of enthalpy, which is easier torelate to the temperature.The TAT technique can be applied to single phase flow or to two-phase flow. For single phase flow, two classicalmodels are obtained within this approximation, namely: Ramey’s steady-state model, and Hasan’s transient model.For the two-phase flow, it can be applied under the assumptions of the Homogeneous Equilibrium Model (HEM),7Separated Flow Model (SFM), Drift Flux Model (DFM), or Two Fluid Model (TFM). The proposed quasi-3D modelis useful to obtain analytic solutions which provide physical insight into the phenomena under scrutiny, including thevalidation of software tools. As can be seen in Equation 104, new terms are obtained should the area of the tubingvary analytically or randomly along the longitudinal direction, as a result of design, deposits, tubing roughness orother.
Acknowledgment
This work is dedicated to Prof. Anatoly A. Barybin of Saint-Petersburg State Electrotechnical University, inmemoriam .The author thanks CNPq for its continued support. He also thanks eng. Maur´ıcio Galasssi and eng. ManoelFeliciano da Silva Jr., both from PETROBRAS Brazilian Oil Company for bringing this problem to his attention.
Symbol nomenclature
A single bar over the quantity indicates transverse area averaging. A double bar over the quantity indicates it is atensor. α – angle between the local normal and the radial direction. α p – coefficient of thermal expansion. β T – compressibility at constant temperature.Φ Q – flux of thermal energy density.Φ q – total flux of kinetic energy density. ρ F – fluid density. ρ G – gas phase density. ρ L – liquid phase density. σ – Cauchy stress tensor. σ ii – normal stresses. τ ij – shear or viscous stresses. µ j – chemical potential for the j -phase. µ JT – Joule-Thomson coefficient. A – tubing cross-section area. C p – heat capacity at constant pressure. C θ – thermal capacity. c a – acoustic speed. c p – specific heat capacity at constant pressure. c w – specific heat of water. e K – total kinetic energy density. F B – probability density function or distribution function or particle density function. f D – Henry Darcy friction factor. f F – John T. Fanning friction factor.8 (cid:126)f ext – external force vector per unit volume. K b – bulk modulus. k – thermal conduction coefficient. G – Gibbs free energy. G N – net generation. H m – enthalpy per unit mass or specific enthalpy. H V – enthalpy per unit volume. H L – liquid fraction or holdup. h – thermal convection coefficient. m – parcel mass. P – total pressure tensor or total flux of momentum density. p – pressure tensor.˙ Q r – rate of thermal energy transfer. R θ – thermal resistance. r ( z ) – tubing cross-section radius as a function of the longitudinal direction. T – temperature. T c – casing temperature. T s – surroundings temperature. T t – tubing temperature. U – internal energy. U θ – coefficient of thermal resistance. (cid:126)u – average parcel velocity or fluid velocity. (cid:126)u G – average gas fraction velocity. (cid:126)u L – average liquid fraction velocity. V – volume. (cid:126)v – fluid parcel velocity. w – mass flow rate. [1] Henry J. Ramey Jr, Wellbore Heat Transmission , Journal of Petroleum Technology
14 (04) , 427-435 (1962).[2] A. R. Hasan, C. S. Kabir, X. Wang,
Development and Application of a Wellbore/Reservoir Simulator for Testing Oil Wells ,SPE Formation Evaluation
12 (03) , 182-188 (1997).[3] Octavio Cazarez-Candia, Mario A. Va´asquez-Cruz,
Prediction of pressure, temperature, and velocity distribution of two-phase flow in oil wells , J. of Petroleum Sci. and Engin., , 195-208 (2005).[4] Bulent Izgec, Transient fluid and heat flow modeling in coupled wellbore/reservoir systems , Texas A& M University, PhDThesis (2008).[5] A. R. Hasan, C. S. Kabir,
Wellbore heat-transfer modeling and applications , Journal of Petroleum Science and Engineering, , 127-136 (2012). [6] A. T. Singhe, J. R. Ursin, J. Henninges, G. Pusch, and L. Ganzer, Modeling of Temperature Effects in CO2 Injection Wells ,Energy Procedia, , 3927-3935 (2013).[7] Effects of tripping and swabbing in drilling and completion operations . Final report. Prepared for The Department of theInterior, Bureau of Safety and Environmental Enforcement Prepared by Bourgoyne Engineering LLC in collaboration withDarryl Andrew Bourgoyne, LLC and Bourgoyne Enterprises, Inc. July 17, 2017.[8] Dan Sui, Vebjørn Haraldstad Lang˚aker,
Downhole temperature Modeling for non-Newtonian fluids in ERD wells , Modeling,Identification and Control, , 131-149 (2018).[9] Fengrui Sun, Yuedong Yao, Mingqiang Chen, Xiangfang Li, Lin Zhao, Ye Meng, Zheng Sun, Tao Zhang, Dong Feng, Performance analysis of superheated steam injection for heavy oil recovery and modeling of wellbore heat efficiency , Energy , 795-804 (2017).[10] Fengrui Sun, Yuedong Yao, Guozhen Li, Xiangfang Li, Chengang Lu, Zhili Chen,
A model for predicting thermophysicalproperties of water at supercritical state in offshore CDTW , Measurement, , 241-251 (2018).[11] B. K. Sinha, M. S. Patel,
Recent developments in high precision quartz and langasite pressure sensors for high temperatureand high pressure applications , IEEE International Frequency Control Symposium (IFCS), pp. 1-13 (2016).[12] Leonardo B. M. Silva, Edval J. P. Santos,
FPGA-based smart sensor implementation with precise frequency to digitalconverter for flow measurement , 6th Southern Programmable Logic Conference, SPL 2010 - Proceedings 5483009, pp.21-26[13] Leonardo B. M. Silva, Edval J. P. Santos,
Modeling quality factor in AT-cut quartz-crystal resonators , Chip in Curitiba2013 - SBMicro 2013: 28th Symposium on Microelectronics Technology and Devices 6676121[14] Leonardo B. M. Silva, Edval J. P. Santos,
Quartz and GaPO4 pressure transducers for high resolution applications in hightemperature: A simulation approach , 2014 29th Symposium on Microelectronics Technology and Devices: Chip in Aracaju,SBMicro 2014 6940116[15] Leonardo B. M. Silva, Edval J. P. Santos,
Modeling high-resolution downhole pressure transducer to achieve semi-distributedmeasurement in oil and gas production wells , Journal of Integrated Circuits and Systems, , 2, 1-9 (2019).[16] Jean-No¨el Jaubert, Fabrice Mutelet, VLE predictions with the Peng-Robinson equation of state and temperature dependent k ij calculated through a group contribution method , Fluid Phase Equilibria Thermo-Fluid Dynamics of Two-Phase Flow , 1sted. New York, Springer Science (2006).[19] Donald A. Drew, Stephen L. Passman,
Theory of Multicomponent Fluids . Springer. Berlin (1999).[20] Richard L. Liboff,
Kinetic theory: classical, quantum and relativistic descriptions , Prentice-Hall (1990).[21] Anatoly A. Barybin, Edval J. P. Santos,
Transverse averaging technique for the depletion capacitance of nonuniformPN-junctions , Semicond. Sci. Technol. Large-signal and high-frequency analysis of nonuniformly doped or shaped pn-junction diodes , J. of Appl. Phys., , 114510 (2011).[23] Stephen Childress,
An introduction to theoretical fluid dynamics , New York University (2008).[24] M. H. Chaudhry,
Applied hydraulics transients , in Chapter 2
Transient-flow equations . 3rd Ed. Springer-Verlag New York(2014).[25] Reijo Kouhia,