3y ago

63 Views

2 Downloads

3.78 MB

86 Pages

Transcription

Solution methods for theIncompressible Navier-Stokes Equations Discretization schemes for the Navier-Stokes equationsPressure-based approachDensity-based approachConvergence accelerationPeriodic FlowsUnsteady FlowsME469B/3/GI1

Background (from ME469A or similar)Navier-Stokes (NS) equationsFinite Volume (FV) discretizationDiscretization of space derivatives (upwind, central, QUICK, etc.)Pressure-velocity coupling issuePressure correction schemes (SIMPLE, SIMPLEC, PISO)Multigrid methodsME469B/3/GI2

NS equationsConservation laws:Rate of change advection diffusion source 0ME469B/3/GI3

NS equationsDifferential form:0The advection term is non-linearThe mass and momentum equations are coupled (via the velocity)The pressure appears only as a source term in the momentum equationNo evolution equation for the pressureThere are four equations and five unknowns (ρ, V, p)ME469B/3/GI4

NS equationsCompressible flows:The mass conservation is a transport equation for density. With an additionalenergy 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 aconstraint on the velocity field; this equation (combined with the momentum) canbe used to derive an equation for the pressureME469B/3/GI5

Finite Volume MethodDiscretize the equations in conservation (integral) formEventually this becomes ME469B/3/GI6

Pressure-based solution of the NS equationThe continuity equation is combined with the momentum and thedivergence-free constraint becomes an elliptic equation for the pressureTo clarify the difficulties related to the treatment of the pressure, wewill define EXPLICIT and IMPLICIT schemes to solve the NS equations:It is assumed that space derivatives in the NS are already discretized:ME469B/3/GI7

Explicit scheme for NS equationsSemi-discrete form of the NSExplicit time integrationThe n 1 velocity field is NOT divergence freeTake the divergence of the momentumElliptic equation for the pressureME469B/3/GI8

Explicit pressure-based scheme for NS equationsVelocity field (divergence free) available at time nCompute HnSolve the Poisson equation for the pressure pnCompute the new velocity field un 1ME469B/3/GI9

Implicit scheme for NS equationsSemi-discrete form of the NSImplicit time integrationTake the divergence of the momentumThe equations are coupled and non-linearME469B/3/GI10

Navier-Stokes EquationsConservation of massConservation of momentumConservation of energyNewtonian fluidIn 3D: 5 equations & 6 unknowns: p, ρ, vi, E(T)Need supplemental information: equation of stateME469B/3/GI11

ApproximationsAlthough the Navier-Stokes equations are considered the appropriateconceptual model for fluid flows they contain 3 major approximations:1. Continuum hypothesis2. Form of the diffusive fluxes3. Equation of stateSimplified conceptual models can be derived introducing additionalassumptions: incompressible flowConservation of mass (continuity)Conservation of momentumDifficulties:Non-linearity, coupling, role of the pressureME469B/3/GI12

A Solution ApproachThe momentum equation can be interpreted as a advection/diffusionequation for the velocity vectorThe mass conservation should be used to derive the pressure taking the divergence of the momentum:A Poisson equation for the pressure is derivedME469B/3/GI13

The Projection MethodImplicit, coupled and non-linearPredicted velocityassumingwe obtainbutand taking the divergencethis is what we wouldlike to enforcecombining (corrector step)ME469B/3/GI14

Alternative View of ProjectionReorganize the NS equations (Uzawa)LU decompositionExact splittingMomentum eqs.Pressure Poisson eq.Velocity correctionME469B/3/GI15

Alternative View of ProjectionExact projection requires the inversion of the LHS of the momentum eq.thus is costly.Approximate projection methods are constructed using two auxiliarymatrices (time-scales)Momentum eqs.Pressure Poisson eq.Velocity correctionThe simplest (conventional) choice isME469B/3/GI16

What about steady state?Solution of the steady-state NS equations is of primary importanceSteady vs. unsteady is another hypothesis that requires formalization Mom. EquationsReference QuantitiesNon dimensional EqnReynolds and Strouhal #sME469B/3/GI17

Implicit scheme for steady NS equationsCompute an intermediate velocity field(eqns are STILL non-linear)Define a velocity and a pressure correctionUsing the definition and combining{{Derive an equation for u’ME469B/3/GI18

Implicit scheme for steady NS equationsTaking the divergence We obtain a Poisson system for the pressure correction Solving it and computing a gradient:So we can updateAnd also the pressure at the next levelME469B/3/GI19

Implicit pressure-based scheme for NS equations (SIMPLE)SIM PLE: Semi-Implicit M ethod for Pressure-Linked EquationsVelocity field (divergence free) available at time nCompute intermediate velocities u*Solve the Poisson equation for the pressure correction p’Neglecting the u*’ termSolve the velocity correction equation for u’Neglecting the u*’ termCompute the new velocity un 1 and pressure pn 1 fieldsME469B/3/GI20

Implicit pressure-based scheme for NS equations (SIMPLEC)SIM PLE: SIM PLE Corrected/ConsistentVelocity field (divergence free) available at time nCompute intermediate velocities u*Solve the Poisson equation for the pressure correction p’Use an approximation to u*’ (neighbor values average u*’ Σ u’)Solve the velocity correction equation for u’Use an approximation to u*’Compute the new velocity un 1 and pressure pn 1 fieldsME469B/3/GI21

Implicit pressure-based scheme for NS equations (PISO)PISO: Pressure Implicit with Splitting OperatorsVelocity field (divergence free) available at time nCompute intermediate velocities u* and p’ as in SIMPLESolve 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 fieldsME469B/3/GI22

SIMPLE, SIMPLEC & PISO - CommentsIn SIMPLE under-relaxation is required due to the neglect of u*’un 1 u* αu u’p pn αp p’There is an optimal relationshipαp 1- αuSIMPLEC and PISO do not need under-relaxationSIMPLEC/PISO allow faster convergence than SIMPLEPISO is useful for irregular cellsME469B/3/GI23

Under-relaxationIs used to increase stability (smoothing)Variable under-relaxationEquation (implicit) under-relaxationME469B/3/GI24

Segregated (pressure based) solver in FLUENTFV discretization for mixed elementsfΩThe quantities at the cell faces can be computed using several different schemesME469B/3/GI25

Discretization of the equationsOptions for the segregated solver in FLUENTDiscretization scheme for convective terms1st order upwind (UD)2nd order upwind (TVD)3rd order upwind (QUICK), only for quad and hexPressure 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 couplingSIMPLESIMPLECPISOME469B/3/GI26

Discretization of the convective termsDetermine the face valueinterpolatedvalueφ(x)φeφP1st Order UpwindφEDepending on the flow direction ONLYVery stable but dissipativePeEFlow directionME469B/3/GI27

Discretization of the convective termsDetermine the face valueinterpolatedvalueφ(x)φPCentral differencing (2nd order)φeSymmetric. Not depending on the flow directionPeφEENot dissipative but dispersive (odd derivatives)Flow directionME469B/3/GI28

Discretization of the convective termsDetermine the face valueφWinterpolatedvalueφ(x)φP2nd order upwindφeDepends on the flow directionLess dissipative than 1st order butnot bounded (extrema preserving)WPeφEEPossibility of using limitersFlow directionME469B/3/GI29

Discretization of the convective termsDetermine the face valueφWinterpolatedvalueφ(x)φPQuick (Quadratic Upwind Interpolationfor Convection Kinetics)φeφEeEFormally 3rd orderDepends on the flow directionWPAs before it is not boundedFlow directionME469B/3/GI30

Evaluation of gradientsGauss GradientLeast SquareGradientLS systemME469B/3/GI31

Solution of the equationφ is one of the velocity component and the convective terms must be linearized:This correspond to a sparse linear system for each velocity componentFluent segregated solver uses:Point Gauss-Seidel techniqueMultigrid accelerationME469B/3/GIbbPb32

GridsME469B/3/GIMultiblock structured - Gambit33

GridsME469B/3/GIHybrid non-conformal - Gambit34

GridsME469B/3/GIHybrid adaptive - non-Gambit35

GridsME469B/3/GIPolyhedral - non-Gambit36

Set-up of problems with FLUENTRead/Import the gridDefine the flow solver optionDefine the fluid propertiesDefine the discretization schemeDefine the boundary conditionDefine initial conditionsDefine convergence monitorsRun the simulationAnalyze the resultsGraphics WindowCommand MenusText WindowME469B/3/GI37

Solver set-upDefine Models SolverDefine Controls SolutionExample: text commands can be used (useful for batch execution)define/models/solver tion-scheme/mom 1solve/set/under-relaxation/mom 0.7 ME469B/3/GI38

Material propertiesDefine MaterialsQuantities are ALWAYS dimensionalME469B/3/GI39

Initial and boundary conditionsSolve Initialize InitializeDefine Boundary ConditionsOnly constant values can be specifiedMore flexibility is allowed via patchingBCs will be discussed case-by-caseME469B/3/GI40

Initial conditions using patchingAdapt Region MarkMark a certain region of the domain(cells are stored in a register)ME469B/3/GISolve Initialize PatchPatch desired values for each variablein the region (register) selected41

Convergence monitorsSolve Monitors ResidualsSolve Monitors SurfaceConvergence history of the equation residuals are stored together with the solutionUser-defined monitors are NOT stored by defaultME469B/3/GI42

PostprocessingDisplay ContoursPlot XY PlotCell-centereddata areComputedThis switchinterpolates theresults on thecell-verticesME469B/3/GI43

Detailed post-processingDefine additional quantitiesDefine plotting lines, planes and surfacesCompute integral/averaged quantitiesDefine Custom Field FunctionME469B/3/GI44

Fluent GUI - SummaryFile: I/OGrid: Modify (translate/scale/etc.), CheckDefine: Models (solver type/multiphase/etc.), Material (fluid properties),Boundary conditionsSolve: Discretization, Initial Condition, Convergence MonitorsAdapt: Grid adaptation, Patch markingSurface: Create zones (postprocessing/monitors)Display: Postprocessing (View/Countors/Streamlines)Plot: XY Plots, ResidualsReport: Summary, IntegralTypical simulationParallel: Load Balancing, MonitorsME469B/3/GI45

Example – Driven cavityClassical test-case forincompressible flow solversProblem set-upMaterial Properties:ρ 1kg/m3µ 0.001kg/msReynolds number:H 1m, Vslip 1m/sRe ρVslipH/µ 1,000Vslip 1Boundary Conditions:Slip wall (u Vslip) on topNo-slip walls the othersHSolver Set-UpSegregated SolverDiscretization:2nd order upwindSIMPLEMultigridV-CycleInitial Conditions:u v p 0Convergence Monitors:Averaged pressure andfriction on the no-slip wallsME469B/3/GI46

Example – Driven cavityThe effect of the meshing schemeQuad-Mapping 1600 cellsTri-Paving 3600 cellsQuad-Paving 1650 cellsEdge size on the boundaries is the sameME469B/3/GI47

Example – Driven cavityThe effect of the meshing scheme – Vorticity ContoursQuad-Mapping 1600 cellsME469B/3/GITri-Paving 3600 cellsQuad-Paving 1650 cells48

Example – Driven cavityThe effect of the meshing scheme – ConvergenceQuad-Mapping 1600 cellsME469B/3/GITri-Paving 3600 cellsQuad-Paving 1650 cells49

Example – Driven cavityThe effect of the meshing schemex-velocity component in the middle of the cavityQuad-Mapping Tri-PavingME469B/3/GIQuad-PavingSymbols corresponds toGhia et al., 198250

Example – Driven cavityGrid Sensitivity – Quad Mapping SchemeVorticity Contours1600 cellsME469B/3/GI6400 cells25600 cells51

Example – Driven cavityGrid Sensitivity – Quad Mapping Schemex-velocity component in the middle of the cavity1600 cellsME469B/3/GI6400 cells25600 cellsSymbols corresponds toGhia et al., 198252

How to verify the accuracy?Define a reference solution (analytical or computed on a very fine grid)Compute the solution on successively refined gridsDefine the error as the deviation of the current solution from the referenceCompute error normsPlot 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 grids2) Definition of the grid sizeME469B/3/GI53

Generation of successively refined grid1) Modify grid dimensions in GAMBIT and regenerate the grid2) Split all the cells in FLUENTAdapt Region AdaptThe region MUST containthe entire domainME469B/3/GIElement shape & metricproperties are preserved54

Driven Cavity - Error evaluationReference solution computed on a 320x320 grid ( 100,000 cells)Reference solution interpolated on coarse mesh to evaluate local errorsQuad-MappingTri-PavingQuad-PavingNote that the triangular grid has more than twice as many grid cellsME469B/3/GI55

Driven Cavity – Accuracy evaluation1stNominalorder accuracyNominal 2ndorder accuracyTri meshing scheme yieldsSlightly higher errors andlower accuracyNote that the definition of Δxis questionable (a change willonly translate the curves notchange the slope)Quad-Mapping Tri-PavingME469B/3/GIError (L2 norm)Quad and Pave meshingschemes yield very similaraccuracy (close to 2nd order)Quad-Paving(N)-1/256

Driven Cavity – Fluent vs. other CFD codesQuad Mapping Scheme (1600 cells)x-velocity component in the middle of the cavityFLUENTME469B/3/GIStarCDNASA INS2DSymbols corresponds toGhia et al., 198257

Techniques for the incompressible NS equationsPressure correction schemesArtificial compressibility approachVorticity-streamfunction formulationDensity-based approachME469B/3/GI58

Techniques for the incompressible NS equationsVorticity-streamfunction approachIt is effectively a change-of-variables; introducing the streamfunction and the vorticityvector the continuity is automatically satisfied and the pressure disappears (if needed thesolution of a Poisson-like equation is still required). It is advantageous in 2D because itrequires the solution of only two PDEs but the treatment of BCs is difficult. In additionin 3D the PDEs to be solved are sixArtificial compressibility approachA time-derivative (of pressure) is added to the continuity equation with the goal oftransforming the incompressible NS into a hyperbolic system and then to apply schemessuitable for compressible flows. The key is the presence of a user-parameter β (relatedto the artificial speed of sound) that determines the speed of convergence to steady stateME469B/3/GI59

Density-based solvers for the NS equationsThe equation are written in compressible form and, for low Mach numbers,the flow is effectively incompressibleThe energy equation is added to link pressureand density through the equation of stateIn compact (vector) form:ME469B/3/GI60

Density-based solvers for the NS equationsStiffness occurs because of the disparity between fluid velocity and speedof sound (infinite in zero-Mach limit)The equations are solved in terms of the primitive variablesNote that the continuity becomes (again) anevolution equation for the pressurewhereME469B/3/GI61

Density-based solvers for the NS equationsThe time derivative is modified (preconditioned) to force all the eigenvaluesto be of the same order (similar to the artificial compressibility approach) whereThe eigenvalues of Γ areME469B/3/GI62

Density-based solvers for the NS equationsLimiting casesCompressible flows(ideal gas):Incompressible flows (ideal gas):All eigenvaluesare comparableIncompressible fluids:ME469B/3/GI63

FLUENT density-based solverExplicit SchemeMultistage Runge-Kutta schemeResidual SmoothingMultigrid accelerationME469B/3/GI64

FLUENT density-based solverImplicit SchemeEuler (one-step) implicit with Newton-type linearizationPoint Gauss-Seidel iterationsMultigrid accelerationME469B/3/GI65

Example – Driven cavityClassical test-case forincompressible flow solversProblem set-upMaterial Properties:ρ 1kg/m3µ 0.001kg/msReynolds number:H 1m, Vslip 1m/sRe ρVslipH/µ 1,000Vslip 1Solver Set-UpCoupled SolverDiscretization:2nd order upwindImplicitMultigridV-CycleBoundary Conditions:Slip wall (u Vslip) on topNo-slip walls the othersHInitial Conditions:u v p 0Convergence Monitors:Averaged pressure andfriction on the no-slip wallsME469B/3/GI66

Example – Driven cavityEffect of the solver - Quad mesh (1600 cells)SegregatedCoupledVorticity ContoursME469B/3/GI67

Example – Driven cavityEffect of the solver - Quad mesh (1600 cells)x-velocity component in the middle of the cavitySegregatedME469B/3/GICoupledSymbols corresponds toGhia et al., 198268

Multigrid accelerationBasic idea: the global error (low-frequency) on a fine grid appears as a localerror (high-frequency) on coarse meshes.Why it is important: linear system solver like Gauss-Seidel are effective inremoving high-frequency errors but VERY slow for global errors. Notethat, on structured, grid line-relaxation (or ADI-type) schemes can be usedto improve the performance of Gauss-Seidel; on unstructured grid similarconcepts are extremely difficult to implement.Convergence Speed: number of iterations on the finest grid required toreach a given level of convergence is roughly independent on the number ofgrid nodes (multigrid convergence)ME469B/3/GI69

Two-grid scheme1.2.3.4.5.α smoothings are performed on the fine grid to reduce the highfrequency components of the errors (pre-smoothing, αS)the residual (error) is transferred to next coarser level (restriction, R)γ iterations are performed on this grid level for the “correction” equationthe problem is transferred back to the fine grid (prolongation, P)β smoothings are performed on the fine grid to remove the highfrequency errors introduced on the coarse mesh (post-smoothing, βS)Parameters to be defined are α, β, γME469B/3/GI70

Multigrid FormalismAfter few sweeps at level hTransfer (restrict) the residualModified system on the coarse gridTransfer (prolong) the solutionCorrectDefinition of the error and residualME469B/3/GI71

Restriction & Prolongation OperatorsFine LevelCoarse LevelME469B/3/GI72

Algebraic MultigridThe coarse levels are generated without the use of any discretization oncoarse levels; in fact no hierarchy of meshes is neededAMG is effectively a solver for linear systems and the restriction andprolongation operators might be viewed as means to modify (group or split)the coefficient matrixFormally:Geometric multigrid should perform better than AMG because nonlinearity of the problem are retained on coarse levels (correction equation)ME469B/3/GI73

Multigrid for unstructured meshesAggregative Coarsening: fine grid cells are collected into a coarse grid elementSelective Coarsening: few fine grid cells are retained on the coarser grids ME469B/3/GI74

Multigrid in FluentLevel Cycling: V, W and F (W V)V-CycleME469B/3/GIW-Cycle75

Multigrid in FluentFlexible CycleRestriction Criteria:A coarser level is invoked as soon as the residual reduction rate is below a certain %Termination Criteria:The corrections are transferred to a finer level as soon as a certain residual level is reachedME469B/3/GI76

Multigrid in FluentME469B/3/GI77

Algebraic Multigrid PerformanceConvergence for the segregated solver1600 cellsME469B/3/GI6400 cells25600 cells78

Algebraic Multigrid PerformanceConvergence for the coupled solver1600 cellsME469B/3/GI6400 cells25600 cells79

Period

Solution methods for the . Compressible flows: The mass conservation is a transport equation for density. With an additional . Define the flow solver option Define the fluid properties Define the discretization scheme Define the boundary condition Define initial conditions

Related Documents: