Vorticity-stream function method and MAC algorithm are adopted to systemically compare the finite volume method (FVM) and finite difference method (FDM) in this paper. Two typical problems—lid-driven flow and natural convection flow in a square cavity—are taken as examples to compare and analyze the calculation performances of FVM and FDM with variant mesh densities, discrete forms, and treatments of boundary condition. It is indicated that FVM is superior to FDM from the perspective of accuracy, stability of convection term, robustness, and calculation efficiency. Particularly ,when the mesh is coarse and taken as 20×20, the results of FDM suffer severe oscillation and even lose physical meaning.
1. Introduction
In the numerical solution of flow and heat transfer problems, the concepts of conservative and nonconservative equations were firstly proposed in [1, 2] in the 1980s. It is noteworthy that, from the perspective of differential unit, the conservative and nonconservative equations are equivalent, which are both the mathematical expression of physical conservation law. Nevertheless, the numerical calculation is implemented on the calculation unit of finite size, for which the two forms of equations are of different characteristics. Practical calculation process shows that the influence of differences between conservative and nonconservative forms on accuracy, stability, and efficiency of numerical calculation is significant.
Conservative governing equation and the corresponding discrete form show some advantages; for example, in the control volume with finite size, only the conservative equation can guarantee that the conservation law of the problem studied is satisfied [3–5]. The result obtained by the conservative equation is of higher accuracy generally. In the calculation of flow problem involving shock wave, the obtained flow field is usually smooth and stable employing the conservation form of governing equation, while using the non-conservation equation might lead to unsatisfactory spatial oscillations in the upstream and downstream regions of the shock wave [6–8]. Moreover, when the conservation equation is used to a body-fitted coordinate system, the conservativeness of governing equation can still be satisfied [9, 10]. More related researches and applications are presented in the literature [11–15]. Hence, the conservation property of discretized equations is desirable in engineering calculation. Based on this, it is of great significance to investigate the influence from conservation property on the calculation performances and provide guidance for the better implementation of it. The present paper compares the calculation performances by analyzing the finite volume method which is conservative and finite difference method which is nonconservative.
Some researchers made comparisons of the two methods. Patankar [16] and West and Fuller [17] pointed out that the results obtained in the situation that grids numbers decrease and mass conservation is rigorously controlled lose its physical meaning. Leonard [18] made a comparison of the accuracy of truncation error of the convection term in FVM and FDM. By theoretical derivation, he indicated that the truncation error of FVM is smaller than that of FDM in second-order central difference and second-order upwind difference. Botte et al. [19] made comparison of the accuracy, mass conservation, and computation time. It was validated by examples that FDM cannot ensure mass conservation while FVM can. Meanwhile, the influence of boundary conditions on the accuracy of algorithm was taken into consideration.
Previous researchers mainly focused on the comparison of accuracy. There are few systemic studies on the stability, robustness, and computation efficiency. This paper compared the calculation performances of FVM and FDM based on previous research. Comprehensive comparison of accuracy, stability of convection term, robustness, and calculation efficiency was performed. The results obtained may enrich the study on the calculation performances of conservative and nonconservative methods, which are expected to provide some reference for the practical calculation.
2. Physical Problem
To make the study possess general significance, it is necessary to analyze the calculation performances of the two methods based on variant problems. Lid-driven flow and natural convection flow in a square cavity which are both typical problems in computational fluid mechanics and numerical heat transfer are solved in this paper to make a systemic analysis on the subject.
As Figure 1(a) shows, in the lid-driven flow the lid moves horizontally at a constant velocity utop while the other three boundaries keep still. The Reynolds number, Re (Reynolds number) adopted in this paper, is 1000.
Schematic of different physical problem in a square cavity.
Lid-Driven Flow
Natural Convection Flow
Figure 1(b) shows the schematic of natural convection flow in a square cavity. The left and right boundaries are the first boundary conditions, between which the temperature on the left boundary is higher. The other two boundaries are adiabatic ones. In calculation, Pr (Prandtl number) is taken as 0.71 and Ra (Rayleigh number) is taken as 10^{6}.
3. Numerical Method
Considering that the calculation approach is influential to the calculation, in order to compare the difference of calculation performances between FVM and FDM, this paper adopts both primitive variable method and nonprimitive variable method to solve the governing equations. For primitive variable method, MAC algorithm is adopted, and for nonprimitive variable method, the generally used vorticity-stream function algorithm is adopted. The conservation property generally refers to the convection term, and thus only the difference of convective term between conservative and nonconservative forms is compared. MAC method is characterized by the direct solving of velocity and pressure. The key point of this method is to obtain a nondivergent velocity field in every time step or iteration step. That is, the velocity field of every time layer or iteration layer should satisfy the continuity equation and pressure, and the velocity is decoupled by solving Poisson equation which is relevant to pressure.
3.1. Governing Equation
In primitive variable method, pressure and velocity are directly used as variables to solve the equations. As to the incompressible fluid, the dimensionless unsteady-state conservative and nonconservative equations can be written as respectively,
(1)∂Φ∂τ+∂(UΦ)∂X+∂(VΦ)∂Y=Γϕ(∂2Φ∂X2+∂2Φ∂Y2)+Sϕ,(2)∂Φ∂τ+U∂Φ∂X+V∂Φ∂Y=Γϕ(∂2Φ∂X2+∂2Φ∂Y2)+Sϕ,
where the values of Φ, Γϕ, and Sϕ for lid-driven flow and natural convection flow are shown in Table 1 and Table 2, respectively. From (1) and (2) we can see that the convection term of conservative equation is given in the form of divergence while that in nonconservative equation is not, where the nondimensional excessive temperature can be described as Θ=(T-Tl)/(Th-Tl) and U, V are defined as nondimensional velocities in separate directions (X direction or Y direction). Governing equations of conservative and nonconservative equations of primitive variable method are given before. In the calculation process, MAC algorithm is adopted to obtain the result in steady state by unsteady calculation over enough time. The result of steady state is analyzed as follows from the angle of primitive variable method. Taking the vorticity-stream function method as example, the governing equations of conservative form and nonconservative form can be written as respectively,
(3)aϕ(∂(UΦ)∂X+∂(VΦ)∂Y)=Γϕ(∂2Φ∂X2+∂2Φ∂Y2)+Sϕ,(4)aϕ(U∂Φ∂X+V∂Φ∂Y)=Γϕ(∂2Φ∂X2+∂2Φ∂Y2)+Sϕ.
As to the vorticity-stream function method, only the lid-driven flow is studied in this paper. Values of Φ, aϕ, Γϕ, and Sϕ are shown in Table 3. where “ω” stands for vorticity and “ψ” stands for stream function.
Parameters of equations for MAC of lid-driven flow.
Name of variables
Φ
Γϕ
Sϕ
X direction
U
1Re
-∂P∂X
Y direction
V
1Re
-∂P∂Y
Parameters of equations for MAC of natural convection flow.
Name of variables
Φ
Γϕ
Sϕ
X direction
U
1
-∂P∂X
Y direction
V
1
RaΘPr-∂P∂Y
Temperature
Θ
1Pr
0
Parameters of equations for vorticity-stream function algorithm of lid-driven flow.
Name of variables
Φ
aϕ
Γϕ
Sϕ
Vorticity equation
ω
1
1Re
0
Stream function equation
ψ
0
1
-ω
where “ω” stands for vorticity, and “ψ” stands for stream function.
3.2. Discretization of Equations
Uniform mesh is adopted in this paper. In order to compare the calculation performances of FVM and FDM with variant discretized schemes, central difference and second-order upwind scheme are adopted. As to FVM, taking e interface on X direction as example, the expressions can be written as,
(5)Φe+=Φe-=ΦP+ΦE2,(6a)Φe+=1.5ΦP-0.5ΦW,Ue≥0(6b)Φe-=1.5ΦE-0.5ΦEE,Ue<0,
where “e” stands for the east interface between two nodes; “E” stands for east node and “W” means west node; and “EE” can be defined as the east node of “E” node.
As to FDM, taking node P as example, the expressions can be written as
(7)∂Φ∂X|P+=∂Φ∂X|P-=ΦE-2ΦP+ΦW2ΔX,(8a)∂Φ∂X|P+=3ΦP-4ΦW+ΦWW2ΔX,UP≥0(8b)∂Φ∂X|P-=-3ΦP+4ΦE-ΦEE2ΔX,UP<0.
3.3. Treatment of the Boundary Condition of Vorticity
Since the influence from boundary condition of vorticity on the solution of equations is significant, in this section, the calculation performances of FVM and FDM with variant treatments of the boundary condition for vorticity are studied.
Three boundary conditions of vorticity are compared and studied. In the boundary as shown in Figure 2, their expressions, Thom equations [20], Woods equation [21], and Jenson equation [22], as well as corresponding accuracy are shown in (9)–(11);
Schematic of the treatment for boundary conditions of vorticity-stream function method.
When MAC is adopted, momentum and energy equations are both solved in explicit scheme while in vorticity-stream function method, the convection term is treated in implicit scheme by deferred correction, which can guarantee the main-diagonal domination of the solving matrix and thus guarantee the calculation stability. As to FVM, taking Φe as example, the expressions of deferred correction are as follows:
(12)whenUe≥0Φe=ΦP+(Φe+-ΦP)*,whenUe<0Φe=ΦE+(Φe--ΦE)*,
where “*” in these equations stands for the former calculation step.
As to FDM, similarly to the treatment of Φe in FVM, the deferred correction for node P can be written as the sum of first-order upwind gradient and correction gradient which is shown as follows:
(13)whenUP≥0UP∂Φ∂X|P+=UP(ΦP-ΦWΔX)+[UP(∂Φ∂X|P+-ΦP-ΦWΔX)]*,(14)whenUP<0UP∂Φ∂X|P-=UP(ΦE-ΦPΔX)+[UP(∂Φ∂X|P--ΦE-ΦPΔX)]*.
The first terms in the right side of (12)–(14) are all first-order upwind term. The second terms are incorporated into the source term as the form of correction term. The discretized equation is solved by G-S solver with underrelaxation, of which the discrete expression is shown as (15) and the discretized equations of FVM and FDM are shown in (16a), (16b), (17a), and (17b), respectively:
(15)aPΦP=aEΦE+aWΦW+aNΦN+aSΦS+b,(16a)aE=max(-Ue,0)ΔX+1ReΔX2,aW=max(Uw,0)ΔX+1ReΔX2,aN=max(-Vn,0)ΔY+1ReΔY2,aS=max(Vs,0)ΔY+1ReΔY2,aP=max(Ue,0)ΔX+max(-Uw,0)ΔX+max(Vn,0)ΔY+max(-Vs,0)ΔY+2ReΔX2+2ReΔY2,(16b)b=-max(Ue,0)[Φe+-ΦPΔX]+max(-Ue,0)[Φe--ΦEΔX]+max(Uw,0)[Φw+-ΦWΔX]-max(-Uw,0)[Φw--ΦPΔX]-max(Vn,0)[Φn+-ΦPΔY]+max(-Vn,0)[Φn--ΦNΔY]+max(Vs,0)[Φs+-ΦSΔY]-max(-Vs,0)[Φs--ΦPΔY],
where aN, aP, aE, aW, and aS are coefficients in discretized equations and b is the source term in discretized equations.
4. Results Comparison and Analysis
In this section, the accuracy, stability of convection term, robustness and calculation efficiency of variant algorithms, schemes, and treatments of boundary condition which have been discussed previously are systemically compared and analyzed based on the calculation results of lid-driven flow and natural convection flow in a square cavity.
4.1. Accuracy
Firstly, the accuracy of FVM and FDM is evaluated by comparing bench mark solution with the results obtained from variant mesh density, discrete schemes, and boundary conditions. Through precision analysis, FVM and FDM are, respectively, adopted with the Taylor expansion commenced at midpoint (P node), and their respective truncation error (TE) is just as shown Table 4 [19]. The truncation errors of FVM are obviously less than the FDM’s.
Truncation error of FVM and FDM governing equations for vorticity-stream function algorithm of lid driven flow in a square cavity.
Convection item
Diffusion item
FVM (format)
ΦeΔX-ΦwΔX
1ΔX∂Φ∂x|e-1ΔX∂Φ∂x|w
FVM (TE)
-ΔX28ΦP′′′-ΔX4128ΦP(5)-⋯
-ΔX224ΦP(4)-13ΔX45760ΦP(6)-⋯
FDM (format)
∂Φ∂X|P
∂2Φ∂X2|P
FDM (TE)
-ΔX26ΦP′′′-ΔX4120ΦP(5)-⋯
-ΔX212ΦP(4)-ΔX4360ΦP(6)-⋯
4.1.1. Influence of Mesh Density
In [4], it is indicated that the conservative property can affect the accuracy while the nonconservative FDM loses its conservative property in the situation when the grid is too small. To make further analysis, the calculation results obtained by FVM and FDM with variant mesh densities are compared. The results of the velocity along the horizontal central axis V and velocity along the vertical central axis U obtained by FDM and FVM with central difference scheme of MAC algorithm are compared and illustrated in Figure 3.
Velocity on the central axis of lid-driven flow obtained by MAC with variant mesh densities.
80×80
40×40
20×20
From Figure 3 we can see that, when the mesh density is taken as 80 × 80, the results obtained by FVM and FDM both agree well with the benchmark solution. When the mesh density is reduced to 40 × 40, the results obtained by FVM are closer to the benchmark solution and if the mesh density is further reduced to 20 × 20, the results obtained by FDM differ from the benchmark solution greatly, in which FVM shows its advantage in accuracy. This phenomenon is in good agreement with the theoretical analysis in [18]; that is, the conservative property of FVM better ensures the accuracy of calculation. Therefore, whether the conservative property can be realized is of much significance to the accuracy of the obtained result.
4.1.2. Influence of Boundary Condition
In the following section, FVM and FDM with variant boundary conditions are compared from the perspective of accuracy by vorticity-stream function method. The results of the velocity along the horizontal central axis V and velocity along the vertical central axis U obtained by FDM and FVM with central difference scheme are compared and illustrated in Figure 4, and the influence of boundary condition of vorticity is also given. From Figure 4 we can see that no matter which treatment of boundary condition is adopted the accuracy of FVM is higher than that of FDM in vorticity-stream function method.
Velocity along the central axis obtained by second-order central difference scheme with two treatments of boundary condition when the mesh density is 40 × 40.
Thom
Woods
4.1.3. Influence of Discrete Scheme
The selection of discrete forms is also influential to the accuracy. Thus the results obtained with variant discrete forms are studied. Taking the vorticity-stream function with boundary condition treated by first-order Thom equation as example, the results of the velocity along the horizontal central axis V and velocity along the vertical central axis U obtained by FDM and FVM with variant discrete schemes are compared and illustrated in Figure 5, from which we can see that, whether second-order central difference or second-order upwind scheme is adopted, the accuracy of FVM is higher than that of FDM, the comparison between (a) and (b) indicates that second-order upwind difference possesses higher accuracy than second-order central difference.
Velocity along the central axis obtained by vorticity-stream function method for lid-driven flow with different schemes when the mesh density is 40 × 40.
Second-order central difference scheme
Second-order upwind difference scheme
4.1.4. Influence of Physical Problem
Finally, two methods in variant physical problems are compared from the perspective of accuracy based on the convection flow in a square cavity. Taking the result obtained by MAC algorithm as example, Figure 6 shows the comparison of results of the velocity along the horizontal central axis V and velocity along the vertical central axis U obtained by FDM and FVM. Temperature fields obtained by FVM and FDM are compared with the grid-independent solution (obtained with grid of 256 × 256) as Figure 7 shows.
Velocity along the central axis of convection flow in a square cavity with variant gird densities.
40 × 40
80 × 80
Comparison of temperature contours of natural convection flow in a square cavity with variant mesh densities (FVM in the left column and FDM in the right column).
40 × 40
80 × 80
From the comparison in Figure 6 and Figure 7 we can see that, based on the results obtained by MAC of natural convection flow in a square cavity, the accuracy of FVM is generally superior to that of FDM. Thus, it is indicated by comparison that, no matter which mesh density, scheme, and treatment of boundary condition is adopted, FVM shows better accuracy than FDM.
4.2. Stability of Convection Term
Stability of mathematical meaning can only ensure that the oscillation of solution is controlled within a range but cannot eliminate the absence of oscillation. The amplification of oscillation may result in divergence. Stability of convection term originates from varied difference schemes for convective item and is independent of the introducing of rounding error in the calculation. This stability is the key point to obtain a solution with physical meaning. It is concluded that all the unstable schemes will result in the oscillation of solution. Stability condition of the convection item (Mesh Peclet number) is affected by the following factors such as mesh density and boundary condition. Taking lid-driven flow as example, the stability of convection term of FVM and FDM with variant mesh densities treatments of boundary condition is compared.
4.2.1. Influence of Mesh Density
Figure 8 shows the comparison of stream lines between FVM and FDM with variant mesh densities. With the grid of 40 × 40, the two methods both can capture the structure of primary vortex and secondary vortex as Figure 8(a) shows. It is seen in Figure 8(b) that, when the mesh of 20 × 20 is adopted, concussion phenomenon appears, which is more serious in FDM, and the location of the vortex center obviously moves upward.
Comparison of stream lines in lid-driven flow with variant mesh densities (FVM in the left column and FDM in the right column).
40 × 40
20 × 20
4.2.2. Influence of Boundary Condition
Similarly to the accuracy analysis, the results obtained by variant treatments for vorticity of wall boundary are compared as Figure 9 shows. Three treatments, Thom, Woods, and Jenson, are analyzed. From Figure 9 we can see that, only when Jenson equation is adopted, the results obtained by FVM and FDM are relatively close, while when Thom and Woods are adopted, FVM shows better results than FDM. The results obtained by FDM suffer severe oscillation and even lack physical meaning. It is clear that, based on the lid-driven flow, stability of convection term of FVM is superior to that of FDM. The influence from the same boundary error (truncation error and rounding error) on difference solving methods is illustrated in Figure 9. With the treatment of boundary condition with the same truncation error or rounding error, the stability of FVM is better than that of FDM, especially in the situation with boundary treatment of (a) and (b) with relatively low precision truncation error; FDM is very unstable and very easy to slip into oscillation.
Stream functions of lid-driven flow under variant boundary treatments with the grid of 20 × 20 (FVM in the left column and FDM in the right column).
Thom
Woods
Jenson
4.3. Robustness
Robustness is a very critical index to evaluate the numerical approaches. In this section, the robustness of FVM and FDM is compared based on the lid-driven flow problem, which is solved by vorticity-stream function. Grids of 20 × 20, 40 × 40, and 80 × 80 are adopted, respectively, in the calculation. Relaxation factor is tentatively taken up from 0.2 to the maximum which can ensure the convergence of calculation. Iterations are correspondingly recorded as well.
It is shown in Figure 10 when the whether the central difference scheme or second-order upwind scheme is adopted, with the refining of calculation grids, the maximum relaxation factor (α) increases in both FVM and FDM, the robustness is enhanced, and the iterations (Niter) correspondingly decrease as well.
Comparison of the robustness of lid-driven flow with variant schemes.
Second-order central difference scheme
Second-order upwind scheme
As Figure 10(a) shows, in the situation that central difference scheme is adopted, with the grid of 20 × 20 and 40 × 40, the maximum relaxation factor of FVM is bigger than that of FVM while with the grid of 80 × 80, the relaxation factor can be taken as 1.0 in both methods. Similarly, as Figure 10(b) shows, generally speaking, when second-order upwind difference scheme is adopted, the maximum relaxation factor is bigger than that of FDM. Based on the result of lid-driven flow, the conservative FVM is superior to the nonconservative FDM from the perspective of numerical robustness.
In the numerical simulation of complex flow field, calculation efficiency is another critical index, on which this paper makes analysis. From Figure 10 we can see that, with the increasing of relaxation factor, the iterations of two methods both decrease generally. What needs to be notified is that the iterations of FVM are smaller than those of FDM.
4.4. Efficiency
In the numerical simulation of complex flow field, calculation efficiency is another critical index, on which this paper makes analysis.
As is shown in Table 5, although the calculation time required by every iteration step in FVM and FDM is almost the same, with the increasing of mesh numbers, calculation time required by FVM slightly decreases. Thus, total calculation time of FVM is much smaller than that of FDM. That is, FVM is much more efficient that FDM. Thus, it is validated that the efficiency of FVM is better than that of FDM in the same calculation condition.
Computation iterations and time of FVM and FDM in second-order central difference scheme for vorticity-stream function algorithm of lid driven flow.
Mesh
α
FVM
FDM
Iterations
Time(s)
Time of each iteration
Iterations
Time(s)
Time of each iteration
20
0.2
4950
73.8
0.015
2100
30.7
0.015
20
0.1
10500
154.4
0.015
4300
61.6
0.014
40
0.6
3250
151.9
0.047
4350
199.3
0.046
40
0.5
4450
203.5
0.046
5800
262.9
0.045
40
0.4
6300
296.3
0.047
7950
369.3
0.046
80
0.8
4950
822.3
0.166
5050
858.4
0.170
80
0.6
8500
1413.4
0.166
9000
1539.5
0.171
80
0.4
17250
2969.4
0.172
17750
3242.0
0.183
5. Conclusions
Based on vorticity-stream function and MAC algorithm, lid-driven flow and natural convection in a square cavity are calculated. The accuracy, stability of convective term, robustness of FVM and FDM with variant mesh densities, discrete forms, and treatments of boundary condition are compared. The following conclusions are drawn.
No matter which algorithm, discrete form, or treatment of boundary condition is adopted, the results obtained by conservative FVM are closer to benchmark solution than those obtained by FDM. The accuracy of FVM is higher than FDM, especially when the mesh density is relatively small.
The results of lid-driven flow indicate that the results obtained by FDM show more serious oscillation and the results obtained even lose physical meaning. Thus, the stability of convection term of FVM is superior than that of FDM.
As to the lid-driven flow, it is indicated that, when vorticity-stream function method is adopted, the relaxation factor of conservative method is bigger than that of the nonconservative one while the iterations are less, which shows that the robustness and calculation efficiency of conservative method are better than those of the nonconservative method.
In the lid-driven flow, with vorticity-stream function method, the calculation time required by every iteration step of FVM is slightly smaller than that of FDM. Together with the smaller total iteration steps, FVM is more efficient than FDM.
Acknowledgment
The study is supported by the National Science Foundation of China (nos. 51134006, 51176204, and 51276198).
AndersonJ. D.Jr.RoacheP. J.TaoW. Q.TaoW. Q.JaluriaY.TorranceK. E.LaxP. D.Weak solutions of nonlinear hyperbolic equations and their numerical computationLaxP.WendroffB.Systems of conservation lawsJamesonA.Transonic potential flow calculations using conservation formProceedings of the Computational fluid dynamics conference1975Hartford, Conn, USA148161VinokurM.Conservation equations of gasdynamics in curvilinear coordinate systemsBridgesT. J.Conservation laws in curvilinear coordinates: a short proof of Vinokur's theorem using differential formsDiPernaR. J.Decay of solutions of hyperbolic systems of conservation laws with a convex extensionBressanA.Hyperbolic systems of conservation lawsHankinR. K. S.The Euler equations for multiphase compressible flow in conservation form-Simulation of shock-bubble interactionsLaiC.BaltzerR. A.SchaffranekR. W.Conservation-form equations of unsteady open-channel flowGaoL. M.LiK. T.LiuB.J. SuGeneral integral conservation form of Navier-Stokes equations and numerical applicationPatankarS.WestA.FullerT.Influence of rib spacing in proton-exchange membrane electrode assembliesLeonardB. P.Comparison of truncation error of finite-difference and finite-volume formulations of convection termsBotteG. G.RitterJ. A.WhiteR. E.Comparison of finite difference and control volume methods for solving differential equationsThomA.The flow past cylinders at low speedsWoodsL. C.A note on the numerical solution of fourth order differential equations19545176182MR0066045JensonV. G.Viscous flow round a sphere at low Reynolds numbers (Re LTHEXA 40)