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

Authors: Xiaoliang Zhu; Yongbin Ge
Journal: Advances in Mathematical Physics
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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–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–13]. Nonlinear degeneration singularity reaction-diffusion equations or convection-reaction-diffusion equations, as a branch of parabolic equations, have been usually employed to handle quenching problems. References [14–16] focused on the quenching 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 [6–18] concentrated on the quenching phenomena of the 1D heat equations. The quenching phenomena of the 2D heat equations are described as follows. References [19–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–15, 19–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–25]. Reference [23] developed a new high-order difference approach for the 2D convection-reaction-diffusion 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 reaction-diffusion 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.


## 2. 2D Problem and Difference Scheme

### 2.1. 2D Nonlinear Convection-Reaction-Diffusion Equation

The typical convection-reaction-diffusion equation of quenching type is written as(1)σx,yut=uxx+uyy+c1x,yux+c2x,yuy+fu,x,y∈Ω,t∈0,T.Its boundary conditions are(2)u0,y,t=uα,y,t=ux,0,t=ux,β,t=0,t∈0,T.Its initial conditions is(3)ux,y,0=u0x,y,x,y∈0,α×0,β.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 ux,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. c1x,y and c2x,y are the convection functions of x and y, and σx,y=x2+y2q/2 is the degeneracy function and degeneration parameter q≥0. The singularity source fu=1/1−uθ is strictly increasing for 0≤u0<1withf0=f0>0,limu→1−fu=∞, and singularity parameter θ is the power of the singular source term 1/1−uθ.With the aid of the intermediate variablesx^ and y^, we replace x^=x/α and y^=y/β into Eq. (1) which can be defined as(4)σαx^,βy^ut=1α2ux^x^+1β2uy^y^+c1αx^,βy^αux^+c2αx^,βy^βuy^+fu,αx^,βy^∈Ω,t∈0,T.For the convenience of expression, we usex and y instead of x^ and y^. Then, the above formulation can be rearranged as(5)σx,yut=1α2uxx+1β2uyy+c1x,yαux+c2x,yβuy+fu,x,y∈Ω,t∈0,T,where σx,y=α2x2+β2y2q/2,Ω∈0,1×0,1. Correspondingly, the boundary conditions of Eq. (5) are(6)u0,y,t=u1,y,t=ux,0,t=ux,1,t=0,t∈0,T.The initial conditions of Eq. (5) is(7)ux,y,0=u0x,y,x,y∈0,1×0,1.In the physical application, we rely on Eqs. (5)–(7) to compute u and ut. Through observing a large number of values of u and ut, we can capture quenching moment, i.e., quenching occurs when u infinitely close to 1 and ut 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 xk,yj can be written as(8)2−αxR−xLc12α2σδx2+c1ασδx+2−βyR−yLc22β2σδy2+c2βσδyuk,j+Ψ1xuxxxk,j+Ψ2xuxxxxk,j+Ψ3xuxxxxxk,j+Ψ1yuyyyk,j+Ψ2yuyyyyk,j+Ψ3yuyyyyyk,j=ut−fσuk,j,l+OxR5+xL5xR+xL+OyR5+yL5yR+yL,where(9)Ψ1x=−2xR−xL+αxRxLc16α2σ,Ψ2x=−2xR2−xRxL+xL2+αxR−xLxRxLc124α2σ,Ψ3x=−2xR−xLxR2+xL2+αxR3xL−xR2xL2+xL3xRc1120α2σ,Ψ1y=−2yR−yL+βyRyLc26β2σ,Ψ2y=−2yR2−yRyL+yL2+βyR−yLyRyLc224β2σ,Ψ3y=−2yR−yLyR2+yL2+βyR3yL−yR2yL2+yL3yRc2120β2σ.We use Taylor series to obtain the derivative expansions ofx−direction: uxxxk,j and uxxxxk,j, and the derivative expansions of y−direction: uyyyk,j and uyyyyk,j from Eq. (5). The four expressions are substituted into the corresponding terms in Eq. (8). OΔ=OxR5+xL5/xR+xL+OyR5+yL5/yR+yL−Ψ2yuyyyyk,j−Ψ3yuyyyyyk,j, which is the truncation error of Eq. (8). Omitting the truncation errors, we consider the situation of point k,j,n for Eq. (8)(10)κ1xδx+κ2xδx2+κ1yδy+κ2yδy2uk,jn+κ3ut+κ4xutx+κ5xutxx+κ4yuty+κ5yutyy+κ6uxy+κ7xuxxy+κ7yuxyy+κ8uxxyyk,jn=−1σ+κ9xδx+κ10xδx2+κ9yδy+κ10yδy2fk,jn+OΔ,where(11)κ1x=c1ασ−αΨ1x−αΨ2xc1c1x−αΨ2xc1xx−β2αΨ1y−βΨ2yc2c1y−β2Ψ2yαc1yy,κ2x=2−αxR−xLc12α2σ−αΨ1x−αΨ2xc1c1−2αΨ2xc1x,κ1y=c2βσ−βΨ1y−βΨ2yc2c2y−βΨ2yc2yy−α2βΨ1x−αΨ2xc1c2x−α2Ψ2xβc2xx,κ2y=2−βyR−yLc22β2σ−βΨ1y−βΨ2yc2c2−2βΨ2yc2y,κ3=α2Ψ1x−αΨ2xc1σx+α2Ψ2xσxx+β2Ψ1y−βΨ2yc2σy+β2Ψ2yσyy−1,κ4x=α2Ψ1x−αΨ2xc1σ+2α2Ψ2xσx,κ4y=β2Ψ1y−βΨ2yc2σ+2β2Ψ2yσy,κ5x=α2Ψ2xσ,κ5y=β2Ψ2yσ,κ6=−α2βΨ1x−αΨ2xc1c2−2α2Ψ2xβc2x−β2αΨ1y−βΨ2yc2c1−2β2Ψ2yαc1y,κ7x=−α2Ψ2xβc2−β2α2Ψ1y−βΨ2yc2,κ7y=−β2Ψ2yαc1−α2β2Ψ1x−αΨ2xc1,κ8=−α2Ψ2xβ2−β2Ψ2yα2,κ9x=α2Ψ1x−αΨ2xc1,κ9y=β2Ψ1y−βΨ2yc2,κ10x=α2Ψ2x,κ10y=β2Ψ2y.The second-order derivative of time from Eq. (8) is utt=1/α2σutxx+1/β2σutyy+c1/ασutx+c2/βσuty+ft/σ. By using the forward difference scheme, we can get(12)utk,jn=uk,jn+1−uk,jnτn−τn2uttk,jn+Oτ2=δt+uk,jn−τn2uttk,jn+Oτ2.Relying on the second-order backward Euler difference scheme, we can get(13)ftk,jn=3fk,jn−4fk,jn−1+fk,jn−22τn.ProvidedOΔ=OΔ+τ2, Eq. (14) can be derived as(14)κ1xδx+κ2xδx2+κ1yδy+κ2yδy2uk,jn+κ3δt+uk,jn+κ4x−κ3τnc12ασutx+κ5x−κ3τn2α2σutxx+κ4y−κ3τnc22βσuty+κ5y−κ3τn2β2σutyyk,jn+κ6uxy+κ7xuxxy+κ7yuxyy+κ8uxxyyk,jn=−1σ+κ9xδx+κ10xδx2+κ9yδy+κ10yδy2fk,jn+3κ34σfk,jn−κ3σfk,jn−1+κ34σfk,jn−2+OΔ.The first-order and the second-order derivate ofu 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 utxk,jn, utxxk,jn, utyk,jn, utyyk,jn, uxyk,jn, uxxyk,jn, uxyyk,jn, and uxxyyk,jn, a linear system is formed as(15)b1yuk,j−1n+1+b1xuk−1,jn+1+b0uk,jn+1+b2xuk+1,jn+1+b2yuk,j+1n+1=d1yuk,j−1n+d1xuk−1,jn+d0uk,jn+d2xuk+1,jn+d2yuk,j+1n+d31uk−1,j−1n+d32uk+1,j−1n+d33uk−1,j+1n+d34uk+1,j+1n+s1yfk,j−1n+s1xfk−1,jn+s0fk,jn+s2xfk+1,jn+s2yfk,j+1n+M1fk,jn−1+M2fk,jn−2,where(16)b1y=−2τn+βc2τnyLκ3−2β2σyLκ4y+4β2σκ5y2β2στnyLyR+yL,b1x=2τn+αc1τnxLκ3−2α2σxLκ4x+4α2σκ5x2α2στnxLxR+xL,b0=κ3τn−2α2σκ5x−τnκ3α2στnxRxL−2β2σκ5y−τnκ3β2στnyRyL,b2x=−2τn−αc1τnxRκ3+2α2σxRκ4x+4α2σκ5x2α2στnxRxR+xL,b2y=−2τn−βc2τnyRκ3+2β2σyRκ4y+4β2σκ5y2β2στnyRyR+yL,d1y=κ1yyR+yL−2κ2yyLyR+yL+−2τn+βc2τnyLκ3−2β2σyLκ4y+4β2σκ5y2β2στnyLyR+yL−2κ7xxRxLyR+yL+4κ8xRxLyLyR+yL,d1x=κ1xxR+xL−2κ2xxLxR+xL+−2τn+αc1τnxLκ3−2α2σxLκ4x+4α2σκ5x2α2στnxLxR+xL−2κ7yxR+xLyRyL+4κ8xLxR+xLyRyL,d0=2κ2xxRxL+2κ2yyRyL+κ3τn−2α2σκ5x−τnκ3α2στnxRxL−2β2σκ5y−τnκ3β2στnyRyL−4κ8xRxLyRyL,d2x=κ1xxR+xL−2κ2xxRxR+xL−2τn−αc1τnxRκ3+2α2σxRκ4x+4α2σκ5x2α2στnxRxR+xL+2κ7yxR+xLyRyL+4κ8xRxR+xLyRyL,d2y=−κ1yyR+yL−2κ2yyRyR+yL+−2τn−βc2τnyRκ3+2β2σyRκ4y+4β2σκ5y2β2στnyRyR+yL+2κ7xxRxLyR+yL+4κ8xRxLyRyR+yL,d31=−κ6xR+xLyR+yL+2κ7xxLxR+xLyR+yL+2κ7yxR+xLyLyR+yL−4κ8xLxR+xLyLyR+yL,d32=κ6xR+xLyR+yL+2κ7xxRxR+xLyR+yL−2κ7yxR+xLyLyR+yL−4κ8xRxR+xLyLyR+yL,d33=κ6xR+xLyR+yL−2κ7xxLxR+xLyR+yL+2κ7yxR+xLyRyR+yL−4κ8xLxR+xLyRyR+yL,d34=−κ6xR+xLyR+yL−2κ7xxRxR+xLyR+yL−2κ7yxR+xLyRyR+yL−4κ8xRxR+xLyRyR+yL,s1y=−κ9yyR+yL+2κ10yyLyR+yL,s1x=−κ9xxR+xL+2κ10xxLxR+xL,s0=3κ3−44σ−2κ10xxRxL−2κ10yyRyL,s2x=κ9xxR+xL+2κ10xxRxR+xL,s2y=κ9yyR+yL+2κ10yyRyR+yL,M1=−κ3σ,M2=κ34σ.Through the prior deducing procedure, it can be seen that the truncation error of the Eq. (15) isOτn2+τnxR2−xRxL+xL2+τnyR2−yRyL+yL2+xR5+xL5/xR+xL+yR5+yL5/yR+yL. When h=xR=xL=yR=yL and τn=h2, it owns spatial accuracy of fourth order and temporal accuracy of second order.


## 2.1. 2D Nonlinear Convection-Reaction-Diffusion Equation

The typical convection-reaction-diffusion equation of quenching type is written as(1)σx,yut=uxx+uyy+c1x,yux+c2x,yuy+fu,x,y∈Ω,t∈0,T.Its boundary conditions are(2)u0,y,t=uα,y,t=ux,0,t=ux,β,t=0,t∈0,T.Its initial conditions is(3)ux,y,0=u0x,y,x,y∈0,α×0,β.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 ux,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. c1x,y and c2x,y are the convection functions of x and y, and σx,y=x2+y2q/2 is the degeneracy function and degeneration parameter q≥0. The singularity source fu=1/1−uθ is strictly increasing for 0≤u0<1withf0=f0>0,limu→1−fu=∞, and singularity parameter θ is the power of the singular source term 1/1−uθ.With the aid of the intermediate variablesx^ and y^, we replace x^=x/α and y^=y/β into Eq. (1) which can be defined as(4)σαx^,βy^ut=1α2ux^x^+1β2uy^y^+c1αx^,βy^αux^+c2αx^,βy^βuy^+fu,αx^,βy^∈Ω,t∈0,T.For the convenience of expression, we usex and y instead of x^ and y^. Then, the above formulation can be rearranged as(5)σx,yut=1α2uxx+1β2uyy+c1x,yαux+c2x,yβuy+fu,x,y∈Ω,t∈0,T,where σx,y=α2x2+β2y2q/2,Ω∈0,1×0,1. Correspondingly, the boundary conditions of Eq. (5) are(6)u0,y,t=u1,y,t=ux,0,t=ux,1,t=0,t∈0,T.The initial conditions of Eq. (5) is(7)ux,y,0=u0x,y,x,y∈0,1×0,1.In the physical application, we rely on Eqs. (5)–(7) to compute u and ut. Through observing a large number of values of u and ut, we can capture quenching moment, i.e., quenching occurs when u infinitely close to 1 and ut 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 xk,yj can be written as(8)2−αxR−xLc12α2σδx2+c1ασδx+2−βyR−yLc22β2σδy2+c2βσδyuk,j+Ψ1xuxxxk,j+Ψ2xuxxxxk,j+Ψ3xuxxxxxk,j+Ψ1yuyyyk,j+Ψ2yuyyyyk,j+Ψ3yuyyyyyk,j=ut−fσuk,j,l+OxR5+xL5xR+xL+OyR5+yL5yR+yL,where(9)Ψ1x=−2xR−xL+αxRxLc16α2σ,Ψ2x=−2xR2−xRxL+xL2+αxR−xLxRxLc124α2σ,Ψ3x=−2xR−xLxR2+xL2+αxR3xL−xR2xL2+xL3xRc1120α2σ,Ψ1y=−2yR−yL+βyRyLc26β2σ,Ψ2y=−2yR2−yRyL+yL2+βyR−yLyRyLc224β2σ,Ψ3y=−2yR−yLyR2+yL2+βyR3yL−yR2yL2+yL3yRc2120β2σ.We use Taylor series to obtain the derivative expansions ofx−direction: uxxxk,j and uxxxxk,j, and the derivative expansions of y−direction: uyyyk,j and uyyyyk,j from Eq. (5). The four expressions are substituted into the corresponding terms in Eq. (8). OΔ=OxR5+xL5/xR+xL+OyR5+yL5/yR+yL−Ψ2yuyyyyk,j−Ψ3yuyyyyyk,j, which is the truncation error of Eq. (8). Omitting the truncation errors, we consider the situation of point k,j,n for Eq. (8)(10)κ1xδx+κ2xδx2+κ1yδy+κ2yδy2uk,jn+κ3ut+κ4xutx+κ5xutxx+κ4yuty+κ5yutyy+κ6uxy+κ7xuxxy+κ7yuxyy+κ8uxxyyk,jn=−1σ+κ9xδx+κ10xδx2+κ9yδy+κ10yδy2fk,jn+OΔ,where(11)κ1x=c1ασ−αΨ1x−αΨ2xc1c1x−αΨ2xc1xx−β2αΨ1y−βΨ2yc2c1y−β2Ψ2yαc1yy,κ2x=2−αxR−xLc12α2σ−αΨ1x−αΨ2xc1c1−2αΨ2xc1x,κ1y=c2βσ−βΨ1y−βΨ2yc2c2y−βΨ2yc2yy−α2βΨ1x−αΨ2xc1c2x−α2Ψ2xβc2xx,κ2y=2−βyR−yLc22β2σ−βΨ1y−βΨ2yc2c2−2βΨ2yc2y,κ3=α2Ψ1x−αΨ2xc1σx+α2Ψ2xσxx+β2Ψ1y−βΨ2yc2σy+β2Ψ2yσyy−1,κ4x=α2Ψ1x−αΨ2xc1σ+2α2Ψ2xσx,κ4y=β2Ψ1y−βΨ2yc2σ+2β2Ψ2yσy,κ5x=α2Ψ2xσ,κ5y=β2Ψ2yσ,κ6=−α2βΨ1x−αΨ2xc1c2−2α2Ψ2xβc2x−β2αΨ1y−βΨ2yc2c1−2β2Ψ2yαc1y,κ7x=−α2Ψ2xβc2−β2α2Ψ1y−βΨ2yc2,κ7y=−β2Ψ2yαc1−α2β2Ψ1x−αΨ2xc1,κ8=−α2Ψ2xβ2−β2Ψ2yα2,κ9x=α2Ψ1x−αΨ2xc1,κ9y=β2Ψ1y−βΨ2yc2,κ10x=α2Ψ2x,κ10y=β2Ψ2y.The second-order derivative of time from Eq. (8) is utt=1/α2σutxx+1/β2σutyy+c1/ασutx+c2/βσuty+ft/σ. By using the forward difference scheme, we can get(12)utk,jn=uk,jn+1−uk,jnτn−τn2uttk,jn+Oτ2=δt+uk,jn−τn2uttk,jn+Oτ2.Relying on the second-order backward Euler difference scheme, we can get(13)ftk,jn=3fk,jn−4fk,jn−1+fk,jn−22τn.ProvidedOΔ=OΔ+τ2, Eq. (14) can be derived as(14)κ1xδx+κ2xδx2+κ1yδy+κ2yδy2uk,jn+κ3δt+uk,jn+κ4x−κ3τnc12ασutx+κ5x−κ3τn2α2σutxx+κ4y−κ3τnc22βσuty+κ5y−κ3τn2β2σutyyk,jn+κ6uxy+κ7xuxxy+κ7yuxyy+κ8uxxyyk,jn=−1σ+κ9xδx+κ10xδx2+κ9yδy+κ10yδy2fk,jn+3κ34σfk,jn−κ3σfk,jn−1+κ34σfk,jn−2+OΔ.The first-order and the second-order derivate ofu 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 utxk,jn, utxxk,jn, utyk,jn, utyyk,jn, uxyk,jn, uxxyk,jn, uxyyk,jn, and uxxyyk,jn, a linear system is formed as(15)b1yuk,j−1n+1+b1xuk−1,jn+1+b0uk,jn+1+b2xuk+1,jn+1+b2yuk,j+1n+1=d1yuk,j−1n+d1xuk−1,jn+d0uk,jn+d2xuk+1,jn+d2yuk,j+1n+d31uk−1,j−1n+d32uk+1,j−1n+d33uk−1,j+1n+d34uk+1,j+1n+s1yfk,j−1n+s1xfk−1,jn+s0fk,jn+s2xfk+1,jn+s2yfk,j+1n+M1fk,jn−1+M2fk,jn−2,where(16)b1y=−2τn+βc2τnyLκ3−2β2σyLκ4y+4β2σκ5y2β2στnyLyR+yL,b1x=2τn+αc1τnxLκ3−2α2σxLκ4x+4α2σκ5x2α2στnxLxR+xL,b0=κ3τn−2α2σκ5x−τnκ3α2στnxRxL−2β2σκ5y−τnκ3β2στnyRyL,b2x=−2τn−αc1τnxRκ3+2α2σxRκ4x+4α2σκ5x2α2στnxRxR+xL,b2y=−2τn−βc2τnyRκ3+2β2σyRκ4y+4β2σκ5y2β2στnyRyR+yL,d1y=κ1yyR+yL−2κ2yyLyR+yL+−2τn+βc2τnyLκ3−2β2σyLκ4y+4β2σκ5y2β2στnyLyR+yL−2κ7xxRxLyR+yL+4κ8xRxLyLyR+yL,d1x=κ1xxR+xL−2κ2xxLxR+xL+−2τn+αc1τnxLκ3−2α2σxLκ4x+4α2σκ5x2α2στnxLxR+xL−2κ7yxR+xLyRyL+4κ8xLxR+xLyRyL,d0=2κ2xxRxL+2κ2yyRyL+κ3τn−2α2σκ5x−τnκ3α2στnxRxL−2β2σκ5y−τnκ3β2στnyRyL−4κ8xRxLyRyL,d2x=κ1xxR+xL−2κ2xxRxR+xL−2τn−αc1τnxRκ3+2α2σxRκ4x+4α2σκ5x2α2στnxRxR+xL+2κ7yxR+xLyRyL+4κ8xRxR+xLyRyL,d2y=−κ1yyR+yL−2κ2yyRyR+yL+−2τn−βc2τnyRκ3+2β2σyRκ4y+4β2σκ5y2β2στnyRyR+yL+2κ7xxRxLyR+yL+4κ8xRxLyRyR+yL,d31=−κ6xR+xLyR+yL+2κ7xxLxR+xLyR+yL+2κ7yxR+xLyLyR+yL−4κ8xLxR+xLyLyR+yL,d32=κ6xR+xLyR+yL+2κ7xxRxR+xLyR+yL−2κ7yxR+xLyLyR+yL−4κ8xRxR+xLyLyR+yL,d33=κ6xR+xLyR+yL−2κ7xxLxR+xLyR+yL+2κ7yxR+xLyRyR+yL−4κ8xLxR+xLyRyR+yL,d34=−κ6xR+xLyR+yL−2κ7xxRxR+xLyR+yL−2κ7yxR+xLyRyR+yL−4κ8xRxR+xLyRyR+yL,s1y=−κ9yyR+yL+2κ10yyLyR+yL,s1x=−κ9xxR+xL+2κ10xxLxR+xL,s0=3κ3−44σ−2κ10xxRxL−2κ10yyRyL,s2x=κ9xxR+xL+2κ10xxRxR+xL,s2y=κ9yyR+yL+2κ10yyRyR+yL,M1=−κ3σ,M2=κ34σ.Through the prior deducing procedure, it can be seen that the truncation error of the Eq. (15) isOτn2+τnxR2−xRxL+xL2+τnyR2−yRyL+yL2+xR5+xL5/xR+xL+yR5+yL5/yR+yL. When h=xR=xL=yR=yL and τn=h2, it owns spatial accuracy of fourth order and temporal accuracy of second order.


## 3. 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(17)τnτn−1=max0<k<K0<j<J1+∂2u/∂t2xk,yj,tn−1/22max0<k<K0<j<J1+∂2u/∂t2xk,yj,tn+1/22,n=1,2,⋯The standard uniform central difference formula is substituted for the∂2u/∂t2 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).


### 3.2. Adaptive Grid Structure in Space

Let a spatial adaptation algorithm be deduced as follows.Wxx,t is taken as the monitor function in x−direction and Wyy,t is taken as the monitor function in y−direction.(18)Wxx,t=1+∑k=0K∂2u/∂t∂xx,t2,0≤x≤1,0≤t<T,(19)Wyy,t=1+∑j=0J∂2u/∂t∂yy,t2,0≤y≤1,0≤t<T.ξkxk,tk=0K and ζjyj,tj=0J point to the set of 2D original spatial mesh points in x−direction and y−direction. The 2D computational area is transformed as xkξk,tk=0K in x−direction and yjζj,tj=0J in y−direction through mesh movement. The continuous function Wxx,t>0 and Wyy,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 x−direction(20)Wxxk+1/2m,tnxk+1m+1−xkm+1−Wxxk−1/2m,tnxkm+1−xk−1m+1=0,k=1,2,⋯,K−1,n≥0.Similarly, we can get the semi-implicit scheme iny−direction(21)Wyyj+1/2m,tnyj+1m+1−yjm+1−Wyyj−1/2m,tnyjm+1−yj−1m+1=0,j=1,2,⋯,J−1,n≥0.


### 3.3. Iterative Adaption Algorithm

This paragraph shows the iterative process of the presented method in the paper.Step 1.
Depending on the spatial meshxkmk=0K, yjmj=0J at the nth time layer, we can obtain the corresponding monitor functions Wxx,t and Wyy,t via Eqs. (18) and (19).Step 2.
The point setxkm+1k=0K refreshed is extracted from the prior grid xkmk=0K via Eq. (20) whereas yjm+1j=0J refreshed is extracted from yjmj=0J by Eq. (21). We use the new point sets to get the replacement for the old point sets iteratively once ∑k=1K−1xkm+1−xkm<er, ∑j=1J−1yjm+1−yjm<er, which er takes 10-5. The last value of m+1 takes Mn.Step 3.
By virtue of the last grid point setsxkMnk=0K, yjMnj=0J, and area ratio algorithm, the solutions in the nth time line uk,jnk=1,j=1K−1,J−1 are all estimated. Combining with τn obtained from Eq. (17), uk,jn+1k=1,j=1K−1,J−1 in the n+1th time line can be calculated by Eq. (15).Step 4.
Repeat steps from the first to the third till quenching occurs, or the solution converges to a steady solution.


## 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(17)τnτn−1=max0<k<K0<j<J1+∂2u/∂t2xk,yj,tn−1/22max0<k<K0<j<J1+∂2u/∂t2xk,yj,tn+1/22,n=1,2,⋯The standard uniform central difference formula is substituted for the∂2u/∂t2 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).


## 3.2. Adaptive Grid Structure in Space

Let a spatial adaptation algorithm be deduced as follows.Wxx,t is taken as the monitor function in x−direction and Wyy,t is taken as the monitor function in y−direction.(18)Wxx,t=1+∑k=0K∂2u/∂t∂xx,t2,0≤x≤1,0≤t<T,(19)Wyy,t=1+∑j=0J∂2u/∂t∂yy,t2,0≤y≤1,0≤t<T.ξkxk,tk=0K and ζjyj,tj=0J point to the set of 2D original spatial mesh points in x−direction and y−direction. The 2D computational area is transformed as xkξk,tk=0K in x−direction and yjζj,tj=0J in y−direction through mesh movement. The continuous function Wxx,t>0 and Wyy,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 x−direction(20)Wxxk+1/2m,tnxk+1m+1−xkm+1−Wxxk−1/2m,tnxkm+1−xk−1m+1=0,k=1,2,⋯,K−1,n≥0.Similarly, we can get the semi-implicit scheme iny−direction(21)Wyyj+1/2m,tnyj+1m+1−yjm+1−Wyyj−1/2m,tnyjm+1−yj−1m+1=0,j=1,2,⋯,J−1,n≥0.


## 3.3. Iterative Adaption Algorithm

This paragraph shows the iterative process of the presented method in the paper.Step 1.
Depending on the spatial meshxkmk=0K, yjmj=0J at the nth time layer, we can obtain the corresponding monitor functions Wxx,t and Wyy,t via Eqs. (18) and (19).Step 2.
The point setxkm+1k=0K refreshed is extracted from the prior grid xkmk=0K via Eq. (20) whereas yjm+1j=0J refreshed is extracted from yjmj=0J by Eq. (21). We use the new point sets to get the replacement for the old point sets iteratively once ∑k=1K−1xkm+1−xkm<er, ∑j=1J−1yjm+1−yjm<er, which er takes 10-5. The last value of m+1 takes Mn.Step 3.
By virtue of the last grid point setsxkMnk=0K, yjMnj=0J, and area ratio algorithm, the solutions in the nth time line uk,jnk=1,j=1K−1,J−1 are all estimated. Combining with τn obtained from Eq. (17), uk,jn+1k=1,j=1K−1,J−1 in the n+1th time line can be calculated by Eq. (15).Step 4.
Repeat steps from the first to the third till quenching occurs, or the solution converges to a steady solution.


## 4. Simulation Demonstration

### 4.1. Numerical Case

This is a common numerical sample utilized to evaluate the performance of Eq. (15) measure by Eq. (22)(22)put=uxx+uyy−pux−puy,x,y∈0,1×0,1,t≥0,(23)ux,y,t=epx−1ep−1+epy−1ep−1+∑q=1∞−1qqπqπ2+p2/4epx−1/2+py−1/2sinqπx+qπye−qπ2/p+p/4t.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: k/K+λx/πsinπk/Kk=0K in x−direction and j/J+λy/πsinπj/Jj=0J in y−direction, where KandJ 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.Table 1
The parameter values for common example.
SchemesKJMax. errorAver. errorCPU timeConv. rateThe scheme in Ref. [26] (uniform)162.59×10−1—1000—325.65×10−2—13442.30645.50×10−3—67663.441283.15×10−4613284.17The scheme in Ref. [26] (nonuniform)161.12×10−1—844—321.22×10−2—16253.34646.84×10−4—111884.251284.42×10−5—1353444.00The scheme in Ref. [27] (uniform)162.4859×10−11.7858×10−247—325.7032×10−21.9148×10−31725.0033645.5382×10−31.4565×10−47023.71661283.1469×10−49.6850×10−629483.91062561.9710×10−56.1742×10−7117943.9714The scheme in Ref. [27] (nonuniform)168.1898×10−38.7725×10−4358—322.9724×10−42.9996×10−528557.9861641.5573×10−51.6521×10−6149144.18241289.3004×10−71.0224×10−7350064.19642561.3364×10−76.4492×10−91659694.7029The proposed scheme (uniform)162.5865×10−11.7571×10−278—325.6511×10−21.9107×10−32974.9902645.5041×10−31.4562×10−412173.71381283.1447×10−49.6804×10−651483.91102561.9357×10−55.9061×10−7209514.0348The proposed scheme (nonuniform)168.1405×10−38.7581×10−4421—322.9716×10−42.9989×10−518577.9849641.5491×10−51.6490×10−669114.18471286.9033×10−78.9949×10−8245854.19642564.1212×10−83.4536×10−91312434.7029With 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 KJ rises. When KJ 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] (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 KJ 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 KJ is 32, 64, 128, and 256, and the tendency becomes more obviously as KJ 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.


### 4.2. Quenching Case without Convection Term

Whenc1x,y=0 and c2x,y=0, with the initial and boundary conditions Eqs. (6) and (7), the original Eq. (5) can be described as(24)α2x2+β2y2q/2ut=1α2uxx+1β2uyy+11−uθ,x,y∈0,1×0,1,t∈0,T.For the situation ofc1x,y=0 and c2x,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=β=10, 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 h0=1/90.Table2 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=β=10, and θ=1 without convection term. The criteria in Table 2 are explained as follows. xmax,ymax refers to the point of quenching location, Tmax refers to the quenching time, maxu refers to the maximum temperature value immediately before quenching happens, and maxut 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 maxu. Tmax of the proposed scheme is 0.5990174307488624, and maxu of the proposed scheme is 0.9901182463562488. But there are greater differences among the three schemes for maxut. maxut 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 ut with regard to spatial variables x, y when the time t is 0.5990174307488624, respectively. From the two 3D surfaces at the time Tmax=0.5990174307488624 that is immediately before quenching occurs, it is found that u⟶1− and ut⟶84− at the quenching location xmax,ymax=0.5,0.5 and the quenching time Tmax=0.5990174307488624.Table 2
Quenching performances of different schemes forq=0, a=β=10, and θ=1 without convection term.
SchemesxmaxymaxTmaxmaxumaxutThe proposed scheme0.5000.5000.59901743074886240.990118246356248883.20947696438735The scheme in [19]——0.587124999937510.990432148.887767The scheme in [20]0.5000.5000.5875540.9992631249.917563Figure 1
(a) The 3D plots ofu immediately before quenching occurs, and (b) the 3D plots of ut immediately before quenching occurs. The parameters are q=0, a=β=10, and θ=1 without convection term.
(a)(b)


### 4.3. Quenching Case with Convection Termb/1+αx+βy

We will investigate what role degeneration parameterq, 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 c1x,y=bx/1+αx+βy and c2x,y=by/1+αx+βy, where bx and by take the constants and convection parameter b=bx=by. When singularity parameter θ is larger than 1, the quenching situation is complex. Therefore, we only consider quenching cases of θ=1. In numerical samples, bx and by 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 h0=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,bx=by=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, ut, x, y, and t.Table 3
Quenching data ofq based on convection term b/1+αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 5.00.042857142857142860.042857142857142863.4284048736271370.9920620807720361125.6265973017968Case 5.1.10. 10.10.64708426291453790.9933649818060033136.4723853886283Case 5.1.20.085714285714285720.085714285714285720.78248470466074350.990381971863921587.67058865539519Case 5.1.30.071428571428571430.071428571428571431.1099426794461640.990847946322951286.47536110256003Case 5.1.40.057142857142857140.057142857142857141.5728604227346150.990278676416439583.87879039070892Case 5.1.50.057142857142857140.057142857142857141.880980950977220.9920140466370957101.1277413288689Case 5.1.60.057142857142857140.057142857142857142.7719285574653510.992005984296346596.8866420271182Case 5.1.70.042857142857142860.042857142857142864.2931560569842690.9908925785572684107.7822215514385Case 5.1.80.042857142857142860.042857142857142865.4853180528869540.9920945867539992126.7837317275917Case 5.1.90.042857142857142860.042857142857142867.1822005120597870.9903718441627313102.1421834048094Case 5.1.100.042857142857142860.042857142857142869.6869022883840930.992048580687276127.8241083054763Table 4
Quenching data ofbx and by based on convection term b/1+αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 5.2.10.17142857142857140.17142857142857148.3790500807861150.99023149627409727.01166111192879Case 5.2.20.12857142857142860.12857142857142866.1446619796208010.990208448695609535.54184621151833Case 5.2.30.085714285714285720.085714285714285724.1332625746227040.990226572804195452.5010571873496Case 5.2.40.057142857142857140.057142857142857143.5406617685147280.990104046454050575.11254034765163Case 5.2.50.042857142857142860.042857142857142863.3893285410551850.992283726259128128.3227508596781Case 5.2.60.042857142857142860.042857142857142863.4494018226441990.9910523423231372108.73967781816Case 5.2.70.028571428571428570.028571428571428573.5783760580867370.990585829058986147.5334079485446Case 5.2.80.028571428571428570.028571428571428574.0188799939399970.9940259443236479246.3627497716133Case 5.2.90.028571428571428570.028571428571428574.3411853450993080.9904731250696808143.5523827878386Case 5.2.100.028571428571428570.028571428571428574.7312223067358770.9969319295625433506.9119438251631Table 5
Quenching data ofα and β based on convection term b/1+αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 5.3.10.200.203.6216245752042090.99171870750864286.83309533401254Case 5.3.20.028571428571428570.028571428571428573.6091533686309490.990718567435303273.0239115821952Case 5.3.30.014285714285714290.014285714285714293.6093906012492290.9912265622789864150.6673183543721Case 5.3.40.014285714285714290.014285714285714293.6092749482193420.9919688690187356165.4554013002607Case 5.3.50.014285714285714290.014285714285714293.3761273432236720.990203949879904790.20671661807094Case 5.3.60.014285714285714290.014285714285714293.3786392565364710.9927000295904764120.3354981361145Case 5.3.70.014285714285714290.014285714285714293.3786681976108720.9923048003313988113.9606848017258Case 5.3.80.014285714285714290.0142857142857142923.381766442034670.990078138653823117.60457567575462Case 5.3.90.014285714285714290.0142857142857142931.09319178344070.990001353838681214.90527975955028Case 5.3.100.014285714285714290.0142857142857142939.496519846580890.99000331782852813.07949247781412Table 6
The maximal values and locations ofu and ut for q=1,bx=by=3,α=β=15 based on convection term b/1+αx+βy.
txymaxx,yumaxx,yut3.4282904954312970.042857142857142860.042857142857142860.98313444981542653.07038382351083.4283426458177350.042857142857142860.042857142857142860.86337297141044968.080632204162383.4283797879896200.042857142857142860.042857142857142860.989320635743985890.305031525598133.4284048736271370.042857142857142860.042857142857142860.9920620807720361125.62659730179683.428420657560073——17.074462901898561.200141619871376×109Except Case 5.0, there are ten quenching cases chosen in Table3. The ten items list as follows: Case 5.1.1 is q=0.1,bx=by=3,α=β=15; Case 5.1.2 is q=0.2,bx=by=3,α=β=15; Case 5.1.3 is q=0.4,bx=by=3,α=β=15; Case 5.1.4 is q=0.6,bx=by=3,α=β=15; Case 5.1.5 is q=0.7,bx=by=3,α=β=15; Case 5.1.6 is q=0.9,bx=by=3,α=β=15; Case1.1.7 is q=1.1,bx=by=3,α=β=15; Case1.1.8 is q=1.2,bx=by=3,α=β=15; Case 5.1.9 is q=1.3,bx=by=3,α=β=15; Case 5.1.10 is q=1.4,bx=by=3,α=β=15. 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 xmax,ymax 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 xmax and ymax do decline as q evolves in the first subdomain and converge to 0.04285714285714286,0.04285714285714286 in the second subdomain. The quenching time Tmax reaches the maximum 9.686902288384093 when q=1.4. Tmax and qare in positive proportion when q in 0.1,1.4.We choose ten items written in Table4 to display the quenching characteristics related to the parameters bx and by. Case 5.2.1 is q=1,bx=by=−5,α=β=15; Case 5.2.2 is q=1,bx=by=−3,α=β=15; Case 5.2.3 isq=1,bx=by=0,α=β=15; Case 5.2.4 is q=1,bx=by=2,α=β=15; Case 5.2.5 is q=1,bx=by=4,α=β=15; Case 5.2.6 is q=1,bx=by=5,α=β=15; Case 5.2.7 is q=1,bx=by=6,α=β=15; Case 5.2.8 is q=1,bx=by=8,α=β=15; Case 5.2.9 is q=1,bx=by=9,α=β=15; Case 5.2.10 is q=1,bx=by=10,α=β=15; bx and by in Eq. (5) are theoretically equivalent, i.e., b=bx=by. A valid scope of b is −5,10 for these cases. Because quenching occurs fastest when b takes 4, b=4 is thought as a special point for quenching time. Quenching time is inversely proportional to b when b∈−5,4 and positively proportional to b when b∈4,10. Quenching location xmax,ymax rises as b increases when b∈−5,5 and takes 0.02857142857142857,0.02857142857142857 when b≥6.Table5 offers seventeen referencing cases. The first ten items are Case 5.3.1 is q=1,bx=by=3,α/β=1,α2=15; Case 5.3.2 is q=1,bx=by=3,α/β=1,α2=1227; Case 5.3.3 is q=1,bx=by=3,α/β=1,α2=1228; Case 5.3.4 is q=1,bx=by=3,α/β=1,α2=1229; Case 5.3.5 is q=1,bx=by=3,α/β=1,α2=2900; Case 5.3.6 is q=1,bx=by=3,α/β=1,α2=2999; Case 5.3.7 is q=1,bx=by=3,α/β=1,α2=3000; Case 5.3.8 is q=1,bx=by=3,α/β=1,α2=8×104; Case 5.3.9 is q=1,bx=by=3,α/β=1,α2=1.1×105; Case 5.3.10 is q=1,bx=by=3,α/β=1,α2=1.43×105. We configure 15≤α2≤1.43×105 and α/β=1 for the prior ten items. Experiments show even if α2 reaches 105, quenching phenomena still happen. Furthermore, if α2 falls in the interval 15,1.43×105, then quenching phenomena will occur normally, and quenching location xmax and ymax do not decrease monotonously as α2 declines. When α2 lies in the area of 1228,1.43×105, xmax,ymax always keeps 0.01428571428571429,0.01428571428571429. In terms of quenching time, there is no obvious linear relationship between Tmax and α2 in a small range of the definition field of α2. If the measurement scale is enlarged in the domain of α2, especially when α2 reaches 2900, Tmax will rise with the increase of α2, i.e., Tmax is 3.376127343223672 when α2 is 2900 and Tmax is 39.49651984658089 when α2 is 1.43×105.We take Case 5.0 to demonstrate specific quenching phenomena. In this case, quenching occurs atT5∗=3.428404873627137, x5∗=0.04285714285714286, and y5∗=0.04285714285714286 with the parameters q=1,bx=by=3,α=β=15. The next context declares more details about quenching. Let us observe the performances of the solution u, the derivative function ut, and their time and space locations of five stages relative to quenching for Case 5.0, which show in 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, ut blows up and reaches 1.200141619871376×109, the position of maxx,yu is (0.5336298640841476×10‐5, 0.2902651144597507), and the position of maxx,yut is (0.5702149685882448×10‐4, 0.5677032305886261×10‐4) when t=3.428420657560073. Figures 2(a) and 2(b) denote the profiles of both solution u and derivative function ut. In the aforementioned figures, the maximum values of u and ut grow from 0 to T5∗. As soon as time variable t arrives at quenching spot x5∗,y5∗ and quenching time T5∗, we have maxu=0.99206208077203618 and maxut=125.6265973017968. By locating the quenching time T5∗, we can draw the plots of maxx,yut and maxx,yutt. Figures 3(a)–3(d) depict the function relationships between u and x, between ut and x, between u and y, and between ut and y, respectively, when t isT5∗. Owing to the identity of x and y in the original equation, Figures 3(a) and 3(c) have the same shape, and Figures 3(b) and 3(d) have the same shape.Figure 2
The graphs of solutionu and ut: (a) maxx,yu as t increases, and (b) maxx,yut as t increases from the initial moment to the occurrence of quenching behavior with the parameters q=1,bx=by=3,α=β=15 based on convection term b/1+αx+βy.
(a)(b)Figure 3
The graphs of (a)u as x increases, (b) ut as x increases, (c) u as yincreases, and (d) ut as y increases when t is T5∗ for convection term b/1+αx+βy.
(a)(b)(c)(d)Figure 4
(a) The variation of spatial stephk as x processes; (b) the variation of spatial step hj as y processes; (c) the curve of the temporal steps τn as time t increases with the parameters q=1,bx=by=3,α=β=15 based on convection term b/1+αx+βy.
(a)(b)(c)Figure4 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 x5∗(y5∗). 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⟶T5∗ 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 T5∗ and lasts through the rest course of calculation. So this quenching phenomenon is caught with the features of maxx,yu⟶1− and maxx,yut⟶126− as t reaches T5∗.There are three groups of 3D profiles of the solutions and their temporal derivatives in Figure5. Figures 5(a) and 5(b) denote the two three-dimensional plots of u and ut when t=0.40, respectively. Figures 5(c) and 5(d) denote the two three-dimensional plots of uandutwhent=3.428342645817735, respectively. Figures 5(e) and 5(f) denote the two three-dimensional plots of u and ut when t=T5∗, respectively. Thought the 3D plots, we can get richer information. The first group is the three-dimensional views from both u and ut at the first 500th temporal layers in the excursion, in which maxx,yu⟶0.2215 and maxx,yut⟶12.9459 when t=0.40. The second group is the three-dimensional views from both u and ut including the penultimate temporal steps before T5∗, in whichmaxx,yu⟶0.9863 and maxx,yut⟶68.0806 when t=3.428342645817735. The third group is the three-dimensional views from both u and ut immediately before quenching happens, in which maxx,yu⟶0.9921 and maxx,yut⟶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 T5∗. 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 T5∗, the maximum ut also increases rapidly and infinitely and then blows up at the next time layer of quenching time.Figure 5
(a) The 3D plot ofu when t=0.40; (b) the 3D plot of ut when t=0.40; (c) the 3D plot of u when t=3.428342645817735; (d) the 3D plot of ut whent=3.428342645817735; (e) the 3D plot of u when t=T5∗; (f) the 3D plot of ut when t=T5∗. All is for Case 5.0.
(a)(b)(c)(d)(e)(f)


### 4.4. Quenching Example with Convection Termb/α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+βywhich is related to the three groups of elements: q and bx and by, α, and β and exemplify the quenching behaviors of the representative Case 6.0. The convection term b/αx+βy should written as c1x,y=bx/αx+βy and c2x,y=by/αx+βy, where bx and by take the constants and convection parameter b=bx=by. For Case 6.0, the initial x−step is h0=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 and defined as q=0.8,bx=by=−2,α=β=11000. Tables 7–10 and Figures 6–9 demonstrate specific quenching information of Case 6.0.Table 7
Quenching data ofq based on convection term b/αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 6.00.01250.01254.0204492345571910.98118665989247131.34049546140411Case 6.1.10.0250.03750.48934004717901080.981554723398659553.96815295250957Case 6.1.20.0250.0251.405465528266710.981067316599417530.97740120847389Case 6.1.30.0250.01251.8068560066123640.983528504653534334.84710326333095Case 6.1.40.01250.01252.3158009120378780.986003207600720548.09280080596852Case 6.1.50.01250.01253.0042833021985410.984760691334703341.42159866823076Case 6.1.60.0250.0255.5195208961405880.980491616477197715.56936805077945Case 6.1.70.0250.0257.61297811818480.980609367578944713.71511944823688Case 6.1.80.0250.02510.863156909568730.980366871397270111.86327696425435Case 6.1.90.0250.02516.290591597143860.980071806660898610.23980977162735Case 6.1.100.0250.02526.338695532458330.98021290675467919.042418195146764Case 6.1.110.0250.02547.981168911737580.9802932569090137.964013508962916Table 8
Quenching data ofbx and by based on convection term b/αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 6.2.10.06250.062510.990932677867840.98044984722915838.583794969514605Case 6.2.20.03750.03757.8096557019533230.980450319388935112.7219405872281Case 6.2.30.03750.03756.2105342225412080.980912286364063413.14389285559599Case 6.2.40.0250.0254.9560041056307880.981837226516544218.9029625958216Case 6.2.50.01250.01252.647996646301790.981096073619703331.29445994657764Case 6.2.60.01250.01259.5596566535946850.981300410148895833.42529207679014Case 6.2.70.0250.02518.16226159738410.980898541828047618.03120629119135Case 6.2.80.03750.02538.775361395170310.980409034462982414.57573131335622Case 6.2.90.03750.037559.072441770307950.980250333413378812.64186937140285Case 6.2.100.03750.037589.056784028483560.980353824767113313.00866087970327Table 9
Quenching data ofα and β based on convection term b/αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 6.3.10.53750.53758.2400227397661650.980523162951724419.33141003501417Case 6.3.20.12290.12296.2407929288174660.980442975903644821.02610806687704Case 6.3.30.06250.06253.4313343512955050.980941810169855521.77295829105509Case 6.3.40.03750.03753.4417379588903710.981637200195769623.0950930660414Case 6.3.50.0250.0253.7505512992212610.980132495794395719.04800266218939Case 6.3.60.0250.0254.0241522731641590.980042038117032517.6646128655794Case 6.3.70.0250.0254.0651551384233620.980594420178104318.01149509967014Case 6.3.80.01250.01254.065585912787060.980706513470496230.77480165240926Case 6.3.90.01250.01253.5201436052809370.980920890316600624.6618068728011Case 6.3.100.01250.012513.442434023300710.98055758605504468.352004051122194Table 10
The maximal values and locations ofu and ut for q=0.8,bx=by=−2,α=β=11000 based on convection term b/αx+βy.
txymaxx,yumaxx,yut3.9469999999996760.01250.01250.81959863941021882.4991221910965863.9999999999996710.01250.01250.89217197319210934.6842992690805284.0169999999996760.01250.01250.948868692550773410.883646217714824.0195956827666690.01250.01250.969284409499635618.78802594909864.0204492345571910.01250.01250.98118665989247131.340495461404114.02068525061823——5.358959032595640×10167.966445413263532×1019Figure 6
The graphs of solutionu and ut: (a) maxx,yu as t increases, and (b) maxx,yut as t increases from the initial time to the occurrence of quench behavior with q=0.8,bx=by=−2,α=β=11000 based on convection term b/αx+βy.
(a)(b)Figure 7
The graphs of (a)u as x increases, (b) ut as x increases, (c) u as y increases, and (d) ut as y increases when t is T5∗ for convection term b/1+αx+βy.
(a)(b)(c)(d)Figure 8
(a) The variation of spatial stephk as space x processes; (b) the variation of spatial step hj as space y processes; (c) the curve of the temporal steps τn as time t increases with the parameters q=0.8,bx=by=−2,α=β=11000 based on convection term b/αx+βy.
(a)(b)(c)Figure 9
(a) The 3D plot ofu when t=T6∗ and (b) the 3D plot of ut when t=T6∗ for Case 6.0.
(a)(b)The effect ofq 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 xmax is not equal to ymax when q takes 0, 0.5 while xmax is equal to ymax in the other cases; xmax,ymax is 0.0125,0.0125 when q takes 0.6, 0.7, and 0.8; xmax,ymax is 0.025,0.025 when q is in 0.9,1.4. Quenching time Tmax increases as q increases in the domain of q. Table 7 gives some of these quenching cases, in which Case 6.1.1 isq=0,bx=by=−2,α=β=11000; Case 6.1.2 is q=0.4,bx=by=−2,α=β=11000; Case 6.1.3 is q=0.5,bx=by=−2,α=β=11000; Case 6.1.4 is q=0.6,bx=by=−2,α=β=11000; Case 6.1.5 is q=0.7,bx=by=−2,α=β=11000; Case 6.1.6 is q=0.9,bx=by=−2,α=β=11000; Case 6.1.7 is q=1,bx=by=−2,α=β=11000; Case 6.1.8 is q=1.1,bx=by=−2,α=β=11000; Case 6.1.9 is q=1.2,bx=by=−2,α=β=11000; Case 6.1.10 is q=1.3,bx=by=−2,α=β=11000; Case 6.1.11 is q=1.4,bx=by=−2,α=β=11000; Case 6.0 is between Case 6.1.5 and Case 6.1.6.The paragraphs describe the relationship betweenbx(by) and quenching features. We write ten cases in Table 8. Case 6.2.1 is q=0.8,bx=by=−16,α=β=11000; Case 6.2.2 is q=0.8,bx=by=−10,α=β=11000; Case 6.2.3 is q=0.8,bx=by=−7,α=β=11000; Case 6.2.4 is q=0.8,bx=by=−5,α=β=11000; Case 6.2.5 is q=0.8,bx=by=1,α=β=11000; Case 6.2.6 is q=0.8,bx=by=12,α=β=11000; Case 6.2.7 is q=0.8,bx=by=21,α=β=11000; Case 6.2.8 is q=0.8,bx=by=40,α=β=11000; Case 6.2.9 is q=0.8,bx=by=57,α=β=11000; Case 6.2.10 is q=0.8,bx=by=80,α=β=11000; Case 6.0 is between Case 6.2.4 and Case 6.2.5. Quenching behavior can be formed when bx(by) takes the values in −16,80, in which the function curve of (b,Tmax) is concave and Tmax is smallest when b=1. When b=−16, b=−13, b∈−10,−7, b∈−5,12, b∈21,39, or b∈57,80, quenching location xmax is equal to ymax whereas when b takes other values in −16,80, xmax is different from ymax.We describe the influence ofα and β on quenching results just depending on quenching data in Table 9. Of course, we do a great quantity of experiments, from which some typical cases chosen in Table 9. We set α/β=1 and α2from 10 to 3.0×106 for Cases 6.3.1-6.3.10. Specifically, Case 6.3.1 is q=0.8,bx=by=−2,α/β=1,α2=10; Case 6.3.2 is q=0.8,bx=by=−2,α/β=1,α2=231; Case 6.3.3 is q=0.8,bx=by=−2,α/β=1,α2=1000; Case 6.3.4 is q=0.8,bx=by=−2,α/β=1,α2=2756; Case 6.3.5 is q=0.8,bx=by=−2,α/β=1,α2=8648; Case 6.3.6 is q=0.8,bx=by=−2,α/β=1,α2=10500; Case 6.3.7 is q=0.8,bx=by=−2,α/β=1,α2=10775; Case 6.3.8 is q=0.8,bx=by=−2,α/β=1,α2=10779; Case 6.3.9 is q=0.8,bx=by=−2,α/β=1,α2=19991; Case 6.3.10 is q=0.8,bx=by=−2,α/β=1,α2=300000; Case 6.0 is between Case 6.3.8 and Case 6.3.9.It does not form quenching status whenα2<10 or α2<3.0×106. After α/β is defined as 1, α2 is evaluated from 10 to 3.0×106. We find that it is more likely to produce quenching behaviors when α2 is between 10 and 3.0×106. By means of continual tests, Cases 6.3.1-6.3.10 represent 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×106. 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, Tmax increases with the increase of α2.Quenching phenomena of Case 6.0 can be considered in the following text. Its quenching time appears atT6∗=4.020449234557191, and its quenching position appears at x6∗,y6∗=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. maxx,yu is 5.358959032595640×1016 at (2.0588245639243×10−4, 9.654095827421034×10−1), and maxx,yut is 7.966445413263532×1019 at (2.080224198836234×10−5, 9.824584434865121×10−1) when quenching occurs. A comprehensive statement of quenching states with the parameters q=0.8,bx=by=−2,α=β=11000 is recorded in the next paragraphs. Although there is slight perturbation in the left boundary of ut, it does not influences on distribution of ut. It is evident that Figures 6(a) and 6(b) give two pairs curves of both maxx,yu and maxx,yut. The first one is for the distribution of the solution maxx,yut as t varies, and the second one is for the distribution of maxx,yutt as t increases. There is small perturbation near the start-up of the adaptive procedure for the left figure of maxx,yut. The four figures in Figure 7 give function relationship physical quantities betweenu, ut, and spatial variables. Figure 7(a) paints the contour of (x,u), Figure 7(b) paints the contour of (x,ut), Figure 7(c) paints the contour of (y,u), and Figure 7(d) paints the contour of (y,ut). Similarly, Figures 7(a) and 7(c) have the same shape, and Figures 7(b) and 7(d) have the same shape.There are three graphics in Figure8, which represent three function relationships between spatial step in x−direction and x, between spatial step in y−direction and y, and between temporal step and t, respectively. Figure 8(a) depicts two curves, in which one is between initial spatial steps hk and x, and another is between quenching spatial steps hk and x. Figure 8(b) also describes two curves, in which one is between initial spatial steps hk and y, and another is between quenching spatial steps hk 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τnvaries as t changes from 0 to T6∗. τ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=T6∗.The three-dimensional surfaces can render more reliable information for both the solutionu and the time derivative ut which display in Figure 9. The group of surfaces is the three-dimensional views of both u and ut when t is equal to T6∗. 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 x6∗,y6∗ is at 0.0125,0.0125 immediately before quenching occurs, the peak value of the solution approaches 1. The three-dimensional 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 maxx,yu⟶1− and maxx,yut⟶32− as t reaches T6∗.


## 4.1. Numerical Case

This is a common numerical sample utilized to evaluate the performance of Eq. (15) measure by Eq. (22)(22)put=uxx+uyy−pux−puy,x,y∈0,1×0,1,t≥0,(23)ux,y,t=epx−1ep−1+epy−1ep−1+∑q=1∞−1qqπqπ2+p2/4epx−1/2+py−1/2sinqπx+qπye−qπ2/p+p/4t.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: k/K+λx/πsinπk/Kk=0K in x−direction and j/J+λy/πsinπj/Jj=0J in y−direction, where KandJ 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.Table 1
The parameter values for common example.
SchemesKJMax. errorAver. errorCPU timeConv. rateThe scheme in Ref. [26] (uniform)162.59×10−1—1000—325.65×10−2—13442.30645.50×10−3—67663.441283.15×10−4613284.17The scheme in Ref. [26] (nonuniform)161.12×10−1—844—321.22×10−2—16253.34646.84×10−4—111884.251284.42×10−5—1353444.00The scheme in Ref. [27] (uniform)162.4859×10−11.7858×10−247—325.7032×10−21.9148×10−31725.0033645.5382×10−31.4565×10−47023.71661283.1469×10−49.6850×10−629483.91062561.9710×10−56.1742×10−7117943.9714The scheme in Ref. [27] (nonuniform)168.1898×10−38.7725×10−4358—322.9724×10−42.9996×10−528557.9861641.5573×10−51.6521×10−6149144.18241289.3004×10−71.0224×10−7350064.19642561.3364×10−76.4492×10−91659694.7029The proposed scheme (uniform)162.5865×10−11.7571×10−278—325.6511×10−21.9107×10−32974.9902645.5041×10−31.4562×10−412173.71381283.1447×10−49.6804×10−651483.91102561.9357×10−55.9061×10−7209514.0348The proposed scheme (nonuniform)168.1405×10−38.7581×10−4421—322.9716×10−42.9989×10−518577.9849641.5491×10−51.6490×10−669114.18471286.9033×10−78.9949×10−8245854.19642564.1212×10−83.4536×10−91312434.7029With 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 KJ rises. When KJ 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] (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 KJ 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 KJ is 32, 64, 128, and 256, and the tendency becomes more obviously as KJ 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.


## 4.2. Quenching Case without Convection Term

Whenc1x,y=0 and c2x,y=0, with the initial and boundary conditions Eqs. (6) and (7), the original Eq. (5) can be described as(24)α2x2+β2y2q/2ut=1α2uxx+1β2uyy+11−uθ,x,y∈0,1×0,1,t∈0,T.For the situation ofc1x,y=0 and c2x,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=β=10, 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 h0=1/90.Table2 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=β=10, and θ=1 without convection term. The criteria in Table 2 are explained as follows. xmax,ymax refers to the point of quenching location, Tmax refers to the quenching time, maxu refers to the maximum temperature value immediately before quenching happens, and maxut 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 maxu. Tmax of the proposed scheme is 0.5990174307488624, and maxu of the proposed scheme is 0.9901182463562488. But there are greater differences among the three schemes for maxut. maxut 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 ut with regard to spatial variables x, y when the time t is 0.5990174307488624, respectively. From the two 3D surfaces at the time Tmax=0.5990174307488624 that is immediately before quenching occurs, it is found that u⟶1− and ut⟶84− at the quenching location xmax,ymax=0.5,0.5 and the quenching time Tmax=0.5990174307488624.Table 2
Quenching performances of different schemes forq=0, a=β=10, and θ=1 without convection term.
SchemesxmaxymaxTmaxmaxumaxutThe proposed scheme0.5000.5000.59901743074886240.990118246356248883.20947696438735The scheme in [19]——0.587124999937510.990432148.887767The scheme in [20]0.5000.5000.5875540.9992631249.917563Figure 1
(a) The 3D plots ofu immediately before quenching occurs, and (b) the 3D plots of ut immediately before quenching occurs. The parameters are q=0, a=β=10, and θ=1 without convection term.
(a)(b)


## 4.3. Quenching Case with Convection Termb/1+αx+βy

We will investigate what role degeneration parameterq, 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 c1x,y=bx/1+αx+βy and c2x,y=by/1+αx+βy, where bx and by take the constants and convection parameter b=bx=by. When singularity parameter θ is larger than 1, the quenching situation is complex. Therefore, we only consider quenching cases of θ=1. In numerical samples, bx and by 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 h0=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,bx=by=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, ut, x, y, and t.Table 3
Quenching data ofq based on convection term b/1+αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 5.00.042857142857142860.042857142857142863.4284048736271370.9920620807720361125.6265973017968Case 5.1.10. 10.10.64708426291453790.9933649818060033136.4723853886283Case 5.1.20.085714285714285720.085714285714285720.78248470466074350.990381971863921587.67058865539519Case 5.1.30.071428571428571430.071428571428571431.1099426794461640.990847946322951286.47536110256003Case 5.1.40.057142857142857140.057142857142857141.5728604227346150.990278676416439583.87879039070892Case 5.1.50.057142857142857140.057142857142857141.880980950977220.9920140466370957101.1277413288689Case 5.1.60.057142857142857140.057142857142857142.7719285574653510.992005984296346596.8866420271182Case 5.1.70.042857142857142860.042857142857142864.2931560569842690.9908925785572684107.7822215514385Case 5.1.80.042857142857142860.042857142857142865.4853180528869540.9920945867539992126.7837317275917Case 5.1.90.042857142857142860.042857142857142867.1822005120597870.9903718441627313102.1421834048094Case 5.1.100.042857142857142860.042857142857142869.6869022883840930.992048580687276127.8241083054763Table 4
Quenching data ofbx and by based on convection term b/1+αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 5.2.10.17142857142857140.17142857142857148.3790500807861150.99023149627409727.01166111192879Case 5.2.20.12857142857142860.12857142857142866.1446619796208010.990208448695609535.54184621151833Case 5.2.30.085714285714285720.085714285714285724.1332625746227040.990226572804195452.5010571873496Case 5.2.40.057142857142857140.057142857142857143.5406617685147280.990104046454050575.11254034765163Case 5.2.50.042857142857142860.042857142857142863.3893285410551850.992283726259128128.3227508596781Case 5.2.60.042857142857142860.042857142857142863.4494018226441990.9910523423231372108.73967781816Case 5.2.70.028571428571428570.028571428571428573.5783760580867370.990585829058986147.5334079485446Case 5.2.80.028571428571428570.028571428571428574.0188799939399970.9940259443236479246.3627497716133Case 5.2.90.028571428571428570.028571428571428574.3411853450993080.9904731250696808143.5523827878386Case 5.2.100.028571428571428570.028571428571428574.7312223067358770.9969319295625433506.9119438251631Table 5
Quenching data ofα and β based on convection term b/1+αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 5.3.10.200.203.6216245752042090.99171870750864286.83309533401254Case 5.3.20.028571428571428570.028571428571428573.6091533686309490.990718567435303273.0239115821952Case 5.3.30.014285714285714290.014285714285714293.6093906012492290.9912265622789864150.6673183543721Case 5.3.40.014285714285714290.014285714285714293.6092749482193420.9919688690187356165.4554013002607Case 5.3.50.014285714285714290.014285714285714293.3761273432236720.990203949879904790.20671661807094Case 5.3.60.014285714285714290.014285714285714293.3786392565364710.9927000295904764120.3354981361145Case 5.3.70.014285714285714290.014285714285714293.3786681976108720.9923048003313988113.9606848017258Case 5.3.80.014285714285714290.0142857142857142923.381766442034670.990078138653823117.60457567575462Case 5.3.90.014285714285714290.0142857142857142931.09319178344070.990001353838681214.90527975955028Case 5.3.100.014285714285714290.0142857142857142939.496519846580890.99000331782852813.07949247781412Table 6
The maximal values and locations ofu and ut for q=1,bx=by=3,α=β=15 based on convection term b/1+αx+βy.
txymaxx,yumaxx,yut3.4282904954312970.042857142857142860.042857142857142860.98313444981542653.07038382351083.4283426458177350.042857142857142860.042857142857142860.86337297141044968.080632204162383.4283797879896200.042857142857142860.042857142857142860.989320635743985890.305031525598133.4284048736271370.042857142857142860.042857142857142860.9920620807720361125.62659730179683.428420657560073——17.074462901898561.200141619871376×109Except Case 5.0, there are ten quenching cases chosen in Table3. The ten items list as follows: Case 5.1.1 is q=0.1,bx=by=3,α=β=15; Case 5.1.2 is q=0.2,bx=by=3,α=β=15; Case 5.1.3 is q=0.4,bx=by=3,α=β=15; Case 5.1.4 is q=0.6,bx=by=3,α=β=15; Case 5.1.5 is q=0.7,bx=by=3,α=β=15; Case 5.1.6 is q=0.9,bx=by=3,α=β=15; Case1.1.7 is q=1.1,bx=by=3,α=β=15; Case1.1.8 is q=1.2,bx=by=3,α=β=15; Case 5.1.9 is q=1.3,bx=by=3,α=β=15; Case 5.1.10 is q=1.4,bx=by=3,α=β=15. 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 xmax,ymax 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 xmax and ymax do decline as q evolves in the first subdomain and converge to 0.04285714285714286,0.04285714285714286 in the second subdomain. The quenching time Tmax reaches the maximum 9.686902288384093 when q=1.4. Tmax and qare in positive proportion when q in 0.1,1.4.We choose ten items written in Table4 to display the quenching characteristics related to the parameters bx and by. Case 5.2.1 is q=1,bx=by=−5,α=β=15; Case 5.2.2 is q=1,bx=by=−3,α=β=15; Case 5.2.3 isq=1,bx=by=0,α=β=15; Case 5.2.4 is q=1,bx=by=2,α=β=15; Case 5.2.5 is q=1,bx=by=4,α=β=15; Case 5.2.6 is q=1,bx=by=5,α=β=15; Case 5.2.7 is q=1,bx=by=6,α=β=15; Case 5.2.8 is q=1,bx=by=8,α=β=15; Case 5.2.9 is q=1,bx=by=9,α=β=15; Case 5.2.10 is q=1,bx=by=10,α=β=15; bx and by in Eq. (5) are theoretically equivalent, i.e., b=bx=by. A valid scope of b is −5,10 for these cases. Because quenching occurs fastest when b takes 4, b=4 is thought as a special point for quenching time. Quenching time is inversely proportional to b when b∈−5,4 and positively proportional to b when b∈4,10. Quenching location xmax,ymax rises as b increases when b∈−5,5 and takes 0.02857142857142857,0.02857142857142857 when b≥6.Table5 offers seventeen referencing cases. The first ten items are Case 5.3.1 is q=1,bx=by=3,α/β=1,α2=15; Case 5.3.2 is q=1,bx=by=3,α/β=1,α2=1227; Case 5.3.3 is q=1,bx=by=3,α/β=1,α2=1228; Case 5.3.4 is q=1,bx=by=3,α/β=1,α2=1229; Case 5.3.5 is q=1,bx=by=3,α/β=1,α2=2900; Case 5.3.6 is q=1,bx=by=3,α/β=1,α2=2999; Case 5.3.7 is q=1,bx=by=3,α/β=1,α2=3000; Case 5.3.8 is q=1,bx=by=3,α/β=1,α2=8×104; Case 5.3.9 is q=1,bx=by=3,α/β=1,α2=1.1×105; Case 5.3.10 is q=1,bx=by=3,α/β=1,α2=1.43×105. We configure 15≤α2≤1.43×105 and α/β=1 for the prior ten items. Experiments show even if α2 reaches 105, quenching phenomena still happen. Furthermore, if α2 falls in the interval 15,1.43×105, then quenching phenomena will occur normally, and quenching location xmax and ymax do not decrease monotonously as α2 declines. When α2 lies in the area of 1228,1.43×105, xmax,ymax always keeps 0.01428571428571429,0.01428571428571429. In terms of quenching time, there is no obvious linear relationship between Tmax and α2 in a small range of the definition field of α2. If the measurement scale is enlarged in the domain of α2, especially when α2 reaches 2900, Tmax will rise with the increase of α2, i.e., Tmax is 3.376127343223672 when α2 is 2900 and Tmax is 39.49651984658089 when α2 is 1.43×105.We take Case 5.0 to demonstrate specific quenching phenomena. In this case, quenching occurs atT5∗=3.428404873627137, x5∗=0.04285714285714286, and y5∗=0.04285714285714286 with the parameters q=1,bx=by=3,α=β=15. The next context declares more details about quenching. Let us observe the performances of the solution u, the derivative function ut, and their time and space locations of five stages relative to quenching for Case 5.0, which show in 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, ut blows up and reaches 1.200141619871376×109, the position of maxx,yu is (0.5336298640841476×10‐5, 0.2902651144597507), and the position of maxx,yut is (0.5702149685882448×10‐4, 0.5677032305886261×10‐4) when t=3.428420657560073. Figures 2(a) and 2(b) denote the profiles of both solution u and derivative function ut. In the aforementioned figures, the maximum values of u and ut grow from 0 to T5∗. As soon as time variable t arrives at quenching spot x5∗,y5∗ and quenching time T5∗, we have maxu=0.99206208077203618 and maxut=125.6265973017968. By locating the quenching time T5∗, we can draw the plots of maxx,yut and maxx,yutt. Figures 3(a)–3(d) depict the function relationships between u and x, between ut and x, between u and y, and between ut and y, respectively, when t isT5∗. Owing to the identity of x and y in the original equation, Figures 3(a) and 3(c) have the same shape, and Figures 3(b) and 3(d) have the same shape.Figure 2
The graphs of solutionu and ut: (a) maxx,yu as t increases, and (b) maxx,yut as t increases from the initial moment to the occurrence of quenching behavior with the parameters q=1,bx=by=3,α=β=15 based on convection term b/1+αx+βy.
(a)(b)Figure 3
The graphs of (a)u as x increases, (b) ut as x increases, (c) u as yincreases, and (d) ut as y increases when t is T5∗ for convection term b/1+αx+βy.
(a)(b)(c)(d)Figure 4
(a) The variation of spatial stephk as x processes; (b) the variation of spatial step hj as y processes; (c) the curve of the temporal steps τn as time t increases with the parameters q=1,bx=by=3,α=β=15 based on convection term b/1+αx+βy.
(a)(b)(c)Figure4 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 x5∗(y5∗). 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⟶T5∗ 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 T5∗ and lasts through the rest course of calculation. So this quenching phenomenon is caught with the features of maxx,yu⟶1− and maxx,yut⟶126− as t reaches T5∗.There are three groups of 3D profiles of the solutions and their temporal derivatives in Figure5. Figures 5(a) and 5(b) denote the two three-dimensional plots of u and ut when t=0.40, respectively. Figures 5(c) and 5(d) denote the two three-dimensional plots of uandutwhent=3.428342645817735, respectively. Figures 5(e) and 5(f) denote the two three-dimensional plots of u and ut when t=T5∗, respectively. Thought the 3D plots, we can get richer information. The first group is the three-dimensional views from both u and ut at the first 500th temporal layers in the excursion, in which maxx,yu⟶0.2215 and maxx,yut⟶12.9459 when t=0.40. The second group is the three-dimensional views from both u and ut including the penultimate temporal steps before T5∗, in whichmaxx,yu⟶0.9863 and maxx,yut⟶68.0806 when t=3.428342645817735. The third group is the three-dimensional views from both u and ut immediately before quenching happens, in which maxx,yu⟶0.9921 and maxx,yut⟶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 T5∗. 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 T5∗, the maximum ut also increases rapidly and infinitely and then blows up at the next time layer of quenching time.Figure 5
(a) The 3D plot ofu when t=0.40; (b) the 3D plot of ut when t=0.40; (c) the 3D plot of u when t=3.428342645817735; (d) the 3D plot of ut whent=3.428342645817735; (e) the 3D plot of u when t=T5∗; (f) the 3D plot of ut when t=T5∗. All is for Case 5.0.
(a)(b)(c)(d)(e)(f)


## 4.4. Quenching Example with Convection Termb/α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+βywhich is related to the three groups of elements: q and bx and by, α, and β and exemplify the quenching behaviors of the representative Case 6.0. The convection term b/αx+βy should written as c1x,y=bx/αx+βy and c2x,y=by/αx+βy, where bx and by take the constants and convection parameter b=bx=by. For Case 6.0, the initial x−step is h0=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 and defined as q=0.8,bx=by=−2,α=β=11000. Tables 7–10 and Figures 6–9 demonstrate specific quenching information of Case 6.0.Table 7
Quenching data ofq based on convection term b/αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 6.00.01250.01254.0204492345571910.98118665989247131.34049546140411Case 6.1.10.0250.03750.48934004717901080.981554723398659553.96815295250957Case 6.1.20.0250.0251.405465528266710.981067316599417530.97740120847389Case 6.1.30.0250.01251.8068560066123640.983528504653534334.84710326333095Case 6.1.40.01250.01252.3158009120378780.986003207600720548.09280080596852Case 6.1.50.01250.01253.0042833021985410.984760691334703341.42159866823076Case 6.1.60.0250.0255.5195208961405880.980491616477197715.56936805077945Case 6.1.70.0250.0257.61297811818480.980609367578944713.71511944823688Case 6.1.80.0250.02510.863156909568730.980366871397270111.86327696425435Case 6.1.90.0250.02516.290591597143860.980071806660898610.23980977162735Case 6.1.100.0250.02526.338695532458330.98021290675467919.042418195146764Case 6.1.110.0250.02547.981168911737580.9802932569090137.964013508962916Table 8
Quenching data ofbx and by based on convection term b/αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 6.2.10.06250.062510.990932677867840.98044984722915838.583794969514605Case 6.2.20.03750.03757.8096557019533230.980450319388935112.7219405872281Case 6.2.30.03750.03756.2105342225412080.980912286364063413.14389285559599Case 6.2.40.0250.0254.9560041056307880.981837226516544218.9029625958216Case 6.2.50.01250.01252.647996646301790.981096073619703331.29445994657764Case 6.2.60.01250.01259.5596566535946850.981300410148895833.42529207679014Case 6.2.70.0250.02518.16226159738410.980898541828047618.03120629119135Case 6.2.80.03750.02538.775361395170310.980409034462982414.57573131335622Case 6.2.90.03750.037559.072441770307950.980250333413378812.64186937140285Case 6.2.100.03750.037589.056784028483560.980353824767113313.00866087970327Table 9
Quenching data ofα and β based on convection term b/αx+βy.
Case namexmaxymaxTmaxmaxumaxutCase 6.3.10.53750.53758.2400227397661650.980523162951724419.33141003501417Case 6.3.20.12290.12296.2407929288174660.980442975903644821.02610806687704Case 6.3.30.06250.06253.4313343512955050.980941810169855521.77295829105509Case 6.3.40.03750.03753.4417379588903710.981637200195769623.0950930660414Case 6.3.50.0250.0253.7505512992212610.980132495794395719.04800266218939Case 6.3.60.0250.0254.0241522731641590.980042038117032517.6646128655794Case 6.3.70.0250.0254.0651551384233620.980594420178104318.01149509967014Case 6.3.80.01250.01254.065585912787060.980706513470496230.77480165240926Case 6.3.90.01250.01253.5201436052809370.980920890316600624.6618068728011Case 6.3.100.01250.012513.442434023300710.98055758605504468.352004051122194Table 10
The maximal values and locations ofu and ut for q=0.8,bx=by=−2,α=β=11000 based on convection term b/αx+βy.
txymaxx,yumaxx,yut3.9469999999996760.01250.01250.81959863941021882.4991221910965863.9999999999996710.01250.01250.89217197319210934.6842992690805284.0169999999996760.01250.01250.948868692550773410.883646217714824.0195956827666690.01250.01250.969284409499635618.78802594909864.0204492345571910.01250.01250.98118665989247131.340495461404114.02068525061823——5.358959032595640×10167.966445413263532×1019Figure 6
The graphs of solutionu and ut: (a) maxx,yu as t increases, and (b) maxx,yut as t increases from the initial time to the occurrence of quench behavior with q=0.8,bx=by=−2,α=β=11000 based on convection term b/αx+βy.
(a)(b)Figure 7
The graphs of (a)u as x increases, (b) ut as x increases, (c) u as y increases, and (d) ut as y increases when t is T5∗ for convection term b/1+αx+βy.
(a)(b)(c)(d)Figure 8
(a) The variation of spatial stephk as space x processes; (b) the variation of spatial step hj as space y processes; (c) the curve of the temporal steps τn as time t increases with the parameters q=0.8,bx=by=−2,α=β=11000 based on convection term b/αx+βy.
(a)(b)(c)Figure 9
(a) The 3D plot ofu when t=T6∗ and (b) the 3D plot of ut when t=T6∗ for Case 6.0.
(a)(b)The effect ofq 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 xmax is not equal to ymax when q takes 0, 0.5 while xmax is equal to ymax in the other cases; xmax,ymax is 0.0125,0.0125 when q takes 0.6, 0.7, and 0.8; xmax,ymax is 0.025,0.025 when q is in 0.9,1.4. Quenching time Tmax increases as q increases in the domain of q. Table 7 gives some of these quenching cases, in which Case 6.1.1 isq=0,bx=by=−2,α=β=11000; Case 6.1.2 is q=0.4,bx=by=−2,α=β=11000; Case 6.1.3 is q=0.5,bx=by=−2,α=β=11000; Case 6.1.4 is q=0.6,bx=by=−2,α=β=11000; Case 6.1.5 is q=0.7,bx=by=−2,α=β=11000; Case 6.1.6 is q=0.9,bx=by=−2,α=β=11000; Case 6.1.7 is q=1,bx=by=−2,α=β=11000; Case 6.1.8 is q=1.1,bx=by=−2,α=β=11000; Case 6.1.9 is q=1.2,bx=by=−2,α=β=11000; Case 6.1.10 is q=1.3,bx=by=−2,α=β=11000; Case 6.1.11 is q=1.4,bx=by=−2,α=β=11000; Case 6.0 is between Case 6.1.5 and Case 6.1.6.The paragraphs describe the relationship betweenbx(by) and quenching features. We write ten cases in Table 8. Case 6.2.1 is q=0.8,bx=by=−16,α=β=11000; Case 6.2.2 is q=0.8,bx=by=−10,α=β=11000; Case 6.2.3 is q=0.8,bx=by=−7,α=β=11000; Case 6.2.4 is q=0.8,bx=by=−5,α=β=11000; Case 6.2.5 is q=0.8,bx=by=1,α=β=11000; Case 6.2.6 is q=0.8,bx=by=12,α=β=11000; Case 6.2.7 is q=0.8,bx=by=21,α=β=11000; Case 6.2.8 is q=0.8,bx=by=40,α=β=11000; Case 6.2.9 is q=0.8,bx=by=57,α=β=11000; Case 6.2.10 is q=0.8,bx=by=80,α=β=11000; Case 6.0 is between Case 6.2.4 and Case 6.2.5. Quenching behavior can be formed when bx(by) takes the values in −16,80, in which the function curve of (b,Tmax) is concave and Tmax is smallest when b=1. When b=−16, b=−13, b∈−10,−7, b∈−5,12, b∈21,39, or b∈57,80, quenching location xmax is equal to ymax whereas when b takes other values in −16,80, xmax is different from ymax.We describe the influence ofα and β on quenching results just depending on quenching data in Table 9. Of course, we do a great quantity of experiments, from which some typical cases chosen in Table 9. We set α/β=1 and α2from 10 to 3.0×106 for Cases 6.3.1-6.3.10. Specifically, Case 6.3.1 is q=0.8,bx=by=−2,α/β=1,α2=10; Case 6.3.2 is q=0.8,bx=by=−2,α/β=1,α2=231; Case 6.3.3 is q=0.8,bx=by=−2,α/β=1,α2=1000; Case 6.3.4 is q=0.8,bx=by=−2,α/β=1,α2=2756; Case 6.3.5 is q=0.8,bx=by=−2,α/β=1,α2=8648; Case 6.3.6 is q=0.8,bx=by=−2,α/β=1,α2=10500; Case 6.3.7 is q=0.8,bx=by=−2,α/β=1,α2=10775; Case 6.3.8 is q=0.8,bx=by=−2,α/β=1,α2=10779; Case 6.3.9 is q=0.8,bx=by=−2,α/β=1,α2=19991; Case 6.3.10 is q=0.8,bx=by=−2,α/β=1,α2=300000; Case 6.0 is between Case 6.3.8 and Case 6.3.9.It does not form quenching status whenα2<10 or α2<3.0×106. After α/β is defined as 1, α2 is evaluated from 10 to 3.0×106. We find that it is more likely to produce quenching behaviors when α2 is between 10 and 3.0×106. By means of continual tests, Cases 6.3.1-6.3.10 represent 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×106. 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, Tmax increases with the increase of α2.Quenching phenomena of Case 6.0 can be considered in the following text. Its quenching time appears atT6∗=4.020449234557191, and its quenching position appears at x6∗,y6∗=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. maxx,yu is 5.358959032595640×1016 at (2.0588245639243×10−4, 9.654095827421034×10−1), and maxx,yut is 7.966445413263532×1019 at (2.080224198836234×10−5, 9.824584434865121×10−1) when quenching occurs. A comprehensive statement of quenching states with the parameters q=0.8,bx=by=−2,α=β=11000 is recorded in the next paragraphs. Although there is slight perturbation in the left boundary of ut, it does not influences on distribution of ut. It is evident that Figures 6(a) and 6(b) give two pairs curves of both maxx,yu and maxx,yut. The first one is for the distribution of the solution maxx,yut as t varies, and the second one is for the distribution of maxx,yutt as t increases. There is small perturbation near the start-up of the adaptive procedure for the left figure of maxx,yut. The four figures in Figure 7 give function relationship physical quantities betweenu, ut, and spatial variables. Figure 7(a) paints the contour of (x,u), Figure 7(b) paints the contour of (x,ut), Figure 7(c) paints the contour of (y,u), and Figure 7(d) paints the contour of (y,ut). Similarly, Figures 7(a) and 7(c) have the same shape, and Figures 7(b) and 7(d) have the same shape.There are three graphics in Figure8, which represent three function relationships between spatial step in x−direction and x, between spatial step in y−direction and y, and between temporal step and t, respectively. Figure 8(a) depicts two curves, in which one is between initial spatial steps hk and x, and another is between quenching spatial steps hk and x. Figure 8(b) also describes two curves, in which one is between initial spatial steps hk and y, and another is between quenching spatial steps hk 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τnvaries as t changes from 0 to T6∗. τ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=T6∗.The three-dimensional surfaces can render more reliable information for both the solutionu and the time derivative ut which display in Figure 9. The group of surfaces is the three-dimensional views of both u and ut when t is equal to T6∗. 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 x6∗,y6∗ is at 0.0125,0.0125 immediately before quenching occurs, the peak value of the solution approaches 1. The three-dimensional 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 maxx,yu⟶1− and maxx,yut⟶32− as t reaches T6∗.


## 5. Conclusion

Relying on the analyses in this paper, we can sort up the relationship between the parametersq, bx(by), 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 high-order difference scheme, which will probably lead to the potentially support to research the next quenching problems.



---

*Source: 3650703-2020-10-29.xml*