Extension of Kershaw diffusion scheme to hexahedral meshes

Numerical solution of the time dependent diffusion equation on non-orthogonal meshes in two spatial dimensions is an essential feature of typical Lagrangian radiation hydrodynamics simulation codes used for designing inertial confinement fusion targets and modeling other high energy density plasmas. Electron thermal conduction is treated with a non-linear diffusion equation and radiative transfer is often modeled using multi-group flux-limited diffusion. For two dimensions, an initially orthogonal mesh can become non-orthogonal due to hydrodynamic flow irregularities and the diffusion equation must be numerically differenced on this non-orthogonal mesh. A novel approach to this problem was reported by Kershaw [1], where he used a variational method to derive the difference operator corresponding to the continuous diffusion operator on a non-orthogonal r–z mesh. Kershaw’s method leads to a nine-point differencing stencil and tractable positive definite matrix solutions for problems of practical significance. It reduces to the standard five-point differencing stencil for orthogonal meshes. For this reason Kershaw’s method is often used as a benchmark for comparison of more recent, higher order methods [2]. In this paper, the discretization scheme developed by Kershaw is extended to three dimensional non-uniform hexahedral x–y–z meshes. As three dimensional radiation hydrodynamics simulations become more commonplace, this three dimensional Kershaw scheme can be a viable approach to solving the diffusion equation. While higher order discretizations exist, the 3D Kershaw method has the benefit of using only zone centered unknowns with a local computational stencil making it relatively easy to implement in Lagrangian codes. Along these lines, the detailed difference equations for the resulting 19 point computational stencil are presented in this paper. While this scheme has limited accuracy, it can still serve as a benchmark for comparison of more elaborate higher order, but more expensive schemes. This extension shares many of the same properties as the original method. The resulting matrix has the benefit of being symmetric positive definite (SPD). However, it is not an ‘‘M-matrix” and therefore negative