[GUIDE] Switching from incompressible to compressible simulation

gabrielfelix · September 2, 2021, 09:11

Hi folks

I'm simulating a propeller in steady-state incompressible flow using simpleFoam with the same boundary conditions as the pimpleFoam or pimpleDyMFoam marine propeller tutorial. I wanted to run my case in compressible flow and I faced some problems to which I could not find a direct guide to overcome them. I saw that there are many posts concerning the same problems I had, and I decided to create this guide, now that I managed to succesfully run my case with rhoSimpleFoam.

I tried to apply the same file configurations of the rhoSimpleFoam aerofoilNACA0012 tutorial in my case, however I had to make other modifications to be able to run the case.

I'm not expert on this compressible configurations, so the solutions I will come up here might not be the most adequate ones, but rather the ones I managed to run my case with after researching in many and many threads.

1. Add the alphat and T files to 0/

Compressible solvers deals with the Navier-Stokes flow energy equation, which is related to the flow temperature. Thus, temperature T and alphat (which I don't know what means) must be added to the boundary conditions.

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
	location    "0";
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 298;

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    inlet
    {
        type            fixedValue;
        value           uniform 298;
	}
    outlet
    {
        type            zeroGradient;
    }
    wall
    {
        type            zeroGradient;
    }
	
	#includeEtc "caseDicts/setConstraintTypes"
}


// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      alphat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 0;

boundaryField
{
     #includeEtc "caseDicts/setConstraintTypes"

    wall
    {
        type            compressible::alphatWallFunction;
        value           uniform 0;
    }

    inlet
    {
        type            calculated;
        value           uniform 0;
    }
	
	outlet
    {
        type            calculated;
        value           uniform 0;
    }
}



// ************************************************************************* //

2. Add the density variable and its value to transportProperties file and create a thermophysicalProperties file

You must alter these files accordingly to your simulation altitude (i.e. rho = f(h) )

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

transportModel  Newtonian;

rho             1.225;

nu              1.8e-05;

// ************************************************************************* //

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     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 // air at room temperature (298 K)
{
    specie
    {
        molWeight   28.96;
    }
    thermodynamics
    {
        Cp          1004.5;
        Hf          0;
    }
    transport
    {
        mu          1.82e-05;
        Pr          0.71;
		Ts			116;
    }
}

// ************************************************************************* //

3. Add new variables to fvSchemes, fvSolution and create fvOptions file

You must add the new state variables of the flow related to energy and temperature.

I used the same configs from the aerofoilNACA0012 tutorial in fvSchemes

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;
	grad(k)         cellLimited Gauss linear 1;
	grad(omega)     cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;

    div(phi,U)      bounded Gauss linearUpwind limited;

    turbulence      bounded Gauss upwind;
    energy          bounded Gauss linearUpwind limited;

    div(phi,k)      $turbulence;
    div(phi,omega)  $turbulence;
	div(phi,epsilon)  $turbulence;

    div(phi,e)      $energy;
    div(phi,K)      $energy;
    div(phi,Ekp)    $energy;

    div(phid,p)     Gauss upwind;
    div((phi|interpolate(rho)),p)  bounded Gauss upwind;

    div(((rho*nuEff)*dev2(T(grad(U)))))    Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear limited corrected 0.33;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited corrected 0.33;
}

wallDist
{
    method meshWave;
}


// ************************************************************************* //

Add the energy equation solver and relaxation factors. It is also necessary to include rho in the relaxation factors. I used the same parameters as the tutorial.

I use potentialFlow to initialize the compressible simulation with an incompressible solution. nNonOrthogonalCorrectors in the SIMPLE section helps the solver when large pressure gradients appear in the flow.

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    "pcorr.*"
    {
        solver          GAMG;
        tolerance       1e-2;
        relTol          0;
        smoother        DICGaussSeidel;
        cacheAgglomeration no;
        maxIter         50;
    }

    p
    {
        $pcorr;
        tolerance       1e-5;
        relTol          0.01;
    }

    pFinal
    {
        $p;
        tolerance       1e-6;
        relTol          0;
    }

    "(U|k|epsilon|e)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0.1;
    }

    "(U|k|epsilon|e)Final"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0;
    }
}

SIMPLE
{
    //correctPhi          no;
    //nOuterCorrectors    2;
    //nCorrectors         1;
    nNonOrthogonalCorrectors 1;
	
	// Set up residual controls. Simulation stops when residual target is reached
	residualControl
	{
		p					1e-6;
		//U					1e-4;
		//"(k|epsilon)"		1e-4;
	}
}

potentialFlow
{
    nNonOrthogonalCorrectors 1;
}

relaxationFactors
{
	fields
	{
		p 0.3 ;
		rho 0.01 ;
	}
	equations
	{
		U 0.7 ;
		e 0.7 ;
		k 0.7 ;
		omega 0.7 ;
	}
}

cache
{
    grad(U);
}


// ************************************************************************* //

Create an fvOptions file to limit the temperature range and avoid getting non-physical values.

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

limitT
{
    type       limitTemperature;
    min        101;
    max        1000;
    selectionMode all;
}

//************************************************************************** //

4. Change the forces or coefficients files if you use them

These files determines how the output forces and coefficients are printed by the solver. In incompressible flows, the rho = inf and the pressure is given by P = p/rho. Now that you have a finite density you must attribute its value to the files.

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/

forces
{
    type          forces;

    libs          (forces);

    writeControl  timeStep;
    timeInterval  1;

    log           yes;

    patches       ("propeller.*");
    rho	   rho;
    log           true;


    CofR          (0 0 0);    // Rotation around centre line of propeller
    pitchAxis     (0 1 0);
}


// ************************************************************************* //

5. Change the dimension of pressure its value in 0/p file

Quoting the user snak:

Quote:

In an incompressible solver, N-S equation is divided by a uniform density rho. This causes (1) the dimension of pressure of [0 2 -2 0 0 0 0] and (2) the kinematic viscosity nu in laplacian term. In an incompressible solver, pressure is assumed to be relative. The atmosphere will be 0 usually.

In a compressible solver, N-S equation is not divided by density. So, the dimension of pressure is [1 -1 -2 0 0 0 0] as usual. The dinamic viscosity mu appears in laplacian term. In a compressible solver, The absolute pressure is must be provided in p file because the value of pressure will be used to calculate other physical properies. The atmosphere will be 1e5 usually.

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 101325;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 101325;
    }

    wall
    {
        type            zeroGradient;
    }
	
	#includeEtc "caseDicts/setConstraintTypes"
}


// ************************************************************************* //

snak · September 3, 2021, 13:39

Hi Gabriel,

Thank you for summarizing the procedure.

Some important points are missing in this guide. The dimension and type of pressure and the physical property of viscosity depend on the type of solver (compressible or incompressible). We have to adjust these as switch the solvers.

In an incompressible solver, N-S equation is divided by a uniform density rho. This causes (1) the dimension of pressure of [0 2 -2 0 0 0 0] and (2) the kinematic viscosity nu in laplacian term. In an incompressible solver, pressure is assumed to be relative. The atmosphere will be 0 usually.

In a compressible solver, N-S equation is not divided by density. So, the dimension of pressure is [1 -1 -2 0 0 0 0] as usual. The dinamic viscosity mu appears in laplacian term. In a compressible solver, The absolute pressure is must be provided in p file because the value of pressure will be used to calculate other physical properies. The atmosphere will be 1e5 usually.

gabrielfelix · September 3, 2021, 13:45

Quote:

Originally Posted by snak

Hi Gabriel,

Thank you for summarizing the procedure.

Some important points are missing in this guide. The dimension and type of pressure and the physical property of viscosity depend on the type of solver (compressible or incompressible). We have to adjust these as switch the solvers.

In an incompressible solver, N-S equation is divided by a uniform density rho. This causes (1) the dimension of pressure of [0 2 -2 0 0 0 0] and (2) the kinematic viscosity nu in laplacian term. In an incompressible solver, pressure is assumed to be relative. The atmosphere will be 0 usually.

In a compressible solver, N-S equation is not divided by density. So, the dimension of pressure is [1 -1 -2 0 0 0 0] as usual. The dinamic viscosity mu appears in laplacian term. In a compressible solver, The absolute pressure is must be provided in p file because the value of pressure will be used to calculate other physical properies. The atmosphere will be 1e5 usually.

Thanks man! I forgot the pressure file changes. I'll add them to the guide. There are so many changes you gotta made that I took some of them for granted.

September 3, 2021, 13:39		#2
snak Senior Member shinji nakagawa Join Date: Mar 2009 Location: Japan Posts: 113 Blog Entries: 1 Rep Power: 18	Hi Gabriel, Thank you for summarizing the procedure. Some important points are missing in this guide. The dimension and type of pressure and the physical property of viscosity depend on the type of solver (compressible or incompressible). We have to adjust these as switch the solvers. In an incompressible solver, N-S equation is divided by a uniform density rho. This causes (1) the dimension of pressure of [0 2 -2 0 0 0 0] and (2) the kinematic viscosity nu in laplacian term. In an incompressible solver, pressure is assumed to be relative. The atmosphere will be 0 usually. In a compressible solver, N-S equation is not divided by density. So, the dimension of pressure is [1 -1 -2 0 0 0 0] as usual. The dinamic viscosity mu appears in laplacian term. In a compressible solver, The absolute pressure is must be provided in p file because the value of pressure will be used to calculate other physical properies. The atmosphere will be 1e5 usually. gabrielfelix and Maxwell-Smoluchowski like this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
compressible vs incompressible flow	Nick R	Main CFD Forum	1	March 15, 2017 12:29
Compressible vs incompressible SGS model!	zhangyan	Main CFD Forum	10	June 4, 2016 05:39
Modify compressible solver to incompressible	mingzhao	OpenFOAM Programming & Development	0	September 18, 2015 18:46
Modify compressible solver to incompressible	mingzhao	OpenFOAM Running, Solving & CFD	0	September 18, 2015 13:50
compressible liquid simulation	kaiser	Main CFD Forum	0	October 13, 2002 05:32