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)

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