FVM Literature
Objective :
To study the theory behind the various interpolation schemes and the flux limiters in case of the Finite Volume Method.
Why the need for FVM?
In advanced CFD approaches with highly unstructured grids, we tend to use a method called as the finite Volume Method ( FVM ). While in the normal discretization of the governing equations we tend to use the finite difference schemes the job becomes rather mundane while doing so for the unstructured grids because of the extra interpolations. For example, consider a 1-D mesh of non-uniform size as is shown below:
As is seen in the figure above, the mesh grid points are not uniformly spaced. This type of situation may arise if we want to capture some phenomenon that is predominant in one region and is very less important in the other. Likewise, we tend to create the mesh in accordance with our need for accuracy. For a physical intuition, one can think of an example of an aeroplane mesh. After long distances, the flow tends to normalize to the free stream flow and hence after long distances, one can create grids which are large in size hence saving unnecessary computation time.
Returning back to our unstructured 1-D grid, we need to apply the Taylor Series approximation at each point in general if we want to discretize the equation via the finite difference method. As the mesh is 1-D it is not so hard to do it, but just imagine doing it for highly unstructured 3-D grids like the ones wrapped around the airfoil or around a car to simulate the flow field. For such cases, the FDM method will not be the best choice. A better alternative is to use the FVM approach.
What is FVM?
FVM stands for the Finite Volume Method. As the name suggests we use the control volume integrated with the governing equation. Consider the following Unstructured mesh and the governing equation for the heat conduction in 1-D with the source term. The dotted line surrounding P represents the control volume and W, E represent the centres of the adjacent control volumes
Governing Equation: `frac(partial)(partialx)(kfrac(partialT)(partialx)) + S = 0`
Now as the governing equation is valid for each and every point in the domain it is valid over an entire control volume too. So hence we can integrate the governing equation for a control volume as:
`int int int _V [frac(partial)(partialx)(kfrac(partialT)(partialx)) + S]dV =0`
As our control volume is a 1-D mesh the volume integral reduces to a 1-D integral and we get:
`int \quad[frac(partial)(partialx)(kfrac(partialT)(partialx)) + S]Adx = 0`
So if we integrate the equation from West to East as is shown in the control volume we get :
`int _w^e[frac(partial)(partialx)(kfrac(partialT)(partialx)) + S]Adx = 0`
Splitting the integrands we get :
`int _w^e[frac(partial)(partialx)(kfrac(partialT)(partialx))]Adx +int_w^e SAdx =`
Solving the equation we get:
`(Akfrac(partialT)(partialx))_w^e +int_w^e SAdx = 0`
Now the term S can be linearised in order to facilitate the easy solution of the equation. Picard's method is the most popular method used for the linearization. For a given control volume P, we start by writing the source term as:
`S = S_C + S_PT_P`
This method of linearizing the source term has an added advantage that if at all the source term present is highly non-linear with several powers of T it can be easily represented in the linear form:
For example consider the case of a source term varying with the temperature as:
`S = 3 - 4T^3` . To linearize it we consider the form as `S = S_C + S_PT_P`.
it is up to us to choose the terms `S_C and S_P`. One way would be to choose
`S_C = 3 and S_P = -(4T^**)^2` where `T^**` is the value from the previous iteration.
Returning back to our governing equation:
`(Akfrac(partialT)(partialx))_w^e +int_w^e SAdx = 0`
`(Akfrac(partialT)(partialx))_e - (Akfrac(partialT)(partialx))_w+int_w^e (S_C + S_PT_P)Adx = 0`
The advantage of writing the source term as a linearized form is that now we can capture the effect of the source term integrated with the Temperature at the control volume.
Simplifying the terms a bit we get:
`Ak_e(frac(partialT)(partialx))_e - Ak_w(frac(partialT)(partialx))_w+int_w^e (S_C )Adx+ int_w^e(S_PT_P)Adx = 0`
We can write:
`Ak_e(frac(partialT)(partialx))_e - Ak_w(frac(partialT)(partialx))_w+S_CADeltax+ (S_PT_P)ADeltax = 0`
The main difficulty lies in the fact of the approximation of the derivative terms like `(frac(partialT)(partialx))_e`. It is imperative to use a robust discretization/interpolation scheme so as to capture the effects nicely.
There are several interpolation schemes available out of which we can choose anyone to work with.
1.) Finite Difference scheme:
This is a simple scheme using the finite difference approximation with the help of Taylor Series:
According to the control volume, the terms will become:
`(frac(partialT)(partialx))_e = frac(T_E-T_P)(Deltax)` and
`(frac(partialT)(partialx))_w = frac(T_P-T_W)(Deltax)` where E and W are the corresponding neighbours.
So we can write the final equation as:
`k_e(frac(T_E-T_P)(Deltax)) - k_w(frac(T_P-T_W)(Deltax))+S_CDeltax+ (S_PT_P)Deltax = 0`
Writing the terms as :
`frac(alpha_e)(Deltax) = a_E` and `frac(alpha_w)(Deltax) = a_W`. So then writing the equation in these terms we can write it as :
`a_ET_E - a_ET_P - a_WT_P + a_WT_W + S_CDeltax + S_PT_PDeltax = 0`
Simplifying the equation we get:
`a_ET_E + a_WT_W + S_CDeltax = a_ET_P + a_WT_P + S_PT_PDeltax`
We can write it as:
`a_PT_P = a_ET_E + a_WT_W + b` in the simpified form for a given control volume and hence this equation can then be solved algebraically to obtain the temperature at the grid point P.
However, this scheme has a downfall. It will not work if by chance `a_P` becomes zero. So we need a more powerful scheme to account for this behaviour. Please note that the scheme was derived via the assumption that the grid sizes are uniform. It can also be derived for the non-uniform grids with the use of appropriate terms as is shown below:
Replacing the derivative terms with accurate enough discretized terms we can modify the equation for the unstructured grids.
2.) Upwind Scheme
This is a new type of scheme in which the value at the interface is taken according to the convection current:
If a given flow field is towards the right then we take the value of the interface as the previous grid point and if the flow field is towards left we take the value at the interface from the right. That is
`phi_e = phi_P\ \text( if positive convection)` and `phi_e = phi_E\ \text( if negative convection)`
This type of scheme is predominant in the fluid flow problems where there are large convection currents in one direction.
3.) Exponential Scheme:
It takes the behaviour inside the CV as an exponential as is obtained by the exact solution for the 1-D convection-diffusion problem and the idea is extended to the higher dimensions, i.e.
`frac( phi - phi_P )(phi_E - phi_P) = frac( e^(kx) - 1 )(e^a - 1)` where `phi` is the value inside a control volume at a distance x and `phi_p and phi_E` are the values of the variable at the adjacent grid points.
This scheme is also used for fluid flow problems involving convection and diffusion.
Similar to this there can be several other schemes like the hybrid, power-law schemes, quadratic interpolation schemes the details of which are out of the scope of this report but is attached in the report below for reference.
So to answer why do we need the interpolation schemes?
We need the interpolation schemes to approximate the property values located at the boundaries of the control volumes which appear in the equations after the integration over the control volume that the FVM can be applied.
Low-order schemes are commonly considered, due to their numerical stability whereas if we need more accuracy higher-order methods are preferred. The preference for the interpolation scheme is crucial in solving the complex coupled non-linear PDE's.
Now to answer the question why do we need Flux Limiters?
Flux limiters are used in the numerical schemes used to solve problems of fluid dynamics, described by highly coupled non-linear Partial Differential Equations. If we tend to use lower order numerical schemes for the solving of the governing equations then we might get a highly oscillatory and unstable solution near the discontinuity if present. Whereas if we tend to use higher-order numerical schemes to capture the phenomenon more accurately we tend to exploit the computational power unnecessarily. The main idea behind the use of the flux limiters is to limit the spatial derivative terms near the discontinuities to physical real values so as to maintain physical consistency. There are several types of flux/slope limiter functions. They come in handy when sharp wave-fronts/peaks are present. For smoothly changing waves/functions, the flux limiters do not operate and the spatial derivatives can be represented normally. It will be called a fool's method to use flux limiters when the user knows that the functions are smooth in the general space.
if put in a sophisticated way, in the implementation we tune the numerical flux of higher and lower order scheme using the flux limiter function in such a way that the resulting scheme gives a high order accuracy in the smooth region of flow and sticks with a lower order of accuracy in the vicinity of socks/discontinuities so as to get a physically consistent solution. These type of methods are predominant in the use of the Bluntnose problems in Aerodynamics where there are shock discontinuities.
It is imperative to 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.
Conclusion:
Finite Volume Method also abbreviated as FVM is a really powerful approach to solve the CFD problems. With the help of proper interpolation schemes and the accurate use of flux limiters in case of any discontinuities present if any, we can solve a wide range of CFD problems. Though it may be difficult to code and implement it in the computer it gives an added advantage of its robustness