CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Conjugated solver for momentum and energy equations system (http://www.cfd-online.com/Forums/openfoam-solving/62746-conjugated-solver-momentum-energy-equations-system.html)

makaveli_lcf March 18, 2009 11:46

Conjugated solver for momentum and energy equations system
 
Hello colleagues!

First of all I would like to express my gratitude to Professor Hrvoje Jasak for developing such powerful CFD tool like OpenFOAM!

My question concerns solution of momentum-energy equations system. In my study case they are conjugated due to drag force source terms for momentum equation. I'm developing transient solver for this task. So it would be rather useful to get an opinion of OpenFOAM experienced users:

1. Is there any built-in tool to solve system of conjugated PDE system?
2. Are there some methods of non-linear PDE solution implemented in OpenFOAM?

Suppose we have system:

dU/dt + F(U,T) = 0
dT/dt + G(U,T) = 0

(d/dt - is partial derivative, specific heat supposed to be constant so we get from enthalpy to temperature)

Adding energy equation (like in http://openfoamwiki.net/index.php/Ho...e_to_icoFoam):

/************************************************** ***********************/
for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;

# include "readPISOControls.H"
# include "CourantNo.H"

/************************************************/
/* MOMENTUM EQUATION */
/************************************************/

// Momentum equation system
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
- fvm::Sp(S_D, U)
);

// Solution of momentum equation
solve(UEqn == - fvc::grad(p) + S_B);

// --- PISO loop

for (int corr=0; corr<nCorr; corr++)
{
volScalarField rUA = 1.0/UEqn.A();

U = rUA*UEqn.H();
// Interpolation using central differencing
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi);

adjustPhi(phi, U, p);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi)
);

pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();

if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}

# include "continuityErrs.H"

// Coorection of flow field
U -= rUA*(fvc::grad(p) - S_B);
U.correctBoundaryConditions();
}

/************************************************/
/* ENERGY EQUATION */
/************************************************/

// Energy equation system
fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(alpha, T)
);

// Solve enqrgy equation
TEqn.solve();

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

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

But now due to S_D (drag force source term) and S_B (buoyancy force source term), which are calculated each time step and both are based on temperature distribution, two equations for U and T are coupled. In previously presented implementation it is only possible to achieve proper steady state at the end of calculations, but not proper transient solution at each time-step. Would it be enough just to add additional loop (like, for example, in solidDisplacementFoam):

/************************************************** ***********************/
for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;


int iCorr = 0;
scalar initialResidual = 0;

do
{

.............

SOLVING MOMENTUM AND ENERGY EQUATIONS
.............


initialResidual = UEqn.solve().initialResidual();

} while (initialResidual > convergenceTolerance && ++iCorr < nCorr);


runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

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

Please correct me if I'm wrong at this point!
Actually it is also point of my interest: when does OpenFOAM switch values to new time step? Is it implemented after execution of "runTime++" operator? To achieve convergence criteria of both equations, old values in ddt operator should be GLOBAL OLD and do not change in internal "do {} while ()" loop.

Thank you! Waiting for your comments!

makaveli_lcf March 27, 2009 04:16

Ok, I found a solution, thank you all for attention.

mahaputra April 13, 2009 11:29

Dear Dr. A. Vakhrushev

which solver you used for simulation?

Regards

Nugroho

makaveli_lcf April 14, 2009 02:16

Quote:

Originally Posted by mahaputra (Post 212718)
Dear Dr. A. Vakhrushev

which solver you used for simulation?

Regards

Nugroho

Hi Nugroho!
It is self-made solver on the basis of icoFoam (for laminar flow ofcourse).


All times are GMT -4. The time now is 13:37.