3. Constitutive laws#
Equation of state for isothermal compressible fluids#
GaPFlow is a compressible fluid solver for confined flows. To close the gap-averaged balance laws, one needs to specify an equation of state in the form \(p(\rho)\).
Below is an overview of fluid models currently implemented GaPFlow:
import numpy as np
import matplotlib.pyplot as plt
from GaPFlow.models.pressure import eos_pressure
def plot_bwr(ax):
props = {'EOS': 'BWR'}
dens = np.linspace(0.2, 0.8, 100)
for temp in [1.5, 2.0, 2.5]:
props['T'] = temp
p = eos_pressure(dens, props)
ax.plot(dens, p, label=rf"$T=${temp}")
ax.legend()
ax.set_xlabel(r'Density')
ax.set_ylabel(r'Pressure')
ax.set_title('Benedict-Webb-Rubin (LJ units)')
def plot_dh(ax):
props = {'EOS': 'DH'}
dens = np.linspace(870., 1050., 100)
p = eos_pressure(dens, props)
ax.plot(dens, p, label='Dowson-Higginson')
ax.set_xlabel(r'Density ($\mathrm{kg\,m^{-3}}$)')
ax.set_ylabel(r'Pressure ($\mathrm{Pa}$)')
ax.set_title('Dowson-Higginson')
def plot_mt(ax):
props = {'EOS': 'MT'}
dens = np.linspace(700., 1200., 100)
p = eos_pressure(dens, props)
ax.plot(dens, p, label='Murnaghan-Tait')
ax.set_xlabel(r'Density ($\mathrm{kg\,m^{-3}}$)')
ax.set_ylabel(r'Pressure ($\mathrm{Pa}$)')
ax.set_title('Murnaghan-Tait')
def plot_pl(ax):
props = {'EOS': 'PL'}
dens = np.linspace(0., 2., 100)
for alpha in [0., 0.5, 1., 1.5]:
props['alpha'] = alpha
p = eos_pressure(dens, props)
ax.plot(dens, p, label=rf'$\alpha=${alpha:.1f}')
ax.set_xlabel(r'Density ($\mathrm{kg\,m^{-3}}$)')
ax.set_ylabel(r'Pressure ($\mathrm{Pa}$)')
ax.set_title('Power law')
ax.legend()
def plot_vdW(ax):
props = {'EOS': 'vdW'}
dens = np.linspace(0., 1000., 100)
for temp in [120., 150., 180., 210.]:
props['T'] = temp
p = eos_pressure(dens, props)
ax.plot(dens, p, label=rf"$T=${temp:.0f} K")
ax.set_xlabel(r'Density ($\mathrm{kg\,m^{-3}}$)')
ax.set_ylabel(r'Pressure ($\mathrm{Pa}$)')
ax.set_title('van der Waals')
ax.legend()
def plot_cubic(ax):
props = {'EOS': 'cubic'}
dens = np.linspace(0.2, 0.6, 100)
p = eos_pressure(dens, props)
ax.plot(dens, p, '--', color='0.0', label='cubic (T=2)')
ax.legend()
if __name__ == "__main__":
sx, sy = plt.rcParams['figure.figsize']
funcs = [plot_bwr, plot_dh, plot_mt, plot_pl, plot_vdW]
n = len(funcs)
nr = n // 2 + n % 2
fig, axes = plt.subplots(nr, 2, figsize=(1.5*sx, nr * sy))
for ax, func in zip(axes.flat, funcs):
func(ax)
plot_cubic(axes[0, 0])
plt.show()
Viscosity models for non-Newtonian fluids#
The viscosity may depend on the pressure and the shear rate. We employ a so-called generalized Newtonian fluid, where we assume that the velocity profiles are still those of a Newtonian fluid, but the viscosity parameter is the non-Newtonian one. Below is an overview of implemented shear-thinning and piezoviscosity models:
import io
import numpy as np
import matplotlib.pyplot as plt
from importlib import resources
from GaPFlow.io import read_yaml_input
from GaPFlow.models.viscosity import shear_thinning_factor, piezoviscosity
eyring = """
properties:
shear: 0.1
bulk: 0.
EOS: DH
thinning:
name: Eyring
"""
carreau = """
properties:
shear: 0.1
bulk: 0.
EOS: DH
thinning:
name: Carreau
"""
barus = """
properties:
shear: 0.1
bulk: 0.
EOS: DH
piezo:
name: Barus
"""
roelands = """
properties:
shear: 0.1
bulk: 0.
EOS: DH
piezo:
name: Roelands
"""
if __name__ == "__main__":
fig, ax = plt.subplots(1, 2, figsize=(1.5*sx, sy))
fig.suptitle('Viscosity models (with default parameters)')
# Shear thinning
shear_rate = np.logspace(0, 11, 100)
for model in [eyring, carreau]:
with io.StringIO(model) as file:
inp = read_yaml_input(file)['properties']
mu0 = inp['shear']
name = inp['thinning']['name']
viscosity = mu0 * shear_thinning_factor(shear_rate, mu0, inp['thinning'])
ax[0].loglog(shear_rate, viscosity, label=name)
ax[0].set_xlabel(r'Shear rate ($\mathrm{s^{-1}}$)')
ax[0].set_ylabel(r'Viscosity ($\mathrm{Pa\,s}$)')
ax[0].legend()
# Piezoviscosity
pressure = np.linspace(1e6, 1e9, 100)
for model in [barus, roelands]:
with io.StringIO(model) as file:
inp = read_yaml_input(file)['properties']
mu0 = inp['shear']
name = inp['piezo']['name']
viscosity = mu0 * piezoviscosity(pressure, mu0, inp['piezo'])
ax[1].loglog(pressure, viscosity, label=name)
ax[1].set_xlabel(r'Pressure ($\mathrm{Pa}$)')
ax[1].set_ylabel(r'Viscosity ($\mathrm{Pa\,s}$)')
ax[1].legend()
plt.show()
*************************************************************
* PROBLEM SETUP *
*************************************************************
- options:
- grid:
- geometry:
- numerics:
- properties:
- shear : 0.1
- bulk : 0.0
- EOS : DH
- rho0 : 877.7007
- P0 : 101325.0
- C1 : 35000000000.0
- C2 : 1.23
- thinning:
- name : Eyring
- tauE : 500000.0
- elastic:
- enabled : False
- gp:
- db:
- md:
*************************************************************
* PROBLEM SETUP COMPLETED *
*************************************************************
*************************************************************
* PROBLEM SETUP *
*************************************************************
- options:
- grid:
- geometry:
- numerics:
- properties:
- shear : 0.1
- bulk : 0.0
- EOS : DH
- rho0 : 877.7007
- P0 : 101325.0
- C1 : 35000000000.0
- C2 : 1.23
- thinning:
- name : Carreau
- mu_inf : 1e-09
- lam : 1e-06
- a : 2.0
- N : 0.6
- elastic:
- enabled : False
- gp:
- db:
- md:
*************************************************************
* PROBLEM SETUP COMPLETED *
*************************************************************
*************************************************************
* PROBLEM SETUP *
*************************************************************
- options:
- grid:
- geometry:
- numerics:
- properties:
- shear : 0.1
- bulk : 0.0
- EOS : DH
- rho0 : 877.7007
- P0 : 101325.0
- C1 : 35000000000.0
- C2 : 1.23
- piezo:
- name : Barus
- aB : 2e-08
- elastic:
- enabled : False
- gp:
- db:
- md:
*************************************************************
* PROBLEM SETUP COMPLETED *
*************************************************************
*************************************************************
* PROBLEM SETUP *
*************************************************************
- options:
- grid:
- geometry:
- numerics:
- properties:
- shear : 0.1
- bulk : 0.0
- EOS : DH
- rho0 : 877.7007
- P0 : 101325.0
- C1 : 35000000000.0
- C2 : 1.23
- piezo:
- name : Roelands
- mu_inf : 0.001
- p_ref : 196000000.0
- z : 0.68
- elastic:
- enabled : False
- gp:
- db:
- md:
*************************************************************
* PROBLEM SETUP COMPLETED *
*************************************************************
