CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   checkerboard pressure and velocity field (http://www.cfd-online.com/Forums/openfoam-programming-development/115524-checkerboard-pressure-velocity-field.html)

pascool April 1, 2013 17:04

checkerboard pressure and velocity field
 
Hi foamers,

I'm computing the periodic flow through a ribbed duct. To achieve the desired flow rate, I coupled simpleFoam and channelFoam as follows:

Code:

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

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

        // --- Pressure-velocity SIMPLE corrector

        // U equation--------------------------------------------------------------------------------------------

      // Momentum predictor

      tmp<fvVectorMatrix> UEqn
      (
          fvm::div(phi, U)
        + turbulence->divDevReff(U)
        ==
          flowDirection*gradP
      );

      UEqn().relax();

      sources.constrain(UEqn());

      solve(UEqn() == -fvc::grad(p));

      // p equation ---------------------------------------------------------------------------------------------
   

          p.boundaryField().updateCoeffs();

          volScalarField rAU(1.0/UEqn().A());
          U = rAU*UEqn().H();
          UEqn.clear();

          phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf();
          adjustPhi(phi, U, p);

          // Non-orthogonal pressure corrector loop
          while (simple.correctNonOrthogonal())
          {
              fvScalarMatrix pEqn
              (
                  fvm::laplacian(rAU, p) == fvc::div(phi)
              );

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

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

          #include "continuityErrs.H"

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

        // Momentum corrector
          U -= rAU*fvc::grad(p);
          U.correctBoundaryConditions();
          sources.correct(U);

          turbulence->correct();

      // Correct driving force for a constant mass flow rate

      // Extract the velocity in the flow direction
      dimensionedScalar magUbarStar =
            (flowDirection & U)().weightedAverage(mesh.V());

      // Calculate the pressure gradient increment needed to
      // adjust the average flow-rate to the correct value
      dimensionedScalar gragPplus =
          (magUbar - magUbarStar)/rAU.weightedAverage(mesh.V());
      U += flowDirection*rAU*gragPplus;

      gradP += gragPplus;

      Info<< "Uncorrected Ubar = " << magUbarStar.value() << tab
          << "pressure gradient = " << gradP.value() << endl;

      runTime.write();

      #include "writeGradP.H"

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

    Info<< "End\n" << endl;

    return 0;
}

The code works, but I get a strange checkerboard pattern.

https://www.dropbox.com/s/43r2go9ohz...eckerBoard.png

Has anybody an idea where it may come from and how I could remedy this issue? I already tried different meshes.

Any help appreciated!!!
Best,
Pascal

pascool April 1, 2013 17:13

1 Attachment(s)
The mentioned checkerBoard pattern:https://www.dropbox.com/s/43r2go9ohz...eckerBoard.pnghttp://www.cfd-online.com/home/pascal/checkerBoard.pnghttp://www.cfd-online.com/home/pascal/checkerBoard.png


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