Finite Volume literature review

Objective

Describe the need of interpolation schemes and flux limiters in FVM.

Finite Volume method

Finite volume method is a numerical technique that transform partial differntial equations representing conservation laws over differential volumes into discrete algebraic equations over finite volumes ( or elements or cells). In a similar fashion to the finite difference or finite element method, the first step in the solution process is the discretization of the geometry domain, which in FVM, is discretized into non-overlapping elements or finite volumes. The partial differential equationa are then discretized/transformed into algebraic equations by integrating them over each discrete element. The system of algebraic equations is then solved to compute the values of the dependent variable for each of the elements.

In finite volume method, some of the terms in the conservation equation are turned into the face fluxes and evaluated at the finite volume faces. Because the flux entering a given volume is identical to that leaving the adjacent volume, the FVM is strictly conservative. This inherent conservayive property of FVM makes it preferred method in CFD. Another important attribute of the FVM is that it can be formulated in the physical space on unstructured polygonal meshes. Finally in FVM it is quite easy to implement a variety of boundary conditions in a noninvasive manner, since the unknown variables are evlaluated at the centroids of the volume elements, not at their boundary face.

These characteristics have made the finite volume quite suitable for the numerical simulation of a variety of applications involving fluid flow and heat and mass transfer, and developments in the method have been closely entwined with advanced in CFD. From a limited potential at inception confined to solving simple physics and geometry over structured grids, the FVM is now capable of dealing with all kinds of complex physics and applications.

similar to other numerical methods developed for the simulation of fluid flow, the finite volume method transforms the set of partial differential equations into a set of linear algebraic equations. The discretization procedure used in the finite volume method is distinctive and involves basic two step. In the first step, the partial differential equations are integrated and transformed into balance equation over an element. This involves changing the surface and volume integrals into discrete algebraic relations over elements and their surfaces using an integration quadrature of a specified order of accuracy. The result is a set of semi-discretized equations. In the second step, interpolation profiles are chosen to approximate the variation of of the variables within the element and relate the surface values of the variables to their cell values and thus transform the algebraic relations into algebraic equations.

Discretization process

• The numerical solution of partial differential consists of finding the values of the dependent variable `phi`at specifed points from which its distribution over the domain of interest constructed. These points are grid elements, or grid nodes and resulting from the discretization of original geometry into a set of non overlapping discrete elements, a process kniwn as meshing.
• The resulting nodes or variables are generally positioned at cell centroids or at vertices depending on the adopted discretization procedure.
• In all methods the focus is on replacing the continuous exact solution of the PDE with discrete values. The distribution `phi`of is hence discretized, and it is appropriate to refer this process of converting the governing equation into a set of algebraic equations for the discrete value of `phi`as the discretization process and the specific method employs to bring about this conversion as the discretization methods.
• The discrete value of `phi`are typically computed by solving a set of algebraic equations relating the values at neighboring grid elements to each other

Step 1: Geometric and physical modeling

• Modeling of physical phenomena is in a way at the heart of the scientific enterprise. A physical phenomena cannot generally be considered as understood unless it can be mathematically formulated, and this formulation tested and validated.
• For our purpose two levels of modeling are performed, one in relation to the gepmetry of the physical domain and a second in relation to the physical phenomena of interest.

Step 2:Domain Discretization

• The geometry discretization of the physical domain results in a mesh on which the conservation equations are eventually solved.

Step3: Equation Discretization

In step 3 the governing PDE , are transformed into a set of algebraic equations, one for each element in the computational domain.

Finite Volume method for one dimensional steady state equation

• Write the governing equation in its integral form
• Apply gauss divergence theorem
• come up wih suitable methods to calculate \"Flux terms\".

We will take the simple example

`frac{partial}{partialx}(alpha frac{partialT}{partialx})+S=0` (1)

Where S is the source

`alpha`is the thermal diffusivity

`frac{partialT}{partialx}`is the temperature gradient along the x-direction.

The first step in the finite volume method is to divide the domain into discrete volumes

The above governing equation is valid at every point in domain, which means it is also valid in entire control volume. If small control volume dv is considered and multiply by equation (1) and integrated it for entire volume then that should be also equal to zero.

`int(frac{partial}{partialx}(alphafrac{partialT}{partialx}+S))dv=0`

`intfrac{partial}{partialx}(alphafrac{partialT}{partialx})dv+intSdv=0`

`intfrac{partial}{partialx}(alphafrac{partialT}{partialx})Adx+intSdv=0`

i.e dv=Adx

`[alpha A frac{partialT}{partialx}]_w^e+barS dv=0`

`barS`is average value of source term (S) which when multiply by dv is going to be equal to `intSdv`

`a_eA frac{partialT}{partialx_e}-a_wAfrac{partialT}{partialx_w}+barSdv=0`

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

`alpha A[frac{T_e-T_p}{Deltax}]-alpha A[frac{T_p-T_w}{Deltax}]+barSdv=0`

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

(`barS`linearization`to``S_u+S_pT_p`)

`alpha A[frac{T_e-T_p}{Deltax}]-alpha A[frac{T_p-T_w}{Deltax}]+barSdv=0`

`[Heat flux]_(out)-[Heat flux]_(In)+Sourceterm=0`

1.Gradient Scheme

When we calculate the temperature gradient `[frac{partialT}{partialx}]_e`at the east face than we need the temp at the neighbouring point P and E and this called an gradient scheme.

2. Interpolation Scheme

Since the results of CFD simulations are typically stored in the cell centers, it is important to interpolate the results from cell centers to the face centers, to obtain the fluxes for the surface integrals in the transport equation. For each term of the transport equation, there is a veriety of discretization/interpolation schemes available.

In general, interpolation needs a flux F through a general face f, and in some cases, one or more parameters. The face value can be evaluated from the values in the neighboring cells using a variety of schemes. The flux satisfy continuity connstraints, which is prerequisite to obtaining the results.

Well known interpolation schemes are

1.First order Upwind scheme

2. Central differencing scheme

3.Quadratic Upwind interpolation for convective Kinematics (QUICK) scheme

4. Hybrid scheme

5. Total variation Diminishing (TVD) scheme.

1.First order upwind scheme

In first order upwind scheme we define `phi` as follows. Here we define two faces, e and w. To obtain flux through face e and w, we need to look its neighbouring values at P/E and W/P respectively. The subscript denote the face at which the face value `phi`or the flux F is located at.

`phi_e=phi_p` if `F_e>0`

`phi_w` is also defined similarly (positive direction is from W to E).

2. Central differencing scheme

Here, we linear interpolation for computing the cell face values.

`phi_e=frac{phi_E+phi_P}{2}` `phi_e=frac{phi_E+phi_P}{2}`

3.Quick Scheme

In Quick scheme 3 point upstream-weighted quadratic interpolation are used for cell face values.

When `F_e>0`, `phi_e=6/8phi_P+3/8phi_E-1/8phi_W`

When `F_w>0`, `phi_w=6/8phi_W+3/8phi_P-1/8phi_(WW)`

Similar expression can be obtained for `F_e<0 and F_w<0`.

4. The Hybrid differencing Scheme

The hybrid differencing scheme of spalding(1972) is based on the combination of central and upwind differencing scheme. The central differencing scheme which is senond order accurate, is employed for small Peclet numbers (`P_e<2`) and the upwind scheme, which is first order accurate but accounts for transportiveness, is employed for large Peclet numbers `P_ege2`. The hybrid differencing scheme uses piecewise formulae based on the local Peclet number to evaluate the net flux through each control volume.

5. The TVD Scheme

Scheme of third-order and above have been developed for the discretization of convective terms with varying degrees of success. Implementation of boundary conditions can be problematic with such higher order schemes. The fact that the QUICK scheme and other higher order schemes can give understood and overshoots has led to the development of second order schemes that avoid these problems. The class TVD scheme has been specially formulated to achieve oscillation-free solutions and has proved to be useful in CFD calaculations. TVD is a property used in the discretisation of equations governing time-dependent gas dynamic problems. More recently, schemes with this property have also become popular in general-purpose CFD solvers.

Flux Limiters

Flux limiters are used in numerical scheme to solve problems in science and engineering, particularly fluid dynamics, that are described by partial differential equation (PDEs). Their main purpose is to avoid the spurious oscillations (wiggles) that would otherwise occur with high-order spatial discretization due to shocks, discontinuities, or steep gradients in the solution domain.They can be used directly on finite differenc schemes for simple applications, such as the 1D advection PDE of eq. However, for more complex systems involving conservation laws, finite volume-resolution schemes are employed that are also designed to avoid entropy violations.

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.

The main idea behind the construction of flux limiter schemes is to limit the spatial derivatives to realistic values, for scientific and engineering problems, this usually means physically realizable values. They are used in high-resolution schemes for solving problems described by PDE\'s and only come into operation when sharp wave fronts are present. For smoothly changing waves, the flux limiters do not operate and the spatial derivatives can be represented by higher-order approximations without introducing non real oscillations. Consider the 1D semidiscrete scheme below

`frac{du_i}{dt}+frac{1}{Deltax_i}[F(u_(i+1/2))-F(u_(i-1/2))]=0`

Where for a finite volume scheme,

`F(u_(i+1/2))` and ` `F(u_(i-1/2))` represent the edge fluxes for the ith cell. Similarly, for a finite difference scheme, they represent flux values on the grid at point `x=x_(i+1/2)` and point `x=x_(i-1/2)`. If these fluxes can be represented by low and high resolution schemes, then a flux limiter can switch between these schemes depending upon the gradients close to the particular cell as follows:

where, `f_(low)=low-resolution flux`, `f_(high)= high resolution flux`, `phi(r)` = flux limiter function, and r represents the ratio of successive gradients on the solution mesh, i.e.,

`r_i= frac{u_i-u_(i-1)}{u_(i+1)-u_i}`

The limiter function is constrained to be greater than or equal to zero, i.e., `rge0`. Therefore, when the limiter is equal to zero (sharp gradient, opposite slopes, or zero gradient), the flux is represented by a low resolution scheme. Similarly, when the limiter is equal to 1 (smooth solution), it is represented by a high resolution scheme. The various limiters listed below have differing switching characteristics and are selected to suit the particular problem and numerical solution scheme. No particular limiter has been found to work well for all problems, and a particular choice is usually made on a trial-and-error basis.