Finite Volume Method¶
- So far we have discussed PDEs and ODEs using the finite difference (FD) method.
- Set a grid of discrete points.
* * * * *
- Approximate the derivatives in the ODE/PDE at those points.
- Solve for function values there.
- Note, in FD, the starting point is the PDE.
Finite Volume (FV) method¶
- Instead of a grid of points, make a grid of volumes.
-------------------------------
| | | | | |
-------------------------------
- Insted of using the PDE (differential conservation law), use the integral conservation law.
- This is a standard balance on control volumes:
- (accumulation) = (input) - (output) + (generation)
- This is a standard balance on control volumes:
3 equivalent approaches.¶
- Apply the above condition directly.
- Integrate the PDE over a given control volume (CV).
- Start from scratch and derive the governing equation from the conservation law (using the Reynolds Transport Theorem.)
Approach 1.¶
- 2-D unsteady heat transfer.
- Energy balance
- 2-D grid
- Control volumes i, j
- Control faces: e, w, t, b.
- Assume properties are uniform throughout a given cell.
- Assume properties are uniform along a cell face.
- (accumulation) = (in) - (out) + (generation).
$$\frac{dE}{dt} = q_w\Delta y + q_b\Delta x - q_t\Delta x - q_e\Delta y + S\Delta x\Delta y.$$
- $E$ (=) J
- $q$ (=) J/m$^2$s
- $S$ (=) J/m$^3$s
- $dE/dt = \Delta x\Delta y \rho c_pdT/dt.$
$$\frac{dT_{i,j}}{dt} = -\frac{1}{\rho c_p}\left(\frac{q_{e,i,j}-q_{w,i,j}}{\Delta x_i}\right) - \frac{1}{\rho c_p}\left(\frac{q_{t,i,j}-q_{b,i,j}}{\Delta y_j}\right) + S_{i,j}.$$
- This equation is usually solved in this form, that is, in terms of fluxes.
- Note, $q_{e,i,j} = q_{w,i+1,j}$.
- Hence, energy is strictly conservered.
- This is an important property that is built-in to FV methods, but not necessarily for FD methods.
- This is especially true for nonuinform grids with variable properties.
- This is an important property that is built-in to FV methods, but not necessarily for FD methods.
Evaluate the fluxes.¶
$$\vec{q} = q_x\vec{i} + q_y\vec{j},$$ $$q_x = -k\frac{\partial T}{\partial x},$$ $$q_e = -k\left.\frac{\partial T}{\partial x}\right|_e = -k\left(\frac{T_{i+1,j}-T_{i,j}}{\Delta x}\right),$$ $$q_w = -k\left.\frac{\partial T}{\partial x}\right|_w = -k\left(\frac{T_{i,j}-T_{i-1,j}}{\Delta x}\right),$$ $$q_t = -k\left.\frac{\partial T}{\partial y}\right|_t = -k\left(\frac{T_{i,j+1}-T_{i,j}}{\Delta y}\right),$$ $$q_b = -k\left.\frac{\partial T}{\partial y}\right|_b = -k\left(\frac{T_{i,j}-T_{i,j-1}}{\Delta y}\right).$$
- If we substitute this into the governing equation, then on a uniform grid we obtain:
$$\frac{dT}{dt} = \frac{\alpha}{\Delta x^2}(T_{i-1,j}-2T_{i,j}+T_{i+1,j}) + \frac{\alpha}{\Delta y^2}(T_{i,j-1}-2T_{i,j}+T_{i,j+1}) + S_{i,j}.$$
- This is the same equation as the FD method (for a uniform grid).
- The the approach is different.
- Nonuniform grids are common.
- Again, we normally solve the equation in terms of fluxes, not with the fluxes substituted into the governing equation.
Approach 2¶
- Integrate the PDE over cell i, j
$$\frac{\partial T}{\partial t} = \alpha\frac{\partial^2T}{\partial x^2} + \alpha\frac{\partial^2T}{\partial y^2} + S,$$
$$ \iint\frac{\partial T}{\partial t}dxdy = \iint\alpha\frac{\partial^2T}{\partial x^2}dxdy + \iint\alpha\frac{\partial^2T}{\partial y^2}dxdy + \iint Sdxdy, $$
- $\partial/\partial t$ commutes with $\iint$.
- Assume $T$, $S$, $\alpha$ are uniform:
$$\iint\frac{\partial T}{\partial t}dxdy = \frac{\partial}{\partial t}\iint Tdxdy = \frac{\partial}{\partial t}\left(T\iint dxdy\right) = \frac{\partial T\Delta x\Delta y}{dt} = \Delta x\Delta y\frac{\partial T}{\partial t}.$$
- Apply similarly for other terms.
- Our equation now becomes
$$\Delta x\Delta y\frac{\partial T}{\partial t} = \alpha\Delta y\int\underbrace{\frac{\partial^2T}{\partial x^2}dx}_{ \int_{x_j}^{x_j+\Delta x_j}\frac{d}{dx}\left(\frac{\partial T}{\partial x}\right)dx \\= \left(\frac{\partial T}{\partial x}\right)_e-\left(\frac{\partial T}{\partial x}\right)_w}+ \alpha\Delta x\int\underbrace{\frac{\partial^2T}{\partial y^2}dy}_{ \int_{y_j}^{y_j+\Delta y_j}\frac{d}{dy}\left(\frac{\partial T}{\partial y}\right)dy \\= \left(\frac{\partial T}{\partial y}\right)_t-\left(\frac{\partial T}{\partial y}\right)_b}+ S\Delta x\Delta y.$$
$$\Delta x\Delta y\frac{\partial T}{\partial t} = \Delta y\alpha\left(\frac{\partial T}{\partial x}\right)_e - \Delta y\alpha\left(\frac{\partial T}{\partial x}\right)_w + \Delta x\alpha\left(\frac{\partial T}{\partial y}\right)_t - \Delta x\alpha\left(\frac{\partial T}{\partial y}\right)_t + S\Delta x\Delta y.$$
(Again, repeating the previous equation) $$\Delta x\Delta y\frac{\partial T}{\partial t} = \Delta y\alpha\left(\frac{\partial T}{\partial x}\right)_e - \Delta y\alpha\left(\frac{\partial T}{\partial x}\right)_w + \Delta x\alpha\left(\frac{\partial T}{\partial y}\right)_t - \Delta x\alpha\left(\frac{\partial T}{\partial y}\right)_t + S\Delta x\Delta y.$$
- Now, using the following, etc., we obtain the same equation as approach 1.
$$\alpha\left(\frac{\partial T}{\partial x}\right)_e = -\frac{1}{\rho c_p}q_e,$$
$$\frac{dT_{i,j}}{dt} = -\frac{1}{\rho c_p}\left(\frac{q_{e,i,j}-q_{w,i,j}}{\Delta x_i}\right) - \frac{1}{\rho c_p}\left(\frac{q_{t,i,j}-q_{b,i,j}}{\Delta y_j}\right) + S_{i,j}.$$
Approach 3¶
- Use the Reynolds Transport Theorem (RTT).
- The RTT is in the form: (generation) = (accumulation) + [(out) - (in)].
- RTT:
$$\frac{dB}{dt} = \frac{d}{dt}\int_{CV}\rho\beta d\mathcal{V} + \int_{CS}\rho\beta\vec{v}\cdot\vec{n}d\mathcal{A}.$$
- Here, $B$ is some extensive property (like total energy), and $\beta=B/m$ is the intensive property (like energy per unit mass).
- Let $B=E$ and $\beta = E/m=e$.
$$\underbrace{\frac{dE}{dt}}_{(1)} = \underbrace{\frac{d}{dt}\int_{CV}\rho e d\mathcal{V}}_{(2)} + \underbrace{\int_{CS}\rho e\vec{v}\cdot\vec{n}d\mathcal{A}}_{(3)}.$$
- In term (3), $\vec{v}=0$, no velocity, so term (3) = 0.
- In term (2),
$$\frac{d}{dt}\int_{CV}\rho ed\mathcal{V} = \frac{d(\rho e\Delta x\Delta y)}{dt} = \Delta x\Delta y \rho c_p\frac{dT}{dt}.$$
In term (1) we apply the conservation law, which is the first law of thermodynamics: $$\frac{dE}{dt} = \dot{Q} + \dot{W} + S\Delta x\Delta y.$$
- There is no work: no velocity, and constant $\rho$.
- For $\dot{Q}$, we have $$\dot{Q} = -\int_{CS}\vec{q}\cdot\vec{n}d\mathcal{A} = -\left[\int_b() + \int_t() + \int_e() + \int_w()\right] = -[q_e\Delta y - q_w\Delta y + q_t\Delta x - q_b\Delta x].$$
Now, substitute the terms back in (and flipping the equation order): $$\Delta x\Delta y \rho c_p\frac{dT}{dt} = -[q_e\Delta y - q_w\Delta y + q_t\Delta x - q_b\Delta x] + S\Delta x\Delta y.$$
Rearranging, this gives the same equation as in the other methods:
$$\frac{dT_{i,j}}{dt} = -\frac{1}{\rho c_p}\left(\frac{q_{e,i,j}-q_{w,i,j}}{\Delta x_i}\right) - \frac{1}{\rho c_p}\left(\frac{q_{t,i,j}-q_{b,i,j}}{\Delta y_j}\right) + S_{i,j}.$$
Practice A¶
- Faces are set midway between grid points:
-----------------------------------------------------
| | | | | |
* | * | * | * | *
| | | | | |
-----------------------------------------------------
- Better accuracy for fluxes (derivatives have symmetric central differences).
- Has half a control volume at the boundaries
- In 2-D, face values are not at the center of faces.
Practice B¶
- Grid points are midway between cell faces:
-------------------------------------------------------
| | | | | | |
| * | * | * | * | * | * |
| | | | | | |
-------------------------------------------------------
- Center value is a better representation of the value in the CV. This is better for evaluating source terms.
- Has a whole control volume at the boundaries