# RD-E: 4300 Perfect Gas Modeling with Polynomial EOS

Polynomial EOS is used to model perfect gas. Pressure or energy can be absolute values or relative. Material LAW6 (/MAT/HYDRO) is used to build material cards for each of these cases.

The purpose of this example is to plot numerical pressure, internal energy, and sound speed for a perfect gas material law. Comparison to theoretical results is made.
Polynomial EOS is often used by Radioss to compute hydrodynamic pressure. It is cubic in compression and linear in expansion.(1)
Where,(2) $E=\frac{{E}_{\mathrm{int}}}{{V}_{0}}$
and (3) $\mu =\frac{\rho }{{\rho }_{0}}-1$
Material LAW6 (/MAT/HYDRO) uses this equation to compute hydrostatic pressure. It is possible to consider absolute values or relative variation (Table 1). This example shows how to build material control cards for each of the following cases:
Table 1. Modeling Formulation for Perfect Gas
Case Mathematical Model Pressure Energy
1 $P\left(\mu ,E\right)$ absolute absolute
2 $\text{Δ}P\left(\mu ,E\right)$ relative absolute
3 $\text{Δ}P\left(\mu ,\text{Δ}E\right)$ relative relative
4 $P\left(\mu ,\text{Δ}E\right)$ absolute relative

A simple test of compression/expansion is made to compare these formulation outputs with theoretical results.

## Options and Keywords Used

• Perfect gas
• Polynomial EOS (/EOS/POLYNOMIAL)
• Absolute / Relative formulations
• Pressure shift
• Hydrodynamic fluid material (/MAT/LAW6 (HYDRO or HYD_VISC))
• Imposed displacement (/IMPDISP)

Nodes on each of the faces are moved with imposed displacement

• Boundary conditions (/ALE/BCS)

Boundary nodes are defined as Lagrangian

Element pressure, density and internal energy density are saved in the Time History file.

## Input Files

The input files used in this example include:
Model 1
Model 2
Model 3
Model 4

## Model Description

This test consists with an elementary volume of perfect gas undergoing spherical expansion and compression.
Initial conditions are:
${P}_{0}$
1e5 Pa
${V}_{0}$
1000 m3
${\rho }_{0}$
1.204 $\left[\frac{\text{kg}}{{\text{m}}^{\text{3}}}\right]$
${\mu }_{0}$
0

The fluid will be assumed to be a perfect gas. Volume is changed in the three directions to consider a pure compression $\left(-1<\mu <0\right)$ followed by an expansion of matter $\left(0<\mu \right)$ (Figure 3).

This test will be modeled with a single ALE element (8 node brick) and polynomial EOS.

Evolutions of pressure, internal energy and sound speed will be compared between numerical output and theoretical results.

### Polynomial EOS

Polynomial EOS is used in material LAW6 (/MAT/HYDRO) to compute hydrodynamic pressure. It is cubic in compression and linear in expansion.(4)

Where,
$P$
Hydrodynamic pressure
(5) $E=\frac{E{}_{\mathrm{int}}}{{V}_{0}}$
and(6) $\mu =\frac{\rho }{{\rho }_{0}}-1$
$\left\{{C}_{i}\right\}\text{ }\text{\hspace{0.17em}}i=0.5$ are called hydrodynamic coefficients and they are input flags. Hypothesis on the material behavior allows determining of these coefficients:

This example is focused only on Perfect Gas modeling.

### Model Method

A single ALE brick element is used. Material is confined inside the element by defining brick nodes as Lagrangian. For each face, displacement is imposed on the four nodes along the normal.

Material LAW6 (/MAT/HYDRO) is used and describes the hydrodynamic viscous fluid material.
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
/MAT/LAW6/mat_ID/unit_ID or /MAT/HYDRO/mat_ID/unit_ID
mat_title
${\rho }_{i}$
$\upsilon$
C0 C1 C2 C3
Pmin Psh
C4 C5 E0

### Pressure Shift

Material LAW6 introduces flag Psh which allows shifting computed pressure in the polynomial equation of state:(7)

Radioss Engine shifts C0 flag and computed pressure $P\left(\mu ,E\right)$ with an offset of -Psh.

### Minimum Pressure

(8) ${P}_{\mathrm{min}}=\underset{\mu \to -1}{\mathrm{lim}}\text{\hspace{0.17em}}P\left(\mu ,E\right)-{P}_{sh}$

The theoretical value is ${P}_{min}=0Pa$ (absolute pressure) with a default value of -1030, to accept a negative value in relative pressure formulation.

This flag has to be manually offset with -Psh.

## Results

### Theoretical Results

The purpose of this section is to plot pressure, internal energy, and sound speed in function of the single parameter $V$ or $\mu$.
1. Pressure
Perfect gas pressure is given by:(9) $PV=\left(\gamma -1\right){E}_{\mathrm{int}}$
Then,(10) $dP\left(V,{E}_{\mathrm{int}}\right)={\frac{\partial P}{\partial V}|}_{{E}_{\mathrm{int}}}dV+{\frac{\partial P}{\partial {E}_{\mathrm{int}}}|}_{V}d{E}_{\mathrm{int}}$
Radioss assumes the hypothesis of an isentropic process to compute the change in internal energy:(11) $d{E}_{\mathrm{int}}=-PdV$
This theory gives the following differential equation:(12) $\frac{dP}{dV}=-\frac{\gamma P}{V}$
This has the form $y\text{'}+\gamma /x=0$ and the general solution is:(13) $y=Cst.{x}^{-\gamma }$
Pressure is also polytropic:(14) $P{V}^{\gamma }=P{}_{0}V{{}_{0}}^{\gamma }$ (15) $P\left(V\right)={P}_{0}{\left(\frac{{V}_{0}}{V}\right)}^{\gamma }$

Here, $\gamma$ is the material constant (ratio of heat capacity). For diatomic gas $\gamma$ =1.4. Air is made mainly of diatomic gas, so set gamma to 1.4 for air.

2. Internal Energy
Equation 9 and Equation 15 lead to the immediate result:(16) ${E}_{\mathrm{int}}\left(V\right)=\frac{{P}_{0}{V}_{0}{}^{\gamma }}{\left(\gamma -1\right){V}^{\gamma -1}}$
3. Sound Speed
Perfect gas sound speed is:(17) $c=\sqrt{\gamma P/\rho }$
Equation 15 gives its expression in term of volume:(18) $c=\frac{\gamma P{}_{0}}{{\rho }_{0}}{\left(\frac{{V}_{0}}{V}\right)}^{\gamma -1}$
The theoretical results are listed in the table below. Pressure, internal energy, and sound speed are expressed both in function of $V$ and $\mu$.
Pressure (Pa) Internal Energy Density (J) Sound Speed (m/s)
PREF(V) PREF($\mu$) $\rho$eREF(V) $\rho$eREF($\mu$) cREF(V) cREF($\mu$)
${P}_{0}{\left(\frac{{V}_{0}}{V}\right)}^{\gamma }$ ${P}_{0}{\left(1+\mu \right)}^{\gamma }$ $\frac{{P}_{0}}{\left(\gamma -1\right)}{\left(\frac{{V}_{0}}{V}\right)}^{\gamma -1}$ $\frac{{P}_{0}}{\left(\gamma -1\right)}{\left(1+\mu \right)}^{\gamma -1}$ $\sqrt{\frac{\gamma {P}_{0}}{{\rho }_{0}}{\left(\frac{{V}_{0}}{V}\right)}^{\gamma -1}}$ $\sqrt{\frac{\gamma {P}_{0}}{{\rho }_{0}}{\left(1+\mu \right)}^{\gamma -1}}$
Corresponding plots are:

### Material Control Cards

Material is supposed to be a perfect gas. The following cases have been investigated:
• Case 1: Both Pressure and Energy are absolute values: $P\left(\mu ,E\right)$
• Case 2: Pressure is relative and Energy is absolute: $\text{Δ}P\left(\mu ,E\right)$
• Case 3: Both Pressure and Energy are relative: $\text{Δ}P\left(\mu ,\text{Δ}E\right)$
• Case 4: Pressure is absolute and Energy is relative: $P\left(\mu ,\text{Δ}E\right)$

### Sound Speed and Time Step

Material law 6 computes sound speed through the usual expression for fluids:(19) ${c}^{2}=\frac{dP}{d\rho }$
It can be written in function of $\mu$:(20) $\mu =\frac{\rho }{{\rho }_{0}}-1⇒\frac{1}{d\rho }=\frac{1}{{\rho }_{0}}\frac{1}{d\mu }$
Then,(21) ${c}^{2}=\frac{1}{{\rho }_{0}}\frac{dP}{d\mu }$
The total differential of $P$ in terms of internal energy $E$ and $\mu$ is:(22) $dP\left(\mu ,E\right)={\frac{\partial P}{\partial \mu }|}_{E}\text{\hspace{0.17em}}d\mu +{\frac{\partial P}{\partial E}|}_{\mu }dE$
In case of an isentropic transformation (reversible and adiabatic), the change of internal energy ${E}_{\mathrm{int}}$ with volume $V$ and pressure $P$ is given by:(23) $d{E}_{\mathrm{int}}=-PdV$
Using relation which links ${E}_{\mathrm{int}}$ and $E$ leads to:(24) $dE=-\frac{P}{{V}_{0}}dV$
$\mu$ can be expressed in terms of volume ratio:(25) $\mu =\frac{{\upsilon }_{0}}{\upsilon }-1$
Its variation in function of the volume change is also:(26) $d\mu =-\frac{{V}_{0}}{{V}^{2}}dV=-\frac{{\left(1+\mu \right)}^{2}}{{V}_{0}}dV$
Change in internal energy per unit volume $E$ is then:(27) $dE=-\frac{P}{{\left(1+\mu \right)}^{2}}d\mu$ (28) $\frac{dP\left(\mu ,E\right)}{d\mu }={\frac{\partial P}{\partial \mu }|}_{E}+\frac{P}{{\left(1+\mu \right)}^{2}}{\frac{\partial P}{\partial E}|}_{\mu }$
Finally, the sound speed is given by:(29) ${c}^{2}=\frac{1}{{\rho }_{0}}{\frac{\partial P}{\partial \mu }|}_{E}+\frac{P}{{\rho }_{0}{\left(1+\mu \right)}^{2}}{\frac{\partial P}{\partial E}|}_{\mu }$
This expression computes the sound speed for a given equation of state $P\left(\mu ,E\right)$. In the case of perfect gas, it was shown that for each type of formulation (absolute or relative), EOS can be written:(30) $P\left(\mu ,E\right)={C}_{0}+{C}_{1}\mu +\left({C}_{4}+{C}_{5}\mu \right)E$
Equation 29 is used to compute sound speed:(31) ${\frac{\partial P}{\partial \mu }|}_{E}={C}_{1}+{C}_{5}E$ (32) ${\frac{\partial P}{\partial E}|}_{\mu }={C}_{4}+{C}_{5}\mu$ (33) ${c}^{2}={\frac{{C}_{1}+{C}_{5}E}{\rho 0}}_{}+\frac{{C}_{0}+{C}_{1}\mu +\left({C}_{4}+{C}_{5}\mu \right)E}{{\rho }_{0}{\left(1+\mu \right)}^{2}}\left({C}_{4}+{C}_{5}\mu \right)$
This calculation is then applied for each of the four cases.
Table 2. Numerical Sound Speed vs. Theoretical Expression
Case C0 C1 C4 C5 c2

Comparison with Theoretical Value
1 0 0 $\gamma -1$ $\gamma -1$ $\frac{\gamma \left(\gamma -1\right)E}{{\rho }_{0}}$ c = cREF
2 0 0 $\gamma -1$ $\gamma -1$ $\frac{\gamma \left(\gamma -1\right)E}{{\rho }_{0}}$ c = cREF
3 ${E}_{0}\left(\gamma -1\right)$ ${E}_{0}\left(\gamma -1\right)$ $\gamma -1$ $\gamma -1$ $\frac{\gamma \left(\gamma -1\right)\left(E+{E}_{0}\right)}{{\rho }_{0}}$ c = cREF
4 ${E}_{0}\left(\gamma -1\right)$ ${E}_{0}\left(\gamma -1\right)$ $\gamma -1$ $\gamma -1$ $\frac{\gamma \left(\gamma -1\right)\left(E+{E}_{0}\right)}{{\rho }_{0}}$ c = cREF

For each of the four formulations, the computed sound speed by Radioss is the same as the theoretical one. Time step and cycle number are also not affected.

### Case 1: Both Pressure and Energy are Absolute Values

1. Equation of State
Equation of state can be written:(34) $P=\left(\gamma -1\right)\frac{{E}_{\mathrm{int}}}{V}=\left(\gamma -1\right)\left(1+\mu \right)\frac{{E}_{\mathrm{int}}}{{V}_{0}}$
with(35) ${{E}_{\mathrm{int}}|}_{t=0}={E}_{0}{V}_{0}=\frac{{P}_{0}{V}_{0}}{\gamma -1}$
Expanding this expression and identifying the polynomial coefficients leads to:(36) $P\left(\mu ,E\right)=\left({C}_{4}+{C}_{5}\mu \right)E$

Where, ${C}_{4}=C{}_{5}\text{\hspace{0.17em}}=\left(\gamma -1\right),\text{ }{E}_{0}=\frac{{P}_{0}}{\gamma -1},\text{ }{P}_{sh}=0$

2. Corresponding Input
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
/MAT/LAW6/mat_ID/unit_ID or /MAT/HYDRO/mat_ID/unit_ID
AbsolutePRESSURE_AbsoluteENERGY
${\rho }_{i}$
$\upsilon$
0 0 0 0
0 0
C4 = $\gamma -1$1 C5 = $\gamma -1$ ${E}_{0}=\frac{{P}_{0}}{\gamma -1}$
3. Output Results
Table 3.
Time History Measure Initial Value Unit
/TH/BRICK (P) P P0 Pressure
/TH (IE) ${E}_{\mathrm{int}}\left(=E\text{\hspace{0.17em}}x\text{\hspace{0.17em}}{V}_{0}\right)$ ${E}_{0}{V}_{0}$ Energy
/TH/BRICK (IE) ${E}_{\mathrm{int}}/V$ ${E}_{0}$ Pressure
4. Comparison with Theoretical Result
Numerical result for perfect gas pressure is given by time history. Element time history (/TH/BRICK) allows displaying it. This result is compared to a theoretical one. Curves are superimposed.
Internal energy can be obtained through two different ways. The first one is internal energy density (${E}_{\mathrm{int}}/V$) recorded by element time history (/TH/BRICK). The second one is the internal energy from the global time history $\left(\sum {}_{elem}\text{\hspace{0.17em}}{E}_{\mathrm{int}}\right)$ because the model is composed of a single element.

### Case 2: Pressure is Relative and Energy is Absolute

1. Equation of State
Equation of state for a perfect gas is:(37) $P\left(\mu ,E\right)=\left(\gamma -1\right)\left(1+\mu \right)\frac{{E}_{\mathrm{int}}}{{V}_{0}}$
Calculating Pressure from a reference one provides relative pressure:(38) $\text{Δ}P=P\left(\mu ,E\right)-{P}_{0}=\left(\gamma -1\right)\left(1+\mu \right)\frac{{E}_{\mathrm{int}}}{{V}_{0}}-{P}_{0}$
Expanding this expression and identifying with polynomial coefficients leads to:(39) $\text{Δ}P\left(\mu ,E\right)=P\left(\mu ,E\right)={P}_{sh}=-{P}_{sh}+\left({C}_{4}+{C}_{5}\mu \right)E$

Where, ${C}_{4}=C{}_{5}\text{\hspace{0.17em}}=\left(\gamma -1\right),\text{ }{E}_{0}=\frac{{P}_{0}}{\gamma -1},\text{ }{P}_{sh}=0$

2. Minimum Pressure(40) $P\left(\mu ,E\right)\ge 0⇒\text{Δ}P\ge -{P}_{0}$

Then, the minimum pressure must be set to a non-zero value ${P}_{\mathrm{min}}=-{P}_{0}$.

3. Corresponding Input
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
/MAT/LAW6/mat_ID/unit_ID or /MAT/HYDRO/mat_ID/unit_ID
RelativePRESSURE_AbsoluteENERGY
${\rho }_{i}$
$\upsilon$
0 0 0 0
-P0 P0
C4 = $\gamma -1$ C5 = $\gamma -1$ ${E}_{0}=\frac{{P}_{0}}{\gamma -1}$
4. Output Result
Time History Measure Initial Value Unit
/TH/BRICK (P) $\text{Δ}P$ 0 Pressure
/TH (IE) ${E}_{\mathrm{int}}\left(=E\text{\hspace{0.17em}}x\text{\hspace{0.17em}}{V}_{0}\right)$ ${E}_{0}{V}_{0}$ Energy
/TH/BRICK (IE) ${E}_{\mathrm{int}}/V$ ${E}_{0}$ Pressure
5. Comparison with Theoretical Result
Element time history (/TH/BRICK) is the pressure relative to Psh. The resulting curve is then shifted with Psh value and starts from 0.
Internal energy can be obtained through two different ways. The first one is internal energy density (${E}_{\mathrm{int}}/V$) recorded by element time history (/TH/BRICK). The second one is the internal energy from the global time history $\left(\sum {}_{elem}\text{\hspace{0.17em}}{E}_{\mathrm{int}}\right)$ because the model is composed of a single element.

### Case 3: Both Pressure and Energy are Relative

1. Equation of State
Equation of state for a perfect gas is:(41) $P\left(\mu ,E\right)=\left(\gamma -1\right)\left(1+\mu \right)\frac{{E}_{\mathrm{int}}}{{V}_{0}}$
Initial internal energy can be introduced:(42) ${E}_{\mathrm{int}}={E}_{\mathrm{int}}+\left({E}_{\mathrm{int}}{}_{|t=0}-{E}_{\mathrm{int}}{}_{|t=0}\right)$
Calculating pressure from a reference one provides:(43) $P\left(\mu ,E\right)-{P}_{0}=\text{Δ}P=\left(\gamma -1\right)\left(1+\mu \right)\left(\text{Δ}E+{E}_{0}\right)-{P}_{0}$
Where,(44) $\text{Δ}E=\frac{{E}_{\mathrm{int}}-{E}_{\mathrm{int}}{}_{|t=0}}{{V}_{0}}\text{\hspace{0.17em}}{E}_{0}={\frac{{E}_{\mathrm{int}}{}_{|t=0}}{{V}_{0}}}_{}$
Expanding this expression and identifying with polynomial coefficients leads to:(45) $\text{Δ}P\left(\mu ,\text{Δ}E\right)=P\left(\mu ,E\right)-{P}_{sh}={C}_{0}-{P}_{sh}+{C}_{1}\mu +\left({C}_{4}+{C}_{5}\mu \right)\text{Δ}E$
Where, ${C}_{0}={C}_{1}={E}_{0}\left(\gamma -1\right)$ ${C}_{4}={C}_{5}=\gamma -1$ $\text{Δ}{E}_{0}=0$ ${P}_{sh}={P}_{0}$
2. Minimum Pressure(46) $P\left(\mu ,E\right)\ge 0⇒\text{Δ}P\ge -{P}_{0}$

The minimum pressure must be set to a non-zero value ${P}_{\mathrm{min}}=-{P}_{0}$.

3. Corresponding Input
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
/MAT/LAW6/mat_ID/unit_ID or /MAT/HYDRO/mat_ID/unit_ID
RelativePRESSURE_RelativeENERGY
${\rho }_{i}$
$\upsilon$
${E}_{0}\left(\gamma -1\right)$ ${E}_{0}\left(\gamma -1\right)$ 0 0
-P0 P0
C4 = $\gamma -1$ C5 = $\gamma -1$ 0
4. Output Results
Time History Measure Initial Value Unit
/TH/BRICK (P) $\text{Δ}P$ 0 Pressure
/TH (IE) $\text{Δ}{E}_{\mathrm{int}}\left(=\text{Δ}E\text{ }x\text{ }{V}_{0}\right)$ 0 Energy
/TH/BRICK (IE) $\text{Δ}{E}_{\mathrm{int}}/V$ 0 Pressure
5. Comparison with Theoretical Result
Element time history (/TH/BRICK) is the pressure relative to Psh. The resulting curve is then shifted with Psh value and starts also from 0.
Internal energy can be obtained through two different ways. The first one is internal energy density (${E}_{\mathrm{int}}/V$) recorded by element time history (/TH/BRICK). The second one is the internal energy from the global time history $\left(\sum {}_{elem}\text{\hspace{0.17em}}{E}_{\mathrm{int}}\right)$ because the model is composed of a single element. This numerical internal energy is relative to its initial value; it is shifted with the ${E}_{0}{V}_{0}$ value from the absolute theoretical one and also starts from 0.

### Case 4: Pressure is Absolute and Energy is Relative

1. Equation of State
Equation of state for a perfect gas is:(47) $P\left(\mu ,E\right)=\left(\gamma -1\right)\left(1+\mu \right)\frac{{E}_{\mathrm{int}}}{{V}_{0}}$
Initial internal energy can be introduced:(48) ${E}_{\mathrm{int}}={E}_{\mathrm{int}}+\left({E}_{\mathrm{int}}{}_{|t=0}-{E}_{\mathrm{int}}{}_{|t=0}\right)$
Which leads to:(49) $P\left(\mu ,E\right)=\left(\gamma -1\right)\left(1+\mu \right)\left({E}_{0}+\text{Δ}E\right)$
Expanding this expression and identifying with polynomial coefficients leads to:(50) $P\left(\mu ,E\right)={C}_{0}+{C}_{1}\mu +\left({C}_{4}+{C}_{5}\mu \right)\text{Δ}E$
Where, ${C}_{0}={C}_{1}={E}_{0}\left(\gamma -1\right)$ ${C}_{4}={C}_{5}=\gamma -1$
2. Corresponding Input
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
/MAT/LAW6/mat_ID/unit_ID or /MAT/HYDRO/mat_ID/unit_ID
AbsolutePRESSURE_RelativeENERGY
${\rho }_{i}$
$\upsilon$
${E}_{0}\left(\gamma -1\right)$ ${E}_{0}\left(\gamma -1\right)$ 0 0
0 0
C4 = $\gamma -1$ C5 = $\gamma -1$ 0
3. Output Results
Time History Measure Initial Value Unit
/TH/BRICK (P) P P0 Pressure
/TH (IE) $\text{Δ}{E}_{\mathrm{int}}\left(=\text{Δ}E\text{ }x\text{ }{V}_{0}\right)$ 0 Energy
/TH/BRICK (IE) $\text{Δ}{E}_{\mathrm{int}}/V$ 0 Pressure
4. Comparison with Theoretical Result
Element time history (/TH/BRICK) gives absolute pressure. This result is compared to a theoretical one. Curves are superimposed.
Internal energy can be obtained through two different ways. The first one is internal energy density ($\text{Δ}{E}_{\mathrm{int}}/V$) recorded by element time history (/TH/BRICK). The second one is the internal energy from the global time history $\left(\sum {}_{elem}\text{\hspace{0.17em}}\text{Δ}{E}_{\mathrm{int}}\right)$ because the model is composed of a single element. This numerical internal energy is relative to its initial value; it is shifted with the ${E}_{0}{V}_{0}$ value from the absolute theoretical one and also starts from 0.