An old FiPy solution to 1D BL

I'm not really sure if the result is correct. The results is indeed incorrect. See the updated code a few cells below that might work slightly better.



In [6]:

    
from fipy import *

u = 1.e-3
L = 100.
nx = 200
dt = 200.
dx = L/nx
muo = 0.002
muw = 0.001
mesh = Grid1D(dx = L/nx, nx = nx)
x = mesh.cellCenters
sw = CellVariable(mesh=mesh, name="saturation", hasOld=True, value = 0.)
sw.setValue(1,where = x<=dx)
sw.constrain(1.,mesh.facesLeft)
#sw.constrain(0., mesh.facesRight)
sw.faceGrad.constrain([0], mesh.facesRight)

eq = TransientTerm(coeff=1) + UpwindConvectionTerm(coeff = u
*(sw**2./muw)/(sw**2./muw+(1-sw)**2./muo) * [[1]]) == 0

sw.constrain(1.,mesh.facesLeft)
#sw.constrain(0., mesh.facesRight)
sw.faceGrad.constrain([0], mesh.facesRight)

steps = 100
viewer = Viewer(vars = sw, datamax=1.1, datamin=-0.1)
for step in range(steps):
    sw.updateOld()
    swres = 1.0e6
    while swres > 1e-5:
        swres = eq.sweep(dt = dt, var = sw)
        print(swres)
        
viewer.plot()

Updates solution on October 8th, 2019



In [4]:

    
from fipy import *

# relperm parameters
swc = 0.1
sor = 0.1
krw0 = 0.3
kro0 = 1.0
nw = 2.0
no = 2.0

# domain and boundaries
u = 1.e-3
L = 100.
nx = 50
dt = 200.
dx = L/nx

# fluid properties
muo = 0.002
muw = 0.001

# define the fractional flow functions
def krw(sw):
    res = krw0*((sw-swc)/(1-swc-sor))**nw
    return res

def kro(sw):
    res = kro0*((1-sw-sor)/(1-swc-sor))**no
    return res

def fw(sw):
    res = krw(sw)/muw/(krw(sw)/muw+kro(sw)/muo)
    return res

import matplotlib.pyplot as plt
import numpy as np

sw_plot = np.linspace(swc, 1-sor, 50)

Visualize the relative permeability and fractional flow curves



In [12]:

    
krw_plot = [krw(sw) for sw in sw_plot]
kro_plot = [kro(sw) for sw in sw_plot]
fw_plot = [fw(sw) for sw in sw_plot]
plt.figure(1)
plt.plot(sw_plot, krw_plot, sw_plot, kro_plot)
plt.show()

plt.figure(2)
plt.plot(sw_plot, fw_plot)
plt.show()



In [13]:

    
# create the grid
mesh = Grid1D(dx = L/nx, nx = nx)
x = mesh.cellCenters

# create the cell variables and boundary conditions
sw = CellVariable(mesh=mesh, name="saturation", hasOld=True, value = 0.)
# sw.setValue(1,where = x<=dx)
sw.constrain(1-sor,mesh.facesLeft)
#sw.constrain(0., mesh.facesRight)
sw.faceGrad.constrain([0], mesh.facesRight)

Some tests on fipy numerix



In [18]:

    
krw_cell = krw(sw)
krw_cell.value
# It works fine









    Out[18]:





array([0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875])



In [17]:









    Out[17]:





array([0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875, 0.0046875,
       0.0046875, 0.0046875])



In [ ]:

    
eq = TransientTerm(coeff=1) + UpwindConvectionTerm(coeff = u
*(sw**2./muw)/(sw**2./muw+(1-sw)**2./muo) * [[1]]) == 0

sw.constrain(1.,mesh.facesLeft)
#sw.constrain(0., mesh.facesRight)
sw.faceGrad.constrain([0], mesh.facesRight)

steps = 100
viewer = Viewer(vars = sw, datamax=1.1, datamin=-0.1)
for step in range(steps):
    sw.updateOld()
    swres = 1.0e6
    while swres > 1e-5:
        swres = eq.sweep(dt = dt, var = sw)
        print(swres)
        
viewer.plot()