CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

The velocity shows checkerboard feature

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 27, 2019, 22:52
Default The velocity shows checkerboard feature
  #1
Member
 
Kai
Join Date: May 2010
Location: Shanghai
Posts: 61
Blog Entries: 1
Rep Power: 16
kaifu is on a distinguished road
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();
        pZones.addResistance(UEqn, tTU.ref());
        trTU = inv(tTU());
        trTU.ref().rename("rAU");

        fvOptions.constrain(UEqn);

        volVectorField gradp(fvc::grad(p));

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

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

        fvOptions.constrain(UEqn);

//        solve(UEqn == -fvc::grad(p));
//        fvOptions.correct(U);
        if (simple.momentumPredictor())  // modify coeffs in rAU
        {
            solve
            (
                UEqn
             ==
                fvc::reconstruct
                (
                    (
            // 2019-09-16
            //          - ghf*fvc::snGrad(rhok)
                      - ghf*fvc::snGrad(rho)
                      - fvc::snGrad(p_rgh)
                    )*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
//    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());

    bool closedVolume = false;

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

    closedVolume = adjustPhi(phiHbyA, U, p_rgh);

    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 = HbyA - trAU()*fvc::grad(p);
    }

    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
//    gradP_rgh = fvc::grad(p_rgh);
//    gradRho =  fvc::grad(rho);
//    driveSource = -gradP_rgh - gh_display*gradRho;

}
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
      + radiation->Sh(thermo, he)
      + 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();
    radiation->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
File Type: jpg Uz.jpg (42.2 KB, 155 views)
File Type: jpg p_rgh.jpg (52.4 KB, 158 views)
File Type: jpg T.jpg (51.6 KB, 159 views)
__________________
Kai
kaifu is offline   Reply With Quote

Reply

Tags
baffles, checkerboarding, natural convection, porous domain

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Problem with Velocity Poisson Equation and Vector Potential Poisson Equation mykkujinu2201 Main CFD Forum 1 August 12, 2017 14:15
Superimposing 2 Velocity fields For Model inputs in Comsol f.black COMSOL 0 August 2, 2017 10:35
atmBoundaryLayerInletVelocity - Velocity Profile not continuous through domain sdfij6354 OpenFOAM Running, Solving & CFD 3 July 26, 2017 17:16
Neumann pressure BC and velocity field Antech Main CFD Forum 0 April 25, 2006 03:15
what the result is negatif pressure at inlet chong chee nan FLUENT 0 December 29, 2001 06:13


All times are GMT -4. The time now is 12:30.