Nodal analysis of flow¶

Marlim3 - Python Tutorials

In this notebook, we will set up a steady-state base case in Marlim3 and perform various simulations, varying some system properties, to illustrate the basic concepts of nodal analysis of flow.

In [1]:

%load_ext autoreload
%autoreload 2

import matplotlib.pyplot as plt
import numpy as np
import marlim3
import copy

%load_ext autoreload %autoreload 2 import matplotlib.pyplot as plt import numpy as np import marlim3 import copy

Introduction¶

One of the main activities of the flow assurance discipline is to predict the pressure drop along the fluid path, from the bottom of the well to the platform. In engineering fluid mechanics, this pressure drop is known as pressure loss. This phenomenon is directly related to flow characteristics in the pipeline, such as fluid properties, line geometry, among other factors, and can be modeled using multiphase flow simulators, such as Marlim3. Since the arrival pressure at the platform is previously defined to meet the demands of the primary processing of the surface plant, the variable of interest becomes the pressure at the bottom of the well. This pressure must be sufficient to overcome the pressure loss and ensure that the specified pressure is reached at the platform.

IPR: Inflow Performance Relationship¶

But what happens before the bottom of the well? The fluid comes from the oil reservoir, where the driving force for flow is the pressure differential between the reservoir pressure (known as static pressure, $P_e$) and the bottom hole pressure ($P_{wf}$). This phenomenon is modeled in the Reservoir discipline and depends on the characteristics of flow in the porous medium of the reservoir.

How, then, does the integration between the Reservoir and Flow Assurance disciplines occur for system modeling? Reservoir typically provides Flow Assurance with a curve that relates the produced flow rate to the static pressure and the bottom hole pressure. This curve, known as IPR (Inflow Performance Relationship), is used as a boundary condition in pipeline flow simulations. Also called the available pressure curve, the IPR relates the bottom hole pressure to the flow rate produced at surface conditions. It is a representation of reservoir productivity and is a function of the characteristics of the porous medium that composes it.

The simplest IPR model is the linear one:

$$Q = \texttt{IP}(P_e-P_{wf}),$$

where $Q$ is the flow rate at surface conditions, $P_e$ is the static pressure of the reservoir, and $P_{wf}$ is the bottom hole pressure.

The $\texttt{IP}$ parameter, known as the productivity index, represents the proportionality constant between the produced flow rate and the pressure differential between the reservoir and the bottom of the well. This differential is interpreted as the driving force for flow in the section in question. The linear model is particularly useful for describing single-phase flows in porous media. However, when pressure is below the saturation pressure and gas begins to be released already in the porous medium, it is recommended to adopt other models, such as the Vogel model, also available in Marlim3.

TPR: Tubing Performance Relationship¶

The IPR curve, provided by the RES discipline, describes the relationship between produced flow rate and pressure differential between the reservoir and the bottom of the well. Analogously, EE uses multiphase flow simulations to generate a curve called TPR (Tubing Performance Relationship), or required pressure curve, which relates the same produced flow rate to the pressure differential between the bottom of the well and the platform. The intersection point between the two curves represents the flow rate and pressure that the system will effectively produce when in equilibrium, in steady state.

This approach is known as nodal analysis, as it determines the state of the system at the intersection point of pressure curves at a specific node in the line — in this case, the bottom of the well.

In this notebook, we will explore how to perform nodal analysis using Marlim3.

Base case setup¶

We will populate a Branch object from Marlim3 for the base case:

In [2]:

base_case = marlim3.Branch()

base_case = marlim3.Branch()

Production fluid¶

We define a black oil type production fluid:

In [3]:

fluid = {
    "id": 0,
    "api": 20,
    "gor": 50,
    "gasDensity": 0.7,
    "bsw": 0.0,
}

base_case.productionFluid = [fluid]

fluid = { "id": 0, "api": 20, "gor": 50, "gasDensity": 0.7, "bsw": 0.0, } base_case.productionFluid = [fluid]

Marlim3 offers other fluid modeling possibilities: flash tables and compositional models.

Formation¶

In Marlim3, it is possible to add a list of thermal information related to rock lithology, which can be used to model heat exchanges between the well and the formation. Below, we include an example of a typical rock with one year of production time. This parameter is relevant because, as the well produces, the surrounding rock is gradually thermally influenced by the flowing fluid.

In [4]:

rock = {
    "id": 0,
    "conductivity": 2.5,    # W/(m.degC)
    "specificHeat": 1e3,    # J/(kg.degC)
    "density": 2.5e3        # kg/m3
}

base_case.initialConfig['formation'] = {
    'properties': [rock],
    'productionTime': 365   # days
}

rock = { "id": 0, "conductivity": 2.5, # W/(m.degC) "specificHeat": 1e3, # J/(kg.degC) "density": 2.5e3 # kg/m3 } base_case.initialConfig['formation'] = { 'properties': [rock], 'productionTime': 365 # days }

Pipeline¶

Materials¶

Below we define two materials to be used to compose the pipeline: steel and an insulating material.

In [5]:

steel = {
    "id": 0,
    "label": "Steel",
    "type": 0,              # solid
    "conductivity": 58,     # W/m.K
    "specificHeat": 480,    # J/kg.K
    "rho": 7850,            # kg/m3
}

completion_fluid = {
    "id": 1,
    "label": "Completion fluid",
    "type": 2               # water
}

cement = {
    "id": 2,
    "label": "Cement",
    "type": 0,              # solid
    "conductivity": 0.6,
    "specificHeat": 1000.0,
    "rho": 500.0
}

insulation = {
    "id": 3,
    "label": "Insulation",
    "type": 0,              # solid
    "conductivity": 0.12,   # W/m.K
    "specificHeat": 1214,   # J/kg.K
    "rho": 510,             # kg/m3
}

base_case.material = [steel, completion_fluid, cement, insulation]

steel = { "id": 0, "label": "Steel", "type": 0, # solid "conductivity": 58, # W/m.K "specificHeat": 480, # J/kg.K "rho": 7850, # kg/m3 } completion_fluid = { "id": 1, "label": "Completion fluid", "type": 2 # water } cement = { "id": 2, "label": "Cement", "type": 0, # solid "conductivity": 0.6, "specificHeat": 1000.0, "rho": 500.0 } insulation = { "id": 3, "label": "Insulation", "type": 0, # solid "conductivity": 0.12, # W/m.K "specificHeat": 1214, # J/kg.K "rho": 510, # kg/m3 } base_case.material = [steel, completion_fluid, cement, insulation]

Cross-sections¶

The production system is composed of several layers that will make up the different cross-sections.

In [6]:

inner_steel_layer = {
    "label": "Inner steel layer",
    "materialId": 0,
    "layerMeasurementType": "THICKNESS",
    "thickness": 0.0254 * 0.25,    # m
}

completion_fluid_layer = {
    "label": "Completion fluid layer",
    "materialId": 1,
    "layerMeasurementType": "THICKNESS",
    "thickness": 0.0254 * 3,       # m
}

casing_layer = {
    "label": "Casing steel layer",
    "materialId": 0,
    "layerMeasurementType": "THICKNESS",
    "thickness": 0.0254 * 0.5,     # m
}

cement_layer = {
    "label": "Cement layer",
    "materialId": 2,
    "layerMeasurementType": "THICKNESS",
    "thickness": 0.0254 * 3,       # m
}

insulation_layer = {
    "label": "Insulation layer",
    "materialId": 3,
    "layerMeasurementType": "THICKNESS",
    "thickness": 0.0254 * 2,       # m
}

inner_steel_layer = { "label": "Inner steel layer", "materialId": 0, "layerMeasurementType": "THICKNESS", "thickness": 0.0254 * 0.25, # m } completion_fluid_layer = { "label": "Completion fluid layer", "materialId": 1, "layerMeasurementType": "THICKNESS", "thickness": 0.0254 * 3, # m } casing_layer = { "label": "Casing steel layer", "materialId": 0, "layerMeasurementType": "THICKNESS", "thickness": 0.0254 * 0.5, # m } cement_layer = { "label": "Cement layer", "materialId": 2, "layerMeasurementType": "THICKNESS", "thickness": 0.0254 * 3, # m } insulation_layer = { "label": "Insulation layer", "materialId": 3, "layerMeasurementType": "THICKNESS", "thickness": 0.0254 * 2, # m }

Two cross-sections are defined in the system: one corresponding to the well, composed only of the steel layer, and another referring to the flowline and riser, formed by two layers — steel and insulating material. Both sections have the same internal diameter ($8$") and the same absolute roughness ($0.183$ mm).

In [7]:

well_section = {
    "id": 0,
    "label": "Well section",
    "innerDiameter": 8 * 0.0254,   # m
    "roughness": 0.183e-3,          # m
    "layers": [inner_steel_layer, completion_fluid_layer, casing_layer, cement_layer]
}

flowline_riser_section = {
    "id": 1,
    "label": "Flowline + Riser section",
    "innerDiameter": 8 * 0.0254,   # m
    "roughness": 0.183e-3,          # m
    "layers": [inner_steel_layer, insulation_layer]
}

well_section = { "id": 0, "label": "Well section", "innerDiameter": 8 * 0.0254, # m "roughness": 0.183e-3, # m "layers": [inner_steel_layer, completion_fluid_layer, casing_layer, cement_layer] } flowline_riser_section = { "id": 1, "label": "Flowline + Riser section", "innerDiameter": 8 * 0.0254, # m "roughness": 0.183e-3, # m "layers": [inner_steel_layer, insulation_layer] }

With the two sections defined, we can add them to the crossSection attribute:

In [8]:

base_case.crossSection = [well_section, flowline_riser_section]

base_case.crossSection = [well_section, flowline_riser_section]

Ducts¶

In Marlim3, the pipeline is divided into sections called ducts. In this case, there will be three portions: well, flowline, and riser. Each of them will have a length of $1000$ m and will be subdivided into 10 cells for discretization.

Well¶

A typical behavior of the geothermal gradient in wells shows a more pronounced temperature decrease in the final portion. To represent this behavior, we define temperatures at three points along the relative length of the duct: $90$°C at the beginning, $50$°C at 80% of the length, and $4$°C at the end, corresponding to the seabed temperature.

In [9]:

n_cells = 10

total_length = 1000  # m

discretization = {
    "numCells": n_cells,
    "length": total_length / n_cells  # m
}

formation_conditions = {
    "measuredPosition": [0, 0.8, 1],
    "ambientTemp": [90, 50, 4],  # degC
}

tubing = {
    "id": 0,
    "label": "Well",
    "crossSectionId": 0,
    "formationId": 0,
    "angle": np.pi / 2,              # rad
    "discretization": [discretization],
    "initialConditions": formation_conditions
}

n_cells = 10 total_length = 1000 # m discretization = { "numCells": n_cells, "length": total_length / n_cells # m } formation_conditions = { "measuredPosition": [0, 0.8, 1], "ambientTemp": [90, 50, 4], # degC } tubing = { "id": 0, "label": "Well", "crossSectionId": 0, "formationId": 0, "angle": np.pi / 2, # rad "discretization": [discretization], "initialConditions": formation_conditions }

Flowline¶

In the flowline, the external environment is specified as seawater. For this case, it is necessary to define the external temperature and external velocity:

In [10]:

n_cells = 10

total_length = 1000  # m

discretization = {
    "numCells": n_cells,
    "length": total_length / n_cells  # m
}

seabed_conditions = {
    "measuredPosition": [0, 1],
    "ambientTemp": [4, 4],    # degC
    "ambientVel": [0.1, 0.1], # m/s
}

flowline = {
    "id": 1,
    "label": "Flowline",
    "crossSectionId": 1,
    "environment": 1,        # seawater
    "angle": 0,                      # rad
    "discretization": [discretization],
    "initialConditions": seabed_conditions
}

n_cells = 10 total_length = 1000 # m discretization = { "numCells": n_cells, "length": total_length / n_cells # m } seabed_conditions = { "measuredPosition": [0, 1], "ambientTemp": [4, 4], # degC "ambientVel": [0.1, 0.1], # m/s } flowline = { "id": 1, "label": "Flowline", "crossSectionId": 1, "environment": 1, # seawater "angle": 0, # rad "discretization": [discretization], "initialConditions": seabed_conditions }

Riser¶

In the riser, the external environment is also seawater, but this time there is a temperature and velocity gradient from the seabed to the surface:

In [11]:

n_cells = 10

total_length = 1000  # m

discretization = {
    "numCells": n_cells,
    "length": total_length / n_cells  # m
}

riser_conditions = {
    "measuredPosition": [0, 1],
    "ambientTemp": [4, 20],   # degC
    "ambientVel": [0.1, 1],   # m/s
}

riser = {
    "id": 2,
    "label": "Riser",
    "crossSectionId": 1,
    "environment": 1,        # seawater
    "angle": np.pi / 2,              # rad
    "discretization": [discretization],
    "initialConditions": riser_conditions
}

n_cells = 10 total_length = 1000 # m discretization = { "numCells": n_cells, "length": total_length / n_cells # m } riser_conditions = { "measuredPosition": [0, 1], "ambientTemp": [4, 20], # degC "ambientVel": [0.1, 1], # m/s } riser = { "id": 2, "label": "Riser", "crossSectionId": 1, "environment": 1, # seawater "angle": np.pi / 2, # rad "discretization": [discretization], "initialConditions": riser_conditions }

With the three ducts defined, we can add them to the productionPipe attribute:

In [12]:

base_case.productionPipe = [tubing, flowline, riser]

base_case.productionPipe = [tubing, flowline, riser]

Visualizing the geometry:

In [13]:

base_case.plot_geometry();

base_case.plot_geometry();

Master 1¶

To indicate the position of the Subsea Christmas Tree (ANM), a masterValve object will be defined. This object represents the valve responsible for controlling production at the ANM. In steady state, masterValve has no effect beyond this indication of the ANM position, so the only property filled in will be measuredLength:

In [14]:

master_valve = {'measuredLength': 1000}  # m
base_case.masterValve = master_valve

master_valve = {'measuredLength': 1000} # m base_case.masterValve = master_valve

In a transient simulation, it would be important to define the valve opening, even if it remained constant over time.

Boundary conditions¶

The boundary conditions (BC) in our case are $1$) IPR upstream and $2$) pressure downstream.

BC1: Upstream IPR¶

BC $1$ must be defined through an ipr object placed right at the beginning of the pipeline (in the case below, $0.1$ m measured length).

In [ ]:

upstream_bc = {
    'id': 0,
    'iprType': 0,               # linear
    'prodFluidId': 0,        # id of the fluid defined earlier
    'measuredLength': 0.1,      # m
    'staticPressureTime': [0],  # s
    'staticPressure': [150],    # kgf/cm2
    'temperaturesTime': [0],    # s
    'temperatures': [90],       # degC
    'ipTime': [0],              # s
    'ip': [200],                # (sm3/d)/(kgf/cm2)
    'iiTime': [0],              # s
    'ii': [200],                # (sm3/d)/(kgf/cm2)
}

base_case.ipr = [upstream_bc]

upstream_bc = { 'id': 0, 'iprType': 0, # linear 'prodFluidId': 0, # id of the fluid defined earlier 'measuredLength': 0.1, # m 'staticPressureTime': [0], # s 'staticPressure': [150], # kgf/cm2 'temperaturesTime': [0], # s 'temperatures': [90], # degC 'ipTime': [0], # s 'ip': [200], # (sm3/d)/(kgf/cm2) 'iiTime': [0], # s 'ii': [200], # (sm3/d)/(kgf/cm2) } base_case.ipr = [upstream_bc]

The injectivity index, ii, is important to define for a possible case of reverse flow (in which the well "drinks").

BC2: Downstream pressure¶

BC $2$ must be defined through the separador field.

In [16]:

downstream_bc = {
    "time": [0],      # s
    "pressure": [10], # kgf/cm2
}

base_case.separator = downstream_bc

downstream_bc = { "time": [0], # s "pressure": [10], # kgf/cm2 } base_case.separator = downstream_bc

Simulation output specification¶

We can choose which output variables (in our case, profiles, or variations along the line) will be reported as simulation results.

Below, we specify pressure, temperature, liquid holdup, flow pattern arrangement, friction and hydrostatic components of pressure drop, and liquid flow rate.

In [17]:

output_vars_list = [
    "pressure",
    "temperature",
    "holdup",
    "flowPattern",
    "frictionPressureGradient",
    "hydrostaticPressureGradient",
    "stdLiqFlowRate",
]

output_vars = {"time": [0]} | {var: True for var in output_vars_list}

base_case.productionProfile = output_vars

output_vars_list = [ "pressure", "temperature", "holdup", "flowPattern", "frictionPressureGradient", "hydrostaticPressureGradient", "stdLiqFlowRate", ] output_vars = {"time": [0]} | {var: True for var in output_vars_list} base_case.productionProfile = output_vars

Running simulations and analyzing results¶

Base simulation¶

Let's run a simulation for the base case:

In [18]:

base_case.simulate()

base_case.simulate()

*******************************************************************************
                                  UFA!!!!!!!!                                  
             Paciencia, nove mulheres nao conseguem gerar uma crianca em um mes. 
                                      Tiao do Linkedin                                  
*******************************************************************************
ARQUIVO DE LOG: simulacao.log

In [19]:

base_case.resultados['productionProfile']

base_case.resultados['productionProfile']

Out[19]:

		Length (m) Boundary F	Length (m) Cell center C	Pressure (kgf/cm2) C	Temperature (C) C	Liquid holdup (-) C	Phase pattern indicator (-) F	Standard dead oil volumetric flow rate (Sm3/d) F	Hydrostatic term (Pa/m) F	Friction term (Pa/m) F	duct id	Elevation (m) F	Elevation (m) C
Tempo (s)
0	0	0	50	134.518559	90.000000	0.927144	1	0.000000	0.000000	0.000000	0	0	50
	1	100	150	126.296684	89.773514	0.931106	1	3096.360754	8000.206249	0.000000	0	100	150
	2	200	250	118.064001	89.520236	0.912150	1	3096.212880	8054.884164	83.910601	0	200	250
	3	300	350	109.959075	89.230407	0.891715	1	3096.204683	7934.112556	86.474855	0	300	350
	4	400	450	101.998195	88.903822	0.869592	1	3096.195357	7797.904344	89.476886	0	400	450
	5	500	550	94.199323	88.540359	0.845657	1	3096.185284	7644.220171	92.996174	0	500	550
	6	600	650	86.581366	88.140003	0.819873	1	3096.174891	7471.665760	97.141651	0	600	650
	7	700	750	79.164421	87.702843	0.792055	1	3096.163510	7279.531730	102.055834	0	700	750
	8	800	850	71.969834	87.229021	0.762140	1	3096.151745	7065.929041	107.900307	0	800	850
	9	900	950	65.019214	86.674747	0.730147	2	3096.140060	6829.990106	114.897091	0	900	950
	10	1000	1050	61.541182	85.996214	0.688905	2	3092.533158	6571.643014	123.336099	1	1000	1000
	11	1100	1150	61.412892	85.756724	0.677731	2	3096.286273	0.000000	126.695559	1	1000	1000
	12	1200	1250	61.285354	85.517850	0.677253	2	3096.287312	0.000000	125.024929	1	1000	1000
	13	1300	1350	61.157616	85.279610	0.676792	2	3096.287283	0.000000	125.222261	1	1000	1000
	14	1400	1450	61.029674	85.041983	0.676314	2	3096.287255	0.000000	125.421448	1	1000	1000
	15	1500	1550	60.901530	84.804966	0.675847	2	3096.287226	0.000000	125.619439	1	1000	1000
	16	1600	1650	60.773179	84.568560	0.675363	2	3096.287197	0.000000	125.824183	1	1000	1000
	17	1700	1750	60.644624	84.332763	0.674883	2	3096.287172	0.000000	126.022026	1	1000	1000
	18	1800	1850	60.515859	84.097574	0.674400	2	3096.287142	0.000000	126.228261	1	1000	1000
	19	1900	1950	60.386886	83.862993	0.673907	2	3096.287118	0.000000	126.430149	1	1000	1000
	20	2000	2050	57.095930	83.629017	0.684961	2	3099.852210	0.000000	126.637922	2	1000	1050
	21	2100	2150	50.760222	83.043170	0.655325	2	3096.147522	6200.138624	134.719073	2	1100	1150
	22	2200	2250	44.705306	82.451140	0.616756	2	3096.102541	5950.305983	149.090162	2	1200	1250
	23	2300	2350	38.979709	81.851804	0.576470	2	3096.093543	5619.660615	165.514754	2	1300	1350
	24	2400	2450	33.597121	81.245606	0.534880	2	3096.087238	5269.129342	186.325866	2	1400	1450
	25	2500	2550	28.566284	80.633042	0.492213	2	3096.082919	4902.484930	213.246289	2	1500	1550
	26	2600	2650	23.889573	80.014641	0.448920	2	3096.082073	4521.986266	248.927632	2	1600	1650
	27	2700	2750	19.559811	79.391027	0.405438	2	3096.085265	4131.980564	297.849180	2	1700	1750
	28	2800	2850	15.556534	78.762941	0.362029	2	3096.092960	3736.728159	368.296713	2	1800	1850
	29	2900	2950	11.836324	78.131230	0.318893	2	3096.108466	3339.010322	476.787730	2	1900	1950

In [20]:

base_case.plot_profiles(indicate_anm=True);

base_case.plot_profiles(indicate_anm=True);

In the profile graphs above, the dashed vertical line represents the position of the ANM. Observe the figure and analyze the behavior of each of the variables. Can you interpret the behavioral changes along the three sections (well, flowline, and riser), each with a length of 1000 meters?

Obtaining TPR curves¶

The solution obtained above corresponds to the equilibrium flow rate of the IPR, at which the available pressure at the bottom of the well equals the pressure necessary to overcome the pressure loss and reach the separator with the pressure specified in BC2.

To generate a TPR curve, it is necessary to remove the boundary condition associated with the IPR (deactivating it through the active boolean) and add a new boundary condition that corresponds to the specified liquid flow rate at the inlet:

In [21]:

base_case.ipr[0]['active'] = False

new_upstream_bc = {
    "id": 0,
    "prodFluidId": 0,          # id of the fluid defined earlier
    "measuredLength": 0.1,     # m
    "time": [0],               # s
    "temperature": [90]        # degC
}

base_case.liquidSource = [new_upstream_bc]

base_case.ipr[0]['active'] = False new_upstream_bc = { "id": 0, "prodFluidId": 0, # id of the fluid defined earlier "measuredLength": 0.1, # m "time": [0], # s "temperature": [90] # degC } base_case.liquidSource = [new_upstream_bc]

Note that we did not define liquidFlowRate in the new boundary condition (BC). It will be added in the for loop below, being configured thirty times, with different values at each iteration. This will result in thirty simulations, each representing a different flow condition. For each flow rate, an associated pressure at the bottom of the well ($P_{wf}$) will be calculated. The combinations of $Q$ (flow rate) and $P_{wf}$ obtained constitute the points that form the TPR curve.

In [22]:

%%time

flow_rates = np.linspace(10, 8000, 30)  # sm3/d
cases = []

# running a simulation for each flow rate
for flow_rate in flow_rates:
    cases.append(copy.deepcopy(base_case))
    cases[-1].initialConfig['classicOutput'] = False
    cases[-1].liquidSource[0]['liquidFlowRate'] = [flow_rate]
    cases[-1].simulate()

# extracting bottom hole pressures from the results
tpr_bottom_pressures = np.array([cases[i].resultados['productionProfile'].iloc[:, 2][0, 0] for i in range(len(cases))])

# re-enabling the IPR BC for the base case
base_case.ipr[0]['active'] = True
base_case.liquidSource[0]['active'] = False

%%time flow_rates = np.linspace(10, 8000, 30) # sm3/d cases = [] # running a simulation for each flow rate for flow_rate in flow_rates: cases.append(copy.deepcopy(base_case)) cases[-1].initialConfig['classicOutput'] = False cases[-1].liquidSource[0]['liquidFlowRate'] = [flow_rate] cases[-1].simulate() # extracting bottom hole pressures from the results tpr_bottom_pressures = np.array([cases[i].resultados['productionProfile'].iloc[:, 2][0, 0] for i in range(len(cases))]) # re-enabling the IPR BC for the base case base_case.ipr[0]['active'] = True base_case.liquidSource[0]['active'] = False

CPU times: total: 609 ms
Wall time: 1min 10s

Plotting IPR and TPR¶

Let's finally visualize the curves!

In [23]:

fig, ax = plt.subplots()

# anonymous function for IPR bottom hole pressure
Pwf_ipr = lambda Pe, Q, IP: Pe - Q / IP

# IPR properties
Pe = base_case.ipr[0]['staticPressure'][0]
IP = base_case.ipr[0]['ip'][0]

# calculating IPR
ipr_bottom_pressures = Pwf_ipr(Pe, flow_rates, IP)

# plotting IPR
ax.plot(flow_rates, ipr_bottom_pressures, label='IPR', marker='.', c='k')

# plotting TPR
ax.plot(flow_rates, tpr_bottom_pressures, label='TPR', marker='o')

ax.set_xlabel('Flow rate (sm3/d)')
ax.set_ylabel('Bottom hole pressure (kgf/cm2)')
ax.legend();

fig, ax = plt.subplots() # anonymous function for IPR bottom hole pressure Pwf_ipr = lambda Pe, Q, IP: Pe - Q / IP # IPR properties Pe = base_case.ipr[0]['staticPressure'][0] IP = base_case.ipr[0]['ip'][0] # calculating IPR ipr_bottom_pressures = Pwf_ipr(Pe, flow_rates, IP) # plotting IPR ax.plot(flow_rates, ipr_bottom_pressures, label='IPR', marker='.', c='k') # plotting TPR ax.plot(flow_rates, tpr_bottom_pressures, label='TPR', marker='o') ax.set_xlabel('Flow rate (sm3/d)') ax.set_ylabel('Bottom hole pressure (kgf/cm2)') ax.legend();

One of the intersection points between the curves represents the flow rate obtained in the base simulation, $Q = 3061$ sm$^3$/d, at which the IPR was used as a boundary condition (BC).

And what about the other intersection point? It corresponds to an unstable equilibrium point. This occurs when the TPR curve crosses the IPR from top to bottom. Naturally, simulators tend to converge to the stable solution, which acts as a numerical attractor of the system.

Finally, one last question about the figure: why is there a change in the direction of growth of the function in the TPR? Initially, bottom hole pressures decrease with increasing flow rate, but from a certain point onwards, they begin to increase. This behavior is related to the multiphase nature of the flow. At low flow rates, gas cannot efficiently entrain the liquid, which results in a significant increase in holdup and, consequently, in the hydrostatic term of pressure loss. As the flow rate increases, holdup decreases, reducing pressure loss until, at a certain point, the increase in pressure loss due to friction becomes predominant, causing pressure loss to increase again with increasing flow rate.

Increasing production¶

Nodal analysis enables the evaluation of different flow scenarios based on available degrees of freedom, allowing the choice of the most appropriate production flow rate for the system. This activity is important for the design of new production systems.

In the following code, we repeat the analysis performed earlier, now considering different gas-oil ratio (GOR) values of the fluid:

In [24]:

%%time

base_case.ipr[0]['active'] = False
base_case.liquidSource[0]['active'] = True

GOR = [25, 50, 100, 200, 300]  # sm3/sm3

flow_rates = np.linspace(10, 8000, 30)  # sm3/d
tpr_bottom_pressures = []

# running a simulation for each GOR and each flow rate
for gor in GOR:
    cases = []
    for flow_rate in flow_rates:
        cases.append(copy.deepcopy(base_case))
        cases[-1].initialConfig['classicOutput'] = False
        cases[-1].productionFluid[0]['gor'] = gor
        cases[-1].liquidSource[0]['liquidFlowRate'] = [flow_rate]
        cases[-1].simulate()
    tpr_bottom_pressures.append(np.array([cases[i].resultados['productionProfile'].iloc[:, 2][0, 0] for i in range(len(cases))]))

base_case.ipr[0]['active'] = True
base_case.liquidSource[0]['active'] = False

%%time base_case.ipr[0]['active'] = False base_case.liquidSource[0]['active'] = True GOR = [25, 50, 100, 200, 300] # sm3/sm3 flow_rates = np.linspace(10, 8000, 30) # sm3/d tpr_bottom_pressures = [] # running a simulation for each GOR and each flow rate for gor in GOR: cases = [] for flow_rate in flow_rates: cases.append(copy.deepcopy(base_case)) cases[-1].initialConfig['classicOutput'] = False cases[-1].productionFluid[0]['gor'] = gor cases[-1].liquidSource[0]['liquidFlowRate'] = [flow_rate] cases[-1].simulate() tpr_bottom_pressures.append(np.array([cases[i].resultados['productionProfile'].iloc[:, 2][0, 0] for i in range(len(cases))])) base_case.ipr[0]['active'] = True base_case.liquidSource[0]['active'] = False

CPU times: total: 2.95 s
Wall time: 5min 51s

In [25]:

fig, ax = plt.subplots()

markers = ['^', 's', 'D', 'p', 'o']

# plotting IPR
ax.plot(flow_rates, ipr_bottom_pressures, label='IPR', marker='.', c='k')

# plotting the various TPR curves
for i in range(len(GOR)):
    ax.plot(flow_rates, tpr_bottom_pressures[i], label=f'TPR - GOR {GOR[i]}', marker=markers[i])

ax.set_xlabel('Flow rate (sm3/d)')
ax.set_ylabel('Bottom hole pressure (kgf/cm2)')
ax.legend();

fig, ax = plt.subplots() markers = ['^', 's', 'D', 'p', 'o'] # plotting IPR ax.plot(flow_rates, ipr_bottom_pressures, label='IPR', marker='.', c='k') # plotting the various TPR curves for i in range(len(GOR)): ax.plot(flow_rates, tpr_bottom_pressures[i], label=f'TPR - GOR {GOR[i]}', marker=markers[i]) ax.set_xlabel('Flow rate (sm3/d)') ax.set_ylabel('Bottom hole pressure (kgf/cm2)') ax.legend();

It is observed that, in the case of GOR equal to 25, the TPR curve does not intersect the IPR curve, which indicates that the reservoir is not capable of providing, at the bottom of the well, the pressure necessary to sustain any flow rate under this condition. In this way, the well is said to be non-flowing and therefore artificial lift methods must be used. On the other hand, for GOR equal to 100, the stable crossing point occurs at a significantly higher flow rate compared to the base case (GOR 50). This behavior can be attributed to the reduction in flow resistance, resulting from the lower average specific mass of the fluid, which reduces the hydrostatic term of pressure loss – the main factor in this scenario.

However, there is a limit to increasing GOR. As GOR increases, the friction term of pressure loss becomes progressively more relevant. Therefore, the productivity gain in equilibrium flow rate between the cases of GOR 100 and GOR 200 is significantly smaller compared to the increment observed between the cases of GOR 50 and GOR 100. In the case of GOR 300, a reversal is observed: the equilibrium flow rate is lower than that recorded in the case of GOR 100.

Even without considering the IPR curve, it is possible to notice that, at high flow rates, the curves tend to cross, resulting in a trend reversal. In this region, higher GORs can result in higher bottom hole pressures (and, consequently, higher pressure losses) for the same flow rate. This behavior contrasts with that observed in the low flow rate region, where the hydrostatic term exerts greater influence and high GORs result in lower pressure loss.

Saving the base case¶

In the command below, we save the base case in a JSON to be able to import it in the next tutorial.

In [26]:

base_case.to_json('base_case1_tutorials')

base_case.to_json('base_case1_tutorials')

Try it yourself!¶

1. Plot new sets of TPR curves like those in the last figure we obtained, varying other parameters such as line diameter, oil density, or water fraction.

2. Choose some simulations that make up the TPR curves (you can access them in the casos list) and visualize the profiles along the flow. Pay attention to the differences from the base case, especially at high flow rates. Compare, for example, the flow patterns and the relative magnitudes of the friction and hydrostatic terms of pressure loss.

In [ ]: