CFD Online Discussion Forums - Solving coupled scalar transport equations in openFoam

Hi, I am trying to solve a series of coupled scalar transport equations i.e. temperature and chemical species transport for reacting flows by modifying icoFoam. I understand from a past forum post (quite a while ago) (https://www.cfd-online.com/Forums/op...equations.html), that an iterative scheme is needed within openFoam.

If anyone could share what is the best practice for doing so, or if there a better way of solving these coupled equations, that would be immensely helpful. The model equations are as follows:
$\frac{\partial T}{\partial t} + \bold{u}\cdot \nabla T = \alpha \nabla^{2}T + \sum_{i}(-\Delta H)r$

$\frac{\partial C_{A}}{\partial t} + \bold{u}\cdot \nabla C_{A} = D\nabla^{2}C_{A} + r$

where

is given by a standard arrhenius expression which couples the two transport equations.:
$r = Aexp(\frac{E_{a}}{RT})C_{A}$

$A, E_{a}$ are the preexponential factor and activation energy and

is the ideal gas constant.

I intend to solve this part outside the PISO loop as the temperature and chemical species do not strongly affect flow and can be neglected. But I am not too sure how to couple these two equations within openFoam. I understand that the following implementation provides a 1 way coupling of the velocity to the temperature equation:

```
// --- PISO loop
while (piso.correct())
{
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);

adjustPhi(phiHbyA, U, p);

// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU);

// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{
// Pressure corrector

fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);

pEqn.setReference(pRefCell, pRefValue);

pEqn.solve();

if (piso.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}
}

#include "continuityErrs.H"

U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
}
// This is where we add the temperature transport equations:
fvScalarMatrix TEqn
(
fvm:: ddt(T)
// The phi bit is a function argument which represents a vol<Type>Field
+ fvm:: div(phi, T)
- fvm:: laplacian(alpha, T)
// We need to add the source term below here:
);

TEqn.solve();
```
Thanks for your help!