[chirag]

Dear Foamers,
I am trying to simulate single phase supersonic flow in an axisymmetric cylinder.
(attached image of the domain)
Mesh(made from blockMesh looks good in CheckMesh option (inlet radius is 2mm, and length from depositor to inlet is 360 mm )

Attachment 44702

Create time
Create polyMesh for time = 0
Time = 0
Mesh stats
points: 321133
internal points: 0
faces: 639485
internal faces: 318353
cells: 159719
faces per cell: 5.99701976596
boundary patches: 7
point zones: 0
face zones: 0
cell zones: 0
Overall number of cells of each type:
hexahedra: 159243
prisms: 476
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
polyhedra: 0
Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).
Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology
INLET 15 31 ok (non-closed singly connected)
OUTLET 41 84 ok (non-closed singly connected)
axi_symm-r 159719 160805 ok (non-closed singly connected)
axi_symm-f 159719 160805 ok (non-closed singly connected)
WALL 1179 2362 ok (non-closed singly connected)
DEPOSITOR 459 919 ok (non-closed singly connected)
defaultFaces 0 0 ok (empty)
Checking geometry...
Overall domain bounding box (0 -0.007077245605 0) (0.81 0.007077245605 0.1620955781)
Mesh (non-empty, non-wedge) directions (1 0 1)
Mesh (non-empty) directions (1 1 1)
Wedge axi_symm-r with angle 2.49999994081 degrees
Wedge axi_symm-f with angle 2.49999994081 degrees
All edges aligned with or perpendicular to non-empty directions.
Boundary openness (1.14422845077e-18 3.11071914992e-14 -1.76617364122e-15) OK.
Max cell openness = 2.47634235044e-16 OK.
Max aspect ratio = 22.5016842093 OK.
Minimum face area = 7.74717731676e-10. Maximum face area = 7.06247921513e-05. Face area magnitudes OK.
Min volume = 1.12277931947e-13. Max volume = 4.99740290618e-08. Total volume = 0.000868728818364. Cell volumes OK.
Mesh non-orthogonality Max: 0.571696746183 average: 0.0101664178209
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.330796465233 OK.
Coupled point location match (average 0) OK.
Mesh OK.
End

Following are my B.C:
T

INLET
{
type totalTemperature;
value uniform 11800;
T0 uniform 12000;
U U;
phi phi;
psi thermosi;
gamma 1.67;


}
OUTLET
{
type zeroGradient;
}

DEPOSITOR
{
type fixedValue;
value uniform 300;
}
WALL
{
type fixedValue;
value uniform 300;
}
axi_symm-r
{
type wedge;
}
axi_symm-f
{
type wedge;
}

defaultFaces
{
type empty;
}

P

U
INLET
{
type zeroGradient;
}
OUTLET
{
type zeroGradient;
}
DEPOSITOR
{
type fixedValue;
value uniform (0 0 0);
}
WALL
{
type fixedValue;
value uniform (0 0 0);
}
axi_symm-r
{
type wedge;
}
axi_symm-f
{
type wedge;
}

defaultFaces
{
type empty;
}

k

INLET
{
type fixedValue;
value uniform 0.1544303;
}
OUTLET
{
type inletOutlet;
inletValue uniform 0.1544303;
value uniform 0.1544303;
}
WALL
{
type compressible::kqRWallFunction;
value uniform 0.1544303;
}
DEPOSITOR
{
type compressible::kqRWallFunction;
value uniform 0.1544303;

epsilon

INLET
{
type fixedValue;
value uniform 67.93;
}
OUTLET
{
type inletOutlet;
inletValue uniform 67.93;
value uniform 67.93;
}
DEPOSITOR
{
type compressible::epsilonWallFunction;
value uniform 67.93;
}
WALL
{
type compressible::epsilonWallFunction;
value uniform 67.93;
}

FVSchemes are as follows:

ddtSchemes
{
default Euler;
}

gradSchemes
{
default leastSquares; //Gauss linear;
}

divSchemes
{
default none;
div((-devRhoReff&U)) Gauss linear;
div((muEff*dev2(T(grad(U))))) Gauss linear;

//momentum equation
div(phiNeg,U) Gauss upwind;//vanLeer;
div(phiPos,U) Gauss upwind;//vanLeer;

//energy equation
div(phiNeg,h) Gauss upwind;//vanLeer;
div(phiPos,h) Gauss upwind;//vanLeer;
div(phiNeg,Ek) Gauss upwind;//vanLeer;
div(phiPos,Ek) Gauss upwind;//vanLeer;

//continuity equation
div(phid_neg,p) Gauss upwind;//vanLeer;
div(phid_pos,p) Gauss upwind;//vanLeer;

div(phi,epsilon) Gauss upwind;//vanLeer;
div(phi,k) Gauss upwind;//vanLeer;
}

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default none;// limitedLinear psi 1.0;//none;

interpolate(rho) linear;
interpolate((rho*U)) linear;

reconstruct(psi) vanLeer;
reconstruct(p) vanLeer;
reconstruct(U) vanLeerV;
reconstruct(Dp) vanLeer;
}

snGradSchemes
{
default limited 0.5;//corrected;
}

fluxRequired
{
default none;
p;

A comparison with experimental and CFD is done.
There is a spike observed near the shock region in rhoCentralFoam and PisoCentralFoam.
I have limited the minimum pressure to 0..0000000000001 pa in Fluent as well as in openfoam

forAll(p,cellI) // added
{
if( p[cellI] < 0.0000000000001) // added
{
p[cellI] = 0.0000000000001;// added
}
}; // added


dimensionedScalar unitp("unitp", dimPressure, scalar(0.000000000000001)); // added
p=max(p,unitp); // added

This is not the case with Fluent and it gives good agreement with literature and experimental data. (attached comparison)
Higher order schemes also results in same results as in upwind.

I am unable to figure out why a discontinuity observed in the results? Is there some modification in the solver required when solving flows near vacuum?

Attachment 44703

Attachment 44704

Attachment 44701