An Adaptive ALE Scheme for Non-Ideal Compressible-Fluid Dynamics over Dynamic Unstructured Meshes

TL;DR

This paper presents an adaptive mesh scheme for simulating non-ideal compressible fluid flows near vapor-liquid saturation, improving robustness and accuracy by avoiding solution interpolation during grid adaptation.

Contribution

It introduces a novel interpolation-free finite-volume adaptive scheme within an ALE framework for NICFD, incorporating non-ideal gas effects via the Peng-Robinson model.

Findings

Successfully captures shock waves and wave interactions in NICFD flows.

Demonstrates robustness and accuracy across various piston movements and initial conditions.

Applicable to 2D and 3D simulations of complex non-ideal flow phenomena.

Abstract

This paper investigates the application of mesh adaptation techniques in the Non-Ideal Compressible Fluid Dynamic (NICFD) regime, a region near the vapor-liquid saturation curve where the flow behavior significantly departs from the ideal gas model, as indicated by a value of the fundamental derivative of gasdynamics less than one. A recent interpolation-free finite-volume adaptive scheme is exploited to modify the grid connectivity in a conservative way, and the governing equations for compressible inviscid flows are solved within the Arbitrary Lagrangian Eulerian framework by including special fictitious fluxes representing volume modifications due to mesh adaptation.The absence of interpolation of the solution to the new grid prevents spurious oscillations that may make the solution of the flow field in the NICFD regime more difficult and less robust.Non-ideal gas effects are taken into account by adopting the polytropic Peng-Robinson thermodynamic model. The numerical results focus on the problem of a piston moving in a tube filled with siloxane $\mathrm{MD_4M}$, a simple configuration which can be the core of experimental research activities aiming at investigating the thermodynamic behavior of NICFD flows. Several numerical tests involving different piston movements and initial states in 2D and 3D assess the capability of the proposed adaption technique to correctly capture compression and expansion waves, as well as the generation and propagation of shock waves, in the NICFD and in the non-classical regime.

Tables

Table 1. Table 1: Test conditions for the oblique shock test, and thermodynamic properties of the working fluid MDM MDM \mathrm{MDM} . The deviation angle due by the wedge is θ = 20 ∘ 𝜃 superscript 20 \theta=20^{\circ} and the oblique shock angle is β = 37.60 ∘ 𝛽 superscript 37.60 \beta=37.60^{\circ}   Gori2017 .

$T_c\left(K\right)$	564.1	$P_c\left(\mathrm{bar}\right)$	14.15
$v_c\left(m^3/\mathrm{kg}\right)$	0.00525	$c_{v_{\infty }}/R$	57.69
	Up-stream		Down-stream
$P/P_c$	0.704		1.160
$v/v_c$	3.000		1.236
$T/T_c$	1.037		1.057
$M$	2.000		2.873
$\Gamma$	0.667		0.279

Table 2. Table 2: Thermodynamic properties for MD 4 ​ M subscript MD 4 M \mathrm{MD_{4}M} . The molecular mass M m subscript 𝑀 m M_{\mathrm{m}} , the boiling temperature T b subscript 𝑇 b T_{\mathrm{b}} , the critical temperature and pressure are those included in the NIST software REFPROP  Refprop . The critical volume is given by the 12-parameter EOS  Colonna2006a implemented in the same software. The isobaric specific heat capacity in the ideal gas limit c p ∞ subscript 𝑐 subscript 𝑝 c_{p_{\infty}} is computed through the polynomial expression given in  Colonna2006a . The acentric factor is approximated by Eq. ( 3 )

$M_m\left(\mathrm{kg}/\mathrm{mol}\right)$	0.45899		$T_b\left(K\right)$	533.9
$T_c\left(K\right)$	653.20		$P_c\left(\mathrm{bar}\right)$	8.7747
$v_c\left(m^3/\mathrm{kmol}\right)$	1.7309		$Z_c$	0.2797
$c_{p_{\infty }}/R$	115.99		$\omega$	$0.7981$

Table 3. Table 3: Initial thermodynamic state for the three tests.

	Test $𝒟$	Test $𝒩ℐ$	Test $𝒩𝒞$
Regime	dilute	non-ideal	non-classical
$T_0/T_c$	1.015	1.015	0.9955
$P_0/P_c$	0.2	0.9	0.9
$v_0/v_c$	16.780	2.374	1.918
$Z_0$	0.9257	0.5886	0.4861
${\Gamma }_0$	0.9306	0.4516	-0.0064
$c_0\left(m/s\right)$	101.9	61.9	45.2

Equations

$$\Gamma =1+\frac{\rho }{c}{\left(\frac{\partial c}{\partial \rho }\right)}_s,$$

$$P(T,v)=\frac{RT}{v-b}-\frac{a{\alpha }_{\omega }(T_{\mathrm{r}}{)}^2}{v(v+b)+b(v-b)},$$

$$a=0.45724\frac{R^2{T}_{\mathrm{c}}^2}{P_{\mathrm{c}}},b=0.0778\frac{R{T}_{\mathrm{c}}}{P_{\mathrm{c}}},$$

$${\alpha }_{\omega }(T_{\mathrm{r}})=1+f_{\omega }\left(1-\sqrt{T_{\mathrm{r}}}\right)$$

$$f_{\omega }=0.37464+1.54226\omega -0.26699{\omega }^2.$$

$$\omega =\frac{3}{7}\left(\frac{T_{\mathrm{b}}}{T_{\mathrm{c}}-T_{\mathrm{b}}}\right){\log }_{10}\left(\frac{P_{\mathrm{c}}}{101325}\right),$$

$$\left(\frac{\partial e(T,v)}{\partial v}\right)=T{\left(\frac{\partial P(T,v)}{\partial T}\right)}_v-P(t,v).$$

$$\varphi (T)={\int }_{T_{\mathrm{ref}}}^T{c}_{v_{\infty }}(\tau )\mathrm{d}\tau ,$$

$$e(T,v)=e_{\mathrm{ref}}+\phi(T)-\dfrac{a}{b}\dfrac{\alpha_{\omega}(T_{\mathrm{r}})}{2\sqrt{2}}\,\Bigl{[}\alpha_{\omega}(T_{\mathrm{r}})\\ +f_{\omega}\sqrt{T_{\mathrm{r}}}\,\Bigr{]}\,\ln\dfrac{1+\sqrt{2}+v/b}{1-\sqrt{2}+v/b}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\mathcal{C}}\mathsf{u}\mathrm{d}\mathbf{x}+{\oint }_{\partial \mathcal{C}}\mathbf{f}(\mathsf{u})\cdot \mathbf{n}d\mathbf{s}=0,$$

$$\mathbf{f}(\mathsf{u})={\left[\mathbf{m},\mathbf{m}\otimes \mathbf{m}/\rho +\Pi (\mathsf{u}){\mathbb{I}}^k,\left[E^{\mathrm{t}}+\Pi (\mathsf{u})\right]\mathbf{m}/\rho \right]}^{\mathrm{T}},$$

$$P=P(e(T,v),v)=P(E,\rho )=P\left(E^{\mathrm{t}}-\frac{\Vert \mathbf{m}{\Vert }^2}{2\rho },\rho \right)=\Pi (\mathsf{u}),$$

$$\frac{{\mathsf{u}}_i^{n+1}-{\mathsf{u}}_i^n}{\Delta t}{V}_i=\sum _{k\in {\mathcal{K}}_{i,\ne }}\Phi ({\mathsf{u}}_i,{\mathsf{u}}_k,\text{contour}[2]black{\eta }_{ik}{)}^{n+1}+{\Phi }^{\partial }({\mathsf{u}}_i,\text{contour}[2]black{\xi }_i{)}^{n+1},$$

$$\text{contour}[2]black{\eta }_{ik}={\int }_{\partial {\mathcal{C}}_{ik}}{\mathbf{n}}_i\mathrm{d}\mathbf{x}\mathrm{and}\text{contour}[2]black{\xi }_i={\int }_{\partial {\mathcal{C}}_{i,\partial }}{\mathbf{n}}_i\mathrm{d}\mathbf{s}$$

$$\Phi =\frac{\mathbf{f}({\mathsf{u}}_i)+\mathbf{f}({\mathsf{u}}_k)}{2}\cdot \text{contour}[2]black{\eta }_{ik}+\frac{1}{2}\mathsf{R}|\tilde{\mathsf{\Lambda }}|\left(\Psi -\mathbb{I}\right)\mathsf{L}({\mathsf{u}}_k-{\mathsf{u}}_i)$$

$$\left(\dfrac{\partial P}{\partial E}\right)_{\!\!\!\!\rho}\bigl{(}E_{i}-E_{k}\bigr{)}+\left(\dfrac{\partial P}{\partial\rho}\right)_{\!\!\!\!E}\bigl{(}\rho_{i}-\rho_{k}\bigr{)}=P_{i}-P_{k}\,.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\mathcal{C}(t)}\mathsf{u}\mathrm{d}\mathbf{x}+{\oint }_{\partial \mathcal{C}(t)}\left[\mathbf{f}(\mathsf{u})-\mathsf{u}\mathbf{v}\right]\cdot \mathbf{n}\mathrm{d}\mathbf{s}=0,$$

$${\nu }_{ik}={\int }_{\partial {\mathcal{C}}_{ik}}\mathbf{v}\cdot {\mathbf{n}}_i\mathrm{d}\mathbf{x}\mathrm{and}{\nu }_i={\int }_{\partial {\mathcal{C}}_{i,\partial }}\mathbf{v}\cdot {\mathbf{n}}_i\mathrm{d}\mathbf{s},$$

$$\left({V_{i}^{n+1}\mathsf{u}_{i}^{n+1}-V_{i}^{n}\mathsf{u}_{i}^{n}}\right)\big{/}{\Delta t}=\\ \sum_{k\in\mathcal{K}_{i,\neq}}\!\!\phi(\mathsf{u}_{i},\mathsf{u}_{k},\contour[2]{black}{$\eta$}_{ik},\nu_{ik})^{n+1}+\phi^{\partial}(\mathsf{u}_{i},\contour[2]{black}{$\xi$}_{i},\nu_{i})^{n+1}\,.$$

$$\Delta V_{ik}^{n+1}=\Delta t{\nu }_{ik}^{n+1}\mathrm{and}\Delta V_{i,\partial }^{n+1}=\Delta t{\nu }_i^{n+1},$$

$$\tilde{\Lambda }=\mathsf{diag}\left\{\min \left(\max \left(\frac{c|{\lambda }_p|}{ϵ},\frac{1}{{\ell }_{\max }^2}\right),\frac{1}{{\ell }_{\min }^2}\right)\right\},$$

$$\Vert \mathbf{w}{\Vert }_{\mathcal{M}}=\sqrt{\mathbf{w}{}^{\mathrm{T}}\mathcal{M}\mathbf{w}}.$$

$$\mathrm{d}(u+c)=\frac{\Gamma }{\rho c}\mathrm{d}P.$$

Full text

††journal: Shock Waves

1 Introduction

Computational Fluid Dynamics (CFD) has undoubtedly become an important prediction, analysis, and design tool in many engineering applications. New methods are continuously developed to better exploit computational resources, to increase solution accuracy, and to extend the range of applicability of CFD to diverse fields. For instance, the numerical simulation of fluid flows within the so-called Non-Ideal Compressible Fluid Dynamics (NICFD) regime still represents some challenges for CFD experts, and it is a very active area of research Vitale2015 ; Gori2017 ; Pini2017 ; Ameli2017 ; Gori2017a ; Head2017 ; Keep2017 .

NICFD deals with flows occurring within the thermodynamic region wherein the fluid thermodynamic behavior significantly departs from the one predicted by the ideal gas model, as occurs for dense vapors, supercritical flows, and compressible two-phase flows. In this region, the attractive and repulsive molecular forces are not negligible and non-ideal gas effects, such as non-monotone variations of the Mach number along isentropic expansions, or non-classical phenomena, such as rarefaction shocks, may occur Colonna2006 ; Harinck2009 ; Nannan2014 .

Reliable simulations of NICFD flows require accurate thermodynamic models that are able to include non-ideal effects, but that are usually characterized by complex mathematical descriptions. The inclusion of such thermodynamic models in CFD software is not a trivial task and leads to a considerable increase of the computational time, especially when the flow states are in the neighborhood of the Vapor-Liquid Equilibrium (VLE) curve ColonnaRebay04 .

In recent years, different efforts have been made to extend standard techniques usually adopted under the perfect, i.e., polytropic ideal, gas assumption to the NICFD regime. As a result, some CFD toolkits presently offer the capability to perform simulations of non-ideal compressible flows, such as for instance the open-source multi-physics software suite SU2, which has been recently equipped with a varied thermodynamic library for pure fluids Vitale2015 , and commercial software such as ANSYS CFD or STAR-CCM+.

In addition, new methods are continuously developed to tackle specific numerical tasks or applications. For example, several contributions have been dedicated to numerical schemes Vinokur1990 ; Abgrall1991 ; Arabi2017 , to efficient evaluation of thermodynamic quantities Colonna2003 ; Cinnella2011 , and, more recently, to shock waves Pantano2017 ; Gori2017 ; Passmann2017 , and to turbulent flows From2017 ; Sciacovelli2017 in the NICFD regime. However, some numerical techniques widely used in ideal gas simulations have not been the subjects of an adequate investigation and assessment in the NICFD regime yet. In particular, in our opinion, an example of these shortcomings concerns mesh adaptation techniques in unsteady simulations.

Mesh adaptation is a valuable strategy for the simulation of flow fields characterized by different spatial scales and for unsteady simulations, where the position of the relevant flow features changes in space and/or time or the computational domain undergoes large deformations. In these situations, the grid spacing can be locally reduced or increased to efficiently optimize solution accuracy, while preventing an excessive growth of the number of grid nodes. Concerning the criteria used to modify the grid, different choices can be made, as for instance integral error indicators Dwight2008 ; Fidkowski2011 , local a-posteriori error estimators based on interpolation error analysis Formaggia2003 ; Coupez2011 or local indicators based on the behavior of the flow solution Kallinderis1989 ; Choi2004 ; Kallinderis2017 . Similarly, various adaptation strategies have been proposed, such as node displacement at fixed connectivity (r-adaptation) Huang2011 , local grid-connectivity modification (r-refinement) Freitag1997 ; Mavriplis2000 ; Dapogny2014 , partial or complete re-meshing, or modification of the polynomial order used to discretize the equations on each element (p-refinement) Babuska1994 ; Dolejsi2015 , which can also be combined together, as hp-refinement. However, standard mesh adaptation techniques require the interpolation of the solution from the old to the adapted grid. As is well known, from a numerical point of view, solution interpolation in compressible flow is an undesirable operation, because it may undermine the accuracy of the solution and the properties of the scheme, such as stability and positivity-preserving properties, and it may introduce spurious oscillations. This latter fact can be particularly dangerous in NICFD simulations, where oscillations may bring the thermodynamic state of the fluid under the VLE curve, where the adopted single-phase Equation of State (EoS) is invalid.

To overcome the difficulties related to the solution interpolation, the inviscid flow solver Flowmesh is currently under-development at the Department of Aerospace Science and Technology of Politecnico di Milano flowmeshWeb . The solver implements a finite-volume scheme for unstructured grids based on an innovative interpretation of the modification of grid elements due to mesh adaptation within the Arbitrary Lagrangian-Eulerian (ALE) framework Guardone2011 ; Isola2015 ; Re2017jcp . Thanks to a series of fictitious continuous deformations of the finite volumes, interpolation of the solution from the initial to the adapted grid is avoided, thus preserving the properties of the fixed-connectivity ALE scheme while guaranteeing conservativeness by construction. In our opinion, the absence in the proposed method of an explicit interpolation step has a favorable impact on NICFD simulations, where oscillations in the thermodynamic fluid state may undermine the robustness of the method to the point of jeopardizing the computation of the solution. In addition, the flow solver is linked to the VThermo library, now included in FluidProp FluidProp , which offers the possibility of evaluating thermodynamic quantities by means of different EoSs, such as the van der Waals EoS VanderWaals1873 , the Peng-Robinson EoS PengRobinson1976 , the Martin-Hou EoS MartinHou1955 , the Redlich-Kwong EoS Redlich1949 , and the Soave–Redlich-Kwong EoS Soave1972 . For what concerns mesh adaptation strategies, node displacement, insertion, deletion, and edge swap are exploited to locally modify the grid spacing, with the support of the automatic triangular and tetrahedral re-meshers of the open-source library Mmg mmgWeb .

A natural application of a numerical tool with the outlined features is, in our opinion, the simulation of a flow inside a tube closed by movable pistons. Different types of piston motions can be enforced to generate various flow fields, which may include shock waves, contact discontinuities, and smooth regions. Indeed, this setup is often the core of experimental test-rig arrangements thanks to its feasibility and flexibility Bryson1965 ; Stalker1966 ; Kewley1974 ; Hannemann2016 . However, the numerical simulation of these flow fields is challenging due to the presence of traveling waves with different spatial scales, and it represents a valuable test to assess the capabilities of the proposed adaptive method.

In essence, this paper describes the assessment of an adaptive Euler solver over 2D and 3D unstructured grids for moving body problems in the NICFD regime. To authors’ knowledge, this is the first application of unsteady mesh adaptation strategies to NICFD problems, and the interpolation-free approach here exploited represents an added value with respect to standard adaptation techniques in this peculiar thermodynamic regime. This work draws on the investigation of mesh adaptation criteria based on flow variables in the NICFD regime for bi-dimensional steady simulations described in Re2017jet .

The paper is organized as follows. Section 2 outlines the main features of the thermodynamic modeling in the NICFD regime and the polytropic Peng-Robinson model used in the simulation of the piston problems. Section 3 sketches the numerical method to solve the unsteady Euler equations over adaptive grids within the ALE framework. Then, the adopted mesh adaptation techniques are described in Sec. 4, along with a summary of the entire computational process. An assessment of the proposed unsteady method in 2D and 3D steady tests is given in Sec. 5. Section 6 describes the results of different piston problems. More specifically, Subsec. 6.1 and 6.2 illustrate an oscillating piston in an infinite-length and in a closed tube in the NICFD regime, starting from two different initial conditions. Then, we show the results for an impulsively started piston in a closed tube. Subsection 6.4 displays the results of two tests about the piston experiencing harmonic motion in the non-classical region. Three-dimensional results are presented in Subsec. 6.5. Finally, the conclusions of the present work are drawn in Sec. 7.

2 Thermodynamic modeling in NICFD regime

In NICFD, a quantitative measure of the deviation of the flow behavior from the ideal gas model is given by the value of the fundamental derivative of gasdynamics $\Gamma$ Thompson1971 ,

[TABLE]

which expresses the isentropic variation of the sound speed, being $\rho$, $c$, $s$ the density, the sound speed and the specific entropy, respectively.

The occurrence of non-ideal gas effects is determined by values of $\Gamma$ less than unity. More specifically, the gasdynamic behavior of fluids characterized by $0<\Gamma<1$ shows qualitative differences with respect to the ideal one. The most peculiar is the non-monotone dependence of the Mach number on the density along isentropic expansions because of the increase of the speed of sound Cramer1991 ; Kluwick1993 . For thermodynamic states at $\Gamma<0$, non-classical gasdynamic phenomena, such as rarefaction shock waves and composite waves, may occur Cramer1984 ; Menikoff1989 , although no experimental evidence of the existence of such exotic flow features in single-phase vapor flows is available. The value of the fundamental derivative varies over the thermodynamic plane and some fluids have a finite region near the VLE curve where $\Gamma<1$ Harinck2009 . However, recent advancements have predicted the possible occurrence of even non-classical phenomena near the critical point also for quite simple molecules Nannan2014 ; Nannan2016 .

Given the inappropriateness for NICFD simulations of the perfect, or so-called Polytropic Ideal Gas (PIG) model—which predicts a constant $\Gamma^{\mathsf{PIG}}=(\gamma+1)/2>1$, where $\gamma=c_{p}/c_{v}>1$ is the constant ratio between the isobaric ($c_{p}$) and isochoric ($c_{v}$) specific heat capacities—a non-ideal thermodynamic model is required to describe the flow behavior. For pure fluids, the complete evaluation of the thermodynamic state can be obtained from a thermodynamic potential expressing directly the fundamental relation, as the Helmholtz free energy, or from two compatible equations of states GuardoneArgrow2005 .

Different modeling approaches have been proposed to describe the thermodynamic behavior of the flow in CFD simulations. For instance, the pressure $P$ can be expressed in terms of temperature $T$ and specific volume $v$ by a cubic equation, generally composed by two terms, one accounting for the attractive intermolecular forces and one for the repulsive covolume effects. The simplest cubic EoS is the van der Waals EoS, but several different pressure EoSs have been proposed PengRobinson1976 ; MartinHou1955 ; Redlich1949 ; Soave1972 . To have a complete thermodynamic description, the pressure EoS has to be complemented by a compatible caloric EoS Callen1985 . A different approach may rely on an accurate multi-parameter EoS expressing the fundamental relation, see for instance Span2003 ; Lemmon2006 . Because of their high non-linearity and complexity, a direct implementation in CFD codes of such EoSs is practically unaffordable, but they can be used in a preceding step to build accurate Look-Up Tables which are then recalled during the simulation to evaluate the thermodynamic state by interpolation Rinaldi2014 ; Pini2015 ; Moraga2017 .

2.1 Peng-Robinson equation of state model

In this work, to model the thermodynamic behavior of the flow we use the pressure Peng-Robinson EoS PengRobinson1976 :

[TABLE]

where $a$ and $b$ are two fluid-specific constants given in terms of the critical pressure $P_{\mathrm{c}}$ and temperature $T_{\mathrm{c}}$,

[TABLE]

and $\alpha_{\omega}$ is a dimensionless function of reduced temperature $T_{\mathrm{r}}=T/T_{\mathrm{c}}$, namely

[TABLE]

where $f_{\omega}$ is a constant characteristic of the substance, defined in terms of the acentric factor $\omega$,

[TABLE]

The acentric factor is a parameter that takes into account the polarity and the lack of spherical shapes of the molecules. In this work, its (approximate) value is computed by means of the Edmister’s equation Poling2001 :

[TABLE]

where $T_{\mathrm{b}}$ is the temperature at the normal boiling point.

Given the pressure EoS, the compatible caloric EoS for the internal energy per unit of mass $e=e(T,v)$ can be derived by integrating the reciprocity relation

[TABLE]

However, the integration with respect to the specific volume allows to determine the compatible EoS only up to an integration constant function of temperature, which is usually connected to the dependence of the energy on temperature in the dilute gas limit. In common practice, it is defined as

[TABLE]

where $c_{v_{\infty}}$ is the isochoric specific heat capacity in the ideal gas limit, that is for $v\to\infty$. For more details, see Ref. GuardoneArgrow2005 .

The caloric EoS compatible to Eq. (2) is

[TABLE]

where $e_{\mathrm{ref}}$ is an arbitrary reference value for the specific energy, used in the integration of Eq. (4). In this work, we assume a polytropic behavior in the dilute gas limit ($v\to\infty$), therefore $\phi(T)$ is a linear function of the temperature, i.e., $c_{v_{\infty}}$ is constant.

A convenient parameter that easily expresses the deviation of the pressure EoS from the ideal gas model is the compressibility factor $Z=Pv/(RT)$, where $R$ is the gas constant. For the PIG model, $Z$ is constant and equal to 1. The Peng-Robinson EoS predicts the same compressibility factor at the critical point for all substances. Indeed, by evaluating Eq. (2) at $P_{\mathrm{c}}=P(T_{\mathrm{c}},v_{\mathrm{c}})$, a cubic equation for $Z_{\mathrm{c}}$ is obtained, which has only one real root $Z_{\mathrm{c}}^{\mathsf{PR}}=0.3214$. As a consequence, this EoS, as all the cubic ones, is not able to correctly compute all three critical point coordinates $P_{\mathrm{c}}$, $T_{\mathrm{c}}$, and $v_{\mathrm{c}}$ GuardoneArgrow2005 .

3 Numerical method

The governing equations for unsteady inviscid compressible flows are provided by the Euler equations. For a control volume $\mathcal{C}$ within the spatial domain $\Omega$ they read

[TABLE]

where $\mathbf{x}\in\mathbb{R}^{k}$ and $t>0$ are as usual the position and the time, $\mathsf{u}=[\rho,\,\mathbf{m},\,E^{\mathrm{t}}]^{\mathrm{T}}$ is the vector of conservative variables (density, momentum density, and total energy density, respectively) and $\mathbf{n}(\mathbf{x},t)$ denotes the outward unit vector normal to the boundary $\partial\mathcal{C}$, over which the position vector is denoted $\mathbf{s}\in\mathbb{R}^{k-1}$. The inviscid flux function $\mathbf{f}(\mathsf{u})$ is defined as

[TABLE]

where $\Pi(\mathsf{u})$ is the pressure function and $\mathbb{I}^{k}$ is the $k\times k$ identity matrix. Suitable initial and boundary conditions have to be specified to complement Eq. (7) LeVeque1992 .

The thermodynamic model completes the Euler equations through the pressure function $\Pi$, which expresses the pressure $P$ as a function of the conservative variables. Indeed, the pressure and caloric EoSs are manipulated so that

[TABLE]

where $E$ is the internal energy density (per unit of volume).

3.1 Spatial and temporal discretization

The discrete form of the Euler equations is now briefly described. Since the focus here is on the influence of the thermodynamic model on the method, only the relevant results of the discretization process are illustrated. A detailed description of the scheme can be found in Guardone2011 ; Isola2015 for the 2D and in Re2017jcp ; RePhd2016 for the 3D cases.

A node-centered edge-based finite-volume scheme is exploited to spatially discretize the governing equations over the computational grid. Thus, the domain is split in a finite number of non-overlapping finite volumes $\mathcal{C}_{i}$ so that $\bigcup_{i}\mathcal{C}_{i}(t)=\Omega(t)$, where $\mathcal{C}_{i}$ is a volume surrounding the grid node $i$. By adopting the backward Euler scheme for time integration, the discretized governing equations can be written as

[TABLE]

where $V_{i}$ is the volume of the cell $\mathcal{C}_{i}$, $\mathsf{u}_{i}$ is the average solution $\mathsf{u}$ over $\mathcal{C}_{i}$, $\mathcal{K}_{i,\neq}$ denotes the set of finite volumes different from $\mathcal{C}_{i}$ that share a portion of their boundary with $\mathcal{C}_{i}$; $\Phi$ and $\Phi^{\partial}$ are suitable integrated numerical fluxes that represent, respectively, the flux across the domain cell interface $\partial\mathcal{C}_{ik}=\partial\mathcal{C}_{i}\cap\partial\mathcal{C}_{k}$ and, if the node $i$ lies on the boundary, across the boundary interface $\partial\mathcal{C}_{i,\partial}=\partial\mathcal{C}_{i}\cap\partial\Omega$. The superscripts $n$ and $n+1$ indicate variables evaluated at time $t^{n}$ and $t^{n+1}$, respectively, while $\Delta t$ is the time step in-between. Moreover, we have introduced the integrated normal vectors:

[TABLE]

with $\mathbf{n}_{i}$ denoting the outward normal with respect to the finite volume $\mathcal{C}_{i}$.

A pseudo-time step method is used to solve the non-linear Eq. (8) and the solution at time $t^{n+1}$ is obtained through an iterative process. The pseudo-time derivative is discretized by means of the backward Euler scheme too, and the unsteady residual, containing the numerical fluxes at the actual pseudo-time step, is approximated by a Taylor expansion. Moreover, according to the defect-correction approach, the exact Jacobian of the integrated fluxes is approximated by the Jacobian of the first order fluxes only, to increase its diagonal dominance Koren1988 . At each pseudo-time step, a Symmetric Gauss Seidel method is used to solve the system of linear equations. Further details about the time integration can be found in IsolaPhd ; Carpentieri2009 .

3.1.1 Integrated numerical fluxes

For the integrated flux $\Phi$, a high-resolution scheme based on the Total Variation Diminishing (TVD) approach Harten1997 is used in the present work to correctly capture shock waves and contact discontinuities within the NICFD regime of interest. More specifically, near discontinuities the second-order centered approximation of the interface flux is replaced by the first order monotonicity-preserving Roe scheme. The switch is controlled by the van Leer limiter. By denoting the limiter $\Psi=\mathrm{diag}\left\{\Psi_{1}\dotsc\Psi_{k+2}\right\}$, the integrated numerical flux across the interface $\partial\mathcal{C}_{ik}$ reads

[TABLE]

where $\widetilde{\mathsf{\Lambda}}$ is the Roe matrix, defined here as the Jacobian of the flux function $\mathbf{f}(\mathsf{u})$ projected along the normal direction \hat{\contour[2]{black}{\eta}}_{ik} and evaluated at the intermediate Roe state $\tilde{\mathsf{u}}=\tilde{\mathsf{u}}(\mathsf{u}_{i},\mathsf{u}_{k})$.

Under the polytropic ideal gas assumption, only two variables are required to define the Roe matrix Roe1981 . Conversely, if a non-ideal thermodynamic model is assumed, the Roe matrix is not uniquely determined by the conservation property Guardone2002 and several approaches were proposed to compute the Roe matrix. Most standard approaches rely on an augmented intermediate state, which includes the pressure derivatives as additional variables and aims at obtaining a quasi-Jacobian form of the Roe matrix Glaister1988 ; Vinokur1990 ; Abgrall1991 ; Cox1994 . However, the average thermodynamic derivatives do not retain their proper significance and this may lead to inconsistencies whenever they are used to compute other thermodynamic quantities, such as the speed of sound Toumi1992 ; Guardone2002 . A different approach is proposed in Guardone2002 , where the Jacobian form of the Roe matrix is enforced. This choice results in an intermediate state that is a one-parameter family of solutions, and the following additional condition is imposed to determine the intermediate density:

[TABLE]

Unfortunately, when using complex EoSs the computation of the intermediate density from the previous relation is not easy and it requires numerical techniques Guardone2007a .

Nevertheless, it has been shown that no relevant differences can be appreciated in the results when using different, simpler approaches Mottura1997 ; Guardone2007a . Therefore, the generalized Roe matrix for non-ideal gas flows is computed here following the simplified approach proposed by Cinnella Cinnella2006 , which consists in selecting the intermediate density as $\widetilde{\rho}=\sqrt{\rho_{i}\,\rho_{k}}$, and the intermediate velocity and total enthalpy through the same definitions used for the PIG model, which are valid also for non-ideal gases as proved in Guardone2002 .

3.1.2 Boundary fluxes

The boundary conditions are imposed in a weak form, i.e., by evaluating the flux in a suitable boundary state $\mathsf{u}^{\partial}_{i}=\mathsf{u}^{\partial}(\mathsf{u}_{i},\mathbf{b})$, function of the solution on the boundary node $\mathsf{u}_{i}$ and of the boundary data $\mathbf{b}$.

In the present work, three types of boundary conditions are imposed: the slip wall, the normal inlet/outlet, and the non-reflecting free surface. The former one is imposed by setting the component of the fluid velocity normal to the boundary equal to the boundary velocity along the normal direction. The second condition models an infinite-length open-end tube, that does not influence the flow inside the domain but preserves only the flow direction parallel to the tube axis, without imposing the sense (inflow or outflow). Hence, this boundary condition simply nullifies the component of the fluid velocity tangential to the open-boundary (or, if need be, sets it equal to the tangential velocity of the boundary). For the latter case, the boundary state is computed via characteristic reconstruction, and the correct number of the physical variables that can be imposed is automatically computed by the eigenvalues analysis Selmin1993 ; Vitale2015 .

Because an implicit pseudo-time method is adopted in this work and a Taylor expansion is used to approximate the residual terms containing the numerical fluxes, the evaluation of the Jacobian matrix at the boundary state and the derivative $\frac{\partial\mathsf{u}^{\partial}_{i}}{\partial\mathsf{u}_{i}}$ are also required IsolaPhd .

3.2 ALE framework for dynamic grids

Unsteady problems often require updating the computational domain to follow the boundary motion or deformation. In such situations, it is common practice to redistribute the displacement occurring at the boundary among the internal grid nodes depending on their shape and volume without modifying the grid connectivity, see for instance Batina1990 ; Venkatakrishnan1996 ; Degand2002 . The same technique can be used to move grid nodes to better describe the local flow features, namely for r-adaptation. As a consequence, on a dynamic grid, the formulation of the flow equations has to be modified to account for the relative motion of the grid with respect to the fluid. The Arbitrary Lagrangian-Eulerian (ALE) formulation is a widely-used and effective strategy to accomplish this task Hirt1974 ; Donea1982 . Within the ALE framework, the governing equations are enforced over control volumes that can move and deform independently of the fluid velocity.

The ALE formulation of the Euler equations reads

[TABLE]

where $\mathbf{v}$ is the velocity of the control volume. The previous equation can be discretized through the same strategy outlined in Sect. 3.1, provided that the contribution of the grid velocity is taken into account in the integrated fluxes. In this regard, two additional quantities labeled as interface velocities are defined as

[TABLE]

and the ALE formulation of Eq. (8) can be written as

[TABLE]

The additional constraint known as the Geometrical Conservation Law (GCL) is beneficial to avoid spurious oscillations and instabilities Formaggia2004 ; Mavriplis2006 ; Etienne2009 . The GCL states that the movement of the computational grid should not affect a uniform solution and it amounts to a suitable computation of the geometrical quantities involved in the grid movement Lesoinne1996 . In this work, the interface velocities defined in Eq. (11) are computed in a GCL-compliant fashion by means of the following relations:

[TABLE]

where $\Delta V_{ik}^{n+1}$ and $\Delta V_{i,\partial}^{n+1}$ are the volumes swept during the time step $\Delta t$ by the interfaces $\partial\mathcal{C}_{ik}$ and $\partial\mathcal{C}_{i,\partial}$, respectively. Assuming that the positions of the grid nodes at the beginning and at the end of the time step are given, the swept volumes can be easily computed by exploiting geometrical relations, as shown in Guardone2011 for triangular grids and in Re2017jcp for tetrahedral ones.

3.2.1 Extension to variable connectivity grids

The ALE scheme presented in the previous subsection has been extended also to variable connectivity grids. Previous works Guardone2011 ; Isola2015 ; Re2017jcp detail how to compute the swept volumes $\Delta V_{ik}^{n+1}$ and $\Delta V_{i,\partial}^{n+1}$ when the local connectivity varies due to mesh adaptation. Connectivity changes are described as fictitious continuous deformations of the involved finite volumes through a series of collapse and expansion operations. More precisely, the following three steps can be defined with the aid of the fictitious dimensionless time $\tau$ (see also Fig. 1 for an example about edge split in 2D).

Collapse, $0\!<\!\tau\!<\!0.5$: all elements involved in the local adaptation collapse into an arbitrary point. 2. 2.

Connectivity change, $\tau\!\!=\!0.5$: when all involved elements reach null volumes, nodes may be inserted or deleted. 3. 3.

Expansion, $0.5\!<\!\tau\!<\!1$: all elements still active (i.e., not deleted at $\tau=0.5$) expand to the final configuration.

The collapse and expansion phases are merely continuous deformations from the initial position to the pivotal point, and from it to the final position. Since the positions of these points are known, the volume swept by the involved interfaces can be easily computed from geometrical considerations. Moreover, the fictitious fluxes generated by these deformations can be expressed within the ALE framework as shown for standard fixed-connectivity deformation and included in Eq. (10). Conversely, the connectivity change at $\tau\!\!=\!0.5$ occurs while the elements have null volume, so no volume is swept by any interface and no flux is exchanged. In this way, we are able to compute the GCL-compliant interface velocities through Eqs. (13) also when local grid topology is modified and to include the volume changes due to mesh adaptation into the standard ALE scheme, without undermining the fixed-connectivity properties and enforcing conservativeness.

The three steps and the core of the procedure are the same for node insertion, node deletion, or edge swap, both in 2D and 3D; only the involved elements depend on each specific grid modification. The three-step procedure is carried out immediately after each connectivity modification and it is instantaneous, i.e., it starts and ends during the same time step $\Delta t$ at which the grid is adapted. Finally, we remark on the local character of the three-step procedure. The interfaces that are not modified by the connectivity change do not take part in three-step procedure, as their contribution would be null since they would sweep the same volume but with different sign during the collapse and the expansion phase Re2017jcp .

4 Mesh adaptation in unsteady NICFD simulations

The accuracy of the numerical solution of the flow equations is strongly influenced by the grid spacing. For simple geometries or in steady simulations, it is possible to generate a priori a mesh that provides sufficiently small numerical errors with respect to the modeling errors. This task becomes considerably more challenging or impossible for complex flow fields characterized by different spatial scales and in unsteady simulations, where the grid spacing has to be related to the behavior of the solution in an efficient way. Solution-dependent adaptive grid techniques represent a valuable tool to deal with this requirement. More specifically, in unsteady simulations, mesh adaptation techniques may be profitably used to tackle two different tasks: to handle large deformations of the boundaries, and to capture the relevant flow features that originate and move through the domain Johnson1999 ; Hassan2000 ; Venkatakrishnan1996 , as outlined in the following subsections.

4.1 Adaptation on flow features

The solution accuracy is here optimized following the equi-distribution principle, whereby the grid spacing is modified to equi-distribute the error over the mesh Borouchaki1997 ; Dolejsi1998 . The first step of this strategy consists in the definition of an error estimator, based upon the flow solution. It appears evident that an effective error indicator will be extremely advantageous to the entire adaptation process. In the NICFD regime, the error estimators based on the Mach number seem to be more efficient than the density-based ones because of the non-monotone dependence of the Mach number on density along isentropic expansions for supersonic flows at $\Gamma<1$ Re2017jet .

A metric map is used to prescribe the size of grid elements during the mesh modification process. If anisotropic mesh adaptation is performed, the shape and the orientation of the elements are also prescribed Borouchaki1997 ; DelPino2011 ; Coupez2011 . Two different procedures are followed to build isotropic or anisotropic metrics. In the former case, a target grid spacing is computed at each grid node by reducing (or augmenting) the actual average size of the edges connected to the node if the error is greater (or less) than a refinement (coarsening) threshold, prescribed in terms of mean and standard deviation of the estimated error over the domain. Conversely, the anisotropic map is built on the basis of the eigenstructure of the Hessian matrix of a certain solution variable FreyAlauzet2005 , as $\mathcal{M}={R}\widetilde{\Lambda}L,$ where $R$ and $L$ are the right and left eigenvectors and

[TABLE]

with $\lambda_{p}$ the $p$-th eigenvalue, $\ell_{\mathrm{max}}$/$\ell_{\mathrm{min}}$ the maximum/ minimum edge length, $c=\frac{k^{2}}{2(k+1)^{2}}$ a constant, and $\epsilon$ a user-defined threshold for the maximum acceptable error, which in all presented computations ranges from $10^{-7}$ to $10^{-6}$. More details about the construction of the metric maps can be found in RePhd2016 ; Re2017aiaa ; Re2017coupled .

The metric $\mathcal{M}(\mathbf{x})$ is a field of symmetric positive matrices $\mathbb{R}^{k\times k}$ that defines a Riemannian structure over $\Omega$. The length of a vector $\mathbf{w}$ in terms of this map is given by

[TABLE]

and the goal of mesh adaptation is to obtain a unit mesh with respect to this metric, i.e., a mesh such that $\lVert\mathbf{w}\lVert_{\mathcal{M}}=1$ for all edges.

Given the initial grid and the metric $\mathcal{M}$, the open-source library Mmg mmgWeb ; Mmg2008 performs automatically a series of different local modifications to make the initial grid as similar as possible to a unit mesh Dapogny2014 . At first, the length of all grid edges is computed according to Eq. (14), then the grid is locally modified both on the interior and on the boundary. Where the edges are too long, the grid spacing is reduced by inserting a new node via element or edge split, or Delaunay triangulation. On the contrary, the nodes connected to edges having $\lVert\mathbf{w}\lVert_{\mathcal{M}}\gg 1$ are removed via edge collapse. Edge swapping and barycentric regularization are also used. In addition, the gradation control technique Borouchaki1998 is exploited to limit the variation in size among adjacent grid elements. All these operations are performed only if they lead to an improvement of the quality of the involved elements, which is defined as

• •

for tetrahedra: \;Q_{m}=\alpha\,V_{m}/\big{(}\,\textstyle{\sum}_{i=1}^{6}\;\ell_{i}^{2}\,\big{)}^{3/2},

• •

for triangles: \;Q_{m}=\ell_{\max}\;\sum_{i=1}^{3}\;\ell_{i}\big{/}(2\;V_{m}),

where $V_{m}$ is the volume (or area in 2D) of the element $m$, $\ell_{i}$ is the edge length, $\ell_{\max}$ is the longest edge in the triangle, and $\alpha$ is a constant introduced to give $Q_{m}=1$ for a regular tetrahedron. It is not possible to specify a maximum or minimum number of grid nodes, but the effect of mesh adaptation can be controlled through the parameter $\epsilon$ in the anisotropic case, and through the refinement/coarsening thresholds in the isotropic case.

4.2 Large boundary displacements

Mesh adaptation here is exploited also to cope with large boundary movements. At each time step, after the boundary nodes are moved to the new positions, the internal nodes are redistributed keeping the grid connectivity fixed. According to the elastic analogy Batina1990 , the largest grid elements account for the major part of the deformation while the smallest ones move almost rigidly. However, when the displacements are large, this strategy may lead to badly-shaped or null-volume elements. In such situations, mesh adaptation is exploited to restore grid quality.

As detailed in Re2017aiaa , the displacement to be imposed at the boundary nodes during the time step (labeled as $\Delta\mathbf{x}_{\mathsf{B}}^{n+1}$) is parameterized by a linear function. If the elastic analogy fails in producing a valid mesh, the displacement is split into smaller portions and the largest portion of $\Delta\mathbf{x}_{\mathsf{B}}^{n+1}$ that leads to a valid mesh is imposed. Then, mesh optimization is carried out aiming only at increasing the quality $Q_{m}$ of all elements, not at modifying the grid spacing. In this specific application of mesh adaptation, the edge swapping technique is particularly effective to restore the element quality, while keeping the number of grid nodes constant Freitag1997 .

4.3 Summary of the computational procedure

To summarize, the main steps of the computational procedure adopted for unsteady adaptive simulations are outlined, with reference to Fig. 2. The generic time step $t^{n}\leq t\leq t^{n+1}$ begins with the grid $\mathcal{K}^{n}$ and the corresponding solution $\mathsf{u}(t^{n},\mathcal{K}^{n})$. Then, the following operations are performed.

Mesh deformation:

the movement $\Delta\mathbf{x}_{\mathsf{B}}^{n+1}$ is imposed at the boundary and the internal nodes are consequently relocated. If the elastic analogy fails, a mesh optimization stage including edge swapping is performed (see Sec. 4.2). The resulting grid is labeled $\mathcal{K}^{n+}$ to highlight that it is intermediate, namely it complies with the new geometry but the solution is the previous one.

Prediction:

the solution $\mathsf{u}^{n+}=\mathsf{u}(t^{n+1},\mathcal{K}^{n+})$ is computed over the new grid. This step prevents a delay between the solution-based mesh adaptation and the actual geometry.

Error estimate:

the metric map $\mathcal{M}$ is generated on the basis of the solution $\mathsf{u}^{n+}$ as summarized in Sec. 4.1.

Mesh adaptation:

the grid $\mathcal{K}^{n+}$ and the metric map are used as inputs to the library Mmg, which locally adapts the grid through node insertions, deletions, relocations, and edge swapping.

• Mesh update:*

the flow solver receives from Mmg all performed local modifications $\Delta\mathcal{K}^{n+}$, that are needed to update the finite-volume discretization and to compute the swept volumes $\Delta V_{ik}$, $\Delta V_{i,\partial}$ by means of the three-step procedure (see Sec. 3.2.1).

Solution:

the final solution $\mathsf{u}^{n+1}$ on the grid $\mathcal{K}^{n+1}$ is computed, using the solution $\mathsf{u}^{n+}$ as initial guess.

Steps three to five can be repeated to perform multiple adaptation cycles, as indicated by the loop s in Fig. 2.

The use of the re-mesher library Mmg mmgWeb ; Mmg2008 alleviates the inherent difficulties of mesh adaptation over unstructured three-dimensional grids, which, for instance, requires to be able to deal with a variable number of elements involved in each local grid modification or with the multiple options that are available for an edge swapping. The choice to exploit an external library to cope with these tasks allows to implement mesh adaptation techniques into the flow solver efficiently, although some modifications to the Mmg library have been required to exploit the three-step procedure of Sec. 3.2.1.

5 Assessment of the unsteady adaptive method for in NICFD regime

This section presents the first validation of the proposed unsteady method for 2D and 3D problems in the NICFD regime. For the validation of the steady version, the reader is referred to Re2017jet . The geometries of the test cases consist in a simple symmetric wedge in 2D and in a portion of a cone in 3D, exposed to a free stream. As it is known, a steady oblique (conical, in 3D) shock wave is generated to deflect the flow such that it is aligned with the downstream wall. This steady problem can be investigated by means of an unsteady approach if we consider the far-field quiet and the wedge (cone, in 3D) moving at the free stream velocity, namely, if we switch from the body to the laboratory reference frame. Obviously, the thermodynamic variables that are independent from the reference frame, such as the pressure, are the same in the steady and unsteady simulation.

5.1 Unsteady simulations of an oblique shock

We selected a specific configuration that highlights the potential non-ideal effects in oblique shock waves. As shown in Gori2017 , when $0<\Gamma<1$, the non-monotone variation of the speed of sound along isentropic expansions can result in an increase of the flow Mach number across oblique shock waves, contrary to what is predicted by the perfect gas model. Here, we reproduce one of the tests described in Gori2017 to compare our results with the analytical and numerical solution presented therein. The working fluid is the siloxane $\mathrm{MDM}$111$\mathrm{MDM}$ is the acronym for the octamethyltrisiloxane, whose chemical formula is $\mathrm{C_{8}H_{24}O_{2}Si_{3}}$, and the test conditions are detailed in Table 1. As in the reference, we adopt the polytropic van der Waals model. Further information about the implementation of this thermodynamic model in the flow solver can be found in Re2017jet .

The domain extends from $x=0$ to $x=1.7$ and from $y=-0.6$ to $y=0.6$ (assuming a unitary dimensional reference length of $1~{}\mathrm{m}$). The leading edge of the wedge is at $x=1.2$, $y=0$. Although the domain is symmetrical with respect to the line $y=0$, it is discretized by an unstructured grid, so the symmetry is lost from the numerical point of view. Inflow boundary conditions for density, velocity, and pressure are enforced on the right end-wall, while the far-field and the outflow are modeled by non-reflecting boundary conditions that, where possible, impose only pressure and density. Starting from an initial coarse grid, initialized everywhere with the up-stream conditions, three adaptation strategies are performed by imposing three different minimal edge sizes: $h_{\mathrm{min}}=\left\{0.0008,0.0004,0.0002\right\}$. In all cases, an isotropic adaptation criterion is built by blending together the Hessian of the pressure and the gradient of the Mach number. Five adaptation cycles (Prediction, Mesh adaptation, Solution) are performed, and all strategies successfully delivered a refinement of the shock region, up to the minimum prescribed edge size. The results obtained in the steady simulations with $h_{\mathrm{min}}=0.0002$ are shown in Fig. 3. Remarkably, the post-shock Mach number is larger than the pre-shock one, thus pointing out the non-ideal nature of the shock. A comparison against the analytical solution is given in Fig. 4: an excellent agreement with the analytical solution is achieved.

The steady solutions are then used to initialize the unsteady simulations. The free stream velocity is subtracted from the solution and a displacement of $\Delta x=-1$ at $M=2$ is imposed to the wedge. With respect to the three different minimal edge sizes, the whole simulation time is divided in $N_{\mathrm{T}}=\left\{1250,2500,5000\right\}$ steps in order to have a Courant number $\mathrm{Co}=1$. The top picture in Fig. 4 shows that the initial pressure profile (obtained in the steady simulation) is conserved well during the unsteady computations. As expected, the size $h_{\mathrm{min}}=0.0002$ allows to achieve the best results, especially in terms of the post-shock state. The pressure is slightly over-predicted if one uses a lower resolution grid. The symmetry of the problem is recovered in the numerical results, since no differences between the profiles extracted from the lower and the upper part can be appreciated. Further simulations are performed with the highest grid resolution at $\mathrm{Co}=\left\{2,4\right\}$, to assess the validity of the adaptive unsteady scheme also at higher Courant numbers. The results, plotted in Figs. 4(bottom) and 5, agree well with the analytical ones, and the shock resolution is deteriorated only to a small degree by the high Courant number.

5.2 Unsteady simulations of a conical shock

A test similar to the one described in the previous subsection is now presented to assess the validity of the proposed method in 3D. The geometry of the test can be viewed as the solid of revolution obtained by rotating half of the domain of the wedge test through 90 degrees about the $x$-axis. The starting wedge geometry used for the rotation is smaller than the one used in Subsec. 5.1, and it extends from $x=0$ to $x=0.7$ and from $y=0$ to $y=0.3$. In the resulting three-dimensional domain, the far-field is represented by a quarter of a cylindrical surface, and the solid body by a quarter of a conical one. The lateral planes at $y=0$ and at $z=0$ are modeled as solid walls to respect the symmetry. For the inflow and the outflow, the same non-reflection conditions used in the oblique shock test are enforced.

An adaptive steady simulation with the free-stream conditions given in Tab. 3 is first performed. Even though the upstream conditions are the same, the flow state downstream of the conical shock differs from the one downstream of the oblique shock, as is well known. No analytical solution of the conical shock waves with van der Waals EoS is available.

The steady results are then used to generate the initial solution for the unsteady computation, which starts from a grid made of 150633 nodes and 882581 elements, with a minimum edge size of $h_{\mathrm{min}}=0.0002$. A displacement of $\Delta x=-0.2$ is imposed to the cone and a Courant number of $\mathrm{Co}=2$ is enforced. The results are shown in Fig. 6. We can observe that the flow field is well conserved during the displacement thanks to the fine grid spacing obtained near the shock as a result of mesh adaptation. Figure 7 confirms the good agreement between the unsteady and the steady pressure fields. Moreover, three planes are sliced at different azimuth angles, and the resulting pressure contours are plotted in the $x$-$R$ plane, where $R=\sqrt{y^{2}+z^{2}}$. An almost perfect overlap between contour plots over the three planes is achieved, which confirms the capability of the proposed approach to correctly recover the cylindrical symmetry of the flow field. Mesh adaptation plays a crucial role to achieve this result, because only a proper refinement of the grid regions where the flow variations are significant can overcome the asymmetry of the unstructured grid.

6 Numerical investigation of unsteady piston problems

The adaptive computational procedure described in the previous sections is here further assessed through the simulations of different piston problems, in 2D and 3D. The piston problem is a standard gasdynamics test which presents diverse numerical challenges. It models the motion of a piston inside a tube with no gap between the tube and the piston itself. If the piston moves towards the right, the fluid lying in the right part of the tube is compressed, the compression wave moves away from the piston and travels rightwards. Correspondingly, in the left part of the tube an expansion wave is generated and travels toward the left. While these nonlinear waves propagate through the tube, their waveform is distorted. If $\Gamma>0$, the compressive waves steepen, possibly forming a classical compression shock, and the rarefaction waves spread out. On the contrary, if $\Gamma<0$, the compressive part of the wave spreads over a longer space as it moves, while the expansion part steepens and it may possibly form a rarefaction shock, the most representative flow structure of the non-classical regime. If the tube is closed, the waves are reflected at the tube ends and interact with the incoming waves, originating a complex flow pattern which contains compression and expansion waves of different intensities, shocks, and contact discontinuities. In addition, the intensity and the type of the perturbation waves can be roughly controlled by the piston motion, thus different tests can be performed. A flow field encompassing simultaneously all these flow features represent a valuable, challenging test to assess the proposed method, since the finite-volume solver is edge-based and the numerical simulations are performed over unstructured grids, so they are inherently multidimensional. Rather, this feature provides an additional opportunity to assess the validity of the proposed method.

The geometry of the reference piston problem is shown in Fig. 8. The tube has different lengths depending on performed tests. For 3D simulations, a square section is assumed, i.e., $L_{\mathrm{P}_{z}}=L_{\mathrm{P}_{y}}$. The tube is filled with the linear siloxane $\mathrm{MD_{4}M}$222$\mathrm{MD_{4}M}$ is the acronym for the tetradecamethylhexasiloxane, whose chemical formula is $\mathrm{C_{1}4H_{42}O_{5}Si_{6}}$, which is a high-molecular complexity fluid that presents a quite extended $\Gamma<1$ region and a small one at $\Gamma<0$. Table 2 summarized the relevant thermo-physical properties of $\mathrm{MD_{4}M}$. As expected, the compressibility factor at the critical point is different from the one predicted by the Peng–Robinson EoS. Indeed, the specific volume computed through the adopted cubic thermodynamic model slightly deviates from the value reported in Tab. 2.

Initially, the fluid inside the tube is at rest (labeled with subscript 0). Three different initial states are used to investigate peculiar gasdynamic behavior of the fluid: the initial conditions for the dilute (Test $\mathcal{D}$) and the NICFD (Test $\mathcal{NI}$) regimes are selected along an isotherm above the critical one, while the initial conditions for the non-classical case (Test $\mathcal{NC}$) are selected along an isotherm below the critical one. The initial states in all tests are detailed in Tab. 3 and they are displayed on the $P-v$ plane, along with the VLE curve, in Fig. 9.

In this work, two different types of motion are simulated: a harmonic motion which transfers to the fluid in contact with the piston a sinusoidal perturbation wave, and an impulsive start, which generates instantaneously a shock wave and a rarefaction fan which travel outwards. In the harmonic motion, the position of the piston is set as $x(t)=x_{0}(t)+A\cos(2\pi\,f\,t)$, where $A$ is the amplitude and $f$ the frequency ($f=30~{}\mathrm{Hz}$ for all tests). In the other case, the piston impulsively acquires a constant positive velocity $V_{\mathrm{P}}=|A\pi\,f|$, that is half of the maximum velocity experienced by the piston during the harmonic motion. The piston motion represents an additional challenge for the CFD software, as it requires dealing with large deformations, especially in the simulation of the impulsive start, where the piston displacement amounts to more than the $40\%$ of the length of the domain.

The harmonic motion is first simulated in a tube with infinite-length to assess the capability of reproducing both the expansion and compressive waves generated by the motion and to observe their distortion. Then, the simulation is repeated considering a closed tube, to assess also the capability to correctly model the interactions with solid boundaries and other waves.

In the following subsections, the exposition and discussion of the results starts from the 2D tests in the NICFD region, with the harmonic motions of the piston in a infinite-length tube and in a finite-length tube, followed by the impulsive start in a closed tube. Then, results concerning the non-classical region are shown for the harmonic motion of the piston in an infinite-length tube, with two different amplitudes. Finally, some 3D results are displayed to show the capability of tackling 3D problems. In all figures of this section, a reference length of $1~{}\mathrm{m}$ is used to make all distances dimensionless.

6.1 2D oscillating piston in an infinite-length tube: dilute and non-ideal conditions

The first test case consists in the harmonic motion of the piston in an infinite-length 2D tube, with an amplitude $A=-0.25~{}\mathrm{m}$. In these simulations, the focus is on the compressive and expansive waves generated by the harmonic motion, so we have limited the influence of the boundary on the flow field by modeling the left and right ends as normal inlets/outlets, and by imposing the tube lengths $L_{\mathrm{L}}=L_{\mathrm{R}}=10~{}\mathrm{m}$. The initial grid is composed of $6998$ triangles and $4502$ nodes and three oscillation periods are simulated.

The piston starts moving towards the right. Its velocity progressively increases during the first quarter of the oscillation period, $0<t<0.25T$, until $V_{\mathrm{P},\mathsf{max}}=-2\pi f\,A$. During this interval, a compressive wave is generated on the right side of the piston and moves towards the right end of the tube. In the next quarter of the period $0.25T<t<0.5T$, the piston still moves towards the right but its speed reduces ($V_{\mathrm{P}}=0$ at $t=0.5T$). This gives rise to an expansion on the right section and a compression on the left part, which continue also in the next quarter of period. At $t=0.75T$ the velocity of the piston is minimum, $V_{\mathrm{P},\mathsf{min}}=2\pi f\,A$ (same magnitude of $V_{\mathrm{P},\mathsf{max}}$ but negative), then, it increases leading to an expansion on the left side and a compression on the right side, similarly to the first quarter. The flow fields at different times during the first oscillation period are shown in Fig. 10 for Test $\mathcal{D}$ and in Fig. 11 for Test $\mathcal{NI}$. The solutions at $t=T$ are shown also in the $P$-$v$ plane in Fig. 9, from which the significant variation of $\Gamma$ in Test $\mathcal{NI}$ can be appreciated.

During the next periods, the flow field evolves in a similar way. Nevertheless, unlike the ones occurring in the right part at $0<t<0.25T$, the following compressions are generated by the largest variation of the piston velocity, i.e., from the minimum to the maximum value, or vice versa. Therefore, these waves may coalesce together generating a shock far from the piston side.

Simulations with $N_{T}=\left\{100,200,300\right\}$ steps per period are carried out to estimate the minimum number that does not impair the solution accuracy. Figures 12 and 13 show the pressure profile obtained with at $t=1.4T$ and $t=1.8T$ for Test $\mathcal{D}$ and Test $\mathcal{NI}$, respectively. For both cases, we can observe that the largest time step leads to a loss of accuracy in proximity of the pressure maxima and minima, while no significant deviations can be noted for the other two time steps. Therefore, $N_{T}=200$ steps per period are used in all following simulations. Although the resulting plots appear as single lines, Figs. 12 and 13 are in fact scatters of all points of the unstructured grids. The one-dimensional behavior of the flow is therefore perfectly recovered in both tests, and this gives the first, partial, assessment of the validity of the proposed scheme. Small symmetry disturbances that can be detected in the contour plots of Figs. 10 and 11 are not present in the pressure profiles, and they do not affect the symmetry of the results. In our opinion, they could be ineffective consequences of the unstructured grids, which do not undermine the validity of the proposed approach.

In these tests, an anisotropic metric is built on the basis of the Hessian of the Mach number and two adaptation cycles are performed every time step. More precisely, the estimator is obtained dividing the Hessian by the gradient to reduce the predominance of the strongest phenomena Re2017jet . The beneficial effects of mesh adaptation are clearly displayed by Figs. 10 and 11, which, in addition to the pressure contours, show the histograms of the number of grid points along the $x$-axis, and two grids. As perturbation waves propagate outwards, mesh adaptation affects a larger portion of the domain, as it can be seen from the grids at $t=0T$ and $t=0.2T$ (grids at the next times are not shown since the node density in some regions becomes so high that mesh elements are indistinguishable at this resolution). Exploiting the one-dimensional character of the flow, the grid density at three subsequent times is represented by the histogram of the number of grid nodes in fixed-width bins. This type of graph allows to highlight the effect of mesh adaptation on a broad scale: the grid points are gathered where the flow field exhibits large variations, while the grid density is reduced in correspondence of almost uniform regions. Grid coarsening is of paramount importance in unsteady simulations, since it allows limiting the number of grid nodes. Mesh adaptation allows following flow variations as, for instance, occurs with the compression in the right part of the tube in Fig. 11: the peak in the histogram around $x=1.15$ at $t=0.6T$ is progressively spread and moved towards the left, following the intensity and position of the compression. Moreover, the grid density reflects the strength of the flow variations, as it can be clearly noticed by looking at the compressive waves and the corresponding histogram peaks in the left part of the tube in both tests. Therefore, the proposed approach is capable of linking the grid spacing with the unsteady flow features.

A minor remark about Fig. 10 and 11 concerns the coarsening of the unperturbed regions of the initial grid, resulting, for example, in the different grid spacing for $x<-1.2$ between the grids at $t=0$ and $t=0.2T$ in Fig. 10. According to our experience, an extremely coarse initial grid may jeopardize the development of important flow features, so it is preferable to start with a barely coarse mesh and exploit mesh coarsening already in the first adaptation step, which is carried out after computing the solution in the Prediction step, as explained in Subsec. 4.3. If, as in these tests, the maximum edge size allowed during the adaptation is different from the one of the initial grid, and the flow field at the beginning of the simulation embeds large uniform regions, the first adaptation step profoundly modifies the initial grid. The regions that remain uniform during the next time steps are not affected by the following mesh adaptations.

Analysis of the flow behavior in dilute regime.

To better analyze the flow behavior, the values of different flow quantities are extracted at fixed locations during the whole simulation time. For Test $\mathcal{D}$, these quantities are shown in Fig. 14. The results obtained with the Polytropic Peng–Robinson (PPR) EoS are compared to the ones obtained with the Polytropic Ideal Gas (PIG) model, computed starting from the same $P_{0}$ and $T_{0}$. Although the initial condition of this test lies in the dilute regime where the flow behavior is expected to differ only quantitatively from the ideal gas model, a first, general observation of Fig. 14 reveals that the differences between the two thermodynamic models are significant. The most notable deviations affect the density and speed of sound profiles. We can interpret these deviations as a measure of the non-ideal gas effects. During the simulation, the value of $\Gamma$ oscillates between $0.83$ and $0.96$, but the PIG model is not able to detect it, as it predicts a constant $\Gamma=1.004$. This difference is extremely important, because it is reflected in the speed of sound. Indeed, reminding Eq. (1), a value of $\Gamma$ constant and close to unity results in a constant value of the speed of sound along an isentropic expansion, as shown in the profile for the PIG model. Conversely, for the PPR model, the shape of the profile of $c$ mirrors the one of $\Gamma$.

The differences on $c$ and $\rho$ affect the propagation of the discontinuities originated on the piston surface. For weak discontinuities, the acoustic theory links the pressure and velocity perturbation ($\delta P$ and $\delta u$, respectively) through the initial acoustic impedance $\rho_{0}c_{0}$, namely $\delta P=\rho_{0}c_{0}\delta u$. This relation suggests that, for the piston problem under investigation here where the only perturbation is the piston movement, if density and sound speed are inaccurate, the same velocity perturbation generates a very different pressure perturbation. For instance, defining the deviation on the pressure as $E_{P}=100(P^{\mathsf{PR}}-P^{\mathsf{IG}})/(P^{\mathsf{PR}}-P_{0})$, the deviations on the maximum and minimum values are

at $\bar{x}_{1}=1.0$:

$\quad E_{P,\mathsf{max}}=51\%,\qquad E_{P,\mathsf{min}}=42\%$,

at $\bar{x}_{2}=3.0$:

$\quad E_{P,\mathsf{max}}=45\%,\qquad E_{P,\mathsf{min}}=29\%$,

which show the failure of the PIG model to describe the flow behavior in this region, despite a value of $\Gamma$ near $1$.

A further difference that can be observed in the profiles concerns the shape. In the PIG model, the perturbation is propagated from the piston surface through the tube with an almost constant shape, while according to the PPR model the initial sinusoidal shape loses its symmetry and the compression part becomes steeper. This can be explained by watching the quantity $u+c$, which represents the wave speed in the right part of the piston. For a one-dimensional homentropic flow, the method of characteristics gives

[TABLE]

Thus, for the PIG model the higher value of $\Gamma$ and the lower value of the acoustic impedance generate much faster waves, which take shorter time to reach the same position. Finally, for both models, being $\Gamma>0$ the compressive part of the wave becomes steeper, until a shock may form, generating a total pressure loss. This phenomenon is captured only by the PPR model, which predicts a faster increase of $u+c$. For this reason, the profiles at the location farther from the piston surface ($\bar{x}_{2}=3.0$) display different peak values, although these differences are limited since the shocks are quite weak.

For the sake of completeness, Fig. 14 reports also the Mach number profile, which embeds all the deviation explained above and clearly displays the loss of the initial sinusoidal shape as the maxima and the minima are not equally spaced.

Analysis of the flow behavior in NICFD regime.

The same quantities analyzed for Test $\mathcal{D}$ are shown in Fig. 15 for Test $\mathcal{NI}$, which takes place in the core of the NICFD regime, as clearly displayed by Fig. 9. The variations of the fundamental derivative $\Gamma$ profile are much more significant than the ones in Test $\mathcal{D}$ and the super-critical regime is also reached, in which the behavior of the flow is more similar to the liquid than the vapor one. A peculiar phenomenon of the NICFD regime can be observed looking at the $c$ profile. According to Eq. (1), the variation of the speed of sound during an isentropic expansion is proportional to $\Gamma-1$. So, in the first part of the expansion where $\Gamma>1$, the speed of sound decreases, as occurs in ideal gas, while it increases in the second part, where $\Gamma<1$. Nevertheless, the increase is not sufficiently strong to result in a non-monotone Mach profile. Watching at the other quantities, a behavior similar to the one described for Test $\mathcal{D}$ can be observed, although the shocks are now stronger, because the lower value of the speed of sound (see $c_{0}$ in Tab. 3) leads, given the same piston velocity, to a higher Mach number. For the same reason, also the wave propagation speed is lower.

The results of the first tests and the previous analysis demonstrate that the proposed method is capable of detecting all expected NICFD phenomena, without introducing spurious oscillations due to mesh adaptation, and it is robust.

6.2 2D oscillating piston in a finite-length tube

In the second type of test, presented here, the piston moves inside a finite-length tube, to simulate the reflection of the waves by solid walls and the interaction between waves with opposite directions. The harmonic motion imposed to the piston is the same one described in the previous section. In Test $\mathcal{D}$, the lengths of the tube are $L_{\mathrm{L}}=10~{}\mathrm{m}$ and $L_{\mathrm{R}}=8~{}\mathrm{m}$ and the initial grid is composed of $4056$ nodes and $6306$ triangles. In Test $\mathcal{NI}$, as the wave propagation speed is lower, the tube is shorter, $L_{\mathrm{L}}=8~{}\mathrm{m}$ and $L_{\mathrm{R}}=5~{}\mathrm{m}$, and the initial grid consists in $2927$ nodes and $4548$ elements.

The details of the reflection of the waves at the right wall are shown in Fig. 16 and 17. When a compressive wave reaches the right wall, it is reflected back and starts moving upstream, in the direction opposite to the flow. Thus, it interacts first with the expansion waves, reducing its intensity, then it comes into contact with the incoming compressive waves. The interaction between compressive waves of opposite family originates two compressions with different strength that move outwards from the interaction point.

Mesh adaptation is crucial to correctly detect such a complex, variable flow pattern, where the a priori determination of the locations of all relevant flow features is not practical. The histograms of the grid nodes along the $x$-axis in Fig. 16 and 17 highlight that the grid density is properly modified according to the position and the intensity of the flow variations. In particular, the histograms in Fig. 17 display at times $t=2.25T$ and $t=2.5$ two peaks corresponding to the incoming compressive waves (left) and to the reflected one (right) which move towards each other; as the waves coalesce at $t=2.75T$ only one peak is present; but at $t=3T$ when the compressive waves go in opposite directions, two refined regions are again generated. So, the present test assesses also the capability of the proposed adaptive method to deal with wave interactions and reflections in the NICFD regime, which is an essential feature to simulate more complex and realistic flow fields.

6.3 2D impulsively started piston

The last 2D test case in the NICFD regime is the impulsively started piston in a closed tube, with sizes $L_{\mathrm{L}}=1.5~{}\mathrm{m}$ and $L_{\mathrm{R}}=2~{}\mathrm{m}$. For brevity, only Test $\mathcal{NI}$ is shown. At time [math], the piston is instantaneously set into motion with a constant velocity $V_{\mathrm{P}}=47.12~{}\mathrm{m/s}$, generating a shock wave in the right part of the tube and an expansion fan in the left part. The simulation is performed until the piston displacement is $\Delta x_{\mathrm{P}}=1.51~{}\mathrm{m}$, so that the shock is reflected twice, once by the right end of the tube and once by the piston face. The initial grid is composed of $789$ nodes and $1574$ elements and the whole simulation time is divided in $N_{T}=1010$ time steps. On the smallest allowed edge ($h_{\mathrm{min}}=0.008~{}\mathrm{m}$), the Courant number is $\mathrm{Co}=1.8$ with respect to the piston velocity, but $\mathrm{Co}=3.1$ with respect to the velocity of the first generated shock.

Figure 18 displays the results of the flow field before the reflections of the shock and the expansion fan, and after the shock reflection on the piston surface. In these pictures the effectiveness of mesh adaptation is evident, especially after the second reflection, when the strongest shock occurs. A more quantitative view of the results is given by Fig. 19, which shows the pressure, Mach, and fundamental derivative profiles at two constant locations, in the left and right part of the tube. From the evolution of the fundamental derivative, we can say that the left part of the tube remains for the whole test in the NICFD region, $0<\Gamma<1$, while the right part ends in the supercritical regime. For the location near the left end of the tube, the results of the open-end case are also shown to highlight the velocity decrease imposed by the no-penetration boundary condition. For the location near the right end, the pressure evolution is computed also analytically according to the shock conditions. A fairly good agreement between the analytical and the numerical result is obtained, although a deviation between the pressure values occurs after the second reflection and becomes larger with the third reflection.

This test proves the capability of the proposed adaptive scheme to deal with sudden generation of expansion fans and shock waves in the NICFD regime, as well as with their interaction with solid boundaries. The specification of the maximum tolerated error during the metric construction allows stronger or weaker grid refinement to be achieved according to the strength of the flow phenomena without requiring any user’s action.

6.4 2D oscillating piston in non-classical region

A couple of tests inside and across the non-classical region, at $\Gamma<0$, are now presented. The piston oscillates harmonically inside an infinite-length tube, but with respect to the previous test cases, the oscillation amplitude is now reduced. More specifically, a small amplitude $A=-0.05~{}\mathrm{m}$ is selected to constrain the fluid state to lie inside or near the region at $\Gamma<0$, and a large amplitude $A=-0.15~{}\mathrm{m}$ is tested to generate stronger perturbations. The evolution of the flow field is shown in the $P$–$v$ plane in Fig. 20. The flow state evolves in close proximity to the VLE curve, being the initial state on an isotherm below the critical one ($T_{0}=0.996T_{c}$).

For the motion characterized by the small amplitude, the evolution of different thermodynamic quantities at two fixed locations is displayed in Fig. 21. Since this motion does not introduce strong perturbations, the flow field on the right side of the tube evolves approximately along an isentrope. During the first part of the compression, the flow state moves more inside the non-classical region (see also Fig. 20), reaching a minimum in the value of the fundamental derivative, but during the final part of the compression, $\Gamma$ slightly increases, as the flow is leaving the core of the non-classical region. From here, the successive expansion crosses the core of the non-classical region and ends in the NICFD region, reaching the maximum value of $\Gamma$ (equal to 0.187 for the location $\bar{x_{1}}$), from which the next compression starts. Comparing to the NICFD results of Sec. 6.1, the most notable differences concern the acoustic impedance and the wave propagation speed. Indeed, re-writing the fundamental derivative as $\Gamma=\frac{1}{c}\left(\frac{\partial\rho c}{\partial\rho}\right)_{s}$, if an isentropic expansion starts in the non-classical region and evolves outwards, the acoustic impedance has a maximum at $\Gamma=0$. In this specific test, since the initial value of $\Gamma$ is close to [math], the maximum of the acoustic impedance is extremely close to the initial value, and $\rho c<\rho_{0}c_{0}$ everywhere. A similar evolution can be observed also for the wave speed, which remains almost everywhere smaller, though close, to the initial value. According to Eq. (15), $\mathrm{d}(u+c)$ is positive in the first part of the expansion ($\mathrm{d}P<0$ and $\Gamma<0$) and negative in the second, and vice versa for the compression. Observing that in this test all pressure maxima occurs at $\Gamma<0$ and the minima at $\Gamma>0$, the wave propagation speed has minima when $\mathrm{d}P=0$ and maxima when $\Gamma=0$. It is remarkable that such an “acoustic” behavior is obtained for finite amplitude waves because the value of $\Gamma$ is close to zero Thompson1972 .

Fig. 22 shows the same quantities for the test in which the large oscillation amplitude is imposed. In this test, the flow state already exits the non-classical region during the first compression and most part of the motion occurs in the NICFD region. Therefore, both non-ideal and non-classical phenomena occur. We can observe the non-monotone variation of speed of sound during the isentropic expansion, as seen for Test $\mathcal{NI}$ in Sec. 6.1. As highlighted for the small-amplitude test, the acoustic impedance has minima and maxima for $\Gamma=0$ (see the vertical dotted lines in Fig. 22), but these are not the absolute extrema, which occur in correspondence of the extrema of $\rho$ and $c$. Similarly, the wave propagation speed exhibits local minima and maxima for $\Gamma=0$, but the absolute ones reflect the extrema of the pressure, since all points at which $\mathrm{d}P=0$ fall in the region at $\Gamma>0$.

Also in these two tests, the proposed mesh adaptation strategy allows to automatically modify the grid spacing to increase the solution accuracy, without introducing oscillations due to solution interpolation. The importance of this feature can be appreciated by looking at Fig. 20: some points of the flow field are very close to the VLE curve and oscillations may cause them to enter the two-phase region, where the speed of sound may become negative and ad hoc numerical tools are required for the thermodynamic modeling.

6.5 3D assessment in the NICFD regime

In this final subsection, we briefly present some 3D results to assess the proposed adaptive scheme also for tetrahedral grids. First, we reproduce an oscillating piston in an infinite-length tube, the Test $\mathcal{NI}$ of Sec. 6.1. The initial grid is made of $27806$ nodes and $97259$ tetrahedra. Given that the physical problem is one-dimensional, we can make a comparison between the results obtained in the 2D and 3D simulations. Figure 23 shows the pressure fields at two different times during the second oscillation period, i.e., at $t=1.4~{}T$ and $t=1.8~{}T$, which are the same shown in Fig. 13. A good matching between the 3D and 2D pressure profiles is obtained. Furthermore, as observed already in the 2D tests, the profiles in Fig. 23 are in fact scatter plots of all grid points, and the method correctly reproduces the physical one-dimensional behavior of the flow assessing the validity of the proposed approach.

In this test, a target isotropic metric is defined by combining an adaptation criterion based on the Hessian of the Mach number and another criterion based on the gradient of the pressure. A portion of the resulting grids at the previous times is shown in Fig. 24 to highlight the crucial role of mesh adaptation also in 3D. By comparing them with the pressure fields shown in Fig. 23, the grids result to be refined near the sharper variations and coarse in quasi-uniform regions.

Lastly, a test case of the impulsively started piston is performed. We simulate the same physical features described in Sec. 6.3 to be able to compare 2D and 3D results. As in the previous test, three-dimensional isotropic mesh adaptation is performed by exploiting a compound error estimator that includes the Hessian of the Mach number and the gradient of the pressure. The pressure field and the computational grid are shown in Fig. 25 before any wave interaction with the solid wall, and in Fig. 26 after the reflection of the expansion fan on the left side of the tube and the second shock reflection, which occurs on the piston right face. The effect of mesh adaptation is clear, especially after the reflections of the expansion fan and of the shock. In particular, we can observe the grid portion close to the shock reflected by the right piston wall, which is the strongest phenomenon in the flow field, and the one next to the left end of the tube, where the rarefaction fan is reflected. Because of the different intensity of the phenomena, a different level of refinement is automatically performed, i.e., a smaller grid spacing is reached near the shock, without the user’s intervention.

Although concise, these results show that the proposed method works well also for three-dimensional simulations. From the displayed contour plots and pressure profiles, we can notice that no spurious oscillations are introduced by the interpolation-free adaptation and the physical mono-dimensional behavior is achieved.

7 Conclusions

This paper presents and discusses the assessment of mesh adaptation in unsteady NICFD simulations, since a detailed investigation of the possible extension of standard adaptation techniques, widely used in ideal gas simulations to deal with geometrically complex moving-body problems or to improve the solution accuracy, to such a peculiar thermodynamic regime was missing. Among the available adaptation techniques, we have selected an interpolation-free adaptive scheme capable of representing the volume changes due to mesh adaptation as fictitious continuous deformations of the finite volumes that compose the domain, and of treating them within the ALE framework. In our opinion, this approach represents an advantage for mesh adaptation in the NICFD regime since it prevents the introduction of spurious oscillations due to interpolation, which may unpredictably modify the thermodynamic state of the flow, moving it under the VLE curve.

The proposed method has been assessed through two sets of numerical tests. First, steady simulations of oblique and conical shocks were performed in the unsteady fashion to assess the validity of the numerical scheme. Then, different problems involving a piston moving into a tube filled with $\mathrm{MD_{4}M}$ were simulated. Two initial conditions have been selected along an isotherm above the critical one: Test $\mathcal{D}$ is representative of the dilute regime (approximately $0.9<\Gamma<1$) and Test $\mathcal{NI}$ is more in the core of the NICFD regime. A third initial state has been chosen near the VLE curve, along an isotherm below the critical one, to perform tests across the non-classical and non-ideal region. The thermodynamic behavior of the flow has been described by the polytropic Peng-Robinson EoS. Different simulations of the harmonic motion and the impulsive starting of the piston were performed, both in an infinite-length and in a closed tube.

The expected NICFD and non-classical behavior was well reproduced by the numerical results and mesh adaptation techniques have been successfully exploited to adapt the computational grid to the unsteady solution. In particular, in all tests, the resulting grid was refined near shocks, rarefaction and expansion waves, and their evolution in the flow field was accurately tracked, even if large deformations were imposed on the domain boundaries. At the same time, the adaptation strategy was capable of removing grid nodes where the solution was smooth, without deteriorating the solution accuracy. Also the interactions between waves and solid walls, or other waves, were correctly treated. No spurious oscillations appeared in the results, despite quite large time steps having been enforced.

In conclusion, the proposed method was assessed to perform adaptive unsteady NICFD simulations. Therefore, in future works it could be successfully exploited to perform more complex unsteady simulations. For instance, it could be useful to devise new experimental test rigs to investigate thermo-physical properties of fluid flows in the NICFD regime, or to improve the design of the components of systems operating in the NICFD regime, as Organic Rankine Cycles or supercritical carbon dioxide power cycles.

Acknowledgements.

This study was partially funded by the European Research Council (Consolidator Grant N. 617603, Project NSHOCK, funded under the FP7-IDEAS-ERC scheme).