Introduction

The purpose of this paper is to present a local support-operators diffusion discretization for arbitrary 3-D hexahedral meshes. We use the standard finite-element definition for hexahedra [2]. The method that we present is a generalization of a similar scheme for 2-D r - z quadrilateral meshes that was developed by Morel, Roberts, and Shashkov [1]. The diffusion equation that we seek to solve can be expressed in the following general form:

$${\frac{{\partial \phi}}{{\partial t}}} \mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$} \cdot$} \mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$} \phi$$ (1)

where t denotes the time variable, $\phi$ denotes a scalar function that we refer to as the intensity, D denotes the diffusion coefficient, and Q denotes the source or driving function. It is sometimes useful to express Eq. (1) in terms of a vector function, $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$}$ , that we refer to as the flux:

$$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \mbox{$\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{\nabla}$}$} \phi$$ (2)

We have taken the terms ``intensity'' and ``flux'' from the radiative transfer literature [3], but we have not explicitly considered the radiative diffusion equation because the subject of this paper relates to essentially any type of diffusion problem.

We define a cell-centered diffusion discretization scheme as one that numerically preserves the integral of Eq. (1) over each spatial cell. In particular, substituting from Eq. (2) into Eq. (1) and integrating that equation over a cell volume, we obtain:

$$\int_{V}^{} {\frac{{\partial \phi}}{{\partial t}}} \oint_{{\partial V}}^{} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{F}$} \mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$} \int_{V}^{}$$ (3)

where V denotes the cell volume, $\partial$V denotes the cell surface, and $\mbox{$\stackrel{^{\mathstrut}\smash{\longrightarrow}}{n}$}$ denotes the outward-directed unit surface normal. Note that we used the divergence theorem to convert the second integral in Eq. (3) from a volume integral to a surface integral. In physical terms, Eq. (3) generally represents a statement of particle or energy conservation over the cell. Thus we can simply state that cell-centered schemes (as we define them) are conservative over each mesh cell.

If one considers only non-orthogonal meshes with material discontinuities, existing vertex-centered diffusion discretizations are generally more advanced than cell-centered discretizations. This is primarily so because of the enormous success of Galerkin finite-element methods [2] and variants of those methods. Nonetheless, there are applications for which cell-centered schemes appear to yield superior accuracy relative to vertex-centered schemes. For instance, when coupling radiation diffusion calculations with cell-centered hydrodynamics calculations, a cell-centered diffusion scheme is highly desirable because it avoids certain difficulties associated with vertex-centered diffusion schemes [4]. Our new scheme has been developed with coupled radiation-diffusion/hydrodynamics applications in mind.

The discretization scheme that we have developed is cell-centered, but it has intensity unknowns at both cell centers and face centers. It can be applied on both structured and unstructured meshes consisting of combinations of arbitrary hexahedra and arbitrary degenerate hexahedra (i.e., wedges, pyramids, and tetrahedra). It yields second-order accurate solutions for the intensities on both smooth and non-smooth meshes even when material discontinuities are present, and it generates a sparse symmetric positive-definite coefficient matrix.

The literature relating to cell-centered diffusion discretization schemes for arbitrary hexahedra is not particularly extensive. One of the earliest relevant papers appeared about ten years ago. In particular, Rose developed a cell-centered hexahedral-mesh discretization scheme for the Laplacian operator.[5] The diffusion operator that we consider degenerates to the Laplacian operator when the diffusion coefficient is everywhere unity. Unlike our scheme, which has only the normal component of the current on each cell face, Rose's scheme has three components of the flux on each cell. Furthermore, the flux is continuous across each cell face in Rose's scheme, whereas only the normal component of the flux is continuous in our scheme. A central aspect of Rose's method is the preservation of an integral expression that is referred to as an energy principle. Our method is actually based upon the preservation of an integral identity. The energy principle used by Rose is not the same as the integral identity that we use, but they are related. In particular, the principle used by Rose can be derived from the diffusion equation together with the integral identity that we use. Rose presented a proof that his hexahedral-mesh method converges with second-order accuracy, but he provided computational results only for a 1 - D version of his method. Arbogast, et al., [6] have recently developed a cell-centered expanded mixed finite-element method for solving the tensor diffusion equation on general meshes (including hexahedral meshes.) Their method has only cell-center intensity unknowns if both the mesh and the diffusion tensor are smooth, but additional face-center intensities are required wherever the mesh or the diffusion tensor is non-smooth. The coefficient matrix generated by their method is always symmetric positive definite (SPD). The method of Arbogast, et al., actually shares some of the best properties of the standard mixed finite-element method and the hybrid mixed finite-element method. Standard mixed finite-element diffusion methods have only cell-center intensities, but this is achieved at the cost of solving a computationally expensive saddle-point linear system. The saddle-point system can be avoided by using the hybrid mixed finite-element approach, which generates a symmetric positive-definite coefficient matrix at the expense of additional face-center unknowns. The method of Arbogast, et al., yields an SPD coefficient matrix like the mixed hybrid method but can sometimes require far fewer unknowns. Although they proved several convergence theorems for their hexahedral-mesh method, Arbogast, et al., provided computational results only for a 2-D version of their method.

Our local support-operators method is similar to hybrid mixed finite-element methods in that it is cell-centered, it has both cell-center and cell-face intensities, and it produces a coefficient matrix that is symmetric positive-definite. However, our scheme is fundamentally a finite-difference technique since basis functions never appear in our formalism. The similarities between our method and mixed hybrid finite element methods suggest that there may be a connection between them, but we cannot presently demonstrate such a connection.

To summarize, the following combination of characteristics appear to be unique to our support-operators diffusion discretization scheme:

• It is a cell-centered discretization for arbitrary hexahedral meshes.
• It gives second-order convergence of the intensity on both smooth and non-smooth meshes both with and without material discontinuities.
• It generates a sparse SPD coefficient matrix.
• It is equivalent to the standard 7-point cell-center diffusion discretization scheme when the mesh is orthogonal.

The remainder of this paper is organized as follows. We first explain the central theme of our local support-operators method, and apply it to an arbitrary hexahedral mesh in Cartesian geometry. We next describe an approximate version of our scheme that we use as a preconditioner in conjunction with a conjugate-gradient solution technique [7]. Finally, computational results are given, followed by a summary and recommendations for future work.