Adaptive High-Order Finite Difference Analysis of 2D Quenching- Type Convection-Reaction-Diffusion Equation

Quenching characteristics based on the two-dimensional (2D) nonlinear unsteady convection-reaction-diffusion equation are creatively researched. The study develops a 2D compact finite difference scheme constructed by using the first and the second central difference operator to approximate the first-order and the second-order spatial derivative, Taylor series expansion rule, and the reminder-correction method to approximate the three-order and the four-order spatial derivative, respectively, and the forward difference scheme to discretize temporal derivative, which brings the accuracy resulted meanwhile. Influences of degenerate parameter, convection parameter, and the length of the rectangle definition domain on quenching behaviors and performances of special quenching cases are discussed and evaluated by using the proposed scheme on the adaptive grid. It is feasible for the paper to offer potential support for further research on quenching problem.


Introduction
As a common class of thermodynamic problems, quenching phenomena research has been widely applied in engineering area. There exist some typical researches including the flow, thermal properties, and its working mechanism of an object in heat associated with quenching phenomena [1,2], the forced convection quench of hot sheets in super-cooled liquid [3], the quench process speed-up, and its quench characteristics in porous interface processing [4,5]. In recent years, more and more scholars concern quenching problem built on parabolic equations [6][7][8][9][10], in which the authors depended on the different parabolic equations to analyze corresponding quench features. Additionally, quenching phenomena based on parabolic equation systems have attracted researchers' attentions [11][12][13]. Nonlinear degeneration singularity reaction-diffusion equations or convectionreaction-diffusion equations, as a branch of parabolic equations, have been usually employed to handle quenching problems. References [14][15][16] focused on the quench-ing phenomena of the nonlinear degeneration singularity reaction-diffusion equations whereas Refs. [17,18] focused on the quenching phenomena of the nonlinear degeneration singularity convection-reaction-diffusion equations. Selcuk discussed quench performances of solution and its time derivative and estimation of quenching time under special conditions [14]. Ge et al. devoted to an adaptive compact difference scheme to solve quench problems [15]. Beauregard analyzed quenching properties of the fractional Kawarada equation by using a new numerical method and proved the proposed method monotonic, nonnegative, and linearly stable [16]. Zhou et al. theoretically investigated quenching characteristics of the nonlinear degeneration singularity convection-reaction-diffusion equation [17]. Zhu and Rui constructed a high-accuracy method with adaptation mechanism for studying quenching behaviors and analyzed the influence of important elements on quenching [18].
References [19][20][21] used numerical method with adaptive algorithm to research degeneration singularity problem of the two-dimensional nonlinear reaction-diffusion equations. A Peaceman-Rachford finite difference scheme based on semiadaptive grids with stability and convergence under certain condition was posed to analyze the quenching behaviors of the 2D reaction-diffusion equation in Ref. [19]. Numerical analyses showed that the method is very efficient and reliable. An adaptive Peaceman-Rachford-Strang splitting method on exponentially evolving meshes was constructed to resolve the 2D problem of quenching type in Ref. [20]. Computational experiments stated that the method is of the satisfactory effectiveness, efficiency, and numerical stability. An exponential splitting method with spatial semidiscretization on arbitrary nonuniform meshes was investigated to explore the 2D degenerate stochastic Kawarada equation in Ref. [21]. Zhu and Ge developed an adaptive ADI strategy to simulate quenching phenomena based on 2D convection-reaction-diffusion equation [22]. Numerical cases illustrated that the method is of the positivity, monotonicity, and numerical stability. The low-accuracy strategies were adopted in the most aforementioned researches [6][7][8][9][10][11][12][13][14][15][19][20][21]. There are only Refs. [15,18,22] delivering the high-accuracy strategies, in which the first two refer to the one-dimensional problem and the last one refers to the two-dimensional problem. It is well known that the development of high-order finite difference scheme of the 2D nonlinear convection-reaction-diffusion equations is maturing [23][24][25]. Reference [23] developed a new highorder difference approach for the 2D convection-reactiondiffusion equation with a small diffusivity, which firstly posed the 1D high-accuracy difference scheme and then extended it to the 2D case a nine-point stencil by using alternating direction algorithm and last resolve the 2D steady incompressible Navier-Stokes equations. In order to simulate groundwater pollution problems, Li et al. put forward a fourth-order compact difference scheme with unconditional stability and the second-order temporal accuracy, the fourth-order spatial accuracy built on the 2D convection-reaction-diffusion model, which was efficient through experiments [24]. Wu and Zhai combined exponential transformation and quadratic interpolation polynomial with Padé approximation to study the 2D time-fractional convection-dominated diffusion equation, which owned the higher accuracy and used alternating direction implicit algorithm to alleviate computational amount [25].
In summary, it is easily found that there are much more low-order strategies than high-order strategies to solve degenerate singularity problems of the 2D reactiondiffusion equation. Especially, high-order difference schemes received few attentions for solving quenching problem of the two-dimensional convection-reaction-diffusion equation. It is the high-precision algorithms that have advantages in dealing with such problems because of its high accuracy and efficiency. So it is a meaningful try to explore the high-order compact difference schemes on adaptation mesh for analyzing quenching problems of the 2D degenerate singular reaction-diffusion equation with convection function. Going down this idea, we extend this study from Ref. [18] and construct a compact difference scheme on adaptive grid for solving the 2D unsteady convection-reaction-diffusion equation to explain the corresponding quenching phenomena. According to Ref. [18], we extend its 1D high-order compact finite difference scheme on adaptive mesh to a 2D strategy and use it to analyze the quenching behaviors of the 2D convection-reaction-diffusion combustion model. There exist three contributions of this paper. Firstly, we represent a new 2D high-order compact difference scheme for solving the corresponding unsteady convection-reaction-diffusion equation and give its accuracy performances. Secondly, we apply the scheme to explain the 2D reaction-diffusion of quenching type with and without convection function, respectively. Thirdly, we investigate a series of quenching characteristics including quenching time, quenching location, and so on, from which we can discover impacts of the parameters q, b, and α (β) on quenching behaviors. 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 introduces adaptive mesh algorithm. Section 4 stimulates some typical numerical samples to explore and explain quenching problems. Section 5 draws the conclusion.

2D Nonlinear Convection-Reaction-Diffusion Equation.
The typical convection-reaction-diffusion equation of quenching type is written as Its boundary conditions are Its initial conditions is This semilinear degenerate problem model involving two spatial dimensions is regarded as Eq. (1) with the boundary and initial conditions of Eqs. (2) and (3). The solution uðx, y , tÞ represents the temperature of the combustion chamber. Ω ∈ ð0, αÞ × ð0, βÞ refers to the smooth domain of the rectangle combustor in which α > 0, β > 0, and combustion chamber sizes α and β are the length of the definition area, ∂Ω, and is its boundary. c 1 ðx, yÞ and c 2 ðx, yÞ are the convection functions of x and y, and σðx, yÞ = ðx 2 + y 2 Þ q/2 is the degeneracy function and degeneration parameter q ≥ 0. The singularity source f ð uÞ = 1/ð1 − uÞ θ is strictly increasing for 0 ≤ u 0 < 1 withf ð0Þ = f 0 > 0, lim u→1 − f ðuÞ = ∞, and singularity parameter θ is the power of the singular source term 1/ð1 − uÞ θ .
With the aid of the intermediate variablesx andŷ, we replacex = x/α andŷ = y/β into Eq. (1) which can be defined as 2 Advances in Mathematical Physics For the convenience of expression, we use x and y instead ofx andŷ. Then, the above formulation can be rearranged as where σðx, yÞ = ðα 2 x 2 + β 2 y 2 Þ q/2 , Ω ∈ ð0, 1Þ × ð0, 1Þ: Correspondingly, the boundary conditions of Eq. (5) are The initial conditions of Eq. (5) is In the physical application, we rely on Eqs. (5)-(7) to compute u and u t . Through observing a large number of values of u and u t , we can capture quenching moment, i.e., quenching occurs when u infinitely close to 1 and u t becomes so huge that its value blows up.
2.2. The Proposed Compact Difference Scheme. According to the idea given in [18], we employ the first and the second central difference operator to approximate the first-order and the second-order spatial derivatives of x − direction and y − direction, respectively, and the forward difference operator to discrete temporal derivative. The proposed high-accuracy finite difference scheme of Eq. (5) is deduced below. After the first-order and the second-order spatial derivatives of x − direction and y − direction are discretized on the nonuniform mesh, respectively, a scheme approximating the 2D unsteady convection-diffusion equation dispersed at point ð x k , y j Þ can be written as where We use Taylor series to obtain the derivative expansions of x − direction: ðu xxx Þ k,j and ðu xxxx Þ k,j , and the derivative expansions of y − direction: ðu yyy Þ k,j and ðu yyyy Þ k,j from Eq. (5). The four expressions are substituted into the corresponding terms in Eq. (8).

Advances in Mathematical Physics
The second-order derivative of time from Eq. (8) is u tt = ð1/α 2 σÞu txx + ð1/β 2 σÞu tyy + ðc 1 /ασÞu tx + ðc 2 /βσÞu ty + ð f t / σÞ. By using the forward difference scheme, we can get Relying on the second-order backward Euler difference scheme, we can get Provided OðΔÞ = OðΔ + τ 2 Þ, Eq. (14) can be derived as The first-order and the second-order derivate of u with regard to space variables x andy are discretized by the central difference scheme. The first-order derivate of u with regard to time variable t is discretized by the forward difference scheme. Subsequently, after the difference approximations are carried out to ðu tx Þ n k,j , ðu txx Þ n k,j , ðu ty Þ n k,j , ðu tyy Þ n k,j , ðu xy Þ n k,j , ðu xxy Þ n k,j , ðu xyy Þ n k,j , and ðu xxyy Þ n k,j , a linear system is formed as where b 1y = −2τ n + βc 2 τ n y L ð Þ κ 3 − 2β 2 σy L κ 4y + 4β 2 σκ 5y Through the prior deducing procedure, it can be seen that the truncation error of the Eq. (15) is- it owns spatial accuracy of fourth order and temporal accuracy of second order.

Adaptive Grid Structure
3.1. Adaptive Grid Structure in Time. The adaptation mesh technique in the temporal and the spatial direction is employed for solving the 2D problem of quenching type. Obviously, the 2D adaptive mesh algorithm can be derived from the 1D case. The arc-length monitor function of the temporal derivative function resulting from equal distribution principle is used to design the temporal and the spatial moving mesh algorithm, respectively. From Ref. [19,20], a self-adaptive grid in the time direction is given as follows. A maximized ratio equation for computing the adaptive temporal interval τ n is The standard uniform central difference formula is substituted for the ∂ 2 u/∂t 2 in Eq. (17). Of course, when the prior temporal step τ n−1 is given, the current temporal step τ n can be calculated through Eq. (17).

Adaptive Grid Structure in Space.
Let a spatial adaptation algorithm be deduced as follows. W x ðx, tÞ is taken as the monitor function in x − direction and W y ðy, tÞ is taken as the monitor function in y − direction.

Advances in Mathematical Physics
fξ k ðx k , tÞg K k=0 and fζ j ðy j , tÞg J j=0 point to the set of 2D original spatial mesh points in x − direction and y − direction. The 2D computational area is transformed as fx k ðξ k , tÞg K k=0 in x − direction and fy j ðζ j , tÞg J j=0 in y − direction through mesh movement. The continuous function W x ðx, tÞ > 0 and W y ðy, tÞ > 0 in the definition area ðx, yÞ ∈ ð0, 1Þ × ð0, 1Þ equally distributes over the 2D mesh refreshed. In fact, the 1D adaptive mesh algorithm individually carries out in each spatial direction. According to [19,20], we can get the semi-implicit scheme in Similarly, we can get the semi-implicit scheme in y − direction 3.3. Iterative Adaption Algorithm. This paragraph shows the iterative process of the presented method in the paper.
Step 1. Depending on the spatial mesh fx ðmÞ k g K k=0 , fy ðmÞ j g J j=0 at the nth time layer, we can obtain the corresponding monitor functions W x ðx, tÞ and W y ðy, tÞ via Eqs. (18) and (19).
Step 2. The point set fx  Step 4. Repeat steps from the first to the third till quenching occurs, or the solution converges to a steady solution.

Simulation Demonstration
4.1. Numerical Case. This is a common numerical sample utilized to evaluate the performance of Eq. (15) measure by Eq.
In the paper, T = 0:5, τ n = 0:01, p = 100. The initial and boundary values can be computed through the exact solution. The tested results display in Table 1, in which there are four criteria containing Max. error, Aver. error, CPU time, and Conv. rate. Max. error means the maximal absolute error between analysis solution and difference solution; Aver. error means the average absolute error; CPU time means the running time of system; Conv. rate means convergence rate between the two interfacing mesh numbers. The nonuniform mesh follows the rule: fk/K + ðλ x /πÞ sin ðπk/KÞg K k=0 in x − direction and fj/J + ðλ y /πÞ sin ðπj/JÞg J j=0 in y − direction, where K and J point to intervals of spatial direction and K = J, whereas λ x and λ y mean telescopic transformation coefficient and λ x = λ y . There exist six difference schemes that are compared on another for the example in Table 1, which are the schemes in Ref. [26] on the uniform and the nonuniform grid, the schemes in Ref. [27] on the uniform and the nonuniform grid, and the proposed schemes on the uniform and the nonuniform grid. The schemes in Ref. [27] and the proposed schemes all chose the five spatial intervals 16, 32, 64, 128, 256, and λ x = λ y = 0:7 and the four aforementioned criteria. Only the schemes in Ref. [26] choose the five spatial intervals 16, 32, 64, 128, and λ x = λ y = 0:4 and the three aforementioned criteria with the exception of Aver. error.
With regard to Max. error and Aver. error, the schemes on the uniform grid are inferior to the schemes on the nonuniform grid; the proposed schemes are superior to the other schemes, and the schemes (nonuniform) in Ref. [27] are superior to the scheme (nonuniform) in Ref. [26], and the tendency becomes more obviously as KðJÞ rises. When KðJ Þ is 128, the maximal absolute error of the scheme in Ref. [27] (uniform) is 3:1469 × 10 −4 , and the maximal absolute error of the scheme in Ref. [27] (nonuniform) is 9:3004 × 10 −7 ; the maximal absolute error of the scheme in Ref. [26] 6 Advances in Mathematical Physics (uniform) is 3:15 × 10 −4 , and the maximal absolute error of the scheme in Ref. [26] (nonuniform) is 4:42 × 10 −5 ; the maximal absolute error of the proposed scheme (uniform) is 3:1447 × 10 −4 , and the maximal absolute error of the proposed scheme (nonuniform) is 6:9033 × 10 −7 . When KðJÞ is 256, the average error of the scheme in Ref. [27] (uniform) is 6:1742 × 10 −7 , and the average error of the scheme in Ref. [27] (nonuniform) is 6:4492 × 10 −9 ; the average error of the proposed scheme (uniform) is 5:9061 × 10 −7 , and the average error of the proposed scheme (nonuniform) is 3:4536 × 10 −9 .
In terms of CPU time, the scheme (nonuniform) in Ref. [27] spends the most time among the schemes when KðJÞ is 32, 64, 128, and 256, and the tendency becomes more obviously as KðJÞ rises; the schemes in Ref. [26] spend much more time than the proposed schemes do. Typically, the running time of the schemes on the uniform grid is less than that of the schemes on the nonuniform grid. It just means that the latter complies with a time-for-space rule, which refers to it may improve the spatial accuracy of an algorithm by adding its running time.

Quenching Case without Convection Term.
When c 1 ðx, yÞ = 0 and c 2 ðx, yÞ = 0, with the initial and boundary conditions Eqs. (6) and (7), the original Eq. (5) can be described as For the situation of c 1 ðx, yÞ = 0 and c 2 ðx, yÞ = 0, we take advantage of Eq. (15) to approximate Eq. (24) and gain quenching case without convection term. We set the three parameters q = 0, a = β = ffiffiffiffiffi 10 p , and θ = 1 to the case which had been illustrated in Refs [19,20]. In this case, the initial temporal step is configured as τ 0 = 4:166677 × 10 −5 , and the initial space step is configured as h 0 = 1/90. Table 2 offers quenching performances of the three different schemes, which are from Refs [19,20] besides the proposed method on adaptive mesh, for q = 0, a = β = ffiffiffiffiffi 10 p , and θ = 1 without convection term. The criteria in Table 2 are explained as follows. ðx max , y max Þ refers to the point of quenching location, T max refers to the quenching time, max u refers to the maximum temperature value immediately before quenching happens, and max u t refers to the maximum variation of temperature with respect to time immediately before quenching happens. For the three schemes, their quenching locations are all ð0:5, 0:5Þ. There are subtle differences among the three schemes for quenching time and max u. T max of the proposed scheme is 0.5990174307488624, and max u of the proposed scheme is 0.9901182463562488. But there are greater differences among the three schemes for max u t . max u t of the proposed scheme is 83.20947696438735 while those of the schemes in Refs [19,20] are 148.887767 and 1249.917563. There is no effect of the differences on quenching research. Figure 1 offers the three-dimensional scenes of u and u t with regard to spatial variables x, y when the time t is 0.5990174307488624, respectively. From the two 3D surfaces at the time T max = 0:5990174307488624 that is immediately before quenching occurs, it is found that u ⟶ 1 − and u t ⟶ 84 − at the quenching location ðx max , y max Þ = ð0:5, 0:5Þ and the quenching time T max = 0:5990174307488624.

Quenching Case with Convection Term
We will investigate what role degeneration parameter q, convection parameter b, and combustion chamber size α(β) in Eqs. (5)-(7) play during the quenching process and enumerate the representative Case 5.0 to show special quenching features when c 1 ðx, yÞ = b x /ð1 + αx + βyÞ and c 2 ðx, yÞ = b y /ð1 + αx + βyÞ, where b x and b y take the constants and convection parameter b = b x = b y . When singularity parameter θ is larger than 1, the quenching situation is complex. Therefore, we only consider quenching cases of θ = 1. In numerical samples, b x and b y range from -200 to 200; α 2 range from 15 to 300000. These are the same to Section 4.4. To better illustrate the problems, we choose some representative data from lots of tests to discuss quenching phenomena in the next paragraphs. For the convection model b/ð1 + αx + βyÞ, we set h 0 = 1/70 as the initial x − (y − ) step and τ 0 = 0:0008 as the initial t − step.
Tables 3-6 offer some typical quenching cases to demonstrate the relationship between the three groups of parameters and quenching phenomena. Case 5.0 recorded as a reference in Table 3, i.e., q = 1, b x = b y = 3, α = β = 15, is compared with other cases in Tables 3-5. Table 6 displays the specific quenching data of Case 5.0. Figures 1-3 show the quenching statuses of Case 5.0 relative to u, u t , x, y, and t.
Except Case 5.0, there are ten quenching cases chosen in Table 3. The ten items list as follows: Although the program may run when q ≥ 1:5, it is hard to form quenching behaviors. If q belongs to the definition domain ½0:1, 1:4, the quenching phenomena will occur. From serial numerical cases, we can easily see the performances of quenching location and time. The quenching location ðx max , y max Þ declines as q evolves and finally converges to ð0:04285714285714286, 0:04285714285714286Þ. There is an intermediate point q = 0:9 dividing the definition domain as ½0:1, 0:9 and ½1:0, 1:4. The quenching location x max and y max do decline as q evolves in the first subdomain and converge to ð0:04285714285714286, 0:04285714285714286Þ in the second subdomain. The quenching time T max reaches the maximum 9.686902288384093 when q = 1:4. T max and qare in positive proportion when q in ½0:1, 1:4.
We take Case 5.0 to demonstrate specific quenching phenomena. In this case, quenching occurs at T 5 * = 3:428404873627137, x 5 * = 0:04285714285714286, and y 5 * = 0:04285714285714286 with the parameters q = 1, b x = b y = 3, α = β = 15. The next context declares more details about quenching. Let us observe the performances of the    Table 6 and Figures 2-4. The first four items reflect the four states before quenching occurs in Table 6. The last item represents the quenching state in which u is more than 1 and reaches 17.07446290189856, u t blows up and reaches 1:200141619871376 × 10 9 , the position of max x,y u is (0:5336298640841476 × 10 -5 , 0.2902651144597507), and the position of max x,y ðu t Þ is (0:5702149685882448 × 10 -4 , 0:5677032305886261 × 10 -4 ) when t = 3:428420657560073.  Figure 4 depicts curve of the adaptive temporal steps as the time variable progresses and graphs of the adaptive spatial steps changing as the space variables change. Figures 4(a) and 4(b) portray the distributions of spatial steps. The plots in Figure 4(a) are similar to those in Figure 4(b), but there are some subtle differences in their specific data. The blue curve depicts the tendency of spatial steps varying before quenching occurs, which underlines x 5 * (y 5 * ). The green line refers to spatial step distribution of uniform mesh with nonadaption, which is extracted from the initial temporal layer. Mesh adaptation is stimulated at the position of red square sign, and τ n becomes gradually smaller and falls rapidly toward its floor whent ⟶ T 5 * in Figure 4(c). The red square mark represents the moment t = 3:425399077950409 and the temporal step τ n = 5:990779506555022 × 10 -4 when the process approximates to quenching time. The time adaption is triggered before the moment T 5 * and lasts through the rest course of calculation. So this quenching phenomenon is caught with the features of max x,y u ⟶ 1 − and max x,y ðu t Þ ⟶ 126 − as t reaches T 5 * .  Figures 5(e) and 5(f) denote the two threedimensional plots of u and u t when t = T 5 * , respectively. Thought the 3D plots, we can get richer information. The first group is the three-dimensional views from both u and u t at the first 500th temporal layers in the excursion, in which max x,y u ⟶ 0:2215 and max x,y ðu t Þ ⟶ 12:9459 when t = 0:40. The second group is the three-dimensional views from both u and u t including the penultimate temporal steps before T 5 * , in whichmax x,y u ⟶ 0:9863 and max x,y ðu t Þ ⟶ 68:0806 when t = 3:428342645817735. The third group is the threedimensional views from both u and u t immediately before quenching happens, in which max x,y u ⟶ 0:9921 and max x,y ðu t Þ ⟶ 125:6266 when t = 3:428404873627137. During the temporal layers' moving forward, the solution u changes smoothly and then almost arrives at the peak value while t approximates to T 5 * . Its peak is the maximal value but before quenching time and location. There is some subtle perturbation at the initial temporal axis of the left boundary for the change rate of the solution. Afterward, the temporal derivative varies smoothly, and their maximums of each time axis increase steadily. While t approaches to the quenching time T 5 * , the maximum u t also increases rapidly and infinitely and then blows up at the next time layer of quenching time.

Quenching Example with Convection Term b/ðαx + βyÞ.
Similarly, we need to only consider θ = 1:0 and conduct a series of simulation experiments to investigate the 2D quenching regularity for the convection term b/ðαx + βyÞ which is related to the three groups of elements: q and b x and b y , α, and β and exemplify the quenching behaviors of the representative Case 6.0. The convection term b/ðαx + βy Þ should written as c 1 ðx, yÞ = b x /ðαx + βyÞ and c 2 ðx, yÞ = b y / ðαx + βyÞ, where b x and b y take the constants and convection parameter b = b x = b y . For Case 6.0, the initial x − step is h 0 = 1/80, and the initial t − step is τ 0 = 0:001. We choose some representative data to record in Tables 7-10 from the simulation results. Case 6.0 is regarded as the reference standard The effect of q on quenching problem can be illustrated by continual experiments. A reasonable quenching behavior will occur when q = 0 or q takes a value in ½0:4, 1:4. When q is 0.1, 0.2, 0.3, or q ≥ 1:5, it does not produce quenching phenomena. There are quenching characteristics for q: quenching location x max is not equal to y max when q takes 0, 0.5 while x max is equal to y max in the other cases; ðx max , y max Þ is ð0:0125, 0:0125Þ when q takes 0.6, 0.7, and 0.8; ð x max , y max Þ is ð0:025, 0:025Þ when q is in ½0:9, 1:4. Quenching time T max increases as q increases in the domain of q. Table 7 gives some of these quenching cases, in which Case     the ten special cases, and we can observe the quenching status as follows. As far as quenching spatial characteristic is concerned, quenching location does not monotonously increase as α 2 increases when α 2 is in the definition domain. For example, quenching location reaches the maximum ð0:5375 , 0:5375Þ when α 2 = 10 and quenching location keeps the coordination point of ð0:0125, 0:0125Þ when α 2 is in ½10780 , 3:0 × 10 6 . Next, quenching temporal characteristic is concerned. In fact, there does not exist strict linear relationship between quenching time and α 2 . Comparatively, quenching time is smaller in the former half than in the latter half of the domain of α 2 . When the measurement scale is enlarged in the domain of α 2 , especially when α 2 reaches 19991, T max increases with the increase of α 2 . Quenching phenomena of Case 6.0 can be considered in the following text. Its quenching time appears at T 6 * = 4:020449234557191, and its quenching position appears at ð x 6 * ,y 6 * Þ = ð0:0125, 0:0125Þ in this situation. Table 10 shows six stages around quenching for Case 6.0. The items from the first to the fifth describe the five stages of representation before the occurrence of quenching, and the last item just notes the quenching moment. max x,y u is 5:358959032595640 × 10 16 at (2:0588245639243 × 10 −4 , 9:654095827421034 × 10 −1 ), and max x,y ðu t Þ is 7:966445413263532 × 10 19 at (2:080224198836234 × 10 −5 , 9:824584434865121 × 10 −1 ) when quenching occurs. A comprehensive statement of quenching states with the parameters q = 0:8, b x = b y = −2, α = β = ffiffiffiffiffiffiffiffiffiffiffi ffi 11000 p is recorded in the next paragraphs. Although there is slight perturbation in the left boundary of u t , it does not influences on distribution of u t . It is evident that Figures 6(a) and 6(b) give two pairs curves of both max x,y u and max x,y ðu t Þ. The first one is for the distribution of the solution max x,y uðtÞ as t varies, and the second one is for the distribution of max x,y u t ðtÞ as t increases. There is small perturbation near the start-up of the adaptive procedure for the left figure of max x,y ðu t Þ. The four figures in Figure 7 give function relationship physical quantities betweenu, u t , and spatial variables. Figure 7(a) paints the contour of (x,u), Figure 7(b) paints the contour of (x,u t ), Figure 7 There are three graphics in Figure 8, which represent three function relationships between spatial step in x − direction and x, between spatial step in y − direction and y,  Figure 9: (a) The 3D plot of u when t = T 6 * and (b) the 3D plot of u t when t = T 6 * for Case 6.0. and between temporal step and t, respectively. Figure 8(a) depicts two curves, in which one is between initial spatial steps h k and x, and another is between quenching spatial steps h k and x. Figure 8(b) also describes two curves, in which one is between initial spatial steps h k and y, and another is between quenching spatial steps h k and y. The green lines with gradient marks reflect the initial spatial step distribution that is uniform regardless of x − direction or y − direction. The red curves covered by blue squares reflect quenching spatial steps that are self-adaptive when quenching occurs. There is a blue curve with a red square marks in Figure 8(c), which portrays temporal stepτ n varies as t changes from 0 to T 6 * . τ n keeps 0.001 when 0 ≤ t < 4:019595682766669, but afterτ n becomes 5:956827669919562 × 10 −4 when t = 4:019595682766669 which signed as a red square,τ n sharply declines up to 3:286760944343214 × 10 −4 when t = T 6 * .
The three-dimensional surfaces can render more reliable information for both the solution u and the time derivative u t which display in Figure 9. The group of surfaces is the threedimensional views of both u and u t when t is equal to T 6 * . It is smooth and steady that the solution u carries forward at the temporal axis in Figure 9(a). Moreover, it is almost at the middle position of each time axis that the case gives birth to the peak value of u. When the position ðx 6 * ,y 6 * Þ is at ð 0:0125, 0:0125Þ immediately before quenching occurs, the peak value of the solution approaches 1. The threedimensional representation of Figure 9(b) reveals the evolution of the temporal derivative with regard to spatial variables. The temporal derivative of the solution varies smoothly and reaches its maximums at the quenching domain ð0:0125, 0:0125Þ. At last, the temporal derivative becomes infinite when quenching occurs at another location (2:080224198836234 × 10 −5 , 9:824584434865121 × 10 −1 ), which is not observed in Figure 9(b). It is obvious for Case 6.0 that quenching occurs based on max x,y u ⟶ 1 − and max x,y ðu t Þ ⟶ 32 − as t reaches T 6 * .

Conclusion
Relying on the analyses in this paper, we can sort up the relationship between the parameters q, b x (b y ), a(β), and quenching behaviors. There exist four aspects concluded as follows. First, the degenerate function acts within the normal range if the degeneracy parameter q is larger than or equal to 0.4 but not more than 1.3, and quenching time increases as q increases. Second, there is a special point b ∧ dividing the definition domain of convection parameter b as the left and the right subdomain, in which quenching phenomena occur normally and have a variety of features. Quenching time decreases as b increases in the left subdomain, and quenching time increases as b increases in the right subdomain. Quenching location either decreases or convergences upon a certain value as bincreases in the domain of b. Third, under the condition of α/β = 1, the influence of α 2 on quenching result is researched. There exists still a special point α ∧ categorizing the definition domain of α 2 as two subdomains, in which quenching phenomena occur normally and display different spatial effects. Specially, with the rise of α 2 , quenching location does not increases monotonously in the left hand of α ∧ and tends to a fixed coordinates in the right hand of α ∧ . Additionally, if we rely on a small scale to observe the domain of α 2 , then there is no linear relationship between quenching time and α 2 . We investigate the definition by virtue of a big scale to find that quenching time also increases when α 2 become larger, especially when α 2 reaches certain value; quenching time increases as α 2 rises. Fourth, it must be a good choice to set θ as 1. Through experiments, it is hard to produce quenching phenomena that can be formed when θ ≥ 2. From the result analyses, it can be discovered that it is meaningful to initially investigate the 2D singularity degeneracy problem of quenching type based on the unsteady convection-reaction-diffusion equation by using the highorder difference scheme, which will probably lead to the potentially support to research the next quenching problems.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.