[hiuluom]

Hi everybody,

When I tried wallHeatFlux the buoyantSimpleFoam OF 2.1.1 in hotroom tutorial it was no problem with thermotype:

thermoType      hPsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>;

pRef            100000;

mixture
{
    specie
    {
        nMoles          1;
        molWeight       28.9;
    }
    thermodynamics
    {
        Cp              1000;
        Hf              0;
    }
    transport
    {
        mu              1.8e-05;
        Pr              0.7;
    }
}

But in sonicFoam with prism tutorial has difference thermo:

thermoType      ePsiThermo<pureMixture<constTransport<specieThermo<eConstThermo<perfectGas>>>>>;

mixture
{
    specie
    {
        nMoles          1;
        molWeight       28.9;
    }
    thermodynamics
    {
        Cv              717.5;
        Hf              0;
    }
    transport
    {
        mu              1.8e-05;
        Pr              0.7;
    }
}

I used the wallHeatFluxRho and wallHeatflux (default) which was provide in forum, OF 2.1.1 and had been a problem:

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 2.1.1-221db2718bbb
Exec   : wallHeatFlux -latestTime
Date   : Jul 25 2015
Time   : 06:53:06
Host   : "CompEng"
PID    : 2979
Case   : /home/huynh/OpenFOAM/huynh-2.1.1/run/prism
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time
Create mesh for time = 0.0003
Time = 0.0003
Selecting thermodynamics package ePsiThermo<pureMixture<constTransport<specieThermo<eConstThermo<perfectGas>>>>>

--> FOAM FATAL ERROR: 
Not implemented
    From function basicThermo::h()
    in file basicThermo/basicThermo.C at line 260.
FOAM aborting
#0  Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#1  Foam::error::abort() in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#2  Foam::basicThermo::h() in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libbasicThermophysicalModels.so"
#3  
 in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/wallHeatFlux"
#4  __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#5  
 in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/wallHeatFlux"
Aborted

therefore, I copy the thermotype from buoyantSimpleFoam to sonicFoam calculation wallHeatFlux again. It was run but the result equal to zero.

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 2.1.1-221db2718bbb
Exec   : wallHeatFluxRho -latestTime
Date   : Jul 25 2015
Time   : 22:05:09
Host   : "CompEng"
PID    : 4590
Case   : /home/huynh/OpenFOAM/huynh-2.1.1/run/prism
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0.0003

Time = 0.0003
Selecting thermodynamics package hPsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>
Reading/calculating face flux field phi

Selecting RAS turbulence model kEpsilon
kEpsilonCoeffs
{
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    C3              -0.33;
    sigmak          1;
    sigmaEps        1.3;
    Prt             1;
}


Wall heat fluxes [W]
prismWall 0

End

I saw the source code wallHeatFlux which was compute by gradient of enthalpy h corresponding hPsiThermo. Otherwise, the ePsiThermo with internal energy e was not described in equation wallHeatFlux. I replace the enthalpy h by internal energy e and compile it unsuccessful.

fvc::interpolate(RASModel->alphaEff())*fvc::snGrad(h)

I would like to calculate wallHeatFlux for sonicFoam but it cannot. Could anybody can show me the solve for this, please?

Thank you so much

[hiuluom]

Hi everybody,

I modified the wallHeatFlux that can be calculated for ePsithermo. Here is the code:
WallHeatFlux.C

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    wallHeatFlux

Description
    Calculates and writes the heat flux for all patches as the boundary field
    of a volScalarField and also prints the integrated flux for all wall
    patches.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "hCombustionThermo.H"
#include "basicPsiThermo.H"
#include "RASModel.H"
#include "turbulenceModel.H"
#include "wallFvPatch.H"

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

int main(int argc, char *argv[])
{
    timeSelector::addOptions();
#   include "addRegionOption.H"
    #include "setRootCase.H"
    #include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
    #include "createNamedMesh.H"

    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;
        mesh.readUpdate();

        #include "createFields.H"

        surfaceScalarField heatFlux
        (
            fvc::interpolate(turbulence->alphaEff())*fvc::snGrad(e)
        );

        const surfaceScalarField::GeometricBoundaryField& patchHeatFlux =
            heatFlux.boundaryField();

        Info<< "\nWall heat fluxes [W]" << endl;
        forAll(patchHeatFlux, patchi)
        {
            if (isA<wallFvPatch>(mesh.boundary()[patchi]))
            {
                Info<< mesh.boundary()[patchi].name()
                    << " "
                    << gSum
                       (
                           mesh.magSf().boundaryField()[patchi]
                          *patchHeatFlux[patchi]
                       )
                    << endl;
            }
        }
        Info<< endl;

        volScalarField wallHeatFlux
        (
            IOobject
            (
                "wallHeatFlux",
                runTime.timeName(),
                mesh
            ),
            mesh,
            dimensionedScalar("wallHeatFlux", heatFlux.dimensions(), 0.0)
        );

        forAll(wallHeatFlux.boundaryField(), patchi)
        {
            wallHeatFlux.boundaryField()[patchi] = patchHeatFlux[patchi];
        }

        wallHeatFlux.write();
    }

    Info<< "End" << endl;

    return 0;
}

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

And createFields.H

    Info<< "Reading thermophysical properties\n" << endl;

    autoPtr<basicPsiThermo> pThermo
    (
        basicPsiThermo::New(mesh)
    );
    basicPsiThermo& thermo = pThermo();

    //volScalarField& p = thermo.p();
    volScalarField& e = thermo.e();
    //const volScalarField& psi = thermo.psi();

    volScalarField rho
    (
        IOobject
        (
            "rho",
            runTime.timeName(),
            mesh
        ),
        thermo.rho()
    );

    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    #include "compressibleCreatePhi.H"


    Info<< "Creating turbulence model\n" << endl;
    autoPtr<compressible::turbulenceModel> turbulence
    (
        compressible::turbulenceModel::New
        (
            rho,
            U,
            phi,
            thermo
        )
    );

Then the result with my nozzle tube in example ran sonicFoam

Create time

Create mesh for time = 0.0090052

Time = 0.0090052
Reading thermophysical properties

Selecting thermodynamics package ePsiThermo<pureMixture<constTransport<specieThermo<eConstThermo<perfectGas>>>>>
Reading field U

Reading/calculating face flux field phi

Creating turbulence model

Selecting turbulence model type RASModel
Selecting RAS turbulence model kOmegaSST
kOmegaSSTCoeffs
{
    alphaK1         0.85034;
    alphaK2         1;
    alphaOmega1     0.5;
    alphaOmega2     0.85616;
    beta1           0.075;
    beta2           0.0828;
    betaStar        0.09;
    gamma1          0.05532;
    gamma2          0.4403;
    a1              0.31;
    c1              10;
    Prt             1;
}


Wall heat fluxes [W]
bottom_heat -104748

But it has been difference with wallHeatFlux default in OF and replace epsiThermo by hpsiThermo.

Create time

Create mesh for time = 0.0090052

Time = 0.0090052
Selecting thermodynamics package hPsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>
Reading/calculating face flux field phi

Selecting RAS turbulence model kOmegaSST
kOmegaSSTCoeffs
{
    alphaK1         0.85034;
    alphaK2         1;
    alphaOmega1     0.5;
    alphaOmega2     0.85616;
    beta1           0.075;
    beta2           0.0828;
    betaStar        0.09;
    gamma1          0.05532;
    gamma2          0.4403;
    a1              0.31;
    c1              10;
    Prt             1;
}


Wall heat fluxes [W]
bottom_heat -146611

End

[hiuluom]

Dear all,

I know the difference result as I mention above. That is wrong at formula calculate heat flux in *.C. I check with paraview this formula after computing heat flux with ePsiThermo

wallHeatFlux*(Cp/Cv)

It will the same result when calculate heat flux with hPsiThermo. And then I modified the code wallHeatFlux.C by

//molecular weight number
		// scalar W(thermo.lookup("molWeight"));

		//Heat capacity
		// scalar Cv(thermo.lookup("Cv"));

        surfaceScalarField heatFlux
        (
            fvc::interpolate(turbulence->alphaEff())*(W + 8314/(Cv*W))*fvc::snGrad(e)
        );

I have error because OF cannot look up molWeight and Cv in thermoPhysicalProperties file

Making dependency list for source file wallHeatFluxCompressible.C
SOURCE=wallHeatFluxCompressible.C ;  g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3  -DNoRepository -ftemplate-depth-100 -I.. -I/opt/OpenFOAM/OpenFOAM-2.1.1/src/turbulenceModels     -I/opt/OpenFOAM/OpenFOAM-2.1.1/src/turbulenceModels/compressible/turbulenceModel -I/opt/OpenFOAM/OpenFOAM-2.1.1/src/thermophysicalModels/specie/lnInclude -I/opt/OpenFOAM/OpenFOAM-2.1.1/src/thermophysicalModels/reactionThermo/lnInclude -I/opt/OpenFOAM/OpenFOAM-2.1.1/src/thermophysicalModels/basic/lnInclude -I/opt/OpenFOAM/OpenFOAM-2.1.1/src/finiteVolume/lnInclude -IlnInclude -I. -I/opt/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-2.1.1/src/OSspecific/POSIX/lnInclude   -fPIC -c $SOURCE -o Make/linux64GccDPOpt/wallHeatFluxCompressible.o
In file included from /opt/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/Field.H:360:0,
                 from /opt/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/labelField.H:39,
                 from /opt/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/primitiveFields.H:37,
                 from /opt/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/pointField.H:36,
                 from /opt/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/edge.H:40,
                 from /opt/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/edgeList.H:32,
                 from /opt/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/primitiveMesh.H:57,
                 from /opt/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/polyMesh.H:44,
                 from /opt/OpenFOAM/OpenFOAM-2.1.1/src/finiteVolume/lnInclude/fvMesh.H:50,
                 from /opt/OpenFOAM/OpenFOAM-2.1.1/src/finiteVolume/lnInclude/fvCFD.H:7,
                 from wallHeatFluxCompressible.C:34:
/opt/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/Field.C: In member function ‘void Foam::Field<Type>::operator=(const Foam::VectorSpace<Form, Cmpt, nCmpt>&)’:
/opt/OpenFOAM/OpenFOAM-2.1.1/src/OpenFOAM/lnInclude/Field.C:680:42: warning: typedef ‘VSType’ locally defined but not used [-Wunused-local-typedefs]
     typedef VectorSpace<Form,Cmpt,nCmpt> VSType;
                                          ^
In file included from /opt/OpenFOAM/OpenFOAM-2.1.1/src/thermophysicalModels/specie/lnInclude/eConstThermo.H:214:0,
                 from wallHeatFluxCompressible.C:36:
/opt/OpenFOAM/OpenFOAM-2.1.1/src/thermophysicalModels/specie/lnInclude/eConstThermoI.H: In member function ‘Foam::scalar Foam::eConstThermo<EquationOfState>::cp(Foam::scalar) const’:
/opt/OpenFOAM/OpenFOAM-2.1.1/src/thermophysicalModels/specie/lnInclude/eConstThermoI.H:108:28: error: ‘specie’ has not been declared
     return Cv_*this->W() + specie::RR;
                            ^
wallHeatFluxCompressible.C: In function ‘int main(int, char**)’:
wallHeatFluxCompressible.C:61:24: error: missing template arguments before ‘.’ token
   scalar W(eConstThermo.lookup("molWeight"));
                        ^
wallHeatFluxCompressible.C:64:25: error: missing template arguments before ‘.’ token
   scalar Cv(eConstThermo.lookup("Cv"));
                         ^
wallHeatFluxCompressible.C:61:10: warning: unused variable ‘W’ [-Wunused-variable]
   scalar W(eConstThermo.lookup("molWeight"));
          ^
wallHeatFluxCompressible.C:64:10: warning: unused variable ‘Cv’ [-Wunused-variable]
   scalar Cv(eConstThermo.lookup("Cv"));
          ^
make: *** [Make/linux64GccDPOpt/wallHeatFluxCompressible.o] Error 1

Someone can help me how to call Cv and molWeight in *.C

Thank you,
Thanh

[hiuluom]

Finally, I compiled successfully wall heat flux for "ePsiThermo" thermal type for sonicFoam or rhoCentralFoam. Here is the code:

for "createFields.h"

Info<< "Reading thermophysical properties\n" << endl;

    autoPtr<basicPsiThermo> pThermo
    (
        basicPsiThermo::New(mesh)
    );
    basicPsiThermo& thermo = pThermo();
	
    //volScalarField& p = thermo.p();
    const volScalarField& e = thermo.e();
	
	const objectRegistry& db = mesh.thisDb();
	
	//const IOdictionary& thermoProps = db.lookupObject<IOdictionary>("thermophysicalProperties");
	
	const dictionary& thermoProps =
        db.lookupObject<IOdictionary>("thermophysicalProperties");
	
	dictionary mixture
		(
			thermoProps.subDict("mixture")
		);
	dictionary specie
		(
			mixture.subDict("specie")
		);
	dictionary thermodynamics
		(
			mixture.subDict("thermodynamics")
		);
	const scalar Cv =
        readScalar(thermodynamics.lookup("Cv"));
		
	const scalar W =
        readScalar(specie.lookup("molWeight"));
	
	//double gamma = (Cv + 8314/(Cv*W));
		
	//const volScalarField& Cv = thermo.Cv();
	//const volScalarField& h = thermo.h();
    //const volScalarField& psi = thermo.psi();

    volScalarField rho
    (
        IOobject
        (
            "rho",
            runTime.timeName(),
            mesh
        ),
        thermo.rho()
    );

    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    #include "compressibleCreatePhi.H"


    Info<< "Creating turbulence model\n" << endl;
    autoPtr<compressible::turbulenceModel> turbulence
    (
        compressible::turbulenceModel::New
        (
            rho,
            U,
            phi,
            thermo
        )
    );

for "wallHeatFLux.C"

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    wallHeatFlux

Description
    Calculates and writes the heat flux for all patches as the boundary field
    of a volScalarField and also prints the integrated flux for all wall
    patches.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "hCombustionThermo.H"
#include "basicPsiThermo.H"
#include "turbulenceModel.H"
#include "wallFvPatch.H"


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

int main(int argc, char *argv[])
{
    timeSelector::addOptions();
#   include "addRegionOption.H"
    #include "setRootCase.H"
    #include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
    #include "createNamedMesh.H"
	
    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;
        mesh.readUpdate();

        #include "createFields.H"
				
        surfaceScalarField heatFlux
        (
            fvc::interpolate(turbulence->alphaEff())*(1 + 8314/(Cv*W))*fvc::snGrad(e)
        );

        const surfaceScalarField::GeometricBoundaryField& patchHeatFlux =
            heatFlux.boundaryField();

        Info<< "\nWall heat fluxes [W]" << endl;
        forAll(patchHeatFlux, patchi)
        {
            if (isA<wallFvPatch>(mesh.boundary()[patchi]))
            {
               Info<< mesh.boundary()[patchi].name()
                    << ": Total "
                    << sum
                       (
                           mesh.magSf().boundaryField()[patchi]
                          *patchHeatFlux[patchi]
                       )
                    << " [W] over "
                    << sum
                       (
                           mesh.magSf().boundaryField()[patchi]
                       )
                    << " [m2] ("
                    << sum
                       (
                           mesh.magSf().boundaryField()[patchi]
                          *patchHeatFlux[patchi]
                       )/
                       sum 
                       (
                           mesh.magSf().boundaryField()[patchi]
                       )
                    << " [W/m2])"
                    << endl;
            }
        }
        Info<< endl;

        volScalarField wallHeatFlux
        (
            IOobject
            (
                "wallHeatFlux",
                runTime.timeName(),
                mesh
            ),
            mesh,
            dimensionedScalar("wallHeatFlux", heatFlux.dimensions(), 0.0)
        );

        forAll(wallHeatFlux.boundaryField(), patchi)
        {
            wallHeatFlux.boundaryField()[patchi] = patchHeatFlux[patchi];
        }

        wallHeatFlux.write();
    }

    Info<< "End" << endl;

    return 0;
}

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

This result have not been difference with "hPsiThermo" thermal type in section 2 above

Create time

Create mesh for time = 0.0090052

Time = 0.0090052
Reading thermophysical properties

Selecting thermodynamics package ePsiThermo<pureMixture<constTransport<specieThermo<eConstThermo<perfectGas>>>>>
Reading field U

Reading/calculating face flux field phi

Creating turbulence model

Selecting turbulence model type RASModel
Selecting RAS turbulence model kOmegaSST
kOmegaSSTCoeffs
{
    alphaK1         0.85034;
    alphaK2         1;
    alphaOmega1     0.5;
    alphaOmega2     0.85616;
    beta1           0.075;
    beta2           0.0828;
    betaStar        0.09;
    gamma1          0.05532;
    gamma2          0.4403;
    a1              0.31;
    c1              10;
    Prt             1;
}


Wall heat fluxes [W]
bottom_heat: Total -146673 [W] over 0.00909676 [m2] (-1.61237e+07 [W/m2])

End

I got much more time to discover this problem because I must do other simulations.

Best regards,
Thanh