A simple and high accurate finite volume scheme for diffusion equations

Abstract

Most of numerical methods for diffusion equations, refer to vertex unknowns directly or indirectly, and their accuracy is ultimately determined by the approximation to vertex unknowns. Based on the “twin-fitting” method, a simple and high accurate treatment for the vertex unknowns is developed and is applied to the nine-point scheme for diffusion problem. Numerical experiments show that the resulting nine-point scheme is high accurate for diffusion problems with discontinuous diffusion coefficients on distorted meshes.

Keywords

Discontinuous diffusion coefficients distorted meshes the “twin-fitting” method the treatment for the vertex unknowns

1 Introduction

Radiation hydrodynamics has been widely applied to energy, national defence and evolvement of celestial bodies, etc. When radiation hydrodynamics is simulated numerically, diffusion calculation usually couples with hydrodynamic calculation and naturally inherits large deformation meshes from the fluid flow simulation. The discontinuous diffusion coefficients are widely distributed over the computational domain, thus discretization of diffusion equations with discontinuous diffusion coefficients on distorted meshes becomes the main task of diffusion calculation.

If time discretization is not considered, the accuracy of discretizing diffusion equation is completely determined by that of numerical fluxes on cell edges. Many works are focused on constructing numerical fluxes and develop various methods. Based on the variational principle, there are related methods such as Kershaw scheme [1], Morel scheme [6] and Support-operator methods [8]. According to interpolation technique, related schemes include the difference scheme on polygonal meshes [5], linearity preserving schemes [7], 21-point scheme [11], Hermeline schemes [3], multi-point flux approximation methods [4] and the “twin-fitting” algorithms [9–10], etc. There are some algorithms [5, 13] relating to preserve the local maximum principle or solution positivity. These methods have made progress in improving their accuracy on distorted meshes, but most lose simplicity and some even increase the computational cost. For the discontinuous diffusion coefficients, most methods bring the low-order accuracy in the interface, but the “twin-fitting” algorithms [9–10] provide the global high-order accuracy.

The approximation to the normal component of fluxes on cell edges relates to the normal direction and the gradient of temperature, only the two cell-centered unknowns defined on the cells sharing an edge in general cannot give high accurate approximation on distorted meshes, and the vertex unknowns must be introduced in the construction of schemes such as scheme [2]. These schemes employ a linear combination of the neighboring cell-centered unknowns to calculate the vertex unknowns, and their accuracy is ultimately determined by the approximation for the vertex unknowns, which has been confirmed by many numerical results and will be further validated in this paper.

There are many works dedicated to improving the accuracy of approximation for the vertex unknowns. The simplest way is to use the arithmetic average weights [2], and another simple way is to make the weights direct proportion to the diffusion coefficients and inverse proportion to the distance between the vertex and its neighboring cell centers. The two methods are often used in practical computation due to their simplicity and fairly good accuracy on smooth meshes, but they lead to poor accuracy on distorted meshes. Recently, there are two kinds of treatments with fairly good accuracy, Taylor method [12] and linearity preserving methods [7]. The method [12] uses Taylor expansion and the requirement of second-order approximation to the vertex unknowns to calculate the weights. The methods [7] apply the linearity preserving criterion to obtain the weights. Numerical results show that the schemes using the methods [7, 12] can attain good accuracy on distorted meshes, but the accuracy lessens rapidly when the mesh is gradually distorted. Based on the “twin-fitting” method, a treatment for the vertex unknowns [10] is established, adopts different linear functions to approximate the temperature on sub-cells sharing a vertex, and shows the same numerical effect as the methods [7, 12].

A new treatment for the vertex unknowns is developed in this paper by preserving the continuity of temperature validly and by using the “twin-fitting” method. The difference between the new treatment and the treatment [10] lies in the choice of sub-cells sharing a vertex such that the new treatment is simple and high accurate.

The rest of the paper is organized as follows. In section 2, the scheme [2] is rederived by the “twin-fitting” method. In section 3, a new treatment for the vertex unknowns on general quadrilateral meshes is given. Numerical experiments in section 4 will demonstrate that the resulting nine-point scheme with the new treatment is simpler and higher accurate than the one with the methods [2, 10], and some conclusions are given in the last section.

2 The finite volume scheme

Consider the diffusion problem as follows, ${\begin{matrix} u_{t} - \nabla \cdot (K (x, y) \nabla u) = R (x, y) & in Ω \subset R^{2}, t ⩾ 0 \\ α u + β K \frac{\partial u}{\partial n} = ψ & on \partial Ω \\ u (x, y, 0) = u_{0} (x, y) \end{matrix}$ (1)

where K ⩾ k₀ = const > 0, R, u₀, α, β, ψ are given functions, Ω a bounded polygonal domain with its boundary denoted by ∂Ω and n the outward unit vector normal to ∂Ω, u a certain material temperature, K the diffusion coefficient and F = - K ∇ u the flux. The continuity of temperature and the normal component continuity of flux are the intrinsic properties of diffusion equation. The initial condition is given by u₀ and as α ≠ 0, β = 0, Dirichlet boundary condition is given, as α = 0, β ≠ 0, Neumann boundary condition is derived, and as αβ ≠ 0, mixed boundary condition can be presented.

Assume Ω is partitioned into a number of nonoverlapping polygonal cells, i.e. the intersection of any two intersected cells is either a common edge or a common vertex. Let $γ = \vec{AB}$ denote a common edge shared by a left cell with center O_L and a right cell with center O_R, see Fig. 1.

Fig. 1

The local mesh of edge AB.

On the left cell to integrate the diffusion equation in (1), and set $u_{L} = \frac{1}{V_{L}} \int_{O_{L}} udxdy$ to be the average value of temperature, where V_L denotes the cell area and n_L the outward unit normal vector of the cell, then there is $\frac{\partial u_{L}}{\partial t} = - \frac{1}{V_{L}} \int_{\partial O_{L}} F \cdot n_{L} ds + R_{L} .$ (2)

If u_L is defined as the unknown on cell center, then the key of discretizing (2) is to approximate ∫_γF · n_ABds on some edge γ.

Two polynomials are employed to approximate the temperature on the two sides of γ for the reason that the gradient of temperature is discontinuous in the case of discontinuous diffusion coefficients, which is the motivation of the “twin-fitting” method. Concretely, p_L (x, y) and p_R (x, y) are composed by area coordinate basis functions on ΔO_LAB and ΔO_RBA respectively and satisfy the given values u_L, u_A, u_B and u_R, u_A, u_B, then are determined by the constraint $F_{AB}^{L} = F_{AB}^{R}$ . Here $F_{AB}^{L} = - K_{L} \nabla p_{L} \cdot R (\vec{AB}),$ $F_{AB}^{R} = - K_{R} \nabla p_{R} \cdot R (\vec{AB}),$

$F_{AB}^{L}$ is the numerical normal component of flux of cell O_L on γ, $F_{AB}^{R}$ is that of cell O_R on γ, and $R (\vec{AB})$ is the normal vector of γ with the length |AB|. So there is $F_{AB}^{L} = F_{AB}^{R} = \tilde{K} [u_{L} - u_{R} + \tilde{D} (u_{B} - u_{A})]$ (3)

where $\begin{matrix} \tilde{K} = \frac{K_{L} K_{R}}{S_{L} K_{R} + S_{R} K_{L}} | AB |^{2} \\ \tilde{D} = \frac{\vec{O_{L} O_{R}} \cdot \vec{AB}}{| AB |^{2}} \end{matrix}$ $\begin{matrix} S_{L} = \vec{O_{L} A} \cdot R (\vec{O_{L} B}), \\ S_{R} = \vec{O_{R} A} \cdot R (\vec{{BO}_{R}}) . \end{matrix}$

If γ is on Dirichlet boundary, then $F_{AB}^{L}$ can be directly solved.

If the time discretization adopts the implicit scheme, then the Equation (2) can be discretized as follows $\frac{u_{L}^{n + 1} - u_{L}^{n}}{Δ t} = - \frac{1}{V_{L}} {\tilde{F}}_{\partial O_{L}}^{n + 1} + R (O_{L}),$ (4) where Δt is the time step, $u_{L}^{n}$ is the value of u_L at the time nΔt, ${\tilde{F}}_{\partial O_{L}}$ denotes the numerical normal component of flux on edges of ∂O_L given by (3). The formula (4) shows that the scheme [2] on polygonal mesh can be rederived by the “twin-fitting” method.

3 The treatments for the vertex unknowns

In this section, a new treatment for the vertex unknowns on quadrilateral meshes is established and can be extended to polygonal mesh.

Consider the vertex shared by four cells, see Fig. 2, assume u₀ be the vertex unknown on point O (0, 0), it is treated as a linear combination of its surrounding four cell-centered unknowns u_i (1 ⩽ i ⩽ 4) respectively defined on the cell-centers O_i (s_i, r_i).

Fig. 2

The local mesh of the vertex.

In order to calculate the vertex unknown, the diffusion coefficients around the vertex are seen as the discontinuous ones when the temperature varies sharply, and four linear functions p_k (x, y) (k = 1, 2, 3, 4) are employed to approximate the temperature on the sub-cells ΔOO_k-1O_k. These functions can be determined by the given temperature values on local stencils and by adopting a minimum problem to preserve the normal component continuities of fluxes on $\vec{{OO}_{k}}$ , that is, $min_{u_{0}} \sum_{k = 1}^{4} S_{k}^{2}, S_{k} = (\nabla p_{k} - \nabla p_{k - 1}) \cdot R (\vec{{OO}_{k}}) .$

This minimum problem is in fact equal to solve the following formula $\frac{\partial}{\partial u_{0}} \sum_{k = 1}^{4} S_{k}^{2} = 0,$

Thus, the vertex unknown u₀ can be expressed as a linear combination of u_i (1 ⩽ i ⩽ 4) as follows, $\begin{matrix} u_{0} = - \frac{1}{\sum_{j = 1}^{4} d_{j 4}^{2}} [(d_{14} d_{11} + d_{24} d_{22} + d_{43} d_{44}) u_{1} \\ + (d_{14} d_{13} + d_{24} d_{21} + d_{34} d_{32}) u_{2} \\ + (d_{24} d_{23} + d_{34} d_{31} + d_{42} d_{44}) u_{3} \\ + (d_{14} d_{12} + d_{34} d_{33} + d_{41} d_{44}) u_{4}] \end{matrix}$ where $d_{k 1} = ct g ∠ O_{k} {OO}_{k - 1} + ct g ∠ O_{k} {OO}_{k + 1},$ $d_{k 2} = - ct g ∠ O_{k} {OO}_{k - 1} - ct g ∠ {OO}_{k} O_{k - 1},$ $d_{k 3} = - ct g ∠ O_{k} {OO}_{k + 1} - ct g ∠ {OO}_{k} O_{k + 1},$ $d_{k 4} = ct g ∠ {OO}_{k} O_{k - 1} + ct g ∠ {OO}_{k} O_{k + 1},$ and this new method is expressed as “tritwin”

4 Numerical experiments

NPS-tritwin denotes scheme (4) with the new treatment, NPS-tritwinv1 scheme (4) with the weights given in [10], and NPS scheme (4) with the arithmetic average weights [2]. Numerical experiments are used to test their performance on distorted meshes of domain Ω = [01] × [0, 1] including smooth, random and Kershaw meshes, see Fig. 3. The discrete L² error Eh and its order of convergence q are used to analyze the results, where h denotes number of cells along x-axis same as that along y-axis and is 12 in Fig. 3.

Fig. 3

Meshes(smooth, random(ex.1), random(ex.2), Kershaw).

Example 1 is to test problem (1) with K = 1 and u (x, y, t) = sin πx sin πy on meshes of Fig. 3.

Example 2 is to test problem (1) with

$\begin{matrix} K = {\begin{matrix} 13 & x \in [0, 0.5) \\ 1 & x \in (0.5, 1] \end{matrix}, u (x, y, t) \\ = {\begin{matrix} sin π x sin π y & x \in [0, 0.5] \\ sin 5 π x sin π y & x \in [0.5, 1] \end{matrix}, \end{matrix}$

on meshes of Fig. 3.

Example 3 is to test the problem in example 1 on meshes of Fig. 4.

Fig. 4

Kershaw meshes with gradual distortion.

Example 4 is to test the problem in example 2 on meshes of Fig. 4.

Table 1 is the result of example 1, and Table 2 the result of example 2. From the two tables, it is shown that the three methods are almost 2nd convergent on smooth meshes, and NPS-tritwinv1 and NPS-tritwin are at least 2nd convergent on random meshes, but NPS has lower order of convergence and even has no convergence. On Kershaw meshes, NPS cannot attain 1st order of convergence, while NPS-tritwinv1 and NPS-tritwin can attain at least 1.5 order of convergence.

Table 1

The result of ex. 1

	Scheme	E12	E24	E48	q
Smooth mesh	NPS	8.47E-03	2.49E-03	6.52E-04	1.93
	NPS-tritwinvl	4.61E-03	1.20E-03	3.03E-04	1.98
	NPS-tritwin	5.72E-03	1.53E-03	3.92E-04	1.96
Random mesh	NPS	7.85E-03	1.09E-02	1.09E-02	0
	NPS-tritwinvl	3.20E-03	8.15E-04	2.02E-04	2.02
	NPS-tritwin	3.25E-03	8.42E-04	2.09E-04	2.01
Kershaw mesh	NPS	7.57E-02	5.23E-02	3.13E-02	0.74
	NPS-tritwinvl	3.16E-02	1.54E-02	5.42E-03	1.51
	NPS-tritwin	3.18E-02	1.16E-02	3.38E-03	1.78

Table 2

The result of ex. 2

	Scheme	E12	E24	E48	q
Smooth mesh	NPS	7.85E-02	1.80E-02	4.41E-03	2.03
	NPS-tritwinvl	7.88E-02	1.80E-02	4.40E-03	2.03
	NPS-tritwin	8.09E-02	1.85E-02	4.52E-03	2.03
Random mesh	NPS	7.92E-02	1.60E-02	9.71E-03	0.72
	NPS-tritwinvl	8.40E-02	1.83E-02	4.03E-03	2.18
	NPS-tritwin	8.39E-02	1.84E-02	4.04E-03	2.18
Kershaw mesh	NPS	8.07E-02	4.18E-02	2.39E-02	0.81
	NPS-tritwinvl	7.42E-02	2.65E-02	8.94E-03	1.57
	NPS-tritwin	6.77E-02	2.13E-02	7.24E-03	1.56

In Figs. 5 and 6, ln(r) measures the distortion of mesh and q-l2 is the order of convergence. From the two figures, it is shown that, when Kershaw meshes are gradually skewed, the convergent order of NPS-tritwin is at least 1.4, that of NPS-tritwinv1 is at least 1.2, and that of NPS lessens rapidly and even attains 0.

Fig. 5

Convergent order of ex. 3.

Fig. 6

Convergent order of ex. 4.

Thus it is obvious that the new treatment is superior to other treatments.

5 Conclusions

A simple and high accurate treatment for the vertex unknowns is established on general quadrilateral meshes by the “twin-fitting” method. Numerical experiments show that this treatment is superior to other treatments due to simplicity and high accuracy, and can adapt well to diffusion problems with discontinuous diffusion coefficients on distorted meshes.

References

Kershaw

D.S.

, Differencing of the diffusion equation in Lagrangian hydrodynamic codes, Journal of Computational Physics39 (1991), 375–395.

D.Y.

and Chen

G.N.

, An introduction to difference methods of parabolic equations, Science Press (in Chinese), 1995.

Hermeline

, Approximation of 2-D and 3-D diffusion operators with variable full tensor coefficients on arbitrary meshes, Computer Methods in Applied Mechanics and Engineering196 (2007), 2497–2526.

Aavatsmark

, An introduction to multipoint flux approximations for quadrilateral grids, Computational Geosciences6 (2002), 405–432.

Droniou

and Potier

C.L.

, Construction and convergence study of schemes preserving the elliptic local maximum principle, SIAM Journal on Numerical Analysis49(2) (2011), 459–490.

Morel

J.E.

, Dendy

J.E.

, Hall

M. L

and White

S.W.

, A cell-centered Lagrangian-mesh diffusion differencing scheme, Journal of Computational Physics103 (1992), 286–299.

J.M.

and Gao

Z.M.

, A nine-point scheme with explicit weights for diffusion equation on distorted meshes, Applied Numerical Mathematics61 (2011), 844–867.

Shashkov

M.J.

, Conservative finite-difference methods on general grids, CRC Press, Boca Raton, FL, 1996.

Song

S.H.

and Wang

S.H.

, High accuracy numerical method for heatconduction equation with discontinuous diffusion coefficient, Acta Mathematicae Applicatae Sinica (Chinese Series)34(2) (2011), 229–239.

10.

Song

S.H.

and Wang

S.H.

, A Simple and High Accurate Algorithm for Diffusion Problem, Acta Mathematicae Applicatae Sinica (Chinese Series)37(2) (2014), 367–378.

11.

Zhukov

V.T.

and Feodoritova

O.B.

, Difference schemes for the heat condunction equation based on local least-squares approximations, Preprint No, 97, IPMathem. Akad. Nauk S. S. S. R., Moscow, 1989.

12.

Liu

X.Z.

, Wang

R.L.

, Lin

and Yu

Y.L.

, Computational Methods and Program for Langrange Radiation Hydrodynamics Based On Unstructured Polygonal Meshes, Inertial Confinement Fusion Annual Report (in Chinese)2 (2014), 92–100.

13.

Sheng

Z.Q.

and Yuan

G.W.

, The finite volume scheme preserving extrenum principle for diffusion equations on polygonal meshes, Journal of Computational Physics230 (2011), 2588–2604.