This tutorial simulates incompressible flow around a cylinder in a 2D rectangular domain.
- Reynolds number: effect of different Re numbers
- Mesh refinement: mesh subdivision
- Domain boundaries: simulation and boundary conditions
- Drag and lift forces: post-processing of drag and lift forces
- Input files: detailed description of case input files
Reynolds number
The case is run first at low Reynolds number and a steady state solution is found.
By increasing velocity the flow becomes unstable and von Karman vortex phenomena are observed in the wake of the cylinder (click here to see the animation )
Mesh refinement
The mesh of the model is structured quadrilateral. DIVISION command is used to refine the mesh in order to obtain more accurate results.
Domain boundaries
The simulation is transient, integrating from 0 to 200 time interval with a prescribed time step of 0.2. Velocity is set to 1, the diameter of the cilynder is 1, density is 1 and dynamic viscosity 0.01, so that the Reynolds number is 100. Units are SI.
Drag and lift forces
As a post process option, forces and moments acting on boundaries are specified to be output in order to calculate drag and lift forces acting on the cylinder surface. Next plots show drag and lift forces due to pressure and viscosity over time.
These plots have been obtained using the Open Source application gnuplot.
-
Type ‘gnuplot’ on a terminal and then:
gnuplot> plot [0:200] [-0.6:0.3] './alya-sets.out' using 1:7 title 'pressure drag force' with lines, './alya-sets.out' using 1:8 title 'pressure lift force' with lines, './alya-sets.out' using 1:4 title 'viscosity drag force' with lines, './alya-sets.out' using 1:5 title 'viscosity lift force' with lines
-
This command produces the plot in the figure by reading the data from columns 1 (time), 4 (F_v_x), 5 (F_v_y), 7 (F_p_x) and 8 (F_p_y) in the file alya-sets.out.
-
Prior to this, alya-sets.out file must be generated with the following command:
/path/Utils/user/alya-sets cylinder-boundary.nsi.set 59 60 61 62
The application alya-sets reads the file cylinder-boundary.nsi.set and does the summation of boundary variables specified in casename.nsi.dat / BOUNDARY_SET, for all elements belonging to boundary sets nº 59, 60, 61 and 62 (cylinder boundary sets in this case). Results are written to file alya-sets.out.
Simulation results show that lift oscillations stabilize at a frequency of 0.14Hz.
Strouhal number (nondimensional frequency) of viscosity and pressure lift oscillations is:
\textnormal{St} = \frac{( \textnormal{f} \cdot \textnormal{d} )}{\textnormal{u}_{\infty}} = 0.14Pressure lift is defined as:
\textnormal{pressure lift} = \frac{\int_{cylinder} {-} pn\mathrm{d}\Gamma}{\rho u_{\infty}^2}where n is the outward normal to the cylinder. Its evolution is plotted in the green curve of the previous figure, showing a clear sinusoidal pattern which is a sign of a sustained vortex shedding process. The simulation gives an amplitude of the pressure lift of 0.26
These results match fairly well with the values reported in the literature:
- G. Houzeaux and J. Principe. A variational subgrid scale model for transient incompressible flows, International Journal of Computational Fluyd Dynamics, Vol. 22, No. 3, March 2008, 135-152.
- C.H.K. Williamson and G.L. Brown. A Series is to represent the Strouhal-Reynolds number relationship of the cylinder wake, J. Fluids Struct. 12,1073 (1998)
Input files
cylinder.dat
$----------------------------------------------------------—
RUN_DATA
ALYA: cylinder $ case name
RUN_TYPE: Preliminary, Frequency=1e6 \ $ initial solution, write restart file frequency
END_RUN_DATA
$----------------------------------------------------------—
PROBLEM_DATA
TIME_COUPLING: Global, Prescribed $ global time step defined here next
TIME_INTERVAL= 0.0, 200.0 $ integration time domain from t=0 to t=200s
TIME_STEP_SIZE= 0.2 $ prescribed integration time step
NUMBER_OF_STEPS= 1e5 $ maximum number of time steps
NASTIN_MODULE: On $ nastin module (Incompressible Navier-Stokes) is used
END_NASTIN_MODULE
PARALL_SERVICE: Off $ parallelization service is not used in this case
END_PARALL_SERVICE
</pre>
END_PROBLEM_DATA
$----------------------------------------------------------—cylinder.dom.dat
$----------------------------------------------------------—
DIMENSIONS
NODAL_POINTS= 1280 $ number of nodes
ELEMENTS= 1200 $ number of elements
SPACE_DIMENSIONS= 2 $ 2D mesh
TYPES_OF_ELEMENTS= 12 $ 4 node quadrilateral elements
BOUNDARIES= 160 $ number of boundary edges
END_DIMENSIONS
$----------------------------------------------------------—
STRATEGY
INTEGRATION_RULE: Open $ open rule is the default
DOMAIN_INTEGRATION_POINTS: 0 $ 0 = automatic, depending on each element type
END_STRATEGY
$----------------------------------------------------------—
GEOMETRY
GROUPS= 20 $ number of groups (for deflation based solvers)
INCLUDE cylinder.geo.dat $ include geometry file (nodes, elements & boundaries)
END_GEOMETRY
$----------------------------------------------------------—
SETS
INCLUDE cylinder.set.dat $ include set file (groups for post-processing)
END_SETS
$----------------------------------------------------------—
BOUNDARY_CONDITIONS, EXTRAPOLATE $ edge BC’s extrapolate to nodes
INCLUDE cylinder.fix.dat $ include boundary conditions file
END_BOUNDARY_CONDITIONS
$----------------------------------------------------------—cylinder.geo.dat
$----------------------------------------------------------—
NODES_PER_ELEMENT
1 4 $ element number, # nodes
2 4
3 4
:
:
1199 4
1200 4
END_NODES_PER_ELEMENT
$----------------------------------------------------------—
ELEMENTS
1 1169 1154 1159 1180 $ element number, node number 1, node number 2, ...
2 1170 1153 1154 1169
3 1171 1152 1153 1170
:
:
1199 791 739 740 792
1200 843 785 739 791
END_ELEMENTS
$----------------------------------------------------------—
COORDINATES
1 1.200000e+01 1.200000e+01 $ node number, coordinate X, coordinate Y
:
:
1279 -5.600000e+00 -1.200000e+01
1280 -6.000000e+00 -1.200000e+01
END_COORDINATES
$----------------------------------------------------------—
BOUNDARIES, ELEMENT
1 790 727 301 $ boundary element #, node #s, element # boundary belongs to
2 727 681 302
3 681 631 303
:
:
159 730 731 1099
160 784 730 1100
END_BOUNDARIES
MATERIALS, NUMBER=0, DEFAULT=1 $ all elements have material 1 END_MATERIALS CHARACTERISTICS END_CHARACTERISTICS
$----------------------------------------------------------—cylinder.set.dat
$----------------------------------------------------------—
ELEMENTS $ element sets definition
1 6 $ element number, set number
2 6
3 6
:
:
:
1199 17
1200 17
END_ELEMENTS
$----------------------------------------------------------—
BOUNDARIES $ boundary sets definition
<pre class="fragment">1 21 $ boundary element number, boundary set number
2 21
3 21
:
:
:
159 62
160 62 $ 1 set for each domain boundary
END_BOUNDARIES
$----------------------------------------------------------—
NODES $ node sets definition (for postprocesing)
389, 270, 113, 611, 220
END_NODES
$----------------------------------------------------------—cylinder.fix.dat
$----------------------------------------------------------—
ON_BOUNDARIES
<pre class="fragment">1 4 $ boundary element number, boundary number
2 4
3 4
:
:
:
159 3
160 3
END_ON_BOUNDARIES
$----------------------------------------------------------—cylinder.ker.dat
$---------------------------------------------------------—
PHYSICAL_PROBLEM
PROPERTIES $ fluid physical properties
DENSITY: CONSTANT, PARAMETERS = 1.0 $ density
VISCOSITY: CONSTANT, PARAMETERS = 0.01 $ dynamic viscosity (μ)
END_PROPERTIES
END_PHYSICAL_PROBLEM
$---------------------------------------------------------—
NUMERICAL_TREATMENT
MESH
DIVISION=0 $ mesh subdivisions (0 = no subdivisions)
END_MESH
ELSEST $ element search strategy
STRATEGY: BIN $ BIN divide mesh into boxes to find elements that host
witness points (suited for structured meshes)
NUMBER: 100,10 $ 2D domain divided in 100 x 10 boxes
END_ELSEST
END_NUMERICAL_TREATMENT
$---------------------------------------------------------—
OUTPUT_&_POST_PROCESS
ON_LAST_MESH
STEPS=10 $ write post-process file every ‘STEPS’ time steps
END_OUTPUT_&_POST_PROCESS
$---------------------------------------------------------—cylinder.nsi.dat
$---------------------------------------------------------—
PHYSICAL_PROBLEM
PROBLEM_DEFINITION
TEMPORAL_DERIVATIVES: On $ transient problem
CONVECTIVE_TERM: On $ off for Stokes flow (negligible for high viscosity fluids)
VISCOUS_TERM: LAPLACIAN $ suitable for constant viscosity fluid assumption
END_PROBLEM_DEFINITION
END_PHYSICAL_PROBLEM
$---------------------------------------------------------—
NUMERICAL_TREATMENT
ELEMENT_LENGTH: Minimum $ element length for Alya to calculate critical time step and stabilization
STABILIZATION: ASGS $ numerical stabilization method
TIME_INTEGRATION: Trapezoidal, ORDER: 2, EULER=20 $ time integration scheme
SAFETY_FACTOR: 100.0 $ multiply global time step: Alya calculates critical time
step as required by explicit solvers, which is suited for transient analysis but makes stationary solutions converge very slowly
STEADY_STATE_TOLER: 1e-12 $ onvergence tolerance for steady state
NORM_OF_CONVERGENCE: LAGGED_ALGEBRAIC_RESIDUAL
MAXIMUM_NUMBER_OF_IT: 1e5 $ max number of inner NS iterations
CONVERGENCE_TOLERANCE: 1e-3 $ convergence tolerance for inner NS loop
ALGORITHM: SCHUR $ NS solution algorithm (uncouples p & V)
END_ALGORITHM
MOMENTUM $ velocity solver
ALGEBRAIC_SOLVER
SOLVER: GMRES, KRYLOV=10 $ solver suited for asymmetric matrix
CONVERGENCE: ITERA=1000, TOLER=1.0e-10, ADAPTIVE, RATIO=0.01 $ max iter #, convergence criteria
ADAPTIVE, RATIO=0.001 means that the loop will end also if the difference of the norm of convergence value changes less than 0.1% in two consecutive iterations
OUTPUT: CONVERGENCE $ solver convergence file (.cso) is generated
PRECONDITIONER: DIAGONAL $ matrix preconditioner type
END_ALGEBRAIC_SOLVER
END_MOMENTUM
CONTINUITY $ pressure solver
ALGEBRAIC_SOLVER
SOLVER: DEFLATED_CG, COARSE: SPARSE $ CG solvers are suited for symmetric matrix
CONVERGENCE: ITERA=1000, TOLER=1.0e-10, ADAPTIVE, RATIO=0.01 $ max iter #, convergence criteria
ADAPTIVE, RATIO=0.001 means that the loop will end also if the difference of the norm of convergence value changes less than 0.1% in two consecutive iterations
OUTPUT: CONVERGENCE $ solver convergence file (.cso) is generated
PRECONDITIONER: DIAGONAL $ matrix preconditioner type
END_ALGEBRAIC_SOLVER
END_CONTINUITY
END_NUMERICAL_TREATMENT
$---------------------------------------------------------—
OUTPUT_&_POST_PROCESS
START_POSTPROCES_AT STEP = 0 $ initial step to post process
POSTPROCESS VELOCITY STEPS = 1 $ post process velocity every 1 step (ker.dat prevails)
POSTPROCESS PRESSURE STEPS = 1 $ post process pressure every 1 step (ker.dat prevails)
POSTPROCESS SCHUR STEPS = 1 $ post process schur every 1 step (ker.dat prevails)
BOUNDARY_SET
FORCE $ obtain pressure & velocity forces acting on boundaries
TORQUE, CENTER=0,0,0 $ obtain torque of forces acting on boundaries, related to CENTER
MASS $ obtain mass flow through boundaries
END_BOUNDARY_SET
END_OUTPUT_&_POST_PROCESS
$---------------------------------------------------------—
BOUNDARY_CONDITIONS
PARAMETERS
INITIAL: COARSE $ uses groups defined in dom.dat to generate an approximate steady initial solution
END_PARAMETERS
CODES, NODES
$ flow around the cylinder $ boundary conditions
1 11 1.000000 0.000000 $ Ux = 1 at inlet (left domain boundary)
2 00 0.000000 0.000000 $ free condition for Ux and Uy (p~0) at the left outlet
3 11 0.000000 0.000000 $ wall condition: Ux = Uy = 0 around the cylinder
4 01 0.000000 0.000000 $ symmetry condition: Ux = free , Uy = 0
1 & 4 11 1.000000 0.000000 $ inlet corners
2 & 4 01 0.000000 0.000000 $ outlet corners
END_CODES $ (note that BC’s apply to boundaries, not to boundary sets)
END_BOUNDARY_CONDITIONS
$---------------------------------------------------------—cylinder.post.alya.dat
$----------------------------------------------------------------—
DATA
FORMAT: visit $ also valid for ParaView
MARK_ELEMENTS: type $ to create automatic layers in post process according to some criterion
ELIMINATE_BOUNDARY_NODES: yes $ for parallel runs, to avoid node duplicity between subdomain
BOUNDARY: ON $ to post process boundary mesh
SUBDOMAINS, ALL $ subdomains (parallel partitions) to post process
END_SUBDOMAINS
END_DATA
$----------------------------------------------------------------—