Modelling Sod's Shock tube Problem and Comparison study between Two different solvers using OpenFOAM
A) INTRODUCTION:
This project focuses on simulating Sod's ShockTube problem. The problem involves simulating the propagation of a normal shockwave inside a shock tube under experimental conditions.
The shockwave is simulated using OpenFOAM using two different solvers: rhoCentralFoam (Density-based) and rhoPimpleFoam [sonicFoam] (Pressure-based). The main reason for selecting the two particular solvers is because usually transport phenomena (such as shockwaves) are approached using density-based and pressure-based solvers.
A1) AIM & OBJECTIVES:
A) To characterize the obtained shockwave
B) Comparison study between both solvers and Grid Independence study.
C) Simulation time comparison of both the solvers
A2) MOTIVATION:
Anything moving faster than the speed of sound could create shockwaves. The study of shockwaves is a very important subject for design purposes. Study of shockwaves characteristics like the momentum, temperature, and strength, energy flux, and position will be crucial in designing the chamber/ structure accordingly. It has wide industrial applications such as in rocket nozzle, wind-tunnels, simulation of explosions and detonations, combustion chambers, etc.
B) BACKGROUND:
B1: SHOCK TUBE
The Shocktube is a long tube that is built with metal which consists of two chambers separated by a diaphragm. The two chambers consist of two gases with a large pressure ratio. One chamber is filled with a gas of higher pressure with is called the Driver Gas and vice-versa is done at the other chamber called the driven gas. Hence the chambers are also termed as Driver and Driven.
The diaphragm is inbuilt with a puncturing mechanism that breaks the diaphragm under predetermined conditions. When the contact is suddenly broken there is a series of pressure waves generated which causes a shockwave.
B2: SHOCKWAVE GENERATION
The breaking of the diaphragm which separates both the driver and driven gas causes a drastic onset of compression pressure waves towards the direction of the low-pressure region due to a sudden gradient in pressure.
The leading pressure wave is called the Shockwave. This shockwave is followed by a series of further pressure waves which raises the temperature and pressure of the driven gas. Also, a sound wave is generated which keeps increasing in velocity. But the leading wave has the maximum velocity.
Simultaneously, a rarefaction wave is set up in the opposite direction of the shockwave which travels towards the driver chamber. This wave is called the Prandtl-Meyer Expansion wave.
The scope of this project is only restricted to the propagation of shockwave from formation until it travels to the end of the shock tube.
B3: GOVERNING EQUATIONS:
For simplicity, the flow is considered to be Inviscid and Laminar. Under Inviscid conditions, the Navier Stokes equation gets reduced to Euler’s equations which can also be applied for flows whose Reynolds number is greater than 1.
The following expressions represent the 1D Euler’s equations of gas dynamics in conservative form.
Where `rho` is density, p is pressure and u is velocity, `m = rho*u` is the momentum and e is the Internal energy given by `e = frac{p}{rho*(gamma-1)}`, gamma is the ratio of specific heats = 1.4
This problem is also a Riemann problem. A Riemann problem is a set of non-linear hyperbolic PDEs which with a piecewise constant Initial Condition. It is this Initial condition that makes this problem unique and is called Sod’s shock tube as it was investigated by Gary A. Sod in 1978. The initial conditions are given by :
B4: REGIONS OF SHOCKWAVE
Shown above is the shockwave shorty after the diaphragm has been punctured. Five different regions will be formed when the diaphragm is burst.
Region 1 & 5 - left and right sides of initial states
Region 2 - Expansion (or) Rarefaction wave
Regions 3 & 4 - Steady states which are separated by a Contact Discontinuity line between 3 and 4 separate fluids of different entropy (but they have the same pressure and velocity)
B5: SOFTWARE AND SOLVERS USED:
The software used for this project is OpenFOAM. OpenFOAM is an open-source CFD software that has been developed by OpenFOAM LTD. It has a wide variety of solvers that can be used to simulate and solve real-time engineering problems.
Shockwave propagation is a transport phenomenon and classified as a high-speed unsteady problem. Solving transport phenomena problems is usually done using a Density-based approach and a Pressure based approach. Hence in this project, both types of solvers are employed to simulate the problem.
The solvers used in this project are 1) Density-based - rhoCentralFoam and 2) Pressure based rhoPimpleFoam [sonicFoam]
rhoCentralFoam
This solver takes into account not only the properties transported by the flow but also the propagation of waves. For this, flux interpolations that take into account the direction of transportation are needed. This solver uses the upwind direction schemes for calculating its volumetric fluxes based on the formula: `alpha = frac{ phiunderset(f+)() }{ phiunderset(f+)() - phiunderset(f-)() }` with the resulting value of `alpha` as 0.5; Hence this solver works on Central upwind differencing Scheme. A flux limiter function `beta(r)` where 'r' is the gradient of successive interpolated variable, is used to switch between the upper-order and the lower-order schemes which provide a reliable solution during shock and discontinuities.
rhoPimpleFoam
Originally sonicFoam was a separate solver on OpenFOAM. Now, the sonicoam has now been merged with the solver rhoPimpleFoam in the latest version of Openfoam. sonicFoam can now be run through rhoPimpleFoam using the transonic mode available in fvSolution file present in the solver's data files.
SonicFoam uses a Non-Iterative method for coupling the pressure and velocity parameters implicitly. This algorithm is called as PIMPLE Algorithm which is a combination of PISO[Pressure Implicit with splitting of operators] and SIMPLE Algorithm [Semi-Implicit method for splitting of Operators] . The splitting of parameters is done and decoupling of pressure from velocity takes place. One other advantage of PISO is that siince it is Implicit, it can handle large values of timesteps. At places of shocks and discontinuties, the solver is expected to produce some oscillations.
B6: ANALYTICAL REFERENCE:
These are the exact characteristic plots obtained from Sod's Shock Tube problem. Parameters are in the non-dimensional form. The reference plots obtained were obtained from the official matlab webiste: http://www.mathworks.com/matlabcentral/fileexchange/46311-sod-shock-tube-problem-solver
D) SOLVING SOD'S SHOCKTUBE PROBLEM
D1) GEOMETRY, INITIAL AND BOUNDARY CONDITIONS:
The geometry of the shocktube is a cuboid with dimenions 1m long X 0.5m wide X 0.2m deep. But essentially the simulation is 2D in nature. The geometry is coded within the blockMeshDict file inside the 'system' folder of the solver. Four different meshes resolutions are done to make comparison between the results with varying mesh size. The created meshes are 100X60, 200X60, 300X60 and 400X60 along the length and height of the geometry respectively.
The front and back walls as seen from the geometry below are given 'EMPTY' boundary condition as the simulation is 2D based. The remaining boundaries are given the value 'WALL'.
The blockMeshDict file is shown below for 300X60 mesh resolution:
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 7
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(-0.5 0 0)
(0.5 0 0)
(0.5 0.5 0)
(-0.5 0.5 0)
(-0.5 0 0.2)
(0.5 0 0.2)
(0.5 0.5 0.2)
(-0.5 0.5 0.2)
);
blocks
(
hex (0 1 2 3 4 5 6 7)(300 60 1) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
sides
{
type wall;
faces
(
(0 4 7 3)
(2 6 5 1)
);
}
topandbottom
{
type wall;
faces
(
(5 6 7 4)
(3 2 1 0)
);
}
empty
{
type empty;
faces
(
(1 5 4 0)
(3 7 6 2)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
The obtained geometry and mesh is shown below:
D2) INITIAL AND BOUNDARY CONDITIONS
The boundary conditions are applied on the geometry in accordance with the Sod's Shock Tube problem.
Both the solvers take Pressure, Temperature and Velocity as Initial conditions. Pressure and velocity is applied directly as a list of Non-uniform field. Temperature is calculated using the Ideal / Perfect gas equation P = ρRT. Hence a non-uniform list of Temperature is obtained.
The Initial conditions files are shown below:
Pressure File:
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 7
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -1 -2 0 0 0 0];
internalField nonuniform List
18000
(
1
1
..150 times
0.1
0.1
..150 times
Above values repeated for 60 times
);
boundaryField
{
sides
{
type zeroGradient;
}
topandbottom
{
type zeroGradient;
}
empty
{
type empty;
}
}
// ************************************************************************* //Temperature File:
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 7
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 1 0 0 0];
internalField nonuniform List
(
0.00348432
0.00348.432
..150 times
0.00278746
0.00278746
..150 times
Above values repeated for 60 times
);
boundaryField
{
sides
{
type zeroGradient;
}
topandbottom
{
type zeroGradient;
}
empty
{
type empty;
}
}
// ************************************************************************* //Velocity File:
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 7
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
sides
{
type slip;
}
topandbottom
{
type slip;
}
empty
{
type empty;
}
}
// ************************************************************************* //
D3) PROGRAM CONTROL AND SIMULATION MODEL
Assuming an initial Courant Number of 0.3, and calculating the timestep based on CFL condition; dx = Length/No of divisions = 1/300 = 0.0034; Assuming an initial shockwave speed of 375m/s [this is the standard speed of shock in the air]; Calculating timestep with these parameters;
We get the timestep as 2.67e-6. Next, the predicted solution is obtained at 0.1 to 0.2 seconds. Hence setting simulation time as 0.2 seconds, the controlDict file is specified.
The controlDict file is shown below:
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 7
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application rhoCentralFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 0.2;
deltaT 2.67e-6;
writeControl adjustableRunTime;
writeInterval 0.005;
cycleWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
adjustTimeStep yes;
maxCo 0.2;
maxDeltaT 0.01;
functions
{
#includeFunc mag(U)
}
// ************************************************************************* //
The fvSchemes and fvSolution files are left unchanged in the problem and are as per rhoCentralFoam default.
The gas that is used in the driver and the driven chamber is assumed to be the same with changes in only the Pressure and Density of the gas according to the chambers. The thermophysical properties of the gas are left unchanged in the rhoCentralFoam solver as is shown below:
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 7
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
thermoType
{
type hePsiThermo;
mixture pureMixture;
transport const;
thermo hConst;
equationOfState perfectGas;
specie specie;
energy sensibleInternalEnergy;
}
mixture
{
specie
{
molWeight 28.96;
}
thermodynamics
{
Cp 1004.5;
Hf 2.544e+06;
}
transport
{
mu 0;
Pr 1;
}
}
// ************************************************************************* //Next, to ensure the simulation follows a laminar model of flow, the same is specified in the Turbulence model file of the solver.
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 7
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType laminar;
// ************************************************************************* //The above steps are repeated for the rhoPimpleFoam solver as well. There is no change in the geometry and Initial conditions file. But the fvSchemes and fvSolution file are as per default of rhoPimpleFoam solver.
E) RESULTS:
E1) CHARACTERIZATION OF SHOCKWAVE
CHARACTERISTIC PLOTS OBTAINED AT TIME = 0.15s USING [rhoCentralFoam]
ANIMATIONS OF PLOTS AND CONTOURS:
PRESSURE:
DENSITY:
VELOCITY:
INFERENCE
a) PRESSURE PLOT
- In the pressure animation, we can see that there exists a rarefaction wave that moves towards the left and a shock front that moves towards the right.
- The rarefaction wave (Expansion wave) reduces the pressure field w.r.t. its propagation in the Driver section and the shockwave increases the pressure field w.r.t. its propagation in the Driven section.
- The pressure increases with the propagation of the shockwave because, as it moves with high speed, the gas particles in the way of the wave do not get time to move away and it gets pushed along with the initial wave.
- The initial wave is followed by a series of subsequent pressure waves.
b) DENSITY PLOT
- In the density animation, all regions of shockwave are observed clearly including the Contact discontinuity which is the region of separation of the two different media of different entropy.
- In the driver chamber, it can be seen that the density reduces w.r.t. rarefaction wave propagation.
- Towards the right, two sharp steps in density is observed, The first step corresponds to the contact discontinuity which separates the high-pressure gas and the low-pressure gas that moves along the direction of the shock; the second one corresponds to the shockfront.
c) VELOCITY PLOT
- When the diaphragm is broken, there is a sudden spike in the velocity towards both sides of the middle section.
- In the velocity plot, it can be observed that there are two shockfronts moving towards the left and right.
- From the midsection, the plot that moves towards the right side is the shockfront of the shockwave and the one that moves towards the left side is the wavefront of the rarefaction wave.
- It can be observed that the speed of the shockfront is much higher than the speed of the expansion front.
E2) GRID INDEPENDANCE STUDY
[WITH RHOCENTRALFOAM]
INFERENCE:
- As the number of mesh grids increases, we get an increase in the accuracy of the plot. This is clearly visible in the plot above.
- The plot with 100 grid points has blunt corners and the sharpness of the corners increases as the mesh resolution increases.
- We can see that plots of mesh resolutions of 300 and 400 almost lie on each other and hence we can say that at this point we have achieved grid independence at mesh resolution of 300X60.
- Further, if we try to increase the number of grid points, there will be very little chance of values but there will be only a drastic increase in the simulation time.
E3) COMPARISON STUDY BETWEEN THE SOLVERS
E3) TIME COMPARISON BETWEEN THE SOLVERS
For mesh resolution of 300X60, the both solvers were run on the same initial conditions:
The simulation time for rhoCentralFoam was ~37 seconds as shown below:
The simulation time for rhoPimpleFoam was ~54 seconds as shown below:
From the above, it can be inferred that RhoCentralFoam is ~45.94% faster than rhoPimpleFoam in this case.
F) INFERENCE:
- From a bird's eye's point of view, both the solvers are able to develop a fairly accurate plot with similar overall trends to the reference but both the solvers are not able to reproduce the squareness near the corners in the plots as compared to the reference.
- It can be observed that rhoCentralFoam is closer to the analytical solution and able to produce sharp characteristic turns. This is due to the fact that rhoCentralFoam uses a Central Upwind differencing scheme with flux limiters that gives accurate values in regions of discontinuities and shockwave.
- There are some oscillations produced in the case of rhoCentralFoam but not in rhoPimpleFoam. This is a point of controversy as rhoCentralFoam is designed to handle shocks and discontinuities.
- This is due to the Numerical diffusion which might have been produced due to the jumping from higher-order schemes to lower-order near the shockfront. As rhoPimplefoam does not use any flux limiters, there are no such effects observed.
- rhoPimpleFoam is less accurate than rhoCentralFoam and the characteristic shapes are observed with blunt turns. This is because the solver runs on the PISO algorithm which was initially designed only for incompressible flows which was then further extended to compressible flows as well. Hence it will not be able to capture the minor details as accurately as rhoCentralFoam.
- Also, rhoPimpleFoam has a higher simulation time due to the extra PISO loop. Rhocentalfoam computes the same problem 45.94% faster.
- Hence, undoubtedly, rhoCentralFoam is a better solver for shock tube applications in OpenFOAM.
G) SUMMARY AND FUTURE SCOPE:
- The formation of the shockwave is characterized and studied using temperature, pressure, density, and velocity plots.
- RhoCentralFoam is more accurate than sonicFoam as rhoCentralFoam is designed for compressible flows which work on Central Upwind Differencing Schemes with the use of Flux Limiters that obtain the numerical values close to the exact values in regions of shocks and discontinuities.
- RhocentralFoam is faster to solve the same set of equations.
- Grid Independence is obtained at a mesh resolution of 300X60 and a further increase in the number of grids will only increase the simulation time.
- Hence it can be concluded by saying that RhoCentralFoam is a better solver for Shocktube simulation as well as any transonic problem with shocks and discontinuities in OpenFoam as it is more accurate and has a lesser simulation time, both of which are essential in any application.
- This study can further be extended to higher pressure ratios between the driver and driven chambers. Also, other inputs such as viscosity, turbulence, and Reynolds number can be varied and these effects can be recorded using rhoCentralFoam.
REFERENCES:
- A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws by Gary A. Sod (1978)
- TEXT: Riemann problem and discontinuities: the shock tube problem
- Greenshields, Christopher J. and Weller, Henry G. and Gasparini, Luca and Reese, Jason M., 2010