Hi everyone,
So I am using v10 to simulate a wind turbine case. I have chosen to model half the domain just to save running time (I really need to do this if possible due to available resources) Initially, I used MRF which ended up working for me, but I now need to use rigid body motion (so that I can use libAcoustics).
I am using the SA turbulence model
My issue is that I am having serious difficulty in just getting a simulation going with the rigid body motion and periodic faces. It will run for like 3 time steps, but then crash on the 3rd time step. The following is an excerpt from running pimpleFOAM (third time step):
Code:
Courant Number mean: 0.011671816 max: 1.5599877
Time = 0.0015s
PIMPLE: Iteration 1
AMI: Creating addressing and weights between 102419 source faces and 101922 target faces
AMI: Patch source sum(weights) min/max/average = 0.38156299, 1.0000005, 0.99971194
AMI: Patch target sum(weights) min/max/average = 0.69604475, 1.000114, 0.99999508
AMI: Creating addressing and weights between 47454 source faces and 47454 target faces
AMI: Patch source sum(weights) min/max/average = 0, 1.0137377, 0.99686404
AMI: Patch target sum(weights) min/max/average = 0, 1.0142746, 0.99638231
................
...................
...............
#0 Foam::error::printStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::divide(Foam::Field<double>&, double const&, Foam::UList<double> const&) at ??:?
#4 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::dimensioned<double> const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) in "/opt/openfoam10/platforms/linux64GccDPInt32Opt/bin/pimpleFoam"
#5 ? in "/opt/openfoam10/platforms/linux64GccDPInt32Opt/bin/pimpleFoam"
#6 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#7 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#8 ? in "/opt/openfoam10/platforms/linux64GccDPInt32Opt/bin/pimpleFoam"
Floating point exception (core dumped)
The AMI weight for the second set of faces has gradually decreased from the first time step until the third, when it reaches 0.
Note that the second set of AMI faces, where the 0 weight issue arises is actually the cyclicAMI faces that separate the rotating and non rotating region. These 2 faces are the exact same. Later I have described the sequence that I used to create these.
As you will see below, when I used moveMesh, the periodic faces also moved with the mesh. Does this indicate that they are not set correctly? I know that with MRF there is a "nonRotatingPatches" option, but it doesn't seem to be a thing in the dynamicMeshDict
What I have tried already:
- I have simplified my case down to just a cylinder in place of a blade to rule out any meshing issues with the blade
- I have generated a 360 degree model with rigid body motion (which appears to be moving fine after running for a while)
- Then I went back to the half domain, with periodic faces, and I have used the moveMesh command, which runs ok (but the entire region including the periodic faces are moving with the rest of the cell zone. (See in the image below)
- I have tried to use cyclicPeriodicAMI, but It does not appear to be implemented in OF10, then I tried cyclicSlip, but was getting incomprehensible error messages
So I am sort of stuck as to where to go from here, and I would appreciate any experience anyone has with this sort of thing. Furthermore, I haven't actually been able to find a single example of a tutorial or case on github that uses both cyclic rotational conditions and dynamic Meshing.
Below I have included
- fvSolution, fvSchemes, controlDict boundary File, boundary conditions, image of a few time steps into moveMesh command (with the top face of the domain removed), SHM dict and checkMesh output, dynamicMeshDict
The sequence I use to create the mesh is:
blockMesh
surfaceFeatures
decomposePar
mpirun -np 4 snappyHexMesh -parallel -overwrite > log.shm
reconstructParMesh -constant
checkMesh -writeSets vtk > log.checkmesh
createBaffles -overwrite
splitBaffles -overwrite
Thanks again, and if anyone wants more information and thinks they could help me, I can upload my entire case to github if needed. Thanks :)
fvSolution
Code:
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-07;
relTol 0;
maxIter 50;
}
pFinal
{
solver PCG;
preconditioner DIC;
tolerance 1e-07;
relTol 0;
maxIter 50;
}
"pcorr.*"
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0;
minIter 3;
}
U
{
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-07;
relTol 0;
maxIter 10;
}
UFinal
{
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-07;
relTol 0;
maxIter 0;
}
k
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
}
epsilon
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
}
nuTilda
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-20;
relTol 0;
maxIter 5;
}
nuTildaFinal
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-20;
relTol 0;
maxIter 5;
}
yPsi
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0;
minIter 2;
}
}
PIMPLE
{
//If enable you will need to add UFinal
momentumPredictor yes;
consistent yes;
//If you set nOuterCorrectors to 1 is equivalent to piso
nOuterCorrectors 3; // for best results, set outer correctors to at least 5 for moving bodies (wolf Dynamics)
nCorrectors 3; // 4 or more for highly transient flows or strongly coupled problems (WD)
nNonOrthogonalCorrectors 2; // at least 1, increase value for bad quality meshes
//pRefCell 0;
//pRefValue 0;
correctPhi yes;
checkMeshCourantNo yes;
}
relaxationFactors
{
fields
{
".*" 0.9;
//p 0.3;
}
equations
{
".*" 0.9;
//U 0.7;
//k 0.7;
//omega 0.7;
}
}
// ************************************************************************* //https://github.com/niall-oneill/open...ualization.png
fvSchemes
Code:
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
div(phi,k) Gauss linear;
div(phi,epsilon) Gauss linear;
div(phi,nuTilda) Gauss linear;
div((nuEff*dev(grad(U).T()))) Gauss linear;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default Gauss linear limited 0.5;;
/*
laplacian(nu,U) Gauss linear corrected;
laplacian(rAU,pcorr) Gauss linear corrected;
laplacian(rAUf,pcorr) Gauss linear corrected;
laplacian(rAU,p) Gauss linear corrected;
laplacian(nuEff,U) Gauss linear limited 0.5;
laplacian((1|A(U)),p) Gauss linear limited 0.5;
laplacian(DkEff,k) Gauss linear limited 0.5;
laplacian(DepsilonEff,epsilon) Gauss linear limited 0.5;
laplacian(DREff,R) Gauss linear limited 0.5;
laplacian(DnuTildaEff,nuTilda) Gauss linear limited 0.5;
laplacian(yPsi) Gauss linear limited 0.5;
laplacian((1|((1|(1|A(U)))-H(1))),p) Gauss linear limited 0.5;
laplacian((1|max(((1|(1|A(U)))-H(1)),(0.1|(1|A(U))))),p) Gauss linear limited 0.5;
*/
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default limited 0.5;
}
fluxRequired
{
default no;
pcorr;
p;
}
wallDist
{
method Poisson;
}
// ************************************************************************* //
controlDict
Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application pimpleFoam;
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 1003;
deltaT 0.0005;
writeControl runTime;
writeInterval 0.01;
purgeWrite 0;
writeFormat binary;
writePrecision 8;
writeCompression uncompressed;
timeFormat general;
timePrecision 6;
//runTimeModifiable yes;
//adjustTimeStep yes;
//maxCo 5; //unstable for 5> Can run up to 15 //(max courant number)
//maxDeltaT 0.01;
// ************************************************************************* //
functions
{
/*
*/
};
checkMesh
Code:
Create time
--> FOAM Warning :
From function static Foam::instantList Foam::timeSelector::select0(Foam::Time&, const Foam::argList&)
in file db/Time/timeSelector.C at line 269
No time specified or available, selecting 'constant'
Create polyMesh for time = constant
Reconstructing and writing vtk representation of all faceSets and cellSets.
Time = 0s
Mesh stats
points: 1281656
faces: 3306166
internal faces: 3066597
cells: 1012712
faces per cell: 6.2927693
boundary patches: 6
point zones: 0
face zones: 1
cell zones: 1
Overall number of cells of each type:
hexahedra: 899961
prisms: 10677
wedges: 0
pyramids: 0
tet wedges: 11
tetrahedra: 0
polyhedra: 102063
Breakdown of polyhedra by number of faces:
faces number of cells
4 81
5 73
6 8183
7 235
8 43
9 84928
12 8122
15 396
18 2
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
outercyl 22402 23515 ok (non-closed singly connected)
inlet 3575 4084 ok (non-closed singly connected)
outlet 3529 4056 ok (non-closed singly connected)
leftperiodic 102419 103147 ok (non-closed singly connected)
rightperiodic 101922 102647 ok (non-closed singly connected)
blade 5722 5762 ok (non-closed singly connected)
Checking geometry...
Overall domain bounding box (-16.549999 -16.859896 -4.8514801e-05) (17.15 25.150302 16.849342)
Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
Mesh has 3 solution (non-empty) directions (1 1 1)
Boundary openness (5.2536918e-16 -1.0812025e-16 5.2740338e-13) OK.
Max cell openness = 3.7686215e-16 OK.
Max aspect ratio = 10.863457 OK.
Minimum face area = 5.622304e-06. Maximum face area = 0.55301834. Face area magnitudes OK.
Min volume = 1.7716102e-06. Max volume = 0.37326841. Total volume = 18720.801. Cell volumes OK.
Mesh non-orthogonality Max: 71.724563 average: 9.9607327
*Number of severely non-orthogonal (> 70 degrees) faces: 1.
Non-orthogonality check OK.
<<Writing 1 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
Max skewness = 1.9453286 OK.
Coupled point location match (average 0) OK.
Mesh OK.
End
boundary
Code:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volScalarField;
object nuSgs;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -1 0 0 0 0];
internalField uniform 1.5e-3;
boundaryField
{
AMI1
{
type cyclicAMI;
// value $internalField;
}
AMI2
{
type cyclicAMI;
// value $internalField;
}
blade
{
type nutUSpaldingWallFunction;
value $internalField;
}
outercyl
{
type slip;
}
///*
"(leftperiodic|rightperiodic)"
{
type cyclicAMI;
}
// */
inlet
{
type calculated;
value $internalField;
}
outlet
{
type calculated;
value $internalField;
}
}
// ************************************************************************* //
p
Code:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 7 0);
boundaryField
{
AMI1
{
type cyclicAMI;
// value $internalField;
}
AMI2
{
type cyclicAMI;
// value $internalField;
}
blade
{
type movingWallVelocity;
value uniform (0 0 0);
}
outercyl
{
type slip;
}
///*
"(leftperiodic|rightperiodic)"
{
type cyclicAMI;
}
//*/
inlet
{
type fixedValue;
value $internalField;
}
outlet
{
type zeroGradient;
}
}
u
Code:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 7 0);
boundaryField
{
AMI1
{
type cyclicAMI;
// value $internalField;
}
AMI2
{
type cyclicAMI;
// value $internalField;
}
blade
{
type movingWallVelocity;
value uniform (0 0 0);
}
outercyl
{
type slip;
}
///*
"(leftperiodic|rightperiodic)"
{
type cyclicAMI;
}
//*/
inlet
{
type fixedValue;
value $internalField;
}
outlet
{
type zeroGradient;
}
}
dynamicMeshDict
Code:
mover
{
type motionSolver;
libs ("libfvMeshMovers.so" "librigidBodyMeshMotion.so");
motionSolver solidBody;
cellZone cell_inner_volume;
solidBodyMotionFunction rotatingMotion;
nonRotatingPatches (leftperiodic rightperiodic);
origin (0.3 0 0);
axis (0 1 0);
omega constant -7.54;
}
createBaffles
Code:
FoamFile
{
format ascii;
class dictionary;
object createBafflesDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Whether to convert internal faces only (so leave boundary faces intact).
// This is only relevant if your face selection type can pick up boundary
// faces.
internalFacesOnly true;
// Baffles to create.
baffles
{
rotating
{
//- Use predefined faceZone to select faces and orientation.
type faceZone;
zoneName face_inner_volume;
//zoneName cell_inner_volume;
//zoneName inner_volume;
patches
{
master
{
//- Master side patch
name AMI1;
type cyclicAMI;
matchTolerance 0.5;// 0.0001
neighbourPatch AMI2;
//transform noOrdering; //v2112 <
transform none; //v10
}
slave
{
//- Slave side patch
name AMI2;
type cyclicAMI;
matchTolerance 0.5;
neighbourPatch AMI1;
//transform noOrdering; //v2112 <
transform none; //v10
}
}
}
}
