# The velocity shows checkerboard feature

 Register Blogs Members List Search Today's Posts Mark Forums Read September 27, 2019, 22:52 The velocity shows checkerboard feature
#1
Member

Kai
Join Date: May 2010
Location: Shanghai
Posts: 61
Blog Entries: 1
Rep Power: 13 Hi Foamers,

I create a new steady state solver based on the standard buoyantSimpleFoam to study natural convection. Part of solid-fluid zones are modelled with porous media. Thus I add some features from the standard solver rhoPorousSimpleFoam. The major equation is coded as,
1. UEqn.H
Code:
```    // Construct the Momentum equation

MRF.correctBoundaryVelocity(U);

tmp<fvVectorMatrix> tUEqn
(
fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
==
fvOptions(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();

UEqn.relax();

// Include the porous media resistance and solve the momentum equation
// either implicit in the tensorial resistance or transport using by
// including the spherical part of the resistance in the momentum diagonal

tmp<volScalarField> trAU;
tmp<volTensorField> trTU;

if (pressureImplicitPorosity)
{
/****
tmp<volTensorField> tTU = tensor(I)*UEqn.A();
trTU = inv(tTU());
trTU.ref().rename("rAU");

fvOptions.constrain(UEqn);

for (int UCorr=0; UCorr<nUCorr; UCorr++)
{
U = trTU() & (UEqn.H() - gradp);
}
U.correctBoundaryConditions();

fvOptions.correct(U);
****/
}
else
{

fvOptions.constrain(UEqn);

//        fvOptions.correct(U);
if (simple.momentumPredictor())  // modify coeffs in rAU
{
solve
(
UEqn
==
fvc::reconstruct
(
(
// 2019-09-16
)*mesh.magSf()
)
);

fvOptions.correct(U);
}

trAU = 1.0/UEqn.A();
trAU.ref().rename("rAU");
}```
2. pEqn.H
Code:
```{
// 1. p --> p_rgh
// 2. comment implicit part

//    const volScalarField& psi = thermo.psi();

rho = thermo.rho();
rho.relax();

//    Info<< "Pre rho max/min : " << max(rho).value() << " " << min(rho).value()
//        << endl;

tmp<volVectorField> tHbyA;
if (pressureImplicitPorosity)
{
//        tHbyA = constrainHbyA(trTU()&UEqn.H(), U, p);
}
else
{
tHbyA = constrainHbyA(trAU()*UEqn.H(), U, p_rgh);
}
volVectorField& HbyA = tHbyA.ref();

tUEqn.clear();

surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*trAU()));
// 2019-09-16

bool closedVolume = false;

surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);

phiHbyA += phig;

// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);

while (simple.correctNonOrthogonal())
{
tmp<fvScalarMatrix> tp_rghEqn;

if (pressureImplicitPorosity)
{
//            tpEqn =
//            (
//                fvm::laplacian(rho*trTU(), p)
//              + fvOptions(psi, p, rho.name())
//             ==
//                fvc::div(phiHbyA)
//            );
}
else
{
tp_rghEqn =
(
fvm::laplacian(rhorAUf, p_rgh)
//              + fvOptions(psi, p, rho.name())
==
fvc::div(phiHbyA)
);
}

fvScalarMatrix& p_rghEqn = tp_rghEqn.ref();

p_rghEqn.setReference
(
pRefCell,
getRefCellValue(p_rgh, pRefCell)
);

p_rghEqn.solve();

if (simple.finalNonOrthogonalIter())
{
// Calculate the conservative fluxes
phi = phiHbyA - p_rghEqn.flux();

// Explicitly relax pressure for momentum corrector
p_rgh.relax();

// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U = HbyA + trAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}

#include "continuityErrs.H"

// 2019-09-16
//    p = p_rgh + rhok*gh;
p = p_rgh + rho*gh;
/****
if (pressureImplicitPorosity)
{
U = HbyA - (trTU() & fvc::grad(p));
}
else
{
}

U.correctBoundaryConditions();
fvOptions.correct(U);
****/
//    pressureControl.limit(p);

// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
/****
if (!thermo.incompressible() && closedVolume)
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
p_rgh = p - rho*gh;
}
****/
rho = thermo.rho();     // depends on p and T
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
<< endl;

// update display var

}```
3.EEqn.H
Code:
```{
volScalarField& he = thermo.he();

/****
dimensionedScalar Cp1
(
"liquidHeatCapacity",
dimEnergy/dimMass/dimTemperature,
thermo.Cp()
);

volScalarField Cp = thermo.Cp();
****/

gammaPhi = gammaf*phi;

volScalarField alphaMixEff
(
"alphaMixEff",
gamma*turbulence->alphaEff() + (1.0-gamma)*lambda2/Cp1
);
// h1 -> fluid
// h2 -> solid
// h2 = Cp2 *T = Cp2 * h1/Cp1
// laplacian(alpha1*gamma, h1) + laplacian(alpha2*(1-gamma), h2)
//=laplacian(alpha1*gamma, h1) + laplacian(alpha2*Cp2/Cp1*(1-gamma), h1)
// alpha = conductivity/heat capacity [kg/m/s]
// Cp  [J/kg/K]
// lambda: conductivity  [kg m /s^3 /K]
// alpha = nu / prandtl
// the following he refers to enthalpy of fluid
fvScalarMatrix EEqn
(
fvm::div(gammaPhi, he)
+ (
he.name() == "e"
? fvc::div(gammaPhi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
: fvc::div(gammaPhi, volScalarField("K", 0.5*magSqr(U)))
)
- fvm::laplacian(alphaMixEff, he)
==
rho*(U&g)*gamma
+ fvOptions(rho, he)
);

fvScalarMatrix EEqnFvO
(
- fvOptions(rho, he)
);

const scalarField& V = mesh.V();
const scalarField& eqnSource = EEqnFvO.source();

scalarField& heSource = volHeatSource;
heSource = eqnSource/V;

//    volHeatSource = EEqnFvO.source()/mesh.V();

//    forAll(mesh.C(), celli)
//    {
//        volHeatSource[celli] = eqnSource[celli]/V[celli];
//    }

EEqn.relax();

fvOptions.constrain(EEqn);

EEqn.solve();

fvOptions.correct(he);

thermo.correct();

Cp1 = thermo.Cp();

// 2019-09-16
//    rhok = rho*(1.0 - beta*(T - TRef));
//    rhok.relax();
}```
In brief, the porous model is implemented in an explicit way.

However I found the velocity along the channel axis shows checkerboard feature, which is located in porous zones. However, there is no checkerboard out of the porous zones. And the temperature is good everywhere.

In details, there are a lot of channel (in which the cross section is square) hang in the space. The fluid will flow into the channel from bottom and flow out the channel from top, as shown in the figure. There are spaces between two channels, and we call it as interchannel. The fluid is flow from top to bottom in interchannels, which is also shown in the figure. Besides, there are holes in the channel wall so that the fluid may flow through the holes aside. The channel wall is modelled as baffles. The major difference between channel and interchannel lies in: 1. the fuild in channels is heated by volumetric source and interchannel no vol heat sources; 2. there are solid structure in channels so the channel is modelled as porous media, whereas interchannel is pure fluid.
Uz: p_rgh: T: I also investigated the reason of checkerboard feature in the forum. The problem should be solved by Rhie-Chow interpolation,
make a "pseudo-staggered" version of the pressure equation in co-located grid, as discussed in,
Differences in solution method for pisoFoam and buoyantBoussinesqPisoFoam
Causes of checkerboarding

I swithed the divSchemes from upwind to bounded Gauss limitedLinear but it doesnt help.

I dont understand why velocity field has checkerboard feature only in porous zones. The source of porous resistance is add as
Code:
`pZones.addResistance(UEqn);`
in my case, I selected Darcy Forchheimer submodel, so the source is added in an explicit way into the matrix in DarcyForchheimerTemplates.C as,
Code:
```        forAll(cells, i)
{
const label celli = cells[i];
const label j = this->fieldIndex(i);
const tensor D = dZones[j];
const tensor F = fZones[j];

AU[celli] += mu[celli]*D + (rho[celli]*mag(U[celli]))*F;
}```
Is there any trouble so that the velocity shows checkerboard feature?
Attached Images Uz.jpg (42.2 KB, 85 views) p_rgh.jpg (52.4 KB, 85 views) T.jpg (51.6 KB, 87 views)
__________________
Kai  Tags baffles, checkerboarding, natural convection, porous domain Thread Tools Search this Thread Show Printable Version Email this Page Search this Thread: Advanced Search Display Modes Linear Mode Switch to Hybrid Mode Switch to Threaded Mode Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules Similar Threads Thread Thread Starter Forum Replies Last Post mykkujinu2201 Main CFD Forum 1 August 12, 2017 14:15 f.black Main CFD Forum 0 August 2, 2017 10:35 sdfij6354 OpenFOAM Running, Solving & CFD 3 July 26, 2017 17:16 Antech Main CFD Forum 0 April 25, 2006 03:15 chong chee nan FLUENT 0 December 29, 2001 06:13

All times are GMT -4. The time now is 09:52.