Numerical solution of elliptic partial differential equations in discontinuous domains using a higher order accurate finite difference scheme

In this paper, we employ a higher order accurate finite difference scheme [1] for numerically solving an elliptic equation with discontinuous coefficients and singular source terms. Here we consider the most general form of elliptic equation which admits discontinuities of the dependent variable, the flux, the convection coefficient and the source term. This new scheme is obtained by coupling a recently developed Higher Order Compact (HOC) methodology with special treatment for the points just next to the points of discontinuity. The scheme is proficient in handling jumps across the interface quite efficiently and the overall order of accuracy of the scheme is at least two. The scheme is used to solve two-dimensional (2D) problems in polar coordinates on a non-uniform space grid. The grid is constructed in such a way that the grid points cluster around the circular interface (points of discontinuity), where the points of discontinuity themselves are not nodes. By wisely choosing the intensity of clustering around the points of discontinuity, which can be easily implemented in actual programming, one can minimize the overall error in computing. Numerous numerical studies on a number of problems are conducted and results are compared with those obtained with immersed interface and other well-known methods in the literature. The numerical results demonstrate that the proposed method provides outstanding results on relatively coarser grids.

Illustrations explain the paper; they do not show study data.

AI-generated illustration

↗ Enlarge

Non-Uniform Polar Grid ConceptClose ×

AI-generated conceptual illustration (A circular two-dimensional mathematical domain with a concentric inner circular interface). It explains the idea and does not show study data.

Open full-size image ↗

Introduction

During the last two decades, Higher Order Compact (HOC) finite difference schemes [2,3] have implanted themselves as highly potential tools in the field of Computational fluid Dynamics to study various aspects of fluid dynamics. These schemes produce very accurate results when the governing equations are smooth throughout the domain, but the scheme is tested for very few problems [2] when there are some discontinuities in the governing equations. These types of problems having discontinuities in the governing equations are generally called Immersed Interface problems which frequently arise in the field of science and engineering. They include biochemical processing, solid mechanics, composite materials, porous media flow, heat transfer, multiphase flow, mining and many others. Consider a two dimensional computational domain Ω split into two (or more) separate sub domains Ω + and Ω -by a lower dimensional interface Γ such that Ω = Ω +  Ω - Γ. The 2D elliptic equation in variable  can be written as

With some specified boundary conditions. Here  is the Dirac Delta function and σ is the value of the point source at the interface Γ.,  and f are piecewise continuous, along with it f may have a delta function singularity. Two physical jump conditions involving the unknown variables and its first derivatives (i.e. flux here) are required to solve this problem numerically. These jump conditions across the interface can be defined as

Where  =  = . the normal derivative of is, is the local normal unit vector to the interface towards the Γ +region; ( , ) → Γ means approaching the interface from Γ + side and vice versa. C And  denotes the jump discontinuity in  and flux along the interface with specific strengths respectively.

In the present paper, the approach similar to Mittal [1,2] is used to develop a numerical method for solving the elliptic equations in the presence of interfaces where both the variable coefficients and the solution itself may be discontinuous. Different approaches are used to discretize the given differential equation at the regular and irregular points. This scheme is at least third order accurate at the regular grid points and exactly second order accurate at the irregular points. Non-uniform space grid is constructed in such a way that the grid points cluster around the point of discontinuity and the points of discontinuity (or the interface) is not a node. Grid points are so clustered that the interval formed by two points on either side of the discontinuity acts almost like one single point. A special scheme is proposed at those points which in conjunction with the regular points produce excellent results when compared with those obtained with other well-known methods .

Mathematical formulation and Discretization procedure

Without loss of generality, we consider the solution domain Ω to be the circle of radius centred at 0. It contains a circular interface Γ of fixed radius = 0.5. We have considered uniform grid spacing along θ-direction and non-uniform grid spacing along -direction in such a way that grid points are clustered around the circular interface. This type of grid has been generated by using the same transformation that has been used by Mittal et al. [1] . Such a typical non uniform grid in polar coordinates along with the circular interface that is used for our calculations is shown in Fig. 1. The grid points are categorized into two types: regular and irregular grid points. Any point which is within one step length of the interface (either in -or θ-direction) is termed as irregular point. All other points are termed as regular points. Consider the 2D elliptic problem (1.1) with some specific boundary conditions. Let us also assume that( , ) ∈ (Ω ) ∪ (Ω ). For the regular points, recently developed higher order compact formulation 6 (HOC) is used. Now, if we consider the grid points ( , ); for each (for m=0,1,2,…,M-1), we shall have a 1D like situation along r-direction and we can develop a 1D like analysis with grid points (for j=0,1,2,…,N-1). Constructing a polar grid has an advantage that it is symmetric in direction. Since the irregular grid points lie on two concentric circles, we can define 1 (2)

↗ Enlarge

Figure 1Close ×

Two dimensional 40x40 mesh in a non-uniform grid along with a magnified view of the grid points clustered around the points of

Open full-size image ↗

Original paper figureSource PDF · page 2 ↗

Where Similarly Taylor's series expansion of ( , )

Hence, for all m=0, 1, 2… M-1, Eq. ( 2) yields

Our mesh structure allows us to let 1 2 , 0 m m d d  (See Fig. 1) which leads to

3) Under similar assumptions, one can show that for all m=0, 1, 2… M-1

And so on, where ( , ), ( , ) are the partial derivatives of ( , ) with respect to r.

Numerical Examples

To validate our formulation, we have solved various interface problems. In this present paper, we consider a problem which has discontinuous coefficients as well as a singular source term which were also used by Zhong [4] to validate their method. In this example, computational domain is a circle centred at (0, 0) containing a circular interface Γ of fixed radius r=0. 5.

According to the exact solution, the jumps in the derivatives across the interface depend upon the parameter b. Therefore, different values of parameter b equal to 0.5 and 1.0 are taken to test the performance of the scheme. Table 2 presents the maximum error norm on successive finer grids for the two values of b. It clearly shows that an overall order of accuracy close to two is obtained for all the successive grids.

—	b = 0.5	—	—	—	—	—	b = 1	—	—	—
Grid	||EN||∞		p	(present)	Zhong [4]	||EN||∞		p	(present)	Zhong [4]
20 × 20	0.001214	1.9967			0.001810	0.000964	1.9940			0.001345
40 × 40	0.000304	2.0355			0.000441	0.000242	2.0289			0.000265
80 × 80	0.00074				0.000117	0.000059				0.000065

Table text extracted from the original PDF; superscript letters are shown inline. The original table image is below.

View original table image

↗ Enlarge

Table 2Close ×

Numerical results showing infinity norm of error (||EN||∞) and order of accuracy (p).

Open full-size image ↗

Original paper tableSource PDF · page 3 ↗

Conclusion

In this paper, we implement an already developed HOC scheme on non-uniform grids for 2D elliptic equations with discontinuous coefficients and singular source term. To validate this formulation, we have considered a problem having discontinuous coefficients as well as a singular source term. Reasonably accurate results are obtained which are comparable with existing numerical results. Since the interface considered in the problem is nothing but a circle, one need not resolve the jumps in other directions resulting in much less errors in computations.

1. Mittal HVR, Kalita JC, Ray RK. A class of finite difference schemes for interface problems with an HOC approach. International Journal of Numerical methods in Fluids. 2016; 82(9):567-606.
2. Mittal HVR. A Class of Higher Order Accurate Schemes for Fluid Interface Problems. PhD thesis, Indian Institute of Technology Mandi, INDIA, 2016.
3. Mittal HVR, Ray RK. A New Direction to Higher Order Compact Scheme for Solving Immersed Interface Problems. Mathematical Sciences International Research Journal, 2015; 4(1):51-54.
4. Zhong X. A new high-order immersed interface method for solving elliptic equations with imbedded interface of discontinuity. Journal of Computational Physics. 2007; 225(1):1066-1096.