Week 9 - FVM Literature Review
Finite volume method is the most commonly used discretization technique in computational fluid dynamics.
It is on the basis of this discretization method most of the CFD tools are developed.
To start with FVM in simple words we can say that, given domain is sub divided into many small domains or specifically many finite volumes and the flow equations are solved inside each of these finite volume.
So broadly we can divide FVM discretization method into 3:
1. Dividing the given domain into number of small finite volumes.
2. Integrating the given flow equation inside each control volume and converting it into discretized form.
3. Solving the algebraic equations obtained from discretization and hence obtaining the distribution of flow parameter over the domian.
Above said steps can we well explained by using an example of 1d steady state heat equation,
1D steady state heat equation without source term (Heat generation) can be given as
`del/(delx)(K*A*(delT)/(delx)) = 0`
above heat equation is used to determine the tempertaure distribution in a rod as shown below:
let the thermal conducivity, K = 1000 w/m*k and Area, A = 10*10^-3 m^2
total length of 0.5 meter is sub divided into 5 finite control volumes each of length 0.1 meter, as shown below
consider finite volume with grid point 2 as shown below
control volume which is under consideration is with a central grid point P, left face of the control volume is given as w and right face as e
adjacent grid points, on the left hand side it is given as W and on the right hand side it is given as E.
above used letters to denote the faces and the grid points are the conventional symbols.
Distance b/w the faces e and w = `deltax`
Distance b/w P and w = `deltax_(w)`
Distance b/w P and e = `deltax_(e)`
since right side of the rod is subjected to high temperature heat flows from right side to left side of the rod,
integarting heat equation , `del/(delx)(K*A*(delT)/(delx)) = 0` within in control volume, i.e, integrating from e to w
`int_e^w del/(delx)(K*A*(delT)/(delx)) = 0`
i.e, `[K*A*(delT)/(delx)]_e^w = 0`
`[K*A*(delT)/(delx)]_w - [K*A*(delT)/(delx)]_e = 0`
now in the above equation `(delT)/(delx)`is determined at e and w using interpolation scheme.
physically heat flux entering the face e must flow continously through face w, there should not be any discontinuity, keeping this physical aspect in mind linear temperature profiles is considered between the grid points and it is as shown below,
this is where the point of flux limiters comes into picture. Instead of assuming a linear temperature profile between the grid points, we can assume constant temperature profiles inside each finite volume as shown below
but if we use this temperature distribution profile at the finite volume faces there is a discontinuity, but physically the heat flux through the faces are continuous and hence above shown constant temperature distrbution profile cannot be used because it is a flux limiter.
as it can be seen from the above image that between the faces e and w the heat flux is continuous
and hence, `[K*A*(delT)/(delx)]_w - [K*A*(delT)/(delx)]_e = 0` can be written as
`(K_(w)*A_(w)*((T_(P)-T_(W))/(deltax_(w)))) - (K_(e)*A_(e)*((T_(E)-T_(P))/(deltax_(e)))) = 0`
multiplying and rearranging we get,
`T_(P)*[(K_(w)*A_(w))/(deltax_(w))+(K_(e)*A_(e))/(deltax_(e))] = T_(W)*[(K_(w)*A_(w))/(deltax_(w))]+T_(E)*[(K_(e)*A_(e))/(deltax_(e))]`
which is expressed as, `T_(P)*(a_(w)+a_(e)) = a_(w)*T_(W) + a_(e)*T_(E)`
where `a_(w) = (K_(w)*A_(w))/(deltax_(w))`
`a_(e) = (K_(e)*A_(e))/(deltax_(e))`
and `a_(p) = a_(e) + a_(w)`
Now replacing `T_(P), T_(W), T_(E)` with `T_(2), T_(1), T_(3)` respectivley we get for node 2 the linear equation as
`T_(2)*a_(p) = a_(w)*T_(1) + a_(e) * T_(3)`
same is being applied for nodes 3,4 and hence we get resulting linear equation for node 3 and 4 as,
for node 3 , `T_(3)*a_(p) = a_(w) *T_(2) + a_(e) *T_(4)`
at node 4, `T_(4)*a_(p) = a_(w) * T_(3) + a_(e) * T_(5)`
Now taking boundary nodes 1 and 5 into consideration, unlike other nodes from 2 to 4 we cannot apply the same for nodes 1 and 5 since these are boundary nodes.
before proceeding into node 1 one point that is to be noted is that, in this example heat generation within the domain is taken as zero i.e,
if heat generation is considered then as per FVM heat generation S can be discretized as
`S = S_(u) + S_(p)*T_(p)`
now coming back to node 1, it can be represented as follow,
Heat flux flows from node 2 to left end A
integrating heat equation from e to A, `int_e^A del/(delx)(K*A*(delT)/(delx)) = 0`
which yields `[K*A*(delT)/(delx)]_e^A = 0`
`[K*A*(delT)/(delx)]_A - [K*A*(delT)/(delx)]_e = 0`
in the above equation `(delT)/(delx)`is determined at A and e using interpolation scheme and same linear approach as that used above is applied
and hence we have
`(K_(A)*A_(A)*(2*((T_(P)-T_(A))/(deltax)))) - (K_(e)*A_(e)*((T_(E)-T_(P))/(deltax_(e)))) = 0`
`T_(P)*[(2*K_(A)*A_(A))/(deltax)+(K_(e)*A_(e))/(deltax_(e))] = T_(A)*[(2*K_(A)*A_(A))/(deltax)]+T_(E)*[(K_(e)*A_(e))/(deltax_(e))]`
In the above equation , `a_(e) = (K_(e)*A_(e))/(deltax_(e))`
`a_(w) = 0`
`S_(u) = (2*T_(A)*K_(A)*A_(A))/(deltax)`
`S_(P) = -(2*K_(A)*A_(A))/(deltax)`
and hence above equation can be written as`a_(p)*T_(P) = a_(e)*T_(E) + S_(u)`
Where `a_(p) = a_(e) + a_(w) -S_(p)`
replacing `T_(P), T_(E)` with `T_(1), T_(2)`respectivley we have
`a_(p)*T_(1) = a_(e)*T_(2) + S_(u)`
same is being applied to node 5, representation of node 5 is as shown below,
Integrating flow heat equation from B to w `int_B^w del/(delx)(K*A*(delT)/(delx)) = 0`
we have `[K*A*(delT)/(delx)]_B^w = 0`
`[K*A*(delT)/(delx)]_w - [K*A*(delT)/(delx)]_B = 0`
in the above equation `(delT)/(delx)`is determined at A and e using interpolation scheme and same linear approach as that used above is applied
and hence we have
`(K_(w)*A_(w)*((T_(P)-T_(W))/(deltax_(w)))) - (K_(B)*A_(B)*(2*((T_(B)-T_(P))/(deltax))) = 0`
`T_(P)*[(2*K_(B)*A_(B))/(deltax)+(K_(w)*A_(w))/(deltax_(w))] = T_(B)*[(2*K_(B)*A_(B))/(deltax)]+T_(W)*[(K_(w)*A_(w))/(deltax_(w))]`
in the above equation, `a_(e) = 0`
`a_(w) = (K_(w)*A_(w))/(deltax_(w)`
`S_(p) = (2*K_(B)*A_(B))/(deltax)`
`S_(u) = ((2*K_(B)*A_(B))/(deltax))*T_(B)`
and hence above equation can be written as `a_(p)*T_(p) = a_(w) *T_(W) + S_(u)`
where `a_(p) = (a_(e)+a_(w)-S_(p))`
replacing `T_(P), T_(W)` with `T_(5), T_(4)` respectivley
we have `a_(p) *T_(5) = a_(w)*T_(4) + S_(u)`
all the `K_(e) = K_(w) = K_(A) = K_(B) = K = 1000 W/(m*k)`
and all `A_(e) = A_(w) = A_(A) =A_(B) = A =10*10^-3`, substituting these values in respective equations we have
for node 2, 3, 4 `a_(p) = 200, a_(w) = 100`
and hence at node 2 we have `200*T_(2) = 100*T_(1) + 100*T_(3)`
at node 3 we have `200*t_(3) = 100*T_(2) + 100*T_(4)`
and at node 4 we have `200*T_(4) = 100*T_(3) +100*T_(5)`
At boundary nodes we have `a_(e) = 100, S_(p) = 200, a_(w) = 100`
hence at node 1 we have `300*T_(1) = 100*T_(2) + 200*T_(A)`
and at node 5 we have `300*T_(5) = T_(4) * 100 + 200*T_(B)`
Now arranging all the linear equations from node 1 to node 5 we have
`300*T_(1) - 100*T_(2) = 200*T_(A)`
`200*T_(2) - 100*T_(1)-100*T_(3) = 0`
`200*T_(3) - 100*T_(2) - 100*T_(4) = 0`
`200*T_(4) - 100*T_(3)-100*T_(5) = 0`
`300*T_(5) - 100*T_(4) = 200*T_(B)`
expressing above linear equations in matrix form and hence solving the equation we will get the temperature distribution across the rod