[JasonG]

I have been trying to get back into using OpenFOAM after being away from it for several months, and I am having some difficulties while using simpleFoam. The case uses the simpleFoam solver set for laminar flow with 1 inlet and 4 outlet boundary conditions.

Attached is a clipped view showing the velocity and pressure magnitudes at the last iteration prior to crashing. The residual plot shows that my p residual was unable to converge at the specified tolerance, and then eventually crashed.

Does the error message below indicate anything obviously wrong with the setup of the model or a boundary condition? I am trying to determine if my mesh has issues, or if one of my boundary conditions is erroneous.


__________________________________________________ ____________________
Time = 232

smoothSolver: Solving for Ux, Initial residual = 0.448189, Final residual = 0.000520345, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.536472, Final residual = 0.00101717, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.638983, Final residual = 0.00085336, No Iterations 2
GAMG: Solving for p, Initial residual = 0.988656, Final residual = 9.96993e+66, No Iterations 100
GAMG: Solving for p, Initial residual = 0.913884, Final residual = 1.33824e+72, No Iterations 100
#0 Foam::error:rintStack(Foam::Ostream&) in "/opt/openfoam170/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::sigFpe::sigFpeHandler(int) in "/opt/openfoam170/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Uninterpreted:
#3 Foam::GAMGSolver::scalingFactor(Foam::Field<double >&, Foam::Field<double> const&, Foam::Field<double> const&, Foam::Field<double> const&) const in "/opt/openfoam170/lib/linuxGccDPOpt/libOpenFOAM.so"
#4 Foam::GAMGSolver::scalingFactor(Foam::Field<double >&, Foam::lduMatrix const&, Foam::Field<double>&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, Foam::Field<double> const&, unsigned char) const in "/opt/openfoam170/lib/linuxGccDPOpt/libOpenFOAM.so"
#5 Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMa trix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const in "/opt/openfoam170/lib/linuxGccDPOpt/libOpenFOAM.so"
#6 Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/opt/openfoam170/lib/linuxGccDPOpt/libOpenFOAM.so"
#7 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/opt/openfoam170/lib/linuxGccDPOpt/libfiniteVolume.so"
#8
in "/opt/openfoam170/applications/bin/linuxGccDPOpt/simpleFoam"
#9 __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#10
in "/opt/openfoam170/applications/bin/linuxGccDPOpt/simpleFoam"
__________________________________________________ _________________

[AlanR]

Jason,

I would look at your pressure boundary conditions - pressure was not converging before the crash. The standard setup for pressure is 'zeroGradient' at the inlet and 'fixedValue' 'uniform 0' at the outlet. For velocity, you should specify the inlet value and 'zeroGradient' or '$internalField' at the outlet. SimpleFoam with turbulence off should converge very quickly with a good mesh and the correct BCs and solvers. If pressure isn't converging on a similar slope to velocity, you don't need to wait for a crash, keep experimenting. In your fvSolution file, you might want to try PGC for pressure instead of GAMG. Basically, go back to the pitzDaily tutorial setup until you get your solution to converge. Then you can start experimenting with different solvers and BCs. You can run the checkMesh utility to see if there are mesh errors, but I think you have BC problems.

Good luck,

Alan

[JasonG]

Alan,

Thank you for the quick response. I believe I have the inlet and outlet BCs for P and U set as mentioned, but I will try "$internalField" instead of "zeroGradient" for the outlet velocity.

I have been tinkering with the settings in the fvSolution file, and set the following:

"relaxationFactors
{
p 0.1; // 0.3 is stable, decrease for bad mesh
U 0.1; // 0.7 is stable, decrease for bad mesh
}
"
Attached is a print screen of the current residuals, but it seems I have some more tinkering left to do.

I did run the checkMesh utility and the only error was 62 "highly skewed" faces.

Thanks again,

Jason

[AlanR]

Jason,

I would stay with the stable relaxation factors. The skewed faces in your mesh could be causing problems. You really want a mesh with no errors. If you used snappyHexMesh, you can adjust the mesh quality controls to generate a mesh with no skewed cells. Have you tried the PGC and PBiGC solvers?

Alan

[JasonG]

Alan,

I think I was originally using the PCG, but tried the GAMG based on someone's suggestion awhile ago (the reason escapes my memory now). I let the last solution run through the night, and the velocity and pressure components stayed around a constant residual of 4E-6 and 1E-3, respectfully.

Looking at the results, it appears I have 3-4 locations where the velocity and pressure gradients are clearly erroneous. Both components are orders of magnitude higher than the rest of the volume in a localized spot, and they are not in a location I would expect any sort of increased pressure/velocity. The rest of the volume does have results within the expected range.

On these next few trials, I will spend some more time on the mesh and try the PCG solver for the pressure. I have not been using snappyHexMesh, instead I have been using Ansys Workbench and exporting the fluent .msh data.

-Jason

[JasonG]

Update:

I have refined the mesh more and added a symmetry plane BC. I still have some skewed faces that I am having difficulty fixing due to the geometry. With the new mesh the solution converged with no anomalies in the pressure/velocity profiles, but the pressure residual is still high.

Now I am running with the fvSolution and fcSchemes copied from the pitzDaily simpleFoam tutorial, and the pressure residual seems to be staying much lower.

[AlanR]

The last residual plot looks a lot better. If that version keeps converging to where the initial residuals are in the 10 -6 range, I'd say you have a reasonable solution (10 -5 is probably OK).

Alan

[JasonG]

All residuals ended up around the 10 -6 range, so I am happy. Getting this model going and using OpenFOAM felt relatively "easy" (thanks in part to this forum). Hopefully in a month some prototype hardware will be ready to test, and I can start getting some test data to compare this to.

I wish I didn't have to wait until April for the next scheduled training .

[jorkolino]

Hello ,

I wish to intervene in the discussion with a problem which I encounter when simulating an extra scalar variable.

dimensionedScalar LMRAsource
    (
    "LMRAsource", dimensionSet(0, 0, -1, 0, 0, 0, 0), -1.0
    ); 

    fvScalarMatrix LMRAEqn
    (
        fvm::div(phi, LMRA)
      - fvm::laplacian(kappaEff, LMRA)
    ==
    LMRAsource
    );

    LMRAEqn.solve();
    LMRAEqn.relax();

I simulate actually one more variable in exactly the same way, and it has excellent results (confirmed by theory) but this one completely outrageous. The only difference between the two is the boundary conditions and the source term (1 vs. -1). For the problematic variable I use precisely the boundary conditions that "pressure" uses, i.e. zero value at outlets (two outlets I have) and zerogradient at inlet. Well, the itteration history of this variable has common with everything but logic. Mainly, the cells right adjacent to outlet get fantastic values in both directions. Why is that? Once again, I simulate two variables with only difference in B.C. After all, how can this zero value at outlet be so devastating for the simulation of one variable? Probably I need to check the coefficient matrix, to see wether it is composed correctly?

[JasonG]

I have been trying to use the simpleFoam solver on another test case that has a predicted Reynolds number around 3,000. Based on the cross-section at the portion of the model that has a high enough velocity to be in transitional or turbulent flow regime, I set approximate values for the Turbulent Kinetic Energy, Turbulent Length Scale, and Turbulent Dissipation coefficients.

The residuals started to drop, and then crashed just after 350 iterations. Is the simpleFoam turbulence model only suitable for higher Re flow, and does the entire model have to be at this sufficiently high Re number? Also, the mesh is comprised of fully tetrahedral elements with the prismatic boundary layer at the boundary surfaces; as I have not been able to get a good hexahedral mesh with my current mesher. Is a hexahedral mesh as detrimental to achieve as it is for a transient LES model?

Thanks.

[khedar]

Hey Jason,
Were you able to figure out the reason for this sudden crash. I too get this kind of behaviour sometimes when the residuals are increasing and all of a sudden it crashes.

[JasonG]

I cannot recall the exact issue I had with that particular model, but in general I would now check for the following:

• Ensure inlet/outlet constraints are sufficiently distanced from regions of adverse pressure/velocity grsdients. When the constraints are not in sufficiently developed flow, the solution can go unstable as the main domain begins resolving itself and is unable to satisfy the boundary condition imposed.
• Check for mesh orthogonality errors via check-mesh. Orthogonality errors can result in accumulation of violations in conservation of energy as the errors build and build at the poor mesh locations.
• View results near the timestep of the run-away instability. This may show areas of poor mesh quality that did not dominate the flow domain until the flow become more resolved in the solution.
• Experiment and tinker. Does the solution crash at certain flowrates, viscosity, RAS models but not others? Try to hone in on these transitions from stabilty, and think about how the mesh and governing equations may be sensitive to those variables.

I by no stretch claim to be a CFD expert, but these are just some suggestions I have inferred from my tinkering over the years. Good luck!

[khedar]

Orthogonality errors can result in accumulation of violations in conservation of energy as the errors build and build at the poor mesh locations.

My mesh some non-orthogonality but the value is around 44(max) and 13 (average). I am using ACMI interfaces though for coupling two regions. I will try to look at these regions for any unphysical results.

When I run the same case in CFX though, I don't get any problem, the residuals stagnate but no crashing occurs. Probably has some thing to with CFX's coupled solver.

I will do some more experiments on this issue.