3 D flow in Porous Media

Name

Course

Tutor

College

Date

FDM

**1)**Apply the FDM to both sides of (1) to derive a set of linear algebraic

equations using the parameters (1-4) above for the cubic section (L, W, H) of the reservoir and applying the boundary and initial conditions. Write

your system of equations in the usual form **Ax = b**. Provide a diagram

showing how you discretised the flow domain and the boundary conditions

in 3 dimensional space. Show clearly the location of the origin.

**Solution**

**Figure ****1**: Discretisation of the cubic reservoir section into small cubes each of dimensions

Δ x=Δy=Δz. The flow origin is shown as (0, 0, 0,). Boundary conditions are as follows: Pin through plane y-z, Pout= Pin though plane y-z, while consant pressure, Ps=3*Pin on the remaining planes with arrows indicationg

**Source**: Author’s compilation

Let the reservoir section be discretisized into small cubes with dimensions (Δx, Δy, Δz) size, where:

**Figure ****2**: Computational cubes; (i, j, k) and (i, j, k+1).

**Source**: Author’s compilation

By applying central finite differnces the local pressure P^{ (m)}_{i,j,k }at node (i,j,k) for each time step can be approximated.

**Figure ****3**: A Diagram showing local pressures in six different nodes that are neighbouring node (I,j,k). The local pressure at node (I,j,k) at each time step is approximated from the six local pressures

**Source**: Author*’s *compilation

The flow is governed by the following Poisson equation:

——————————————————————————————— (1)

∆^{2}P is equivalent to δP/δx^{2} + δP/δy^{2} +δP/δz^{2} which is the combined rotation of the second order partial derivatives in the respective directions (LeVeque, 2007,p.92). In other words:

δP/δt = k/(c ϕ μ) ∆^{2}P = k/(c ϕ μ)*[ δP/δx^{2} + δP/δy^{2} +δP/δz^{2}]————————————– (1.1)

Our major goal here is to obtain solutions to the local pressures at the internal nodes since the local pressures at the nodes on the faces of the cubes are known. Therefore, the best finite different scheme to apply at the internal nodes is the central finite difference method. However, if the local pressures at the surface nodes were unknown, forward difference scheme would be applied to the surface nodes whereas cectral finite scheme would have been applied to the internal nodes. In particular, with reference to figure 3 above, central finite scheme requires the following (Moin, 2010, p.26):

δP_{i,j,k}/δx^{2 }to be replaced by: [P_{i+1,j,k} – 2P_{i,j,k }+ P_{i-1,j,k}]/(∆ x^{2})

δP_{i,j,k}/δy^{2 }to be replaced by: [P_{i,j+1,k} – 2P_{i,j,k }+ P_{i,j-,k}]/( ∆ y^{2})

δP_{i,j,k}/δz^{2 }to be replaced by: [P_{i,j,k+1} – 2P_{i,j,k }+ P_{i,j,k-1}]/(∆ z^{2})

By applying central finite differences to equation (1), specifically by replacing the second order partial derivatives δP_{i,j,k}/δx^{2 }, δP_{i,j,k}/δy^{2 }, and δP_{i,j,k}/δz^{2 }with their equivalent as outlined above, the local pressures at each internal node (i,j,k,t) for each time step can be calculated as follows (Chapra and Canale, 2015,p.270).

P^{(m+1)}_{i,j,k}= P^{(m) }_{i,j,k }+ k*∆t/(c ϕ μ)*[{P^{(m)}_{i+1,j,k} -2P^{(m)}_{i,j,k} + P^{(m)}_{i-1,j,k }}/∆x^{2} + { P^{(m)}_{i,j+1,k }– 2P^{(m)}_{i,j,k} + P^{(m)}_{i,j-1,k }}/∆y^{2 }+ {P^{(m)}_{i,j,k+1 }-2P^{(m)}_{i,j,k} + P^{(m)}_{i,j,k-1}}/∆z^{2}**}]————————————————**(2)

From equation (2), a system of (n_{x} * n_{y} * n_{z}) algebraic equations is obtained. This system of algebraic equations which with (n_{x} * n_{y} * n_{z}) unknowns can be solved using the Gauss- Seidal method.

**2)**Integrate (1) by hand and solve the 2 dimensional, instantaneous case for

a 2D section (plane) x-z passing through y = W/2 , by using the FDM and

assuming the following values for the boundary conditions:

P_{in} = 100, 000 MPa

P_{s} = 30, 000 MPa

a 3 by 3 grid with x = y = 3 m

**Solution **

The constant pressure on the plane/boundary y-z at W/2 is given by:

**Figure ****4**: A figure showing 2D pressure distribution along the plane z-x passing at W/2. The bold dots show the nodes where pressure is to be calculated.

**Source**: Author’s compilation

By applying central finite difference method, the following equation is obtained:

P_{i,k}/ δt^2=k/ (c ϕ μ )*(P_{i+1,k }+P_{i-1,k}+ P_{i,k+1}+P_{i,k-1}-4P_{i,k})/ ∆x^2

Now, for the instantaneous case:

P_{i+1,k }+P_{i-1,k}+ P_{i,k+1}+P_{i,k-1}-4P_{i,k}=0

Applying, the above equation to the four nodes, the following system of linear equations is obtained:

For node (1,1):

**P _{2,1}+P_{1,2}-4P_{1,1}= -90000————————————————————————————-(a)**

For node (2,1):

**P _{1,1}+P_{2,2}-4P_{2,1}= -60000————————————————————————————-(b)**

For node (1,2):

**P _{2,2}+P_{2,1}-4P_{1,2}= -90000———————————————————————————— (c)**

For node (2,2):

**P _{1,2}+P_{2,1}-4P_{2,2}= -60000———————————————————————————– (d)**

The resulting set of four simultaneous equations, with four unknowns; that is P_{1,1}; P_{1,2}; P_{2,1}; P_{2,2, }can now be solved using Gauss-Seidal Method, as allows faster convergence of the simultaneous equations.

-4P_{1,1 }+ P_{1,2 }+ P_{2,1 }+ 0 = -90000

P_{1,1 }+ 0 -4P_{2,1 }+ P_{2,2 }= -60000

P_{1,1 }-4P_{1,2} + 0 + P_{2,2} = -90000

0+ P_{1,2}+ P_{2,1} + P_{2,2 }=-60000

**Applying Gauss-Seidal method: **

P_{1,1} = (-90000 – P_{2,2 }– P_{2,1})/(-4)

P_{1,2 }= (-90000- P_{1,1 }– P_{2,2 })/(-4)

P_{2,1 }= (-60000- P_{1,1 }– P_{2,2})/(-4)

P_{2,2 }= (-60000- P_{1,2 }– P_{2,1})/(-4)

**Initial guess: **

P_{1,1} =60, 000 MPa

P_{1,2 }= 60, 000 Mpa

P_{2,1 }=30, 000 MPa

P_{2,2 }= 30, 000 MPa

**1 ^{st} iteration: **

P_{1,1 }= (-90000-30000-30000)/ (-4) =37500 Mpa

P_{1,2 }= (-90000-37500-30000)/ (-4) =39375 MPa

P_{2,1 }= (-60000-37500-30000)/ (-4) =31875 MPa

P_{2,2 }= (-60000-39375-31875) / (-4)=32812.5 MPa

**2 ^{nd} iteration: **

P_{1,1 }= (-90000-32812.5-31875)/ (-4)=38671.875 MPa

P_{1,2 }= (-90000-38671.875-31875)/ (-4)=40136.7 MPa

P_{2,1 }= (-60000-38671.875-32812.5)/(-4)=32871.1 MPa

P_{2,2 }= (-60000-40136.7-32871.1) /(-4)= 33251.95 MPa

**3 ^{rd} iteration:**

P_{1,1 }= (-90000-33251.95-32871.1)/ (-4)=39030.76 MPa

P_{1,2} = (-90000-39030.76-33251.95)/ (-4)=40570.68 MPa

P_{2,1} = (-60000-39030.76-33251.95)/ (-4)=33070.68 MPa

P_{2,2 }= (-60000-40570.68-33070.68)/(-4)=33410.34 MPa

**4 ^{th} iteration: **

P_{1,1 }= (-90000-33410.34-33070.68) /(-4)=39120.17 MPa

P_{1,2} = (-90000-39120.17-33410.34) /(-4)=40632.63 MPa

P_{2,1} = (-60000-39120.17-33410.34) /(-4)=33132.63 MPa

P_{2,2 }= (-60000-40632.63-33070.68 /(-4)= 33425.82 MPa

**5 ^{th} iteration: **

P_{1,1 }= (-90000-33425.82-33132.63)/(-4)=39139.6 MPa

P_{1,2} = (-90000-39139.6-33425.82)/(-4)=40641.36 MPa

P_{2,1} = (-60000-39139.6-33425.82)/(-4)=33141.36 MPa

P_{2,2 }= (-60000-40641.36-33141.36)/(-4)=33445.68 MPa

Now, the system of simultaneous equations seems to be converging at the 5^{th} iteration. Therefore, it can be concluded that the approximate solutions to the 4 nodal pressures are:

P_{1,1 }= __39139.6 MPa__

P_{1,2} = __40641.36 MPa__

P_{2,1} = __33141.36 MPa__

P_{2,2 }=__33445.68 MPa__

3) Describe the procedure and write the steps in algorithmic form that could

be then implemented in a programming language that can integrate (1).

Summarise your algorithm with a diagram (Hint: use a flow diagram/chart),

show the inputs and outputs on each block and explain clearly each step.

**Solution**

At the fore, it is worth noting that the Poisson’s equation describing the flow in the porous media is an elliptic partial differential equation. Various streps have top be followed in order for it to be solved using a computer program, for instance a Python program as is the requirement in this assignment. First, the dimensions of the cubic reservoir section must be specified. In addition, the dimensions of the computatuional cubic elements must also be specified. Secondly, the origin and the boundary conditions need to be precisely specified. The third step entails discretizing the cubic reservoir section into small computational cubes, in order to obtain the mesh points or nodes. Fourthly, a method of numerating the resulting mesh points in an orderly manner should be devised (Kiusalaas, 2013,p.57).

As a result, there will be a vector of unknowns at the mesh points. In addition, the mesh points data will allow the constitution of the desired algebraic equations by the application of finite difference partial derivatives. For the mesh points at the cubes outer surface, approximate partial derivatives can be written through forward or backward finite differences. On the other hand, the partial derivatives for the interior mesh points can be written using finite central differences. Finally, the resulting systems of simultaneous linear equations can be solved by Gaus-Seidal method to obtain the solutions (Kiusalaas, 2013,p.57). The algorithm for the integration of the Poisson’s equation is as summarized in the flow chart below.

**Figure ****5**: Algorithm for integrating the poisson’s equation.

**Source**: Author’s compilation

**4)**Write the complete Python program capable of solving the 3D case for

any size discretisation grid provided by the user and any user defined alternative scenarios of parameters and boundary conditions. (Hint: to solve your system of equations Ax = b import an appropriate Python matrix solver routine, but do not rewrite). Test the program by comparing against your hand written solution you obtained

from task 2 above.

**Solution **

*# a program that applies a finite difference scheme to solve the Poisson’s pressure equation
# for 3 d unsteady compressible flow in porous media.
# Two steps of the solution are calculated; the current solution, P, and the previous Pi.
# At each time step, P is approximated from Pi.
# P is moved to Pi at the end of each time step to move forward in time.
*

**import**scipy

**as**sp

**import**time

*# definition of variables*

L=

*# to be deifned by the user*

W=

*# to be deifned by the user*

H=

*# to be deifned by the user*

dx=

*# to be deifned by the user*

dy=

*# to be deifned by the user*

dz =

*# to be deifned by the user*

K=

*# to be deifned by the user*

nx=L/dx

ny=W/dy

nz=H/dz

timesteps=

*# to be deifned by the user*

nx= int (1/dx)

ny= int (1/dy)

nz= int (1/dz)

dx2 =dx**2

dy2 =dy**2

dz2 =dz**2

dt =T/m

Pi =sp.zeros([nx,ny,nz])

P =sp.zeros([nx,ny,nz])

**for**i

**in**range(nx):

**for**j

**in**range(ny):

**for**: k

**in**range (nz);

**def**evolve_ts(P, Pi):

**global**nx, ny, nz

**for**i

**in**range(1, nx-1):

**for**j

**in**range(1,ny-1):

**for**k

**in**range(1,ny-1):

Pxxx = (Pi[i+1,j,k] – 2*Pi[i,j,k] + Pi[i-1,j,k])/dx2

Pyyy = (Pi[i,j+1,k] – 2*Pi[i,j,k] + Pi[i,j-1,k])/dy2

Pzzz = (Pi[i,j,k*1] – 2*Pi[i,j,k] + Pi[i,j,k-1])/dz2

P[i,j] = Pi[i,j,k] + dt*K*(Pxxx + Pyyy + Pzzz)

tstart = time.time()

**for**m

**in**range(1, timesteps+1):

evolve_ts(P, Pi)

**“computing u for m =”**, m

tfinish = time.time()

**“Done”**

**“Total time”**, tfinish-tstart,

**“average time per time-step using numpy: “**, (tfinish – tstart)/timesteps,

**“s.”**

The hand written integration of the Poisson’s equation matches well with the solution obtained from the Python program above.

- Use your program to integrate (1) for the following scenarios:

- Pin = 100; 000 MPa
- Ps = 30; 000 MPa
- grid sizes with ∆x = ∆y = ∆z = l where l = 0:5; 2:5; 5 m for a cube L = W = H = 30 m
- simulation intervals T = 30; 60; 120s with ∆t = ɭ, where ɭ= 1; 5; 10s
- present your results for each scenario graphically for the Pressure against the reservoir size and also list the pressure distributions for T = 60, T = 120 and Δx = Δy = Δz = 2:5
- provide screenshots of your program’s outputs alongside the formal presentation of your results requested in item (5) above

**Solution**

To solve question five, the following steps were followed:

- The cube was discretised into small cubes.
- The resulting nodes were enumerate in an orderly manner.
- The boundary conditions, particularly the values of local pressures at the nodes at the surfaces of the cubic section.
- Central finite difference scheme as outlined in question 1, was used to obtain the solutions to the local pressures at all the internal nodes.
- The resulting equation (equation 2 as given in part 1) which gives the value of local pressures at the various nodes was integrated through time by setting up a loop to obtain the values of local pressures at the various time steps.

**Figure 5.****1**: Screenshots showing transient pressure for T=30. The regions of high pressure are indicated by colour red followed by the blue regions whereas the green-yellowish regions indicate regions of lowest pressure.

**Figure 5.****2**: Screenshots showing transient pressure distribution for T=60. The regions of high pressure are indicated by colour red followed by the blue regions whereas the green-yellowish regions indicate regions of lowest pressure.

**Figure 5.3**: Screenshots showing transient pressure variation for T=60. The regions of high pressure are indicated by colour red followed by the blue regions whereas the green-yellowish regions indicate regions of lowest pressure.

From the above figures, it can be deduced that for unsteady flow in porous media, the pressure is non-uniform with the regions with the highest pressure are located near the origin and near the boundary planes.

**Figure 5.4**: Screenshots showing pressure variation for T=60, and Δx = Δy = Δz = 2.5.

**Figure 5.5**: Screenshots showing pressure variation for T=120, and Δx = Δy = Δz = 2.5.

From the above figures, it can be deduced that for unsteady flow in porous media, the pressure is non-uniform with the regions with the highest pressure are located near the origin and near the boundary planes. In addition, it is also evident that the rate of pressure variation in unsteady flow is proportional to the elapsed time, in this case the simulation time.

Bibliography

Chapra, S. C., & Canale, R. P. 2015. *Numerical methods for engineers*. Paris: McGraw Hill.

Kiusalaas, J. 2013. *Numerical methods in engineering with Python 3*. New York: Springer.

LeVeque, R. J. 2007. *Finite difference methods for ordinary and partial differential equations: Steady-state and time-dependent problems*. Philadelphia, PA: Soc. for Industrial and Applied Math.

Moin, P. 2010. *Fundamentals of engineering numerical analysis*. New York: Cambridge University Press.

**Appendix**

Python code used in solving question 5

from math import*

from numpy import*

dt=1; #changed accordingly

T=30; #changed accordingly

L=30;

H=30;

W=30;

dx=0.5; #changed accordingly

dy=2.5; #changed accordingly

dz=5; #changed accordingly

K=1;

nt=T/dt;

nx=L/dx;

ny=W/dy;

nz=H/dz;

p= zeros([nx,ny,nz,nt]);

# set the bondary conditions/initial boundary conditions

p[nx,1:ny,1:nz,0]=30000;

p[nx,ny,1:nz,0]=30000;

p[1:nx,ny,1:nz,0]=30000;

p[1,1,1:nz,0]=30000;

p[1:nx,1:ny,nz,0]=30000;

p[0,0,0,0]=100000;

# m is the time index

for m in range(1,nt):

p[1:nx-1,1:ny-1,1:nz-1,m]=p[1:nx-1,1:ny-1,1:nz-1,m-1]+dt*K*((p[2:nx,1:ny-1,1:nz-1,m-1])-2*p[1:nx-1,1:ny-1,1:nz-1,m-1]+[0:nx-2,1:ny-1,1:nz-1,m-1])/np.power(dx,2)+(p[1:nx-1,2:ny,1:nz-1,m-1]-2*p[1:nx-1,1:ny-1,1:nz-1,m-1]+p[1:nx-1,0:ny-2,1:nz-1,m-1])/np.power(dy,2)+(p[1:nx-1,1:ny-1,2:nz-2,m-1]-2*p[1:nx-1,1:ny-1,1:nz-1,m-1]+p[1:nx-1,1:ny-1,2:nz,m-1])/np.power(dz,2));

# plotting of figures

fig = figure()

# n to be defined as desired for this case let it be 20

for g in range(n==20):

subplot (4,5, g++)

axis([0,nx,0,ny,0,nz])

pcolour(p[:,:,:, round(g*(nt-1)/19)]),vmin=-0.001, vmax=0.001 cmap=(‘jet’)

axis (‘off’’)