[kl1812]

Hi everyone,

I have been using OpenFOAM for a few months as part of my Masters project. Currently I am trying to simulate heave and surge of a truncated cylinder in still water, using interDyMFoam. I've been having some trouble with my simulations, they seem to crash as the time step becomes very small. I am using OpenFOAM 2.3.0, and I'm running the case in parallel on a cluster. I have created the mesh in gmsh and from the dropbox link below are screenshots of the domain from gmsh and ParaView. The total domain height is 1.3 m and the water depth is 1 m.

https://www.dropbox.com/sh/i4izf7xz1...5RrafMXNa?dl=0

I will now list the following files:
1. checkMesh log file
2. fvSolution
3. fvSchemes
4. controlDict
5. dynamicMeshDict
6. pointDisplacement
7. Lines where run goes wrong

1. checkMesh log

/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.3.0 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 2.3.0
Exec : checkMesh
Date : May 02 2016
Time : 12:14:24
Host : "login-2"
PID : 8762
Case : /export111/work/kl1812/CASEFOLDER/CYLDYN_D009_T1_A015_TR06_heave
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = -3

Time = -3

Mesh stats
points: 1372338
faces: 4016912
internal faces: 3918088
cells: 1322500
faces per cell: 6
boundary patches: 7
point zones: 0
face zones: 0
cell zones: 1

Overall number of cells of each type:
hexahedra: 1322500
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
polyhedra: 0

Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology
bottom 10522 10818 ok (non-closed singly connected)
symmetric 34820 35331 ok (non-closed singly connected)
inlet 5330 5502 ok (non-closed singly connected)
back 28340 28689 ok (non-closed singly connected)
outlet 5330 5502 ok (non-closed singly connected)
cylinder 4464 4593 ok (non-closed singly connected)
atmosphere 10018 10320 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (-3 0 -1) (3 1 0.3)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (-1.3031048e-15 -1.6108113e-14 -3.6396876e-15) OK.
Max cell openness = 2.9759405e-16 OK.
Max aspect ratio = 9.6552261 OK.
Minimum face area = 3.5702381e-06. Maximum face area = 0.00075188365. Face area magnitudes OK.
Min volume = 3.5702381e-08. Max volume = 7.5188365e-06. Total volume = 7.7971397. Cell volumes OK.
Mesh non-orthogonality Max: 42.803886 average: 6.7349899
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.57632634 OK.
Coupled point location match (average 0) OK.

Mesh OK.

End

2. fvSolution:

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.3.0 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
"alpha.water.*"
{
nAlphaCorr 2;
nAlphaSubCycles 1;
alphaOuterCorrectors yes;
cAlpha 1;

MULESCorr yes;
nLimiterIter 3;

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

pcorr
{
solver PCG;
preconditioner DIC;
tolerance 1e-5;
relTol 0;
}

pcorrFinal
{
$pcorr;
relTol 0;
}

p_rgh
{
solver PCG;
preconditioner DIC;
tolerance 1e-07;
relTol 0.05;
}

p_rghFinal
{
$p_rgh;
relTol 0;
}

"(U|k|epsilon).*"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-06;
relTol 0;
}

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

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

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

}

PIMPLE
{
momentumPredictor no;
nOuterCorrectors 1;
nCorrectors 3;
nNonOrthogonalCorrectors 4;
}

relaxationFactors
{

equations
{
U 0.9;
p 0.9;
p_rgh 0.9;
".*" 1;
}
}


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

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

ddtSchemes
{
default Euler;
}

gradSchemes
{
// default cellLimited fourth 1;
default Gauss linear;
grad(U) Gauss linear;
grad(alpha1) Gauss linear;
}

divSchemes
{
// default Gauss linear;
div(rhoPhi,U) Gauss limitedLinearV 1;
// div(phi,gamma) Gauss vanLeer;
// div(rho*phi,U) Gauss MUSCL;
div((muEff*dev(T(grad(U))))) Gauss linear;
div(phi,alpha) Gauss MUSCL;
div(phirb,alpha) Gauss interfaceCompression;
// div(rhoPhi,U) Gauss limitedLinearV 1;
}

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p_rgh;
pcorr;
alpha.water;
}

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

4. controlDict
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object controlDict;
}

startFrom startTime;

startTime -3;

stopAt endTime;

endTime 3;

deltaT 0.0001;

writeControl adjustableRunTime;

writeInterval 0.1;

purgeWrite 0;

writeFormat binary;

writePrecision 8;

writeCompression uncompressed;

timeFormat general;

timePrecision 8;

runTimeModifiable no;

adjustTimeStep yes;

maxCo 0.50;

maxAlphaCo 0.5;

maxDeltaT 0.01;

functions
{

forces
{
type forces;
functionObjectLibs ("libforces.so");
outputControl timeStep;
outputInterval 2;
patches (cylinder);
pName p;
UName U;
log true;
rhoName rho;
rhoInf 1000;
nuInf 1e-6;
CofR ( 0 0 -1.2 );
}

forcesDynamics
{
type forces;
functionObjectLibs ("libforces.so");
outputControl timeStep;
outputInterval 2;
patches (cylinder);
pName p_rgh;
UName U;
log true;
rhoName rho;
rhoInf 1000;
nuInf 1e-6;
CofR ( 0 0 -1.2 );
}

patches
{
type surfaces;
functionObjectLibs ("libsampling.so");
outputControl outputTime;
surfaceFormat vtk;
fields (p p_rgh);
interpolationScheme cellPointFace;
surfaces
(
cylinder
{
type patch;
patches (cylinder);
triangulate true;
interpolate true;
}
);

}

}


5. dynamicMeshDict
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object motionProperties;
}

dynamicFvMesh dynamicMotionSolverFvMesh;
motionSolverLibs ("libfvMotionSolvers.so");
solver displacementLaplacian;

displacementLaplacianCoeffs
{

diffusivity inverseDistance 1(cylinder);
}


6. pointDisplacement
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5-dev |
| \\ / A nd | Revision: 1740 |
| \\/ M anipulation | Web: http://www.OpenFOAM.org |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class pointVectorField;
object pointDisplacement;
}

dimensions [0 1 0 0 0 0 0];

internalField uniform (0 0 0);

boundaryField
{

inlet
{
type fixedValue;
value uniform (0 0 0);
}

outlet
{
type fixedValue;
value uniform (0 0 0);
}

symmetric
{
type symmetryPlane;
}

atmosphere
{
type fixedValue;
value uniform (0 0 0);
}

bottom
{
type fixedValue;
value uniform (0 0 0);
}

back
{
type fixedValue;
value uniform (0 0 0);
}

cylinder
{
type oscillatingDisplacement;
omega 6.283185;
amplitude (0 0 -0.15); \\ for heave motions
\\ amplitude (-0.15 0 0); \\ for surge motions
value uniform (0 0 0);

}
}


7. Error message: (for heave motions)

Interface Courant Number mean: 4.0286563e+21 max: 4.4729917e+25
Courant Number mean: 4.5707886e+25 max: 2.0822537e+33
deltaT = 5.0499871e-98
--> FOAM Warning :
From function Time:perator++()
in file db/Time/Time.C at line 1055
Increased the timePrecision from 11 to 12 to distinguish between timeNames at time -2.4862885
Time = -2.48628851541

GAMG: Solving for cellDisplacementx, Initial residual = 0, Final residual = 0, No Iterations 0
GAMG: Solving for cellDisplacementy, Initial residual = 0, Final residual = 0, No Iterations 0
GAMG: Solving for cellDisplacementz, Initial residual = 0.26870788, Final residual = 4.7402301e-06, No Iterations 6
Execution time for mesh.update() = 0.89 s
DICPCG: Solving for pcorr, Initial residual = 1, Final residual = 9.0295617e-06, No Iterations 370
DICPCG: Solving for pcorr, Initial residual = 0.77634439, Final residual = 9.6291055e-06, No Iterations 328
DICPCG: Solving for pcorr, Initial residual = 0.82315231, Final residual = 8.7830744e-06, No Iterations 339
DICPCG: Solving for pcorr, Initial residual = 0.86288272, Final residual = 9.4384348e-06, No Iterations 302
DICPCG: Solving for pcorr, Initial residual = 0.90947341, Final residual = 8.9003081e-06, No Iterations 322
time step continuity errors : sum local = 3.8561497e-09, global = 5.4902406e-12, cumulative = -1.4254648e+23
smoothSolver: Solving for alpha.water, Initial residual = 2.7574382e-06, Final residual = 3.8891055e-10, No Iterations 1
Phase-1 volume fraction = 0.76926425 Min(alpha1) = -6.1210092e-06 Max(alpha1) = 1.0021505
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.76926425 Min(alpha1) = -0.31261314 Max(alpha1) = 1.0021505
DICPCG: Solving for p_rgh, Initial residual = 1, Final residual = 0.049260868, No Iterations 192
DICPCG: Solving for p_rgh, Initial residual = 0.12477516, Final residual = 0.0042661957, No Iterations 308
DICPCG: Solving for p_rgh, Initial residual = 0.83326361, Final residual = 0.041521828, No Iterations 328
DICPCG: Solving for p_rgh, Initial residual = 0.91849496, Final residual = 0.045002832, No Iterations 347
DICPCG: Solving for p_rgh, Initial residual = 0.91435506, Final residual = 0.035006983, No Iterations 338
time step continuity errors : sum local = 2.3111299, global = 0.47075552, cumulative = -1.4254648e+23
[1] #0 Foam::error:rintStack(Foam::Ostream&) in "/apps/OpenFOAM/2.3.0//OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[1] #1 Foam::sigFpe::sigHandler(int) in "/apps/OpenFOAM/2.3.0//OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
[1] #2 in "/lib64/libc.so.6"
[1] #3 double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) in "/apps/OpenFOAM/2.3.0//OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"

...etc


For heave motions, the model runs for 0.5 s before crashing, due to the very small delta T. Running checkMesh through these timesteps gives me a non-orthogonality max of around 80 which is quite high, but the output is still mesh OK. I have tried changing small things in the fvSolution file and controlDict but results are still the same. I also added the under-relaxation factors in the fvSolution file to see if that helped, but it hasn't.

For surge motions, the model will not even run for 0.1 s, and crashes very quickly due to the small timestep again. I was able to run one surge case for a small amplitude of 0.02 m, but on running checkMesh the mesh failed some of the checks at certain timesteps. I have then also tried using solidBodyMotionFvMesh instead without success, I'm not quite sure whether it is best to use this one or dynamicMotionSolverFvMesh.

Any help is greatly appreciated. I don't have a lot of experience with OpenFOAM so if anyone has experienced similar problems or is running similar cases please get in touch. Do let me know if you need more documentation on the cases I have been running.

Kind regards,
Karoline

[mo_na]

Hi Karoline,

I experienced something similar and in my case the simulation actually crashed due to non-orthogonal cells. Maybe you can try to adapt your controlDict to write the last time close to when your simulation crashed and see if the non-orthogonality gets higher, if yes you can use the following command after running checkMesh:

foamToVTK -faceSet nonOrthogonalFaces -time 0.5

Of course you would need to adjust the time and check if thats the right name for the faceset, I am just guessing.
Another idea is to run the Co-utility which calculates the courant number and works just by typing Co in your case folder. This way you can check where you have the biggest problems with the courant number.

Hope that might help you somehow.
Cheers,
Mona

[kl1812]

Hi Mona,

Thanks for your help. You were right, my simulations were crashing due to non-orthogonal cells. I had to use a different mesh from the one I posted above in order to get the simulations to work. I also found using the commands moveMesh and checkMesh before running the simulations very useful, as I could check that the mesh deformation was okay first without having to run the full simulation!

Many thanks,
Karoline