Week 9 - FVM Literature Review
Objective
To describe the differences between finite difference method (FDM) and finite volume method (FVM) and to describe the need for interpolation schemes and flux limiters in FVM
Finite Difference Method
It is a numerical method used in CFD to solve the Navier Stokes equation. It makes use of discrete points around the grid point at which the solution is to be calculated.
The grid mesh is divided into a finite number of points, and the solution of the Navier Stokes equation at a particular grid point is based on Taylor series approximations of values at neighbouring points.
A few well known finite difference schemes are iterative methods like Gauss Seidel Jacobi and SOR.
Drawbacks of Finite Difference method
1. It cannot be directly used to a body whose grid points are not aligned properly along the x,y directions. These may occur due to folds,or cuts in the structure
2. Therefore, the properties at the adjacent grid points cannot be exactly calculated. Rather, a correction factor needs to be introduced for approximating the difference. This method is known as Variable Transformation Method, which creates a grid using a different variable space for approximation, and then it is transferred back to the original grid.
3. Variable Transformation is a complex and time consuming process, which may be close to impossible in case of 2D or 3D simulations which involve a moving grid.
Therefore, the Finite Volume Method is preferred for this.
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 FVM, volume integrals in a partial differential equation are converted to surface integrals, using the divergence theorem. These terms are then evaluated as fluxes at the surface of each finite volume.
Since the flux entering a given volume is identical to that leaving the adjacent volume, these methods are conservative. Another advantage of FVM is that it can be easily formulated to allow for unstructured meshes.
In FVM, equations are written in an entire volume, instead of at a single point. No assumptions are made over the shape of the volume. Inside the considered volume, the flow is conserved. This is known as Finite Volume Representation. The governing equation is in the integral form.
Steps in the FVM
1. Obtain the governing equation in integral form.
2. Apply Gauss divergence theorem
3. A discretized form which can be solved - Flux terms
Example
Consider a 1D steady state heat conduction equation
`del /(del x) (alpha. (delT)/(delx)) + S = 0`
S is the source term
`alpha` is thermal diffusivity (function of temperature)
`(delT)/(delx)` is the temperature gradient along x direction
Considering finite volumes
Governing equation - `del /(del x) (alpha. (delT)/(delx)) + 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 multiplied with the governing equation and integrated over the entire volume, it gives the integral form of the governing equation
`int(del/(delx)(alpha.(delT)/(delx))+S) dv = 0`
`int(del/(delx)(alpha.(delT)/(delx))dv+ intS dv = 0`
`because dv = A .dx`
`intdel/(delx)(alpha.(delT)/(delx))A.dx+ intS dv = 0`
`[alphaA((delT)/(delx))]_w^e +overlineS dv = 0`
`overlineS` is the average value of the source term (S), which when multiplied by dv is equal to `intSdv`
`alpha_eA(delT)/(delx_e)-alpha_wA(delT)/(delx_w) +overlineS dv = 0`
Assuming `alpha` as a constant
`alpha A[(T_e-T_p)/(Deltax)]-alphaA[(T_p-T_w)/(Deltax)] + overlineSdv =0`
The terms inside the square braces are the heat flux terms. The source term needs to be linearized, if it is a function of temperature
`(overline S underline("linearization") to S_u+S_pT_p)`
`alpha A[(T_e-T_p)/(Deltax)]-alphaA[(T_p-T_w)/(Deltax)] + overlineSdv =0`
`[Heat flux]_("out")-[Heat flux]_("in")+Source term =0`
The heat flow in a particular volume is conserved.
Different schemes are used to solve the equation
1. Gradient scheme
2. Interpolation schemes
1. Gradient scheme
When we calculate the temperature gradient in the east face, we use temperature from the neighbouring points. This is known as the gradient scheme.
2. Interpolation schemes
In FVM formulation, evaluation of surface and volume integrals may sometimes require 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, a simple 1D convection and diffusion is described.
`d/(dx)(rho u phi)=d/(dx)(Gamma(dphi)/(dx))`
where `phi` is the solution variable (scalar), `rho` is the density, `u` is the velocity and `Gamma`, the diffusion coefficient.
The above figure shows the three neighbouring control volumes, used for discretization.
Integration of the governing differential equation over the control volume around grid point P yields
`(rho u phi)_e - (rho u phi)_w = (Gamma(dphi)/(dx))_e - (Gamma(dphi)/(dx))_w`
The evaluation of the LHS of the above equation requires the value of `phi` over the faces `e` and `w`. These values are estimated by employing certain interpolation schemes.
1. Upwind interpolation scheme
The upwind interpolation scheme takes into account the direction of fluid flow, and based on it, the value of solution variable at the face is approximated by its value at the upstream computational node.
`underset(e)f = { (underset(p)f, if underset(e)((v.n))>0),(underset(E)f, if underset(e)((v.n))<0):}`
It is equivalent to applying a backward or forward finite difference approximation. It is first order accurate.
The upwind relation is given below.
`underset(e)Gamma^(n um) =underset(e)((rho*u))*(Deltax)/2`
2. Linear interpolation
To approximate the value of the variable at control volume face center by using two nearest computational nodes. It is second order accurate.
It is equal to central difference schemes of the first order derivative.
`underset(e)f = underset(E)f*underset(e)(lambda)+underset(p)f*(1-underset(e)lambda)`
3. Quadratic upwind interpolation
It is third order accurate, but is prone to oscillations and is derived from polynomial fittings. The interpolation on a uniform cartesian grid is given by
`underset(e)(phi)= 6/8*underset(U)(phi) +3/8*underset(D)(phi)-1/8*underset(UU)(phi)`
where D, U and UU denote the downstream, the first upstream and second upstream nodes respectively.
In regions with strong gradients, overshoots and undershoots can lead to stability problems.
Flux limiters
They operate functionally as a switch and tune the numerical fluxes obtained from lower order schemes and higher order schemes such that a lower order (stable) scheme is employed near the discontinuity and a higher order (accurate) scheme away from it.
If `F_e`represents flux at face e calculated using the high resolution scheme, `L_e`represents flux calculated at e using lower order scheme and `H_e` represents flux calculated using higher order scheme.
`psi(r)` represents the flux limiter function such that
`F_e = L_e + psi(r)(H_e-L_e)`
The flux limiter function depends upon the smoothness parameter r defined as
`r_P = (phiP-phiW)/(phiE-phiP)`
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 and lower order schemes for numerical representation of flux, thereby keeping the scheme stable and accurate throughout.
The flux limiters are also referred to as slope limiters, since they 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 the system fluxes and slope limiter is used, when the limiter acts on system states (like pressure, velocity etc)
Consider 1D semi discrete scheme
`(du_i)/(dt)+1/(Deltax_i)[F(u_(i+1/2)) - F(u_(i-1/2))] = 0`
`F(u_(i+1/2)) & F(u_(i-1/2))` are edge fluxes for the "i"th cell
`F(u_(i+1/2)) = f_(i+1/2)^(low) - phi(ri)(f_(i+1/2)^(low)-f_(i+1/2)^(high))`
`F(u_(i-1/2)) = f_(i-1/2)^(low) - phi(ri)(f_(i-1/2)^(low)-f_(i-1/2)^(high))`
`f^(high)`- low precision, high resolution flux
`f^(low)` - high precision, low resolution flux
`phi(r)` - flux limiter function
`r_i = (u_i-u_(i-1))/(u_(i+1)-u_i)`
r represents the ratio of successive gradients on the solution mesh.
`phi(r)ge 0`
This is the limiter function case
If the limiter function = 0, there will be a low resolution and sharp gradient.
If the limiter function > 0, there will be a high resolution and positive gradient.