Reservoir simulation code (1D-horizontal-single phase) FORTRAN

With this code you can practice with examples of reservoir simulation excersices with the follows considerations:

One dimension
Horizontal
One phase (oil)
Slightly compressible fluid
4 different types of boundary conditions
Implicit formulation
Peaceman well model
One well (specified production rate or specified flowing pressure)

You can change the number of cells, rock and fluid properties, boundaries conditions and well especification.

Just copy the code in a new file and try it

       

Program ARS_fortran_code
 implicit none
! variables
 real::dx,dy,dz,poro,Cr,k,Bo,vis,Cf,rw,s,qo,pwf,cell
 real::Gpe,Gpw,PE,PW,QE,QW,dt,time,beta,alpha,TT
 real::vol,Ax,TE,TW,Acum,cellP,FG,req,AB,AC,Po
 real, dimension(:), allocatable ::aa,bb,cc,dd,rr,Pt,Ptdt
 integer::BCE,BCW,i,WT,nx
 nx=5   !number of cells

 allocate(aa(nx))
 allocate(bb(nx))
 allocate(cc(nx))
 allocate(dd(nx))
 allocate(rr(nx))
 allocate(Pt(nx))
 allocate(Ptdt(nx))

 beta = 0.001127  !convertion factor
    alpha = 5.614  !convertion factor
 
 dx= 1000 !ft, dimension in x direction
 dy= 1000 !ft, dimension in y direction
 dz= 75 !ft, dimension in z direction
 
 poro=0.18 !porosity, fraction
 Cr=0 !rock Compressibility psi^-1
 Cf=0.0000035 !fluid Compressibility  psi^-1
 k=15 !permeability, md
 Bo=1 !Formation volume factor, RB/STB
 vis=10 !viscosity,cp 

 
 
 !Well information
 rw=3.5 !in, well radius
 s=0 !skin
 cellP=4 !well cell
 !WT = 1 Specified production rate
 !WT = 2 Specified flowing pressure
 !WT = 3 no well
 WT=1
 if(WT==1) then
 qo=-150 !well rate, STB/D
 end if

 if(WT==2) then
 pwf=4000 !Well Flowing pressure, psi
 end if

 if(WT==3) then
 qo=0 !Well rate, STB/D
 end if


 

 dt=1 !time step, days
 TT = 365 !total simulation time

 !Initial condition
 Po=6000 !psi


  

 Ax = dy * dz !Faces area
    vol = dx * dy * dz !Cells volume
    TE = ((beta * Ax * k) / (vis * Bo * dx))  !  East transmisibility
    TW = ((beta * Ax * k) / (vis * Bo * dx))  !  West transmisibility


 do i=1,nx
 Pt(i) = Po
 aa(i) = 0
 bb(i) = 0
 cc(i) = 0
 dd(i) = 0
 end do

 Acum = ((vol * poro * (Cf + Cr)) / (alpha * Bo * dt))
    time = dt

 !Boundary condition
 ! 1 = "No-flow"
 ! 2 = "Specified pressure gradient"
 ! 3 = "Specified flow rate"
 ! 4 = "Specified pressure boundary"

 BCE=1 !East boundary contidion
 BCW=1 !West boundary condition

  !Boundary condition
 !West
If (BCW == 2) Then
GPw=-0.2 !psi/ft
End If

If (BCW == 3) Then
QW=10 !STB/D
End If

If (BCW == 4) Then
PW=6000 !psi
End If

!East

If (BCE == 2) Then
GPe=-0.2 !psi/ft
End If

If (BCE == 3) Then
QE=10 !STB/D
End If

If (BCE == 4) Then
PE=6000 !psi
End If
  

 write(*,300)
 !Time loop
 DO WHILE ( time <= TT ) 

 do i=1,nx

 aa(i) = TW
 bb(i) = -(TW + TE + Acum)
 cc(i) = TE
 dd(i) = -Acum * Pt(i)
 end do
   !Boundary condition
 !West
If (BCW == 1) Then
bb(1) = -(TE + Acum)
End If

If (BCW == 2) Then
bb(1) = -(TE + Acum)
dd(1) = -GPw * TW * dx - Acum * Pt(1)
End If

If (BCW == 3) Then
bb(1) = -(TE + Acum)
dd(1) = -QW - Acum * Pt(1)
End If

If (BCW == 4) Then
bb(1) = -(TW * 2 + TE + Acum)
dd(1) = -TW * 2 * PW - Acum * Pt(1)
End If

!East

If (BCE == 1) Then
bb(nx) = -(TW + Acum)
End If

If (BCE == 2) Then
bb(nx) = -(TW + Acum)
dd(nx) = -TE * GPe * dx - Acum * Pt(nx)
End If

If (BCE == 3) Then
bb(nx) = -(TW + Acum)
dd(nx) = -QE - Acum * Pt(nx)
End If

If (BCE == 4) Then
bb(nx) = -(TW * 2 + TE + Acum)
dd(nx) = -TW * 2 * PE - Acum * Pt(nx)
End If

!Well

If (WT==1) Then
dd(cellP) = dd(cellP) - qo
AB = dd(cellP)
req = 0.14 * ((dx) **2 + (dy) ** 2) ** 0.5
FG = (2 * 3.1416 * beta * k * dz) / (Log(req / (rw / 12)) + s)
End If

If (WT==2) Then
AB = dd(cellP)
req = 0.14 * ((dx) ** 2 + (dy) ** 2) ** 0.5
FG = (2 * 3.1416 * beta * k * dz) / (Log(req / (rw / 12)) + s)
dd(cellP) = AB - (FG / (Bo * vis)) * pwf
bb(cellP) = bb(cellP) - (FG / (Bo * vis))
End If

call thomas(aa,bb,cc,dd,Ptdt,nx)

 do i = 1,nx

    Pt(i) = Ptdt(i)
    end do
If (WT==3) Then
qo = 0
pwf = 0
End If

If (WT==1)Then
pwf = qo / (FG / (Bo * vis)) + Pt(cellP)
End If

If (WT==2) Then
qo = -(FG / (Bo * vis)) * (Pt(cellP) - pwf)
End If
WRITE(*,400)Time,Ptdt(1),Ptdt(2),Ptdt(3),Ptdt(4),Ptdt(5),qo,pwf
time = time + dt

end do
300 FORMAT(5X,'time',5X,'Cell 1',5X,'Cell2',5X,'Cell 3',5X,'Cell 4',5X,'Cell 5',5X,'qo,'5X,'Pwf')
400 FORMAT(F8.0,2X,F8.3,2X,F8.3,2X,F8.3,2X,F8.3,2X,F8.3,2X,F8.3,2X,F8.3,2X,F8.3)

end program

!*****************************************************

SUBROUTINE thomas(a,B,c,d,x,n)
 ! Subroutine to solve a tridiagonal system
 ! a = subdiagonal vector
 ! B = diagonal vector
 ! c = superdiagonal vector
 ! d = right hand side vector
 ! x = solution vector
 ! n = number of diagonal vector elements 
 
 dimension a(*),B(*),c(*),d(*),x(*)
 integer n
 !Forward elimination
 do i=2,n
 B(i)=B(i)-a(i)*c(i-1)/B(i-1);
 d(i)=d(i)-a(i)*d(i-1)/B(i-1);
 end do
 !Back substitution
 x(n)=d(n)/B(n);
 do i=n-1,1,-1
 x(i)=(d(i)-c(i)*x(i+1))/B(i);
 end do
 end subroutine thomas

        

We can compare result with books examples (Ertekin et. al.) and understand better reservoir simulation.

if you have any question of how to use the code just ask us.

references

Ertekin, Turgay, Abou-Kassem, Jamal H. & King, Gregory R.(2001)Basic Applied Reservoir Simulation

Abou-Kassem, Jamal H., Farouq Ali, S. M. & Rafiq Islam, M. (2006) Petroleum Reservoir Simulation A basic approach