Week 9 - FVM Literature Review

Objective -

To describe the need for interpolation schemes and flux limiters in FVM.

Basics

In computational fluid dynamics there is 3 numerical method-based solvers to solve the governing Navier strokes equations these are:

1. Finite difference method
2. Finite element method
3. Finite volume method

In finite difference method:

The grid is linear and aligned so that all grid points lie in a straight line. The main thing is that the properties are considered to concentrate at a single point.

In finite volume method:

The method is applied over a volume such that the governing equation is valid over the entire volume, as the properties are assumed to be concentrated at the geometric center of the volume concerned.

In finite element method:

The properties are valid over an elemental area. the properties assumed to be concentrated at the area center.

Drawbacks of Finite Difference Method:

1. we cannot use FDM directly on a body whose grid points are not aligned properly along the x, and y directions. these may occur due to complex folds, and cuts in the structure.
2. as it is not aligned we cannot calculate the properties at the next adjacent point exactly, but we can make some correction factor and approximate it, this method is called as \' The Variable Transformation Method\' which basically seeks to create a grid using different variable space which looks very similar to a cartesian grid.
3. variable transformation method gives an approx solution of the property but it is time-consuming and complex.
4. for complex profiles, this method is less preferred.

That's why the Finite Volume Method (FVM) is used

Finite Volume Method

The finite volume method (FVM) is a method for representing and evaluating partial differential equations in the form of algebraic equations. In the finite volume method, volume integrals in a partial differential equation that contain a divergence term are converted to surface integrals, using the divergence theorem. These terms are then evaluated as fluxes at the surfaces of each finite volume. Because the flux entering a given volume is identical to that leaving the adjacent volume, these methods are conservative. Another advantage of the finite volume method is that it is easily formulated to allow for unstructured meshes. The method is used in many computational fluid dynamics packages. \"Finite volume\" refers to the small volume surrounding each node point on a mesh.

In FVM instead of writing the equation in a single point, equations are written in an entire volume. It means for a considered volume no assumptions are made over the shape of the volume. Inside the particular volume, the flow is conserved, flow quantities are conserved, energy is conserved, mass is conserved, momentum is conserved, and so on. This is called as \' Finite Volume Representation\' which begins with \' Integral form of governing equation\'.

Steps:-

1. once we have this integral form of an equation.
2. Apply Gauss divergence theorem
3. We get a discretized form which can be solved using the computer (\"Flux Terms\").

consider a steady state heat conduction equation with heat generation:

∂∂x(α.∂T∂x)+S=0∂∂x(α.∂T∂x)+S=0

S is the source term

ααis thermal diffusivity term

∂T∂x∂T∂x is the temperature gradient along the x direction.

Governing equation - ∂∂x(α.∂T∂x)+S=0∂∂x(α.∂T∂x)+S=0 is valid at every point in the domain, which means it is also valid in the entire control volume.

If a small control volume (dv) is considered and it is multiplied with the governing equation and integrated it for the entire volume then that should be also equal to zero i.e the \' integral form of governing equation\'.

Taking the conservation equation validated at a single point and if it is validated at single point then it is validated at the entire volume and in this entire volume properties are conserved.

Multiplying the equation by \'dv\' and integrating it -

∫(∂∂x(α.∂T∂x)+S)dv=0∫(∂∂x(α.∂T∂x)+S)dv=0

∫(∂∂x(α.∂T∂x)dv+∫Sdv=0∫(∂∂x(α.∂T∂x)dv+∫Sdv=0

∵dv=A.dx∵dv=A.dx

∫∂∂x(α.∂T∂x)A.dx+∫Sdv=0∫∂∂x(α.∂T∂x)A.dx+∫Sdv=0

[αA(∂T∂x)]ew+¯¯¯Sdv=0[αA(∂T∂x)]we+S¯dv=0

¯¯¯SS¯ is the average value of source term (S) which when multiplied by dv is going to be equal to ∫Sdv∫Sdv.

αeA∂T∂xe−αwA∂T∂xw+¯¯¯Sdv=0αeA∂T∂xe-αwA∂T∂xw+S¯dv=0

If we assume αα as a constant and not a function of temperature,

αA[Te−TpΔx]−αA[Tp−TwΔx]+¯¯¯Sdv=0αA[Te-TpΔx]-αA[Tp-TwΔx]+S¯dv=0

The source term which is to be added in a cell, if it is a function of temperature we have to linearize it. If ααis a function of temperature we have to linearize this also.

(¯¯¯Slinearization\–––––––––––––––→Su+SpTp)(S¯linearization\̲→Su+SpTp)

αA[Te−TpΔx]−αA[Tp−TwΔx]+¯¯¯Sdv=0αA[Te-TpΔx]-αA[Tp-TwΔx]+S¯dv=0

[Heatflux]out\−[Heatflux]in\+Sourceterm=0[Heatflux]out\-[Heatflux]in\+Sourceterm=0

The amount of heat that is removed from a particular box is conserved.

Different schemes are used to solve the equation:

1. Gradient scheme
2. Interpolation scheme

1. Gradient Scheme

From the above figure when we calculate the temperature gradient in the east face, we basically use temperature from neighboring points trying to evaluate or numerically represent the gradient of temperature is called as Gradient scheme.

2. Interpolation scheme

In the finite volume formulation, evaluation of surface and volume integrals sometimes requires the values of unknown variables at locations other than the computational nodes of the control volume. These values are obtained by employing certain interpolation schemes.

As an example let us consider the simple case of steady one-dimensional convection and diffusion.

ddx(ρuϕ)=ddx(Γdϕdx)ddx(ρuϕ)=ddx(Γdϕdx)

where ϕϕis our solution variable (a scalar), ρρ is the density, u is the velocity and ΓΓ is the diffusion coefficient.

The above figure shows the three neighboring control volumes, which we are going to use for deriving the discretization equation.

Integration of the governing differential equation over the control volume around grid point P yields:

(ρuϕ)e−(ρuϕ)w=(Γdϕdx)e−(Γdϕdx)w(ρuϕ)e-(ρuϕ)w=(Γdϕdx)e-(Γdϕdx)w

Evaluation of the left-hand side of the above equation requires the values of ϕϕ over faces e and w. These values are estimated by employing certain interpolation schemes.

Some interpolation schemes used in CFD are as follows:

1. Upwind interpolation scheme
2. Linear interpolation scheme
3. Quadratic Upwind Interpolation for Convective Kinematics (QUICK) scheme
4. Hybrid Scheme
5. Total Variation Diminishing (TVD) scheme

Explaining the Upwind interpolation scheme for an example

The Upwind interpolation scheme takes into account the direction of fluid flow and based on it, the value of the solution variable at the face is approximated by its value at the upstream computational node. For instance, if the fluid is moving from west to east (refer to figure), then:

ϕe≈ϕPϕe≈ϕP

And if the fluid is moving from east to west (refer figure), then:

ϕe≈ϕEϕe≈ϕE

In addition to these Gradient scheme and Interpolation schemes there are Flux Limiters.

Flux Limiters

While calculating the variable gradients they tend to toggle and are not bounded to a particular value and keep jumping, this occurs due to discontinuity due to shock occurring, having a flux limiter will help to bound the gradient value thereby giving a stable solution.

Numerical approximation using higher-order spatial discretization schemes results in better accuracy but exhibits spurious oscillations near discontinuities such as shocks. On the other hand, lower order spatial discretization schemes remain stable near discontinuities but dampen out the distributions because of their dissipative nature. Hence in order to capture discontinuities, high-resolution schemes are desired - which are highly accurate and stable throughout the domain.

Flux limiters operate functionally as a switch and tune the numerical fluxes obtained from lower order scheme and higher order scheme such that near the discontinuity lower order (stable) scheme is employed and away from it higher order (accurate) scheme is employed.

Suppose represents flux at face e calculated using a high-resolution scheme, LeLe represents flux at e calculated using a lower order scheme and HeHe represents flux at e calculated using a higher order scheme. Then ψ(r)ψ(r) represents the flux limiter function such that:

Fe=Le+ψ(r)(He−Le)Fe=Le+ψ(r)(He-Le)

Flux limiter function depends upon the smoothness parameter r which is defined as:

rP=ϕP−ϕWϕE−ϕPrP=ϕP-ϕWϕE-ϕP

The flux limiter function is constrained to be greater than or equal to zero. It assumes a value of 0 near discontinuities and a value of 1 away from them. Hence it switches between higher order scheme and a lower order scheme for numerical representation of flux, thereby keeping the scheme stable and accurate throughout.

Note that flux limiters are also referred to as slope limiters because they both have the same mathematical form, and both have the effect of limiting the solution gradient near shocks or discontinuities. In general, the term flux limiter is used when the limiter acts on system fluxes and slope limiter is used when the limiter acts on system states (like pressure, velocity etc.).

consider 1D- semi-discrete scheme

duidt+1Δxi[F(ui+12)−F(ui−12)]=0duidt+1Δxi[F(ui+12)-F(ui-12)]=0

F(ui+12)&F(ui−12)F(ui+12)&F(ui-12) these are edge fluxes for the ith cell

F(ui+12)=flowi+12−ϕ(ri)(flowi+12−fhighi+12)F(ui+12)=fi+12low-ϕ(ri)(fi+12low-fi+12high)

F(ui−12)=flowi−12−ϕ(ri)(flowi−12−fhighi−12)F(ui-12)=fi-12low-ϕ(ri)(fi-12low-fi-12high)

fhighfhigh - low precision, high-resolution flux

flow- high precision, low-resolution flux

ϕ(r)ϕ(r)- flux limiter function

ri=ui−ui−1ui+1−uiri=ui-ui-1ui+1-ui

\'r\' represents the ratio of successive gradients on the solution mesh.

ϕ(r)≥0ϕ(r)≥0

this is the limiter function case

if the limiter function = 0 , then there will be a low resolution and sharp gradient.

if limiter function > 0 , then there will be a high resolution and +ve gradient.