[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 )

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;
}

dimensions [1 -1 -2 0 0 0 0];

internalField uniform 100;

boundaryField
{
/* INLET
{
type fixedValue;
value uniform 20000;
}

OUTLET
{
type fixedValue;
value uniform 100;
}
*/
INLET
{
type totalPressure;
p0 uniform 20000;
U U;
phi phi;
rho rho;
psi none;
gamma 1.67;
value uniform 20000;
}

OUTLET
{
type waveTransmissive;
field p;
phi phi;
rho rho;
psi thermosi;
gamma 1.67;
fieldInf 100;
lInf 0.0001;
value uniform 100;
}

WALL
{
type zeroGradient;
}

DEPOSITOR
{
type zeroGradient;
}

axi_symm-r
{
type wedge;
}

axi_symm-f
{
type wedge;
}

defaultFaces
{
type empty;
}
}

U

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.
This is not the case with Fluent and it gives good agreement with literature and experimental data.
I am unable to figure out why a discontinuity observed in the results?

comparison-1.png

meshmain.jpg

Vmag-OF1.jpg

pisocent_foam.zip

[LukaD]

Hey there,

Can you also give the boundary conditions for pressure? It looks like you put those of velocity but not pressure. Although I suppose you put totalPressure at the inlet and zeroGradient at the outlet, which would be the correct b.c.

The first thing that comes to mind is: do you know for sure that the simulation is converged? Yours for some reason makes me think it is not fully converged yet. There are quite some threads here how to recognize a properly converged solution.

[chirag]

Hi LukeD,
I guess It has converged I have simulated it for really long time. My maximum velocity is 3000 m/s . so the particle sweeps through domain atleast 30 times.
Thanks in advance!
Chirag

[chirag]

Hi LukeD,
https://drive.google.com/file/d/0B5_...ew?usp=sharing
I have shared the results in google drive too along with all the files
Can you please go through them ?
Thanks in advance!
Chirag

[chirag]

Originally Posted by LukaD

Hey there,

Can you also give the boundary conditions for pressure? It looks like you put those of velocity but not pressure. Although I suppose you put totalPressure at the inlet and zeroGradient at the outlet, which would be the correct b.c.

The first thing that comes to mind is: do you know for sure that the simulation is converged? Yours for some reason makes me think it is not fully converged yet. There are quite some threads here how to recognize a properly converged solution.

Hi LukeD,
I have to maintain the chamber at T and P of(100 Pa, 300K). High energy gas (20000 Pa, 12000 k) enters the chamber and results in underexpanded flow through nozzle.
I am using totalPressure at inlet and waveTransmissive pressure B.C at outlet.
I have used similar B.C in fluent and results seem to be in very good agreement with exp and dsmc.
Thank you,
Chirag

[mkraposhin]

Hi,

can you add a small converging region at the inlet?

[mkraposhin]

Also, i would suggest to try running case with lower inlet/outlet pressure ratio

[chirag]

Hello Sir!
I will get it done and share it within few hours!
Thanks!
Chirag

[chirag]

Originally Posted by mkraposhin

Hi,

can you add a small converging region at the inlet?

I am trying to validate this paper by by selezneva et. al "Stationary supersonic plasma expansion: continuum fluid mechanics versus direct simulation Monte Carlo method" and I am trying to use same operating conditions (as I did in fluent for verifying the results). I will change the straight region with convergent one at inlet and fire it. I will also reduce the pressure ratio as advised.
Thanks!
http://iopscience.iop.org/article/10.../35/12/312/pdf

[mkraposhin]

Hi, chirag!
After some steps i was able to reach solution that is very close to fluent, see attached figures for Ma and T distribution along axis.

I done next steps:

1) I changed from upwind to Minmod scheme, i reverted to 3 PISO correctors. I think, that for this case it is possible to use less dissipative schemes, like vanLeer, vanAlbdada, filteredLinear2, etc

2) I switched off momentum predictor at system/fvSolution

3) Also, i changed this case for steady-state simulation (LTS), which is available for OpenFOAM 3.0.0, (see attached) tgz file. For steady-state simulation i used different rDeltaTSmoothingCoeff at different iterations:
0.001 for starting iterstions (0 - 5000)
0.01 for iterations (5001 - 30000)
0.02 for iterations > 30000

4) I changed inlet conditions for k, epsilon, U and p:
it is better to calculate k and epsilon at inlet using empirical formula depending on current velocity:
k:

    INLET
    {
        type            turbulentIntensityKineticEnergyInlet;
        intensity       0.01;
        value           uniform 0.1544303;
    }

epsilon:

    INLET
    {
        type            turbulentMixingLengthDissipationRateInlet;
        mixingLength    0.00036;
        value           uniform 67.93;
    }

mixing length is 9% of nozzle diameter

velocity must be calculated from current pressure with uniform profile, that's why i used pressureInletUniformVelocity:

    INLET
    {
        type            pressureInletUniformVelocity;
        value           uniform (0 0 0);
    }

And, pressure at inlet must be increased with time (or iterations), also it is important to set gamma to 0.9 ( if gamma < 1, it is not used) and psi to 'thermosi' - this will switch B.C. to more stable implementation for compressible flow:

    INLET
    {
        type            timeVaryingTotalPressure;
	value           uniform 100;
	p0              table
	21
	(
	    (0 100)
	    (10 200)
	    (20 300)
	    (30 400)
	    (50 500)
	    (60 600)
	    (70 700)
	    (80 800)
	    (90 900)
	    (100 1000)
	    (200 2000)
	    (300 2000)
	    (400 3000)
	    (500 3000)
	    (600 5000)
	    (700 5000)
	    (800 7500)
	    (900 7500)
	    (1000 10000)
	    (2000 10000)
	    (3000 10000)
	);
        U               U;
	phi             phi;
	psi		thermo:psi;
	gamma           0.9;
    }

This B.C. is not present in standard OpenFOAM distribution, to use it you have to download our library libcompressibleTools from github: https://github.com/unicfdlab/libcompressibleTools

On the outlet, i imposed special B.C. for presssure, which is also present only in libcompressibleTools - subsonicSupersonicPressureOutlet

    OUTLET
    {
        type            subsonicSupersonicPressureOutlet;
        value           $internalField;
        p0              $internalField;
        U               U;
        phi             phi;
        psi		thermo:psi;
        gamma           1.67;
        refValue	$internalField;
        refGradient	uniform 0;
        valueFraction	uniform 1;
    }

This B.C. swithes between fixedValue and zeroGradient depending on local Mach number.

As i understand from the article, authors used power law for viscousity and conductivity, for our case this leads to change in enthalpy and velocity diffusion coefficients more than 10 times, and it is also very important. So, i made a simple library, which implements linear approximation to power law properties (see attached file)

Also, authors told that they imposed 12000K and 10000Pa at the outlet of nozzle. It means that if you will set fixedValue for temperature at nozzle wall, you will get a strong heat flux, and you will loose a lot of thermal energy in the nozzle. So, i separated nozzle wall from other walls and used zeroGradient for temperature on nozzle surface.

At last, i made nozzle to be slightly convering, to prevent flow from reaching Ma > 1, see blockMeshDict in attached case

What next steps must be done to get closer to more physical solution:
1) We need to implement temperature-dependent Cp - because now it is constant, but it can vary significantly for such case
2) We need to implement power law for viscosity and thermal conductivity
3) We need to implement limiter for lower boundary of pressure - there is still a small region with p ~ 1 and this can lead to diverging solution. This can be done with fvOptions, but i need add fvOptions to current implementation of pisoCentralFoam
4) Also, i would propose to make outlet of nozzle diverging, this may also help to avoid low pressure region

[mkraposhin]

And a temperature field (case didn't converged to steady-state)

[chirag]

Hello Sir,
Thank you very much for looking into the problem.
There is no gradient of V & T observed in the results u have arrived as expected from experiments /simulation.

The author has inlet b.c as 20000 Pa and 12000K. She studied for outlet flow:
1) P = 100 Pa, 300k
2) P = 20 Pa, 300K

The pressure will sometimes go below 0 Pa and needs bounded as you said ( I had applied that before in Fluent). I had limited its value to 1e-8. The region where you have observed the low pressure is where fluent predicts (please see the attached images). We have good agreement with the experiments and the N-S simulation which she had performed. I had applied varying viscosity and specific heat in fluent for getting good agreement.
I am trying to get the same in openFoam.
Its really very kind of you to look into it.
Thanks,

[chirag]

hi,
I tried to install the new boundary condition library. I observed following error:
"cannot find -lfftw3"
Can you please gupdate the repository with the library file fftw3
Thanks!

[chirag]

Originally Posted by mkraposhin

Hi, chirag!
After some steps i was able to reach solution that is very close to fluent, see attached figures for Ma and T distribution along axis.

I done next steps:

1) I changed from upwind to Minmod scheme, i reverted to 3 PISO correctors. I think, that for this case it is possible to use less dissipative schemes, like vanLeer, vanAlbdada, filteredLinear2, etc

2) I switched off momentum predictor at system/fvSolution

3) Also, i changed this case for steady-state simulation (LTS), which is available for OpenFOAM 3.0.0, (see attached) tgz file. For steady-state simulation i used different rDeltaTSmoothingCoeff at different iterations:
0.001 for starting iterstions (0 - 5000)
0.01 for iterations (5001 - 30000)
0.02 for iterations > 30000

4) I changed inlet conditions for k, epsilon, U and p:
it is better to calculate k and epsilon at inlet using empirical formula depending on current velocity:
k:

    INLET
    {
        type            turbulentIntensityKineticEnergyInlet;
        intensity       0.01;
        value           uniform 0.1544303;
    }

epsilon:

    INLET
    {
        type            turbulentMixingLengthDissipationRateInlet;
        mixingLength    0.00036;
        value           uniform 67.93;
    }

mixing length is 9% of nozzle diameter

velocity must be calculated from current pressure with uniform profile, that's why i used pressureInletUniformVelocity:

    INLET
    {
        type            pressureInletUniformVelocity;
        value           uniform (0 0 0);
    }

And, pressure at inlet must be increased with time (or iterations), also it is important to set gamma to 0.9 ( if gamma < 1, it is not used) and psi to 'thermosi' - this will switch B.C. to more stable implementation for compressible flow:

    INLET
    {
        type            timeVaryingTotalPressure;
	value           uniform 100;
	p0              table
	21
	(
	    (0 100)
	    (10 200)
	    (20 300)
	    (30 400)
	    (50 500)
	    (60 600)
	    (70 700)
	    (80 800)
	    (90 900)
	    (100 1000)
	    (200 2000)
	    (300 2000)
	    (400 3000)
	    (500 3000)
	    (600 5000)
	    (700 5000)
	    (800 7500)
	    (900 7500)
	    (1000 10000)
	    (2000 10000)
	    (3000 10000)
	);
        U               U;
	phi             phi;
	psi		thermo:psi;
	gamma           0.9;
    }

This B.C. is not present in standard OpenFOAM distribution, to use it you have to download our library libcompressibleTools from github: https://github.com/unicfdlab/libcompressibleTools

On the outlet, i imposed special B.C. for presssure, which is also present only in libcompressibleTools - subsonicSupersonicPressureOutlet

    OUTLET
    {
        type            subsonicSupersonicPressureOutlet;
        value           $internalField;
        p0              $internalField;
        U               U;
        phi             phi;
        psi		thermo:psi;
        gamma           1.67;
        refValue	$internalField;
        refGradient	uniform 0;
        valueFraction	uniform 1;
    }

This B.C. swithes between fixedValue and zeroGradient depending on local Mach number.

As i understand from the article, authors used power law for viscousity and conductivity, for our case this leads to change in enthalpy and velocity diffusion coefficients more than 10 times, and it is also very important. So, i made a simple library, which implements linear approximation to power law properties (see attached file)

Also, authors told that they imposed 12000K and 10000Pa at the outlet of nozzle. It means that if you will set fixedValue for temperature at nozzle wall, you will get a strong heat flux, and you will loose a lot of thermal energy in the nozzle. So, i separated nozzle wall from other walls and used zeroGradient for temperature on nozzle surface.

At last, i made nozzle to be slightly convering, to prevent flow from reaching Ma > 1, see blockMeshDict in attached case

What next steps must be done to get closer to more physical solution:
1) We need to implement temperature-dependent Cp - because now it is constant, but it can vary significantly for such case
2) We need to implement power law for viscosity and thermal conductivity
3) We need to implement limiter for lower boundary of pressure - there is still a small region with p ~ 1 and this can lead to diverging solution. This can be done with fvOptions, but i need add fvOptions to current implementation of pisoCentralFoam
4) Also, i would propose to make outlet of nozzle diverging, this may also help to avoid low pressure region

hi,
I tried to install the new boundary condition library. I observed following error:
"cannot find -lfftw3"
Can you please gupdate the repository with the library file fftw3
Thanks!

[mkraposhin]

Hi,
You must compile library with script ./makeLib.sh

[mkraposhin]

If you will wait, i'm going to update library with thermophysical properties for your case

[mkraposhin]

Can you post here expression for Cp?

[chirag]

Hello Sir,
Please find the properties of argon needed for simulating this flow
Viscosity:
mu = mu0 (T/To)^0.72
where mu0 = 2.125e-5 kg/m.s
To = 273.11 K
specitic heat is function of temperature and linear variation has been calculated using the data provided in Thermal Plasmas Fundamentals and Applications
Volume 1 by Maher I. Boulos et. al

Cp = 6e-17 *T^5 - 4e-13*T^4 - 3e-9*T^3 +3e-5T^2 - 0.0597*T + 551.74
I am also attaching the table in text format if linear input of viscosity is required instead of expression
Thanks!!

[mkraposhin]

Hi,

I made it. It is interesting, what results you will obtain.

Don't forget to downoald rescent version of pisoCentralFoam and compressibleTools, case for simulation in the attachment

I used polynomial expression for Cp, because it is already done in OpenFOAM. In the future, it will be better to implement tabular properties, because it is much faster.

minimum value for pressure is set at system/fvOptions

[chirag]

Originally Posted by mkraposhin

Hi,

I made it. It is interesting, what results you will obtain.

Don't forget to downoald rescent version of pisoCentralFoam and compressibleTools, case for simulation in the attachment

I used polynomial expression for Cp, because it is already done in OpenFOAM. In the future, it will be better to implement tabular properties, because it is much faster.

minimum value for pressure is set at system/fvOptions

Hello Sir,
Thank you very much!
I will install 3.0.0 as I was using 2.3.1. I will share the results soon!
Regards,
Chirag