In this benchmark problem we describe the influence of a heat flux on a twophase system (aqueous phase gas phase), also commonly known as the heat pipe effect. The heat pipe effect describes heat transport in a porous medium caused by convection due to capillary forces. Detailed description of the heatpipe prblem can be found in (Helmig et al., 1997).
This benchmark case considers the heat pipe effect for a twophase system in a horizontal, onedimensional column. As illustrated in the figure below, a constant heat flux $q$ is applied at one end of the column, which is sufficient to heat the water at this end of the column above boiling temperature. As the water vaporizes, the vapour moves towards the other end of the column, condensing to water again as it moves through cooler regions and, hence, giving up the latent heat of vaporization. The resulting nonlinear water saturation profile along the horizontal column causes a capillary pressure gradient to arise, leading to a flux of the liquid phase pointing back towards the heat source: the mass of the vapour moving away from the heat source is equal to the mass of the condensate moving back towards the heat source.
When a heat pipe evolves from a singlephase liquid region as a consequence of constant heat injection, the vapour that is created at the heat source, displaces the water first. As a result, the latent heat of evaporation is given up by the vapour at the condensation front, and a mass flux is formed by the condensate pointing away from the heat source towards the completely cooled end of the column (here liquid exists in a singlephase region). The capillary forces which arise due to the nonlinear saturation profile along the column are strong enough to transport the entire condensate back to the heat source, as long as the heat source is strong enough to ensure constant vaporization and liquid saturation reduction at the heat source. Like this, a closed loop of mass cycle is created within the heat pipe, and the heat pipe reaches its maximum length when all the water is evaporized at the heat source and the evaporation front is far enough from the heat source so that the temperature is just enough for the vaporization. At this point the heat pipe will not propagate further and remains stationary.
While heat conduction is essential in the singlephase regions in front and behind the heat pipe for the heat transfer, the temperature gradient between the two ends of the heat pipe is relatively small, and therefore, within the heat pipe itself, heat conduction has no primary importance compared to convection.
This benchmark case assumes that the gas phase contains an additional component  air. The presence of air in the system results in the following considerations that are taken into account:
For verification purposes, the numerical solution of this benchmark problem is compared to a semianalytical solution. The semianalytical solution is briefly introduced in the following sections.
The analytical solution of the heat pipe problem solves a coupled system of first order differential equations for pressure, temperature, saturation and mole fraction derived by (Udell and Fitch, 1985). The here presented solution is slightly modified by (Helmig et al., 1997), as not all original assumptions proposed by (Udell and Fitch, 1985) are satisfied in this benchmark, thus some generalization are applied and will be mentioned alongside the original considerations:
\begin{align} \mu_g=x_g^a\mu_g^a+x_g^w\mu_g^w \end{align}
\begin{align} \rho_g=\frac{p_g}{(RT)(x_g^aM^a+x_g^wM^w)} \end{align}
The following material parameters and heat flow value are considered as input parameters.
import math
import numpy as np
q = 100.0 # heat injection [W/m²]
K = 1e12 # permeability [m²]
phi = 0.4 # porosity []
p_ref = 101325 # reference pressure [Pa]
T_ref = 373.15 # reference temperature [K]
lambda_G = 0.2 # thermal conductivity of gas phase [W/mK]
lambda_L = 0.5 # thermal conductivity of liquid phase [w/mK]
lambda_S = 1.0 # thermal conductivity of solid matrix [w/mK]
dh_evap = 2258000 # latent heat of evaporation [J/kg]
D_pm = 2.6e6 # binary diffusion coefficient [m²/s]
rho_L = 1000.0 # density of liquid phase [kg/m³]
MW = 0.018016 # molecular weight of water component [kg/mol]
MC = 0.028949 # molecular weight of air component [kg/mol]
R = 8.3144621 # universal gas constant
mu_L = 2.938e4 # dynamic viscosity of liquid phase [Pa.s]
muA_G = 2.194e5 # dynamic viscosity of air component in gas phase [Pa.s]
muW_G = 1.227e5 # dynamic viscosity of water component in gas phase [Pa.s]
s_LRes = 0.0 # residual saturation of liquid phase []
s_GRes = 0.0 # residual saturation of gas phase []
k_rG_min = 1e5 # used for normalization of BC model
k_rL_min = 1e5 # used for normalization of BC model
p_thr_BC = 5.0e3 # entry pressure for BrooksCorey model [Pa]
exp_BC = 3.0 # Corey exponent for BrooksCorey model []
The system of balance equations from the previous chapter has to be closed by a variety of constitutive relations, which in turn specify the necessary material properties and express the thermodynamic equilibrium between the constituents of liquid and gas phases. We apply the same relations in the analytical solution and in the numerical model as well.
Due to the existence of multiple phases within the pore space, the movement of a ﬂuid phase is obstructed by the presence of the other phase. In multiphase ﬂow applications, this effect is usually realised by introducing relative permeabilities as functions of saturation which calculate the effective permeability of each phase as described in the extended Darcy law. Additionally if the present phases are immiscible, one also needs to consider the arising capillary effects by introducing capillary pressure accounting for the difference of phase pressures.
The BrooksCorey formulation accounts for residual saturations of both ﬂuid phases and links effective saturation to capillary pressure. Additionally, it also models the relative permeability of both phases as a function of effective liquid saturation.
def capillary_pressure(sL_eff):
return p_thr_BC * (sL_eff ** (1.0 / exp_BC))
def capillary_pressure_derivative(sL_eff):
return p_thr_BC / exp_BC * (sL_eff ** ((exp_BC + 1.0) / exp_BC))
def saturation_effective(p_c):
return (p_c / p_thr_BC) ** (exp_BC)
def relative_permeability_gas(sL_eff):
return max(
k_rG_min, ((1.0  sL_eff) ** 2) * (1  (sL_eff ** ((2.0 + exp_BC) / exp_BC)))
)
def relative_permeability_liquid(sL_eff):
return max(k_rL_min, sL_eff ** ((2.0 + 3 * exp_BC) / exp_BC))
Determining the composition of the gas phase, we assume that the sum of all constituents’ partial pressures accounts for the entire gas phase pressure (Dalton’s law). The partial pressure of water vapour is derived from the true vapour pressure that accounts for the impact of capillary effects as well: due to wettability and capillarity, the interfaces prevailing in porous media are not flat, but rather curving. Above curved interfaces, the vapour pressure may change depending on the direction of curvature. The KelvinLaplace equation accounts for this and expresses the true vapour pressure as a function of capillary pressure and the saturation vapour pressure of pure water.
The saturation vapour pressure is determined by the approximate ClausiusClapeyron equation.
def vapour_pressure(p_sat, p_G, p_c, xA_G, T):
return p_sat * math.exp((p_c  xA_G * p_G) * MW / rho_L / R / T)
def saturation_vapour_pressure(T):
return p_ref * math.exp((1.0 / T_ref  1.0 / T) * dh_evap * MW / R)
def partial_pressure_vapour(p_G, p_c, xA_G, T):
p_sat = saturation_vapour_pressure(T)
return vapour_pressure(p_sat, p_G, p_c, xA_G, T)
The gas phase density can be determined using an appropriate equation of state for binary mixtures, such as the PengRobinson equation of state. For simplicity, we use the thermal equation of state of ideal gases in this work, which gives sufﬁcient results at high temperatures and low pressures.
Dynamic viscosity of composite phases is obtained from the simple mixing rule.
def molar_mass_gas_phase(xA_G):
return xA_G * MC + (1  xA_G) * MW
def density_gas_phase(p_G, xA_G, T):
M = molar_mass_gas_phase(xA_G)
return p_G * M / (R * T)
def viscosity_gas_phase(xA_G):
return xA_G * muA_G + (1.0  xA_G) * muW_G
def kinematic_viscosity_gas_phase(p_G, xA_G, T):
mu_G = viscosity_gas_phase(xA_G)
rho_G = density_gas_phase(p_G, xA_G, T)
return mu_G / rho_G
The macroscopic diffusion constant (diffusivity) is determined using a simple correction accounting for porosity and gas saturation on the microscopic diffusion constant ($D_{pm}^G$).
def diffusivity(sL_eff):
return phi * (1.0  sL_eff) * D_pm
When heat conduction occurs over multiple phases, a mixing rule can describe averaged heat conduction if local thermal equilibrium is assumed. In this case, we apply a very simple model (upper Wiener bound, Wiener 1912) to ﬁnd an effective heat conductivity by averaging individual phase conductivities by volume fraction.
def thermal_conductivity(sL_eff):
sL = sL_eff * (1.0  s_GRes  s_LRes) + s_LRes
phi_G = (1.0  sL) * phi
phi_L = sL * phi
phi_S = 1.0  phi
return lambda_G * phi_G + lambda_L * phi_L + lambda_S * phi_S
The solution for the 4 primary variables (liquid phase saturation, gas phase pressure, mole fraction of air in the gas phase and temperature) is obtained from the numerical integration of a coupled system of first order differential equations.
\begin{align} \frac{\partial z}{\partial S_{L,eff}}=\frac{\frac{dp_c}{dS_{L,eff}}}{\theta \omega \frac{1}{k_{rG}}( \frac{1}{1x_G^a} + \frac{\nu_L}{\nu_G} \frac{k_{rG}}{k_{rL}} )} \end{align}
\begin{align} \frac{\partial p_G}{\partial S_{L,eff}}=\frac{\frac{dp_c}{dS_{L,eff}}}{1 + \frac{\nu_L}{\nu_G} \frac{k_{rG}}{k_{rL}} ({1x_G^a})} \end{align}
\begin{align} \frac{\partial x_G^a}{\partial S_{L,eff}}=\frac{\frac{dp_c}{dS_{L,eff}} \frac{K}{\nu_G \rho_G D_{pm}} \frac{1}{1x_G^a}}{\frac{1}{k_{rG}} ( \frac{1}{1x_G^a} + \frac{\nu_L}{\nu_G} \frac{k_{rG}}{k_{rL}} )} \end{align}
\begin{align} \frac{\partial T}{\partial S_{L,eff}}=\frac{\frac{dp_c}{dS_{L,eff}} \frac{1\theta}{\theta} \frac{dh_{evap} K}{\nu_g D_{pm}}}{\frac{1}{k_{rG}} ( \frac{1}{1x_G^a} + \frac{\nu_L}{\nu_G} \frac{k_{rG}}{k_{rL}} )} \end{align}
The differentail equations given by (Udell and Fitch, 1985) that are shown above are derived on the basis of mass and energy balance. It should be noted that the equation system is integrated over the effective saturation instead of the spatial coordinate $z$. The above given equations are reached when:
\begin{align} \frac{dz}{dS_{L,eff}}=\frac{dp_c/dS_{L,eff}}{\partial p_c / \partial z} \end{align}
\begin{align} \frac{\partial p_G}{\partial S_{L,eff}}=\frac{\partial p_G}{\partial z} \frac{\partial z}{\partial S_{L,eff}} \end{align}
The parameter $\eta$ in the energy balance equation represents the ratio of the heat flux caused by convection to the total heat flux $q$:
\begin{align} \eta = 1+\frac{\lambda \frac{\partial T}{\partial z}}{q} \end{align}
With some algebraic transformation, $\eta$ can also be explicitly expressed by
\begin{align} \alpha = 1+ \frac{p_c  x_G^a p_G}{dh_{evap} \rho_L} \end{align}
\begin{align} \delta = \frac{\rho_L dh_{evap}^2 K \alpha}{\lambda \nu_G T} \end{align}
\begin{align} \xi = \frac{1}{k_{rG}} (1+ \frac{\rho_L R T}{p_G M^w} \frac{1}{1x_G^a})+ \frac{\nu_L}{\nu_G k_{rL}} \end{align}
\begin{align} \zeta = \frac{K \rho_L R T}{\rho_G \nu_G D_{pm} M^w} \frac{x_G^a}{1x_G^a} (\frac{p_G M^w}{\rho_L R T} + \frac{1}{1x_G^a}) \end{align}
\begin{align} \eta = \frac{\delta}{\delta \xi \zeta} \end{align}
Including the reformulated $\eta$ and grouping some of the parameters under $\gamma$ and $\omega$, the 4 governing equation can be calculated in the following form:
\begin{align} \gamma = \frac{1}{k_{rG}} ( \frac{1}{1x_G^a} + \frac{\nu_L}{\nu_G} \frac{k_{rG}}{k_{rL}} ) \end{align}
\begin{align} \omega = \frac{q \nu_G}{dh_{evap} K} \end{align}
\begin{align} \frac{\partial z}{\partial S_{L,eff}} = \frac{\frac{dp_c}{dS_{L,eff}}}{\eta \omega \gamma} \end{align}
\begin{align} \frac{\partial p_G}{\partial S_{L,eff}} = \frac{\frac{dp_c}{dS_{L,eff}}}{\gamma k_{rG} (1x_G^a)} \end{align}
\begin{align} \frac{\partial x_G^a}{\partial S_{L,eff}} = \frac{\frac{dp_c}{dS_{L,eff}} \frac{K}{\nu_G \rho_G D_{pm}} \frac{1}{1x_G^a}}{\gamma} \end{align}
\begin{align} \frac{\partial T}{\partial S_{L,eff}} = \frac{\frac{dp_c}{dS_{L,eff}} \frac{1\theta}{\theta} \frac{dh_{evap} K}{\nu_g D_{pm}}}{\gamma} \end{align}
The full derivation of the analytical solution can be found in (Helmig et al., 1997).
# Parameter grouping
def alpha_(sL_eff, p_G, xA_G):
p_c = capillary_pressure(sL_eff)
return 1.0 + ((p_c  xA_G * p_G) / (dh_evap * rho_L))
def delta_(sL_eff, p_G, xA_G, T):
alpha = alpha_(sL_eff, p_G, xA_G)
nu_G = kinematic_viscosity_gas_phase(p_G, xA_G, T)
th_cond = thermal_conductivity(sL_eff)
return (rho_L * (dh_evap**2) * K * alpha) / (th_cond * nu_G * T)
def xi_(sL_eff, p_G, xA_G, T):
k_rG = relative_permeability_gas(sL_eff)
k_rL = relative_permeability_liquid(sL_eff)
nu_G = kinematic_viscosity_gas_phase(p_G, xA_G, T)
nu_L = mu_L / rho_L
return ((1.0 + ((rho_L * R * T) / (p_G * MW * (1.0  xA_G)))) / k_rG) + (
nu_L / nu_G
) / k_rL
def zeta_(sL_eff, p_G, xA_G, T):
D = diffusivity(sL_eff)
mu_G = viscosity_gas_phase(xA_G)
a = K * rho_L * R * T * xA_G
b = mu_G * D * MW * (1.0  xA_G)
c = p_G * MW / (rho_L * R * T) + 1.0 / (1.0  xA_G)
return (a / b) * c
def eta_(sL_eff, p_G, xA_G, T):
delta = delta_(sL_eff, p_G, xA_G, T)
xi = xi_(sL_eff, p_G, xA_G, T)
zeta = zeta_(sL_eff, p_G, xA_G, T)
return delta / (delta + xi + zeta)
def gamma_(sL_eff, p_G, xA_G, T):
k_rG = relative_permeability_gas(sL_eff)
k_rL = relative_permeability_liquid(sL_eff)
nu_G = kinematic_viscosity_gas_phase(p_G, xA_G, T)
nu_L = mu_L / rho_L
return 1.0 / k_rG * ((1.0 / (1.0  xA_G)) + (nu_L / nu_G) * (k_rG / k_rL))
# Differential equations
# Spatial variable (1D) derivative
def dz_dsL_eff(sL_eff, p_G, xA_G, T):
dpC_dsL_eff = capillary_pressure_derivative(sL_eff)
eta = eta_(sL_eff, p_G, xA_G, T)
gamma = gamma_(sL_eff, p_G, xA_G, T)
nu_G = kinematic_viscosity_gas_phase(p_G, xA_G, T)
omega = (q * nu_G) / (dh_evap * K)
return dpC_dsL_eff / (eta * omega * gamma)
# Gasphase pressure derivative
def dp_G_dsL_eff(sL_eff, p_G, xA_G, T):
dpC_dsL_eff = capillary_pressure_derivative(sL_eff)
gamma = gamma_(sL_eff, p_G, xA_G, T)
k_rG = relative_permeability_gas(sL_eff)
return dpC_dsL_eff / (gamma * k_rG * (1.0  xA_G))
# Mole fraction of air component in gas phase derivative
def dxA_G_dsL_eff(sL_eff, p_G, xA_G, T):
dpC_dsL_eff = capillary_pressure_derivative(sL_eff)
gamma = gamma_(sL_eff, p_G, xA_G, T)
mu_G = viscosity_gas_phase(xA_G)
D = diffusivity(sL_eff)
return dpC_dsL_eff * K / (mu_G * D) * xA_G / (1.0  xA_G) / gamma
# Temperature derivative
def dT_dsL_eff(sL_eff, p_G, xA_G, T):
dpC_dsL_eff = capillary_pressure_derivative(sL_eff)
eta = eta_(sL_eff, p_G, xA_G, T)
gamma = gamma_(sL_eff, p_G, xA_G, T)
nu_G = kinematic_viscosity_gas_phase(p_G, xA_G, T)
th_cond = thermal_conductivity(sL_eff)
return dpC_dsL_eff * (1.0  eta) / eta * dh_evap / (nu_G * th_cond) * K / gamma
To obtain the analytical solution, the above introduced four coupled differential equation can now be integrated with the wellknow Forward Euler method.
To get a unique solution, the following boundary and initial conditions are considered:
On the lefthand side ($z = 0$)  the cool end: $S_L = 1$ , $p_G = p_{G,i}$ , $x_G^a = x_{G,i}^a$ and $T = T_i$.
On the righthand side  the hot end of the the heat pipe, a constant heatflux is considered: $q = 100 W/m^2$.
As initial conditions we chose: $S_L = 1$ , $p_G = 101325 Pa$ , $x_G^a = 0.25$ (computed according to Dalton’s law using the true vapour pressure assuming $p_{c,i} = 5001 Pa$) and $T = 365 K$.
# Define righthandsides of the coupled system of first order derivative equations
def dy_dsL_eff(_y, sL_eff, p_G, xA_G, T):
dydsL = np.zeros(4)
dydsL[0] = dz_dsL_eff(sL_eff, p_G, xA_G, T)
dydsL[1] = dp_G_dsL_eff(sL_eff, p_G, xA_G, T)
dydsL[2] = dxA_G_dsL_eff(sL_eff, p_G, xA_G, T)
dydsL[3] = dT_dsL_eff(sL_eff, p_G, xA_G, T)
return dydsL
# Numerical integration  Forward Euler method
# to estimate the integrals of the coupled equation system
def step_Euler(y, sL_eff, dsL_eff, p_G, xA_G, T):
return y + dsL_eff * dy_dsL_eff(y, sL_eff, p_G, xA_G, T)
def full_Euler(dsL_eff, y0, sL_eff_low, sL_eff_high):
max_steps = int(abs((sL_eff_low  sL_eff_high) / dsL_eff))
sL_eff_list = np.linspace(sL_eff_low, sL_eff_high, max_steps + 1)
M = np.zeros(
(4, max_steps + 1)
) # Solution matrix containing the 4 primary variable
M[:, 0] = y0
for i in range(max_steps):
p_G = M[1, i]
xA_G = M[2, i]
T = M[3, i]
if (dz_dsL_eff(sL_eff_list[i], p_G, xA_G, T) * dsL_eff) < 1.0:
M[:, i + 1] = step_Euler(M[:, i], sL_eff_list[i], dsL_eff, p_G, xA_G, T)
else:
M[:, i + 1] = np.nan
return M, sL_eff_list
# initial condition
z_0 = 0
p_G0 = 101325
p_c0 = 5001
T_0 = 365
xA_G0 = 1  partial_pressure_vapour(p_G0, p_c0, 0, T_0) / p_G0
y0 = np.array([z_0, p_G0, xA_G0, T_0])
# initial effective saturation
sL_eff_0 = saturation_effective(p_c0)
# integration boundaries and saturation step size
sL_eff_low = sL_eff_0
sL_eff_high = 10e16
n_dsL_eff = 10**2
dsL_eff = (sL_eff_high  sL_eff_low) / n_dsL_eff
# execute analytical solution
M, sL_eff_list = full_Euler(dsL_eff, y0, sL_eff_low, sL_eff_high)
In the following section, the previously analitically solved heatpipe problem is solved using numerical methods (FEM in OpenGeoSys). The numerical solution takes the exact same input parameters, and considers the same  modified  assumptions that are mentioned in the analytical solution.
Additionally, it is necessary to introduce a discretized spacial domain for the numerical simulation. We chose a 1 dimensional domain with length of 1 m discretized into 200 equally spaced identical elements.
The numerical problem considers the same consitutive relationships and identical boundary and initial conditions. We use the TH2M model of OGS to solve the coupled partial differential equations describing the system behavior. Detailed description of the numerical model can be found in (Grunwald et al., 2022).
import os
from pathlib import Path
from ogs6py import ogs
out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
out_dir.mkdir(parents=True)
prj_file = "heat_pipe_rough.prj"
model = ogs.OGS(INPUT_FILE=prj_file, PROJECT_FILE=prj_file)
model.run_model(logfile=f"{out_dir}/out.txt", args=f"o {out_dir}")
OGS finished with project file heat_pipe_rough.prj.
Execution took 93.63611078262329 s
Project file written to output.
To compare the results produced by the analytical solution and those by the numerical compution by OpenGeoSys, the four primary variables are plotted along the 1D domain ($z$): liquid phase saturation ($S_{L,eff}$), air molar fraction in the gas phase ($x_G^a$), temperature ($T$) and gas phase pressure ($p_G$).
# Import OGS simulation results
import pyvista as pv
pv.set_plot_theme("document")
pv.set_jupyter_backend("static")
pvd_file = f"{out_dir}/results_heatpipe_rough.pvd"
reader = pv.get_reader(pvd_file)
reader.set_active_time_value(1.0e7) # set reader to simulation endtime
mesh = reader.read()[0]
# Define line along mesh and extract data along line for plotting
pt1 = (0, 0.0025, 0)
pt2 = (1, 0.0025, 0)
xaxis = pv.Line(pt1, pt2, resolution=2)
line_mesh = mesh.slice_along_line(xaxis)
x_num = line_mesh.points[:, 0] # x coordinates of each point
S_num = line_mesh.point_data["saturation"]
xA_G_num = line_mesh.point_data["xnCG"]
p_G_num = line_mesh.point_data["gas_pressure"]
T_num = line_mesh.point_data["temperature"]
# Resampling dataset via linear interpolation for error calculation
S_num_interp = np.interp(M[0, :], x_num, S_num)
xA_G_num_interp = np.interp(M[0, :], x_num, xA_G_num)
p_G_num_interp = np.interp(M[0, :], x_num, p_G_num)
T_num_interp = np.interp(M[0, :], x_num, T_num)
As one can see from the figures below, the numerical results are in really good agreement with the analytical solution. To better understand and visualize the deviation, we also perform a quick error analysis by simply calculating the difference (absolute and relative error) between the analytical and the numerical solution.
import matplotlib.pyplot as plt
plt.rcParams["lines.linewidth"] = 2.0
plt.rcParams["lines.color"] = "black"
plt.rcParams["legend.frameon"] = True
plt.rcParams["font.family"] = "serif"
plt.rcParams["legend.fontsize"] = 14
plt.rcParams["font.size"] = 14
plt.rcParams["axes.axisbelow"] = True
plt.rcParams["figure.figsize"] = (16, 6)
fig1, (ax1, ax2, ax3) = plt.subplots(1, 3)
ax1.plot(M[0, :], sL_eff_list, "kx", label=r"$S_{L,eff}$ analytical")
ax1.plot(M[0, :], M[2, :], "kx", label=r"$x_G^a$ analytical")
ax1.plot(x_num, S_num, "b", label=r"$S_{L,eff}$ numerical")
ax1.plot(x_num, xA_G_num, "g", label=r"$x_G^a$ numerical")
ax1.set_xlabel(r"$z$ / m")
ax1.set_ylabel(r"$S_{L,eff}$ and $x_G^a$ / ")
ax1.legend()
ax1.grid(True)
ax1.set_xlim(0, 1)
ax1.set_ylim(0, 1)
ax2.plot(M[0, :], S_num_interp  sL_eff_list, "b", label=r"$\Delta S_{L,eff}$")
ax2.plot(M[0, :], xA_G_num_interp  M[2, :], "g", label=r"$\Delta x_G^a$")
ax2.set_xlabel(r"$z$ / m")
ax2.set_ylabel(r"Absolute error / ")
ax2.legend()
ax2.grid(True)
ax2.set_xlim(0, 1)
ax2.set_ylim(0.001, 0.02)
relError_S_w = np.zeros(len(M[0, :]))
relError_xA_G = np.zeros(len(M[0, :]))
for i in range(len(M[0, :])):
if (sL_eff_list[i]) >= 0.001:
relError_S_w[i] = (S_num_interp[i]  sL_eff_list[i]) / sL_eff_list[i]
else:
relError_S_w[i] = np.nan
if (M[2, i]) >= 0.01:
relError_xA_G[i] = (xA_G_num_interp[i]  M[2, i]) / M[2, i]
else:
relError_xA_G[i] = np.nan
ax3.plot(M[0, :], relError_S_w, "b", label=r"$\Delta S_{L,eff}/S_{L,effanalytical}$")
ax3.plot(M[0, :], relError_xA_G, "g", label=r"$\Delta x_G^a/x_{Ganalytical}^a$")
ax3.set_xlabel(r"$z$ / m")
ax3.set_ylabel(r"Relative error / ")
ax3.set_xlim(0, 1)
ax3.legend()
ax3.grid(True)
fig1.tight_layout()
plt.show()
fig2, (ax1, ax2, ax3) = plt.subplots(1, 3)
ax1.plot(M[0, :], M[3, :], "kx", label=r"$T$ analytical")
ax1.plot(x_num, T_num, "r", label=r"$T$ numerical")
ax1.set_xlabel(r"$z$ / m")
ax1.set_ylabel(r"$T$ / K")
ax1.legend()
ax1.grid(True)
ax1.set_xlim(0, 1)
ax1.set_ylim(365, 375)
ax2.plot(M[0, :], T_num_interp + M[3, :], "r", label=r"$\Delta T$")
ax2.set_xlabel(r"$z$ / m")
ax2.set_ylabel(r"Absolute error / K")
ax2.legend()
ax2.grid(True)
ax2.set_xlim(0, 1)
ax2.set_ylim(0, 0.12)
ax3.plot(
M[0, :],
(T_num_interp + M[3, :]) / M[3, :],
"r",
label=r"$\Delta T/T_{analytical}$",
)
ax3.set_xlabel(r"$z$ / m")
ax3.set_ylabel(r"Relative error / ")
ax3.set_xlim(0, 1)
ax3.set_ylim(0, 0.0003)
ax3.legend()
ax3.grid(True)
fig2.tight_layout()
plt.show()
fig3, (ax1, ax2, ax3) = plt.subplots(1, 3)
ax1.plot(M[0, :], M[1, :] / 1000, "kx", label=r"$p_G$ analytical")
ax1.plot(x_num, p_G_num / 1000, "c", label=r"$p_G$ numerical")
ax1.set_xlabel(r"$z$ / m")
ax1.set_ylabel(r"$p_G$ / kPa")
ax1.legend()
ax1.grid(True)
ax1.set_xlim(0, 1)
ax1.set_ylim(100, 106)
ax2.plot(M[0, :], p_G_num_interp + M[1, :], "c", label=r"$\Delta p_G$")
ax2.set_xlabel(r"$z$ / m")
ax2.set_ylabel(r"Absolute error / Pa")
ax2.legend()
ax2.grid(True)
ax2.set_xlim(0, 1)
ax2.set_ylim(0, 30)
ax3.plot(
M[0, :],
(p_G_num_interp + M[1, :]) / M[1, :],
"c",
label=r"$\Delta p_G/p_{Ganalytical}$",
)
ax3.set_xlabel(r"$z$ / m")
ax3.set_ylabel(r"Relative error / ")
ax3.set_xlim(0, 1)
ax3.set_ylim(0, 0.0004)
ax3.legend()
ax3.grid(True)
fig3.tight_layout()
plt.show()
[1] Helmig, R. et al. (1997). Multiphase ﬂow and transport processes in the subsurface: a contribution to the modeling of hydrosystems. SpringerVerlag.
[2] Udell, K. and Fitch, J. (1985) Heat and mass transfer in capillary porous media considering evaporation, condensation, and noncondensible gas effects. 23rd ASME/AIChE national heat transfer conference, Denver, pp. 103110.
[3] Wiener, O. (1912). “Abhandl. math.phys”. In: Kl. Königl. Sächsischen Gesell 32, p. 509.
[4] Grunwald, N. et al. (2022). Nonisothermal twophase ﬂow in deformable porous media: Systematic opensource implementation and veriﬁcation procedure. Geomech. Geophys. Geoenerg. Georesour., 8, 107 https://link.springer.com/article/10.1007/s40948022003942
This article was written by Kata Kurgyis. If you are missing something or you find an error please let us know.
Generated with Hugo 0.122.0
in CI job 447067

Last revision: September 14, 2022