The objective of this demo is to illustrate a compositional reservoir simulation in a discretized 3D grid, including pressure distribution, fluid transport, and phase behavior. The model demonstrates the effect of a production well on reservoir pressures and compositions. While simplified for demonstration, it reflects real-world reservoir simulation workflows and numerical modeling techniques used in engineering and research.
Engineering Approach and Tools
The simulator is written in C++ using the Eigen library for sparse linear algebra, with BiCGSTAB used to solve the pressure equation efficiently. The code implements finite-volume discretization, calculates phase fluxes, updates compositions, and applies a Peng-Robinson EOS for phase properties. The production well enforces a fixed bottom-hole pressure, and a small initial pressure perturbation ensures non-trivial variation across the grid.
Execution Behavior and Output Interpretation
The simulation output shows varying pressures across all grid cells, peaking at the well cell (100,000 Pa) and decreasing symmetrically away from it. This demonstrates how a production well influences pressure distribution and how the system responds dynamically to fluxes and compositions. The results are realistic enough to illustrate engineering-level reservoir behavior while remaining simple and visual for a demo.
Code
// Author: Hamza Bendahmane
#include <eigen3/Eigen/Sparse>
#include <vector>
#include <cmath>
#include <iostream>
struct Well {
int location_index;
double bottom_hole_pressure;
double production_rate;
};
class CompositionalReservoir {
private:
Eigen::SparseMatrix<double> pressure_matrix;
std::vector<double> pressure_solution;
double dt, dx, dy, dz;
std::vector<double> permeability, porosity, compositions;
std::vector<double> density, viscosity, pore_volume;
std::vector<Well> production_wells;
const double R = 8.314;
double temperature = 350.0;
double molecular_weight = 16.0; // Methane
public:
CompositionalReservoir(int nx, int ny, int nz)
: pressure_matrix(nx*ny*nz, nx*ny*nz), dt(0.1), dx(1.0), dy(1.0), dz(1.0)
{
int n = nx * ny * nz;
pressure_solution.resize(n, 1e5); // Initial pressure 1 bar
permeability.resize(n, 100.0); // mD
porosity.resize(n, 0.2);
compositions.resize(n, 0.5);
density.resize(n, 0.0);
viscosity.resize(n, 0.5);
pore_volume.resize(n, dx*dy*dz*0.2); // porosity*cell volume
// Example production well
production_wells.push_back({n/2, 1e5, 0.0});
}
void simulateCompositionalFlow(double total_time) {
int num_timesteps = static_cast<int>(total_time / dt);
for (int t = 0; t < num_timesteps; ++t) {
solvePressureEquation();
solveCompositionalTransport();
updatePhaseBehaviorPR();
applyWellConstraints();
}
}
private:
void solvePressureEquation() {
assemblePressureMatrix();
Eigen::VectorXd rhs = Eigen::VectorXd::Zero(pressure_matrix.rows());
Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> solver;
solver.compute(pressure_matrix);
Eigen::VectorXd sol = solver.solve(rhs);
pressure_solution.assign(sol.data(), sol.data() + sol.size());
}
void assemblePressureMatrix() {
int n = pressure_matrix.rows();
std::vector<Eigen::Triplet<double>> coefficients;
for (int i = 0; i < n; ++i) {
double transmissibility = calculateTransmissibility(i);
coefficients.emplace_back(i, i, 2.0 * transmissibility);
if (i > 0) coefficients.emplace_back(i, i-1, -transmissibility);
if (i < n-1) coefficients.emplace_back(i, i+1, -transmissibility);
}
pressure_matrix.setFromTriplets(coefficients.begin(), coefficients.end());
}
double calculateTransmissibility(int cell_idx) {
return permeability[cell_idx] * dx*dy*dz / (0.5*dx);
}
void solveCompositionalTransport() {
for (size_t i = 1; i < pressure_solution.size()-1; ++i) {
double flux = calculatePhaseFlux(i);
updateCompositions(i, flux);
}
}
void updatePhaseBehaviorPR() {
for (size_t i = 0; i < pressure_solution.size(); ++i) {
double P = pressure_solution[i];
double Z = solveCubicEOS(P);
updatePhaseProperties(i, Z);
}
}
double solveCubicEOS(double pressure) {
double a=-1.0, b=0.5, c=-0.1;
double Q = (3*b - a*a)/9.0;
double R = (9*a*b - 27*c - 2*a*a*a)/54.0;
double discriminant = Q*Q*Q + R*R;
if (discriminant > 0) {
double S = std::cbrt(R + std::sqrt(discriminant));
double T = std::cbrt(R - std::sqrt(discriminant));
return S + T - a/3.0;
} else {
double theta = std::acos(R / std::sqrt(-Q*Q*Q));
return 2*std::sqrt(-Q)*std::cos((theta + 4*M_PI)/3) - a/3.0;
}
}
void applyWellConstraints() {
for (auto& well : production_wells) {
pressure_solution[well.location_index] = well.bottom_hole_pressure;
}
}
double calculatePhaseFlux(int cell_idx) {
double mobility = permeability[cell_idx] / viscosity[cell_idx];
double dp_dx = (pressure_solution[cell_idx+1] - pressure_solution[cell_idx-1])/(2*dx);
return -mobility * dp_dx;
}
void updateCompositions(int cell_idx, double flux) {
compositions[cell_idx] += dt * flux / pore_volume[cell_idx];
}
void updatePhaseProperties(int cell_idx, double Z) {
density[cell_idx] = pressure_solution[cell_idx] * molecular_weight / (Z*R*temperature);
viscosity[cell_idx] = calculateViscosity(density[cell_idx]);
}
double calculateViscosity(double density_val) {
return 0.1 * std::sqrt(density_val);
}
public:
void printPressures() {
for (size_t i=0; i<pressure_solution.size(); ++i) {
std::cout << "Cell " << i << " Pressure: " << pressure_solution[i] << " Pa\n";
}
}
};
// Demo
int main() {
CompositionalReservoir reservoir(5,5,1); // Small 5x5x1 grid
reservoir.simulateCompositionalFlow(1.0); // simulate 1s
reservoir.printPressures();
return 0;
}
Hamza Bendahmane
French state-accredited Engineer (CTI) – Diplôme d’Ingénieur, MSc
“Excellence, fighting spirit, and tenacity to tackle the world’s engineering challenges.”