ipython markdown

http://daringfireball.net/projects/markdown/basics

http://notes.zatvia.com/notes-md.html

Four tildes to set verbatim on and off.

Circular cylinder howto

This page describes one of my efforts to simulate the flow around a circular cylinder with openFoam. Comments are most welcome. This is not all my own work; I've used several articles from foamwiki and others but this is a work in progress. I'll credit eventually...

The process here is:
1) Specify the mesh.
2) Set flow / material parameters.
3) Solution control:
3a) Set turbulence modelling approach.
3b) Specify solution schemes / types.
3c) Specify solution parameters.
3d) Specify software settings (save rate, start time etc)
4) Create the mesh.
5) Solve.

The turbulent analysis benefits from a laminar pre-solution to set the initial field for the turbulent analysis and progressively complex turbulence models benefit from progressively more accurate initial conditions. For this reason, a Reynolds Stress analysis may benefit from an initial laminar solution, then say a k-epsilon solution to generate the initial field. That's another post, anyway. This post runs from scratch with turbulence.

1) Create the mesh: Convergence studies suggested little change was occurring with further refinement to the boundary layer mesh. The file,

constant/polyMesh/blockMeshDict

FoamFile{version 2.0; format ascii; class dictionary; object blockMeshDict;}

convertToMeters 0.05;

// Vertex dimensions are in diameters!

// space between `vertices' and `(' is required.
vertices (  
   (0.5 0 0)                // cylinder wall
   (1 0 0)                  // b/l
   (20 0 0)                 // exit
   (20 0.707107 0)
   (0.707107 0.707107 0)
   (0.353553 0.353553 0)
   (20 10 0)
   (0.707107 10 0)
   (0 10 0)
   (0 1 0)
   (0 0.5 0)
   (0.5 0 0.5)
   (1 0 0.5)
   (20 0 0.5)               // downstream exit at 2
   (20 0.707107 0.5)
   (0.707107 0.707107 0.5)
   (0.353553 0.353553 0.5)
   (20 10 0.5)
   (0.707107 10 0.5)
   (0 10 0.5)
   (0 1 0.5)

   (0 0.5 0.5)

   (-0.5 0 0)
   (-1 0 0)
   (-10 0 0)
   (-10 0.707107 0)
   (-0.707107 0.707107 0)
   (-0.353553 0.353553 0)
   (-10 10 0)
   (-0.707107 10 0)
   (-0.5 0 0.5)
   (-1 0 0.5)
   (-10 0 0.5)
   (-10 0.707107 0.5)
   (-0.707107 0.707107 0.5)
   (-0.353553 0.353553 0.5)
   (-10 10 0.5)
   (-0.707107 10 0.5)
   );


// Defined anti-clockwise on front face, then back face.
blocks(
   // cylinder wall, 0 to 45 deg.  radial, circumferential.
   hex(5 4 9 10 16 15 20 21) (80 10 1) simpleGrading (10 1 1) 

   // cylinder wall, 45 to 90 deg. radial, circumferential.
   hex(0 1 4 5 11 12 15 16) (80 10 1) simpleGrading (10 1 1)  

   // downstream, behind cylinder.
   hex(1 2 3 4 12 13 14 15) (160 10 1) simpleGrading (15 1 1) 
   
   // downstream, upper 1/4
   hex(4 3 6 7 15 14 17 18) (160 60 1) simpleGrading (15 14 1)

   // above cylinder, 45 to 90 deg. 
   hex(9 4 7 8 20 15 18 19) (10 60 1) simpleGrading (1 14 1)  

   // cylinder wall, 90 to 135 deg. radial, circumferential
   hex(10 9 26 27 21 20 34 35)  (80 10 1) simpleGrading (10 1 1) 

   // cylinder wall, 135 to 180 deg. radial, circ.
   hex(27 26 23 22 35 34 31 30) (80 10 1) simpleGrading (10 1 1) 

   // above cylinder, 90 to 135 deg.
   hex(23 26 25 24 31 34 33 32) (10 60 1) simpleGrading (1 14 1) 

   // upstream, upper 1/4
   hex(26 29 28 25 34 37 36 33) (60 60 1) simpleGrading (14 14 1) 

   // directly upstream of cylinder. 
   hex(26 9 8 29 34 20 19 37)   (10 60 1) simpleGrading (1 14 1)  
   );

edges(
   arc 0 5 (0.469846 0.17101 0)
   arc 5 10 (0.17101 0.469846 0)
   arc 1 4 (0.939693 0.34202 0)
   arc 4 9 (0.34202 0.939693 0)
   arc 11 16 (0.469846 0.17101 0.5)
   arc 16 21 (0.17101 0.469846 0.5)
   arc 12 15 (0.939693 0.34202 0.5)
   arc 15 20 (0.34202 0.939693 0.5)

   arc 22 27 (-0.469846 0.17101 0)
   arc 27 10 (-0.17101 0.469846 0)
   arc 23 26 (-0.939693 0.34202 0)
   arc 26 9 (-0.34202 0.939693 0)
   arc 30 35 (-0.469846 0.17101 0.5)
   arc 35 21 (-0.17101 0.469846 0.5)
   arc 31 34 (-0.939693 0.34202 0.5)
   arc 34 20 (-0.34202 0.939693 0.5) );

// patches ( type name (nodes)
// These are GEOMETRY effect, not flow etc.
// So flow inlet is a patch in this file,
//  while physical walls are wall.
patches(
   patch inlet(( 28 25 33 36 )
               ( 25 24 32 33 ))

   patch outlet((6 3 14 17)
                (3 2 13 14))

   symmetry symmetry((24 23 31 32)
                     (23 22 30 31)
                     (0 1 12 11)
                     (1 2 13 12))

   wall top((28 29 37 36)
            (29 8 19 37)
            (8 7 18 19)
            (7 6 17 18))

   wall wall((22 27 35 30)
             (27 10 21 35)
             (10 5 16 21)
             (5 0 11 16))

   empty emptyf((28 29 26 25)
                (29 8 9 26)
                (26 9 10 27)
                (26 27 22 23)
                (25 26 23 24)
                (8 7 4 9)
                (7 6 3 4)
                (9 4 5 10)
                (4 3 2 1)
                (4 1 0 5))

   empty emptyb((36 37 34 33)
                (37 19 20 34)
                (34 20 21 35)
                (34 35 30 31)
                (33 34 31 32)
                (19 18 15 20)
                (18 17 14 15)
                (20 15 16 21)
                (15 14 13 12)
                (15 12 11 16)));
mergePatchPairs();

This file creates some skewed cells; these will require

nNonOrthogonalCorrectors = 4

or greater. This was found experimentally. This parameter may effect the number of iterations of the pressure solver, per `time-step'.

2) Specify the fluid properties as kinematic viscosity. Specify nu where nu = mu / rho with mu, the viscosity and rho the density. mu = 1.78E-5 kg/(m s) = 1.78E-5 Pa s at 15 degrees centigrade. Air density is 1.225 kg/m3 at this temperature. nu therefore is, 1.78E-5 kg/(ms) / 1.225 kg/m3 = 1.45E-5 m2 s-1. (recall that viscosity, mu sets the proportionality between the force required to move two parallel, massless plates with fluid between them; and the plate area, their relative velocity and separation. F = mu A u/y (well, du/dy). Its units are prescribed by 1/mu = A u / (y F) -> 1/(m2 m s-1 m-1 N-1) = m-2 s N = Pa s). Of similar interest for this study, the Reynolds number Re = u D rho / mu = u D / nu with velocity, u; and diameter, D. For a 0.05m cylinder in 10m/s flow, the Reynolds number is 10m/s 0.05m 1.225 kg m-3 / (1.78E-5 kg m-1 s-1) = 34400. This cylinder will display transition in shear layer behaviour (Zdravkovitch) so the boundary layer will be laminar when it separates and separation will occur around 80 degrees.

constant/transportProperties

FoamFile{version 2.0; format ascii; class dictionary;object transportProperties;}
transportModel Newtonian;
nu nu [0 2 -1 0 0 0 0] 1.45E-5;

3) Solution control. 3a) Set the turbulence treatment for RAS. Create

constant/RASProperties

Set as laminar for now (turbulence off):

FoamFile{version 2.0; format ascii; class dictionary; object RASProperties;}
RASModel LaunderGibsonRSTM;
turbulence off;

3b) Solution types / schemes are specified in

 system/fvSchemes 

FoamFile{version 2.0; format ascii; class dictionary; location "system"; object fvSchemes;}

snGradSchemes{
    // limited 0.7; for bad meshes, corrected for good meshes.
    // must fit to laplacianSchemes
    default         limited 0.7;                                
    }

laplacianSchemes{
   default none;  // If something is missing, tell me.

   // for good meshes
   //laplacian(nuEff,U) Gauss linear corrected;    
   // for bad meshes  
   laplacian(nuEff,U) Gauss linear limited 0.7;    

   // for good meshes
   //laplacian((1|A(U)),p) Gauss linear corrected; 
   // for bad meshes 
   laplacian((1|A(U)),p) Gauss linear limited 0.7; 
   
   // for turbulent flow involving k only:
   // laplacian(DkEff,k) Gauss linear corrected;

   // for turbulent flow involving epsilon only:
   laplacian(DepsilonEff,epsilon) Gauss linear corrected;

   // for turbulent flow including Reynolds stress only:
   // Needed for LRR. corrected or limited.
   laplacian(DREff,R) Gauss linear corrected; 

   // for turbulent flow requiring nuTilda only, ie SA   
   // laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
 
   // These are from $FOAM_TUTORIALS/incompressible/simpleFoam/airFoil2D
   // SA not kE
   // laplacian(nuEff,U) Gauss linear corrected; // got it
   // laplacian((1|A(U)),p) Gauss linear corrected; // got it

   // Why not limited? nuTilda for SA only.
   // laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; 

   // didnt have it.   limited?
   // laplacian(1,p)  Gauss linear corrected;  
   }

fluxRequired{
    default no; 
    p; 
    //phi;
    }

gradSchemes{
   // Gauss linear or leastSquares
   default   Gauss linear;  
   
   grad(p)   faceLimited leastSquares 0 1;   
   
   grad(U)   Gauss linear;             
   }

interpolationSchemes{
   default  linear;     // linear is fine
   U               linear;

   // new tet mesh interpolator. Use with reconCentral in divSchemes.
   // This is not available in vanilla openFoam yet.  There is an extended
   // version.  I'll try compile it in another article.
   // I've read articles where people stated it gave little benefit,
   // but none that stated it helped.
   // interpolate(U) reconCentral phi cellLimited leastSquare 1.0; 
   }

divSchemes{
   default none;
   
   // stable: Gauss upwind
   div(phi,U)      Gauss upwind;               
   // faster linearUpwind Gauss
   // div(phi,U)      Gauss linearUpwind Gauss; 
   
   // new, tets mesh setting:
   // div(phi,U)      Gauss reconCentral cellLimited leastSquares 1.0;    

   // for turbulent flow only
   // div(phi,k)      Gauss upwind;             

   // for turbulent flow only
   div(phi,epsilon) Gauss upwind;              

   // for turbulent flow only
   div(phi,R)      Gauss upwind;               

   // for turbulent flow only
   div(R)          Gauss linear;               
   
   // for SA turbulent flow only.
   // div(phi,nuTilda) Gauss upwind;           

   // this scheme is not used for steady and laminar flow, but
   // it is read by the solver.  (I haven't tested this)
   div((nuEff*dev(grad(U).T()))) Gauss linear; 
                                               
   // not advised, taken from motorbike case
   div((nuEff*dev(T(grad(U))))) Gauss linear; 

   // from airFoil2D:
   // default         none;  // same
   // div(phi,U)      Gauss linearUpwind grad(U);
   
   // set above but this I think is the fast version.  SA only
   // div(phi,nuTilda) Gauss linearUpwind grad(nuTilda);  
   
   // I think this is the same as the term above.
   // div((nuEff*dev(T(grad(U))))) Gauss linear; 
   }

ddtSchemes{
   default steadyState;
   // ddt(epsilon)
   }

3c) Specify the solution parameters. Create

system/fvSolution

FoamFile{version 2.0; format ascii; class dictionary; location "system"; object fvSolution;}

SIMPLE{
   // 0 is fine for perfect mesh; increase to 1 or 2 for non orthogonal mesh
   nNonOrthogonalCorrectors 2;         

   // convergence criteria steady state
   convergence         1e-7;           
   
   // ref cell is 0 in the airfoil tutorial
   pRefCell               0;           
   
   // ref value is 0 in the airfoil case.
   pRefValue              0;           
   residualControl{
      p 1E-5; 
      U 1E-5; 
      k 1E-5; 
      epsilon 1E-5;
      }                               
   }

relaxationFactors{
   // 0.3 is stable, decrease for bad mesh
   p        0.1;                

   // 0.7 is stable, decrease for bad mesh
   U        0.3;                

   // 0.7 is stable, decrease for bad mesh
   //k        0.3;                

   // Not sure yet about R relaxation.
   //R        0.1;

   // Not sure yet about epsilon relaxation.                
   //epsilon  0.3;
   }

solvers{
   // linear equation system solver for U
   U{
      // solver type
      solver          smoothSolver;
      
      // smoother type
      smoother        GaussSeidel;
   
      // solver finishes if either absolute or relative 
      // tolerance is reached. airFoil2D uses 1E-8
      tolerance       1e-06;          
      relTol          0.01;           

      // smoothSolver setting.
      //: airFoil uses 2 sweeps with 0.1 relTol.
      nSweeps         1;
      maxIter         100;
      }

   // linear equation system solver for p
   p{                                 
      // very efficient multigrid solver
      solver          GAMG;           

      // solver finishes if either absolute or relative 
      // tolerance is reached.  
       // airFoil2D uses relTol 0.1, tol 1E-6
      tolerance       1e-07;          
      relTol          0.001;

      // a minimum number of iterations
      minIter         3;              
      
      // limitation of iterions number
      maxIter         100;            

      // setting for GAMG. airFoil2D uses GaussSeidel.
      smoother        DIC;            

      // 1 for p, set to 0 for all other!  airfoil2D uses 0
      nPreSweeps      1;              

      // 2 is fine.  airfoil2D uses 2
      nPostSweeps     2;              

      // 2 is fine
      nFinestSweeps   2;              
      
      scaleCorrection true;
      directSolveCoarsestLevel false;
      // on unless dynamic mesh refinement is used.
      cacheAgglomeration on;          

      // 500 or sqrt(cell count)
      nCellsInCoarsestLevel 1000;

      agglomerator    faceAreaPair;  
      mergeLevels     1;
      }

   //epsilon{solver PBiCG; preconditioner DILU; tolerance 1E-9; relTol 0;}
   //R{solver PBiCG; preconditioner DILU; tolerance 1E-9; relTol 0;}
   }

3d) Create the control dictionary: file system/controlDict

FoamFile{version 2.0; format ascii; class dictionary; object controlDict;}
application simpleFoam;
startFrom latestTime; // latestTime / startTime
startTime 0;
stopAt endTime;  endTime 100000; 
deltaT 1;
writeControl timeStep; writeInterval 300; purgeWrite 0;
writeFormat ascii; writePrecision 6; writeCompression uncompressed;
timeFormat general;timePrecision 6; 
runTimeModifiable yes;

4) Create the mesh by typing,

$ blockMesh

5) Solve with a single processor by typing,

 $ nohup simpleFoam > log &

or in parallel with

$ mpirun -np 4 simpleFoam -parallel > log &

X11 connection rejected because of wrong authentication.

rm -rf ~/.Xauthorit*

gmsh sinusoidal mesh density

Point (1) = {0, 0, 0};  
Point (2) = {1, 0, 0};  
Point (3) = {2, 0, 0};  
Point (5) = {1, 1, 0};  
Point (9) = {1, -1, 0}; 

Line(1) = {1, 2}; 
Line(2) = {2, 3}; 
Line(3) = {3, 5}; 
Line(5) = {5, 1}; 
Line(7) = {5, 2}; 

Line Loop(8) = {1, -7, 5}; Plane Surface (9) = {8}; 
Line Loop(10) = {7, 2, 3}; Plane Surface (11) = {10}; 

Line(12) = {1, 9}; 
Line(13) = {9, 2}; 
Line(14) = {9, 3};

Line Loop(15) = {1, -13, -12}; Plane Surface(16) = {15}; 
Line Loop(17) = {2, -14, 13}; Plane Surface(18) = {17}; 

// Sets the element size in the left region (9) 
Field[1] = MathEval; Field[1].F = "0.011+0.009*Sin(x*12)"; 

Field[6] = Restrict; Field[6].IField = 1; 
Field[6].FacesList = {16}; 
Background Field = 6;

http://engits.eu/wiki/index.php/Tutorials/Version1.2/Unstructured_Grids_for_OpenFOAM_With_Blender_and_enGrid

gmsh min function

minab = (a < b) ? a : b;

ssh security measures

After a failed login attempt, refuse attempts to log in from that ip address for one minute:

$ iptables -A INPUT -p tcp -m state --syn --state NEW --dport 22 -m limit --limit 1/minute --limit-burst 1 -j ACCEPT

$ iptables -A INPUT -p tcp -m state --syn --state NEW --dport 22 -j DROP

Does this need to be re-specified on each restart?

rotating log files

from http://ubuntuforums.org/showthread.php?t=676837

/var/log/news/* {
daily
rotate 7
olddir /var/log/news/oldnews
missingok
postrotate
kill -1 `cat /var/run/syslogd.pid`
endscript
compress
}


I bet this can be updated to rotating archival.....