CFD Online Discussion Forums - Include friction to heat the fluid

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Include friction to heat the fluid (https://www.cfd-online.com/Forums/openfoam-programming-development/98155-include-friction-heat-fluid.html)

Include friction to heat the fluid

Hi all,

i am using OF since 2010 but i never programmed sth. in the code so i am a beginner and need help.

I wanna simulate an arc of a circle. The fluid is water and one wall is rotating with the following bc:

Code:

    

    wallOut

    {

        type            rotatingWallVelocity;

        origin          (0 0 0);

        axis            (0 0 1);

        omega         constant 40;

    }

the other wall is fixed (0 0 0).
Couse of friction the fluid should heated (maybe not very much but a little bit).

The question is how i realize that? I have no experience and no idea!
I am using buoyantBoussinesqPimpleFoam and looked into the temperature file. I think there is something missing that calculates the friction coused of pressure forces or something like that. Is that correct?

Maybe someone could help me.
Thx in advance.

Tobi

Hi Tobi,

I don't believe that that is a compressible solver. You would need to include some thermodynamic equation of state in order to do any pressure-temperature coupling.

Hi,
you are right. That is a incompressible Solver. Is it possible to implement it in that solver? I am a beginner in implementing things to OpenFOAM and its not easy for me to do that. :(

If you are more interested in heat due to friction you should be able to add a a source term for viscous work to the energy equation ("T" transport) that is in the solver you are currently using.

Quote:

Originally Posted by kmooney (Post 350223)

If you are more interested in heat due to friction you should be able to add a a source term for viscous work to the energy equation ("T" transport) that is in the solver you are currently using.

Okay thx. Do you know if there 's a solver in OpenFOAM which the heat due to firction is implemented?

Thx for your answers

Hello all,

i'm very interested in the creation of heat due to friction, because it's part of my diploma thesis, but i have no idea how to start in this case.

I'm using the simpleFoam solver and i already implemented the scalar temperatur like it is done in die "Add temperatur to icoFoam" - tutorial.
This works for the heat circulation, but not for creating heat with a rotating cylinder.

Now i hope to find an answer on how to add the term of viscous heat creation to my solver.

Perhaps someone could help Tobi and me in this case.
Tanks a lot.

Sebastian

If I understood your problem right, you might find this thread interesting:
http://www.cfd-online.com/Forums/ope...-enthalpy.html

Hi and thanks for the link.

Sounds really good, but i did't manage to integrate your terms in a simpleFoam solver.
I was trying, but i'm new to this, so perhaps you can help me a little bit?

Thanks a lot!
Sebastian

I think you have all the mathematics (and the implementation as well) in the other thread. However in the attachement you find an example for the implementation of viscous heating in pisoFoam, from there it is straightforward to implement the same changes into simpleFoam.
Let me know if you have problems compiling or managing the code.
Regards.

Thanks for your reply! I will try it right away!

Quote:

Originally Posted by cosimobianchini (Post 356371)

You may also find it useful to examine the rhoEpsilonEff() method available in the compressible turbulence models.

Hi all,

while compiling i got some little errors. I think it's easy for you to get the problem. But i have tried many variants and i think it's a problem of the included basic-files.

Here is, how i've integrated the code to simpleFoam:

simpleFoamT.C

Code:

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

========= |

\\ / F ield | OpenFOAM: The Open Source CFD Toolbox

\\ / O peration |

\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.

\\/ 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

simpleFoamT

Description

Steady-state solver for incompressible, turbulent flow with additional equation for temperatur

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

#include "fvCFD.H"

#include "singlePhaseTransportModel.H"

#include "RASModel.H"

#include "simpleControl.H"

 

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

int main(int argc, char *argv[])

{

#include "setRootCase.H"

#include "createTime.H"

#include "createMesh.H"

#include "createFields.H"

#include "initContinuityErrs.H"

#include "readTimeControls.H"

#include "CourantNo.H"

#include "setInitialDeltaT.H"

 

simpleControl simple(mesh);

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

 

Info<< "\nStarting time loop\n" << endl;

while (simple.loop())

{

Info<< "Time = " << runTime.timeName() << nl << endl;

 

p.storePrevIter();

 

#include "readTimeControls.H"

#include "simpleControl.H"

#include "CourantNo.H"

#include "setDeltaT.H"

 

// --- Pressure-velocity SIMPLE corrector

 

{

#include "UEqn.H"

#include "pEqn.H"

}

 

turbulence->correct();

 

#include "TEqn.H"

 

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

<< " ClockTime = " << runTime.elapsedClockTime() << " s"

<< nl << endl;

}

 

Info<< "End\n" << endl;

return 0;

}

 

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

UEqn.H:

Code:

// Momentum predictor

tmp<fvVectorMatrix> UEqn

(

fvm::div(phi, U)

+ turbulence->divDevReff(U)

);

UEqn().relax();

solve(UEqn() == -fvc::grad(p));

pEqn.H:

Code:

{

    p.boundaryField().updateCoeffs();

    volScalarField rAU(1.0/UEqn().A());

    U = rAU*UEqn().H();

    UEqn.clear();

    phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf();

    adjustPhi(phi, U, p);

    // Non-orthogonal pressure corrector loop

    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)

    {

        fvScalarMatrix pEqn

        (

            fvm::laplacian(rAU, p) == fvc::div(phi)

        );

        pEqn.setReference(pRefCell, pRefValue);

        pEqn.solve();

        if (nonOrth == simple.nNonOrthCorr())

        {

            phi -= pEqn.flux();

        }

    }

    #include "continuityErrs.H"

    // Explicitly relax pressure for momentum corrector

    p.relax();

    // Momentum corrector

    U -= rAU*fvc::grad(p);

    U.correctBoundaryConditions();

}

TEqn.H:

Code:

 DpDt = fvc::DDt(surfaceScalarField("phiU",phi),p);

 

        fvScalarMatrix TEqn

        (

              fvm::ddt(Cp,T)

            + fvm::div(phi*Cp,T)

            - fvm::laplacian((turbulence->nut()/Prt + turbulence->nu()/Pr)*Cp,T)

            ==

              DpDt

        );

        bool VH(mesh.solutionDict().found("VH"));

        if(VH)

        {

            TEqn +=  (turbulence->R() && fvc::grad(U));       

        }

        TEqn.solve();

createFields.H:

Code:

    Info<< "Reading field p\n" << endl;

    volScalarField p

    (

        IOobject

        (

            "p",

            runTime.timeName(),

            mesh,

            IOobject::MUST_READ,

            IOobject::AUTO_WRITE

        ),

        mesh

    );

    Info<< "Reading field U\n" << endl;

    volVectorField U

    (

        IOobject

        (

            "U",

            runTime.timeName(),

            mesh,

            IOobject::MUST_READ,

            IOobject::AUTO_WRITE

        ),

        mesh

    );

    #include "createPhi.H"

 

    label pRefCell = 0;

    scalar pRefValue = 0.0;

    setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);

    singlePhaseTransportModel laminarTransport(U, phi);

    autoPtr<incompressible::RASModel> turbulence

    (

        incompressible::RASModel::New(U, phi, laminarTransport)

    );

IOdictionary transportPropDict

    (

        IOobject

        (

            "transportProperties",

            U.time().constant(),

            U.db(),

            IOobject::MUST_READ,

            IOobject::NO_WRITE

        )

    );

 

    dimensionedScalar Cp(transportPropDict.lookup("Cp"));

    dimensionedScalar Pr(transportPropDict.lookup("Pr"));

    dimensionedScalar Prt(transportPropDict.lookup("Prt"));

 

    Info<< "Creating field DpDt\n" << endl;

    volScalarField DpDt =    fvc::DDt(surfaceScalarField("phiU", phi), p);

But while compiling following error appiers:

Quote:

Making dependency list for source file simpleFoamT.C
SOURCE=simpleFoamT.C ; g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -O3 -DNoRepository -ftemplate-depth-100 -I/na/share/openfoam/OpenFOAM-2.0.x/src/turbulenceModels -I/na/share/openfoam/OpenFOAM-2.0.x/src/turbulenceModels/incompressible/RAS/RASModel -I/na/share/openfoam/OpenFOAM-2.0.x/src/transportModels -I/na/share/openfoam/OpenFOAM-2.0.x/src/transportModels/incompressible/singlePhaseTransportModel -I/na/share/openfoam/OpenFOAM-2.0.x/src/finiteVolume/lnInclude -IlnInclude -I. -I/na/share/openfoam/OpenFOAM-2.0.x/src/OpenFOAM/lnInclude -I/na/share/openfoam/OpenFOAM-2.0.x/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linux64Gcc45DPOpt/simpleFoamT.o
In file included from simpleFoamT.C:78:0:
TEqn.H: In function ‘int main(int, char**)’:
TEqn.H:7:28: error: no matching function for call to ‘ddt(Foam::dimensionedScalar&, <unresolved overloaded function type>)’
TEqn.H:8:32: error: no matching function for call to ‘div(Foam::tmp<Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> >, <unresolved overloaded function type>)’
TEqn.H:9:80: error: no matching function for call to ‘laplacian(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >, <unresolved overloaded function type>)’
/na/share/openfoam/OpenFOAM-2.0.x/src/finiteVolume/lnInclude/readTimeControls.H:38:8: warning: unused variable ‘maxDeltaT’
make: *** [Make/linux64Gcc45DPOpt/simpleFoamT.o] Fehler 1

What do i have to change to get it running? Please help! :)

Thankfully and with best regards!

Sebastian

It seems you did not create the volScalarField T.
Furthermore consider that simpleFoam is a steady-state solver so all the ddt terms actually are not defined.

Quote:

Originally Posted by cosimobianchini (Post 356633)

It seems you did not create the volScalarField T.
Furthermore consider that simpleFoam is a steady-state solver so all the ddt terms actually are not defined.

Damn, you are right. I did net impliment the T volScalarFiled to creatField.H. Thanks a lot!

After creating T-volField the compiler got no errors.

But how can i define those ddt terms? Sorry i'm new to this.

Well the ddt (partial derivative in time) should be zero for a steady state problem, for the DDt (total derivative in time) instead you can a take a look at hEqn.H in rhoSimpleFoam. Basically the TEqn should look exactly the same apart from the unknown and the viscous term of course.
By the way if you are running with timeScheme steadyState the transient operators should be consistent anyway.
Regards

Hi again,

still have no good results...i have tried it with following TEqn.H file:

Code:

 //DpDt = fvc::DDt(surfaceScalarField("phiU",phi),p);

 

        fvScalarMatrix TEqn

        (

       ////fvm::ddt(Cp,T)

       //+ 

  //  fvm::div(phi*Cp,T)

         //- fvm::laplacian((turbulence->nut()/Prt + turbulence->nu()/Pr)*Cp,T)

 

          fvm::div(phi*Cp, T)

        - fvm::Sp(fvc::div(phi)*Cp, T)

        - fvm::laplacian((turbulence->nut()/Prt + turbulence->nu()/Pr)*Cp,T)

        ==

        - fvc::div(phi, 0.5*magSqr(U), "div(phi,k)")

 

     //DpDt

 

);

        bool VH(mesh.solutionDict().found("VH"));

        if(VH)

        {

            TEqn +=  (turbulence->R() && fvc::grad(U));       

        }

        TEqn.solve(); 

/*

  fvScalarMatrix TEqn

  (

      fvm::div(phi*Cp, T)

    - fvm::Sp(fvc::div(phi)*Cp, T)

    - fvm::laplacian((turbulence->alphaEff())*Cp, T)

   ==

     - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")

       );

*/

Everything compiled und also converged, but the results of temperatur are very very low. somewhere at the third position after decimal point.

Now i have startet a calculation with this TEqn.H:

Code:

 DpDt = fvc::DDt(surfaceScalarField("phiU",phi),p);

 

        fvScalarMatrix TEqn

        (

       ////fvm::ddt(Cp,T)

       //+ 

  //  fvm::div(phi*Cp,T)

         //- fvm::laplacian((turbulence->nut()/Prt + turbulence->nu()/Pr)*Cp,T)

 

          fvm::div(phi*Cp, T)

        - fvm::Sp(fvc::div(phi)*Cp, T)

        - fvm::laplacian((turbulence->nut()/Prt + turbulence->nu()/Pr)*Cp,T)

        ==

     //- fvc::div(phi, 0.5*magSqr(U), "div(phi,k)")

 

        DpDt

 

        );

        bool VH(mesh.solutionDict().found("VH"));

        if(VH)

        {

            TEqn +=  (turbulence->R() && fvc::grad(U));       

        }

        TEqn.solve(); 

/*

  fvScalarMatrix TEqn

  (

      fvm::div(phi*Cp, T)

    - fvm::Sp(fvc::div(phi)*Cp, T)

    - fvm::laplacian((turbulence->alphaEff())*Cp, T)

   ==

     - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")

       );

*/

I'm very excited which results i get now.
edit: After the first 100 iterations, the residuals were under 0.01.... but the Temperatur was exactly the some as it came from the inlet. so no heating has taken place.

What do you think about it? Perhaps you can help me one more time?? :)
Don't know how to return the favour, but thanks a lot!

Greetings

Are you sure the bool VH is set to true and you have sufficient velocity gradients inside your domain (i.e. a no slip wall)?
Viscous heating is generally low but not zero.

Well, in my case, i expect a temperature rise of about 0.7 kelvin.
This is, what i got from StarCCM+ for the same case.

Can you explain me, where i have to set the bool VH to true? Perhaps this is the answer to the question.

in my case i have a rotating cylinder (4000rpm) with nutkwallfunction for example in nut-file.
so i think i should have enought velocity gradients?

i also tried a total energy equation i've found in this forum. with that, i got a temperatur rise of about 1.4 K. Not to bad at all, but a little bit to high - perhaps this is because of the different implimentation of the turbulence models in starccm and openfoam.
i almost have the same problem when calculating the torque of the cylinder. in laminar case they are nearly the same. but not in turbulence.

i want to try your solver too, but unfortunately it wont work.

thanks a lot for your help!

greetings

bool VH(mesh.solutionDict().found("VH"));

tells you that the solver is checking for the word VH in the fvSolution dictionary.
You can write VH; in that dictionary or recompile the solver setting
bool VH(true);

Concerning the other energy equation, if you are mentioning the one discussed in this thread http://www.cfd-online.com/Forums/openfoam-solving/59468-full-energy-equation-enthalpy.html, it is actually very close to the one you are trying (once VH is set to true of course).

I'm not familiar with starCCM+, is it solving energy equation in terms of thermal (static) or total enthalpy. There is a sensible difference in case of high acceleration due to the different discretization of the source terms including viscous heating.

Hi,

you have been right. After including "VH;" in de fvSolution - file, heat is produced. Unfortunately the results are at about 2.8 K, and not the expected 0.7K like in StarCCM+.

Edit: I've started a new calculation in CCM+ with the same case and the same BC and now i have a temperatur rise of 2.08K. Still rising! hope it will stop at 2.8K ;)! Don't know why there is a difference between this and the first one.
But i will have a look at the first one, when this is ready.

At the Moment, i'm searching for the Equation, that is used by CCM+ and
I've found out, that it uses the Total Energy Equation.

Now, i'll try to figure out all the parts belonging to the Equation of CCM+, to get an overview of the difference between the two solvers.

I hope i can give you an overview soon.

Tanks for your help!

Greetings