<jupyter_start><jupyter_text>**Unit 5 Finite-Difference Approximations to Linear-Flow Problems** Explicit Formulation (*Forward-Difference Approximation*)![Fig 5 29 Explicit Method](https://user-images.githubusercontent.com/51282928/76113357-94c97880-6016-11ea-8a3f-7658d9208cf9.PNG)$$p_{i}^{n+1} = p^i_n + (\frac{\alpha_c B_l \Delta t}{V_b \phi c_l})_i \cdot q_{lsc_i} + (\frac{\alpha_c B_l \Delta t}{V_b \phi c_l})_i \cdot [T_{lx_{i+1/2}}^{n} p_{i+1}^{n}-(T_{lx_{i+1/2}}^{n} + T_{lx_{i-1/2}}^{n})p_{i}^{n}+T_{lx_{i-1/2}}p_{i-1}^{n}]$$In the following codes, $(\frac{\alpha_c B_l \Delta t}{V_b \phi c_l})_i$ will be written as `factor`Transmissibility of coupling cells, $T_{lx_{i\pm1/2}}^{n}$, written as `T_min_half` or `T_plus_half`<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
!git clone https://github.com/yohanesnuwara/reservoir-engineering<jupyter_output>fatal: destination path 'reservoir-engineering' already exists and is not an empty directory.<jupyter_text>Example 5.8 Five reservoir gridblocks with no flow at both boundaries![Example 5 8](https://user-images.githubusercontent.com/51282928/75949840-09e36380-5eda-11ea-824c-dd5c9ac6f6d5.PNG)  For timestep 10 days<jupyter_code>"Task. determine pressure distribution during first year of production with timestep 10 days"

# known
p_initial = 6000 # initial pressure, in psia
delta_x = 1000 # ft
delta_y = 1000
delta_z = 75
ngrid = 5
grid_loc = 4 # grid location where production well is located
B = 1 # phase FVF, assumed constant over pressure, rb/stb
c = 3.5E-06 # phase compressibility, psi^-1
k_x = 15 # perm in x direction, md
poro = 0.18
mu = 10 # phase viscosity, cp

delta_t = 10 # days
qsc = -150 # minus means production, stb/d

# conversion
k_x = 15 * 1E-03 # darcy to mD

# calculate factor
Vb = delta_x * delta_y * delta_z
alpha = 5.615 # volume conversion factor, is a constant
factor = (alpha * B * delta_t) / (Vb * poro * c)
factor<jupyter_output><empty_output><jupyter_text>Because of uniform gridblocks, the equation for `tr_plus_coupling` and `tr_min_coupling` become:$$T_{lx_{i+1/2}}^{n} = T_{lx_{i-1/2}}^{n} = (\beta_c \frac{A_x k_x}{\mu_l B_l \Delta x})_{i+1/2}^{n} = (\beta_c \frac{A_x k_x}{\mu_l B_l \Delta x})_{i-1/2}^{n}$$<jupyter_code># calculate transmissibility of coupling cells
beta = 1.127 # transmissibility conversion factor, is a constant
Ax = delta_y * delta_z
T_plus_half = beta * ((Ax * k_x) / (mu * B * delta_x))
T_min_half = T_plus_half
T_min_half, T_plus_half<jupyter_output><empty_output><jupyter_text>There are 5 grids (`grid_1`, `grid_2`, `grid_3`, `grid_4`, `grid_5`), so the transmissibility coefficients will be: $$T_{lx_{i+1/2}}^{n}=(T_{lx_{1+1/2}}^{n}, T_{lx_{2+1/2}}^{n}, T_{lx_{3+1/2}}^{n}, T_{lx_{4+1/2}}^{n}, T_{lx_{5+1/2}}^{n})$$and$$T_{lx_{i-1/2}}^{n}=(T_{lx_{1/2}}^{n}, T_{lx_{1+1/2}}^{n}, T_{lx_{2+1/2}}^{n}, T_{lx_{3+1/2}}^{n}, T_{lx_{4+1/2}}^{n})$$The values of $T_{lx_{5+1/2}}^{n}=0$ and $T_{lx_{1/2}}^{n}=0$ (the edge of gridblocks), so the values of each are:$$T_{lx_{i+1/2}}^{n}=(0.127, 0.127, 0.127, 0.127, 0)$$and$$T_{lx_{i-1/2}}^{n}=(0, 0.127, 0.127, 0.127, 0.127)$$<jupyter_code>q = np.full(ngrid-1, T_min_half)
Ti_plus_halves = np.append(q, [0])
print(Ti_plus_halves)
p = np.full(ngrid-1, T_plus_half)
Ti_min_halves = np.append([0], p)
print(Ti_min_halves)
print("At grid 1, the coupling transmissibility coeffs are:", Ti_min_halves[0], "for T_min_half and:", Ti_plus_halves[0], "for T_plus_half.")
print("At grid 3, the coupling transmissibility coeffs are:", Ti_min_halves[2], "for T_min_half and:", Ti_plus_halves[2], "for T_plus_half.")
print("At grid 5, the coupling transmissibility coeffs are:", Ti_min_halves[4], "for T_min_half and:", Ti_plus_halves[4], "for T_plus_half.")<jupyter_output>[0.1267875 0.1267875 0.1267875 0.1267875 0.       ]
[0.        0.1267875 0.1267875 0.1267875 0.1267875]
At grid 1, the coupling transmissibility coeffs are: 0.0 for T_min_half and: 0.1267875 for T_plus_half.
At grid 3, the coupling transmissibility coeffs are: 0.1267875 for T_min_half and: 0.1267875 for T_plus_half.
At grid 5, the coupling transmissibility coeffs are: 0.1267875 for T_min_half and: 0.0 for T_plus_half.<jupyter_text>Requires array for $q_{sc}$<jupyter_code>qsc = [0, 0, 0, -150, 0] # production well in grid 4<jupyter_output><empty_output><jupyter_text>Calculate $p_{i}^{n+1}$ for each grid in each time<jupyter_code>pi = np.full(ngrid, p_initial) # array of pressure in each grid [6000, 6000, 6000, 6000, 6000]
time = np.arange(15, 370, delta_t)

pi_arr = []
min_arr = []
plus_arr = []

for j in range(len(time)):
  pnew_arr = []

  minus = pi[0]
  plus = pi[-1]
  minus_arr = []
  plusus_arr = []
  for i, obj in enumerate(pi):
    if i > 0:
      minus = pi[i-1]
    if i < (len(pi) - 1):
      plus = pi[i+1]

    pnew = pi[i] + (factor * ((Ti_plus_halves[i] * plus) - ((Ti_plus_halves[i] + Ti_min_halves[i]) * pi[i]) + (Ti_min_halves[i] * minus))) + (factor * qsc[i])
    pnew_arr.append(float(pnew))

    minus_arr.append(float(minus))
    plusus_arr.append(float(plus))

  pi = pnew_arr

  min_arr.append(minus_arr)
  plus_arr.append(plusus_arr)
  pi_arr.append(pi)

df = pd.DataFrame.from_records(pi_arr)
df = pd.DataFrame(pd.np.column_stack([time, df]), columns=['time', 'grid 1', 'grid 2', 'grid 3', 'grid 4', 'grid 5'])
df<jupyter_output><empty_output><jupyter_text>For timestep 15 days<jupyter_code>delta_t = 15 # days

# calculate factor
Vb = delta_x * delta_y * delta_z
alpha = 5.615 # volume conversion factor, is a constant
factor = (alpha * B * delta_t) / (Vb * poro * c)

pi = np.full(ngrid, 6000) # array of pressure in each grid [6000, 6000, 6000, 6000, 6000]
time = np.arange(15, 370, delta_t)

pi_arr = []
min_arr = []
plus_arr = []

for j in range(len(time)):
  pnew_arr = []

  minus = pi[0]
  plus = pi[-1]
  minus_arr = []
  plusus_arr = []
  for i, obj in enumerate(pi):
    if i > 0:
      minus = pi[i-1]
    if i < (len(pi) - 1):
      plus = pi[i+1]

    pnew = pi[i] + (factor * ((Ti_plus_halves[i] * plus) - ((Ti_plus_halves[i] + Ti_min_halves[i]) * pi[i]) + (Ti_min_halves[i] * minus))) + (factor * qsc[i])
    pnew_arr.append(float(pnew))

    minus_arr.append(float(minus))
    plusus_arr.append(float(plus))

  pi = pnew_arr

  min_arr.append(minus_arr)
  plus_arr.append(plusus_arr)
  pi_arr.append(pi)

df = pd.DataFrame.from_records(pi_arr)
df = pd.DataFrame(pd.np.column_stack([time, df]), columns=['time', 'grid 1', 'grid 2', 'grid 3', 'grid 4', 'grid 5'])
df<jupyter_output><empty_output><jupyter_text>Example 5.9 Similar to Ex 5.8, but **constant pressure at left boundary of reservoir**![Example 5 9](https://user-images.githubusercontent.com/51282928/76109912-081bbc00-6010-11ea-9f00-023e9d954ac6.PNG)The only modification is in the line:```pnew = pi[i] + (factor * ((Ti_plus_halves[i] * plus) - ((Ti_plus_halves[i] + Ti_min_halves[i]) * pi[i]) + (Ti_min_halves[i] * minus))) + (factor * qsc[i])pnew_arr.append(float(pnew))pnew_arr[0] = pi[0]```Where `pnew_arr[0] = pi[0]` means the pressure in the first gridblock is kept constant, equals to initial pressure `pi` 6000 psi.  For timestep 15 days<jupyter_code>pi = np.full(ngrid, p_initial) # array of pressure in each grid [6000, 6000, 6000, 6000, 6000]
time = np.arange(15, 370, delta_t)

pi_arr = []
min_arr = []
plus_arr = []

for j in range(len(time)):
  pnew_arr = []

  plus = pi[-1]
  minus_arr = []
  plusus_arr = []
  for i, obj in enumerate(pi):
    if i > 0:
      minus = pi[i-1]
    if i < (len(pi) - 1):
      plus = pi[i+1]

    pnew = pi[i] + (factor * ((Ti_plus_halves[i] * plus) - ((Ti_plus_halves[i] + Ti_min_halves[i]) * pi[i]) + (Ti_min_halves[i] * minus))) + (factor * qsc[i])
    pnew_arr.append(float(pnew))
    pnew_arr[0] = pi[0]

    minus_arr.append(float(minus))
    plusus_arr.append(float(plus))

  pi = pnew_arr

  min_arr.append(minus_arr)
  plus_arr.append(plusus_arr)
  pi_arr.append(pi)

df = pd.DataFrame.from_records(pi_arr)
df = pd.DataFrame(pd.np.column_stack([time, df]), columns=['time', 'grid 1', 'grid 2', 'grid 3', 'grid 4', 'grid 5'])
df<jupyter_output><empty_output><jupyter_text>Example 5.10 Halving the grid spacing of Ex 5.9 (more refined grids)![Example 5 10](https://user-images.githubusercontent.com/51282928/76109917-0c47d980-6010-11ea-95bf-87492c658368.PNG)  For timestep 15 days<jupyter_code># known
p_initial = 6000 # initial pressure, in psia
delta_x = 500 # ft, half the previous one
delta_y = 1000
delta_z = 75
ngrid = 10 # now twice the number of previous gridblocks
grid_loc = 4 # grid location where production well is located
B = 1 # phase FVF, assumed constant over pressure, rb/stb
c = 3.5E-06 # phase compressibility, psi^-1
k_x = 15 # perm in x direction, md
poro = 0.18
mu = 10 # phase viscosity, cp

delta_t = 15 # days
qsc = -150 # minus means production, stb/d

# conversion
k_x = 15 * 1E-03 # darcy to mD

# calculate factor
Vb = delta_x * delta_y * delta_z
alpha = 5.615 # volume conversion factor, is a constant
factor = (alpha * B * delta_t) / (Vb * poro * c)
factor
# calculate transmissibility of coupling cells
beta = 1.127 # transmissibility conversion factor, is a constant
Ax = delta_y * delta_z
T_plus_half = beta * ((Ax * k_x) / (mu * B * delta_x))
T_min_half = T_plus_half
T_min_half, T_plus_half
q = np.full(ngrid-1, T_min_half)
Ti_plus_halves = np.append(q, [0])
print(Ti_plus_halves)
p = np.full(ngrid-1, T_plus_half)
Ti_min_halves = np.append([0], p)
print(Ti_min_halves)
print("At grid 1, the coupling transmissibility coeffs are:", Ti_min_halves[0], "for T_min_half and:", Ti_plus_halves[0], "for T_plus_half.")
print("At grid 3, the coupling transmissibility coeffs are:", Ti_min_halves[2], "for T_min_half and:", Ti_plus_halves[2], "for T_plus_half.")
print("At grid 5, the coupling transmissibility coeffs are:", Ti_min_halves[4], "for T_min_half and:", Ti_plus_halves[4], "for T_plus_half.")
qsc = [0, 0, 0, 0, 0, 0, -75, -75, 0, 0] # production well in grid 7 and 8
pi = np.full(ngrid, p_initial) # array of pressure in each grid [6000, 6000, 6000, 6000, 6000]
time = np.arange(15, 370, delta_t)

pi_arr = []
min_arr = []
plus_arr = []

for j in range(len(time)):
  pnew_arr = []

  plus = pi[-1]
  minus_arr = []
  plusus_arr = []
  for i, obj in enumerate(pi):
    if i > 0:
      minus = pi[i-1]
    if i < (len(pi) - 1):
      plus = pi[i+1]

    pnew = pi[i] + (factor * ((Ti_plus_halves[i] * plus) - ((Ti_plus_halves[i] + Ti_min_halves[i]) * pi[i]) + (Ti_min_halves[i] * minus))) + (factor * qsc[i])
    pnew_arr.append(float(pnew))
    pnew_arr[0] = pi[0]

    minus_arr.append(float(minus))
    plusus_arr.append(float(plus))

  pi = pnew_arr

  min_arr.append(minus_arr)
  plus_arr.append(plusus_arr)
  pi_arr.append(pi)

df = pd.DataFrame.from_records(pi_arr)
df = pd.DataFrame(pd.np.column_stack([time, df]), columns=['time', 'grid 1', 'grid 2', 'grid 3', 'grid 4', 'grid 5', 'grid 6', 'grid 7', 'grid 8', 'grid 9', 'grid 10'])
df<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1643535427319211, 'fraction_numerical': 0.06464572680788896, 'mean_word_length': 4.675129533678756, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 39, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '47336', 'license': 'MIT', 'n_tokens_mistral': 4993, 'n_tokens_neox': 4445, 'n_words': 1082, 'path': 'Reservoir Simulation Ertekin/Unit 5 Finite-Difference Approximation to Linear-Flow Equations/fd1d_linearflow_explicit.ipynb', 'repo_name': 'datasolver/reservoir-engineering', 'size': '10952'}