[MikeWorth]

I'm running an axisymmetric interfoam model, with a moving, submerged boundary (via AMIs). A couple of seconds into the simulation I start seeing high velocities near the moving surface, which then lead to low values of alpha somehow. This looks to be a numerical error of some kind as the only boundaries nearby are wedge and wall type, and there is no air (i.e. low alpha) anywhere near this point.

I'm not exactly sure how the symmetry implicitly works with the wedges, but it seems fishy that it's happening right here - could something be slightly off with exact mesh positioning causing the maths to go wonky?

Here is what it looks like as the interface begins to form - the white outline is an alpha=0.5 contour. Velocities are an order of magnitude faster than the surface is moving:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1912                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;

    div(rhoPhi,U)   Gauss linearUpwind grad(U);

    div(phi,alpha)  Gauss vanLeer;
    div(phirb,alpha) Gauss interfaceCompression;

    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

oversetInterpolation
{
    method         cellVolumeWeight;
}

fluxRequired
{
    default         no;
    pcorr           ;
    p               ;
}

oversetInterpolationSuppressed
{
    grad(p_rgh);
    surfaceIntegrate(phiHbyA);
}

// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1912                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    cellDisplacement
    {
        solver          GAMG;
        tolerance       1e-5;
        relTol          0;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    "alpha.latex.*"
    {
        nAlphaCorr      3;
        nAlphaSubCycles 2;
        cAlpha          1;
        icAlpha         0;
        scAlpha         0;

        MULESCorr       yes;
        nLimiterIter    5;
        alphaApplyPrevCorr  no;

        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0;
    }


    p_rgh
    {
        solver          GAMG;
        smoother        DIC;

        tolerance       1e-8;
        relTol          0.01;

    }

    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    pcorr
    {
        $p;
        solver          PCG;
        preconditioner  DIC;
    }

    pcorrFinal
    {
        $pcorr;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0.01;
        maxIter         200;
        minIter         1;
    }

    UFinal
    {
        $U;
        relTol          0;
    }

    cellMotionUy
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-08;
        relTol          0;
    }

}

PIMPLE
{
    momentumPredictor   no;
    correctPhi          no;
    nOuterCorrectors    2;
    nCorrectors         3;
    nNonOrthogonalCorrectors 0;

    ddtCorr             true;

    //Not needed if we have a fixed pressure BC
    //pRefPoint           (1e-3 0.001 0.00);
    //pRefValue           0.0;
}

relaxationFactors
{
   fields
   {
   }
   equations
   {
       ".*"             1;
   }
}

// ************************************************************************* //

Thanks,
Mike

[kk415]

Hi Mike,

What is the CFL number you are using? Maybe you can try reducing your dt.

[Ship Designer]

Hello Mike,

What gets my attention in the code you posted is fluxRequired { p; } and pcorr { $p; }. Are you sure these shouldn't be set to p_rgh instead? Since there's no definition of p in your fvSolution file, I'm not sure what settings apply to pcorr. Additionally you may want to try setting the relTol for pcorr to e.g. 0.1 and see if that makes a difference.

Good luck, Claudio

[MikeWorth]

Thanks for the suggestions.

MaxCo is set to 0.25, although when things are running well the timestep is limited via the Brackbill capillary time constraint such that Co stays <0.1. When non-physical high velocities start appearing it does get pushed up to the 0.25 limit. I can try nudging that down to see if it makes a difference.

The fluxrequired and pcorr lines are unchanged from the tutorial (overInterDyMFoam/twoSimpleRotors) I originally copied many moons ago... I've changed a hell of a lot since then, so I'll try tweaking them as you suggest and see if that does the trick

Thanks,
Mike

[Ship Designer]

Hi Mike,

I've found the overInterDyMFoam tutorial you've mentioned and I've seen the p file in the 0 folder there. As you mentioned, you are using interFoam for your case however. I've noticed that as opposed to the p field interFoam normally puts out, the one in the tutorial has dimensions [0 2 -2 0 0 0 0] instead of [1 -1 -2 0 0 0 0]. I tried to find in the source code of overInterDyMFoam whether p is actually a required field but couldn't readily find any indication for it. I tried giving my plain (non-mesh overset) interFoam a faulty p file in the 0 folder to provoke an error but it seems to be ignored, even after including fluxRequired { p; } in fvSchemes. I tried with both Foundation and ESI versions of OpenFOAM. I also couldn't find any interFoam cases in the tutorials that use a p input field. Therefore I presume that interFoam doesn't read p fields as input at all. Do you really need p as input in your case?

In your case pcorr uses the settings of p, but the big question is, what are they and thus what are the actual settings for pcorr in your case?

You mention high velocities in the low alpha cells. What does p_rgh look like there? Are these values present only in one wedge boundary or in the coupled one as well, which they should, don't they?

One last question, what application did you use to make the mesh? Your viscous layers look very neat! Thanks.

Cheers, Claudio

[MikeWorth]

The p does seem to be a vestigial appendage leftover from the old tutorial. I don't define boundary or initial conditions for p in my 0/, although I do see p in the generated output (I can't remember if this is normal or not for interFoam).

Since starting this work I've switched from using overset to using AMI based movement - any leftover overset stuff is (hopefully) not actually used. I'm using v1912.

The fact that this isn't set up is puzzling, as I've been successfully running variants of the model for quite a while now, until finding a couple of variants that are unstable. By using "foamDictionary system/fvSolution" I can see how it gets interpreted, and the answer is that it is empty! The pcorr block looks like this:

    pcorr
    {
        solver          PCG;
        preconditioner  DIC;
    }

Swapping $p out for $p_rgh gives this instead:

    pcorr
    {
        solver          PCG;
        smoother        DIC;
        tolerance       1e-08;
        relTol          0.01;
        preconditioner  DIC;
    }

I'm not an expert on the subtleties of these configurations, but something strikes me as much better than nothing!

The mesh was made using snappyHexMesh, but in order to get what I wanted I ran it multiple times with only the layers part of the dictionary enabled with a series of different settings. This allowed me to insert a thickness distribution with gradation between the square background and an intermediate scale, a bunch of constant thickness cells at this intermediate scale, followed by a further gradation down to a fine boundary layer. This is based on an expected property of the interface which I need to capture with intermediate scale cells.

Other things I've been playing with are tightening up the mesh tolerances, by moving any points that have been warped off the 2D plane or away from the axis by SHM back into place. I've also increased the writePrecision during meshing to ensure that faces remain properly planar (previously I got warnings that they were off by 1e-15).

It's now very unstable - the issue seems to be at the actual interface. It could be a coincidence, but looks to be related to a cell splitting by SHM - "checkMesh -allGeometry" often moans about concave cells in these sorts of places, could this be the issue?



Thanks,
Mike

[MikeWorth]

I've been looking closer at the differences between versions that work and those that don't, and the SHM multi-layering approach is looking like a smoking gun. Somehow if I strip things back to only 2 sets of layers it still falls over, while working OK (at least initially) with it all done in one go.

The 'normal' approach like this works:

addLayersControls
{
    relativeSizes false;

    layers
    {
        "movingpart.*"
        {
            nSurfaceLayers 14;//Calculated using https://openfoamwiki.net/index.php/Scripts/blockMesh_grading_calculation
        }
    }

    // Expansion factor for layer mesh
    expansionRatio 1.258498951;

    firstLayerThickness 5e-6;

But if I try to make almost identical layers in two steps like this it doesn't:

Outer layers (first)

addLayersControls
{
    relativeSizes false;

    layers 
    {
        "movingpart.*"
        {
            //Calculated using https://openfoamwiki.net/index.php/Scripts/blockMesh_grading_calculation
            nSurfaceLayers 10;
        }
    }

    // Expansion factor for layer mesh
    expansionRatio 1.258498951;

    firstLayerThickness 12.50e-6;

Inner (second)

addLayersControls
{
    relativeSizes false;

    layers
    {
        "movingpart.*"
        {
            //Calculated using https://openfoamwiki.net/index.php/Scripts/blockMesh_grading_calculation
            nSurfaceLayers 5;
        }
    }

    // Expansion factor for layer mesh
    expansionRatio 1.258498951;

    firstLayerThickness 5e-6;

It seems odd that such a subtle change in the mesh completely ruins everything. Is there something subtle going on in the way that SHM pushes the old layers back that might be making a substantial difference? Nothing jumps out at me from checkMesh (although I might not be looking in the right place)

Thanks,
Mike

[Ship Designer]

Hello Mike,

I don't define boundary or initial conditions for p in my 0/, although I do see p in the generated output (I can't remember if this is normal or not for interFoam).

Yes it is! interFoam writes a p field as output which is defined as:

p == p_rgh + rho*gh;

The fact that this isn't set up is puzzling, as I've been successfully running variants of the model for quite a while now, until finding a couple of variants that are unstable.

This is possible, if for the other cases, whatever the pcorr settings, they were acceptable. Maybe the cases where it doesn't work are inherently less stable and only then it becomes an issue. It's probable that there are some standard settings for pcorr, but I don't know what they are. At least by seeing them now, you could try changing them and see the effect.


If the mesh non-orthogonality is okay, I don't think that concavity is a problem either, unless checkMesh explicitly reports a mesh error for that. I rather believe that the pcorr settings are very tight and I would suggest loosening them up. Try reducing tolerance to 1e-6 or 1e-5 and relTol to 0.1. See if it makes any difference at all or if the velocity buildup occurs later. I've had some stability problems in the past, mostly p_rgh buildups near the phase interface in places of error-free mesh and far away from the geometry and only this helped where nothing else would.

Normally the pcorr subdict contains a subdict for the preconditioner with its own tolerance settings, something like this:

"pcorr.*"  {
          solver  PCG;
  preconditioner  {
    preconditioner  GAMG;
          smoother  GaussSeidel;
         tolerance  1e-6;
            relTol  0;
  };
       tolerance  1e-6;
          relTol  0.1;
         minIter  1;
 };

I suggest to include that in your fvSolution and set the tolerance the same as for pcorr.

For the points off an AMI or symmetry plane, unless checkMesh or the interFoam complain, I would ignore that and try the above first, since high writePrecision settings bloat your case size. If you still want to correct this, since you use v1912 you could try use cfMesh's utility improveSymmetryPlanes that corrects exactly points being off-plane. Copy your polyMesh to a separate cfMesh case and let the utility work on it. Make sure the mesh is in ascii format, otherwise cfMesh will abort.

I find your viscous layer creation strategy in two passes very original! I take it sHM is unable to create the viscous layers you need in one step only? Can you please confirm whether the layers generated with the one-pass and two-passes method are ultimately the same? If not, the problem might simply occur due to the different geometry i.e. having more layers. If after the creation of the layers the mesh has no errors, then it's not a mesh quality or sHM issue I would say, but rather the flow that slightly changes due to different cells, giving magically rise to the high velocity. Is p_rgh high as well where this occurs?

Good luck and have a nice weekend, Claudio

[MikeWorth]

Thanks for the suggestions. I've tried adding that snipped to my fvSolution and adjusted the tolerance and relTol of the pcorr block. The regex also matches, but I'm guessing that the explicit match takes priority for that? Doesn't seem to help much. The problem area now seems to be firmly at the liquid interface, with high velocity but no visible spikes in p or p_rgh. Is this therefore indicative of spurious currents from the interfacial surface tension applying wrong?

The SHM approaches should give notionally the same mesh. having loaded them both into paraview over the top of each other in different colours, this seems broadly true. There are slight differences where the warping is subtly different, but they're basically the same. I've also diffed their checkMesh output and it's very similar. Bizarrely the problem seems to only occur on the double-layered version, so it seems that some very subtle change here is breaking things. Do any of the differences in checkmesh below look fishy to you?

--- ../checkMesh.SHMx1	2021-04-19 09:15:56.378157603 +0000
+++ ../checkMesh.SHMx2	2021-04-19 09:53:57.548388806 +0000
@@ -1,142 +1,142 @@
 /*---------------------------------------------------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  v1912                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 Build  : _f3950763fe-20191219 OPENFOAM=1912
 Arch   : "LSB;label=32;scalar=64"
 Exec   : checkMesh -allGeometry
 Date   : Apr 19 2021
-Time   : 09:15:54
+Time   : 09:53:56
 Host   : ip-172-31-8-169
-PID    : 5238
+PID    : 20125
 I/O    : uncollated
 Case   : /home/ubuntu/casename
 nProcs : 1
 trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
 fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
 allowSystemOperations : Allowing user-supplied system call operations
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 Create time
 
 Using #calc at line 69 in file "/home/ubuntu/casename/system/controlDict"
 Using #codeStream with "/home/ubuntu/casename/dynamicCode/platforms/linux64GccDPInt32Opt/lib/libcodeStream_90540247e08e92844c35514f8198f20ebc47281f.so"
 Using #calc at line 70 in file "/home/ubuntu/casename/system/controlDict"
 Using #codeStream with "/home/ubuntu/casename/dynamicCode/platforms/linux64GccDPInt32Opt/lib/libcodeStream_fae589c14ad73a95fb40be0cd9ce4cbb1179cdeb.so"
 Using #calc at line 72 in file "/home/ubuntu/casename/system/controlDict"
 Using #calc at line 73 in file "/home/ubuntu/casename/system/controlDict"
 Using #codeStream with "/home/ubuntu/casename/dynamicCode/platforms/linux64GccDPInt32Opt/lib/libcodeStream_c5283eb0bd727a9daf6cb00bc145aac9a2439283.so"
 Create mesh for time = 0
 
 Enabling all geometry checks.
 
 Time = 0
 
 Mesh stats
     points:           187080
     internal points:  0
     faces:            365251
     internal faces:   179277
     cells:            90642
     faces per cell:   6.00745791134
     boundary patches: 10
     point zones:      0
     face zones:       4
     cell zones:       4
 
 Overall number of cells of each type:
     hexahedra:     88815
     prisms:        626
     wedges:        0
     pyramids:      0
     tet wedges:    0
     tetrahedra:    0
     polyhedra:     1201
     Breakdown of polyhedra by number of faces:
         faces   number of cells
             6   6
             7   1088
             8   107
 
 Checking topology...
     Boundary definition OK.
     Cell to face addressing OK.
     Point usage OK.
     Upper triangular ordering OK.
     Face vertices OK.
    *Number of regions: 2
     The mesh has multiple regions which are not connected by any face.
   <<Writing region information to "0/cellToRegion"
   <<Writing region 0 with 7200 cells to cellSet region0
   <<Writing region 1 with 83442 cells to cellSet region1
 
 Checking patch topology for multiply connected surfaces...
                    Patch    Faces   Points                  Surface topology Bounding box
                 outerAMI      600     1202  ok (non-closed singly connected) (0.0249762055395 0 -0.00109048468413) (0.0249762055395 0.6 0.00109048468413)
                      top       37       77  ok (non-closed singly connected) (1.0597848687e-75 0.6 -0.0016357270262) (0.0374643083093 0.6 0.0016357270262)
                   bottom       37       77  ok (non-closed singly connected) (0 0 -0.0016357270262) (0.0374643083093 0 0.0016357270262)
                     back    90642    93755  ok (non-closed singly connected) (-4.83868159767e-40 0 -0.0016357270262) (0.0374643083093 0.6 1.00830515897e-43)
                    front    90642    93755  ok (non-closed singly connected) (-4.83868159767e-40 0 -2.11261400988e-41) (0.0374643083093 0.6 0.0016357270262)
                 innerAMI     1348     2698  ok (non-closed singly connected) (0.0249762055395 0 -0.00109048468413) (0.0249762055395 0.6 0.00109048468413)
                  thing      441      883  ok (non-closed singly connected) (0 0.3 -0.000785148938775) (0.0179828672143 0.340000003576 0.000785148938775)
           thing_wetted     1627     3255  ok (non-closed singly connected) (0 0.10000000149 -0.000959626516838) (0.0219790607557 0.3 0.000959626516838)
                 bathWall      325      652  ok (non-closed singly connected) (0.0374643083093 0 -0.0016357270262) (0.0374643083093 0.325 0.0016357270262)
           atmosphereSide      275      552  ok (non-closed singly connected) (0.0374643083093 0.325 -0.0016357270262) (0.0374643083093 0.6 0.0016357270262)
 
 Checking faceZone topology for multiply connected surfaces...
                 FaceZone    Faces   Points                  Surface topology Bounding box
                      top       25       51  ok (non-closed singly connected) (1.0597848687e-75 0.6 -0.00109048468413) (0.0249762055395 0.6 0.00109048468413)
-         centralTopFaces      102      205  ok (non-closed singly connected) (0 0.342870999563 -0.00109048468413) (0.0249762055395 0.343124069787 0.00109048468413)
-      centralBottomFaces      102      205  ok (non-closed singly connected) (0 0.0968728834598 -0.00109048468413) (0.0249762055395 0.0971234991189 0.00109048468413)
+         centralTopFaces      102      205  ok (non-closed singly connected) (0 0.342870810231 -0.00109048468413) (0.0249762055395 0.343123893846 0.00109048468413)
+      centralBottomFaces      102      205  ok (non-closed singly connected) (0 0.0968732213925 -0.00109048468413) (0.0249762055395 0.0971243427523 0.00109048468413)
                   bottom       25       51  ok (non-closed singly connected) (0 0 -0.00109048468413) (0.0249762055395 0 0.00109048468413)
 
 Checking basic cellZone addressing...
                 CellZone        Cells       Points       VolumeBoundingBox
              background         7200        15626 2.04271272065e-05 (0.0249762055395 0 -0.0016357270262) (0.0374643083093 0.6 0.0016357270262)
-          topBlockCells         7033        14532 6.99977663975e-06 (-4.83868159767e-40 0.342870999563 -0.00109048468413) (0.0249762055395 0.6 0.00109048468413)
-      centralBlockCells        73514       151236 3.07632624824e-06 (0 0.0968728834598 -0.00109048468413) (0.0249762055395 0.343124069787 0.00109048468413)
-       bottomBlockCells         2895         6096 2.64138591764e-06 (0 0 -0.00109048468413) (0.0249762055395 0.0971234991189 0.00109048468413)
+          topBlockCells         7032        14530 6.99971599829e-06 (-4.83868159767e-40 0.342870810231 -0.00109048468413) (0.0249762055395 0.6 0.00109048468413)
+      centralBlockCells        73515       151238 3.07638209319e-06 (0 0.0968732213925 -0.00109048468413) (0.0249762055395 0.343123893846 0.00109048468413)
+       bottomBlockCells         2895         6096 2.64139071415e-06 (0 0 -0.00109048468413) (0.0249762055395 0.0971243427523 0.00109048468413)
 
 Checking geometry...
     Overall domain bounding box (-4.83868159767e-40 0 -0.0016357270262) (0.0374643083093 0.6 0.0016357270262)
     Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
     Mesh has 3 solution (non-empty) directions (1 1 1)
-    Wedge back with angle 2.50000000138 degrees
-    Wedge front with angle 2.50000000138 degrees
+    Wedge back with angle 2.50000000133 degrees
+    Wedge front with angle 2.50000000133 degrees
     All edges aligned with or perpendicular to non-empty directions.
-    Boundary openness (1.64702991369e-14 -8.71103408019e-19 2.62100096573e-14) OK.
+    Boundary openness (1.85672186208e-14 -8.50625678685e-19 6.18299811514e-15) OK.
     Max cell openness = 1.53490106268e-15 OK.
     Max aspect ratio = 25.0289312037 OK.
     Minimum face area = 5.44039160854e-11. Maximum face area = 3.35723618435e-06.  Face area magnitudes OK.
     Min volume = 3.38966568599e-15. Max volume = 3.35723618435e-09.  Total volume = 3.31446160121e-05.  Cell volumes OK.
-    Mesh non-orthogonality Max: 25.3896599773 average: 2.69630443086
+    Mesh non-orthogonality Max: 25.386544632 average: 2.69326805544
     Non-orthogonality check OK.
     Face pyramids OK.
-    Max skewness = 0.465880892795 OK.
+    Max skewness = 0.465999415964 OK.
     Coupled point location match (average 0) OK.
     Face tets OK.
     Min/max edge length = 4.99999951583e-06 0.0032714540524 OK.
   <<Writing 488 near (closer than 6.01177408764e-07 apart) points to set nearPoints
     All angles in faces OK.
     Face flatness (1 = flat, 0 = butterfly) : min = 1  average = 1
     All face flatness OK.
-    Cell determinant (wellposedness) : minimum: 0.00809562991112 average: 0.396592299529
+    Cell determinant (wellposedness) : minimum: 0.00809562991112 average: 0.396372492196
     Cell determinant check OK.
- ***Concave cells (using face planes) found, number of cells: 360
-  <<Writing 360 concave cells to set concaveCells
-    Face interpolation weight : minimum: 0.294523453407 average: 0.483772637659
+ ***Concave cells (using face planes) found, number of cells: 358
+  <<Writing 358 concave cells to set concaveCells
+    Face interpolation weight : minimum: 0.294516579567 average: 0.483787294555
     Face interpolation weight check OK.
-    Face volume ratio : minimum: 0.12468161947 average: 0.932602925846
+    Face volume ratio : minimum: 0.124681118055 average: 0.932634914103
     Face volume ratio check OK.
 Calculating AMI weights between owner patch: outerAMI and neighbour patch: innerAMI
 AMI: Creating addressing and weights between 600 source faces and 1348 target faces
 AMI: Patch source sum(weights) min:0.999999999996 max:1 average:1
 AMI: Patch target sum(weights) min:0.999996339997 max:1.0000039496 average:1.00000000146
 
 Failed 1 mesh checks.
 
 End