This benchmark is available as a Jupyter notebook: TH2M/H2/mcWhorter/mcWhorter.ipynb.
McWhorter and Sunada propose an analytical solution to the twophase flow equation. A onedimensional problem was considered which describes the flow of two incompressible, immiscible fluids through a porous medium, where the wetting phase (water) displaces the nonwetting fluid (air or oil) in the horizontal direction (without the influence of gravity).
A detailed semianalytical solution and a convenient tool for calculating the solution for different material parameters can be found here.
Property  Symbol  Value  Unit 

Porosity  \(\phi\)  0.15  1 
Intrinsic permeability  \(K\)  \[1.0\cdot 10^{10}\]  \(m^2\) 
Residual saturation of the wetting phase  \[s_\mathrm{L}^{res}\]  0.02  1 
Residual saturation of the nonwetting phase  \[s_\mathrm{G}^{res}\]  0.001  1 
Dynamic viscosity of the wetting phase  \(\mu_\mathrm{L}\)  \[1.0\cdot 10^{3}\]  Pa s 
Dynamic viscosity of the nonwetting pha  \(\mu_\mathrm{G}\)  \[5.0\cdot 10^{3}\]  Pa s 
Brooks and Corey model parameter: entry pressure  \(p_b\)  5000  Pa 
Brooks and Corey model parameter: pore size distribution index  \(\lambda\)  3.0  1 
Property  Symbol  Value  Unit 

Initial saturation  \[s_\mathrm{L}(t=0)\]  0.05  1 
Injection boundary saturation  \[s_\mathrm{L}(x=0)\]  0.8  1 
The exact solution is not yet calculated in this notebook, instead the online tool by Radek Fučík is used. This tool calculates the solution and outputs the results with arbitrary accuracy as CSV files, which are plotted below.
import numpy as np
# Import analytical solution from a CSV file
exact = np.loadtxt('data/ref_solution_saturation.csv', delimiter=",")
# Zeroth column is location, first column is saturation
import os
out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')
if not os.path.exists(out_dir):
os.makedirs(out_dir)
from ogs6py.ogs import OGS
# run OGS
model=OGS(PROJECT_FILE="mcWhorter_h2.prj")
model.run_model(logfile=f"{out_dir}/out.txt", args=f"o {out_dir}")
OGS finished with project file mcWhorter_h2.prj.
Execution took 98.09925818443298 s
import vtuIO
# read OGS results from PVD file
pvdfile = vtuIO.PVDIO(f"{out_dir}/result_McWhorter_H2.pvd", dim=2, interpolation_backend="vtk")
import numpy as np
# The sampling routine requires a line of points in space
x_axis=[(i,0,0) for i in exact[:,0]]
# Only the last timestep is written (together with the initial condition).
# Thus the second element of the time vector (index = 1) is the to be sampled
time = pvdfile.timesteps[1]
# The numerical solution is sampled at the same supporting points as the analytical solution
sL_num = pvdfile.read_set_data(time, "saturation", pointsetarray=x_axis);
# Absolute and relative errors
err_abs = exact[:,1]  sL_num
err_rel = err_abs / exact[:,1]
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (14, 4)
fig1, (ax1, ax2) = plt.subplots(1, 2)
fig1.suptitle(r"Liquid saturation and errors at t="+str(time)+" seconds")
# Saturation vs. time
ax1.set_xlabel("$x$ / m", fontsize=12)
ax1.set_ylabel("$s_\mathrm{L}$", fontsize=12)
ax1.plot(exact[:,0], sL_num,'b', label=r"$s_\mathrm{L}$ numerical")
ax1.plot(exact[:,0], exact[:,1],'g', label=r"$s_\mathrm{L}$ exact")
lns2 = ax2.plot(exact[:,0], err_abs,'b', linewidth=2,
linestyle="", label=r"absolute error")
ax2.set_xlabel("$x$ / m", fontsize=12)
ax2.set_ylabel("$\epsilon_\mathrm{abs}$", fontsize=12)
ax3 = ax2.twinx()
lns3 = ax3.plot(exact[:,0], err_rel,'g', linewidth=2,
linestyle="", label=r"relative error")
ax3.set_ylabel("$\epsilon_\mathrm{rel}$", fontsize=12)
ax1.legend()
lns = lns2+lns3
labs = [l.get_label() for l in lns]
ax2.legend(lns, labs, loc=0)
fig1.savefig(f"{out_dir}/mcWhorter.pdf")
The numerical approximation fits the analytical solution very well in the whole area. Only in the area of the saturation front are deviations recognisable. These deviations are mainly due to the grid resolution and are well known in multiphase simulations. The error can be reduced almost arbitrarily by lowering the size of the grid elements.
This article was written by Norbert Grunwald. If you are missing something or you find an error please let us know.
Generated with Hugo 0.101.0
in CI job 262272

Last revision: October 19, 2022