CFD Online Discussion Forums - energy equation problem in porousSimpleFoam

Hi dear foamers
I wrote a TEqn for my solver.
the formula is something like thisA. Bejan Convection Heat Transfer,2nd edition, pages 571,575,576)

-dp/dx = [(mu/K) (u)] + [rho*beta*U^2] (darcy-forcheheimer)

and the energy equation is something like this

rho(fluid)*Cp(fluid)[sigma dT/dt + UdT/dx]=k d2T/dx2 + q'" + (mu/K) u^2

k and sigma are:

k= (porosity*kf) + (1- porosity)*ks

Sigma = [porosity(rho*cp)fluid+(1-porosity)*(rho)solid]/(rho*cp)f

for my case there is no heat generation and also I ignored the last term which shows viscous dissipation.
omitting Q"' and (mu/K) u^2, we have :

rho(fluid)*Cp(fluid)[sigma dT/dt + UdT/dx]=k d2T/dx2

dividing both sides of equation to rho(fluid)*Cp(fluid) we have:

[sigma dT/dt + UdT/dx]=alpha d2T/dx2

alpha= k / rho(fluid)*Cp(fluid)

I have all the properties of solid and fluid and porosity,so I can easily define alpha,sigma.I wrote my TEqn.H file like this

// Solve the Energy equation
{
fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(alpha, T)
);
TEqn.solve();

}

I modified createFields.H file either.I defined alpha and sigma in transportProperties file and ask the code to read them from there.

Info<< "Reading transportProperties\n" << endl;

IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);

dimensionedScalar sigma
(
transportProperties.lookup("sigma")
);
dimensionedScalar alpha
(
transportProperties.lookup("alpha")
);

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

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

and finally my code is something like this,I used include TEqn.H out of simple loop.

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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
porousSimpleFoam

Description
Steady-state solver for incompressible, turbulent flow with
implicit or explicit porosity treatment and support for multiple reference
frames (MRF)

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

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "RASModel.H"
#include "simpleControl.H"
#include "IOMRFZoneList.H"
#include "IOporosityModelList.H"
#include "fvIOoptionList.H"

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

int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"

simpleControl simple(mesh);

#include "createFields.H"
#include "createFvOptions.H"
#include "createZones.H"
#include "initContinuityErrs.H"

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

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

while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;

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

// ************************************************** *********************** //
I defined fvscheme and fvsolution file for T equation but I did not understand why this code give me wrong data.

would you please help me in this case,I'm really confused
every suggestion would be greatly appreciated.
if anyone has an experienced in this case please send me his/her email.
this error happened after running the code

Time = 0.014

GAMG: Solving for p, Initial residual = 0.964789, Final residual = 0.0417207, No Iterations 17
time step continuity errors : sum local = 2.99708e+12, global = 2.90211e+11, cumulative = 2.90211e+11
DILUPBiCG: Solving for T, Initial residual = 0.994907, Final residual = 0.0487249, No Iterations 22
ExecutionTime = 5.1 s ClockTime = 5 s

Time = 0.015

#0 Foam::error:rintStack(Foam::Ostream&) in "/opt/openfoam220/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler(int) in "/opt/openfoam220/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2 in "/lib/x86_64-linux-gnu/libc.so.6"
#3 at tensorField.C:0
#4
in "/home/vaio/OpenFOAM/vaio-2.2.0/platforms/linux64GccDPOpt/bin/newPorousSimpleFoam"
#5
in "/home/vaio/OpenFOAM/vaio-2.2.0/platforms/linux64GccDPOpt/bin/newPorousSimpleFoam"
#6 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#7
in "/home/vaio/OpenFOAM/vaio-2.2.0/platforms/linux64GccDPOpt/bin/newPorousSimpleFoam"
Floating point exception (core dumped)
vaio@vaio-SVE14122CVW:~/OpenFOAM/vaio-2.2.0/run/m$

best regards
Mohsen