Advanced Computational Models Fluent_notes

Advanced Computational Models Grid Adaptation Non-conformal Interfaces Moving Boundaries Deforming Boundaries…

ME469B/4/GI

1

Grid Adaptivity

The computational grid can be refined and/or coarsened based on geometrical and numerical solution data

Useful for: Capture flow features in details Increase resolution in near-wall regions Improve grid quality …

ME469B/4/GI

2

Grid Adaptivity Example Computational Domain Flow in a complex passage

A uniform triangular grid is likely to be inappropriate to capture all the feature of the flow ME469B/4/GI

3

Adaptation Process

Definition of the adaptation function (based on geometrical and/or solution data) Selection of the cells to be refined or coarsened (marking or tagging) Selection of grid refinement/coarsening scheme Adaptation Interpolation of the previous solution onto the new grid (automatic)

ME469B/4/GI

4

Adaptation Functions

Geometrical Region Boundary Volume

Solution based Isovalue Gradient y+

ME469B/4/GI

5

Region Adaptation This is the same option we used for “global” grid refinement to study the grid convergence of the solutions Adapt Æ Region Select a region shape

Input the geometrical definition of the region ME469B/4/GI

6

Boundary Adaptation Adapt Æ Boundary Select a boundary of the computational domain

Three options: cell distance normal distance volume distance

ME469B/4/GI

7

Boundary Adaptation Adaptation based on a cell’s distance from the selected boundary measured in number of cells. (1 means only the cells attached to the boundary) Adaptation based on a cell’s normal distance from the selected boundary

Adaptation based on a target boundary volume (specify a target volume and a growth factor a)

Vcell = Vtarget e ad This adaptation attempts to generate boundary-layer type grid

ME469B/4/GI

8

Volume Adaptation Adapt Æ Volume Based on cell volume (area in 2D):

Two options: magnitude (threshold values) change (neighbor change)

ME469B/4/GI

Allow to compute the range in your grid

9

Isovalue Adaptation Adapt Æ IsoValue All field values are available (including equation residuals, customized and gridrelated functions)

Method: specify the function specify a range specify the inside/outside option ME469B/4/GI

10

Gradient Adaptation As before all field values are available (including equation residuals, customized and grid-related functions)

Adapt Æ Gradient

Method: specify the function specify a range Note that the option is refine/coarsen Cells above the Refine Threshold are Refined Cells below the Coarsen Threshold are Coarsened ME469B/4/GI

11

y+ Adaptation Useful for turbulent flow simulations

Adapt Æ y+

Method: specify the wall boundaries specify the range (as before)

Note y* is a friendlier version of y+ y+ = r yp ut /m ME469B/4/GI

y* = r Cm1/4 kp1/2 yp /m 12

Mark (Tag) and display Before performing the adaptation you can mark the cell selected (based on a certain adaptation function) and display them

Adapt Æ Region Æ Control Æ Display

Only the selected cells will be displayed

ME469B/4/GI

13

Mark or Adapt? Mark allows to evaluate the effect of a refining/coarsening procedure without actually changing the grid

Marked cells are saved in registers that can be combined or manipulated before Adapting

adapt

mark

Adapt is usually used only at the end when all The desired cells are marked

ME469B/4/GI

14

Managing Registers Marked cells are saved in registers

Adapt Æ <method> Æ Manage

Adapt registers: collection of cells Mask registers: binary tagging to all cells in active and inactive

Operations on the registers are: Union: combine two adapt registers (or more) Intersection: combine adapt and mask registers Change Type: convert a adapt in mask and viceversa Delete: eliminate a register Exchange: swap an adapt register (coarse Æ refine) Invert: swap a mask register (active Æ inactive) Limit: apply the adaptations limits to a register Fill: mark for coarsening all the cells not in the register ME469B/4/GI

15

Adaptation controls Set limits on the adaptation procedure (i.e. maximum number of final cells)

Adapt Æ <method> Æ Manage Æ Control

Select the adaptation scheme Conformal Hanging Nodes

ME469B/4/GI

16

Hanging-Node and Conformal Adaptation

Cell to be refined

ME469B/4/GI

17

Hanging-Node Adaptation

father

Selected cells (father) are refined homothetically (kids) Nodes are added on the edges of the father cell and connectivity information are generated to link the kids to the father neighbors

kids

Memory penalty associated to the additional connectivity Information required and to the presence of an inactive father cell Neighboring cells are not allowed to differ more than one-level of refinement (because of inaccuracy related to large volume variations) Coarsening can be only performed on previously refined regions. Kids are deleted and the father cell becomes active ME469B/4/GI

18

Conformal Adaptation Selected cells are refined by splitting the longest edge It is inherently conservative because the cell connectivity is not modified No memory penalty (the old cells are deleted and the only the new are stored) Low quality meshes can be improved (refinement not homothetic) Coarsening can be applied everywhere and corresponds to a local remeshing Can only be applied to triangular (tetrahedral) grids In is NOT as local as the hanging-node approach

ME469B/4/GI

19

Grid Adaptivity Guidelines

Surface mesh should be fine enough to capture all the essential geometrical features of the model (especially for high curvature surfaces) The initial mesh should be fine enough to capture the overall features of the flow A reasonably converged solution must be obtained before adapting the grid Suitable flow-adaptation criteria are crucial to obtain increased resolution of selected region (i.e. velocity gradients are better than pressure gradient in incompressible flows and high values of turbulent quantities are relevant for turbulent flows)

ME469B/4/GI

20

Boundary Grid Adaptivity Poor boundary resolution cannot be improved via grid refinement

Original geometry ME469B/4/GI

21

Non-conformal grids Grid generation can be simplified in certain problems by meshing various components separately. The grids have to be coupled using a non-conformal interface

Quadrilateral grid

Triangular grid ME469B/4/GI

22

Non-conformal grid interfaces The matching conditions between the two faces have to be defined Coupled interface: the interface is actually an internal face (no bc) Periodic interface: the matching allows to specify pressure gradients Define Æ Grid Interface

ME469B/4/GI

23

Non-conformal grid interfaces - Example Channel flow Periodic Boundaries

Non-conformal interface ME469B/4/GI

Velocity contours (cell values) 24

Non-conformal Interface Guidelines

Grid interface can be of any shape (2D and 3D) but the surfaces to be mached MUST be based on the same geometry especially for highly curved surfaces Grid resolution at the two sides of the interfaces can be different; accuracy (and fluxes conservation) degrades for highly different mesh size Non-conformal interfaces are binary connectivity between two (and only two) zones.

ME469B/4/GI

25

Moving Zones Several industrial applications involve flow through a domain which contains a moving component (propellers, turbines, etc.)

stationary

moving

Time ME469B/4/GI

26

Advanced Computational Models

• Grid Adaptivity • Non-Conformal Grid Interfaces • Moving Zones

ME469B/4/GI

27

Modeling Moving Zones

Different approaches can be followed: 1) 2) 3) 4) 5) 6)

Single Reference Frame Model (SRFM) Multiple Reference Frame Model (MRFM) Mixing Plane Model (MPM) Sliding Mesh Model (SMM) Mesh Deformation Remeshing

ME469B/4/GI

28

Single Reference Frame Model Simplest model available; entire computational domain is referred to a rotating reference frame (domain moves with the reference frame) The equations are rewritten in the moving reference frame and a Coriolis acceleration term appears as a source for the momentum balance Boundaries that move with the frame can assume any shape BUT boundaries that are stationary MUST be surface of revolution

ME469B/4/GI

29

Single Reference Frame Set-Up Define Æ Models Æ Solver Velocity formulation Absolute: the unknowns are the field variables in the stationary reference frame Relative: the unknowns are the field variables in the moving reference frame Note: the relative is only available in the Segregated solver and it is usually faster

ME469B/4/GI

30

Single Reference Frame Boundary Conditions Define Æ Boundary Condition Æ Fluid

Define the axis of rotation Select moving reference frame Define the angular speed

ME469B/4/GI

Æ Wall

Stationary surfaces: zero absolute rotational speed Moving surfaces: zero rotational speed relative to the adjacent zone 31

Multiple Reference Frame Model

Many rotating machinery problems involve stationary components that cannot be represented as surface of revolution or move at different velocity The extension is to consider separate reference frames for each component The set-up is similar to the SRFM MRFM ignores the relative motions of subdomains and the Coriolis body-forces are local to each region (no equivalence between stationary and moving reference frame) Suitable for problems where the interaction between rotating and non-rotating components is small ME469B/4/GI

32

Sliding Mesh Model Like the MRFM the domain is divided into moving and stationary components Unlike the MRFM the mesh in each subdomain moves with respect to one another and the problem is inherently unsteady The equations are solved in the stationary reference frame and the meshes are moved at each time step (no approximations to the governing equations) The sliding mesh interface is defined as a non-conformal interface

ME469B/4/GI

33

SMM Example

time

ME469B/4/GI

Velocity distribution

34

Grid Generation

• Geometry definition (simple shapes, CAD import) • Grid generation algorithms • GAMBIT • Grid quality and improvement • Automation

ME469B/2/GI

1

Grid generation package: GAMBIT

Special Geometry Grid BC Tools

Geometry Tools

Graphics Window

Volume Tools

Visualization Tools

Text Window

The User Manual is available from the class Web Site ME469B/2/GI

2

1

Geometrical types - Topology

• Vertex • Edge (2 or more vertices) • Face (3 or more edges) • Volume (3 or more faces)

Bottom-up approach: generate low dimensional entities and build on top of them higher dimensional entities

ME469B/2/GI

3

Vertex: • Input Coordinates…

Edge: • Line segment (connect 2 vertices) • Circular arc • Quadratic functions • NURBS: Non-Uniform Rational B-Splines (connect N vertices)

Topologically any edge is ALWAYS a connection between 2 vertices (additional vertices used to build the geometry are NOT part of the edge) ME469B/2/GI

4

2

Face: • Rectangular • Circular •… • Sweep (translation or rotation of an edge) • Wireframe (connecting 3 or more edges)

ME469B/2/GI

5

Volume: • Cuboid • Sphere • Cone • Pyramids •… • Sweep (translation or rotation of a face) • Wireframe (connecting 3 or more faces)

ME469B/2/GI

6

3

Manipulate Geometry - Boolean Operations: Volume 1

• Unite • Substract • Intersect Volume 2

it generates the intersection edge

ME469B/2/GI

7

Manipulate Geometry - Blend

smooth sharp edges ME469B/2/GI

8

4

Create Entities - Faces

ME469B/2/GI

9

Create Entities - Faces

Some entities generated using “primitives” have fewer lower topological entities Example: Cube volume: 6 faces, 12 edges, 8 vertices Cylinder volume: 3 faces, 2 edges, 2 vertices ME469B/2/GI

10

5

Manipulate Geometry – Create Entities Create a vertex on a face Parametric representation of the face

ME469B/2/GI

11

Manipulate Geometry – Create Entities Create two edges by splitting an edge Parametric representation of the edge

ME469B/2/GI

12

6

Import Geometries • Realistic geometries are TOO complicated to be generated from “simple” shapes • Engineering design is based on CAD systems

Translation between CAD and CFD system is a major bottleneck • Gambit is based on ACIS geometrical libraries • ACIS (Andy, Charles & Ian’s System) is the most widely used 3D modeling software technology (http://www.spatial.com) • It can also import: • STEP (STandard for Exchange of Product model data; ISO standard) • IGES (Initial Graphics Exchange Specification; ANSI standard) • STL (STereo Lithography; Rapid Prototyping Standard) • …. ME469B/2/GI

13

Clean-Up a CAD Model

• Eliminate components not exposed to the flow • Eliminate duplicated entities • Eliminate small details • Water-proofing the surfaces • Rebuild geometrical connectivity between parts

ME469B/2/GI

14

7

Example: Helicopter Rotor

The rotor-shaft connection is VERY complicated ME469B/2/GI

15

Geometrical entities

The geometry consists of 10 components 3 blades, 1 shaft support, 6 connectors ME469B/2/GI

16

8

Example: Import IGES Model

IGES export is available from every CAD system IGES models are a collection of “untrimmed” edges and faces

Imported geometry consists of: 0 volumes ~250 faces ~1100 edges ~1000 vertices ME469B/2/GI

17

Example: Import STEP Model

STEP export is available from many CAD system STEP models are a collection of parts or components

Imported geometry consists of: 10 volumes ~190 faces ~450 edges ~300 vertices ME469B/2/GI

18

9

Components “exploded”

connector-1

support

blade

connector-2 For aerodynamic analysis the details of the rotor-shaft connection are not important. The geometry of the blades MUST be preserved ME469B/2/GI

19

Geometry simplification

Connectors eliminated

Support generated as a “simple” cylinder with blended side

ME469B/2/GI

20

10

Geometry simplification

Blade edge cleaned and sealed

Blade-support connector is a “simple” cuboid

ME469B/2/GI

21

Clean Geometry

This example is available on the class Web site ME469B/2/GI

22

11

Grid Generation

• Geometry definition (simple shapes, CAD import) • Grid generation algorithms • GAMBIT • Grid quality and improvement • Automation

ME469B/2/GI

23

Grid generation techniques

• Structured grids • Ordered set of (locally orthogonal) lines • Several Techniques can be used to Map a computational domain into a physical domain: Transfinite Interpolation, Morphing, PDE Based, etc. • The grid lines are curved to fit the shape of the boundaries • Unstructured grids • Unorganized collection of polygons (polyhedron) • Three main techniques are available to generate automatically triangles (tetrahedra): Delaunay triangulation, Advancing front, OCTREE • Paving for automatic generation of quads in 2D

ME469B/2/GI

24

12

Structured Grids: Mapping

Side B

Side A1

Physical Domain

ME469B/2/GI

Side B

k j

j=NJ=1 k=NK Side C

Side A1 Side A2

Side D

eC

Side A2

Sid

j=1=NJ

k=1 Side D

Computational Domain

25

Unstructured Grids: Triangulations • Delaunay • Empty circle principle: any node must not be contained within the circumcircle (circle passing through the vertices of a triangle) on any triangle within the mesh • Automatic triangulation of random set of nodes • Nodes are inserted locally in a triangulation and triangles are redefined locally to satisfy the Delaunay criterion (available mathematical tools) + Inherent grid quality + Elegant mathematical basis - Boundary integrity

ME469B/2/GI

26

13

Unstructured Grids: Triangulations • Advancing front • Triangles are built inward from the boundary surfaces • The last layer of elements constitutes the active front • An optimal location for a new nodes is generated for each segment on the front; the new node is generated by checking all existing nodes and this new optimal location • Intersection checks are required to avoid front overlap + Surface grid preserved + Specialized layers near surfaces - Computationally complex - Low quality

ME469B/2/GI

27

Unstructured Grids: Triangulations • OCTREE • Squares containing the boundaries are recursively subdivided until desired resolution is obtained • Irregular cells (or triangulation) are generated near the surface where square intersect the boundary + Requires least of surface representation + Highly automated - Cannot match surface grid - Low quality near surfaces

ME469B/2/GI

28

14

Unstructured Grids: Paving • Advancing front technique based on quads (instead of triangles) • Only in 2D

Triangulation ME469B/2/GI

Paving 29

Unstructured Grids: Coopering • 2D mesh sweeping • Only for cylindrical volumes • unstructured surface mesh is generated on surface A (source face) • structured grids are generated on cylindrical surfaces C & D • mesh on surface A is sweeped in the volume to generate the full 3D m esh

ME469B/2/GI

30

15

Unstructured Grids: 3D elements

ß Standard Elements:

Hex

ß ß ß ß ß

Tet

Pyramid

Wedge

Hex: Maximum Volume Covered per Edge Size Hex: Maximum Ratio Nodes/Elements Hex/Wedges: Clustering at Solid Wall with High Quality Elements Tets: Automatic Meshing of Extremely Complicated Regions Pyramids/Wedges: Transition Between Tets & Hex

ME469B/2/GI

31

Unstructured Grids: Hex or Tets? We NEED Hex-Based Meshing because: ß Equiangular Tets are NOT Good for Thin Volumes ß Too Many Elements for Reasonable Resolutions (estimated >2M grid Points in conical-annular Swirler)

35K Total Elements

ME469B/2/GI

Cross-Section of the Swirler

410K Total Elements

32

16

Speed

Robustness

Quality & Control

Complex Geometry

Mesh sizes

• Structured gridding (mapping)

+

+

+

-

+

• Unstructured triangulation (2D/3D)

-

+

+/-

+

-

-

-

+/-

+

+

+

-

+/-

-

+

What is available in GAMBIT

• Unstructured paving (2D) • Unstructured coopering (3D)

All GAMBIT meshes are exported as unstructured collection of (mixed) elements ME469B/2/GI

33

Grid generation – 1D - Edges

Straightforward Select number of points Select distribution of points Edge direction is defined from 1st to 2nd vertex Clustering toward one side is defined accordingly

ME469B/2/GI

34

17

Grid generation – 2D - Faces

Easy Select number of points • use predefined edge meshes • use uniform spacing Select meshing scheme • constraints on the edge meshing for mapping and paving schemes

ME469B/2/GI

35

Grid generation – 2D - Faces

It is possible to force the cell element type at face-vertices elements forced to be triangles

Mixing element-type is one of the main advantages of unstructured mesh technology

ME469B/2/GI

36

18

Grid generation – 2D - Mesh-patching options

Matching interface

Overlapping interface ME469B/2/GI

Non-conformal interface

Mixed-element interface 37

Grid generation – 3D - Volumes

Not so easy Select number of points • use predefined face meshes • use uniform spacing Select meshing scheme • constraints on the face meshing for mapping and cooper schemes

ME469B/2/GI

38

19

3D Grid generation – Advanced Cooper technique

“Creative” way of coopering: multisurface to multisurface sweep ME469B/2/GI

39

Grid generation – Sizing functions Instead of the bottom-up approach (1D to 3D) grid generation Sizing functions can be specified to mesh volumes directly

ME469B/2/GI

40

20

Grid generation – Clustering points Sizing functions can be used effectively to define the size of the cells BUT they cannot provide directional control (anisotropy) One option is to build (grow) elements from the boundaries and to form “viscous” layers

ME469B/2/GI

41

Grid generation – Boundary Layers

ME469B/2/GI

42

21

Example – meshing a circle Mapping

Boundary Layer Paving

Triangulation

Boundary Layer Paving

Paving

Multiblock mapping

Boundary Layer Multiblock Paving

Boundary Layer Transition Multiblock Paving

Circle defined as segments ME469B/2/GI

43

Mesh linking

Edges, faces and volume meshes can be linked Define corresponding entities and ALSO reference entities Needed to “enforce” coincident grids on different entities (i.e. for periodicity bc)

ME469B/2/GI

44

22

Virtual Geometry GAMBIT operates on two different type of entities REAL: with corresponding geometrical and topological characteristics VIRTUAL: defined only with reference to REAL entities REAL entities are what we used and described so far; VIRTUAL are used to SIMPLIFY, CLEAN UP, DECOMPOSE real entities Note that some geometry tools cannot be applied to virtual entities (boolean operation, volume blending, creation of volumes by sweeping faces, etc.)

Real entities can be transformed in virtual but NOT viceversa ME469B/2/GI

45

Virtual Geometry Clean-Up Example of edge connecting operation

REAL ME469B/2/GI

VIRTUAL 46

23

Grid quality Quality measures are NOT absolute but should be considered in connection with solution schemes The final accuracy of a procedure is ALWAYS a function of the grid quality

Several geometrical measures can be defined: • • •

Depending on the size of the elements Depending on the shape of the elements Depending on relative dimensions of neighboring elements

ME469B/2/GI

47

Quality measures available in GAMBIT

ME469B/2/GI

48

24

Examine meshes • Define mesh element to examine • Define a cutting plane • Define the quality measure

Aspect ratio

Equiangular skew

ME469B/2/GI

49

Mesh improvement • Smoothing operators are applied to redistribute nodes

Aspect ratio

Equiangular skew

Typically, improvement in 3D meshes are based on improved 2D meshes ME469B/2/GI

50

25

Automation

• GAMBIT saves a “journal” file with the commands issues during a session • Journal file are ASCII editable files • Commands are quasi-English and easy-to-use • They are useful to trace-back sessions to find errors • They can be made general by introducing User Defined Parameters

ME469B/2/GI

51

GAMBIT Journal file

The command: Volume create width 1 depth 1 height 1 offset 0 0 0 brick

Generates a cube of size 1 centered at 0 0 0 On the other hand the sequence $W = 2.3 $D = 1.5 $H = 4 Volume create width $W depth $D height $H offset 0 0 0 brick

Generates a cuboid of User-Specified Size centered at 0 0 0

ME469B/2/GI

52

26

Example of Journal file

ME469B/2/GI

/ -------------------------------------------------------/ CYCLONE GRID GENERATION / ME269B - Spring 2002 / -------------------------------------------------------/ / R1 = External radius of Cyclone / R2 = Gas Outlet Pipe (External) / R3 = Gas Outlet Pipe (Internal) / RB = Particles Outlet (Bottom) / H1 = Height of the Cylindrical Part of Cyclone / H2 = Height of the Conical Part / HE = Depth of the Outlet Channel into the Cyclone / / inletl = Length (x) of the Gas Inlet Channel / inleta = Height (z) of the Gas Inlet Channel / inletb = Span (y) of the Gas Inlet Channel / outletl = Length (z) of the outlet (gas) pipe / / cellsize = Average size of the cells / / Remark: z-axis is the Cyclone Axis / -------------------------------------------------------/ / Input Quantities / $R1 = 1.555 $R2 = 0.45 $R3 = 0.4 $RB = 0.75 $H1 = 4.5 $H2 = 5 $HE = 3.6 $inletl = 2 $inleta = 1.03 $inletb = 0.7 $outletl = 7.2 $cellsize = 0.18 / /

53

Grid generation research ß Structured grids: automatic generation of mappable subdomains NLR 2000-366 Report (PDF available from class web site)

ß Unstructured tetrahedral grids: anisotropic Delaunay schemes Shimada et al. “High quality anisotropic tetrahedral mesh generation via ellipsoidal bubble packing” (PDF available)

ME469B/2/GI

54

27

Grid generation research ß Unstructured hex-dominant grids: OCTREE based SAMM – Computational Dynamics Ltd. Hexpress – Numeca International Inc.

ß Unstructured purely hexahedral grids: Whisker-Weaving CUBIT – Sandia National Lab. (PS report available)

ME469B/2/GI

55

Grid generation – Links and References

• Links • Mesh generation and grid generation on the Web http://www-users.informatik.rwth-aachen.de/~roberts/meshgeneration.html • Meshing research corner http://www.andrew.cmu.edu/user/sowen/mesh.html • General CFD: Topic mesh generation http://www.cfd-online.com

• References: • Handbook of grid generation. Thompson, Soni, Weatherill, CRC Press • Numerical Grid Generation: Foundation & Applications. Thompson, Warsi, Mastin. North Holland Press ME469B/2/GI

56

28

Example #1 Geometry Modeling Objective: generation of a swirl-generator model in GAMBIT • •

Four rectangular ducts inject tangentially in an expanding pipe. Select the dimensions appropriately

The duct offset DO is the main design variable: with fixed mass flow the swirl generated is a function of DO

DO

Use the “cyclone.jou” example as a baseline ME469B/2/GI

The shape of the pipe/duct intersection is a function of 57 the duct offset

Example #2 Grid Generation

Objective: generation of a grid for a swirl-generator model in GAMBIT • Generate an all-tetrahedra grid • Generate an all-hexahedra grid • Compare the quality of the grid especially in the channel-pipe intersection

ME469B/2/GI

58

29

Solution methods for the Incompressible Navier-Stokes Equations

l l l l l l

Discretization schemes for the Navier-Stokes equations Pressure-based approach Density-based approach Convergence acceleration Periodic Flows Unsteady Flows

ME469B/3/GI

1

Background (from ME469A or similar)

Navier-Stokes (NS) equations Finite Volume (FV) discretization Discretization of space derivatives (upwind, central, QUICK, etc.) Pressure-velocity coupling issue Pressure correction schemes (SIMPLE, SIMPLEC, PISO) Multigrid methods

ME469B/3/GI

2

1

NS equations Conservation laws: Rate of change

+

advection

+ diffusion = source

=0

The advection term is non-linear The mass and momentum equations are coupled (via the velocity) The pressure appears only as a source term in the momentum equation No evolution equation for the pressure There are four equations and five unknowns (r, V, p) ME469B/3/GI

3

NS equations

Compressible flows: The mass conservation is a transport equation for density. With an additional energy equation p can be specified from a thermodynamic relation (ideal gas law) Incompressible flows: Density variation are not linked to the pressure. The mass conservation is a constraint on the velocity field; this equation (combined with the momentum) can be used to derive an equation for the pressure

ME469B/3/GI

4

2

Pressure-based solution of the NS equation The continuity equation is combined with the momentum and the divergence-free constraint becomes an elliptic equation for the pressure

To clarify the difficulties related to the treatment of the pressure, we will define EXPLICIT and IMPLICIT schemes to solve the NS equations:

It is assumed that space derivatives in the NS are already discretized:

ME469B/3/GI

5

Explicit scheme for NS equations

Semi-discrete form of the NS Explicit time integration The n+1 velocity field is NOT divergence free Take the divergence of the momentum Elliptic equation for the pressure

ME469B/3/GI

6

3

Explicit pressure-based scheme for NS equations

Velocity field (divergence free) available at time n

Compute Hn

Solve the Poisson equation for the pressure pn

Compute the new velocity field un+1

ME469B/3/GI

7

Implicit scheme for NS equations

Semi-discrete form of the NS Implicit time integration Take the divergence of the momentum

The equations are coupled and non-linear

ME469B/3/GI

8

4

Implicit scheme for NS equations Compute an intermediate velocity field (eqns are STILL non-linear)

{

Define a velocity and a pressure correction

Using the definition and combining

{

Derive an equation for u’

Enforce divergence-free condition at n+1

ME469B/3/GI

9

Implicit pressure-based scheme for NS equations (SIMPLE) SIMPLE: Semi-Implicit Method for Pressure-Linked Equations

Velocity field (divergence free) available at time n Compute intermediate velocities u* Solve the Poisson equation for the pressure correction p’ Neglecting the u*’ term

Solve the velocity correction equation for u’ Neglecting the u*’ term

Compute the new velocity un+1 and pressure pn+1 fields

ME469B/3/GI

10

5

Implicit pressure-based scheme for NS equations (SIMPLEC) SIMPLE: SIMPLE Corrected/Consistent

Velocity field (divergence free) available at time n Compute intermediate velocities u* Solve the Poisson equation for the pressure correction p’ Use an approximation to u*’ (neighbor values average u*’ ~ S u’)

Solve the velocity correction equation for u’ Use an approximation to u*’

Compute the new velocity un+1 and pressure pn+1 fields

ME469B/3/GI

11

Implicit pressure-based scheme for NS equations (PISO) PISO: Pressure Implicit with Splitting Operators

Velocity field (divergence free) available at time n Compute intermediate velocities u* and p’ as in SIMPLE

Solve the Poisson equation for the pressure correction p(m+1)’ u*’ is obtained from u m’

Solve the velocity correction equation for u(m+1)’ u*’ is obtained from u m’

Compute the new velocity un+1 and pressure pn+1 fields ME469B/3/GI

12

6

SIMPLE, SIMPLEC & PISO - Comments In SIMPLE under-relaxation is required due to the neglect of u*’ un+1 = u* + au u’

p = pn + ap p’

There is an optimal relationship

ap =1- au

SIMPLEC and PISO do not need under-relaxation SIMPLEC/PISO allow faster convergence than SIMPLE PISO is useful for irregular cells

ME469B/3/GI

13

NS equations Conservation laws: Rate of change

+

advection

+ diffusion = source

=0

ME469B/3/GI

14

7

Segregated solver in FLUENT FV discretization for mixed elements

f W

The quantities at the cell faces can be computed using several different schemes ME469B/3/GI

15

Discretization of the equations Options for the segregated solver in FLUENT Discretization scheme for convective terms 1st order upwind (UD) 2nd order upwind (TVD) 3rd order upwind (QUICK), only for quad and hex Pressure interpolation scheme (pressure at the cell-faces) linear (linear between cell neighbors) second-order (similar to the TVD scheme for momentum) PRESTO (mimicking the staggered-variable arrangement) Pressure-velocity coupling SIMPLE SIMPLEC PISO ME469B/3/GI

16

8

Solution of the equation f is one of the velocity component and the convective terms must be linearized:

This correspond to a sparse linear system for each velocity component

Fluent segregated solver uses: Point Gauss-Seidel technique Multigrid acceleration

b b

ME469B/3/GI

P

b

17

Set-up of problems with FLUENT Read/Import the grid Define the flow solver option Define the fluid properties Define the discretization scheme Define the boundary condition Define initial conditions Define convergence monitors Run the simulation Analyze the results

Graphics Window

Command Menus

Text Window ME469B/3/GI

18

9

Solver set-up Define Æ Models Æ Solver

Define Æ Controls Æ Solution

Example: text commands can be used (useful for batch execution) define/models/solver segregated define/models/steady

solve/set/discretization-scheme/mom 1 solve/set/under-relaxation/mom 0.7

…

… ME469B/3/GI

19

Material properties Define Æ Materials

Quantities are ALWAYS dimensional ME469B/3/GI

20

10

Initial and boundary conditions Solve Æ Initialize Æ Initialize

Only constant values can be specified More flexibility is allowed via patching

Define Æ Boundary Conditions

BCs will be discussed case-by-case

ME469B/3/GI

21

Initial conditions using patching Adapt Æ Region Æ Mark

Mark a certain region of the domain (cells are stored in a register)

ME469B/3/GI

Solve Æ Initialize Æ Patch

Patch desired values for each variable in the region (register) selected

22

11

Convergence monitors Solve Æ Monitors Æ Residuals

Solve Æ Monitors Æ Surface

Convergence history of the equation residuals are stored together with the solution User-defined monitors are NOT stored by default ME469B/3/GI

23

Postprocessing Display Æ Contours

Plot Æ XY Plot

Cell-centered data are Computed This switch interpolates the results on the cell-vertices

ME469B/3/GI

24

12

Detailed post-processing Define additional quantities Define plotting lines, planes and surfaces Compute integral/averaged quantities Define Æ Custom Field Function

ME469B/3/GI

25

Fluent GUI - Summary File: I/O Grid: Modify (translate/scale/etc.), Check Define: Models (solver type/multiphase/etc.), Material (fluid properties), Boundary conditions Solve: Discretization, Initial Condition, Convergence Monitors Adapt: Grid adaptation, Patch marking Surface: Create zones (postprocessing/monitors) Display: Postprocessing (View/Countors/Streamlines) Plot: XY Plots, Residuals Report: Summary, Integral Typical simulation Parallel: Load Balancing, Monitors

ME469B/3/GI

26

13

Example – Driven cavity Classical test-case for incompressible flow solvers

Problem set-up

Solver Set-Up

Material Properties: r = 1kg/m3 m = 0.001kg/ms

Segregated Solver Discretization: 2nd order upwind SIMPLE

Reynolds number: H = 1m, Vslip= 1m/s Re = rVslipH/m = 1,000

Vslip=1

Boundary Conditions: Slip wall (u = Vslip) on top No-slip walls the others

H

Multigrid V-Cycle

Initial Conditions: u=v=p=0 Convergence Monitors: Averaged pressure and friction on the no-slip walls

ME469B/3/GI

27

Example – Driven cavity

The effect of the meshing scheme

Quad-Mapping 1600 cells

Tri-Paving 3600 cells

Quad-Paving 1650 cells

Edge size on the boundaries is the same ME469B/3/GI

28

14

Example – Driven cavity

The effect of the meshing scheme – Vorticity Contours

Quad-Mapping 1600 cells

Tri-Paving 3600 cells

Quad-Paving 1650 cells

ME469B/3/GI

29

Example – Driven cavity The effect of the meshing scheme – Convergence

Quad-Mapping 1600 cells ME469B/3/GI

Tri-Paving 3600 cells

Quad-Paving 1650 cells 30

15

Example – Driven cavity The effect of the meshing scheme x-velocity component in the middle of the cavity

Quad-Mapping Tri-Paving

Quad-Paving

ME469B/3/GI

Symbols corresponds to Ghia et al., 1982 31

Example – Driven cavity Grid Sensitivity – Quad Mapping Scheme

Vorticity Contours

1600 cells ME469B/3/GI

6400 cells

25600 cells 32

16

Example – Driven cavity Grid Sensitivity – Quad Mapping Scheme x-velocity component in the middle of the cavity

1600 cells

6400 cells

25600 cells

ME469B/3/GI

Symbols corresponds to Ghia et al., 1982 33

How to verify the accuracy?

Define a reference solution (analytical or computed on a very fine grid) Compute the solution on successively refined grids Define the error as the deviation of the current solution from the reference Compute error norms Plot norms vs. grid size (the slope of the curve gives the order of accuracy) Problems with unstructured grids: 1) Generation of a suitable succession of grids 2) Definition of the grid size

ME469B/3/GI

34

17

Generation of successively refined grid 1) Modify grid dimensions in GAMBIT and regenerate the grid 2) Split all the cells in FLUENT Adapt Æ Region Æ Adapt

The region MUST contain the entire domain

Element shape & metric properties are preserved

ME469B/3/GI

35

Driven Cavity - Error evaluation Reference solution computed on a 320x320 grid (~100,000 cells) Reference solution interpolated on coarse mesh to evaluate local errors

Quad-Mapping

Tri-Paving

Quad-Paving

Note that the triangular grid has more than twice as many grid cells ME469B/3/GI

36

18

Driven Cavity – Accuracy evaluation

Nominal 1st order accuracy

Nominal 2nd order accuracy

Tri meshing scheme yields Slightly higher errors and lower accuracy Note that the definition of Dx is questionable (a change will only translate the curves not change the slope)

Quad-Mapping Tri-Paving

Error (L2 norm)

Quad and Pave meshing schemes yield very similar accuracy (close to 2nd order)

(N)-1/2

Quad-Paving

ME469B/3/GI

37

Driven Cavity – Fluent vs. other CFD codes Quad Mapping Scheme (1600 cells) x-velocity component in the middle of the cavity

FLUENT ME469B/3/GI

StarCD

NASA INS2D

Symbols corresponds to Ghia et al., 1982 38

19

Techniques for the incompressible NS equations

Pressure correction schemes Fractional step methods Artificial compressibility approach Vorticity-streamfunction formulation Density-based approach

ME469B/3/GI

39

Techniques for the incompressible NS equations Fractional-step approach Similar in concept to pressure-correction schemes, but the pressure (or a related quantity) is only used in the final correction step; it is equivalent to an approximate factorization

Vorticity-streamfunction approach It is effectively a change-of-variables; introducing the streamfunction and the vorticity vector the continuity is automatically satisfied and the pressure disappears (if needed the solution of a Poisson-like equation is still required). It is advantageous in 2D because it requires the solution of only two PDEs but the treatment of BCs is difficult. In addition in 3D the PDEs to be solved are six

Artificial compressibility approach A time-derivative (of pressure) is added to the continuity equation with the goal of transforming the incompressible NS into a hyperbolic system and then to apply schemes suitable for compressible flows. The key is the presence of a user-parameter b (related to the artificial speed of sound) that determines the speed of convergence to steady state ME469B/3/GI

40

20

Density-based solvers for the NS equations The equation are written in compressible form and, for low Mach numbers, the flow is effectively incompressible

The energy equation is added to link pressure and density through the equation of state

In compact (vector) form:

ME469B/3/GI

41

Density-based solvers for the NS equations Stiffness occurs because of the disparity between fluid velocity and speed of sound (infinite in zero-Mach limit) The equations are solved in terms of the primitive variables

Note that the continuity becomes (again) an evolution equation for the pressure

The time derivative is modified (preconditioned) to force all the eigenvalues to be of the same order (similar to the artificial compressibility approach)

Several choices are available for G ME469B/3/GI

42

21

FLUENT density-based solver

Explicit Scheme Multistage Runge-Kutta scheme Residual Smoothing Multigrid acceleration Implicit Scheme Euler (one-step) implicit Newton-type linearization Point Gauss-Seidel iterations Multigrid acceleration

ME469B/3/GI

43

Example – Driven cavity Classical test-case for incompressible flow solvers

Problem set-up Material Properties: r = 1kg/m3 m = 0.001kg/ms Reynolds number: H = 1m, Vslip= 1m/s Re = rVslipH/m = 1,000

Vslip=1

Solver Set-Up Coupled Solver Discretization: 2nd order upwind Implicit Multigrid V-Cycle

Boundary Conditions: Slip wall (u = Vslip) on top No-slip walls the others

H

Initial Conditions: u=v=p=0 Convergence Monitors: Averaged pressure and friction on the no-slip walls

ME469B/3/GI

44

22

Example – Driven cavity Effect of the solver - Quad mesh (1600 cells) Segregated

Coupled

Vorticity Contours ME469B/3/GI

45

Example – Driven cavity Effect of the solver - Quad mesh (1600 cells) x-velocity component in the middle of the cavity

Segregated ME469B/3/GI

Coupled

Symbols corresponds to Ghia et al., 1982 46

23

Multigrid acceleration

Basic idea: the global error (low-frequency) on a fine grid appears as a local error (high-frequency) on coarse meshes. Why it is important: linear system solver like Gauss-Seidel are effective in removing high-frequency errors but VERY slow for global errors. Note that, on structured, grid line-relaxation (or ADI-type) schemes can be used to improve the performance of Gauss-Seidel; on unstructured grid similar concepts are extremely difficult to implement. Convergence Speed: number of iterations on the finest grid required to reach a given level of convergence is roughly independent on the number of grid nodes (multigrid convergence)

ME469B/3/GI

47

Two-grid scheme

1. 2. 3. 4. 5.

a smoothings are performed on the fine grid to reduce the highfrequency components of the errors (pre-smoothing, aS) the residual (error) is transferred to next coarser level (restriction, R) g iterations are performed on this grid level for the “correction” equation the problem is transferred back to the fine grid (prolongation, P) b smoothings are performed on the fine grid to remove the highfrequency errors introduced on the coarse mesh (post-smoothing, bS)

Parameters to be defined are a, b, g

ME469B/3/GI

48

24

Restriction & Prolongation Operators

Fine Level Coarse Level ME469B/3/GI

49

Algebraic Multigrid

The coarse levels are generated without the use of any discretization on coarse levels; in fact no hierarchy of meshes is needed Geometric multigrid should perform better than AMG because non-linearity of the problem are retained on coarse levels (correction equation) AMG is effectively a solver for linear systems and the restriction and prolongation operators might be viewed as means to modify (group or split) the coefficient matrix The same multigrid options (cycle type and a, b, g smoothing) are available

ME469B/3/GI

50

25

Driven Cavity – Accuracy evaluation

Nominal 1st order accuracy

Nominal 2nd order accuracy

Tri meshing scheme yields Slightly higher errors and lower accuracy Note that the definition of Dx is questionable (a change will only translate the curves not change the slope)

Quad-Mapping Tri-Paving

Error (L2 norm)

Quad and Pave meshing schemes yield very similar accuracy (close to 2nd order)

Quad-Paving

ME469B/3/GI

(N)-1/2 51

Multigrid for unstructured meshes GEOMETRIC MULTIGRID ALGEBRAIC MULTIGRID

ME469B/3/GI

52

26

Algebraic Multigrid Performance Convergence for the segregated solver

1600 cells

6400 cells

25600 cells

ME469B/3/GI

53

Algebraic Multigrid Performance Convergence for the coupled solver

1600 cells ME469B/3/GI

6400 cells

25600 cells 54

27

Periodic Flows Geometrical periodicity

Periodicity simply corresponds to matching conditions on the two boundaries

The velocity field is periodic BUT the pressure field is not. The pressure gradient drives the flow and is periodic. A pressure JUMP condition on the boundary must be specified

ME469B/3/GI

55

Periodic Flows – Set-Up In the segregated solver periodicity can be imposed by fixing either the mass flow or the pressure drop In the coupled solver periodicity is enforced by fixing the pressure drop Define Æ Periodic Conditions

Segregated solver ME469B/3/GI

Define Æ Boundary Conditions

Coupled Solver 56

28

Periodic Flow Example – 2D channel An analytical solution of the NavierStokes equations (Poiseuille flow) can be derived:

h

y

Solution in the form u=u(y) The pressure drop balances the viscous drag on the walls Navier-Stokes equations Velocity distribution in the channel Averaged velocity ME469B/3/GI

57

Periodic Flow Example – 2D channel Problem set-up

Solver Set-Up

Material Properties: r = 1kg/m3 m = 0.1kg/ms

Coupled Solver

Reynolds number: h = 2m, Vave= 1m/s Re = rVsliph/m = 20

h

Periodic boundaries

Discretization: 2nd order upwind SIMPLE Multigrid V-Cycle

Boundary Conditions: Periodicity Dp=0.3 No-slip top/bottom walls Initial Conditions: u = 1; v = p = 0 Exact solution: Vave = 1

ME469B/3/GI

58

29

Periodic Flow Example – 2D channel x-velocity distribution in the channel

Cell-centered values showed (no interpolation)

Quad-Mapping

Tri-Paving

ME469B/3/GI

59

Periodic Flow Example – 2D channel The error in this case CAN be computed with reference to the exact solution In this case the computed averaged velocity error is plotted

Nominal 1st order accuracy Nominal 2nd order accuracy

(N)-1/2 Quad-Mapping Tri-Paving ME469B/3/GI

This test-case is available on the class web site

60

30

Unsteady flows The algorithms we used so far are time-marching: From an initial condition they iterate until a steady-state is reached The time-evolution of the solution is NOT accurate Typical Implicit Time-Accurate Scheme

Unsteady transport equation 1st order time integration 2nd order time integration

ME469B/3/GI

61

Unsteady flows – Dual time stepping At each time step we build a steady-state-like problem and we use an iterative algorithm

Outer iteration (time-step)

ME469B/3/GI

Inner iteration steady-state pseudo time-step

62

31

Unsteady flows – Set Up Define Æ Models Æ Solver

Solve Æ Iterate Outer iteration

Inner iteration

ME469B/3/GI

63

Unsteady Flow – Impulsive start-up of a plate Again an analytical solution of the NavierStokes equations can be derived:

y

Vwall

Solution in the form u=u(y,t) The only force acting is the viscous drag on the wall Navier-Stokes equations Velocity distribution Wall shear stress ME469B/3/GI

64

32

Unsteady Flow – Impulsive start-up of a plate Problem set-up periodic boundaries

Solver Set-Up

Material Properties: r = 1kg/m3 m = 0.1kg/ms

Segregated Solver Discretization: 2nd order upwind SIMPLE

Reynolds number: Re = r VwallL/m Vwall = 5.605 L = mVwall/twall

H

Multigrid V-Cycle

Boundary Conditions: Slip wall (u = Vwall) on bottom No-slip wall (top) Periodicity Dp=0 Initial Conditions: u=v=p=0

Vwall

Exact Solution twall = 1 @ t = 1

ME469B/3/GI

H/L ~ 10

65

Unsteady Flow – Impulsive start-up of a plate

Time discretization 1st order

2nd order

Influence of the BCs ME469B/3/GI

Nominal 1st order accuracy

Nominal 2nd order accuracy

Error (L2 norm)

The error CAN be computed with reference to the exact solution In this case the computed wall shear stress is plotted

Dt

This test-case is available on the class web site

66

33

Overview of commercial CFD codes About 30 packages. Three major general-purpose products (covering ~50% of the market): FLUENT, StarCD, CFX

Grid Type

Pressure Based

Density Based

Multigrid

System Solver

Discretization

FLUENT

Unstructured Mixed

SIMPLE SIMPLEC PISO

Coupled Implicit Preconditioned

Algebraic Geometric

Gauss-Seidel

UD/TVD QUICK

StarCD

Unstructured Mixed

SIMPLE SIMPISO PISO

-

-

Conjugate Gradient

UD/TVD QUICK CD

Unstructured Mixed

SIMPLE -

Algebraic Coupled

ILU

UD/TVD QUICK CD

CFX

ME469B/3/GI

67

34

Advanced Physical Models

• Heat Transfer • Buoyancy • Combustion and reaction modeling • Multiphase flows • Solidification and melting

ME469B/5/GI

1

Heat Transfer Thermal analysis are crucial in many industrial applications

Turbulence is enhanced is internal turbine cooling passages to improve heat transfer Roughness elements (ribs) are placed in the channels

ME469B/5/GI

2

Heat Transfer Modeling An energy equation must be solved together with the momentum and the continuity equations For incompressible flows the energy equation is decoupled from the others (r is NOT a function of the temperature) For laminar flows the energy equation can be solved directly; for turbulent flows after Reynolds-averaging the equation contains an unclosed correlation: Momentum

Energy

ME469B/5/GI

3

Heat Transfer Modeling

Momentum

Energy

The Prandtl number is the measure of the momentum diffusivity vs. the thermal diffusivity Pr=cp m/k Pr is order 1 for gases (typically 0.7 for air) Prt is an additional parameter in the turbulence model (typically 0.9) ME469B/5/GI

4

Set-Up for Heat Transfer Calculations Define Æ Models Æ Energy Activate the energy equation

Define Æ Materials

Specify material properties

ME469B/5/GI

5

Wall thermal boundary conditions The options are: 1) Fixed temperature 2) Fixed thermal flux (temperature gradient)

wall thickness computational domain

ME469B/5/GI

BC

6

Flow-thermal simulations For the energy equations all the numerical options (discretization, underrelaxation, etc.) are available

For incompressible fluids the temperature and momentum equations are decoupled

ME469B/5/GI

7

Periodic flows Many heat-transfer devices are characterized by geometrically periodic configurations (ribbed passages) Temperature behaves like the pressure: it varies in the streamwise direction but its variation (gradient) is periodic The energy equations can be rewritten in terms of a scaled temperature:

And the modified energy equation can be solved with periodic BC

ME469B/5/GI

8

Example of Heat Transfer Calculations Ribbed Passages

Wall heat transfer (temperature gradients) are strongly connected to wall friction coefficients and therefore to turbulence modeling ME469B/5/GI

9

Example: Ribbed Channel Flow h h

H

L

Periodic boundaries

Problem set-up

Solver Set-Up

Material Properties: r = 1kg/m3 m = 0.0001kg/ms Cp = 1000 J/Kg/oK k = 0.142 W/m oK

Segregated Solver

Reynolds number: h = 1m, L=10m, H=L Reh = rUbh/m = 10,000

Multigrid V-Cycle

Discretization: 2nd order upwind SIMPLE

Boundary Conditions: Periodicity m=rUbH=10Kg/s No-slip walls Initial Conditions: u = 1; v = p = 0

ME469B/5/GI

Turbulence model: k-e

10

Grids in Ribbed Channel Flows

Grid points are clustered at the walls and in the shear layers

Unstructured gridding allows to separate the bottom and top BLs having different resolutions in the streamwise direction ME469B/5/GI

11

Heat Transfer Predictions

x/e ME469B/5/GI

12

Conjugate Heat Transfer In many cases the correct prediction of the thermal field in a device requires the inclusion of conduction effects in solids Conjugate simulations are referred to coupled fluid-solid temperature calculations

ME469B/5/GI

13

Set-Up Conjugate Heat Transfer We need to specify two zones (fluid and solid) in the grid generation And then specify the material properties Define Æ Materials Æ Fluid

ME469B/5/GI

Define Æ Materials Æ Solid

14

Wall thermal boundary conditions The boundary between the two zones is ALWAYS a wall and a shadow zone is created automatically by Fluent

wall thickness

shadow wall thickness Solid zone

Fluid zone

coupling condition

ME469B/5/GI

15

Conjugate Heat Transfer The temperature field is affected by the treatment of the rib walls

The rib is insulated

The rib is heated from the base

The rib is uniformly heated

ME469B/5/GI

16

Effect of the Thermal Wall Bc

ME469B/5/GI

17

Effect of the Thermal Wall Bc

ME469B/5/GI

18

Comments on Thermo-fluid simulations

In incompressible flows the energy equation is decoupled from the momentum equations and can be solved a posteriori with the velocity field frozen Additional modeling is involved for the solution of thermal equation in the RANS context (therefore additional approximations and errors) Wall quantities (temperature and heat flux) are very sensitive to the modeling of near-wall turbulence Conjugate heat transfer (coupled fluid/solid) are often necessary to describe accurately a physical device

ME469B/5/GI

19

User Programming

• What are User Defined Functions • Introduction to C • Set-Up C User Routines in Fluent • Implementing a Custom Turbulence Model • Programming in other CFD Commercial Codes

Acknowledgement: this handout is partially based on Fluent training material ME469B/6/GI

1

Introduction to UDF Programming Why programming in commercial codes? • The codes are general-purpose but cannot anticipate all needs • New (physical) models can be developed in a user-friendly environment • Large number of problems (test-cases) can be addressed with the same implementation

What is a the User Defined Function (UDF)? • C (Fluent) or FORTRAN (StarCD, CFX) routines programmed by the user linked to the solver to perform certain operations: • initialization • special boundary condition (i.e. space or time dependent) • material properties • source terms • reaction rates • postprocessing and reporting • debugging • …. ME469B/6/GI

2

1

A Simple UDF Example Specify a Parabolic Velocity Profile at the Inlet Goal: The UDF (inlet_parab) set the values of the x-velocity component at the cell faces of the inlet boundary zone

Inlet zone

1) 2) 3) 4)

Determine the cell-faces belonging to the inlet zone Loop on all those faces Determine the coordinate of the face centroid Specify the x-velocity component ME469B/6/GI

3

A Simple UDF Example Parabolic Velocity Profile:

function inlet_parab Definitions Loop over all the inlet cell faces Evaluate the face centroid coordinates The function return the velocity

ME469B/6/GI

4

2

A Simple UDF Example Compile/interpret the UDF: Define Æ User Defined

Attach the profile to the inlet zone Define Æ Boundary Condition

Equivalent to attach the profile from a separate simulation ME469B/6/GI

5

A Simple UDF Example Solve the equations Solve Æ Iterate

Final Result

ME469B/6/GI

6

3

C Programming Typical C function:

ME469B/6/GI

7

C vs. FORTRAN

ME469B/6/GI

8

4

Basic Programming Rules Statements MUST end with a semicolon

Æ

;

Comments are placed anywhere in the program between Æ /* …… */ Statements are grouped by curly brackets

Æ { …… }

Variables defined within the body functions are local Variables defined outside the body functions are global and can be used by all the functions that follows Variables MUST be ALWAYS defined explicitly Integer/real/double functions have to return a integer/real/double value C is case sensitive! ME469B/6/GI

9

Basic Statements Arithmetic expressions in C look like FORTRAN a = y + (i-b)/4;

Functions which return values can be used in assignment a = mycompute(y,i,b);

Functions which do not return values can be called directly mystuff(a,y,i,b);

Mathematical and various other default functions are available a = pow(b,2);

ME469B/6/GI

/* pow(x,y) returns x raised to y */

10

5

Data Types and Arrays Arrays of variables are defined using the notation int a[10]; /*define an array of ten integers */ Note the brackets are square [] not round ()!! Note that the arrays ALWAYS start from 0

Standard C types are: Integer, Real, Double, Char C allows to define additional types typedef struct list{int id, float x[3], int id_next};

ME469B/6/GI

11

Pointers A pointer is a variable which contain the address of another variable Pointers are defined as: int *a_p; int a;

/* pointer to an integer variable */ /* an integer variable */

Set-up a pointer a = 1; a_p = &a;

/* &a return the address of variable a */

*a_p returns the content of the address pointed by a_p

ME469B/6/GI

12

6

Operators Arithmetic Operators

Logical Operators

ME469B/6/GI

13

Conditional Statements if and if-else

ME469B/6/GI

Example

14

7

Loop Procedure for loops

Example

ME469B/6/GI

15

C Preprocessor It handles Macro substitutions #define A B #define my_f(x,y) x+y*3-5

The preprocessor replaces A with B It also handles file inclusion #include “udf.h” #include “math.h”

The files to be included MUST reside in the currect directory

ME469B/6/GI

16

8

Programming in C

Of course much more than just this….

Additional information are: http://www.cs.cf.ac.uk/Dave/C/CE.html Plenty of books: “The C Programming Language”, Kernighan & Ritchie, 1988

ME469B/6/GI

17

UDF in Commercial CFD Codes Commercial CFD codes allow the development of User Defined Functions for various applications. Anyway, the core of the software is closed. UDF must be compiled and linked to the main code Most codes provide macros and additional tools to facilitate the use of UDFs In Fluent there are two options: Interpreted The code is executed on a “line-by-line” basis at run time + does not need a separate compiler (completely automatic) - slows down the execution and uses more memory - somewhat limited in scope Compiled A library of UDF is compiled and linked to the main code Overcomes all the disadvantages reported above ME469B/6/GI

18

9

Interpret the UDFs Define Æ User Defined Æ Interpreted

Display of code translation in assembly (and eventual compiling errors)

Default stack size might be too small for large arrays! ME469B/6/GI

19

Compile the UDFs

Define Æ User Defined Æ Compiled

The library MUST be precompiled in a directory tree ME469B/6/GI

20

10

Directory tree for compiled UDFs Machine dependent irix6r10 ultra hp700 …

my_library

src

Makefile makefile

lnx86

my_udf.c

3d

2d makefile

my_udf.c

libudf.so makefile

my_udf.c

ME469B/6/GI

libudf.so 21

Makefiles for UDFs In the directory /usr/local/Fluent.Inc/fluent6/src

There are two files Makefile.udf makefile.udf

to be copied in the directory my_library to be copied in the directory my_library/src

The first one does not require modifications. In the second one two macros MUST be modified SOURCE = my_udf.c FLUENT_INC = /usr/local/Fluent.Inc

ME469B/6/GI

22

11

UDFs in FLUENT Boundary Conditions Initial Conditions Adjust Function Source Terms Material Properties Execute on Demand User Defined Scalars User Defined Memory

Available MACROS

Pointers to threads Geometry Macros Cell and Face Variables Arithmetic and Trigonometric Functions

Programming Tools

ME469B/6/GI

23

Boundary Conditions The inlet profile specification is based on the macro DEFINE_PROFILE DEFINE_PROFILE can be used for wall boundary conditions as well to impose temperature, velocity, shear stress, etc. It is possible to specify a constant value, a position-dependent or time-dependent values and to link the values on the boundary to the values in the internal cells Internal Values

x

x

x

BC

Note that the BCs are applied to the faces of the cells (face centroids) ME469B/6/GI

24

12

Initial Conditions The initial condition specification can be performed using the macro DEFINE_INIT Space-dependent conditions might be imposed The routine is executed once at the beginning of the solution process It is attached in the UDFs hooks Define Æ User Defined Æ Function Hooks

ME469B/6/GI

25

Initial Conditions Note that the difference between the DEFINE_PROFILE and DEFINE_INIT is that the former performs a loop on the face-centroids belonging to a certain zone whereas the latter loops on the cell-centroids Example:

ME469B/6/GI

26

13

Adjust Function Similar to the DEFINE_INIT but it is executed every iteration: DEFINE_ADJUST Can be used for postprocessing, clipping, etc. It is attached in the UDFs hooks Define Æ User Defined Æ Function Hooks

ME469B/6/GI

27

Source Terms

To add a source term in the equations: DEFINE_SOURCE Can be used for: Continuity Momentum (component by component) Turbulent quantities Energy It is different from the previous macros because it works on a cell-by-cell basis (no need to perform loops!)

ME469B/6/GI

28

14

Source Terms Define Æ Boundary Conditions Æ Fluid

Activate source terms

Attach the UDF

ME469B/6/GI

29

Source Terms

In FLUENT the source terms are written in the form S (f) = A + B f where f is the dependent variable A is treated explicitly and B f is treated implicitly In the DEFINE_SOURCE both terms A and B have to be specified

ME469B/6/GI

30

15

Material Properties

To define the properties of the materials : DEFINE_PROPERTIES Can be used for: Viscosity Density Thermal Conductivity …. As for the source terms works on a cell-by-cell basis

ME469B/6/GI

31

Material Properties Define Æ Material Property All the available UDFs Are shown in the menu

ME469B/6/GI

32

16

Execute on Demand To perform operations instantaneously : EXECUTE_ON_DEMAND It is executed when activated by the user Can be used for: Debugging Postprocessing …. Define Æ User Defined Æ Execute on Demand

ME469B/6/GI

33

User Defined Scalar In addition to the Continuity and Momentum Equations (and Energy) generic transport equations can be solved (up to 50!) To include scalars in the calculations 1) Define the number of scalars Define Æ User Defined Æ User Defined Scalars 2) Define the diffusion and source terms in the generic equation

ME469B/6/GI

34

17

User Defined Scalar When UDS are defined in the material property panel the diffusivity of the scalars MUST be defined Define Æ Material Property

Note that the default diffusivity for the scalars is constant and equal to 1! ME469B/6/GI

35

User Defined Scalar Boundary conditions for the scalars can be defined as Constant Value Constant Gradient User Defined Define Æ Boundary Conditions Æ Wall

ME469B/6/GI

36

18

User Defined Scalar The Source terms for the scalars are set using the DEFINE_SOURCE macro Introduced before The scalar equations can be further customized 1) Convective terms can be modified using the macro DEFINE_UDS_FLUX 2) Unsteady term can be modified using the macro DEFINE_UDS_UNSTEADY

ME469B/6/GI

37

User Defined Memory The User Defined Scalars are used to solve additional equations eventually coupled to the mean flow equations Sometimes it is useful to define temporary field variables to store and retrieve values not directly available UDM are defined:

Define Æ User Defined Æ User Defined Memory

More efficient storage compared to User Defined Scalars

ME469B/6/GI

38

19

I/O in UDF User Defined Scalars and User Defined Memory are AUTOMATICALLY Stored in the Fluent Data files Additional I/O (from and to the Fluent Case and Data files) can be Accomplished using the macro DEFINE_RW_FILE

Define Æ User Defined Æ Function Hooks

ME469B/6/GI

39

Detailed Programming The macros introduced before are interpreted/compiled and attached to the various Fluent panels The detailed programming is based on additional macros that allow to loop on cells to retrieve field variables, etc. Loop Macros Geometry Macros Field Variable Macros Control Macros Arithmetic and Trigonometric Macros Before looking at the Macros, the Fluent Data structure in introduced ME469B/6/GI

40

20

Data Structure It is based on a hierarchy of structures Threads (zones) are collection of cells or faces Domain is a collections of threads

Domain Thread

Thread

Thread

Cell

Face

Face

ME469B/6/GI

41

Threads Every zone is associated to a single ID (available in the boundary conditions menu)

zone ID Given the ID of a thread the pointer can be retrieved as: Thread *tf = Lookup_Thread(domain, ID);

Given the thread pointer the ID of the zone can be retrieved as: ID = THREAD_ID(thread); ME469B/6/GI

42

21

Loop Macros

ME469B/6/GI

43

Cell-face connectivity When looping on faces the surrounding cells can be accessed using the macros: cell0_thread = F_C0_THREAD(face,thread) cell1_thread = F_C1_THREAD(face,thread) cell0_id = F_C0(face,thread) cell1_id = F_C1(face,thread)

x

cell0

cell1

face

When looping on cells adjacent faces and cells can be accessed using the macros: for(nf=0;nf
Number of faces of each cell is unknown 44

22

Geometry Macros

Many more available; refer to the FLUENT UDF Manual ME469B/6/GI

45

Field Variables Macros

ME469B/6/GI

46

23

Field Variables Macros

Many more available; refer to the FLUENT UDF Manual ME469B/6/GI

47

Control Macros

Many more available; refer to the FLUENT UDF Manual ME469B/6/GI

48

24

An example of User Defined Programming Low-Re k-e model Development of a custom turbulence model can be accomplished using the UDFs. k-e model with damping functions formulation:

Transport Equations

ME469B/6/GI

49

An example of User Defined Programming Eddy Viscosity P= S=

Low-Re k-e model

Turbulent Kinetic Energy Production Strain Rate Magnitude

Damping functions

y

Turbulent Reynolds number definitions Wall boundary conditions

ME469B/6/GI

50

25

Required UDF Routines

Source Terms

DEFINE_SOURCE(k_source, thread, eqn) DEFINE_SOURCE(d_source, thread, eqn)

Diffusivity

DEFINE_PROPERTY(ke_diffusivity, domain)

Boundary Condition

DEFINE_PROFILE(wall_d_bc, domain)

Adjust Routine

DEFINE_ADJUST(turb_adjust, domain)

ME469B/6/GI

51

Required field variables

Density

C_R(cell,thread)

Molecular viscosity

C_MU(cell,thread)

Eddy viscosity

C_MU_T(cell,thread)

Strain Rate Magnitude

Strain_rate(cell,thread)

Wall distance

C_WALL_DIST(cell,thread)

ME469B/6/GI

52

26

UDF Header #include "udf.h" /* Turbulence model constants */ #define C_MU 0.09 #define SIG_TKE 1.0 #define SIG_TDR 1.3 #define C1_D 1.44 #define C2_D 1.92 /* User-defined scalars */ enum { TKE, TDR, N_REQUIRED_UDS };

Required: includes all Fluent macros

Constant definitions (global)

Assign a number to each scalar

ME469B/6/GI

53

Damping Functions These are defined on a cell-by-cell basis /* Reynolds number definitions */ real Re_y(cell_t c, Thread *t) { return C_R(c,t)*sqrt(C_UDSI(c,t,TKE))*C_WALL_DIST(c,t)/C_MU_L(c,t);} real Re_t(cell_t c, Thread *t) { return C_R(c,t)*SQR(C_UDSI(c,t,TKE))/C_MU_L(c,t)/C_UDSI(c,t,TDR);} /* Damping Functions */ real f_mu(cell_t c, Thread *t) { return tanh(0.008*Re_y(c,t))*(1.+4/pow(Re_t(c,t),0.75));} real f_1(cell_t c, Thread *t) { return 1.;} real f_2(cell_t c, Thread *t) { return (1.-2/9*exp(-Re_t(c,t)*Re_t(c,t)/36))*(1.-exp(-Re_y(c,t)/12));}

ME469B/6/GI

54

27

Source Term Routines The production term in the k equation is: e is obtained from the definition of the eddy viscosity to increase the coupling between the equations and to define an implicit term

DEFINE_SOURCE(k_source, c, t, dS, eqn) { real G_k; G_k = C_MU_T(c,t)*SQR(Strainrate_Mag(c,t)); dS[eqn] = -2.*C_R(c,t)*C_R(c,t)*C_MU*f_mu(c,t)*C_UDSI(c,t,TKE)/C_MU_T(c,t); return G_k - C_R(c,t)*C_R(c,t)*C_MU*f_mu(c,t)*SQR(C_UDSI(c,t,TKE))/C_MU_T(c,t); }

ME469B/6/GI

55

Source Term Routines The production term in the e equation is: It contains already both k and e. No need for manipulations! DEFINE_SOURCE(d_source, c, t, dS, eqn) { real G_k; G_k = C_MU_T(c,t)*SQR(Strainrate_Mag(c,t)); dS[eqn] = C1_D*f_1(c,t)*G_k/C_UDSI(c,t,TKE) - 2.*C2_D*f_2(c,t)*C_R(c,t)*C_UDSI(c,t,TDR)/C_UDSI(c,t,TKE); return C1_D*f_1(c,t)*G_k*C_UDSI(c,t,TDR)/C_UDSI(c,t,TKE) - C2_D*f_2(c,t)*C_R(c,t)*SQR(C_UDSI(c,t,TDR))/C_UDSI(c,t,TKE); }

ME469B/6/GI

56

28

Diffusivity The diffusion terms in the scalar equations are set-up together DEFINE_DIFFUSIVITY(ke_diffusivity, c, t, eqn) { real diff; real mu = C_MU_L(c,t); real mu_t = C_R(c,t)*C_MU*f_mu(c,t)*SQR(C_UDSI(c,t,TKE))/C_UDSI(c,t,TDR); switch (eqn) { case TKE: diff = mu_t/SIG_TKE + mu; break; case TDR: diff = mu_t/SIG_TDR + mu; break; default: diff = mu_t + mu; } return diff;

But each equation can have a different value

}

ME469B/6/GI

57

Eddy Viscosity The eddy viscosity is set in the adjust routine (called at the beginning of each Iteration) and it is used in the mean flow and in the scalar equations

DEFINE_ADJUST(turb_adjust, domain) { Thread *t; cell_t c; /* Set the turbulent viscosity */ thread_loop_c (t, domain) if (FLUID_THREAD_P(t)) { begin_c_loop(c,t) { C_MU_T(c,t) = C_R(c,t)*C_MU*f_mu(c,t)*SQR(C_UDSI(c,t,TKE))/C_UDSI(c,t,TDR); } end_c_loop(c,t) } }

ME469B/6/GI

58

29

Wall Boundary Conditions Only the boundary condition for e is complicated because it requires the value of the derivative of k (the square root of k) DEFINE_PROFILE(wall_d_bc, t, position) { face_t f; cell_t c0; Thread *t0 = t->t0; /* t0 is cell thread */ real xw[ND_ND], xc[ND_ND], dx[ND_ND], rootk, dy, drootkdy;

The derivative is rootk/dy begin_f_loop(f,t) { rootk is the sqrt of k in c0 = F_C0(f,t); the adjacent cell center rootk = sqrt(MAX(C_UDSI(c0,t0,TKE), 0.)); dy is the distance between F_CENTROID(xw,f,t); cell and face center C_CENTROID(xc,c0,t0); NV_VV(dx, =, xc, -, xw); dy = ND_MAG(dx[0], dx[1], dx[2]); drootkdy = rootk/dy; F_PROFILE(f,t,position) = 2.*C_MU_L(c0,t0)/C_R(c0,t0)*drootkdy*drootkdy; } end_f_loop(f,t) } ME469B/6/GI

59

Set-Up the Problem

Set-Up a case using the standard k-e Define two scalars (TKE, TDR) Compile and Attach the UDFs Deactivate the Equations for k and e Initialize and solve the RANS + TKE and TDR

ME469B/6/GI

60

30

Attach the UDF

Source terms

Boundary conditions

ME469B/6/GI

61

Open UDF Library

After compiling the library It can be opened in Fluent

Output in the text window Are the functions available

ME469B/6/GI

62

31

Perform the Simulation The convergence is for the mean flow + the two scalars (turbulent quantities deactivated!)

For postprocessing purposes the UDS have to be used instead of turbulent quantities ME469B/6/GI

63

Additional Information

Many additional macros are available to implement different physical models combustion models particles-based models …. It is formally possible to develop additional numerical methods flux discretizations variable reconstruction and clipping …. Information are available in the FLUENT UDF manual (class web site)

ME469B/6/GI

64

32

UDF in other commercial CFD Codes

CFX

CFX (v4) is a structured code; the data structure is much simpler because the field variables are accessed directly. It uses 1D arrays where the quantities are stored in succession CFX UDFs are written in FORTRAN For example: U(I,J,K) Æ U(IJK) where IJK = (K-1)*NI*NJ+(J-1)*NI There are no macros, but examples of subroutines to perform the customization USRBC USRSRC USRDIF …. ME469B/6/GI

UDF in other commercial CFD Codes

65

Star-CD

StarCD is an unstructured code similar to Fluent in term of data organization But similar to CFX for the organization of the UDF. StarCD has threads (as Fluent) and the UDF work normally on a cell-by-cell basis StarCD UDFs are written in FORTRAN There are no macros, but examples of subroutines to perform the customization BCDEFW SORSCA DIFFUS ….

ME469B/6/GI

66

33

Simulation of Turbulent Flows

• From the Navier-Stokes to the RANS equations • Turbulence modeling • k-e model(s) • Near-wall turbulence modeling • Examples and guidelines

ME469B/3/GI

1

Navier-Stokes equations The Navier-Stokes equations (for an incompressible fluid) in an adimensional form contain one parameter: the Reynolds number:

Re = r Vref Lref / m it measures the relative importance of convection and diffusion mechanisms

What happens when we increase the Reynolds number? ME469B/3/GI

2

1

Reynolds Number Effect 350K < Re

Turbulent Separation Chaotic

200 < Re < 350K

Laminar Separation/Turbulent Wake Periodic

40 < Re < 200

Laminar Separated Periodic

5 < Re < 40

Laminar Separated Steady

Re < 5

ME469B/3/GI

Laminar Attached Steady Experimental Observations

Re 3

Laminar vs. Turbulent Flow

Laminar Flow

The flow is dominated by the object shape and dimension (large scale) Easy to compute ME469B/3/GI

Turbulent Flow

The flow is dominated by the object shape and dimension (large scale) and by the motion and evolution of small eddies (small scales) Challenging to compute 4

2

Why turbulent flows are challenging? Unsteady aperiodic motion Fluid properties exhibit quasi-random spatial variations (3D) Strong dependence from initial conditions Contain a wide range of scales (eddies)

The implication is that the turbulent simulations MUST be always three-dimensional, time accurate with extremely fine grids

ME469B/3/GI

5

Direct Numerical Simulation The objective is to solve the time-dependent NS equations resolving ALL the scale (eddies) for a sufficient time interval so that the fluid properties reach a statistical equilibrium

Grid requirement: N ~ (Ret)9/4 ~ 1x107 for Ret = 800 Time step requirement: Dt ~ (Ret)-1/2 ~ 1x10-5 for Ret = 800 ME469B/3/GI

y+ = r yp ut /m

ut = (tw/r)1/2

6

3

Typical DNS channel flow simulations

ME469B/3/GI

Instantaneous Flow Field – Vortical Structures

7

DNS of Channel Flow: Averaged Results

ME469B/3/GI

Laminar and Turbulent Velocity profiles

8

4

DNS using Commercial CFD Codes?

Channel flow @ Ret = 395 ME469B/3/GI

9

Beyond DNS DNS is possible but only for low Reynolds number flows (and simple geometries) The (time and space) details provided by DNS are NOT required for design purposes Time averaged quantities are appropriate for engineering purposes Large scale resolution (not to the level of the smallest eddies) is enough for applications

Can we extract time-average and large-scale quantities at a reasonable computational cost?

ME469B/3/GI

10

5

Flow over a Backstep

Istantaneous high velocity regions

low velocity regions

Time averaged

Length of the recirculation region is of engineering interest ME469B/3/GI

11

Velocity & Velocity Fluctuations is a function of time with rapid fluctuations (turbulence) By re-defining it as the time average and the fluctuations

(x1)

(x2)

t

Velocity measurements in the turbulent wake of a cylinder ME469B/3/GI

12

6

Reynolds-Averaged Navier-Stokes Equations

Define Reynolds-averaged quantities

Substitute and average:

Closure problem ME469B/3/GI

13

Unsteady RANS? Every turbulent flow is unsteady BUT not all the unsteadiness is turbulence! RANS averaging based on time average can be applied only to “statistically” steady flows. What if flow has a large scale periodicity (vortex shedding)?

We can define the Reynolds-Averaging procedure in terms of Ensemble Average:

ME469B/3/GI

14

7

Turbulence modeling Define the Reynolds stresses in terms on known (averaged) quantities 1) Boussinesq hypothesis - simple relationship between Reynolds stresses and velocity gradients through the eddy viscosity (similar to molecular viscosity) - isotropic (eddy viscosity is a scalar!) 2) Reynolds stress transport models - equations derived directly manipulating the NS equations - still contain unknown (undetermined) quantities - no assumption of isotropy - very complicated and expensive to solve 3) Non-linear Eddy viscosity models (Algebraic Reynolds stress) 4) Model directly the divergence of the Reynods Stresses ME469B/3/GI

15

Eddy viscosity models Boussinesq relationship:

with:

Guidelines for defining the eddy viscosity: 1) Dimensional arguments - units are [m2/s] - define 2 out of three scales: velocity, length, time

2) Physical arguments - asymptotic analysis - consistency with experimental findings

3) Numerical arguments - simple and easy to compute ME469B/3/GI

16

8

Classification of eddy viscosity models

The various models (about 200) are classified in terms of number of transport equations solved in addition to the RANS equations: 1) 2) 3) 4) 5)

zero-equation/algebraic models: Mixing Length, Cebeci-Smith, Baldwin-Lomax, etc one-equation models: Wolfstein, Baldwin-Barth, Spalart-Allmaras, k-model, etc two-equation models: k-e, k-w, k-t, k-L, etc. three-equation models: k-e-A four-equation models: v2-f model

ME469B/3/GI

Zero-equation model

17

Prandtl Mixing Length

From dimensional arguments and analogy with molecular transport Definition of L is different for each problem (boundary layes, mixing layers, etc.) Eddy viscosity is zero if the velocity gradients are zero No “history” effect; purely local L can be made “universal” using ad hoc functions of distance from the walls, pressure gradients, etc. ME469B/3/GI

18

9

One-equation model

k-model

An equation from k can be derived directly from the NS equations (using the definition)

k1/2 is assumed to be the velocity scale it still requires a length scale L as before to define the eddy viscosity 4 out of 7 terms in the k equation require further assumptions Production is computed using the Boussinesq approximation Dissipation is modeled (using dimensionality arguments) as k3/2/L Turbulent transport and pressure diffusion are modeled together:

ME469B/3/GI

One-equation model

19

k-model

The final form of the model is:

The ONLY advantage with respect to zero-equation models is the inclusion of the history effects. Modern one-equation models abandoned the k-equation and are based on a ad-hoc Transport equation for the eddy viscosity directly. Spalart-Allmaras model:

ME469B/3/GI

20

10

Two-equation model

k-f family

The main drawback of the k one-equation model is the incomplete representation of the two scales required to build the eddy viscosity; two-equation models attempt to represent both scales independently. • All models use the transport equation for the turbulent kinetic energy k • Several transport variables are used e: turbulence dissipation rate L: turbulent length scale w: inverse of turbulent time scale w2 g y t

ME469B/3/GI

21

k-e model The k equation is the same as before:

The e equation can be obtained from the NS equations but it contains several undetermined quantities; it is therefore derived “mimicking” the k equation

The eddy viscosity is obtained as: There are 5 free constants

ME469B/3/GI

22

11

Determining the constants? The constants can be determined studying simple flows: 1.

Decaying homogeneous isotropic turbulence

2.

Homogeneous shear flow

3.

The Logarithmic Layer

4.

…

Or by comparison with experimental data Standard k-e refers to a certain choice of the constants (Launder & Sharma 1972)

ME469B/3/GI

23

Structure of the Turbulent Boundary Layer Universal Law (velocity profile)

At High Reynolds number the viscous dominated layer is so thin that it is very difficult to resolve it ME469B/3/GI

24

12

Wall Function Approach (High-Re k-e) The laminar sublayer is NOT resolved First grid point is assumed to be in the logarithmic layer (y+>11) and the velocity is assumed to be described by: u+ = (1/k)ln(Ey+) A slip condition (u ≠ 0) is imposed at the wall (imposed shear stress) k boundary condition is usually imposed as a zero-gradient. e is obtained by equilibrium condition (Pk=e ) If first grid point is too close (viscous layer) then the velocity is: u+ = y+ ME469B/3/GI

25

Near Wall Region Modeling

From a physical point of view: It is important because solid walls are the main source of vorticity and turbulence (local extrema of turbulent kinetic energy, large variations of turbulence dissipation, etc.) In engineering applications: Wall quantities (velocity gradients, pressure, etc.) are very important in several applications Flow separation and reattachment are strongly dependent on a correct prediction of the development of turbulence near walls ME469B/3/GI

26

13

Damping Functions Approach (Low-Re k-e) The equations are integrated to the wall WITHOUT assuming an universal law for the velocity profile and an equilibrium conditions for k and e The problem is that the model predicts the wrong behavior for k and e near the solid walls (from DNS and experimental observation) The equations are modified using algebraic functions to “damp” certain terms:

The damping functions are designed to correct the behavior of the eddy viscosity ME469B/3/GI

27

Damping Functions We need to specify 5 constants plus 3 functions:

MANY choices are available (about 30 different formulations!). Fluent has 6 different versions Classic model is the Launder and Sharma model:

Others might be function of the distance from the wall, the pressure gradient, etc. ME469B/3/GI

28

14

Two-Layer Approach The computational domain is divided in two regions: viscosity affected (near wall) and fully turbulent core (outer region) Two different models are used: the complete k-e model for the outer region and a simplified model (typically a one-equation k-based model) for the near-wall The separation between the two regions is defined in terms of a distance from the wall (y+~30) The main assumption is related to the definition of e which is based on

ME469B/3/GI

29

Near-Wall Treatments for k-e models

Approach

Physics

Grid requirement

Wall functions

-

+

+

+/-

Two layer

+/-

+/-

+

+/-

Damping functions

+/-

-

-

+/-

Numerics

Accuracy

Summary and Comparison

ME469B/3/GI

30

15

Example: Turbulent Channel Flow

h

L

Periodic boundaries

Problem set-up

Solver Set-Up

Material Properties: r = 1kg/m3 m = 0.0017kg/ms

Segregated Solver

Reynolds number: h = 2m, L=1m Ret = rUth/m = 590 Boundary Conditions: Periodicity Dp/L=1=Ut No-slip top/bottom walls Initial Conditions: u = 1; v = p = 0

Discretization: 2nd order upwind SIMPLE Multigrid V-Cycle

Grid SAME GRID USED FOR THE LAMINAR FLOW @ Re=20

Turbulence model: k-e with wall functions ME469B/3/GI

31

Turbulent Channel Flow Ret = 590

Velocity and k profiles ME469B/3/GI

32

16

Grid Convergence?

Turbulence Channel Flow

Ret = 590

Velocity and k profiles ME469B/3/GI

33

From High-Re to Low-Re k-e Wall clustering y+=30 y+~1

Wall boundary condition dk/dy=0 k=0

Grid ME469B/3/GI

Turbulent kinetic energy

Ret = 590

34

17

Low-Re k-e model

Turbulence Channel Flow

Ret = 590

Velocity and k profiles ME469B/3/GI

35

Pros & Cons of k-e model + Simple + Affordable + Reasonably accurate for wide variety of flows (without separation) + History effects - Overly diffusive - Cannot predict different flows with the same set of constants (universality) - Source terms are stiff numerically - Not accurate in the region close to no-slip walls where k and e exhibit large peaks (DNS and experimental observations) - Near wall treatment A lot of variants have been introduced to overcome these problems One (or more) constants become coefficients varying with S, distance from the walls, pressure gradient, etc.: RNG k-e, realizable k-e… ME469B/3/GI

36

18

Alternatives/Improvements to k-e models The k-w model was developed from the realization that most of the problems experienced by k-e-type models are due to the modeling of the e equation which is neither accurate or easy to solve (e has a local extrema close to the wall) Mathematically this is equivalent to a change of variables w~e/k The v2-f model is based on the argument that k/e is the correct turbulent time scale in the flow (close to the wall and in the outer region) but k is not the appropriate turbulent velocity scale An additional equation for the correct velocity scale v2 (independent from k) has to be solved. Moreover, the damping effect produced from the presence of the wall is NOT local (as assumed in the damping function approach) but must be accounted for globally using an elliptic equation ME469B/3/GI

37

k-w model

Turbulence Channel Flow

Ret = 590

Velocity and k profiles ME469B/3/GI

38

19

Reynolds Stress Models Attempt to model directly the “new” terms appearing in the RANS equations Mathematically is expensive because we have 6 additional equations:

More importantly ONLY the production term are exact, everything else MUST be modeled

RSM are extremely stiff numerically due to the coupling between the equations and suffer of the same near-wall problems of the k-e ME469B/3/GI

39

Models available in Fluent Define Æ Models Æ Viscous

One-equation model Two-equation model k-e (3 versions + 3 wall treatments) k-w (2 versions)

Reynolds Stress model

Note that the coefficients might be adjusted!!! ME469B/3/GI

40

20

Models in Fluent hidden behind the GUI Some turbulence models are NOT directly available in the GUI!!!

Mixing Length define/model/viscous/mixing-length y

Two-equation model k-e Low-Re (6 versions)

define/model/viscous/turbulence-expert/low-re-ke y define/model/viscous/turbulence-expert/low-re-ke-index 2 ME469B/3/GI

41

Comparison of Models

Turbulence Channel Flow

Ret = 590

Velocity and k profiles ME469B/3/GI

42

21

Channel Flow Summary

Wall functions are accurate if the first grid point is in the logarithmic layer Grid Convergence Study with wall functions approach FAILS Damping Functions (and Two-Layer approaches) are accurate for the velocity profiles But the turbulent kinetic energy peak is underpredicted k-w model is a viable alternative to k-e and has less sensitivity to the grid clustering SA model and v2-f model are equivalent in capturing the velocity profile v2-f model is accurate in predicting the peak of turbulence kinetic energy near the wall

ME469B/3/GI

43

Inlet boundary conditions for turbulent quantities At inlet boundary conditions additional quantities have to be specified in turbulent flows depending on the turbulence model selected. Typically there are three options: 1) k-e 2)

Turbulence Intensity and Turbulence Length Scale

3)

Turbulence Intensity and Turbulent Viscosity

ME469B/3/GI

44

22

Example: Asymmetric Diffuser Problem set-up

Solver Set-Up

Material Properties: r = 1kg/m3 m = 0.0001kg/ms

Initial Conditions: u = 1; v = p = 0

Reynolds number: h = 2m, Re = rUh/m = 20,000

Turbulence model: Two-equation models

Boundary Conditions: Inlet profiles available from experiments No-slip top/bottom walls

Segregated Solver Discretization: 2nd order upwind SIMPLE Multigrid V-Cycle

Outflow Inflow

ME469B/3/GI

45

Flow in Asymmetric Diffuser k-e Low-Re

Experiments indicated the presence of a large recirculation region v2-f

k-e models with damping function do NOT capture it

k-e Low-Re

v2-f ME469B/3/GI

46

23

Flow in Asymmetric Diffuser

ME469B/3/GI

Skin friction on the bottom wall

47

Flow in Asymmetric Diffuser Streamwise Velocity

ME469B/3/GI

Turbulence Kinetic Energy

48

24

k-e with wall functions

Series of grids generated with different clustering at wall

ME469B/3/GI

No separation is captured!

49

k-e with wall functions In “complex” configuration it is impossible to generated a grid with the first grid point in the logarithmic layer…..

…in addition, for complicated flows with recirculation the Universal Law is inaccurate ME469B/3/GI

50

25

Example: NLR Two Component Airfoil Problem set-up Material Properties: r = 1kg/m3 m = 3.98E-7Kg/ms Reynolds number: h = 1m, Re = rUh/m = 2,512,600

Solver Set-Up Initial Conditions: u = cosa; v = sina, p = 0

Turbulence model: SA and k-e models

Boundary Conditions: Constant velocity at AOA=a No-slip walls

Segregated Solver Discretization: 2nd order upwind SIMPLE Multigrid V-Cycle

a ME469B/3/GI

51

Computational Grid

Due to the high Reynolds number resolution of the boundary layers requires extreme clustering

ME469B/3/GI

52

26

How this grid is generated?

Multiblock Structured Grid ME469B/3/GI

53

Pressure distributions

Pressure Distribution at Low and High Angle of Attack ME469B/3/GI

54

27

Why k-e fails? Streamwise Velocity

Turbulent kinetic energy

Spurious production of k in the stagnation regions Fix: Use of a production limiter: define/model/viscous/turbulence-expert/kato-launder-model y

ME469B/3/GI

55

Guidelines Simulations of turbulence flows require “decisions” based on:

physics accuracy resources Turbulence model

Computational grid

CFD Process

1) Flow Physics to characterize the flow features (turbulence, high gradients, etc.) 2) Computational requirements to evaluate the grids resolution required for a certain accuracy 3) Project Requirement to evaluate the need for sophisticated turbulence models

Simulation ME469B/3/GI

56

28

Guidelines

Modeling procedure: • Determine relevant Reynolds number to estimate if the flow is turbulent • Select a turbulence model option and a near-wall treatment • Estimate the physical dimension of the first grid point off the wall (y+) • Generate the grid • Perform the simulation • “Reality” check (experiments, literature, model consistency)

ME469B/3/GI

57

Estimating y+ Definition: y+ = r yp ut /m y+ = r yp ut /m ut = Ue (cf/2)1/2 To estimate ut we can use “classic” relationships for the skin friction: • •

Flat Plate: Pipe:

cf/2 ~ 0.052 (Rex)-0.142 cf/2 ~ 0.046 (Rex)-0.2

Note that there are different laws for different Re number ranges

ME469B/3/GI

58

29

“Final” Guidelines

For k-e simulations: Two-layer is preferable over wall-functions (grid dependence + accuracy) Realizable k-e or Kato&Launder limiter have to be used For k-w simulations: SST is usually better than standard Grid should be clustered at wall SA is usually a good compromise between accuracy and simplicity. It also has very good convergence properties (as compared to the two-equation models) Reynolds stress model is expensive and it require a good initial guess (typically a k-e-type simulation)

ME469B/3/GI

59

Summary of turbulence models in Commercial Codes

Zero equation

One equation

Two equation

RSM

Non-Linear Models

Custom

FLUENT

y

y

y

y

n

y

StarCD

n

n

y

n

y

y

CFX

y

y

y

y

n

y

ME469B/3/GI

60

30