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

Authors: Xiaoliang Zhu; Yongbin Ge
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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


---

## Body

## 1. Introduction

Nonlinear reaction-diffusion equation with a singular or near-singular source term has been widely applied in ion conduct polarization theory [1], computational fluid dynamics [2], electromagnetism [3], material research [4], ecology [5], thermology [6, 7], and so on. Different from the linear reaction-diffusion equation, it has a chance to produce a solution with singularity. Blow-up and quenching are considered as singularity of a solution. At a limited time, the former indicates that its solution tends to infinite, and the latter indicates that its temporal derivative tends to infinite, but its solution is restricted in a certain scope [8, 9]. Researchers often use the nonlinear singular degenerate reaction-diffusion model to simulate ignition and burning and blow-up of the fuel, which can describe burning transformation from a stable status to an unstable status in a combustion chamber and from which quenching may arise [10]. To a large extent, it is to cognize features and change laws of quenching precisely that decides the design and optimization of the combustion chamber. We pay attention to the quenching problem of the degenerative singular reaction-diffusion equation with a convection term in this paper.Let us recall the development of the singular reaction-diffusion equation of the quenching type. Since Kawarada concluded that the time differential quotient of the solution leans toward infinity when the solution approximates to one indefinitely for this equation [1], more and more scholars have focused on nonlinear singular problems. The academic system of researching the quenching phenomena of the nonlinear reaction-diffusion equation has been gradually established, which ranges from one dimension and two dimension to three dimension [11–17]. More one-dimensional problems are investigated rather than higher-dimensional ones. However, there are still some quenching principles of the two-dimensional degenerate singular reaction-diffusion equations discovered [6–10, 12–15]. Specially, they qualitatively depicted quenching or nonquenching under certain conditions, asymptotic status of quenching time, the spatial definition domain size, i.e., critical length, resulting in quenching behaviors, and so on. Sheng et al. have been engaged in research in this field. Especially, Padgett and Sheng investigated the degenerate stochastic Kawarada equations from 1D to 2D and 3D by combining the finite difference schemes with arbitrary meshes [11, 16, 17]. The low-accuracy concepts are exploited in the above resolutions. Ge et al. firstly rendered a high-order compact difference scheme based on an adaptive mesh to study the quenching principle of the one-dimensional singular degenerate reaction-diffusion equation in 2018 [18].Let us recall the development of the reaction-diffusion equation with a convection term of the quenching type again. Compared to the quenching equation mentioned in the prior paragraph, less attention is paid to the convection-reaction-diffusion equation with quenching singularity, especially to the 2D and 3D problems [19]. So far, it is difficult to find research on the 2D quenching problems. The authors in [19] discussed quenching phenomena under the conditions of different convection terms and singularity indexes: asymptotic status of quenching time, the spatial definition domain size, i.e., critical length, and so on, for the one-dimension singular equation. Sheng and Cheng studied quenching behaviors of the same type equation by using the fully adaptive difference scheme of low-order accuracy in [20]. Sheng and Khliq investigated the relationship between the critical value of the degenerative term that was not a function of spatial points and both quenching time and convergence time for the 1D singular degenerate problem in [21]. Zhou et al. theoretically depicted the critical length and the quenching spatial point for the convection-reaction-diffusion equation without the degenerate term [22]. Research studies above belong to the low-order difference strategy. Recently, a high-order adaptive difference scheme is firstly used to study the influence of degenerate function, convection function, nonlinear source function, and spatial definition length on quench behaviors for the one-dimensional convection-reaction-diffusion equation, respectively [23]. It provides the basis for adopting the high-precision algorithm to analyse the quenching states of the 2D convection-reaction-diffusion equation.Depending on the previous research studies, we discover that there is an important significance for researching this two-dimensional convection-reaction-diffusion problem of the quenching type. It is well known that the finite element scheme and the finite difference scheme are popular numerical methods to solve the convection-reaction-diffusion equation [24–29]. The improved schemes based on the finite element are utilized to resolve the two-dimensional convection-reaction-diffusion problems [24, 25]. In the past years, there have emerged many finite difference schemes for solving the two-dimensional convection-reaction-diffusion problem [26–29]. Ge and Cao built up the transformation-free high-order compact difference algorithm proposed by Kalita et al. [26] on the multigrid to solve the 2D steady convection-diffusion equation [27]. Huang proposed an ADI scheme to solve the 2D unsteady convection-reaction-diffusion equation, which is built by using the first and the second central difference operators to approximate the first-order and the second-order spatial derivative, respectively, and the Crank–Nicolson (CN) method to the discrete temporal derivative [28]. Li et al. represented a compact finite difference scheme of unconditional stability and the fourth order for solving the groundwater pollutant problem based on the 2D nonlinear unsteady convection-diffusion equation [29]. Of course, there exists other schemes for solving the 2D problem. A Fourier spectral discretization method was posted for explaining the high-dimensional fractional reaction-diffusion equations by using two separate mathematical techniques [30]. Pindza and Owolabi developed the Fourier spectral method in space and the exponential integrator scheme in time to solve the 2D and 3D fractional reaction-diffusion equations [31]. Nonlinear high-dimensional convection-diffusion problems were resolved by a novel mesh-less local Petrov–Galerkin method proposed in [32]. Different from the above research studies, Owolabi and Atangana [33] studied the high-dimensional fractional reaction-diffusion systems of spatial interactions of three components’ species. The dynamical system was considered to be of asymptotical stability both in local and global style by analyzing the main equation theoretically. Furthermore, a theory certifying the existence and permanence of the three species was put forward.In general, a high-order finite difference scheme has a greater advantage comparing with other methods. An adaptive ADI difference scheme based on the two-dimensional unsteady reaction-diffusion equation with convection function is used to analyze quenching phenomena of the 2D combustion model through serial cases and investigate relationship between burner shape and quenching features. Due to significant practical meaning of quenching behaviors, the important quenching information is beneficial to engineering applications such as fuel burning and industry manufacturing. The compact scheme with the adaptive algorithm is employed to compute a series of important information including quenching time and quenching location. There are five parts in the paper. Section1 introduces the theme of this study. Section 2 describes carefully the proposed scheme of high-order accuracy. Section 3 proves adaptive mesh theory. Section 4 stimulates some typical numerical samples to explore and explain quenching problems. Section 5 draws the conclusion.


## 2. Problem and Difference Discretization

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

This semilinear degenerate problem model involving two spatial dimensions is regarded as equation (1) with the boundary and initialization conditions of equation (2) and (3). The solution ux,y,t represents the temperature of the combustion chamber. Ω∈0,α×0,β refers to the smooth rectangle domain of the combustor, 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 q≥0. The singularity source fu=1/1−uθ is strictly increasing for 0≤u0<1 with f0=f0>0,limu⟶1−fu=∞. The typical convection-reaction-diffusion equation of the 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 initialization condition is(3)ux,y,0=u0x,y,x,y∈0,α×0,β.With the aid of the intermediate variablesx^ and y^, we replace x^=x/α and y^=y/β in equation (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σαx^,βy^ut=1/α2ux^x^+1/β2uy^y^+c1αx^,βy^/αux^c2αx^,βy^/βuy^+fu,αx^,βy^∈Ω,t∈0,T. 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 and Ω∈0,1×0,1. Correspondingly, the boundary conditions of equation (5) are(6)u0,y,t=u1,y,t=ux,0,t=ux,1,t=0,t∈0,T.The initialization condition of equation (5) is(7)ux,y,0=u0x,y,x,y∈0,1×0,1.In the physical application, we rely on equations (5) and (7) to compute u and ut. Through observing a large number of values of u and ut, we can capture the quenching moment, i.e., quenching occurs when u is infinitely close to 1 and ut and becomes so huge that its value blows up.


### 2.2. Nonuniform ADI Difference Scheme

According to Taylor formula [28], the high-accuracy ADI difference scheme of equation (5) is deduced below. The steady convection-diffusion equation derived from the 1D unsteady convection-diffusion equation in the x−direction is considered as 1/α2σxuxx+c1x/ασxux=gx where g is a function of x. By virtue of the Taylors’ series expansions, we get expressions of uxk and uxxk, including the central difference operators of the first and the second spatial derivatives, discretized on the nonuniform mesh, respectively. The two expressions are substituted into the above steady convection-diffusion equation. So, a 1D steady convection-diffusion equation dispersed at point xk can be written as(8)1α2σk−c1kxR−xL2ασkδx2uk+c1kασkδxuk+ψ1xux3k+ψ2xux4k+ψ3xux5k=gk+OxR5+xL5xR+xL,where(9)ψ1x=−2xR−xL+αc1kxRxL6α2σk,ψ2x=−2xR2−xRxL+xL2+αc1kxR−xLxRxL24α2σk,ψ3x=−2xR3+xRxL2−xR2xL−xL3xR−xL−αc1k−xR3xL+xR2xL2−xRxL3120α2σk.After discretizingux3k and ux4k in the spatial direction, equation (8) is reorganized as(10)Cx−1Pxuk=gk+Cx−1ψ3xux5k+OxR5+xL5xR+xL,where(11)Px=P1xδx+P2xδx2,P1x=1ασk+P9xδx−αψ2xδx2c1k,P9x=−αψ1x+α2ψ2xc1k,P2x=1α2σk+−xR−xL2ασk+P9x−2αψ2xδxc1k,Cx=C0x+C1xδx+C2xδx2,C0x=1+αP9xδx−α2ψ2xδx2σk,C1x=αP9x−2α2ψ2xδxσk,C2x=−α2ψ2xσk.Similarly, from the one-dimensional steady convection-diffusion equation in they−direction: 1/β2σyuyy+c2y/βσyuy=gy, and the equation in the y−direction can be derived as(12)Cy−1Pyuj=gj+Cy−1ψ3yuy5j+OyR5+yL5yR+yL,where(13)Py=P1yδy+P2yδy2,P1y=1βσj+P9yδy−βψ2yδy2c2j,P9y=−βψ1y+β2ψ2yc2j,P2y=1β2σj+−yR−yL2βσj+P9y−2βψ2yδyc2j,Cy=C0y+C1yδy+C2yδy2,C0y=1+βP9yδy−β2ψ2yδy2σj,C1y=βP9y−2β2ψ2yδyσj,C2y=−β2ψ2yσj,ψ1y=−2yR−yL+βyRyLc2j6β2σj,ψ2y=−2yR2−yRyL+yL2+βyR−yLyRyLc2j24β2σj,ψ3y=−2yR3+yRyL2−yR2yL−yL3−β−yR3yL+yR2yL2−yRyL3c2j120β2σj.Adding the two hands of equation (10) and equation (12), we can get a 2D steady convection-diffusion equation:(14)Cx−1Pxuk,j+Cy−1Pyuk,j=gk,j+Cx−1ψ3xux5k,j+Cy−1ψ3yuy5k,j+OxR5+xL5xR+xL+OyR5+yL5yR+yL.We can mark the truncation errorCx−1ψ3xux5k,j+Cy−1ψ3yuy5k,j+OxR5+xL5/xR+xL+OyR5+yL5/yR+yL as OΔ and get(15)CyPx+CxPyuk,j=CxCygk,j+OΔ.Subsequently, after substitutingut−fu/σx,yk,j for gk,j in equation (14), the situation of point n+1/2 for equation (15) is considered as(16)CxCyutk,jn+1/2=CyPx+CxPyuk,jn+1/2+CxCyfuσx,yk,jn+1/2+OΔ.The Crank–Nicolson method is used to discretize the temporal term of equation (16), and an equation of point n+1 can be formed as(17)CxCy−τn2CyPx+CxPyuk,jn+1=CxCy+τn2CyPx+CxPyuk,jn+τnCxCy2σx,yfk,jn+1+fk,jn+OτnΔ+Oτn3,where τn is the time step and OτnΔ+Oτn3 is the truncation error. When τn2/4PyPxuk,jn+1 is added to the left side of equation (18) and τn2/4PyPxuk,jn+1 is added to the right and the truncation error is deleted, a group of linear system is organized as(18)Cx−τn2PxCy−τn2Pyuk,jn+1=Cx+τn2PxCy+τn2Pyuk,jn+τnCxCy2σx,yfk,jn+1+fk,jn.By virtue of the intermediate variableuk,j∗, equation (18) is reorganized as(19)Cx−τn2Pxuk,j∗=Cx+τn2PxCy+τn2Pyuk,jn+τnCxCy2σx,yfk,jn+1+fk,jn,(20)Cy−τn2Pyuk,jn+1=uk,j∗.After equation (19) is computed and arranged again, a linear system is formed as(21)B∗jU∗j=U∧jn+Fjn,n+1,j=1,2,3,…,J−1.In equation (21),(22)B∗j=b1,2x∗jb1,3x∗jb2,1x∗jb2,2x∗jb2,3x∗j⋯⋯⋯bK−2,1x∗jbK−2,2x∗jbK−2,3x∗jbK−1,1x∗jbK−1,2x∗j,U∗j=u0,j∗u1,j∗u2,j∗⋯uk,j∗⋯uK−2,j∗uK−1,j∗uK,j∗,U∧jn=u1,j∧u2,j∧⋯uk,j∧⋯uK−2,j∧uK−1,j∧,Fjn,n+1=F1,jF2,j⋯Fk,j⋯FK−2,jFK−1,j.The elements ofB∗j are(23)bk,1x∗j=2C2x−τn/2P2x−xLC1x−τn/2P1xxLxR+xL,bk,2x∗j=−2C2x−τn/2P2x+C0xxRxLxRxL,bk,3x∗j=2C2x−τn/2P2x+xRC1x−τn/2P1xxRxR+xL.The elementuk,j∧ of U∧jn is(24)uk,j∧=dk,1x∗jdk,1y∗juk−1,j−1n+dk,2y∗juk−1,jn+dk,3y∗juk−1,j+1n+dk,2x∗jdk,1y∗juk,j−1n+dk,2y∗juk,jn+dk,3y∗juk,j+1n+dk,3x∗jdk,1y∗juk+1,j−1n+dk,2y∗juk+1,jn+dk,3y∗juk+1,j+1n,where(25)dk,1x∗j=2C2x+τn/2P2x−xLC1x+τn/2P1xxLxR+xL,dk,2x∗j=−2C2x+τn/2P2x+C0xxRxLxRxL,dk,3x∗j=2C2x+τn/2P2x+xRC1x+τn/2P1xxRxR+xL,dk,1y∗j=2C2y+τn/2P2y−yLC1y+τn/2P1yyLyR+yL,dk,2y∗j=−2C2y+τn/2P2y+C0yyRyLyRyL,dk,3y∗j=2C2y+τn/2P2y+yRC1y+τn/2P1yyRyR+yL.The elementFk,j of Fjn,n+1 is(26)Fk,j=τn2σsk,1x∗jsk,1y∗jfk−1,j−1n+fk−1,j−1n+1+sk,2y∗jfk−1,jn+fk−1,jn+1+sk,3y∗jfk−1,j+1n+fk−1,j+1n+1+τn2σsk,2x∗jsk,1y∗jfk,j−1n+fk,j−1n+1+sk,2y∗jfk,jn+fk,jn+1+sk,3y∗jfk,j+1n+fk,j+1n+1+τn2σsk,3x∗jsk,1y∗jfk+1,j−1n+fk+1,j−1n+1+sk+1,2y∗jfk+1,jn+fk+1,jn+1+sk,3y∗jfk+1,j+1n+fk+1,j+1n+1,where(27)sk,1x∗j=2C2x−xLC1xxLxR+xL,sk,2x∗j=−2C2x+C0xxRxLxRxL,sk,3x∗j=2C2x+xRC1xxRxR+xL,sk,1y∗j=2C2y−yLC1yyLyR+yL,sk,2y∗j=−2C2y+C0yyRyLyRyL,sk,3y∗j=2C2y+yRC1yyRyR+yL.Similarly, from equation (20), another linear system is written as(28)BUn+1=U∗,where(29)B=b1,2yb1,3yb2,1yb2,2yb2,3y⋯⋯⋯bJ−2,1ybJ−2,2ybJ−2,3ybJ−1,1ybJ−1,2y,in which(30)bj,1y=2C2y−τn/2P2y−yLC1y−τn/2P1yyLyR+yL,bj,2y=−2C2y−τn/2P2y+C0yyRyLyRyL,bj,3y=2C2y−τn/2P2y+yRC1y−τn/2P1yyRyR+yL.It can be observed that eachuk,jn+1 on the n+1th time layer is obtained, according to uk,jn on the prior one, by computing the tridiagonal linear equations twice. We adopt the catch-up method to solve the linear equations. We rely on equation (21) to get the medial variable uk,j∗, which represents the x−direction. Subsequently, we depend on uk,j∗ and equation (28) to gain uk,jn+1, which carries information in the y−direction and so on, until the solutions on the last time axis are worked out. The accuracy analysis of equation (18) is performed, and its truncation error is(31)Oτn2+xR2+xL2xR−xL+yR2+yL2yR−yL+xR5+xL5xR+xL+yR5+yL5yR+yL.Note that equations (19) and (20) have the second-order accuracy in time and the third- or fourth-order accuracy in space. The scheme of uniform owns spatial accuracy of the fourth order, whereas its nonuniform scheme owns spatial accuracy of the third order.


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

This semilinear degenerate problem model involving two spatial dimensions is regarded as equation (1) with the boundary and initialization conditions of equation (2) and (3). The solution ux,y,t represents the temperature of the combustion chamber. Ω∈0,α×0,β refers to the smooth rectangle domain of the combustor, 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 q≥0. The singularity source fu=1/1−uθ is strictly increasing for 0≤u0<1 with f0=f0>0,limu⟶1−fu=∞. The typical convection-reaction-diffusion equation of the 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 initialization condition is(3)ux,y,0=u0x,y,x,y∈0,α×0,β.With the aid of the intermediate variablesx^ and y^, we replace x^=x/α and y^=y/β in equation (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σαx^,βy^ut=1/α2ux^x^+1/β2uy^y^+c1αx^,βy^/αux^c2αx^,βy^/βuy^+fu,αx^,βy^∈Ω,t∈0,T. 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 and Ω∈0,1×0,1. Correspondingly, the boundary conditions of equation (5) are(6)u0,y,t=u1,y,t=ux,0,t=ux,1,t=0,t∈0,T.The initialization condition of equation (5) is(7)ux,y,0=u0x,y,x,y∈0,1×0,1.In the physical application, we rely on equations (5) and (7) to compute u and ut. Through observing a large number of values of u and ut, we can capture the quenching moment, i.e., quenching occurs when u is infinitely close to 1 and ut and becomes so huge that its value blows up.


## 2.2. Nonuniform ADI Difference Scheme

According to Taylor formula [28], the high-accuracy ADI difference scheme of equation (5) is deduced below. The steady convection-diffusion equation derived from the 1D unsteady convection-diffusion equation in the x−direction is considered as 1/α2σxuxx+c1x/ασxux=gx where g is a function of x. By virtue of the Taylors’ series expansions, we get expressions of uxk and uxxk, including the central difference operators of the first and the second spatial derivatives, discretized on the nonuniform mesh, respectively. The two expressions are substituted into the above steady convection-diffusion equation. So, a 1D steady convection-diffusion equation dispersed at point xk can be written as(8)1α2σk−c1kxR−xL2ασkδx2uk+c1kασkδxuk+ψ1xux3k+ψ2xux4k+ψ3xux5k=gk+OxR5+xL5xR+xL,where(9)ψ1x=−2xR−xL+αc1kxRxL6α2σk,ψ2x=−2xR2−xRxL+xL2+αc1kxR−xLxRxL24α2σk,ψ3x=−2xR3+xRxL2−xR2xL−xL3xR−xL−αc1k−xR3xL+xR2xL2−xRxL3120α2σk.After discretizingux3k and ux4k in the spatial direction, equation (8) is reorganized as(10)Cx−1Pxuk=gk+Cx−1ψ3xux5k+OxR5+xL5xR+xL,where(11)Px=P1xδx+P2xδx2,P1x=1ασk+P9xδx−αψ2xδx2c1k,P9x=−αψ1x+α2ψ2xc1k,P2x=1α2σk+−xR−xL2ασk+P9x−2αψ2xδxc1k,Cx=C0x+C1xδx+C2xδx2,C0x=1+αP9xδx−α2ψ2xδx2σk,C1x=αP9x−2α2ψ2xδxσk,C2x=−α2ψ2xσk.Similarly, from the one-dimensional steady convection-diffusion equation in they−direction: 1/β2σyuyy+c2y/βσyuy=gy, and the equation in the y−direction can be derived as(12)Cy−1Pyuj=gj+Cy−1ψ3yuy5j+OyR5+yL5yR+yL,where(13)Py=P1yδy+P2yδy2,P1y=1βσj+P9yδy−βψ2yδy2c2j,P9y=−βψ1y+β2ψ2yc2j,P2y=1β2σj+−yR−yL2βσj+P9y−2βψ2yδyc2j,Cy=C0y+C1yδy+C2yδy2,C0y=1+βP9yδy−β2ψ2yδy2σj,C1y=βP9y−2β2ψ2yδyσj,C2y=−β2ψ2yσj,ψ1y=−2yR−yL+βyRyLc2j6β2σj,ψ2y=−2yR2−yRyL+yL2+βyR−yLyRyLc2j24β2σj,ψ3y=−2yR3+yRyL2−yR2yL−yL3−β−yR3yL+yR2yL2−yRyL3c2j120β2σj.Adding the two hands of equation (10) and equation (12), we can get a 2D steady convection-diffusion equation:(14)Cx−1Pxuk,j+Cy−1Pyuk,j=gk,j+Cx−1ψ3xux5k,j+Cy−1ψ3yuy5k,j+OxR5+xL5xR+xL+OyR5+yL5yR+yL.We can mark the truncation errorCx−1ψ3xux5k,j+Cy−1ψ3yuy5k,j+OxR5+xL5/xR+xL+OyR5+yL5/yR+yL as OΔ and get(15)CyPx+CxPyuk,j=CxCygk,j+OΔ.Subsequently, after substitutingut−fu/σx,yk,j for gk,j in equation (14), the situation of point n+1/2 for equation (15) is considered as(16)CxCyutk,jn+1/2=CyPx+CxPyuk,jn+1/2+CxCyfuσx,yk,jn+1/2+OΔ.The Crank–Nicolson method is used to discretize the temporal term of equation (16), and an equation of point n+1 can be formed as(17)CxCy−τn2CyPx+CxPyuk,jn+1=CxCy+τn2CyPx+CxPyuk,jn+τnCxCy2σx,yfk,jn+1+fk,jn+OτnΔ+Oτn3,where τn is the time step and OτnΔ+Oτn3 is the truncation error. When τn2/4PyPxuk,jn+1 is added to the left side of equation (18) and τn2/4PyPxuk,jn+1 is added to the right and the truncation error is deleted, a group of linear system is organized as(18)Cx−τn2PxCy−τn2Pyuk,jn+1=Cx+τn2PxCy+τn2Pyuk,jn+τnCxCy2σx,yfk,jn+1+fk,jn.By virtue of the intermediate variableuk,j∗, equation (18) is reorganized as(19)Cx−τn2Pxuk,j∗=Cx+τn2PxCy+τn2Pyuk,jn+τnCxCy2σx,yfk,jn+1+fk,jn,(20)Cy−τn2Pyuk,jn+1=uk,j∗.After equation (19) is computed and arranged again, a linear system is formed as(21)B∗jU∗j=U∧jn+Fjn,n+1,j=1,2,3,…,J−1.In equation (21),(22)B∗j=b1,2x∗jb1,3x∗jb2,1x∗jb2,2x∗jb2,3x∗j⋯⋯⋯bK−2,1x∗jbK−2,2x∗jbK−2,3x∗jbK−1,1x∗jbK−1,2x∗j,U∗j=u0,j∗u1,j∗u2,j∗⋯uk,j∗⋯uK−2,j∗uK−1,j∗uK,j∗,U∧jn=u1,j∧u2,j∧⋯uk,j∧⋯uK−2,j∧uK−1,j∧,Fjn,n+1=F1,jF2,j⋯Fk,j⋯FK−2,jFK−1,j.The elements ofB∗j are(23)bk,1x∗j=2C2x−τn/2P2x−xLC1x−τn/2P1xxLxR+xL,bk,2x∗j=−2C2x−τn/2P2x+C0xxRxLxRxL,bk,3x∗j=2C2x−τn/2P2x+xRC1x−τn/2P1xxRxR+xL.The elementuk,j∧ of U∧jn is(24)uk,j∧=dk,1x∗jdk,1y∗juk−1,j−1n+dk,2y∗juk−1,jn+dk,3y∗juk−1,j+1n+dk,2x∗jdk,1y∗juk,j−1n+dk,2y∗juk,jn+dk,3y∗juk,j+1n+dk,3x∗jdk,1y∗juk+1,j−1n+dk,2y∗juk+1,jn+dk,3y∗juk+1,j+1n,where(25)dk,1x∗j=2C2x+τn/2P2x−xLC1x+τn/2P1xxLxR+xL,dk,2x∗j=−2C2x+τn/2P2x+C0xxRxLxRxL,dk,3x∗j=2C2x+τn/2P2x+xRC1x+τn/2P1xxRxR+xL,dk,1y∗j=2C2y+τn/2P2y−yLC1y+τn/2P1yyLyR+yL,dk,2y∗j=−2C2y+τn/2P2y+C0yyRyLyRyL,dk,3y∗j=2C2y+τn/2P2y+yRC1y+τn/2P1yyRyR+yL.The elementFk,j of Fjn,n+1 is(26)Fk,j=τn2σsk,1x∗jsk,1y∗jfk−1,j−1n+fk−1,j−1n+1+sk,2y∗jfk−1,jn+fk−1,jn+1+sk,3y∗jfk−1,j+1n+fk−1,j+1n+1+τn2σsk,2x∗jsk,1y∗jfk,j−1n+fk,j−1n+1+sk,2y∗jfk,jn+fk,jn+1+sk,3y∗jfk,j+1n+fk,j+1n+1+τn2σsk,3x∗jsk,1y∗jfk+1,j−1n+fk+1,j−1n+1+sk+1,2y∗jfk+1,jn+fk+1,jn+1+sk,3y∗jfk+1,j+1n+fk+1,j+1n+1,where(27)sk,1x∗j=2C2x−xLC1xxLxR+xL,sk,2x∗j=−2C2x+C0xxRxLxRxL,sk,3x∗j=2C2x+xRC1xxRxR+xL,sk,1y∗j=2C2y−yLC1yyLyR+yL,sk,2y∗j=−2C2y+C0yyRyLyRyL,sk,3y∗j=2C2y+yRC1yyRyR+yL.Similarly, from equation (20), another linear system is written as(28)BUn+1=U∗,where(29)B=b1,2yb1,3yb2,1yb2,2yb2,3y⋯⋯⋯bJ−2,1ybJ−2,2ybJ−2,3ybJ−1,1ybJ−1,2y,in which(30)bj,1y=2C2y−τn/2P2y−yLC1y−τn/2P1yyLyR+yL,bj,2y=−2C2y−τn/2P2y+C0yyRyLyRyL,bj,3y=2C2y−τn/2P2y+yRC1y−τn/2P1yyRyR+yL.It can be observed that eachuk,jn+1 on the n+1th time layer is obtained, according to uk,jn on the prior one, by computing the tridiagonal linear equations twice. We adopt the catch-up method to solve the linear equations. We rely on equation (21) to get the medial variable uk,j∗, which represents the x−direction. Subsequently, we depend on uk,j∗ and equation (28) to gain uk,jn+1, which carries information in the y−direction and so on, until the solutions on the last time axis are worked out. The accuracy analysis of equation (18) is performed, and its truncation error is(31)Oτn2+xR2+xL2xR−xL+yR2+yL2yR−yL+xR5+xL5xR+xL+yR5+yL5yR+yL.Note that equations (19) and (20) have the second-order accuracy in time and the third- or fourth-order accuracy in space. The scheme of uniform owns spatial accuracy of the fourth order, whereas its nonuniform scheme owns spatial accuracy of the third order.


## 3. Adaptive Mesh Method

### 3.1. Temporal Adaptive Mesh Method

Comparing with solutionu, the temporal derivative ut is more sensible. Therefore, like the 1D adaptive method, the 2D adaptive mesh method in the time direction is based on the arc-length monitor function of the temporal derivative function ut. We can extend the 1D adaptive temporal grid schemes in Reference [23] to the 2D case:(32)τnτn−1=max0<k<K0<j<J1+uttxk,yj,tn−1/22max0<k<K0<j<J1+uttxk,yj,tn+1/22,n=1,2,….Afterutt in equation (32) is discretized through the central difference formulas to calculate the maximum ratio, the adaptive time step τn is easily obtained by equation (32) and τn−1 in the prior time layer.


### 3.2. Spatial Adaptive Mesh Method

Similarly, we can utilize the arc-length monitor function in the space direction based onut and extend from the 1D adaptive spatial algorithm to the 2D problem by relying on Reference [23]. Therefore, Wxx,t representing the monitor function in the x−direction and Wyy,t representing the monitor function in the y−direction can be written as(33)Wxx,t=1+∑k=0Kutxx,t2,0≤x≤1,0≤t<T,(34)Wyy,t=1+∑j=0Jutyy,t2,0≤y≤1,0≤t<T.The two spatial monitor functions are positive in their definition area.ξkxk,tk=0K and ζjyj,tj=0J represent the set of 2D original spatial mesh points. After an adaptive improvement is complemented, the grid in the 2D computational area is updated as xkξk,tk=0K in the x−direction and yjζj,tj=0J in the y−direction. It is easily observed that the 1D adaptive grid scheme acts on in the x−direction and y−direction, respectively, to produce new grid points.We take the equidistributing mesh in they−direction as an example below. From [23], we know that the subareas of Wyy,t in intervals of yj−1t,yjtj=1,2,…,J are equal. Obviously, when Wy≡1, the scope 0,1 is uniform. It can be expressed as(35)∫y0y1Wyy,tdy=∫y1y2Wyy,tdy=⋯=∫yJ−1yJWyy,tdy=1J∫01Wyy,tdy=ςyt,where ςyt is the average arc-length. Equation (35) is equal to(36)∫yj−1yjWyy,tdy=∫ζj−1ζjWyyζ,t,ty′ζ,tdζ,j=1,2,…,J.After the second derivation of the aforementioned formula, the following equation can be obtained:(37)∂∂ζWyyjm,t∂∂ζyjm+1=0,where m represents the iterative number, yjmj=0J represents the spatial grid node set before iteration, and yjm+1j=0J represents the space grid node set after iteration. Equation (37) can be solved by the difference scheme:(38)Wyyj+1/2m,τnyj+1m+1−yjm+1−Wyyj−1/2m,τnyjm+1−yj−1m+1=0,j=1,2,…,J−1,n≥0,which approximates to the former. Similarly, the semi-implicit scheme in the x−direction is defined as(39)Wxxk+1/2m,τnxk+1m+1−xkm+1−Wxxk−1/2m,τnxkm+1−xk−1m+1=0,k=1,2,…,K−1,n≥0.


### 3.3. Iterative Adaption Method

The adaptive iterative procedure of the presented method is carried forward as follows. Firstly, the spatial monitor functionsWxx,t and Wyy,t can be obtained through equations (33) and (34) based on the grids xkmk=0K and yjmj=0J in the nth time layer. Secondly, depending on equations (38) and (39), we can gain xkm+1k=0K derived from xkmk=0K in the x−direction and yjm+1j=0J derived from yjmj=0J in the y−direction. The new points substitute for the old points iteratively when ∑k=1K−1xkm+1−xkm<er and ∑j=1J−1yjm+1−yjm<er are reached for 0<er≤1. mmax is the terminal value of m+1. Thirdly, from the value built on the last grid point sets xkmmaxk=0K and yjmmaxj=0J, uk,jnk=1,j=1K−1,J−1 is computed by using the method of area ratio. Thus, we can get uk,jn+1k=1,j=1K−1,J−1 in the n+1th time layer depending on equation (18). Lastly, repeat the aforementioned produce till quenching occurs or the solution converges to a steady solution.


## 3.1. Temporal Adaptive Mesh Method

Comparing with solutionu, the temporal derivative ut is more sensible. Therefore, like the 1D adaptive method, the 2D adaptive mesh method in the time direction is based on the arc-length monitor function of the temporal derivative function ut. We can extend the 1D adaptive temporal grid schemes in Reference [23] to the 2D case:(32)τnτn−1=max0<k<K0<j<J1+uttxk,yj,tn−1/22max0<k<K0<j<J1+uttxk,yj,tn+1/22,n=1,2,….Afterutt in equation (32) is discretized through the central difference formulas to calculate the maximum ratio, the adaptive time step τn is easily obtained by equation (32) and τn−1 in the prior time layer.


## 3.2. Spatial Adaptive Mesh Method

Similarly, we can utilize the arc-length monitor function in the space direction based onut and extend from the 1D adaptive spatial algorithm to the 2D problem by relying on Reference [23]. Therefore, Wxx,t representing the monitor function in the x−direction and Wyy,t representing the monitor function in the y−direction can be written as(33)Wxx,t=1+∑k=0Kutxx,t2,0≤x≤1,0≤t<T,(34)Wyy,t=1+∑j=0Jutyy,t2,0≤y≤1,0≤t<T.The two spatial monitor functions are positive in their definition area.ξkxk,tk=0K and ζjyj,tj=0J represent the set of 2D original spatial mesh points. After an adaptive improvement is complemented, the grid in the 2D computational area is updated as xkξk,tk=0K in the x−direction and yjζj,tj=0J in the y−direction. It is easily observed that the 1D adaptive grid scheme acts on in the x−direction and y−direction, respectively, to produce new grid points.We take the equidistributing mesh in they−direction as an example below. From [23], we know that the subareas of Wyy,t in intervals of yj−1t,yjtj=1,2,…,J are equal. Obviously, when Wy≡1, the scope 0,1 is uniform. It can be expressed as(35)∫y0y1Wyy,tdy=∫y1y2Wyy,tdy=⋯=∫yJ−1yJWyy,tdy=1J∫01Wyy,tdy=ςyt,where ςyt is the average arc-length. Equation (35) is equal to(36)∫yj−1yjWyy,tdy=∫ζj−1ζjWyyζ,t,ty′ζ,tdζ,j=1,2,…,J.After the second derivation of the aforementioned formula, the following equation can be obtained:(37)∂∂ζWyyjm,t∂∂ζyjm+1=0,where m represents the iterative number, yjmj=0J represents the spatial grid node set before iteration, and yjm+1j=0J represents the space grid node set after iteration. Equation (37) can be solved by the difference scheme:(38)Wyyj+1/2m,τnyj+1m+1−yjm+1−Wyyj−1/2m,τnyjm+1−yj−1m+1=0,j=1,2,…,J−1,n≥0,which approximates to the former. Similarly, the semi-implicit scheme in the x−direction is defined as(39)Wxxk+1/2m,τnxk+1m+1−xkm+1−Wxxk−1/2m,τnxkm+1−xk−1m+1=0,k=1,2,…,K−1,n≥0.


## 3.3. Iterative Adaption Method

The adaptive iterative procedure of the presented method is carried forward as follows. Firstly, the spatial monitor functionsWxx,t and Wyy,t can be obtained through equations (33) and (34) based on the grids xkmk=0K and yjmj=0J in the nth time layer. Secondly, depending on equations (38) and (39), we can gain xkm+1k=0K derived from xkmk=0K in the x−direction and yjm+1j=0J derived from yjmj=0J in the y−direction. The new points substitute for the old points iteratively when ∑k=1K−1xkm+1−xkm<er and ∑j=1J−1yjm+1−yjm<er are reached for 0<er≤1. mmax is the terminal value of m+1. Thirdly, from the value built on the last grid point sets xkmmaxk=0K and yjmmaxj=0J, uk,jnk=1,j=1K−1,J−1 is computed by using the method of area ratio. Thus, we can get uk,jn+1k=1,j=1K−1,J−1 in the n+1th time layer depending on equation (18). Lastly, repeat the aforementioned produce till quenching occurs or the solution converges to a steady solution.


## 4. Numerical Illustration

### 4.1. Common Example

It is not a quenching case. It is just for demonstrating the accuracy and convergence of the proposed difference scheme. We have solved the following unsteady convection-diffusion equation with coefficient whose exact solution is offered:(40)put=uxx+uyy−pux−puy,x,y∈0,1×0,1,t≥0,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.The initial and boundary value conditions are derived from the exact solution. There are four metrics for the data test, which are maximum error, average error, CPU time, and convergence rate. We adopt these parametersT=0.5,τn=0.01,andp=100 to fulfill the two difference schemes: the proposed scheme based on the uniform mesh and the proposed scheme based on the nonuniform mesh. k/K+λx/πsinπk/Kk=0K and j/J+λy/πsinπj/Jj=0J represent the nonuniform grid, where KandJ point to intervals of the spatial direction and K=J, whereas λxandλy point to the mesh transformation coefficient and λx=λy.Table1 offers these measuring data. It displays the accuracies and running time of our scheme (uniform and nonuniform). We select five spatial intervals: 16, 32, 64, 128, and 256 as test cases, which are given to each of the above two meshes (i.e., uniform and nonuniform). For the nonuniform mesh, the different values of λx impact these precision criteria. In the paper, λx takes 0.7. The maximum errors of the first scheme are higher than those of the second one, and the average errors of the latter are superior to those of the former. Furthermore, the larger the number of spatial intervals is, the higher the accuracies of the scheme based on the nonuniform grid are. Its accuracies arrive at 1.3364 × 10−7 in the maximum error and 6.4492 × 10−9 in the average error when K=J=256. In terms of computing time of the two schemes whose measurement unit is millisecond, the first scheme is lower than the second one. Obviously, by comparing with the uniform mesh, the difference scheme based on the nonuniform mesh spends more time to run. Lastly, the conv. rate values are observed. The tendency of the conv. rate from the former scheme is similar to that from the latter scheme. We take the latter as the example. The four values are positive, which means the accuracy increases as K rises when K∈16,256. The conv. rate value is maximum when K transfers from 16 to 32, which means the accuracy promotes greatly in the stage. The latter three values are close, which means the accuracy promotes gently in the three stages.Table 1
The parameter values for common example.
Difference schemesKJMax. errorAver. errorCPU timeConv. rateThe proposed scheme (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.1743 × 10−7117943.9714The proposed scheme (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.01432561.3364 × 10−76.4492 × 10−91659693.9867


### 4.2. Quenching Case without Convection Term

This section will describe a quenching status without the convection term, i.e.,c1x,y=0 and c2x,y=0. The mathematical model can be expressed as(41)α2x2+β2y2q/2ut=1α2uxx+1β2uyy+11−uθ,x,y∈0,1×0,1,t∈0,T.Its initial and boundary conditions are just as in equations (6) and (7).After equations (19) and (20) are used to approximate equation (41), we can get some quenching behaviors. For the 2D nonlinear degeneration singularity equation without the convection term, we demonstrate the quenching case of q=0 that has been described in [13, 14]. In Table 2 and Figure 1, q=0, a=β=10, and θ=1 which are the same for the three approaches, but the initial space step h0 and the initial temporal step τ0 are different among them, where we take h0=1/50 and τ0=4.166677×10−3. In Table 2, the five parameters represent the same senses, in which xmax,ymax represents the point of quenching location, tmax represents the quenching time, maxurepresents the maximum temperature value immediately before quenching happens, and maxut means the maximum variation of temperature with respect to time immediately before quenching happens. From Table 2, we can see that the quenching location 0.5,0.5 of the referred scheme is identical with that of the scheme in [14], and tmax=0.5958348109999995 and maxu=0.998562875689732 of the first one are close to those of the last two. Obviously, maxut is different among the various methods. From Figure 1, it is found that u⟶1− and ut⟶690−at the location xmax,ymax=0.5,0.5 and the quenching time tmax=0.5958348109999995.Table 2
Quenching performances of different schemes forq=0, a=β=10, and θ=1 without the convection term.
SchemestmaxxmaxymaxmaxumaxutThe proposed scheme0.59583481099999950.5000.5000.998562875689732689.2724781649539The scheme in [13]0.58712499993751——0.990432148.887767The scheme in [14]0.5875540.5000.5000.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 the convection term.
(a)(b)


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

The standard 2D degenerate semilinear quenching problem with the initial boundary condition lists in equations (5)–(7) is extended from the equation in one-space dimension [23]. The first convection term b/1+αx+βy can be written as two functions c1x,y=bx/1+αx+βy and c2x,y=by/1+αx+βy, where bx and by take the constants. The four groups of variables: q, bx and by, α and β, and θ in equations (5)–(7) affect the quenching characteristics of the 2D problem, whose meaning has been stated in Section 2. There follows a special case for b/1+αx+βy. The initial time step is considered as τ0=0.0003, and the initial space step is considered as h0=1/55, in which the x-direction is identical with the y-direction. The four groups of parameters are configured as q=0.5,α=β=200,bx=by=10,andθ=1 which is defined as Case 7.0. This will lead to a quenching phenomenon. Quenching location x7∗,y7∗ is 0.03636363636363636,0.03636363636363636; the quenching time T7∗ is 4.487776579638927; maxx,yu is 0.9992409498317957; maxx,yut is 1539.153786381783. When t is 4.486155872847072, the temporal step starts to become smaller and smaller until quenching occurs. We choose six rows of representative data given in Table 3 which poins to features relative to x, y, t, u, and ut before quenching occurs of Case 7.0. maxx,yu and maxx,yut refer to maximal u and maximal ut in each time layer before quenching occurs; x and y refer to the coordinate point of synchronously producing maxx,yu and maxx,yut in each time layer before quenching occurs; t refers to the moment at which each time layer is located before quenching occurs. x,y first appears at the point x7∗,y7∗ when t is 4.487088240486234 which is the first row in Table 3 and maintains at this point until quenching occurs. When t is T7∗ which is the immediate time before quenching happens, we can obtain important characteristics. When t is 4.487777138788559, maxx,yut blows up and reaches 1.051076326841211 × 109.Table 3
Time, location,maxx,yu, and maxx,yut in the six temporal layers which contain quenching occurrence for Case 7.0.
txymaxx,yumaxx,yut4.4870882404862340.036363636363636360.036363636363636360.979448569828961552.184431339763224.4874898808006030.036363636363636360.036363636363636360.985461462054538975.528997590273484.4876538240986460.036363636363636360.036363636363636360.9897483333199262109.07447404344314.4877707022550420.036363636363636360.036363636363636360.9972318816440022418.00237744504054.4877765796389270.036363636363636360.036363636363636360.99924094983179571539.1537863817834.487777138788559———1.051076326841211 × 109Figures2–7 provide seven groups of graphs to illustrate more particular quenching statuses. Figure 2 shows the 2D plots of maxx,yu and maxx,yut changing as t rises from the initial moment until quenching occurs. Figure 2(a) corresponds to u, and Figure 2(b) corresponds to ut. Figure 3 sketches the 2D contours of u and ut regarding x immediately before quenching happens. Figure 3(a) corresponds to u, and Figure 3(b) corresponds to ut. Figure 4 draws the 2D outlines of u and ut regarding y immediately before quenching happens. Figure 4(a) corresponds to u, and Figure 4(b) corresponds to ut. The three 2D plots in Figure 5 describe the adaptive mechanism before quenching occurs. The first one reflects spatial adaptation, and the last two reflect temporal adaptation. Figure 5(a) describes the trend of the temporal step over time. To distinguish more clearly the variation, we intercept the curves of the period from 4.470000000001407 to 4.487777138788559. The red sign on the blue curve in Figure 5(a) denotes the initial temporal adaptive moment (i.e., 4.486155872847072). The temporal step is uniform before the adaptive moment but becomes nonuniform when t varies from 4.486155872847072 to 4.487777138788559. The green sign on the blue curve in Figure 5(a) denotes the initial moment (i.e., 4.487088240486234) of producing the quenching location x7∗,y7∗. The green lines in Figures 5(b) and 5(c) describe the uniform spatial steps in the x−direction and the y−direction, respectively. The spatial adaptation occurs at 4.487777138788559 after which the spatial step becomes nonuniform. The red curves with blue spots in Figures 5(b) and 5(c) describe, respectively, the nonuniform spatial steps in the x−direction and the y−direction at T7∗. Although the two curves are very similar, there is still subtle difference between their data. Figures 6–7 provide the 3D surfaces of u and ut with respect to spatial variables x and y at certain special time. Figure 6(a) records the 3D surfaces about x, y, and u, and Figure 6(b) records the 3D surfaces about x, y, and ut when t=4.487770702255042. Figure 7(a) records the 3D surfaces about x, y, and u, and Figure 7(b) records the 3D surfaces about x, y, and ut when T7∗.Figure 2
The graphs of (a)maxx,yu as t increases and (b) maxx,yut as t increases from the initial moment to the occurrence of quenching behavior for the convection term b/1+αx+βy.
(a)(b)Figure 3
The graphs of (a)u as x increases and (b) ut as x increases when t is T7∗ for the convection term b/1+αx+βy.
(a)(b)Figure 4
The graphs of (a)u as y increases and (b) ut as y increases when t is T7∗ for the convection term b/1+αx+βy.
(a)(b)Figure 5
(a) The curve of the temporal stepsτn as time t increases, (b) the variation of spatial steps hk as x processes, and (c) the variation of spatial steps hk as y processes and for the convection term b/1+αx+βy.
(a)(b)(c)Figure 6
(a) The 3D plot ofu and (b) the 3D plot of ut when t=4.487770702255042 for the convection term b/1+αx+βy.
(a)(b)Figure 7
(a) The 3D plot ofu and (b) the 3D plot of ut when t=T7∗ for the convection term b/1+αx+βy.
(a)(b)


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

The second convection termb/αx+βy can be written as the two convection functions: c1x,y=bx/αx+βy and c2x,y=by/αx+βy, where bx and by take the constants. Comparing with the first convection term b/1+αx+βy, it is hard to conduct quenching behavior for the second convection term b/αx+βy. We make a following configuration to gain a quenching status. The initial x−step (y−step) is h0=1/40, and the initial t−step is τ0=0.0005. The four groups of important parameters are set as q=1.1,α=β=20,bx=by=3,andθ=1 which is defined as Case 8.0. The parameter statuses of the case are recorded in Table 4. For the case, the initial time of producing the quenching location is earlier than the temporal adaptive time. The former time is 5.596499999999618, which is the first row in Table 4, the latter one is 5.603982012140827, which is the fifth row in Table 4, and the quenching location x8∗,y8∗ is 0.025,0.025. Moreover, it is so transient from the front time of quenching location appearing to the terminative time of quenching location appearing that there are 18 time layers between the two moments. The latter quenching time is T8∗=5.604262757909247, maxx,yu is 0.9835341519876587, and maxx,yut is 73.15674241637485 immediately before quenching occurs, which is the sixth row in Table 4. maxx,yut gets the value 6.637517818582074 × 1011 and blows up when t is 5.604272047535358, which is the last row in Table 4.Table 4
Time, location,maxx,yu, and maxx,yut in the last seven temporal layers which contains quenching occurrence for Case 8.0.
txymaxx,yumaxx,yut5.5964999999996180.0250.0250.965613248035384527.204925805149065.5984999999996170.0250.0250.966848091050629728.766902187783995.5999999999996160.0250.0250.967807952885317130.064978162417135.6019999999996150.0250.0250.96913921620335732.000736898491995.6039820121408270.0250.0250.970836300971370934.727044766461555.6042627579092470.0250.0250.983534151987658773.156742416374855.604272047535358———6.637517818582074 × 1011This paragraph supplies five sets of figures (i.e., Figures8–12) to more vividly depict the quenching information of the case. The 2D curves in Figure 8 describe the relationship between maxx,yu and t and the relationship between maxx,yut and t. The 2D curves in Figure 9 describe the relationship between u and x and the relationship between ut and x. The 2D curves in Figure 10 describe the relationship between u and y and the relationship between ut and y. In Figures 8–10, (a) corresponds to u and (b) corresponds to ut. Figure 11(a) draws the tendency of the temporal step varying with temporal variable t; Figure 11(b) draws the tendency of the initial spatial step and the adaptive step varying with spatial variable x; and Figure 11(c) draws the tendency of the initial spatial step and the adaptive step varying with spatial variable y. For Figure 11(a), to distinguish more clearly the variation, we intercept the curves of the period from 5.580499999999627 to 5.604272047535358; and the red spot on the blue curve expresses the earliest temporal adaptive time (i.e., 5.603982012140827). maxx,yut in the time step becomes 1.185615675885128 × 1012 when t is 5.603982012140827. The temporal step is uniform before the adaptive moment but becomes nonuniform since then. The green spot on the blue curve expresses the initial moment of producing the quenching location x8∗,y8∗, which is 5.604262757909247. In Figures 11(b) and 11(c), the green lines indicate the uniform initial spatial steps while the red curves with blue signs indicate the nonuniform adaptive spatial steps. Lastly, we provide two 3D surfaces in Figure 12 to state quenching statuses of u and ut when t=T8∗. Figure 12(a) denotes the relationship between u, x, and y, and Figure 12(b) denotes the relationship between ut, x, and y. Figure 12(a) records the 3D surfaces about x, y, and u, and Figure 12(b) records the 3D surfaces about x, y, and ut when t=T8∗. There is a narrowband projecting the rear plane of the 3D axis box in Figure 12(b). This is because some larger values of maxx,yut appear at the positions, in which the x−axis and the y−axis are simultaneously near the origin of the coordinates. From the data source of the 3D surface, we can observe some larger maxx,yut varying from 1 to 27, whereas maxx,yut at other positions are positive numbers less than 1, which can be seen via Figures 9(b) and 10(b). Obviously, the larger maxx,yut are all less than 73.15674241637485, which is reasonable for a quenching case.Figure 8
The graphs of (a)maxx,yu as t increases and (b) maxx,yut as t increases from the initial moment to the occurrence of quenching behavior for the convection term b/αx+βy.
(a)(b)Figure 9
The graphs of (a)u as x increases and (b) ut as x increases when t is T8∗ for the convection term b/αx+βy.
(a)(b)Figure 10
The graphs of (a)u as y increases and (b) ut as y increases when t is T8∗ for the convection term b/αx+βy.
(a)(b)Figure 11
(a) The curve of the temporal stepsτn as time t increases, (b) the variation of spatial steps hk as x processes, and (c) the variation of spatial steps hk as y processes and for the convection term b/αx+βy.
(a)(b)(c)Figure 12
(a) The 3D plot ofu and (b) the 3D plot of ut when t=T8∗ for the convection term b/αx+βy.
(a)(b)


### 4.5. The Effect of Burner Size on Quenching

The relationship between combustion the chamber shape and quenching features is discussed in this paragraph, which shows how the ratio of length to width of the burner affects the quenching behavior. In terms ofcx,y=b/1+αx+βy, we also consider serial values of α/β taken in 0.7,0.97 besides α=β corresponding to Case 7.0 in Table 5 because it does not produce quenching behavior when α/β is smaller than 0.7. We choose ten values: 0.7, 0.72, 0.82, 0.85, 0.88, 0.89, 0.9, 0.92, 0.95, and 0.97 for Case 7.0.1–Case 7.0.10 in Table 5. Specially, Case 7.0.1: α/β=0.7, Case 7.0.2: α/β=0.72, Case 7.0.3: α/β=0.82, Case 7.0.4: α/β=0.85, Case 7.0.5: α/β=0.88, Case 7.0.6: α/β=0.89, Case 7.0.7: α/β=0.9, Case 7.0.8: α/β=0.92, Case 7.0.9: α/β=0.95, and Case 7.0.10: α/β=0.97. The quenching data associated with different α/β based on q=0.5,α=200,bx=by=10,θ=1, τ0=0.0003, and h0=1/55 for the convection term b/1+αx+βy are recorded in Table 5. From Table 5, we can discover that quenching time tmax that is 4.517239999869761 is the largest when α/β=0.85 and quenching location xmax=ymax that is 0.03636363636363636 when 0.89≤α/β≤1.Table 5
Quenching time, quenching location,maxu, and maxut relative to different α/β when cx,y=b/1+αx+βy and q=0.5,α=200,bx=by=10,andθ=1.
Case nametmaxxmaxymaxmaxumaxutCase 7.04.4877765796389270.036363636363636360.036363636363636360.99924094983179571539.153786381783Case 7.0.14.2796229436678040.036363636363636360.018181818181818180.99926817235246561710.18659619492Case 7.0.24.3318311716037690.036363636363636360.018181818181818180.99919038601226571552.306642881425Case 7.0.34.4965841665480550.036363636363636360.018181818181818180.99919878527678541599.733801339674Case 7.0.44.5172399998697610.036363636363636360.018181818181818180.99932564573152071911.903484162842Case 7.0.54.514654812100120.036363636363636360.018181818181818180.99922345142780731453.952639953373Case 7.0.64.5130262860982670.036363636363636360.036363636363636360.99959385434728312793.85091710218Case 7.0.74.5112082758536160.036363636363636360.036363636363636360.99976333747722144813.578815453525Case 7.0.84.507259249814290.036363636363636360.036363636363636360.99936721522188381807.556556363968Case 7.0.94.5005082887814940.036363636363636360.036363636363636360.99967514335532653556.703280522448Case 7.0.104.4956643737398980.036363636363636360.036363636363636360.99923969308800791524.670120016285Similarly, we investigate quenching behavior variation with the changes ofα/β for convection function b/αx+βy. From a number of tests, we select eleven meaningful cases to be given in Table 6. Case 8.0: α/β=1, Case 8.0.1: α/β=0.15, Case 8.0.2: α/β=0.25, Case 8.0.3: α/β=0.3, Case 8.0.4: α/β=0.4, Case 8.0.5: α/β=0.5, Case 8.0.6: α/β=0.6, Case 8.0.7: α/β=0.7, Case 8.0.8: α/β=0.8, Case 8.0.9: α/β=0.85, and Case 8.0.10: α/β=0.95. The initial t−step is τ0=0.0005, and the initial x−step (y−step) is h0=1/40. q, α, bxby, and θ are set as cq=1.1,α=20,bx=by=3,andθ=1. It is easily seen that quenching occurs, and quenching location xmax=ymax is 0.025 when α/β≥0.15 and quenching time tmax is the smallest when α/β=0.375 from Table 6.Table 6
Quenching time, quenching location,maxu, and maxut relative to different α/β when cx,y=b/αx+βy and q=1.1,α=20,bx=by=3,andθ=1.
Case nametmaxxmaxymaxmaxumaxutCase 8.05.6042627579092470.0250.0250.983534151987658773.15674241637485Case 8.0.148.532751622336880.0250.0250.980067370148566211.64671094996052Case 8.0.229.172570437063890.0250.0250.980617614577940320.45608182578629Case 8.0.325.293710733071440.0250.0250.981360578040261426.01159161851154Case 8.0.424.182158212581830.0250.0250.981979446412260934.70682386174796Case 8.0.524.677549002639620.0250.0250.981092748693986839.65241020973182Case 8.0.646.469330577447670.0250.0250.982587000051935751.20986827535954Case 8.0.714.790355402763990.0250.0250.980480825232632448.88835666486023Case 8.0.815.630820529256960.0250.0250.981019432588903854.83303601502359Case 8.0.922.864815498774960.0250.0250.980643384714718556.65280512266988Case 8.0.105.1602778087048190.0250.0250.982476989886559566.75815635770577Figures13 and 14 display the relationships between tmax, xmax, ymax, and α/β, respectively, in which, the red curves in Figures 13(a) and 14(a) reflect quenching time tmax variations as α/β rises; the blue curves in Figures 13(b) and 14(b) reflect quenching location xmax variations as α/β rises; the green curves in Figures 13(b) and 14(b) reflect quenching location ymax variations as α/β rises. Figure 13 corresponds to Table 5, and Figure 14 corresponds to Table 6. The highest value marked by the green sign of the red contour in Figure 13(a) is the quenching time 4.517239999869761 of Case 7.0.4. Quenching time tmax drops monotonously as α/β increases when 0.85<α/β≤1.0, and it promotes monotonously as α/β increases when 0.7≤α/β≤0.85. The shape of the red contour in Figure 14(a) is different from that of the red contour in Figure 13(a). Its lowest green point is just the quenching time 5.160277808704819 of Case 8.0.10, and α/β is 0.95. The contour shows the fluctuation of the points up and down. For the relationship between xmax and α/β, we rely on the tendencies of the red dots on the blue curves in Figures 13(b) and 14(b), respectively, to see that xmax keeps a same value in each case. Depending on the trend of the purple dots on the green curves in Figures 13(b) and 14(b), respectively, we can find that ymax keeps a smaller value than xmax when 0.7≤α/β≤0.89 and becomes the same value as xmax when 0.89≤α/β≤1 for Figure 13(b) and ymax is equal to xmax for Figure 14(b). Meanwhile, it is clearly observed that there is only the blue contour with ten red points extending to the end in Figure 14(b).Figure 13
(a) The plot oftmax as α/β increases and (b) the plots of xmax and ymax as α/β increases for the convection term b/1+αx+βy.
(a)(b)Figure 14
(a) The plot oftmax as α/β increases and (b) the plots of xmax and ymax as α/β increases for the convection term b/αx+βy.
(a)(b)Finally, we can conclude that quenching problems based on different convection functions require different combustion chamber sizes. It can lead to quenching behaviors to properly take the ratio ofα and β when α/β≥0.7 for convection function b/1+αx+βy. But, it can result in quenching phenomena to properly take the ratio of α and β when α/β≥0.15 for the convection term b/αx+βy. By comparing, it is easier to form quenching behaviors for the latter than the former. Generally, we can consider it will produce quenching phenomena when α/β takes a value not less than 0.7. There must be a suitable ratio that will spend the most quenching time. When the length and the width of the combustion chamber are close (i.e., 0.89≤α/β≤1), the horizontal coordinate of quenching location is equivalent to its vertical coordinate.


## 4.1. Common Example

It is not a quenching case. It is just for demonstrating the accuracy and convergence of the proposed difference scheme. We have solved the following unsteady convection-diffusion equation with coefficient whose exact solution is offered:(40)put=uxx+uyy−pux−puy,x,y∈0,1×0,1,t≥0,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.The initial and boundary value conditions are derived from the exact solution. There are four metrics for the data test, which are maximum error, average error, CPU time, and convergence rate. We adopt these parametersT=0.5,τn=0.01,andp=100 to fulfill the two difference schemes: the proposed scheme based on the uniform mesh and the proposed scheme based on the nonuniform mesh. k/K+λx/πsinπk/Kk=0K and j/J+λy/πsinπj/Jj=0J represent the nonuniform grid, where KandJ point to intervals of the spatial direction and K=J, whereas λxandλy point to the mesh transformation coefficient and λx=λy.Table1 offers these measuring data. It displays the accuracies and running time of our scheme (uniform and nonuniform). We select five spatial intervals: 16, 32, 64, 128, and 256 as test cases, which are given to each of the above two meshes (i.e., uniform and nonuniform). For the nonuniform mesh, the different values of λx impact these precision criteria. In the paper, λx takes 0.7. The maximum errors of the first scheme are higher than those of the second one, and the average errors of the latter are superior to those of the former. Furthermore, the larger the number of spatial intervals is, the higher the accuracies of the scheme based on the nonuniform grid are. Its accuracies arrive at 1.3364 × 10−7 in the maximum error and 6.4492 × 10−9 in the average error when K=J=256. In terms of computing time of the two schemes whose measurement unit is millisecond, the first scheme is lower than the second one. Obviously, by comparing with the uniform mesh, the difference scheme based on the nonuniform mesh spends more time to run. Lastly, the conv. rate values are observed. The tendency of the conv. rate from the former scheme is similar to that from the latter scheme. We take the latter as the example. The four values are positive, which means the accuracy increases as K rises when K∈16,256. The conv. rate value is maximum when K transfers from 16 to 32, which means the accuracy promotes greatly in the stage. The latter three values are close, which means the accuracy promotes gently in the three stages.Table 1
The parameter values for common example.
Difference schemesKJMax. errorAver. errorCPU timeConv. rateThe proposed scheme (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.1743 × 10−7117943.9714The proposed scheme (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.01432561.3364 × 10−76.4492 × 10−91659693.9867


## 4.2. Quenching Case without Convection Term

This section will describe a quenching status without the convection term, i.e.,c1x,y=0 and c2x,y=0. The mathematical model can be expressed as(41)α2x2+β2y2q/2ut=1α2uxx+1β2uyy+11−uθ,x,y∈0,1×0,1,t∈0,T.Its initial and boundary conditions are just as in equations (6) and (7).After equations (19) and (20) are used to approximate equation (41), we can get some quenching behaviors. For the 2D nonlinear degeneration singularity equation without the convection term, we demonstrate the quenching case of q=0 that has been described in [13, 14]. In Table 2 and Figure 1, q=0, a=β=10, and θ=1 which are the same for the three approaches, but the initial space step h0 and the initial temporal step τ0 are different among them, where we take h0=1/50 and τ0=4.166677×10−3. In Table 2, the five parameters represent the same senses, in which xmax,ymax represents the point of quenching location, tmax represents the quenching time, maxurepresents the maximum temperature value immediately before quenching happens, and maxut means the maximum variation of temperature with respect to time immediately before quenching happens. From Table 2, we can see that the quenching location 0.5,0.5 of the referred scheme is identical with that of the scheme in [14], and tmax=0.5958348109999995 and maxu=0.998562875689732 of the first one are close to those of the last two. Obviously, maxut is different among the various methods. From Figure 1, it is found that u⟶1− and ut⟶690−at the location xmax,ymax=0.5,0.5 and the quenching time tmax=0.5958348109999995.Table 2
Quenching performances of different schemes forq=0, a=β=10, and θ=1 without the convection term.
SchemestmaxxmaxymaxmaxumaxutThe proposed scheme0.59583481099999950.5000.5000.998562875689732689.2724781649539The scheme in [13]0.58712499993751——0.990432148.887767The scheme in [14]0.5875540.5000.5000.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 the convection term.
(a)(b)


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

The standard 2D degenerate semilinear quenching problem with the initial boundary condition lists in equations (5)–(7) is extended from the equation in one-space dimension [23]. The first convection term b/1+αx+βy can be written as two functions c1x,y=bx/1+αx+βy and c2x,y=by/1+αx+βy, where bx and by take the constants. The four groups of variables: q, bx and by, α and β, and θ in equations (5)–(7) affect the quenching characteristics of the 2D problem, whose meaning has been stated in Section 2. There follows a special case for b/1+αx+βy. The initial time step is considered as τ0=0.0003, and the initial space step is considered as h0=1/55, in which the x-direction is identical with the y-direction. The four groups of parameters are configured as q=0.5,α=β=200,bx=by=10,andθ=1 which is defined as Case 7.0. This will lead to a quenching phenomenon. Quenching location x7∗,y7∗ is 0.03636363636363636,0.03636363636363636; the quenching time T7∗ is 4.487776579638927; maxx,yu is 0.9992409498317957; maxx,yut is 1539.153786381783. When t is 4.486155872847072, the temporal step starts to become smaller and smaller until quenching occurs. We choose six rows of representative data given in Table 3 which poins to features relative to x, y, t, u, and ut before quenching occurs of Case 7.0. maxx,yu and maxx,yut refer to maximal u and maximal ut in each time layer before quenching occurs; x and y refer to the coordinate point of synchronously producing maxx,yu and maxx,yut in each time layer before quenching occurs; t refers to the moment at which each time layer is located before quenching occurs. x,y first appears at the point x7∗,y7∗ when t is 4.487088240486234 which is the first row in Table 3 and maintains at this point until quenching occurs. When t is T7∗ which is the immediate time before quenching happens, we can obtain important characteristics. When t is 4.487777138788559, maxx,yut blows up and reaches 1.051076326841211 × 109.Table 3
Time, location,maxx,yu, and maxx,yut in the six temporal layers which contain quenching occurrence for Case 7.0.
txymaxx,yumaxx,yut4.4870882404862340.036363636363636360.036363636363636360.979448569828961552.184431339763224.4874898808006030.036363636363636360.036363636363636360.985461462054538975.528997590273484.4876538240986460.036363636363636360.036363636363636360.9897483333199262109.07447404344314.4877707022550420.036363636363636360.036363636363636360.9972318816440022418.00237744504054.4877765796389270.036363636363636360.036363636363636360.99924094983179571539.1537863817834.487777138788559———1.051076326841211 × 109Figures2–7 provide seven groups of graphs to illustrate more particular quenching statuses. Figure 2 shows the 2D plots of maxx,yu and maxx,yut changing as t rises from the initial moment until quenching occurs. Figure 2(a) corresponds to u, and Figure 2(b) corresponds to ut. Figure 3 sketches the 2D contours of u and ut regarding x immediately before quenching happens. Figure 3(a) corresponds to u, and Figure 3(b) corresponds to ut. Figure 4 draws the 2D outlines of u and ut regarding y immediately before quenching happens. Figure 4(a) corresponds to u, and Figure 4(b) corresponds to ut. The three 2D plots in Figure 5 describe the adaptive mechanism before quenching occurs. The first one reflects spatial adaptation, and the last two reflect temporal adaptation. Figure 5(a) describes the trend of the temporal step over time. To distinguish more clearly the variation, we intercept the curves of the period from 4.470000000001407 to 4.487777138788559. The red sign on the blue curve in Figure 5(a) denotes the initial temporal adaptive moment (i.e., 4.486155872847072). The temporal step is uniform before the adaptive moment but becomes nonuniform when t varies from 4.486155872847072 to 4.487777138788559. The green sign on the blue curve in Figure 5(a) denotes the initial moment (i.e., 4.487088240486234) of producing the quenching location x7∗,y7∗. The green lines in Figures 5(b) and 5(c) describe the uniform spatial steps in the x−direction and the y−direction, respectively. The spatial adaptation occurs at 4.487777138788559 after which the spatial step becomes nonuniform. The red curves with blue spots in Figures 5(b) and 5(c) describe, respectively, the nonuniform spatial steps in the x−direction and the y−direction at T7∗. Although the two curves are very similar, there is still subtle difference between their data. Figures 6–7 provide the 3D surfaces of u and ut with respect to spatial variables x and y at certain special time. Figure 6(a) records the 3D surfaces about x, y, and u, and Figure 6(b) records the 3D surfaces about x, y, and ut when t=4.487770702255042. Figure 7(a) records the 3D surfaces about x, y, and u, and Figure 7(b) records the 3D surfaces about x, y, and ut when T7∗.Figure 2
The graphs of (a)maxx,yu as t increases and (b) maxx,yut as t increases from the initial moment to the occurrence of quenching behavior for the convection term b/1+αx+βy.
(a)(b)Figure 3
The graphs of (a)u as x increases and (b) ut as x increases when t is T7∗ for the convection term b/1+αx+βy.
(a)(b)Figure 4
The graphs of (a)u as y increases and (b) ut as y increases when t is T7∗ for the convection term b/1+αx+βy.
(a)(b)Figure 5
(a) The curve of the temporal stepsτn as time t increases, (b) the variation of spatial steps hk as x processes, and (c) the variation of spatial steps hk as y processes and for the convection term b/1+αx+βy.
(a)(b)(c)Figure 6
(a) The 3D plot ofu and (b) the 3D plot of ut when t=4.487770702255042 for the convection term b/1+αx+βy.
(a)(b)Figure 7
(a) The 3D plot ofu and (b) the 3D plot of ut when t=T7∗ for the convection term b/1+αx+βy.
(a)(b)


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

The second convection termb/αx+βy can be written as the two convection functions: c1x,y=bx/αx+βy and c2x,y=by/αx+βy, where bx and by take the constants. Comparing with the first convection term b/1+αx+βy, it is hard to conduct quenching behavior for the second convection term b/αx+βy. We make a following configuration to gain a quenching status. The initial x−step (y−step) is h0=1/40, and the initial t−step is τ0=0.0005. The four groups of important parameters are set as q=1.1,α=β=20,bx=by=3,andθ=1 which is defined as Case 8.0. The parameter statuses of the case are recorded in Table 4. For the case, the initial time of producing the quenching location is earlier than the temporal adaptive time. The former time is 5.596499999999618, which is the first row in Table 4, the latter one is 5.603982012140827, which is the fifth row in Table 4, and the quenching location x8∗,y8∗ is 0.025,0.025. Moreover, it is so transient from the front time of quenching location appearing to the terminative time of quenching location appearing that there are 18 time layers between the two moments. The latter quenching time is T8∗=5.604262757909247, maxx,yu is 0.9835341519876587, and maxx,yut is 73.15674241637485 immediately before quenching occurs, which is the sixth row in Table 4. maxx,yut gets the value 6.637517818582074 × 1011 and blows up when t is 5.604272047535358, which is the last row in Table 4.Table 4
Time, location,maxx,yu, and maxx,yut in the last seven temporal layers which contains quenching occurrence for Case 8.0.
txymaxx,yumaxx,yut5.5964999999996180.0250.0250.965613248035384527.204925805149065.5984999999996170.0250.0250.966848091050629728.766902187783995.5999999999996160.0250.0250.967807952885317130.064978162417135.6019999999996150.0250.0250.96913921620335732.000736898491995.6039820121408270.0250.0250.970836300971370934.727044766461555.6042627579092470.0250.0250.983534151987658773.156742416374855.604272047535358———6.637517818582074 × 1011This paragraph supplies five sets of figures (i.e., Figures8–12) to more vividly depict the quenching information of the case. The 2D curves in Figure 8 describe the relationship between maxx,yu and t and the relationship between maxx,yut and t. The 2D curves in Figure 9 describe the relationship between u and x and the relationship between ut and x. The 2D curves in Figure 10 describe the relationship between u and y and the relationship between ut and y. In Figures 8–10, (a) corresponds to u and (b) corresponds to ut. Figure 11(a) draws the tendency of the temporal step varying with temporal variable t; Figure 11(b) draws the tendency of the initial spatial step and the adaptive step varying with spatial variable x; and Figure 11(c) draws the tendency of the initial spatial step and the adaptive step varying with spatial variable y. For Figure 11(a), to distinguish more clearly the variation, we intercept the curves of the period from 5.580499999999627 to 5.604272047535358; and the red spot on the blue curve expresses the earliest temporal adaptive time (i.e., 5.603982012140827). maxx,yut in the time step becomes 1.185615675885128 × 1012 when t is 5.603982012140827. The temporal step is uniform before the adaptive moment but becomes nonuniform since then. The green spot on the blue curve expresses the initial moment of producing the quenching location x8∗,y8∗, which is 5.604262757909247. In Figures 11(b) and 11(c), the green lines indicate the uniform initial spatial steps while the red curves with blue signs indicate the nonuniform adaptive spatial steps. Lastly, we provide two 3D surfaces in Figure 12 to state quenching statuses of u and ut when t=T8∗. Figure 12(a) denotes the relationship between u, x, and y, and Figure 12(b) denotes the relationship between ut, x, and y. Figure 12(a) records the 3D surfaces about x, y, and u, and Figure 12(b) records the 3D surfaces about x, y, and ut when t=T8∗. There is a narrowband projecting the rear plane of the 3D axis box in Figure 12(b). This is because some larger values of maxx,yut appear at the positions, in which the x−axis and the y−axis are simultaneously near the origin of the coordinates. From the data source of the 3D surface, we can observe some larger maxx,yut varying from 1 to 27, whereas maxx,yut at other positions are positive numbers less than 1, which can be seen via Figures 9(b) and 10(b). Obviously, the larger maxx,yut are all less than 73.15674241637485, which is reasonable for a quenching case.Figure 8
The graphs of (a)maxx,yu as t increases and (b) maxx,yut as t increases from the initial moment to the occurrence of quenching behavior for the convection term b/αx+βy.
(a)(b)Figure 9
The graphs of (a)u as x increases and (b) ut as x increases when t is T8∗ for the convection term b/αx+βy.
(a)(b)Figure 10
The graphs of (a)u as y increases and (b) ut as y increases when t is T8∗ for the convection term b/αx+βy.
(a)(b)Figure 11
(a) The curve of the temporal stepsτn as time t increases, (b) the variation of spatial steps hk as x processes, and (c) the variation of spatial steps hk as y processes and for the convection term b/αx+βy.
(a)(b)(c)Figure 12
(a) The 3D plot ofu and (b) the 3D plot of ut when t=T8∗ for the convection term b/αx+βy.
(a)(b)


## 4.5. The Effect of Burner Size on Quenching

The relationship between combustion the chamber shape and quenching features is discussed in this paragraph, which shows how the ratio of length to width of the burner affects the quenching behavior. In terms ofcx,y=b/1+αx+βy, we also consider serial values of α/β taken in 0.7,0.97 besides α=β corresponding to Case 7.0 in Table 5 because it does not produce quenching behavior when α/β is smaller than 0.7. We choose ten values: 0.7, 0.72, 0.82, 0.85, 0.88, 0.89, 0.9, 0.92, 0.95, and 0.97 for Case 7.0.1–Case 7.0.10 in Table 5. Specially, Case 7.0.1: α/β=0.7, Case 7.0.2: α/β=0.72, Case 7.0.3: α/β=0.82, Case 7.0.4: α/β=0.85, Case 7.0.5: α/β=0.88, Case 7.0.6: α/β=0.89, Case 7.0.7: α/β=0.9, Case 7.0.8: α/β=0.92, Case 7.0.9: α/β=0.95, and Case 7.0.10: α/β=0.97. The quenching data associated with different α/β based on q=0.5,α=200,bx=by=10,θ=1, τ0=0.0003, and h0=1/55 for the convection term b/1+αx+βy are recorded in Table 5. From Table 5, we can discover that quenching time tmax that is 4.517239999869761 is the largest when α/β=0.85 and quenching location xmax=ymax that is 0.03636363636363636 when 0.89≤α/β≤1.Table 5
Quenching time, quenching location,maxu, and maxut relative to different α/β when cx,y=b/1+αx+βy and q=0.5,α=200,bx=by=10,andθ=1.
Case nametmaxxmaxymaxmaxumaxutCase 7.04.4877765796389270.036363636363636360.036363636363636360.99924094983179571539.153786381783Case 7.0.14.2796229436678040.036363636363636360.018181818181818180.99926817235246561710.18659619492Case 7.0.24.3318311716037690.036363636363636360.018181818181818180.99919038601226571552.306642881425Case 7.0.34.4965841665480550.036363636363636360.018181818181818180.99919878527678541599.733801339674Case 7.0.44.5172399998697610.036363636363636360.018181818181818180.99932564573152071911.903484162842Case 7.0.54.514654812100120.036363636363636360.018181818181818180.99922345142780731453.952639953373Case 7.0.64.5130262860982670.036363636363636360.036363636363636360.99959385434728312793.85091710218Case 7.0.74.5112082758536160.036363636363636360.036363636363636360.99976333747722144813.578815453525Case 7.0.84.507259249814290.036363636363636360.036363636363636360.99936721522188381807.556556363968Case 7.0.94.5005082887814940.036363636363636360.036363636363636360.99967514335532653556.703280522448Case 7.0.104.4956643737398980.036363636363636360.036363636363636360.99923969308800791524.670120016285Similarly, we investigate quenching behavior variation with the changes ofα/β for convection function b/αx+βy. From a number of tests, we select eleven meaningful cases to be given in Table 6. Case 8.0: α/β=1, Case 8.0.1: α/β=0.15, Case 8.0.2: α/β=0.25, Case 8.0.3: α/β=0.3, Case 8.0.4: α/β=0.4, Case 8.0.5: α/β=0.5, Case 8.0.6: α/β=0.6, Case 8.0.7: α/β=0.7, Case 8.0.8: α/β=0.8, Case 8.0.9: α/β=0.85, and Case 8.0.10: α/β=0.95. The initial t−step is τ0=0.0005, and the initial x−step (y−step) is h0=1/40. q, α, bxby, and θ are set as cq=1.1,α=20,bx=by=3,andθ=1. It is easily seen that quenching occurs, and quenching location xmax=ymax is 0.025 when α/β≥0.15 and quenching time tmax is the smallest when α/β=0.375 from Table 6.Table 6
Quenching time, quenching location,maxu, and maxut relative to different α/β when cx,y=b/αx+βy and q=1.1,α=20,bx=by=3,andθ=1.
Case nametmaxxmaxymaxmaxumaxutCase 8.05.6042627579092470.0250.0250.983534151987658773.15674241637485Case 8.0.148.532751622336880.0250.0250.980067370148566211.64671094996052Case 8.0.229.172570437063890.0250.0250.980617614577940320.45608182578629Case 8.0.325.293710733071440.0250.0250.981360578040261426.01159161851154Case 8.0.424.182158212581830.0250.0250.981979446412260934.70682386174796Case 8.0.524.677549002639620.0250.0250.981092748693986839.65241020973182Case 8.0.646.469330577447670.0250.0250.982587000051935751.20986827535954Case 8.0.714.790355402763990.0250.0250.980480825232632448.88835666486023Case 8.0.815.630820529256960.0250.0250.981019432588903854.83303601502359Case 8.0.922.864815498774960.0250.0250.980643384714718556.65280512266988Case 8.0.105.1602778087048190.0250.0250.982476989886559566.75815635770577Figures13 and 14 display the relationships between tmax, xmax, ymax, and α/β, respectively, in which, the red curves in Figures 13(a) and 14(a) reflect quenching time tmax variations as α/β rises; the blue curves in Figures 13(b) and 14(b) reflect quenching location xmax variations as α/β rises; the green curves in Figures 13(b) and 14(b) reflect quenching location ymax variations as α/β rises. Figure 13 corresponds to Table 5, and Figure 14 corresponds to Table 6. The highest value marked by the green sign of the red contour in Figure 13(a) is the quenching time 4.517239999869761 of Case 7.0.4. Quenching time tmax drops monotonously as α/β increases when 0.85<α/β≤1.0, and it promotes monotonously as α/β increases when 0.7≤α/β≤0.85. The shape of the red contour in Figure 14(a) is different from that of the red contour in Figure 13(a). Its lowest green point is just the quenching time 5.160277808704819 of Case 8.0.10, and α/β is 0.95. The contour shows the fluctuation of the points up and down. For the relationship between xmax and α/β, we rely on the tendencies of the red dots on the blue curves in Figures 13(b) and 14(b), respectively, to see that xmax keeps a same value in each case. Depending on the trend of the purple dots on the green curves in Figures 13(b) and 14(b), respectively, we can find that ymax keeps a smaller value than xmax when 0.7≤α/β≤0.89 and becomes the same value as xmax when 0.89≤α/β≤1 for Figure 13(b) and ymax is equal to xmax for Figure 14(b). Meanwhile, it is clearly observed that there is only the blue contour with ten red points extending to the end in Figure 14(b).Figure 13
(a) The plot oftmax as α/β increases and (b) the plots of xmax and ymax as α/β increases for the convection term b/1+αx+βy.
(a)(b)Figure 14
(a) The plot oftmax as α/β increases and (b) the plots of xmax and ymax as α/β increases for the convection term b/αx+βy.
(a)(b)Finally, we can conclude that quenching problems based on different convection functions require different combustion chamber sizes. It can lead to quenching behaviors to properly take the ratio ofα and β when α/β≥0.7 for convection function b/1+αx+βy. But, it can result in quenching phenomena to properly take the ratio of α and β when α/β≥0.15 for the convection term b/αx+βy. By comparing, it is easier to form quenching behaviors for the latter than the former. Generally, we can consider it will produce quenching phenomena when α/β takes a value not less than 0.7. There must be a suitable ratio that will spend the most quenching time. When the length and the width of the combustion chamber are close (i.e., 0.89≤α/β≤1), the horizontal coordinate of quenching location is equivalent to its vertical coordinate.


## 5. Conclusions

Nonlinear singular degenerate reaction-diffusion equation plays more and more important role in the quenching area such as diffusion of thermal and energy, combustion theory, aerospace engineering, and biomedical treatment. We put forward a compact difference scheme and built an adaptation algorithm to solve the 2D convection-reaction-diffusion equation of the quenching type. Depending on theoretical and numerical experiments above, we can get some conclusions as follows. As a valid and powerful mathematical tool, the high-order compact difference scheme combining with the proper adaptive procedure may offer elegant manners to resolve the semilinear singular degenerate partial differential equation. It is very efficient to employ the two-layer self-adaptive grid principle for solving the quenching problem. The self-adaptive temporal step makes time adaptation activated when matching the critical condition. It not only saves computational cost, but also impacts on specific moment of quenching. Similarly, it is important to choose spatial adaptation parameters, which may decide where the temporal derivative will blow-up, even whether quenching occurs or not. The study provides the feasibility of analyzing and resolving different engineering applications of the quenching problem in different finite difference schemes on an adaptive mesh. We will conduct further research in the future, which includes the 3D reaction-diffusion equations with the convection function of the quenching type.



---

*Source: 8161804-2020-06-15.xml*