Multiphysics Learning & Networking

Introduction

One major difference between the FVM and the FDM is that the FVM is based on the integral form of the governing equations instead of the differential form. In the FVM, this discretization is conducted over each control volume, which can ensure that the flux entering the volume is identical to that leaving the volume.

Integral Form of the Governing Equation

The following equation is the general differential form of conservation problems. \[\underbrace{\frac{\partial u}{\partial t} }_{{\rm Accumulation}}\underbrace{+\nabla \cdot \left(uv\right)}_{{\rm Advection}}=\underbrace{-\nabla \cdot \left(\alpha \nabla u\right)}_{{\rm Diffusion\; }\left({\rm Conduction}\right)}+\underbrace{g}_{{\rm Source}}\] The diffusion term in the equation also contains other mechanisms with similar mathematical formulations such as dispersion. The FVM requires that the above equation is satisfied over the control volume $V_{P} $ around the centroid of the volume, leading to the following integral form of the governing equation: \[\begin{array}{l} {\int _{t}^{t+\Delta t}\left[\frac{\partial }{\partial t} \int _{V_{P} }^{}udV +\int _{V_{P} }\nabla \cdot \left(uv\right)dV \right]dt } \\ {=\int _{t}^{t+\Delta t}\left[\int _{V_{P} }\nabla \cdot \left(\alpha \nabla u\right)dV +\int _{V_{P} }gdV \right]dt } \end{array}\] or we have \[\frac{\partial }{\partial t} \int _{V_{P} }^{}udV +\int _{V_{P} }\nabla \cdot \left(uv\right)dV =\int _{V_{P} }\nabla \cdot \left(\alpha \nabla u\right)dV +\int _{V_{P} }gdV\] Applying the general form of the Green theorem, we reformat the above equation into \[\frac{\partial }{\partial t} \int _{V_{P} }^{}udV +\int _{\partial V}uv\cdot dS =\int _{\partial V}\alpha \nabla u\cdot dS +\int _{V_{P} }gdV ,\] where $S$ is the vectorial surface area of the control volume's boundary $\partial V$. Let us first assume the variation of the dependent variable is linear over a control volume, which can ensure a second-order accuracy: \[u\left(r\right)=u_{{\rm C}} +\left(r-r_{{\rm C}} \right)\cdot \left(\nabla u\right)_{{\rm C}} ,\] where the conserved quantity $u$ is a function of the position vector $r$ and the subscript ${\rm C}$ represents the centroid of the control volume under consideration. In the above equation, the volume integral in the first term on the left-hand side of the equation can be reformulated as \[\int _{V_{P} }^{}udV =\int _{V_{P} }\left[u_{{\rm C}} +\left(r-r_{{\rm C}} \right)\cdot \left(\nabla u\right)_{{\rm C}} \right] dV=u_{{\rm C}} V_{{\rm C}}\] The other terms can be transformed in a similar way, leading to the following equation: \[u_{{\rm C}} V_{{\rm C}} =\int _{t}^{t+\Delta t}\left[\sum \left(\nabla u\right)_{{\rm c}} \cdot \alpha _{{\rm c}} S +g_{{\rm C}} V_{{\rm C}} -\sum vu_{{\rm c}} \cdot S \right] dt,\] where the lower-case letter "c" in the subscripts indicates the centroid of the surface around the control volume. The summation is conducted over all the boundary segments, such as all the faces of a hexagon cell in an unstructured grid. To numerically deal with the integration, the above equation can be further changed into \[u_{{\rm C}} ^{t_{2} } =u_{{\rm C}} ^{t_{1} } +\frac{\delta t}{V_{{\rm C}} } \left[\sum \alpha _{{\rm c}} S\cdot \left(\nabla u\right)_{{\rm c}} +g_{{\rm C}} V_{{\rm C}} -\sum \alpha _{{\rm c}} S\cdot vu_{{\rm c}} ^{t_{1} } \right],\] where $t_{2} $ and $t_{1} $ are two consecutive time points for temporal integration, for which $t_{2} -t_{1} =\delta t$, and all the terms on the right-hand side of the equation are based on $u_{{\rm C}} ^{t_{1} } $ or the value(s) of $u_{{\rm C}} $ from even earlier time steps.

Example

Let us first revisit the general parabolic equation with the typical boundary conditions: \[\frac{\partial u}{\partial t} =\frac{\partial ^{2} u}{\partial x^{2} } ,\] \[\left. u\right|_{x=0} =a,\] \[\left. \frac{\partial u}{\partial x} \right|_{x=L} =\beta ,\] \[\left. u\right|_{t=0} =u^{0}\] In order to obtain an algebraic equation which is comparable to that in the FDM, we need to discretize the domain in the following way. The 1D domain with a length of $L$ is divided into $I-1$ portions with a uniform mesh size of $h$: there are two portions with a length of $h$/2 on the two ends (boundary) and $I-3$ portions in the interior of the domain. Therefore, we have $I=\frac{L}{h} +1$. The physical process to be simulated lasts a time of $Nk$. The general form of the integral equation including an accumulation and a diffusion term is as follows, \[u_{{\rm C}} ^{t_{2} } =u_{{\rm C}} ^{t_{1} } +\frac{\Delta t}{V_{{\rm C}} } \left[\sum \alpha _{{\rm c}} S\cdot \left(\nabla u\right)_{{\rm c}} \right]\] This is a 1D problem, so there is no surface area between different control volumes (or cells). Let us use a unit area instead, thus $\left|S\right|=1$. $S$ has a positive direction if it is on the right side of the cell and negative on the left. The centroid of the cell is the middle point of the cell. Accordingly, the above integral equation can be transformed into \[u^{n+1} =u^{n} +\frac{k}{h} \left[\left(\frac{\partial u}{\partial x} \right)_{{\rm right}} -\left(\frac{\partial u}{\partial x} \right)_{{\rm left}} \right]\] The subscript "${\rm C}$" is removed as we use the function value in the middle of a cell $u_{{\rm C}} $ to represent the $u$ of that cell. $u^{n+1} $ and $u^{n} $ are the values at the next and current time steps, respectively. $\left(\frac{\partial u}{\partial x} \right)_{{\rm right}} $ and $\left(\frac{\partial u}{\partial x} \right)_{{\rm left}} $ are the derivative on the right and left ends of the cell, respectively, whose signs in the equation are determined by the direction of $S$ and $\left(\nabla u\right)_{{\rm c}} $. The centroid of the boundary surface does not require special treatments as the boundary segments are two infinitesimally small points in this 1D domain.

Without considering the boundary conditions, the governing equation can be directly discretized as

Point 1: $u_{1}^{n+1} =u_{1}^{n} +\frac{k}{h} \left(\frac{u_{2}^{n} -u_{1}^{n} }{h} -\frac{u_{1}^{n} -u_{0}^{n} }{h} \right)$

Point 2: $u_{2}^{n+1} =u_{2}^{n} +\frac{k}{h} \left(\frac{u_{3}^{n} -u_{2}^{n} }{h} -\frac{u_{2}^{n} -u_{1}^{n} }{h} \right)$

$\dots $

Point $i$: $u_{i}^{n+1} =u_{i}^{n} +\frac{k}{h} \left(\frac{u_{i+1}^{n} -u_{i}^{n} }{h} -\frac{u_{i}^{n} -u_{i-1}^{n} }{h} \right)$

$ \dots $

Point $I$: $u_{I}^{n+1} =u_{i}^{n} +\frac{k}{h} \left(\frac{u_{I+1}^{n} -u_{I}^{n} }{h} -\frac{u_{I}^{n} -u_{I-1}^{n} }{h} \right)$

The boundary conditions are implemented by changing the discretized equations in the boundary cells, i.e., Point 1 and Point $I$, with the following equations:

$u_{0}^{n} =a$ , at Point 1

$\left. \left(\frac{\partial u}{\partial x} \right)\right|_{{\rm right}} =\beta $ , at Point $I$