FiPy 1D two-phase flow in porous mediaq, 11 October, 2019

Different approaches:

Coupled
Sequential
...



In [71]:

    
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
k = 1e-12 # m^2
phi = 0.4
u = 1.e-5
p0 = 100e5 # Pa
Lx = 100.
Ly = 10.
nx = 100
ny = 10
dx = Lx/nx
dy = Ly/ny

# 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 dkrw(sw):
    res = krw0*nw/(1-swc-sor)*((sw-swc)/(1-swc-sor))**(nw-1)
    return res


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

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

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

def dfw(sw):
    res = (dkrw(sw)/muw*kro(sw)/muo-krw(sw)/muw*dkro(sw)/muo)/(krw(sw)/muw+kro(sw)/muo)**2
    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 [72]:

    
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 [73]:

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

# create the cell variables and boundary conditions
sw = CellVariable(mesh=mesh, name="saturation", hasOld=True, value = swc)
p = CellVariable(mesh=mesh, name="pressure", hasOld=True, value = p0)
# sw.setValue(1,where = x<=dx)
sw.constrain(1.0,mesh.facesLeft)
#sw.constrain(0., mesh.facesRight)
sw.faceGrad.constrain([0], mesh.facesRight)
p.constrain(p0, mesh.facesRight)
p.constrain(1.05*p0, mesh.facesLeft)

Equations

$$\varphi \frac{\partial S_w}{\partial t}+u \frac{\partial f_w}{\partial x}=0$$

or $$\varphi \frac{\partial S_w}{\partial t}+\nabla.\left( u \frac{\partial f_w}{\partial S_w} S_w\right)+ \nabla. \left( u f_w-u\frac{\partial f_w}{\partial S_w} S_{w0} \right)=0$$



In [74]:

    
eq_p = DiffusionTerm(var=p, coeff=-k*(krw(sw.faceValue)/muw+kro(sw.faceValue)/muo))+ \
UpwindConvectionTerm(var=sw, coeff=k*(dkrw(sw.faceValue)/muw+dkro(sw.faceValue)/muo)*p.faceGrad)- \
(k*(dkrw(sw.faceValue)/muw+dkro(sw.faceValue)/muo)*sw.faceValue*p.faceGrad).divergence == 0

eq_sw = TransientTerm(coeff=phi, var=sw) + \
DiffusionTerm(var=p, coeff=-k*krw(sw.faceValue)/muw)+ \
UpwindConvectionTerm(var=sw, coeff=-k*dkrw(sw.faceValue)/muw*p.faceGrad)- \
(-k*dkrw(sw.faceValue)/muw*p.faceGrad*sw.faceValue).divergence == 0

eq = eq_p & eq_sw
steps = 200
dt0 = 500.
dt = dt0
t_end = steps*dt0
t = 0.0
viewer = Viewer(vars = sw, datamax=1.1, datamin=-0.1)
while t<t_end:
    swres = 1.0e6
    loop_count = 0
    while True:
        swres_new = eq.sweep(dt = dt)
        loop_count+=1
        if loop_count==1:
            sw_res = swres_new
        if swres_new>sw_res or loop_count>5:
            dt = dt/3
            continue
        swres=swres_new
        print(swres)
        if swres_new<1e-5:
            sw.updateOld()
            p.updateOld()
            t+=dt
            dt = dt0
            break
        
# Note: try to use the appleyard method; the overflow is a result of wrong rel-perm values        
viewer.plot()









    



0.0013423184598508084






    



/home/ali/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/tools/numerix.py:659: RuntimeWarning: overflow encountered in square
  return sqrt(add.reduce(arr**2))
/home/ali/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/tools/numerix.py:659: RuntimeWarning: overflow encountered in reduce
  return sqrt(add.reduce(arr**2))






    



---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-74-d407c58fe62a> in <module>
     19     loop_count = 0
     20     while True:
---> 21         swres_new = eq.sweep(dt = dt)
     22         loop_count+=1
     23         if loop_count==1:

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/term.py in sweep(self, var, solver, boundaryConditions, dt, underRelaxation, residualFn, cacheResidual, cacheError)
    230             self.residualVector = None
    231 
--> 232         solver._solve()
    233 
    234         return residual

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/solvers/scipy/scipySolver.py in _solve(self)
     24              raise Exception("SciPy solvers cannot be used with multiple processors")
     25 
---> 26          self.var[:] = numerix.reshape(self._solve_(self.matrix, self.var.ravel(), numerix.array(self.RHSvector)), self.var.shape)

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/solvers/scipy/linearLUSolver.py in _solve_(self, L, x, b)
     32                                             relax=1,
     33                                             panel_size=10,
---> 34                                             permc_spec=3)
     35 
     36         error0 = numerix.sqrt(numerix.sum((L * x - b)**2))

~/miniconda3/envs/myfipy/lib/python3.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py in splu(A, permc_spec, diag_pivot_thresh, relax, panel_size, options)
    309         _options.update(options)
    310     return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
--> 311                           ilu=False, options=_options)
    312 
    313 

RuntimeError: Factor is exactly singular

Analytical solution



In [ ]:

    
import fractional_flow as ff
xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R_oil = ff.frac_flow_wf(muw=muw, muo=muo, ut=u, phi=1.0, \
  k=1e-12, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0, \
  nw=nw, sw0=swc, sw_inj=1.0, L=Lx, pv_inj=5.0)



In [ ]:

    
plt.figure()
plt.plot(xt_prf, sw_prf)
plt.plot(x.value.squeeze()/(steps*dt), sw.value)
plt.show()



In [ ]:

    
?eq.sweep