Confined compression of a cube

In this test, the thermodynamic relationships between gas pressure, temperature and density are tested. For that, a cube-shaped domain consisting of an ideal gas is compressed by 50% of its initial volume over a short period of time, starting from the top surface. The boundaries of that domain are impermeable to fluid flow, therefore gas pressure and density must increase as a result of the decreasing volume. Since energy flow across the boundaries is also prevented, this compression is an adiabatic change of state. No frictional losses are taken into account, so the process can be reversed at any time and the entropy in the overall system remains constant.

Analytical solution

Density evolution

The mass balance for such a system can be found by simplifying the mass balance (eq. 44 in Grunwald et al., 2022). With $\phi=s_\mathrm{G}=\alpha_\mathrm{B}=1$ and $\mathrm{A}^\zeta_\alpha=\mathrm{J}^\zeta_\alpha=0$ one obtains

With volume strain $e=\text{div}\left(\mathbf{u}_\text{S}\right)'_\text{S}$ we write

which can be integrated

so that we find

or

Gas pressure evolution

The evolution of gas pressure can be found from energy balance equation (eq. 51 in the paper): Starting from

with $u_\text{S}=h_\text{S}$ and $u_\text{G}=h_\text{G}-\frac{p_\text{GR}}{\rho_\text{GR}}$ and with $h_\alpha=c_{p,\alpha}T$. Considering that $\phi_\alpha=\text{const}$ and $c_{p,\text{S}}=0$, we find

Assuming ideal gas behaviour, we can write $\frac{p_\text{GR}}{\rho_\text{GR}}=\frac{M}{RT}$ and find

With $c_{v,\text{G}}=\text{const}$ it follows

Density and temperature cancel out so the final equation can be written as

with the adiabatic index $\kappa = \frac{c_v}{c_p}$. The equation can be integrated to find the solution for the pressure evolution

Temperature evolution

The temperature evolution follows Poissons equations for isentropic processes

import os
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import ogstools as ot

out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
    out_dir.mkdir(parents=True)

p_0 = 1e6
T_0 = 270
R = 8.3144621
M = 0.01
c_p = 1000

c_v = c_p - R / M
rho_0 = p_0 * M / R / T_0
kappa = c_p / c_v


def volume_strain(t):
    return -t / 100


def gas_density(t):
    return rho_0 * np.exp(-volume_strain(t))


def gas_pressure(t):
    return p_0 * np.exp(-kappa * volume_strain(t))


def temperature(t):
    return gas_pressure(t) * M / R / gas_density(t)

Numerical solution

# run OGS
cube_compression = ot.Project(
    input_file="compression_gas.prj", output_file="compression_gas.prj"
)
cube_compression.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir}")

OGS finished with project file compression_gas.prj.
Execution took 0.31472015380859375 s
Project file written to output.

# read PVD file
ms = ot.MeshSeries(f"{out_dir}/result_compression_gas.pvd")
time = ms.timevalues()


def plot_results_errors(var: ot.variables.Scalar, ref_vals: np.ndarray):
    "Plot numerical results, analytical solution and errors."
    num_vals = ms.probe([0.0, 1.0, 0.0], var.data_name)[:, 0]
    abs_err = ref_vals - num_vals
    rel_err = abs_err / ref_vals
    assert np.all(np.abs(rel_err) <= 0.002)

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
    ax3 = ax2.twinx()
    ax1.set_xlabel(r"$t$ / s")
    ax2.set_xlabel(r"$t$ / s")
    ax1.set_ylabel(var.get_label())
    ax2.set_ylabel("absolute error / " + var.data_unit)
    ax3.set_ylabel("relative error")

    ax1.plot(time, var.transform(num_vals), "kx", label="numerical")
    ax1.plot(time, var.transform(ref_vals), "b", label="analytical")
    lns_abs = ax2.plot(time, abs_err, "b", label=r"absolute")
    lns_rel = ax3.plot(time, rel_err, "g", label=r"relative")
    ax1.legend()
    lns = lns_abs + lns_rel
    ax2.legend(lns, [label.get_label() for label in lns])

    fig.tight_layout()

Gas density evolution

rho_g = ot.variables.Scalar("gas_density", r"kg/$m^3$", r"kg/$m^3$")
plot_results_errors(rho_g, gas_density(ms.timevalues()))

  0%|          | 0/11 [00:00<?, ?it/s]

Gas pressure evolution

p_g = ot.variables.Scalar("gas_pressure_interpolated", "Pa", "MPa")
plot_results_errors(p_g, gas_pressure(ms.timevalues()))

  0%|          | 0/11 [00:00<?, ?it/s]

Temperature evolution

plot_results_errors(ot.variables.temperature, temperature(ms.timevalues()))

  0%|          | 0/11 [00:00<?, ?it/s]