# Weird AMI Courant Number

 Register Blogs Members List Search Today's Posts Mark Forums Read

 March 16, 2016, 14:02 Weird AMI Courant Number #1 New Member     Vyssion Join Date: Mar 2016 Posts: 12 Rep Power: 10 Hello, I am trying to simulate a fan rotating within a fixed duct just floating in space. As it is currently set up, the fan is 1m in diameter and sits in the centre of a 4x4x4m cube domain. I am using snappyHexMesh and pimpleDyMFoam to try and simulate this rotation, and seem to have it working for the most part... until about 50 ish iterations in, when the time step begins to approach zero (i.e. 1e-13) and the maximum courant number increases above the defined maximum I have set. (I have said "weird" in the title because the average courant number that I have is on the order of 1e-05 give or take)... At this point, the simulation hits a floating point error and terminates. I am experienced in commercial packages for CFD, but have never used OpenFOAM prior to a fortnight ago, and so am looking for some help to get it running smoothly. I should add that I also have tried a simple MRF setup (which I have managed to converge) however, what I dont like about this option is that I have to state the Cp and Ct values within the MRFProperties script when I am looking to measure the "lift" generated just by rotating some .stl geometry. Here is the error I receive: Code: ```Courant Number mean: 9.258629e-06 max: 1.998997 deltaT = 1.771031e-06 Time = 5.24984e-05 solidBodyMotionFunctions::rotatingMotion::transformation(): Time = 5.249842e-05 transformation: ((0 0 0) (0.9999862 (0 -0.005249818 0))) PIMPLE: iteration 1 smoothSolver: Solving for Ux, Initial residual = 0.0005662521, Final residual = 6.517957e-13, No Iterations 2 smoothSolver: Solving for Uy, Initial residual = 5.593663e-05, Final residual = 3.318731e-14, No Iterations 2 smoothSolver: Solving for Uz, Initial residual = 0.001630861, Final residual = 1.829808e-12, No Iterations 2 GAMG: Solving for p, Initial residual = 0.8991979, Final residual = 0.008066495, No Iterations 1 [5] #0 Foam::error::printStack(Foam::Ostream&) at ??:? [5] #1 Foam::sigFpe::sigHandler(int) at ??:? [5] #2 ? in "/lib/x86_64-linux-gnu/libc.so.6" [5] #3 Foam::GAMGSolver::scale(Foam::Field&, Foam::Field&, Foam::lduMatrix const&, Foam::FieldField const&, Foam::UPtrList const&, Foam::Field const&, unsigned char) const at ??:? [5] #4 Foam::GAMGSolver::Vcycle(Foam::PtrList const&, Foam::Field&, Foam::Field const&, Foam::Field&, Foam::Field&, Foam::Field&, Foam::Field&, Foam::Field&, Foam::PtrList >&, Foam::PtrList >&, unsigned char) const at ??:? [5] #5 Foam::GAMGSolver::solve(Foam::Field&, Foam::Field const&, unsigned char) const at ??:? [5] #6 Foam::fvMatrix::solveSegregated(Foam::dictionary const&) at ??:? [5] #7 Foam::fvMatrix::solve(Foam::dictionary const&) at ??:? [5] #8 ? at ??:? [5] #9 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6" [5] #10 ? at ??:? [vmware-vm:07270] *** Process received signal *** [vmware-vm:07270] Signal: Floating point exception (8) [vmware-vm:07270] Signal code: (-6) [vmware-vm:07270] Failing at address: 0x3e800001c66 [vmware-vm:07270] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x352f0) [0x7fa0de9382f0] [vmware-vm:07270] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x37) [0x7fa0de938267] [vmware-vm:07270] [ 2] /lib/x86_64-linux-gnu/libc.so.6(+0x352f0) [0x7fa0de9382f0] [vmware-vm:07270] [ 3] /opt/openfoam30/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver5scaleERNS_5FieldIdEES3_RKNS_9lduMatrixERKNS_10FieldFieldIS1_dEERKNS_8UPtrListIKNS_17lduInterfaceFieldEEERKS2_h+0xa2) [0x7fa0dfb901b2] [vmware-vm:07270] [ 4] /opt/openfoam30/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver6VcycleERKNS_7PtrListINS_9lduMatrix8smootherEEERNS_5FieldIdEERKS8_S9_S9_S9_S9_S9_RNS1_IS8_EESD_h+0x783) [0x7fa0dfb93073] [vmware-vm:07270] [ 5] /opt/openfoam30/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver5solveERNS_5FieldIdEERKS2_h+0x671) [0x7fa0dfb95d21] [vmware-vm:07270] [ 6] /opt/openfoam30/platforms/linux64GccDPInt32Opt/lib/libfiniteVolume.so(_ZN4Foam8fvMatrixIdE15solveSegregatedERKNS_10dictionaryE+0x16c) [0x7fa0e257e1ec] [vmware-vm:07270] [ 7] pimpleDyMFoam(_ZN4Foam8fvMatrixIdE5solveERKNS_10dictionaryE+0x199) [0x4575c9] [vmware-vm:07270] [ 8] pimpleDyMFoam() [0x41ee71] [vmware-vm:07270] [ 9] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf0) [0x7fa0de923a40] [vmware-vm:07270] [10] pimpleDyMFoam() [0x422529] [vmware-vm:07270] *** End of error message *** -------------------------------------------------------------------------- mpirun noticed that process rank 5 with PID 7270 on node vmware-vm exited on signal 8 (Floating point exception).``` Here is my snappyHexMesh.dict Code: ```// Which of the steps to run castellatedMesh true; snap true; addLayers true; // GEOMETRY - DEFINITION OF ALL SURFACES geometry { fan_blades.stl // Filename { type triSurfaceMesh; // Type name fan_blades; // Name to be referred to as regions { fan_blades // Name of Solid which is inside the .stl { name fan_blades; } } } duct.stl // Filename { type triSurfaceMesh; // Type name duct; // Name to be referred to as regions { duct // Name of Solid which is inside the .stl { name duct; } } } refinementCylOuter // User Defined Region { type searchableCylinder; point1 (0 -0.15 0); // Base point2 (0 0.3 0); // Top radius 0.65; } refinementCylinder // User Defined Region { type searchableCylinder; point1 (0 -0.1 0); // Base point2 (0 0.25 0); // Top radius 0.5; } rotatingCylinder // User Defined Region { type searchableCylinder; point1 (0 -0.05 0); point2 (0 0.15 0); radius 0.5; name rotatingCylinder; } }; // SETTINGS FOR castellatedMeshControls GENERATION castellatedMeshControls { // REFINEMENT PARAMETERS // If local number of cells is >= maxLocalCells on any processor switches from // refinement followed by balancing (current method) to (weighted) balancing // before refinement. maxLocalCells 200000; // Overall cell limit (approximately). Refinement will stop immediately // upon reaching this number. // NOTE: that this is the number of cells before removing the part which // is not 'visible' from the keepPoint. The final number of cells might // actually be a lot less. maxGlobalCells 15000000; // The surface refinement loop might spend lots of iterations refining just a // few cells. This setting will cause refinement to stop if <= minimumRefine // are selected for refinement. Note: it will at least do one iteration // (unless the number of cells to refine is 0) minRefinementCells 5; // Allow a certain level of imbalance during refining // (since balancing is quite expensive) // Expressed as fraction of perfect balance (= overall number of cells / // nProcs). 0=balance always. maxLoadUnbalance 0.10; // Number of buffer layers between different levels. // 1 means normal 2:1 refinement restriction, larger means slower refinement. nCellsBetweenLevels 3; // EXPLICIT FEATURE EDGE REFINEMENT // Specifies a level for any cell intersected by its edges. // This is a featureEdgeMesh, read from constant/triSurface for now. features ( { file "fan_blades.eMesh"; // File Containing Mesh level 9; // Level of Refinement } { file "duct.eMesh"; // File Containing Mesh level 8; // Level of Refinement } ); // SURFACE BASED REFINEMENT // Specifies two levels for every surface. The first is the minimum level, // every cell intersecting a surface gets refined up to the minimum level. // The second level is the maximum level. Cells that 'see' multiple // intersections where the intersections make an // angle > resolveFeatureAngle get refined up to the maximum level. // DO NOT PUT THE FILE NAME HERE!!!!!! // Inside the .stl file there will be a solid defined with a given name // THAT is the string that you type here... refinementSurfaces { fan_blades { level (6 8); // ( ) Surface Refinement patchInfo // Optional specification of patch type { // No constraint types (e.g. symm) allowed. type wall; } } duct { level (5 8); patchInfo { type wall; } } rotatingCylinder { level (6 8); faceZone rotatingCylinder; cellZone rotatingCylinder; cellZoneInside inside; } } // RESOLVE SHARP ANGLES // Cells above angle use the max ref level, cells below use min resolveFeatureAngle 30; // REGION-WISE REFINEMENT // Specifies refinement level for cells in relation to a surface. One of // three modes // 1. distance 'levels' specifies per distance to the surface the // wanted refinement level. The distances need to be // specified in descending order. // 2. inside 'levels' is only one entry and only the level is used. // All cells inside the surface get refined up to the level. // The surface needs to be closed for this to be possible. // 3. outside Same as inside, but cells outside. refinementRegions // See "geometry" sub-Dictionary for { // User defined refinement regions refinementCylOuter { mode inside; // inside, outside, distance? levels ((1e15 6)); // FORMAT ( ) // NOTE: for "inside" and "outside", // is ignored, but MUST be specified!! } refinementCylinder { mode inside; levels ((1e15 7)); } fan_blades { mode distance; // Refinement level levels ((0.05 8)); // within of .stl } duct { mode distance; levels ((0.02 7)); } rotatingCylinder { mode inside; levels ((1e15 7)); } } // MESH SELECTION // After refinement patches get added for all refinementSurfaces and // all cells intersecting the surfaces get put into these patches. The // section reachable from the locationInMesh is kept. // NOTE: This point should never be on a face, always inside a cell, // even after refinement. locationInMesh (3.9999 3.9999 3.9999); // Whether any faceZones (as specified in the refinementSurfaces) // are only on the boundary of corresponding cellZones or also allow // free-standing zone faces. Not used if there are no faceZones. allowFreeStandingZoneFaces true; } // SETTINGS FOR THE SNAPPING PROCESS snapControls { // Number of patch smoothing iterations before finding correspondence to surface nSmoothPatch 3; // Relative distance for points to be attracted by surface feature point // or edge. True distance is this factor times local // maximum edge length. // Tolerance*cellsize on surface = distance to look for closest cell // distance relative to local cell size to look for adjacent cells tolerance 4.0; // Number of mesh displacement relaxation iterations. // Deforms entire mesh to make good quality iterations nSolveIter 30; // Maximum number of snapping relaxation iterations. Should stop before upon // reaching a correct mesh. Relaxing sharp intersection makes smooth. nRelaxIter 5; // FEATURE SNAPPING // Number of feature edge snapping iterations. Leave out altogether to disable. // TRY WITH THIS DISABLED?! // Do not use here since mesh resolution too low and baffles present nFeatureSnapIter 10; // Detect (geometric only) features by sampling surface (default=false) implicitFeatureSnap true; // Use castellatedMeshControls::features (default = true) explicitFeatureSnap true; // Detect points on multiple surfaces (only for explicitFeatureSnap) multiRegionFeatureSnap true; } // SETTINGS FOR THE LAYER ADDITION PROCESS addLayersControls { // Are the thickness parameters below relative to the undistorted // size of the refined cell outside layer (true) or absolute sizes (false). // Relative to local cell size at surface relativeSizes true; layers { fan_blades { nSurfaceLayers 5; } duct { nSurfaceLayers 5; } } // Expansion factor for layer mesh expansionRatio 1.2; // Wanted thickness of final added cell layer. If multiple layers // is the thickness of the layer furthest away from the wall. // Relative to undistorted size of cell outside layer. // See relativeSizes parameter. // If "relativeSizes" parameter is false you need to specify in meters the // finallayerthickness, with it on we just use a factor of the local cell size // Per final patch (so not geometry!) the layer information finalLayerThickness 0.2; // Minimum thickness of cell layer. If for any reason layer cannot be above // minThickness do not add layer. Relative to undistorted size of cell outside // layer. Not enforced but seeks to meet this. minThickness 0.01; // If points get not extruded do nGrow layers of connected faces that are // also not grown. This helps convergence of the layer addition process // close to features. // Note: changed(corrected) w.r.t 17x! (didn't do anything in 17x) nGrow 0; // Sometimes cells don't extrude away from surface due to cells meeting at a // sharp angle this forces the cells to grow a layer away from the surface. // ADVANCED SETTINGS // When not to extrude surface. 0 is flat surface, 90 is when two faces // are perpendicular featureAngle 180; // At non-patched sides allow mesh to slip if extrusion direction makes // angle larger than slipFeatureAngle. slipFeatureAngle 359; // Maximum number of snapping relaxation iterations. Should stop // before upon reaching a correct mesh. nRelaxIter 5; // Number of smoothing iterations of surface normals nSmoothSurfaceNormals 3; // Number of smoothing iterations of interior mesh movement direction nSmoothNormals 3; // Smooth layer thickness over surface patches nSmoothThickness 10; // Stop layer growth on highly warped cells maxFaceThicknessRatio 0.3; // Reduce layer growth where ratio thickness to medial distance is large maxThicknessToMedialRatio 0.5; // Angle used to pick up medial axis points // NOTE: changed(corrected) w.r.t 17x! 90 degrees corresponds to 130 in 17x. minMedianAxisAngle 80; // Create buffer region for new layer terminations nBufferCellsNoExtrude 0; // Overall max number of layer addition iterations. The mesher will exit // if it reaches this number of iterations; possibly with an illegal mesh nLayerIter 30; } // GENERIC MESH QUALITY SETTINGS // At any undoable phase these determine where to undo. meshQualityControls { #include "meshQualityDict" // ADVANCED QUALITY CONTROLS // Number of error distribution iterations nSmoothScale 4; // Amount to scale back displacement at error points errorReduction 0.75; } // ADVANCED snappyHexMesh PARAMETERS // Flags for optional output // 0 Only write final meshes // 1 Write intermediate meshes // 2 Write volScalarField with cellLevel for postprocessing // 4 Write current intersections as .obj files debug 0; writeFlags ( scalarLevels layerSets layerFields // Write volScalarField for layer coverage ); // MERGE TOLERANCE // Is a fraction of overall bounding box of initial mesh. // NOTE: the write tolerance needs to be higher than this. mergeTolerance 1.001e-6; //``` Here is my fvSchemes script: Code: ```ddtSchemes { default Euler; } gradSchemes { default Gauss linear; grad(p) Gauss linear; grad(U) cellLimited Gauss linear 1; } divSchemes { default none; div(phi,U) Gauss linearUpwind grad(U); div(phi,k) Gauss upwind; div(phi,omega) Gauss upwind; div((nuEff*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear limited corrected 0.33; // ???? } interpolationSchemes { default linear; } snGradSchemes { default limited corrected 0.33; // ???? } wallDist { method meshWave; }``` Here is my fvSolution script: Code: ```solvers { pcorr { solver GAMG; tolerance 1e-03; relTol 0.0; smoother DICGaussSeidel; nPreSweeps 0; nPostSweeps 2; cacheAgglomeration false; // ???? nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; } p { \$pcorr; tolerance 1e-06; relTol 0.1; } pFinal { \$p; tolerance 1e-08; relTol 0.1; } U { solver smoothSolver; smoother symGaussSeidel; nSweeps 2; tolerance 1e-08; relTol 0.1; } UFinal { \$U; } k { solver smoothSolver; smoother symGaussSeidel; nSweeps 2; tolerance 1e-08; relTol 0.1; } kFinal { \$k; } omega { solver smoothSolver; smoother symGaussSeidel; nSweeps 2; tolerance 1e-08; relTol 0.1; } omegaFinal { \$omega; } } PIMPLE { correctPhi no; nOuterCorrectors 2; nCorrectors 1; nNonOrthogonalCorrectors 1; } relaxationFactors { fields { p 0.5; } equations { U 0.5; k 0.5; omega 0.5; } } cache { grad(U); }``` Here is my controlDict script: Code: ```application pimpleDyMFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 0.1; deltaT 1e-5; writeControl adjustableRunTime; adjustTimeStep true; writeInterval 0.001; purgeWrite 0; ////- For testing with moveDynamicMesh //deltaT 0.01; //writeControl timeStep; //writeInterval 1; writeFormat ascii; writePrecision 7; writeCompression compressed; timeFormat general; timePrecision 6; runTimeModifiable true; maxCo 2; functions { #include "readFields" #include "Q" // #include "surfaces" #include "forces" }``` Now it seems that the issue is arrising for the Pressure script within the "0" directory and so I have also included my "p" file and "omega" since I am using the k-omega SST turbulence model: Initial Conditions (referenced in the "p" and "omega" scripts (yes I know that my initial omega was set to 500 -- I was getting immediate divergence with it set to 0.1 and read that setting it much higher would then fix itself as it iterated): Code: ```turbulentKE 0.1; // kInlet Value turbulentOmega 500.0; //``` "p" file: Code: ```dimensions [0 2 -2 0 0 0 0]; internalField uniform 0; boundaryField { top // Face 5 { type zeroGradient; } bottom // Face 4 { type fixedValue; value uniform 0; } fan_blades { type zeroGradient; } duct { type zeroGradient; } #include "caseDicts/setConstraintTypes" #include "include/frontBackLeftRightPatches" }``` And "omega" file: Code: ```#include "include/initialConditions" dimensions [0 0 -1 0 0 0 0]; internalField uniform \$turbulentOmega; boundaryField { //- Set patchGroups for constraint patches // #includeEtc "caseDicts/setConstraintTypes" top { type fixedValue; value uniform \$turbulentOmega; } bottom { type inletOutlet; inletValue uniform \$turbulentOmega; value uniform \$turbulentOmega; } fan_blades { type omegaWallFunction; value \$internalField; } duct { type omegaWallFunction; value uniform \$turbulentOmega; } AMI1 { type cyclicAMI; } AMI2 { type cyclicAMI; } #include "include/frontBackLeftRightPatches" }``` Any help at all would be welcome!! Thanks __________________ Optimism, pessimism, f*ck that: We're going make it happen. As God is my bloody witness, I'm hell-bent on making it work. - Elon Musk in 2008, after three unsuccessful Falcon launches.

 March 18, 2016, 03:53 SOLVED - Though next error is present... #2 New Member     Vyssion Join Date: Mar 2016 Posts: 12 Rep Power: 10 I looked into the code more closely and found the issue. In my snappyHexMesh script, I forgot to define the faceType of the rotatingCylinder at all. Once I set that to "boundary", the rest of it executed fine. However... I am now at the point of attempting to "moveDynamicMesh -checkAMI" and am finding that I am getting aspect ratios on the order of >1e+94, in addition to a lot of "bad" (collectively") cells including negative cell volumes. I have heard that this issue arises due to the way the meshing traverses the AMI boundary and that one possible way to solve this issue is to mesh inside the AMI first, then mesh the rest of the domain and use a "combineMeshes" function (which makes sense seeing as the meshes should end perfectly at the boundaries). This may also explain why I was able to get about 50 iterations out of the pimpleDyMFoam before things blew up. I will try meshing separately and combining and let you know how it goes. One final thing though; before I defined the faceType of the AMI as a boundary, I was getting fatal errors stating that "cannot find patchField entry for " for p, U, etc. Once I defined it as a boundary, these errors disappeared... mostly. I now receive this error: Code: ```Create time Create mesh for time = 0 Selecting dynamicFvMesh solidBodyMotionFvMesh Selecting solid-body motion function rotatingMotion Applying solid body motion to cellZone rotatingCylinder PIMPLE: no residual control data found. Calculations will employ 2 corrector loops Reading field p Reading field U Reading/calculating face flux field phi Selecting incompressible transport model Newtonian Selecting turbulence model type RAS Selecting RAS turbulence model kOmegaSST Selecting patchDistMethod meshWave --> FOAM FATAL IO ERROR: Cannot find patchField entry for AMI1 file: /home/osboxes/OpenFOAM/osboxes-3.0.1/run/base_model_bestMesh_AMI/0/omega.boundaryField from line 30 to line 27. From function GeometricField::GeometricBoundaryField::readField(const DimensionedField&, const dictionary&) in file /home/openfoam/OpenFOAM/OpenFOAM-3.0.1/src/OpenFOAM/lnInclude/GeometricBoundaryField.C at line 209. FOAM exiting``` Any suggestions as to why this is? I will keep poking around. __________________ Optimism, pessimism, f*ck that: We're going make it happen. As God is my bloody witness, I'm hell-bent on making it work. - Elon Musk in 2008, after three unsuccessful Falcon launches.

 March 23, 2016, 04:35 Update #3 New Member     Vyssion Join Date: Mar 2016 Posts: 12 Rep Power: 10 Figured I would keep updating this so that anyone in the future who has the same issue can have some ideas to try. I have found that when I performed the "moveDynamicMesh -checkAMI" operation, I was getting meshing errors build up after about 4 time steps (of 1e-05). Poking around, it seems to be the way that snappyHexMesh works that is the issue. Because I defined the AMI as a boundary, it will mesh through it without any consideration that it will be rotating. Hence, when it does eventually rotate, the cells are dragged along with it. Playing around with some of the settings in snappyHexMesh, (such as keeping the min-max cell sizes at the interface the same, along with any refinements which overlap the AMI zone being limited in the same way), I have managed to get it to get it to run until about 9 time steps before the -checkAMI blows up and I can simulate through this until 0.00018 seconds. Considering I need about 0.032 seconds for a full revolution, I am a ways off!! I need to look into meshing separately and then merging them... Problem is, I have no clue of the nitty-gritty which I would need to do in order to set this up... Has anyone ever done this before or knows of a tutorial that would help?? P.S. I am wondering if part of the issue is that my fan blades touch the duct itself (Im not interested in capturing tip clearance shear interactions at this point) and so right at the AMI boundary, there are zero cells at some points... Could this be an issue?? hogsonik likes this. __________________ Optimism, pessimism, f*ck that: We're going make it happen. As God is my bloody witness, I'm hell-bent on making it work. - Elon Musk in 2008, after three unsuccessful Falcon launches.

 April 13, 2016, 02:31 #4 New Member   Join Date: Mar 2016 Posts: 3 Rep Power: 10 I ran the similar case(Re=1.4*10^5) with the standard k-epsilon turbulence model and I got a higher St=0.25~0.26. I'm confused about y+. I think that the min y+ should be higher than 30 when using wallfunction to make the first point near the cylinder wall locate in the logarithmic layer. And I realize that y+ is different with that when calculated in initial situation. Which y+ should be based on to decide the distance of the first point to the wall. Can you give me some opinions? Thanks a lot!

 Tags ami, courant, cyclic, floating point, mrf