
[Sponsors] 
March 4, 2012, 16:36 
Include friction to heat the fluid

#1 
Senior Member
Tobias Holzmann
Join Date: Oct 2010
Location: Leoben (Austria)
Posts: 1,113
Blog Entries: 6
Rep Power: 19 
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; } 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 

March 13, 2012, 14:15 

#2 
Senior Member
Kyle Mooney
Join Date: Jul 2009
Location: Amherst, MA USA  San Diego, CA USA
Posts: 290
Rep Power: 9 
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 pressuretemperature coupling. 

March 19, 2012, 08:18 

#3 
Senior Member
Tobias Holzmann
Join Date: Oct 2010
Location: Leoben (Austria)
Posts: 1,113
Blog Entries: 6
Rep Power: 19 
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. 

March 19, 2012, 11:18 

#4 
Senior Member
Kyle Mooney
Join Date: Jul 2009
Location: Amherst, MA USA  San Diego, CA USA
Posts: 290
Rep Power: 9 
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.


March 19, 2012, 13:05 

#5  
Senior Member
Tobias Holzmann
Join Date: Oct 2010
Location: Leoben (Austria)
Posts: 1,113
Blog Entries: 6
Rep Power: 19 
Quote:
Thx for your answers 

April 16, 2012, 08:27 

#6 
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 5 
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 

April 19, 2012, 04:23 

#7 
Member

If I understood your problem right, you might find this thread interesting:
Full energy equation for enthalpy
__________________
Cosimo Bianchini Energy Engineering Department "S. Stecco" University of Florence Via di S.Marta, 3 50139 Florence  ITALY Tel: +39 055 4796575 Fax: +39 055 4796342 Mob: +39 320 9460153 email: cosimo.bianchini@htc.de.unifi.it URL: www.htc.de.unifi.it 

April 23, 2012, 05:10 

#8 
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 5 
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 

April 23, 2012, 07:12 

#9 
Member

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.
__________________
Cosimo Bianchini Energy Engineering Department "S. Stecco" University of Florence Via di S.Marta, 3 50139 Florence  ITALY Tel: +39 055 4796575 Fax: +39 055 4796342 Mob: +39 320 9460153 email: cosimo.bianchini@htc.de.unifi.it URL: www.htc.de.unifi.it 

April 23, 2012, 07:37 

#10 
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 5 
Thanks for your reply! I will try it right away!


April 24, 2012, 02:46 

#11  
Senior Member
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 778
Rep Power: 18 
Quote:
You may also find it useful to examine the rhoEpsilonEff() method available in the compressible turbulence models. 

April 24, 2012, 05:41 

#12  
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 5 
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 basicfiles. 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) 20042011 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 Steadystate 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" //  Pressurevelocity 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)); 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); // Nonorthogonal 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(); } 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(); 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); Quote:
Thankfully and with best regards! Sebastian 

April 24, 2012, 05:58 

#13 
Member

It seems you did not create the volScalarField T.
Furthermore consider that simpleFoam is a steadystate solver so all the ddt terms actually are not defined.
__________________
Cosimo Bianchini Energy Engineering Department "S. Stecco" University of Florence Via di S.Marta, 3 50139 Florence  ITALY Tel: +39 055 4796575 Fax: +39 055 4796342 Mob: +39 320 9460153 email: cosimo.bianchini@htc.de.unifi.it URL: www.htc.de.unifi.it 

April 24, 2012, 06:16 

#14  
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 5 
Quote:
Damn, you are right. I did net impliment the T volScalarFiled to creatField.H. Thanks a lot! After creating TvolField the compiler got no errors. But how can i define those ddt terms? Sorry i'm new to this. 

April 24, 2012, 06:28 

#15 
Member

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
__________________
Cosimo Bianchini Energy Engineering Department "S. Stecco" University of Florence Via di S.Marta, 3 50139 Florence  ITALY Tel: +39 055 4796575 Fax: +39 055 4796342 Mob: +39 320 9460153 email: cosimo.bianchini@htc.de.unifi.it URL: www.htc.de.unifi.it 

May 2, 2012, 05:19 

#16 
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 5 
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)") ); */ 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)") ); */ 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 Last edited by SiCaRiUs; May 2, 2012 at 06:47. 

May 3, 2012, 03:45 

#17 
Member

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.
__________________
Cosimo Bianchini Energy Engineering Department "S. Stecco" University of Florence Via di S.Marta, 3 50139 Florence  ITALY Tel: +39 055 4796575 Fax: +39 055 4796342 Mob: +39 320 9460153 email: cosimo.bianchini@htc.de.unifi.it URL: www.htc.de.unifi.it 

May 3, 2012, 08:54 

#18 
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 5 
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 nutfile. 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 

May 3, 2012, 11:25 

#19 
Member

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.cfdonline.com/Forums/openfoamsolving/59468fullenergyequationenthalpy.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.
__________________
Cosimo Bianchini Energy Engineering Department "S. Stecco" University of Florence Via di S.Marta, 3 50139 Florence  ITALY Tel: +39 055 4796575 Fax: +39 055 4796342 Mob: +39 320 9460153 email: cosimo.bianchini@htc.de.unifi.it URL: www.htc.de.unifi.it 

May 7, 2012, 06:28 

#20 
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 5 
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 Last edited by SiCaRiUs; May 7, 2012 at 08:24. 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Water subcooled boiling  Attesz  CFX  7  January 5, 2013 04:32 
heat conducting solid within a fluid medium?  sieginc.  STARCCM+  1  February 19, 2012 12:48 
Solid / Fluid Heat Transfer  Koranten  FLUENT  3  March 19, 2011 08:21 
Simulation of a single bubble with a VOFmethod  Suzzn  CFX  18  October 2, 2009 04:18 
Conjugate Heat Transfer between fluid and solid  Li Yang  Main CFD Forum  8  March 27, 2004 12:05 