Reservoir Fluid Flow Automation

Automating PVT Analysis with EOS

In modern reservoir engineering, relying on static PVT tables from lab reports is a bottleneck. We need dynamic, predictive models that tell us how fluids will behave under changing reservoir conditions. This is where automating Pressure-Volume-Temperature (PVT) analysis with Equations of State (EOS) comes in. By implementing these models in MATLAB, we can build powerful tools that predict fluid phase behaviour on the fly.

The goal is to create scripts that take raw fluid composition data and directly output the key fluid property tables needed for reservoir simulation. This process hinges on translating complex thermodynamic equations into efficient, reusable code.

Implementing EOS Models in MATLAB

The workhorses of the petroleum industry for phase behaviour are cubic Equations of State. We'll focus on two of the most common: the Soave-Redlich-Kwong (SRK) model and the Peng-Robinson (PR) model. While their thermodynamic underpinnings are deep, their implementation in code is quite straightforward. They provide a direct mathematical link between pressure, temperature, and volume for a given fluid composition.

P = \frac{RT}{V_m - b} - \frac{a(T)}{V_m(V_m + b) + b(V_m - b)}

To implement this in MATLAB, we create a function that takes temperature, pressure, and fluid component properties (like critical temperature, critical pressure, and acentric factor) as inputs. The function then calculates the parameters $a$ and $b$ for the mixture and solves the EOS to find the molar volume $V_m$ and compressibility factor $Z$ .

function [Z] = solve_peng_robinson(P, T, fluid_props)
    % fluid_props is a struct containing Tc, Pc, omega, etc. for each component

    % 1. Calculate mixture parameters 'a_mix' and 'b_mix'
    %    based on composition and binary interaction parameters.
    [a_mix, b_mix] = calculate_mixing_rules(T, fluid_props);

    % 2. Form the cubic equation for the compressibility factor Z
    R = 8.314; % Gas constant
    A = a_mix * P / (R^2 * T^2);
    B = b_mix * P / (R * T);
    
    % Z^3 - (1-B)Z^2 + (A-2B-3B^2)Z - (AB-B^2-B^3) = 0
    coeffs = [1, -(1 - B), (A - 2*B - 3*B^2), -(A*B - B^2 - B^3)];

    % 3. Solve the cubic equation for roots
    roots_Z = roots(coeffs);

    % 4. Select the correct physical root for Z
    %    For vapor: largest real root. For liquid: smallest real root.
    Z = select_physical_root(roots_Z);

end

Vectorised Calculations for Efficiency

A reservoir simulation might involve thousands of grid blocks, each with its own pressure and temperature. Calculating fluid properties one by one using a for loop would be incredibly slow. This is where MATLAB's strength in vectorisation becomes essential. By writing our functions to accept arrays of pressures and temperatures, we can perform calculations for the entire reservoir simultaneously.

Instead of passing single values for P and T, we pass vectors. All the intermediate calculations for parameters like A and B in our EOS function become array operations. MATLAB is highly optimised for this, leading to a massive performance increase. This allows us to rapidly generate property tables or update fluid states across a whole model in response to production changes.

% Example of vectorised PVT calculation

% Define a range of pressures for our analysis
pressure_vec = 1000:100:5000; % Pressure in psi

% Define a constant reservoir temperature
reservoir_T = 250; % Temperature in F

% Pre-allocate a vector for the results
Z_factor_vec = zeros(size(pressure_vec));

% Call the EOS function with the entire pressure vector
% The function must be written to handle vector inputs
Z_factor_vec = solve_peng_robinson_vectorised(pressure_vec, reservoir_T, fluid_props);

% Now we have Z-factors for the entire pressure range in a single variable
plot(pressure_vec, Z_factor_vec);
xlabel('Pressure (psi)');
ylabel('Compressibility Factor (Z)');
title('Z-Factor vs. Pressure at 250 F');

Automating Flash Calculations

The most fundamental task in phase behaviour prediction is the Its job is to determine, for a fluid of a known composition at a specific pressure and temperature, how much of it will be in the vapour phase and how much will be in the liquid phase. It also tells us the composition of each phase.

This is an iterative process. We make an initial guess for the split between vapour and liquid, calculate the properties of each potential phase using our EOS function, and check if the system is in equilibrium. If not, we adjust our guess and repeat. The Rachford-Rice equation is commonly used to solve the material balance at the heart of this loop.

\sum_{i=1}^{N} \frac{z_i (K_i - 1)}{1 + V(K_i - 1)} = 0

In MATLAB, we can implement this as a while loop that continues until the Rachford-Rice function result is close enough to zero. Inside the loop, we update our K-values based on the fugacities of the liquid and vapour phases, which we get from our EOS solver.

% Simplified Flash Calculation Loop

V = 0.5; % Initial guess for vapor fraction
K = initial_K_values_guess(P, T, fluid_props); % Wilson's correlation

iteration = 0;
max_iter = 100;

while iteration < max_iter
    % Solve Rachford-Rice equation for V
    % This often requires a numerical solver like fzero
    V = solve_rachford_rice(z, K);

    % Calculate liquid (x) and vapor (y) compositions
    x = z ./ (1 + V * (K - 1));
    y = K .* x;

    % Use EOS to find fugacities of liquid and vapor phases
    phi_L = calculate_fugacity(P, T, x, 'liquid');
    phi_V = calculate_fugacity(P, T, y, 'vapor');

    % Check for convergence (fugacity_liquid == fugacity_vapor)
    residual = sum((phi_L - phi_V).^2);
    if residual < 1e-10
        break; % Converged!
    end

    % Update K-values for next iteration
    K = K .* (phi_L ./ phi_V);
    
    iteration = iteration + 1;
end

Scripting Reservoir Fluid Properties

With our EOS and flash calculation tools in place, we can now compute the properties that reservoir engineers use every day. The Oil Formation Volume Factor ( $B_o$ ) and Gas Formation Volume Factor ( $B_g$ ) are two of the most critical. They relate the volume of a fluid at reservoir conditions to its volume at standard surface conditions.

$B_o$ is the ratio of oil volume at reservoir conditions to the volume of the same oil once it reaches the stock tank. $B_g$ is the ratio of gas volume at reservoir conditions to its volume at standard conditions.

To calculate these, we use our automated workflow. For a given reservoir pressure, we first perform a flash calculation to find the amounts and compositions of the oil and gas phases. Then, we take the resulting oil and flash it again down to standard (surface) conditions. The ratio of the initial oil volume to the final stock tank oil volume gives us $B_o$ . A similar process yields $B_g$ .

function [Bo, Bg] = calculate_fvf(P_res, T_res, fluid_props)
    % Standard conditions
    P_std = 14.7; % psi
    T_std = 60;   % F

    % 1. Flash at reservoir conditions
    [V_res, x_res, y_res] = flash_calc(P_res, T_res, fluid_props.z, fluid_props);
    L_res = 1 - V_res; % Moles of liquid (oil)

    % 2. Get molar volume of oil at reservoir conditions using EOS
    Z_oil_res = solve_eos(P_res, T_res, x_res, 'liquid');
    vol_oil_res = Z_oil_res * R * (T_res + 460) / P_res;
    
    % We need total volume of the liquid part
    total_vol_oil_res = L_res * vol_oil_res;

    % 3. Flash the reservoir liquid (x_res) to standard conditions
    [V_std, x_std, ~] = flash_calc(P_std, T_std, x_res, fluid_props);
    L_std = 1 - V_std; % Moles of stock tank oil

    % 4. Get molar volume of stock tank oil
    Z_oil_std = solve_eos(P_std, T_std, x_std, 'liquid');
    vol_oil_std = Z_oil_std * R * (T_std + 460) / P_std;
    
    % We need total volume of this final liquid
    total_vol_sto = L_std * vol_oil_std;
    
    % 5. Calculate Bo (rb/stb)
    % We must account for the fact that we started with L_res moles of oil
    % and ended up with L_std moles from that initial amount.
    Bo = total_vol_oil_res / (L_res * vol_oil_std);
    
    % Bg calculation would follow a similar logic for the gas phase.
    Bg = ...; % Left as an exercise

end

By wrapping these functions into a single script, we build a robust workflow. We can now take any fluid sample, input its composition, and programmatically generate all the necessary PVT tables for professional reservoir simulation software, bridging the gap between raw data and actionable reservoir insights.

Quiz Questions 1/5

What is the primary purpose of a flash calculation in PVT analysis?

Quiz Questions 2/5

Why is vectorisation a crucial technique when implementing Equations of State (EOS) in MATLAB for large-scale reservoir simulation?

Advanced MATLAB for Subsurface Engineering Applications