Introduction

Pressure buildup in subsea oil and gas equipment has remained as a challenging area of research for the past several decades. A major part of the early works were primarily focused on pressure buildup in simple cylindrical geometries such as pipes or tubes. One of the early works on this area of pressure buildup is from Sudhir (1984), which concentrates on relief valves for pressure release in trapped volumes. The paper by Sudhir also proposes an analytical method to determine the pressure buildup. Ramirez et al. (2010) performed studies on pressure rise in trapped volumes with rigid walls, and in their paper they emphasized that the heating of a trapped liquid is a well-recognized hazard scenario and the minor source of over-pressure can result in a catastrophic accident. The two primary strategies identified b Ramirez et al. (2010) for preventing over-pressure hazard are pressure containment and pressure relief.

Pressure rise in the Annulus of the casing string in subsea Christmas trees (XT) also known as “annular pressure buildup” abbreviated as APB has remained as an area of major concern for subsea oil and gas developments for the past several decades. The trapped pressure can be disastrous if proper mitigation techniques are not used. According to Hasan et al. (2010) managing APB for sustaining well production is extremely crucial in subsea wells (Hasan et al. 2010). Hasan et al. (2010) also proposed models for understanding and mitigating APB during continuous production. The intrinsic idea was to mimic the physical process with minimal input parameters to estimate pressure buildup in the annuli. Mathias et al. (2014) proposed analytical solutions to estimate the trapped pressure during carbon dioxide injection into gas reservoirs. Ai et al. (2015) presented a new model to calculate the critical production rate of high pressure and high temperature gas wells. The model of Ai et al. was applied to two gas wells in China.

Wang et al. (2018) in their work on APB tried to account for the variability of fluid properties. Wang et al. claim that their model improves the APB prediction. Bellary et al. (2013), emphasized in their paper that the annular pressure buildup (APB) remains an important consideration for many wells, not just subsea wells. The work of Bellarby et al. proposes a methodology for analyzing APB issues applicable to both subsea and surface oil and gas wells. They even presented a case to disprove a common belief that the fluid external to the sleeves equalizes with the reservoir over the time frame of the stimulation operation which prevents over-pressurization. Zhang and Wang (2016) proposed a methodology to determine the pressure buildup in single and multilayer annuli. Pattillo et al. (2006) in their work provide an analysis of a 16″ casing failure in the Gulf of Mexico (GoM). According to the Patillo et al., it was the first instance of an APB failure during drilling.

In the last few decades, a number of researchers have focused on the mitigation strategies for APB. Zhang et al. (2016) proposed completion techniques that could mitigate APB to a certain extent. They proposed several active and passive methods of APB mitigation. Udaya et al. (2017), in their studies, concluded that it is extremely difficult to quantify the reliability of a mitigation strategy or its corresponding interaction with other parts of the wellbore. Udaya et al. used a quantitative approach to design APB mitigation strategies. The use of computational methods to assess alternative scenarios and the corresponding mitigation strategies during the lifetime of the well has been discussed in detail. Dong and Chen (2017) investigated five different types of control measures for trapped annular pressure (TAP) or annular pressure buildup (APB), which contained increasing casing strengths, eliminating the trapped annular, releasing trapped pressure, balancing the expansion volume and blocking heat transfer. The advantages and disadvantages of a dozen control technologies were also analyzed and compared systematically, such as working principles, critical equipment, applications and key technical indicators.

Guan et al. (2016) recommended thermally insulated pipes as an alternative to mitigate annular pressure buildup. They also proposed a model to calculate APB with insulated pipes and a process to design insulated pipes to mitigate annular pressure. Zhou et al. (2015) conducted experiments and determined the influencing factors of APB. Based on experimental results, they proposed pressure management by injecting nitrogen in casing annulus. Experiments conducted by Zhou et al. show that isothermal compressibility factor and thermal expansion coefficient are the key influencing factors for annulus pressure buildup.

A detailed review of literature reveals that the majority of the works on annular pressure buildup are more focused on simplified analytical solutions. The application of a judicious combination of computational fluid dynamics (CFD) coupled with structural finite element analysis (FEA) to predict structural damage is unavailable till now. The present study is an attempt to propose a methodology for more realistic prediction of structural damage due to thermal expansion of trapped liquids in subsea equipment.

Even though there exist a vast collection of literature on APB and the corresponding mitigation strategies in subsea wells, a full-fledged analysis methodology is not shown explicitly in any publication. A combination of computational fluid dynamics (CFD) and structural finite element analysis (FEA) to predict pressure rise due to the thermal expansion of trapped liquids and the corresponding structural damage/failure is still not paid adequate attention. Moreover, the variation of thermophysical properties with respect to temperature and pressure has not been bestowed adequate attention in the vast majority of the available literature. The analytical equations deduced from thermodynamic equations are quite simplified and are incapable of taking the variation of thermophysical properties into account. A judicious combination of CFD and structural FEA can handle the variation of the thermophysical properties of the trapped liquids, account for geometric variation and can calculate the structural damage due to the pressure buildup.

Advantages and disadvantages/limitations of the present study

The present study splits the complex case of a fluid structure interaction (FSI) problem into a combination of CFD and static structural analyses. The current state of the art on this particular topic of pressure buildup is limited to simple analytical solutions, and the results obtained are mostly restricted to constant fluid properties, i.e., the variation of thermophysical properties is not taken into account. The present study can account for the variation of fluid properties as a function of both pressure and temperature and could provide deeper insights on the flow physics in the entrapped volume. The literature available on this topic is either focused on the thermal solution of pressure buildup or the structural problems due to pressure buildup. The present study bridges the gap between the two classes of problem widely investigated on pressure buildup, namely the thermal and structural analysis, and is more realistic and provides a conservative prediction of pressure buildup and the subsequent structural damage associated with it. The computational costs associated with such a segregated analysis is a lot less, as compared to a full-fledged FSI analysis.

The only limitation/disadvantage of the present study is that the thermal and structural analyses are not fully coupled, i.e., it does not perform a two-way FSI analysis to predict the pressure buildup. A complete FSI study is computationally expensive and is not suitable for projects with stiff deadlines.

The current paper is organized as follows: firstly, the problem is setup, as presented in section two. Section two presents the geometry considered for the analysis, mesh generation, material properties, and the initial and boundary conditions applied for the CFD analysis. Section three presents the results obtained from the steady-state and transient CFD analysis along with the structural analysis to predict structural damage.

Problem setup

CFD analysis work flow

The various steps in a CFD analysis can be summarized as preprocessing, solver/solution and post-processing stages, and the work flow is presented in Fig. 1. The CAD geometry received for analysis is cleaned and made analysis ready by simplification and removal of unwanted features such as fillets, chamfers, threads and slivers. The clean geometry is then meshed and supplied with initial and boundary conditions. The analysis is then set up with the respective models and run in a desktop computer or a super computer based on the number of elements considered suitable for the analysis. The variables of interest are regularly monitored to ensure convergence. The results obtained are then subject to post-processing to retrieve the required data for further analysis or conclusion of the project.

Fig. 1
figure 1

Work flow of the analysis

Geometry considered for the analysis

A subsea Christmas tree (XT) along with the hydraulic lines shown in Fig. 2 is considered for the present analysis. The XT consists of 5″ gate valves to control the main production flow path/bore and 2″ gate valves for regulating the annulus fluid. The 5″ production line has a production master valve (PMV), production swab valve (PSV), production wing valve (PWV) and a cross over valve (XOV). The 2″ annulus line consists of annulus master valve (AMV), annulus swab valve (ASV) and annulus wing valve. The XOV connects the 5″ production line to the 2″ annulus line. In case of a valve malfunction the fluid can be routed through the annulus line or vice versa. A major portion of the XT is insulated to prevent the heat dissipation to the seawater. The temperature of the production fluid has to be maintained above the hydrate formation temperature (HFT). Hydrates are formed at high pressure and relatively lower temperatures; this varies based on the composition of the production fluid. The hydrate inhibitors such as methanol and mono-ethylene glycol (MEG) are circulated during workover operations to inhibit the formation of hydrates. The circulation of hydrate inhibitors is quite expensive, and hence, the temperature in the tree must be maintained above the HFT. The XT’s and manifolds also have a number of chemical injection valves in critical regions to circulate the various inhibitors such as the paraffin/wax inhibitor, the hydrate inhibitor and the asphaltene inhibitor. This varies from field to field and is based on the various fluids produced, the flow assurance group along with the subject matter experts determine the suitability of the various inhibitors in a field.

Fig. 2
figure 2

Geometry of the Christmas tree and the hydraulic lines considered for the analysis

The prime intention of the present study is to capture the pressure buildup in small volumes. The subsea Christmas tree (XT) is a massive structure as compared to the thin control lines. The temperature distribution in the entire XT has to be calculated prior to the determination of the pressure buildup in the thin control lines. In order to determine the temperature distribution, one might be more inclined toward CFD than thermal FEA, as the heat transfer coefficients evolve as part of the CFD solution. In the case of thermal FEA the solution is as good as the heat transfer coefficients calculated using empirical correlations found in the literature and heat transfer textbooks. In the present study the temperature distribution through the entire XT is determined using a three-dimensional, steady-state, incompressible, and turbulent analysis. The seawater surrounding the structure is also modeled, and the natural convection is modeled using Boussinesq approximation. The use of Boussinesq approximation is justified as the temperature gradient between the structure and the seawater is not too high. The CFD analysis of the entire XT has been carried out with 40 million elements. The number of bodies in the XT is too many, and hence, generating a hexahedral mesh is quite impossible for the present study. One could even use a polyhedral mesh for such an analysis.

The seawater domain considered around the XT is thrice the size of the tree in the horizontal and transverse directions and twice the height of the XT in the vertical direction. The geometry of the XT along with the seawater considered for the steady-state CFD analysis is shown in Figs. 3 and 4. The domain must be sufficiently large enough to obtain proper numerical convergence. The boundary effects can be avoided by locating the boundaries far from the region of interest. The subsea current broadly referred to as “oceanic current” varies from location to location and depth of installation and hence is difficult to account for in such analysis. A conservative approach to model such problems is to consider the fluid as stagnant, as it would yield the maximum temperature in the XT. In brief, the flow around the XT is assumed to be predominantly due to free or natural convection. In some cases, it has been a common practice to assume an oceanic current of 0.3 m/s (approx 0.1 ft/s), as explained earlier the maximum temperature inside the XT will be mainly due to natural convection.

Fig. 3
figure 3

Top view—geometry of the Christmas tree surrounded by seawater

Fig. 4
figure 4

Isometric view of the XT inside the seawater domain

Mesh generation

The mesh generated for the XT is shown in Fig. 5. The mesh is mostly tetrahedral; in some primitive volumes, a combination of tetrahedral and hexahedral elements is used. A mesh sensitivity study has been carried at every stage of this analyses process to ensure the solution is truly grid independent. The mesh sensitivity study has not been included in this paper for the sake of brevity. One could even argue that a mesh independent study would be impossible for such complex geometries due to the thousands of parts present; it is still possible to perform analysis with at least three sets of meshes and concluding that the solution does not vary much. As a rule of thumb a variation of less than 1% is considered acceptable for such complex computational domains. A sectional view of the mesh used for the present study is presented in Fig. 6. It is advisable to attach an inflation layer outside the surface of the XT establishing contact with the seawater to capture the steep temperature gradients close to the wall. A detailed explanation of the turbulence model and the effect of mesh on the final solution is included in the** result section for the benefit of the readers.

Fig. 5
figure 5

Mesh generated for the XT

Fig. 6
figure 6

Sectional view of the mesh

As it can be seen from the sectional view of the mesh (Fig. 6), the mesh is densely packed closer to the main bore (5″ production line). The temperature gradients are too steep in the regions establishing close contact with the main bore. The only source of heat in the present analysis is the production fluid which is flowing at a very high pressure and temperature. It must be ensured that there exist at least 4 elements in thin slivers/sections. The mesh generated for the entire XT is more than 40 million elements, and hence, it is impossible to selectively refine the thin and long control lines without compromising on the quality of the results. In order to circumvent this problem, the results from the steady-state analysis are provided as initial and boundary conditions for the transient analysis performed to calculate the pressure rise/buildup.

Sixteen million elements have been packed in the thin and long control lines that can be split into the XT and tubing hanger control lines. These lines penetrate all the way to the downhole safety valve/surface-controlled subsea safety valve (DHSV/SCSSV). These hydraulic lines supply hydraulic fluid and electrical power for downhole operation. The mesh considered for the tree and that for tubing hanger lines are shown in Figs. 7 and 8. The mesh is quite fine, and inflation layers have been attached to ensure a grid-independent solution.

Fig. 7
figure 7

Mesh generated for the control lines

Fig. 8
figure 8

Mesh generated for the tubing hanger lines

Material properties

The subsea XT has thousands of metal parts, and it is extremely difficult to identify each and every part within the XT for this analysis. The prime objective of the present study is to propose a numerical methodology to predict the structural damage due to the pressure buildup. and hence. the intricate details about the various materials and their corresponding properties have been intentionally avoided. The XT comprises mainly of parts made out of steel and nickel-based alloys; one could conservatively use the following thermophysical properties to obtain a reasonable prediction of temperature distribution within the XT.

  • Density = 8000 kg/m3

  • Thermal conductivity = 42 W/m K

  • Specific heat capacity = 475 J/kg K

The thermophysical properties for some of the fluids used in the present study are listed in Table 1.

Table 1 Thermophysical properties of liquids

Initial conditions

Dirichlet (isothermal) boundary condition has been assumed for the entire production fluid volume for the steady-state analysis. The production fluid volume is shown in Fig. 9. As it can be seen, there exist slivers and thin volumes within the gate valves that makes it difficult to generate a hexahedral mesh and also to capture the velocity gradients. In this case, the fluid is trapped and hence there exists no flow. If the fluid is assumed to be at a non-uniform temperature, i.e., if there exists a temperature gradient within the production fluid volume, it could eventually lead to Rayleigh–Benard convection. The fluid domain should then be modeled with the Navier–Stokes equations. It is conservative to assume the fluid at the maximum operating temperature as that could give us the maximum temperatures at the control lines of interest. The production fluid is considered to be at 300 °F for all the cases in the present stidy.

Fig. 9
figure 9

Production fluid at isothermal condition at 300 °F

Boundary conditions

The bottom face of the seawater domain is considered as wall, and all the other faces are pressure outlet boundaries. The walls/surfaces of the XT establishing contact with the seawater are considered as coupled walls. The main assumptions of the analysis are listed below:

Assumptions for the analysis

  • The numerical investigations are performed for a steady-state temperature distribution.

  • The pressure increase calculated using thermodynamic calculations was using fixed properties and thereby there is a degree of uncertainty.

  • The present numerical investigations do not take into account the effect of thermal expansion of the metal bodies.

Results and discussion

Steady-state CFD analysis

A steady-state, three-dimensional, turbulent, and incompressible CFD analysis is performed on the entire XT with the production fluid at 300 °F. The seawater outside is assumed to be at a constant temperature of 40 °F. The convection on the seawater side is assumed to be due to natural/free convection, and Boussinesq approximation is used for modeling the buoyancy effects. The temperature contours at the mid-plane along the X and Y axes of the XT are presented in Figs. 10 and 11.

Fig. 10
figure 10

Temperature contours at the mid-plane along the X axis

Fig. 11
figure 11

Temperature contours at the mid-plane along the Y axis

As it can be seen in Figs. 10 and 11, the maximum temperature is at regions closer to the production fluid. The top portion of the XT can retain heat as it is insulated and the heat is not dissipated to the surroundings. The bottom portion of the XT has a couple of hydraulic connectors that are mostly uninsulated and hence dissipate heat rapidly to the seawater. The connectors are left uninsulated for operational and installation purposes. The flow lines from the XT are also insulated to prevent hydrate formation. The temperature scale has been adjusted to present a clear picture of the plume formation outside the XT. This can be seen in Figs. 12 and 13. The plume formation can be clearly seen at a temperature scale of 40–45 °F. The entire XT is at a temperature greater than 45 °F and hence is shown as red in Figs. 12 and 13.

Fig. 12
figure 12

Temperature contours at the mid-plane along the X axis (temperature scale has been adjusted)

Fig. 13
figure 13

Temperature contours at the mid-plane along the Y axis (temperature scale has been adjusted)

The velocity contours at the mid plane of the XT are shown in Fig. 14. The present analysis uses Boussinesq approximation and the validity of the model is restricted to a minimal change in density to the original density of the fluid which in turn happens at a very low temperature gradient.

Fig. 14
figure 14

Velocity contours at the mid-plane of the XT

The temperature contours in the production fluid and the control lines are shown in Fig. 15. A close-up view of the control lines is shown in Fig. 16. The temperature of the control fluid far from the production fluid is at a lower temperature due to the cooling effect from the seawater. The heat dissipates rapidly in metal regions establishing direct contact with the seawater.

Fig. 15
figure 15

Temperature contours in the production fluid and the control lines

Fig. 16
figure 16

A close-up view of the temperature contours in the control lines

The maximum temperature is extracted for the tree and the tubing hanger lines separately for the transient analysis. Figure 17 shows the maximum temperature at the tree line as 280.6 °F, and Fig. 18 presents the maximum temperature at the tubing hanger (TH) line. The maximum temperature at the TH line is 272.6 °F. These values correspond to the maximum values shown in the legend of the corresponding figures.

Fig. 17
figure 17

Maximum temperature at tree lines

Fig. 18
figure 18

Maximum temperature at the tubing hanger lines

Three cases have been considered for the present study. The maximum and minimum temperatures retrieved from the tree and tubing hanger lines are supplied as inputs (initial and boundary conditions) for the transient analysis to determine the pressure rise in the trapped volumes. The tree and tubing hanger lines are the trapped volumes in the present analyses. The maximum and minimum temperatures for the cases are given in Table 2.

Table 2 Maximum and minimum temperatures from the tree and TH lines

Thermodynamics calculations for pressure buildup

A simple equation from thermodynamics which can be derived from the Maxwell’s equations can serve as a good starting point for the pressure rise calculations. The derivation of the analytical pressure rise equation assuming the walls to be rigid is presented as “Appendix.”

The pressure rise calculated for the three case under investigation is presented in Table 3.

Table 3 Pressure rise calculated using the analytical expression

Transient analysis

The analytical expression deduced from Maxwell’s equations is not capable of taking the variation of thermophysical properties into account. CFD can serve as a valuable tool while accounting for the variation of properties. Such as an exercise does not exist in the literature, and a comparison of CFD results with standard analytical would serve as a good starting point. The trapped hydraulic fluid considered for the present study is Transaqua HT, and the properties of Transaqua HT are given below:

Properties of Transaqua HT

  • Density

    $$\rho = 1384.947978 + 4.11292 \times 10^{ - 7} P - 18.916878T^{0.5} \;{\text{in}}\,{\text{kg}}/{\text{m}}^{{3}}$$
  • Viscosity

    $$\mu = - 12.8482192 + 1.06027 \times 10^{ - 8} P + 36666.111/T^{1.5} \;{\text{in kg}}/{\text{m}}\,{\text{s}}$$
  • Thermal conductivity

    $$k = 0.42\;{\text{in W}}/{\text{m}}\,{\text{K}}$$
  • Specific heat capacity at constant pressure

    $$C_{p} = 3257\;{\text{in}}\,{\text{ J}}/{\text{kg}}\,{\text{K}}$$

The temperature contours extracted from the steady-state analysis are supplied as initial and boundary conditions for the transient analysis. The tree and the tubing hanger lines are simulated using a transient analysis with the shear stress transport k-ω model. The time step considered for the present analysis must be sufficiently fine enough to capture the pressure rise. An adaptive time stepping scheme has been used for the present analysis. The QUICK (quadratic upwinding with interpolation of convective kinetics) scheme has been used for the momentum and turbulence equations. The kω model utilized for the present analysis is y+ insensitive. The mesh used for the present analysis is too fine and can capture the steep gradients near the walls adequately. A mesh sensitivity analysis has been performed for the present study to ensure that the solution is truly mesh independent. The details of the mesh sensitivity study are not presented for the sake of brevity. The velocity contours within the tree and the TH lines are presented in Figs. 19 and 20. The heat pipe effect can be seen in the velocity contours, as the fluid parcels circulates within the trapped volume. The temperature contours within the tree and TH lines are presented in Figs. 21 and 22. The transient analysis for all the three cases have been run for 60 min, and the pressure buildup has been continuously monitored. The pressure rises to the maximum value within a duration of 5 min and thereafter remains constant at a fixed value.

Fig. 19
figure 19

Velocity contours within the tree lines

Fig. 20
figure 20

Velocity contours within the TH lines

Fig. 21
figure 21

Temperature contours within the tree lines

Fig. 22
figure 22

Temperature contours within the TH lines

The maximum pressure rise calculated from the transient CFD analysis and a comparison of the obtained results with the thermodynamic equation are presented in Table 4.

Table 4 A comparison of the results obtained from CFD analysis with thermodynamic calculations

Summary of results

Maximum temperature and the corresponding pressure increase observed in the trapped fluid for the various scenarios are as shown below:

  • Case 1 280.6 °F (411.2 K)—absolute tree line pressure 40,846.6 psi.

  • Case 2 280.6 °F (411.2 K)—absolute tree line pressure 29,200.4 psi.

  • Case 3 272.6 °F (406.8 K)—absolute TH line pressure 40,845.2 psi.

The maximum difference of approximately 6% between the thermodynamics and CFD calculations can be attributed to the variation of properties with respect to the local distribution of pressure and temperature within the tree and TH lines.

Structural analysis

The pressure calculated from the transient CFD analysis is used to determine the structural damage in the XT. The tree lines are within the composite valve block (CVB) which is a thick metal body and can withstand local failure to a certain extent, whereas the tubing hanger (TH) that hosts the TH lines is the most vulnerable for a structural damage. The extent of structural damage helps the Engineers decide whether a mitigation mechanism for the trapped volume is required or not. The TH body considered for the structural analysis is shown in Fig. 23. The holes in the TH body are for hydraulic and electric lines to the downhole equipment.

Fig. 23
figure 23

Geometry of the TH body

As the geometry is quite irregular, a tetrahedral mesh has been generated with midside nodes. The mesh generated for the TH body is shown in Fig. 24. The maximum pressure obtained from the transient analysis is imposed as a boundary condition for the structural analysis.

Fig. 24
figure 24

Mesh generated for the TH body

The pressure boundary condition applied inside the ports/holes is shown in Fig. 25. A maximum pressure of 40,847 psi has been applied for the present analysis. The top and bottom faces are fixed to obtain conservative numbers.

Fig. 25
figure 25

Pressure applied within the TH port/holes

The maximum directional deformation within the TH ports/holes containing the TH lines is less than 0.01″ along the X and the Y axes. The deformation is quite minimal, and hence, no mitigation mechanism/device to vent the pressure is required for the present operational conditions. The maximum deformation along the X axis is presented in Fig. 26. This localized failure or damage may not affect the integrity of the system in anyways. The proposed methodology can be extended to any geometry or operational conditions of interest.

Fig. 26
figure 26

Maximum directional deformation along X axis

Summary and conclusions

It is possible to numerically predict the maximum pressure rise due to the thermal expansion of trapped liquids in subsea oil and gas equipment and subsequently capture the structural damage caused by the trapped pressure. The availability of high-end desktop computers and super-computers have made the numerical analysis of problems with multiple degrees of freedom a reality. A class of problem that could immensely benefit due to the growth in CPU and GPU capabilities is the fluid structure interaction (FSI) problem in subsea oil and gas equipment. The size of the equipment used in subsea oil and gas fields is quite large and hence requires a lot of computational time and cost to analyze.

  • The present work has given a methodology to perform a one-way coupled CFD and FEA analysis to predict the structural damage in a subsea oil and gas equipment. The coupling of the CFD results to the structural FEA code seems to yield good results. Hence, the prediction of structural damage due to pressure buildup can be effectively simulated using the methodology discussed in this study.

  • The present study is a conservative approach to determine the structural damage and a more accurate way to calculate the damage will be to perform a full-fledged fluid structure interaction with complete coupling between the CFD and FEA codes. Nevertheless, the present study itself is sufficient from an application point of view.

  • An FSI study would be computationally expensive if performed for very fine time steps, and hence, it is advisable to use a loosely coupled one-way or a two-way interaction study for such problems.

  • The present analysis is confined to a single-phase liquid and the presence of gas bubbles or dissolved gases can potentially serve as a cushion, reducing the magnitude of the pressure buildup.

  • As it is shown from the present study, the structural damage can be restricted to just a thin metal layer close to the trapped volumes and hence may not affect the integrity of the entire system or equipment to any notable extent.

  • The proposed methodology can be extended to thin slivers and gaps in critical areas of interest, especially regions closer to the seals and gaskets, as this could potentially compromise the integrity of the system through continuous or intermittent leakage.

  • The numerical experiment performed in the present study could serve as an effective substitute for expensive laboratory testing.