Research Article | Open Access

Volume 2020 |Article ID 8161804 | https://doi.org/10.1155/2020/8161804

Xiaoliang Zhu, Yongbin Ge, "Adaptive ADI Numerical Analysis of 2D Quenching-Type Reaction: Diffusion Equation with Convection Term", Mathematical Problems in Engineering, vol. 2020, Article ID 8161804, 19 pages, 2020. https://doi.org/10.1155/2020/8161804

# Adaptive ADI Numerical Analysis of 2D Quenching-Type Reaction: Diffusion Equation with Convection Term

Accepted29 Apr 2020
Published15 Jun 2020

#### Abstract

An adaptive high-order difference solution about a 2D nonlinear degenerate singular reaction-diffusion equation with a convection term is initially proposed in the paper. After the first and the second central difference operator approximating the first-order and the second-order spatial derivative, respectively, the higher-order spatial derivatives are discretized by applying the Taylor series rule and the temporal derivative is discretized by using the Crank–Nicolson (CN) difference scheme. An alternating direction implicit (ADI) scheme with a nonuniform grid is built in this way. Meanwhile, accuracy analysis declares the second order in time and the fourth order in space under certain conditions. Sequentially, the high-order scheme is performed on an adaptive mesh to demonstrate quenching behaviors of the singular parabolic equation and analyse the influence of combustion chamber size on quenching. The paper displays rationally that the proposed scheme is practicable for solving the 2D quenching-type problem.

#### 1. Introduction

Nonlinear reaction-diffusion equation with a singular or near-singular source term has been widely applied in ion conduct polarization theory , computational fluid dynamics , electromagnetism , material research , ecology , thermology [6, 7], and so on. Different from the linear reaction-diffusion equation, it has a chance to produce a solution with singularity. Blow-up and quenching are considered as singularity of a solution. At a limited time, the former indicates that its solution tends to infinite, and the latter indicates that its temporal derivative tends to infinite, but its solution is restricted in a certain scope [8, 9]. Researchers often use the nonlinear singular degenerate reaction-diffusion model to simulate ignition and burning and blow-up of the fuel, which can describe burning transformation from a stable status to an unstable status in a combustion chamber and from which quenching may arise . To a large extent, it is to cognize features and change laws of quenching precisely that decides the design and optimization of the combustion chamber. We pay attention to the quenching problem of the degenerative singular reaction-diffusion equation with a convection term in this paper.

Let us recall the development of the singular reaction-diffusion equation of the quenching type. Since Kawarada concluded that the time differential quotient of the solution leans toward infinity when the solution approximates to one indefinitely for this equation , more and more scholars have focused on nonlinear singular problems. The academic system of researching the quenching phenomena of the nonlinear reaction-diffusion equation has been gradually established, which ranges from one dimension and two dimension to three dimension . More one-dimensional problems are investigated rather than higher-dimensional ones. However, there are still some quenching principles of the two-dimensional degenerate singular reaction-diffusion equations discovered [610, 1215]. Specially, they qualitatively depicted quenching or nonquenching under certain conditions, asymptotic status of quenching time, the spatial definition domain size, i.e., critical length, resulting in quenching behaviors, and so on. Sheng et al. have been engaged in research in this field. Especially, Padgett and Sheng investigated the degenerate stochastic Kawarada equations from 1D to 2D and 3D by combining the finite difference schemes with arbitrary meshes [11, 16, 17]. The low-accuracy concepts are exploited in the above resolutions. Ge et al. firstly rendered a high-order compact difference scheme based on an adaptive mesh to study the quenching principle of the one-dimensional singular degenerate reaction-diffusion equation in 2018 .

Let us recall the development of the reaction-diffusion equation with a convection term of the quenching type again. Compared to the quenching equation mentioned in the prior paragraph, less attention is paid to the convection-reaction-diffusion equation with quenching singularity, especially to the 2D and 3D problems . So far, it is difficult to find research on the 2D quenching problems. The authors in  discussed quenching phenomena under the conditions of different convection terms and singularity indexes: asymptotic status of quenching time, the spatial definition domain size, i.e., critical length, and so on, for the one-dimension singular equation. Sheng and Cheng studied quenching behaviors of the same type equation by using the fully adaptive difference scheme of low-order accuracy in . Sheng and Khliq investigated the relationship between the critical value of the degenerative term that was not a function of spatial points and both quenching time and convergence time for the 1D singular degenerate problem in . Zhou et al. theoretically depicted the critical length and the quenching spatial point for the convection-reaction-diffusion equation without the degenerate term . Research studies above belong to the low-order difference strategy. Recently, a high-order adaptive difference scheme is firstly used to study the influence of degenerate function, convection function, nonlinear source function, and spatial definition length on quench behaviors for the one-dimensional convection-reaction-diffusion equation, respectively . It provides the basis for adopting the high-precision algorithm to analyse the quenching states of the 2D convection-reaction-diffusion equation.

Depending on the previous research studies, we discover that there is an important significance for researching this two-dimensional convection-reaction-diffusion problem of the quenching type. It is well known that the finite element scheme and the finite difference scheme are popular numerical methods to solve the convection-reaction-diffusion equation . The improved schemes based on the finite element are utilized to resolve the two-dimensional convection-reaction-diffusion problems [24, 25]. In the past years, there have emerged many finite difference schemes for solving the two-dimensional convection-reaction-diffusion problem . Ge and Cao built up the transformation-free high-order compact difference algorithm proposed by Kalita et al.  on the multigrid to solve the 2D steady convection-diffusion equation . Huang proposed an ADI scheme to solve the 2D unsteady convection-reaction-diffusion equation, which is built by using the first and the second central difference operators to approximate the first-order and the second-order spatial derivative, respectively, and the Crank–Nicolson (CN) method to the discrete temporal derivative . Li et al. represented a compact finite difference scheme of unconditional stability and the fourth order for solving the groundwater pollutant problem based on the 2D nonlinear unsteady convection-diffusion equation . Of course, there exists other schemes for solving the 2D problem. A Fourier spectral discretization method was posted for explaining the high-dimensional fractional reaction-diffusion equations by using two separate mathematical techniques . Pindza and Owolabi developed the Fourier spectral method in space and the exponential integrator scheme in time to solve the 2D and 3D fractional reaction-diffusion equations . Nonlinear high-dimensional convection-diffusion problems were resolved by a novel mesh-less local Petrov–Galerkin method proposed in . Different from the above research studies, Owolabi and Atangana  studied the high-dimensional fractional reaction-diffusion systems of spatial interactions of three components’ species. The dynamical system was considered to be of asymptotical stability both in local and global style by analyzing the main equation theoretically. Furthermore, a theory certifying the existence and permanence of the three species was put forward.

In general, a high-order finite difference scheme has a greater advantage comparing with other methods. An adaptive ADI difference scheme based on the two-dimensional unsteady reaction-diffusion equation with convection function is used to analyze quenching phenomena of the 2D combustion model through serial cases and investigate relationship between burner shape and quenching features. Due to significant practical meaning of quenching behaviors, the important quenching information is beneficial to engineering applications such as fuel burning and industry manufacturing. The compact scheme with the adaptive algorithm is employed to compute a series of important information including quenching time and quenching location. There are five parts in the paper. Section 1 introduces the theme of this study. Section 2 describes carefully the proposed scheme of high-order accuracy. Section 3 proves adaptive mesh theory. Section 4 stimulates some typical numerical samples to explore and explain quenching problems. Section 5 draws the conclusion.

#### 2. Problem and Difference Discretization

##### 2.1. Semilinear Reaction-Diffusion Equation with Convection Term

This semilinear degenerate problem model involving two spatial dimensions is regarded as equation (1) with the boundary and initialization conditions of equation (2) and (3). The solution represents the temperature of the combustion chamber. refers to the smooth rectangle domain of the combustor, and is its boundary. and are the convection functions of and , and is the degeneracy function and . The singularity source is strictly increasing for with . The typical convection-reaction-diffusion equation of the quenching type is written as

Its boundary conditions are

Its initialization condition is

With the aid of the intermediate variables and , we replace and in equation (1) which can be defined as

For the convenience of expression, we use and instead of and . Then, the above formulation can be rearranged aswhere and Correspondingly, the boundary conditions of equation (5) are

The initialization condition of equation (5) is

In the physical application, we rely on equations (5) and (7) to compute and . Through observing a large number of values of and , we can capture the quenching moment, i.e., quenching occurs when is infinitely close to 1 and and becomes so huge that its value blows up.

##### 2.2. Nonuniform ADI Difference Scheme

According to Taylor formula , the high-accuracy ADI difference scheme of equation (5) is deduced below. The steady convection-diffusion equation derived from the 1D unsteady convection-diffusion equation in the direction is considered as where is a function of . By virtue of the Taylors’ series expansions, we get expressions of and , including the central difference operators of the first and the second spatial derivatives, discretized on the nonuniform mesh, respectively. The two expressions are substituted into the above steady convection-diffusion equation. So, a 1D steady convection-diffusion equation dispersed at point can be written aswhere

After discretizing and in the spatial direction, equation (8) is reorganized aswhere

Similarly, from the one-dimensional steady convection-diffusion equation in the direction: , and the equation in the direction can be derived aswhere

Adding the two hands of equation (10) and equation (12), we can get a 2D steady convection-diffusion equation:

We can mark the truncation error as and get

Subsequently, after substituting for in equation (14), the situation of point for equation (15) is considered as

The Crank–Nicolson method is used to discretize the temporal term of equation (16), and an equation of point can be formed aswhere is the time step and is the truncation error. When is added to the left side of equation (18) and is added to the right and the truncation error is deleted, a group of linear system is organized as

By virtue of the intermediate variable , equation (18) is reorganized as

After equation (19) is computed and arranged again, a linear system is formed as

In equation (21),

The elements of are

The element of iswhere

The element of iswhere

Similarly, from equation (20), another linear system is written aswherein which

It can be observed that each on the th time layer is obtained, according to on the prior one, by computing the tridiagonal linear equations twice. We adopt the catch-up method to solve the linear equations. We rely on equation (21) to get the medial variable , which represents the direction. Subsequently, we depend on and equation (28) to gain , which carries information in the direction and so on, until the solutions on the last time axis are worked out. The accuracy analysis of equation (18) is performed, and its truncation error is

Note that equations (19) and (20) have the second-order accuracy in time and the third- or fourth-order accuracy in space. The scheme of uniform owns spatial accuracy of the fourth order, whereas its nonuniform scheme owns spatial accuracy of the third order.

##### 3.1. Temporal Adaptive Mesh Method

Comparing with solution , the temporal derivative is more sensible. Therefore, like the 1D adaptive method, the 2D adaptive mesh method in the time direction is based on the arc-length monitor function of the temporal derivative function . We can extend the 1D adaptive temporal grid schemes in Reference  to the 2D case:

After in equation (32) is discretized through the central difference formulas to calculate the maximum ratio, the adaptive time step is easily obtained by equation (32) and in the prior time layer.

##### 3.2. Spatial Adaptive Mesh Method

Similarly, we can utilize the arc-length monitor function in the space direction based on and extend from the 1D adaptive spatial algorithm to the 2D problem by relying on Reference . Therefore, representing the monitor function in the direction and representing the monitor function in the direction can be written as

The two spatial monitor functions are positive in their definition area. and represent the set of 2D original spatial mesh points. After an adaptive improvement is complemented, the grid in the 2D computational area is updated as in the direction and in the direction. It is easily observed that the 1D adaptive grid scheme acts on in the direction and direction, respectively, to produce new grid points.

We take the equidistributing mesh in the direction as an example below. From , we know that the subareas of in intervals of are equal. Obviously, when , the scope is uniform. It can be expressed aswhere is the average arc-length. Equation (35) is equal to

After the second derivation of the aforementioned formula, the following equation can be obtained:where represents the iterative number, represents the spatial grid node set before iteration, and represents the space grid node set after iteration. Equation (37) can be solved by the difference scheme:which approximates to the former. Similarly, the semi-implicit scheme in the direction is defined as

The adaptive iterative procedure of the presented method is carried forward as follows. Firstly, the spatial monitor functions and can be obtained through equations (33) and (34) based on the grids and in the th time layer. Secondly, depending on equations (38) and (39), we can gain derived from in the direction and derived from in the direction. The new points substitute for the old points iteratively when and are reached for . is the terminal value of . Thirdly, from the value built on the last grid point sets and , is computed by using the method of area ratio. Thus, we can get in the th time layer depending on equation (18). Lastly, repeat the aforementioned produce till quenching occurs or the solution converges to a steady solution.

#### 4. Numerical Illustration

##### 4.1. Common Example

It is not a quenching case. It is just for demonstrating the accuracy and convergence of the proposed difference scheme. We have solved the following unsteady convection-diffusion equation with coefficient whose exact solution is offered:

The initial and boundary value conditions are derived from the exact solution. There are four metrics for the data test, which are maximum error, average error, CPU time, and convergence rate. We adopt these parameters to fulfill the two difference schemes: the proposed scheme based on the uniform mesh and the proposed scheme based on the nonuniform mesh. and represent the nonuniform grid, where point to intervals of the spatial direction and , whereas point to the mesh transformation coefficient and .

Table 1 offers these measuring data. It displays the accuracies and running time of our scheme (uniform and nonuniform). We select five spatial intervals: 16, 32, 64, 128, and 256 as test cases, which are given to each of the above two meshes (i.e., uniform and nonuniform). For the nonuniform mesh, the different values of impact these precision criteria. In the paper, takes 0.7. The maximum errors of the first scheme are higher than those of the second one, and the average errors of the latter are superior to those of the former. Furthermore, the larger the number of spatial intervals is, the higher the accuracies of the scheme based on the nonuniform grid are. Its accuracies arrive at 1.3364 × 10−7 in the maximum error and 6.4492 × 10−9 in the average error when . In terms of computing time of the two schemes whose measurement unit is millisecond, the first scheme is lower than the second one. Obviously, by comparing with the uniform mesh, the difference scheme based on the nonuniform mesh spends more time to run. Lastly, the conv. rate values are observed. The tendency of the conv. rate from the former scheme is similar to that from the latter scheme. We take the latter as the example. The four values are positive, which means the accuracy increases as rises when . The conv. rate value is maximum when transfers from 16 to 32, which means the accuracy promotes greatly in the stage. The latter three values are close, which means the accuracy promotes gently in the three stages.

 Difference schemes Max. error Aver. error CPU time Conv. rate The proposed scheme (uniform) 16 2.4859 × 10−1 1.7858 × 10−2 47 — 32 5.7032 × 10−2 1.9148 × 10−3 172 5.0033 64 5.5382 × 10−3 1.4565 × 10−4 702 3.7166 128 3.1469 × 10−4 9.6850 × 10−6 2948 3.9106 256 1.9710 × 10−5 6.1743 × 10−7 11794 3.9714 The proposed scheme (nonuniform) 16 8.1898 × 10−3 8.7725 × 10−4 358 — 32 2.9724 × 10−4 2.9996 × 10−5 2855 7.9861 64 1.5573 × 10−5 1.6521 × 10−6 14914 4.1824 128 9.3004 × 10−7 1.0224 × 10−7 35006 4.0143 256 1.3364 × 10−7 6.4492 × 10−9 165969 3.9867
##### 4.2. Quenching Case without Convection Term

This section will describe a quenching status without the convection term, i.e., and . The mathematical model can be expressed as

Its initial and boundary conditions are just as in equations (6) and (7).

After equations (19) and (20) are used to approximate equation (41), we can get some quenching behaviors. For the 2D nonlinear degeneration singularity equation without the convection term, we demonstrate the quenching case of that has been described in [13, 14]. In Table 2 and Figure 1, , , and which are the same for the three approaches, but the initial space step and the initial temporal step are different among them, where we take and . In Table 2, the five parameters represent the same senses, in which represents the point of quenching location, represents the quenching time, represents the maximum temperature value immediately before quenching happens, and means the maximum variation of temperature with respect to time immediately before quenching happens. From Table 2, we can see that the quenching location of the referred scheme is identical with that of the scheme in , and and of the first one are close to those of the last two. Obviously, is different among the various methods. From Figure 1, it is found that and at the location and the quenching time .

 Schemes The proposed scheme 0.5958348109999995 0.500 0.500 0.998562875689732 689.2724781649539 The scheme in  0.58712499993751 — — 0.990432 148.887767 The scheme in  0.587554 0.500 0.500 0.999263 1249.917563
##### 4.3. Quenching Case with Convection Term

The standard 2D degenerate semilinear quenching problem with the initial boundary condition lists in equations (5)–(7) is extended from the equation in one-space dimension . The first convection term can be written as two functions and , where and take the constants. The four groups of variables: , and , and , and in equations (5)–(7) affect the quenching characteristics of the 2D problem, whose meaning has been stated in Section 2. There follows a special case for . The initial time step is considered as , and the initial space step is considered as , in which the -direction is identical with the -direction. The four groups of parameters are configured as which is defined as Case 7.0. This will lead to a quenching phenomenon. Quenching location is ; the quenching time is 4.487776579638927; is 0.9992409498317957; is 1539.153786381783. When is 4.486155872847072, the temporal step starts to become smaller and smaller until quenching occurs. We choose six rows of representative data given in Table 3 which poins to features relative to , , , , and before quenching occurs of Case 7.0. and refer to maximal and maximal in each time layer before quenching occurs; and refer to the coordinate point of synchronously producing and in each time layer before quenching occurs; refers to the moment at which each time layer is located before quenching occurs. first appears at the point when is 4.487088240486234 which is the first row in Table 3 and maintains at this point until quenching occurs. When is which is the immediate time before quenching happens, we can obtain important characteristics. When is 4.487777138788559, blows up and reaches 1.051076326841211 × 109.

 4.487088240486234 0.03636363636363636 0.03636363636363636 0.9794485698289615 52.18443133976322 4.487489880800603 0.03636363636363636 0.03636363636363636 0.9854614620545389 75.52899759027348 4.487653824098646 0.03636363636363636 0.03636363636363636 0.9897483333199262 109.0744740434431 4.487770702255042 0.03636363636363636 0.03636363636363636 0.9972318816440022 418.0023774450405 4.487776579638927 0.03636363636363636 0.03636363636363636 0.9992409498317957 1539.153786381783 4.487777138788559 — — — 1.051076326841211 × 109

Figures 27 provide seven groups of graphs to illustrate more particular quenching statuses. Figure 2 shows the 2D plots of and changing as rises from the initial moment until quenching occurs. Figure 2(a) corresponds to , and Figure 2(b) corresponds to . Figure 3 sketches the 2D contours of and regarding immediately before quenching happens. Figure 3(a) corresponds to , and Figure 3(b) corresponds to . Figure 4 draws the 2D outlines of and regarding immediately before quenching happens. Figure 4(a) corresponds to , and Figure 4(b) corresponds to . The three 2D plots in Figure 5 describe the adaptive mechanism before quenching occurs. The first one reflects spatial adaptation, and the last two reflect temporal adaptation. Figure 5(a) describes the trend of the temporal step over time. To distinguish more clearly the variation, we intercept the curves of the period from 4.470000000001407 to 4.487777138788559. The red sign on the blue curve in Figure 5(a) denotes the initial temporal adaptive moment (i.e., 4.486155872847072). The temporal step is uniform before the adaptive moment but becomes nonuniform when varies from 4.486155872847072 to 4.487777138788559. The green sign on the blue curve in Figure 5(a) denotes the initial moment (i.e., 4.487088240486234) of producing the quenching location . The green lines in Figures 5(b) and 5(c) describe the uniform spatial steps in the direction and the direction, respectively. The spatial adaptation occurs at 4.487777138788559 after which the spatial step becomes nonuniform. The red curves with blue spots in Figures 5(b) and 5(c) describe, respectively, the nonuniform spatial steps in the direction and the direction at . Although the two curves are very similar, there is still subtle difference between their data. Figures 67 provide the 3D surfaces of and with respect to spatial variables and at certain special time. Figure 6(a) records the 3D surfaces about , , and , and Figure 6(b) records the 3D surfaces about , , and when . Figure 7(a) records the 3D surfaces about , , and , and Figure 7(b) records the 3D surfaces about , , and when .

##### 4.4. Quenching Case with Convection Term

The second convection term can be written as the two convection functions: and , where and take the constants. Comparing with the first convection term , it is hard to conduct quenching behavior for the second convection term . We make a following configuration to gain a quenching status. The initial step (step) is , and the initial step is . The four groups of important parameters are set as which is defined as Case 8.0. The parameter statuses of the case are recorded in Table 4. For the case, the initial time of producing the quenching location is earlier than the temporal adaptive time. The former time is 5.596499999999618, which is the first row in Table 4, the latter one is 5.603982012140827, which is the fifth row in Table 4, and the quenching location is . Moreover, it is so transient from the front time of quenching location appearing to the terminative time of quenching location appearing that there are 18 time layers between the two moments. The latter quenching time is , is 0.9835341519876587, and is 73.15674241637485 immediately before quenching occurs, which is the sixth row in Table 4. gets the value 6.637517818582074 × 1011 and blows up when is 5.604272047535358, which is the last row in Table 4.

 5.596499999999618 0.025 0.025 0.9656132480353845 27.20492580514906 5.598499999999617 0.025 0.025 0.9668480910506297 28.76690218778399 5.599999999999616 0.025 0.025 0.9678079528853171 30.06497816241713 5.601999999999615 0.025 0.025 0.969139216203357 32.00073689849199 5.603982012140827 0.025 0.025 0.9708363009713709 34.72704476646155 5.604262757909247 0.025 0.025 0.9835341519876587 73.15674241637485 5.604272047535358 — — — 6.637517818582074 × 1011

This paragraph supplies five sets of figures (i.e., Figures 812) to more vividly depict the quenching information of the case. The 2D curves in Figure 8 describe the relationship between and and the relationship between and . The 2D curves in Figure 9 describe the relationship between and and the relationship between and . The 2D curves in Figure 10 describe the relationship between and and the relationship between and . In Figures 810, (a) corresponds to and (b) corresponds to . Figure 11(a) draws the tendency of the temporal step varying with temporal variable ; Figure 11(b) draws the tendency of the initial spatial step and the adaptive step varying with spatial variable ; and Figure 11(c) draws the tendency of the initial spatial step and the adaptive step varying with spatial variable . For Figure 11(a), to distinguish more clearly the variation, we intercept the curves of the period from 5.580499999999627 to 5.604272047535358; and the red spot on the blue curve expresses the earliest temporal adaptive time (i.e., 5.603982012140827). in the time step becomes 1.185615675885128 × 1012 when is 5.603982012140827. The temporal step is uniform before the adaptive moment but becomes nonuniform since then. The green spot on the blue curve expresses the initial moment of producing the quenching location , which is 5.604262757909247. In Figures 11(b) and 11(c), the green lines indicate the uniform initial spatial steps while the red curves with blue signs indicate the nonuniform adaptive spatial steps. Lastly, we provide two 3D surfaces in Figure 12 to state quenching statuses of and when . Figure 12(a) denotes the relationship between , , and , and Figure 12(b) denotes the relationship between , , and . Figure 12(a) records the 3D surfaces about , , and , and Figure 12(b) records the 3D surfaces about , , and when . There is a narrowband projecting the rear plane of the 3D axis box in Figure 12(b). This is because some larger values of