# Second-Order Semi-Discretized Schemes for Solving Stochastic Quenching Models on Arbitrary Spatial Grids

Authors: Nina Garcia-Montoya; Julienne Kabre; Jorge E. Macías-Díaz; Qin Sheng
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

Reaction-diffusion-advection equations provide precise interpretations for many important phenomena in complex interactions between natural and artificial systems. This paper studies second-order semi-discretizations for the numerical solution of reaction-diffusion-advection equations modeling quenching types of singularities occurring in numerous applications. Our investigations particularly focus at cases where nonuniform spatial grids are utilized. Detailed derivations and analysis are accomplished. Easy-to-use and highly effective second-order schemes are acquired. Computational experiments are presented to illustrate our results as well as to demonstrate the viability and capability of the new methods for solving singular quenching problems on arbitrary grid platforms.


---

## Body

## 1. Introduction

Nonlinear reaction-diffusion-advection equations have been playing an important role in interactions between natural and artificial systems. The partial differential equations provide precisely mathematical interpretations for numerous natural phenomena, such as diffusion of heat and energy, burning or quenching of fuels, energy concentrations, and transformations in different environments. Nonlinear partial differential equations also deliver popular computational tools to cell biology and cancer treatment plans [1–4].Consider anm-dimensional heat engine, m≥1. Let u be the temperature distribution inside of its combustion chamber and u0 be the initial temperature distribution due to the sparks. Assume that combustion occurs at the unit temperature. Further, let 0<σ≤1,t>0, be the reciprocal chamber size index. The use of variable σ reflects the change of chamber size, which is typical when different chemical fuels are selected. Then, an idealized combustion model can be comprised as a nonlinear reaction-diffusion-advection initial-boundary value problem:(1)ut=σ2∇a∇u+fuXu,x∈D,t>0,(2)ux,t=0,x∈∂D,t>0,(3)ux,=0=u0x,x∈D∪∂D,where D⊂ℝm, ∂D is its boundary, ∇ is the m-dimensional gradient vector, the coefficient a=ax∈C1D is positive, 0≤u0≪1, and Xu is the stochastic indicator for u<1 [5, 6].The nonlinear source functionfu>0 is strictly increasing with respect to 0≤u<1 and(4)f0>0,limu⟶1−fu=∞.A typical example of such reaction term is(5)fu=11−uκ,for0≤u<1,where κ>0 is the internal combustion index utilized [3, 5–7].The solutionu of (1)–(4) is said to quench if there exists a finite Tσ such that(6)sup0<x<1utx,t⟶∞,ast⟶Tσ−.Such a valueTσ is called the quenching time which is a possible combustion-reaction time [7, 8]. It has been shown that a necessary condition for quenching to occur is(7)max0≤x≤1ux,t⟶1−,ast⟶Tσ−.The study of singular reaction-diffusion-advection equations of the form (1) and (5) can be traced back to Hale and Kawarada’s pioneering work [9, 10]. It was observed that when m=κ=1 and a≡1, there exists a critical valueσ∗≈0.65340 such that the solution of (1)–(3) and (5) quenches if σ<σ∗. Verifications of the existence and uniqueness to more general quenching models can be found in numerous recent publications including [1, 5, 9, 11].To solve the nonlinear differential equation problem (1)–(3) numerically is an interesting, yet challenging, multitask since we often do not know if the solution of (1)–(3) will quench or not until a proper solution procedure is taking place. In other words, in addition to approximating the solution u, a numerical method must be capable of evaluating simultaneously correct critical value σ∗, quenching time Tσ∗, quenching location x∗, and extremely sensitive and possibly unbounded derivative function ut [5, 6, 11]. This can be difficult since values of maxx∈D¯ut remain bounded and well-behaved until time t reaches a certain neighborhood of T∗, if it exists. Furthermore, we may also observe from properties (6) and (7) that, while maxx∈D¯ut grows exponentially, or faster than exponentially, in the aforementioned neighborhoods, the solution u itself continuously grows but stays bounded throughout the computation until a quenching suddenly erupts [7, 9, 10]. These coexisted and very distinct characters make it extremely hard to design an effective scheme.Due to their extremely important applications in nature and societies, numerous investigations, in both theory and computations, have been carried out for nonlinear quenching problems including (1)–(3) in recent years. Most of the existing algorithms are constructed either based on the reduced problems, that is, the stationary problems by removing the variable t, or by using fixed mesh steps (cf. [4, 7, 8, 12–14] and references therein). In these procedures, critical values, quenching times, solution u, and the rate-of-change function ut are often approximated incorrectly. Computational procedures are less efficient and less reliable, in particular when multidimensional modeling equations are considered. Although adaptive strategies have been particularly in favor due to their great flexibility and geometric accuracy in capturing quenching singularities involved, orders of accuracies of the spatial discretization of existing adaptive methods have been limited to one due to the unpredictability of the nonuniform grids generated by adaptations [6, 8, 12].On the other hand, effective new approaches in high-order approximations via finite difference, spectral, or finite element methods, have been introduced for solving other nonlinear partial differential equations [2, 15, 16]. Investigations have also been extended to fractional order partial differential equation problems [16, 17]. Since it is the high dimension that causes a tremendous increase in the computational cost, modern splitting techniques capable of converting higher dimensional problems into sets of single dimensional subproblems have reached a new height [13, 18]. Rigorous numerical analysis has also been given on split adaptations [18]. These have motivated our study of highly applicable and effective higher-order finite difference methods.Note that problems (1)–(3) can be numerically treated through proper second-order stable dimensional splitting [18, 19]. The decomposition may effectively reduce the total number of operations from On3m to On3, where m is the number of dimensions and n is the number of internal mesh points used [7, 18, 20]. Thus, in this paper, we shall focus on second-order semi-discretization for a one-dimensional modeling equation with m=1 and a≡1. In the circumstance, problems (1)–(3) can be simplified to(8)ut=σ2uxx+fuXu,0<x<1,t>0,(9)u0,t=u1,t=0,t>0,(10)ux,0=u0x,0≤x≤1.We consider an arbitrary partition of the spatial domain0,1, that is, grids(11)Ωn=0=x0,x1,x2,…,xk−1,xk,xk+1,…,xn,xn+1=1,for which 0<hk=xk+1−xk≪1,k=0,1,…,n. Denote Hn=h1,h2,…,hn+1. Apparently, such sets Ωn,Hn are t-dependent when they are used for solving (8)–(10). However, due to the feature of semi-discretization, we prefer dropping time location indicators for the simplicity of notations. We further assume that in general hk+1≠hk,k=0,1,…,n. At any spatial location xk, we adopt usual notations ϕxk,t=ϕkt or ϕxk,t=ϕk without confusion.


## 2. Semi-Discretization Approximations

At the time levelt>0, we set(12)hk=maxhk−3,hk−2,hk−1,hk,hk+1,hk+2,for any available index k. Without loss of generality, we drop the superindex for simplicity. We have expansions of the spatial derivatives(13)uk″=bk,−1uk−1+bk,0uk+bk,1uk+1+bk,2uk+2+bk,3uk+3+Oh2,(14)uk″=ak,−2uk−2+ak,−1uk−1+ak,0uk+ak,1uk+1+ak,2uk+2+Oh2,(15)uk″=ck,−3uk−3+ck,−2uk−2+ck,−1uk−1+ck,0uk+ck,1uk+1+Oh2.Strategy 1.
From the above, we must have(16)u1″=b1,−1u0+b1,0u1+b1,1u2+b1,2u3+b1,3u4+Oh2,uk″=ak,−2uk−2+ak,−1uk−1+ak,0uk+ak,1uk+1+ak,2uk+2+Oh2,k=2,…,n−1,un″=cn,−3un−3+cn,−2un−2+cn,−1un−1+cn,0un+cn,1un+1+Oh2.
Note thatu0=un+1=0 due to the homogeneous condition (9). Substituting the above into (8)–(10), we acquire that(17)ut=σ2Au+b+Oh2,u0=u0,where(18)A=b1,0b1,1b1,2b1,3a2,−1a2,0a2,1a2,2a3,−2a3,−1a3,0a3,1a3,2a4,−2a4,−1a4,0a4,1a4,2⋱⋱⋱⋱⋱an−2,−2an−2,−1an−2,0an−2,1an−2,2an−1,−2an−1,−1an−1,0an−1,1cn,−3cn,−2cn,−1cn,0∈ℝn×n,u=u1u2u3⋮⋮un∈ℝn,b=fu1Xu1fu2Xu2fu3Xu3⋮⋮funXun∈ℝn.
The band matrixA is large in size since n is usually large. But it is relatively sparse and can be handled conveniently by many existing software packages such as “sparse” subroutines in MATLAB. Dropping the truncation error term in (19), we obtain a second-order semi-discretization scheme for solving (8)–(10):(19)ut=σ2Au+b,(20)u0=u0.
Therefore, the solution of (8)–(10) can be readily approximated numerically through the solution of the nonlinear system of ordinary differential equations (19) and (20).Strategy 2.
The coefficient matrixA in (19) can be extremely stiff particularly due to the use of one-sided finite difference formulas (13) and (15). To avoid using the formulas, we may assume first that, in addition to the boundary values u0=un+1=0, we also have values u−1,un+2 calculated via some other methods that are at least of second order in accuracy. In this way, we may only use (14) for our second derivative approximations. The fully centralized formulas may hopefully ease the stiffness [19]. This yields our semi-discretized new approximation(21)ut=σ2Bu+c+Oh2,u0=u0,where(22)B=a1,0a1,1a1,2a2,−1a2,0a2,1a2,2a3,−2a3,−1a3,0a3,1a3,2a4,−2a4,−1a4,0a4,1a4,2⋱⋱⋱⋱⋱an−2,−2an−2,−1an−2,0an−2,1an−2,2an−1,−2an−1,−1an−1,0an−1,1an,−2an,−1an,0∈ℝn×n,u=u1u2u3⋮⋮un∈ℝn,c=fu1Xu1−a1,−2u−1fu2Xu2fu3Xu3⋮⋮funXun−an,2un+2∈ℝn.
Dropping the truncation error term from the differential system, we obtain another second-order semi-discretization scheme for solving (8)–(10):(23)ut=σ2Bu+c,(24)u0=u0.
Therefore, the solution of (8)–(10) can be approximated numerically through solving the nonlinear system of ordinary differential equations (23) and (24), given that values of u−1, un+2 can indeed be calculated successfully.


## 3. Parameter Determinations

Now, let us investigate the issue and see if valuesu−1,un+2 can be computed through some direct expansions. To this end, we have(25)uk+1=uk+hkuk′+hk22!uk″+hk33!uk‴+hk44!uk4+hk55!uk5+hk66!uk6+Ohk7.We observe readily that(26)u−1=u0−h−1u0′+h−122!u0″−h−133!u0‴+h−144!u04−h−155!u05+h−166!u06+Oh−17,un+2=un+1+hn+1un+1′+hn+122!un+1″+hn+133!un+1‴+hn+144!un+14+hn+155!un+15+hn+166!un+16+Ohn+17.For the simplicity in calculations, we may seth−1=h0,hn+1=hn. Thus,(27)u−1=u0−h0u0′+h022!u0″−h033!u0‴+h044!u04−h055!u05+h066!u06+Oh07,un+2=un+hnun′+hn22!un″+hn33!un‴+hn44!un4+hn55!un5+hn66!un6+Ohn7.To use the above two formulations for a quadratic order approximation, we need to balance coefficients related tou0′,u0″ and un′,un″ at internal spatial mesh points. Recall (13) and (15) and set(28)u0′=β1,0u0+β1,1u1+β1,2u2+Oh2,(29)u0″=β2,0u0+β2,1u1+β2,2u2+β2,3u3+β2,4u4+Oh2,(30)un′=γ1,n−2un−2+γ1,n−1un−1+γ1,nun+Oh2,(31)un″=γ2,n−4un−4+γ2,n−3un−3+γ2,n−2un−2+γ2,n−1un−1+γ2,nun+Oh2.Recall (25). To find the coefficients β1,0,β1,1,β1,2, we let(32)u0′=β1,0u0+β1,1u0+h0u0′+h022!u0″+h033!u0′′′+h044!u04+h055!u05+h066!u06+⋯+β1,2u0+h0+h1u0′+h0+h122!u0″+h0+h133!u0′+h0+h144!u04+h0+h155!u05+h0+h166!u06+⋯+Oh2.It follows therefore(33)β1,0+β1,1+β1,2=0,β1,1h0+β1,2h0+h1=1,β1,1h02+β1,2h0+h12=0.It turns out that(34)β1,0=−β1,1−β1,2,1=β1,1h0+β1,2h0+h1,β1,1=−h0+h12h02β1,2.From the last two equations,(35)−h0+h12h02β1,2h0+β1,2h0+h1=h0+h1−h0+h12h0β1,2=1.Hence,(36)β1,1=h0+h1h0h1,β1,2=−h0h0+h1h1.Consequently,(37)β1,0=−h0+h1h0h1+h0h0+h1h1=1h1h02−h0+h12h0h0+h1=−2h0+h1h0h0+h1.Similarly, to find coefficientsγ1,n−2,γ1,n−1,γ1,n, we consider(38)un′=γ1,n−2un−hn−2+hn−1un′+hn−2+hn−122!un″−hn−2+hn−133!un″′+hn−2+hn−144!un4−hn−2+hn−155!un5+hn−2+hn−166!un6+⋯+γ1,n−1un−hn−1u0′+hn−122!un″−hn−133!un″′+hn−144!un4−hn−155!un5+hn−166!un6+⋯+γ1,nun+Oh2.Balancing the coefficients, we have(39)γ1,n−2+γ1,n−1+γ1,n=0,−γ1,n−2hn−2+hn−1−γ1,n−1hn−1=1,γ1,n−2hn−2+hn−12+γ1,n−1hn−12=0.It follows subsequently that(40)γ1,n=−γ1,n−2−γ1,n−1,−1=γ1,n−2hn−2+hn−1+γ1,n−1hn−1,γ1,n−1=−hn−2+hn−12hn−12γ1,n−2.From the last two equations,(41)γ1,n−2hn−2+hn−1−hn−2+hn−12hn−12γ1,n−2hn−1=−hn−2+hn−1hn−2hn−1γ1,n−2=−1.Hence,(42)γ1,n−2=hn−1hn−2+hn−1hn−2,γ1,n−1=−hn−2+hn−1hn−2hn−1,and furthermore,(43)γ1,n=−hn−1hn−2+hn−1hn−2+hn−2+hn−1hn−2hn−1=hn−2+2hn−1hn−2+hn−1hn−1.Now, recall the expansion (25). To determine coefficients β2,0,β2,1,β2,2,β2,3,β2,4 in (29), we have(44)u0″=β2,0u0+β2,1u0+h0u0′+h022!u0″+h033!u0″′+h044!u04+h055!u05+h066!u06+⋯+β2,2u0+h0+h1u0′+h0+h122!u0″+h0+h133!u0″′+h0+h144!u04+h0+h155!u05+h0+h166!u06+⋯+β2,3u0+h0+h1+h2u0′+h0+h1+h222!u0″+h0+h1+h233!u0″′+h0+h1+h244!u04+h0+h1+h255!u05+h0+h1+h266!u06+⋯+β2,4u0+h0+h1+h2+h3u0′+h0+h1+h2+h322!u0″+h0+h1+h2+h333!u0″′+h0+h1+h2+h344!u04+h0+h1+h2+h355!u05+h0+h1+h2+h366!u06+⋯+Oh2.From the above, we obtain the following linear system:(45)β2,0+β2,1+β2,2+β2,3+β2,4=0,β2,1h0+β2,2h0+h1+β2,3h0+h1+h2+β2,4h0+h1+h2+h3=0,β2,1h02+β2,2h0+h12+β2,3h0+h1+h22+β2,4h0+h1+h2+h32=2,β2,1h03+β2,2h0+h13+β2,3h0+h1+h23+β2,4h0+h1+h2+h33=0,β2,1h04+β2,2h0+h14+β2,3h0+h1+h24+β2,4h0+h1+h2+h34=0.The system can be conveniently solved digitally via any existing solution package. By the same token, for evaluating coefficientsγ2,n−4,γ2,n−3,γ2,n−2,γ2,n−1,γ2,n, in (31), we find that(46)un″=+γ2,n−4un−hn−4+hn−3+hn−2+hn−1un′+hn−4+hn−3+hn−2+hn−122!un″−hn−4+hn−3+hn−2+hn−133!un′′′+hn−4+hn−3+hn−2+hn−144!un4−hn−4+hn−3+hn−2+hn−155!un5+hn−4+hn−3+hn−2+hn−166!un6+⋯+γ2,n−3un−hn−3+hn−2+hn−1un′+hn−3+hn−2+hn−122!un″−hn−3+hn−2+hn−133!un″′+hn−3+hn−2+hn−144!un4−hn−3+hn−2+hn−155!un5+hn−3+hn−2+hn−166!un6+⋯+γ2,n−2un−hn−2+hn−1un′+hn−2+hn−122!un″−hn−2+hn−133!un′′′+hn−2+hn−144!un4−hn−2+hn−155!un5+hn−2+hn−166!un6+⋯+γ2,n−1un−hn−1u0′+hn−122!un″−hn−133!un″′+hn−144!un4−hn−155!un5+hn−166!un6+⋯+γ2,nun+Oh2.The above expansion leads to the following non-homogeneous linear system:(47)γ2,n−4+γ2,n−3+γ2,n−2+γ2,n−1+γ2,n=0,γ2,n−4hn−4+hn−3+hn−2+hn−1+γ2,n−3hn−3+hn−2+hn−1+γ2,n−2hn−2+hn−1+γ2,n−1hn−1=0,γ2,n−4hn−4+hn−3+hn−2+hn−12+γ2,n−3hn−3+hn−2+hn−12+γ2,n−2hn−2+hn−12+γ2,n−1hn−12=2,γ2,n−4hn−4+hn−3+hn−2+hn−13+γ2,n−3hn−3+hn−2+hn−13+γ2,n−2hn−2+hn−13+γ2,n−1hn−13=0,γ2,n−4hn−4+hn−3+hn−2+hn−14+γ2,n−3hn−3+hn−2+hn−14+γ2,n−2hn−2+hn−14+γ2,n−1hn−14=0.Now, recall (14). To determine coefficients ak−2,ak−1,ak,ak+1,ak+2, we have(48)uk″=ak−2+ak−1+ak+ak+1+ak+2uk+−hk−2+hk−1ak−2−hk−1ak−1+hkak+1+hk+hk+1ak+2uk′+hk−2+hk−122!ak−2+hk−122!ak−1+hk22!ak+1+hk+hk+122!ak+2uk″+−hk−2+hk−133!ak−2−hk−133!ak−1+hk33!ak+1+hk+hk+133!ak+2uk″′+hk−2+hk−144!ak−2+hk−144!ak−1+hk44!ak+1+hk+hk+144!ak+2uk4+−hk−2+hk−155!ak−2−hk−155!ak−1+hk55!ak+1+hk+hk+155!ak+2uk5+hk−2+hk−166!ak−2+hk−166!ak−1+hk66!ak+1+hk+hk+166!ak+2uk6+Oh7.Therefore, to guarantee overallsecond-order approximations of the required function values for k=1,2,…,n, we must ask that(49)ak−2+ak−1+ak+ak+1+ak+2=0,−hk−2+hk−1ak−2−hk−1ak−1+hkak+1+hk+hk+1ak+2=0,hk−2+hk−122!ak−2+hk−122!ak−1+hk22!ak+1+hk+hk+122!ak+2=1,−hk−2+hk−133!ak−2−hk−133!ak−1+hk33!ak+1+hk+hk+133!ak+2=0,hk−2+hk−144!ak−2+hk−144!ak−1+hk44!ak+1+hk+hk+144!ak+2=0,which can be compressed into a matrix form(50)Mka=f,where(51)Mk=11111−hk−2−hk−1−hk−10hkhk+hk+1hk−2+hk−12hk−120hk2hk+hk+12−hk−2+hk−13−hk−130hk3hk+hk+13hk−2+hk−14hk−140hk4hk+hk+14,a=ak−2ak−1akak+1ak+2,f=00200,k=1,2,…,n.Furthermore, for (13), we observe that(52)uk″=bk−1uk−hk−1uk′+hk−122!uk″−hk−133!uk′′′+hk−144!uk4−hk−155!uk5+hk−166!uk6+bkuk+bk+1uk+hkuk′+hk22!uk″+hk33!uk′′′+hk44!uk4+hk55!uk5+hk66!uk6+bk+2uk+hk+hk+1uk′+hk+hk+122!uk″+hk+hk+133!uk′′′+hk+hk+144!uk4+hk+hk+155!uk5+hk+hk+166!uk6+bk+3uk+hk+hk+1+hk+2uk′+hk+hk+1+hk+222!uk″+hk+hk+1+hk+233!uk′′′+hk+hk+1+hk+244!uk4+hk+hk+1+hk+255!uk5+hk+hk+1+hk+266!uk6+Oh7=bk−1+bk+bk+1+bk+2+bk+3uk+−hk−1bk−1+hkbk+1+hk+hk+1bk+2+hk+hk+1+hk+2bk+3uk′+hk−122!bk−1+hk22!bk+1+hk+hk+122!bk+2+hk+hk+1+hk+222!bk+3uk″+−hk−133!bk−1+hk33!bk+1+hk+hk+133!bk+2+hk+hk+1+hk+233!bk+3uk″′+hk−144!bk−1+hk44!bk+1+hk+hk+144!bk+2+hk+hk+1+hk+244!bk+3uk4+−hk−155!bk−1+hk55!bk+1+hk+hk+155!bk+2+hk+hk+1+hk+255!bk+3uk5+hk−166!bk−1+hk66!bk+1+hk+hk+166!bk+2+hk+hk+1+hk+266!bk+3uk6+Oh7.Therefore, the result is as follows:(53)bk−1+bk+bk+1+bk+2+bk+3=0,−hk−1bk−1+hkbk+1+hk+hk+1bk+2+hk+hk+1+hk+2bk+3=0,hk−12bk−1+hk2bk+1+hk+hk+12bk+2+hk+hk+1+hk+22bk+3=2,−hk−13bk−1+hk3bk+1+hk+hk+13bk+2+hk+hk+1+hk+23bk+3=0,hk−14bk−1+hk4bk+1+hk+hk+14bk+2+hk+hk+1+hk+24bk+3=0.It can be rewritten as(54)Mlb=g,where(55)Ml=11111−hk−10hkhk+hk+1hk+hk+1+hk+2hk−120hk2hk+hk+12hk+hk+1+hk+22−hk−130hk3hk+hk+13hk+hk+1+hk+23hk−140hk4hk+hk+14hk+hk+1+hk+24,b=bk−1bkbk+1bk+2bk+3,g=00200.Now, for (15), we have(56)uk′=ck−3+ck−2+ck−1+ck+ck+1uk+−hk−3+hk−2+hk−1ck−3−hk−2+hk−1ck−2−hk−1ck−1+hkck+1uk′+hk−3+hk−2+hk−122!ck−3+hk−2+hk−122!ck−2+hk−122!ck−1+hk22!ck+1uk″+−hk−3+hk−2+hk−133!ck−3−hk−2+hk−133!ck−2−hk−133!ck−1+hk33!ck+1uk‴+hk−3+hk−2+hk−144!ck−3+hk−2+hk−144!ck−2+hk−144!ck−1+hk44!ck+1uk4+hk−3+hk−2+hk−155!ck−3−hk−2+hk−155!ck−2−hk−155!ck−1+hk55!ck+1uk5+hk−3+hk−2+hk−166!ck−3+hk−2+hk−166!ck−2+hk−166!ck−1+hk66!ck+1uk6+Oh7.This gives us the linear system(57)Mrc=r,where(58)Mr=11111−hk−3−hk−2−hk−1−hk−2−hk−1−hk−10hkhk−3+hk−2+hk−12hk−2+hk−12hk−120hk2−hk−3+hk−2+hk−13−hk−2+hk−13−hk−130hk3hk−3+hk−2+hk−14hk−2+hk−14hk−140hk4,c=ck−3ck−2ck−1ckck+1,r=00200.Based on the above investigations, we acquire the following result.Theorem 1.
If matricesMk,Ml, and Mr are nonsingular based on the selection of Hn, then the coefficient B in (23) is uniquely determined.We note that properties, such as the spectrums, of matricesMk,Ml,Mr and A,B still remain to be studied based on arbitrary sets of Hn at each temporal level of t. Bear in mind that nonuniform mesh steps in Hn are determined through particular adaptive procedure employed, such as the arc-length monitoring functions investigated in [7, 13, 17, 19].


## 4. Simulation Experiments

Letm=2 and κ=1. We consider two typical two-dimensional reaction-diffusion-advection modeling problems of (1)–(3) associated with (5).Example 1.
We are interested in natural quenching situations [21] with σ=5 and(59)ax,y=exp−10x−122+y−122,0≤x,y≤1,u0x,y=1401−cos2πx1−cos4πy,0≤x,y≤1.
We further setXu≡1 for the simplicity in illustrations. This problem can be viewed as a simplification of the internal combustion when dual ignition sparks are utilized [2, 3]. The partial differential equation problem is first split to one-dimensional subproblems via exponential splitting and then Strategies 1 and 2 are applied, respectively [8, 14, 18].
The semi-discretization mesh density controller is chosen to ben=200, while the temporal discretization parameter τ=λh is fixed, where λ=1/10>0 is the Courant-Friedrichs-Lewy constant, and h=maxkhk is taken in an appropriate temporal level, since temporal adaptations are not a goal of the current study. Standard exponentially evolving grid adaption is adopted [6, 12]. Cubic splines are used for passing data on each temporal level since spatial grids are moved. To see more clearly solution profiles, surface plots are only shown on internal mesh points without the fixed homogeneous boundary data (2).
It is observed that, due to the strong smoothness effect of the diffusion operator, the initial temperature distributionu0 due to the dual ignition quickly converges to a smooth surface as time increases. The surface further forms a symmetric distribution with respect to the center of the spatial domain. Effects of the two sparks become hardly noticeable. The phenomenon is demonstrated in Figure 1 at the second, 10,000th, 20,000th, and 30,000th temporal levels, respectively. It has also been found that the numerical solution is monotonically increasing pointwise as t increases.
Another extremely interesting physical function to monitor in quenching-combustion situations is the temporal derivative, or rate of change function,ut [7, 21]. The evolution of ut accelerates rapidly as time t is approaching the quenching time T∗, which is approximately 69727τ≈0.47624479. To see more details, we show the rate of change function ut in four temporal levels: T=69600τ,69700τ,69720τ, and 69726τ under the same scale in Figure 2. A rapid increase of the peak values can be observed. These consistently match the latest observations reported in [7, 21].
To see more precisely profiles ofu and ut, we plot evolutional trajectories of the maximal values and mean values of them for T∈0,T∗ in Figure 3. It can be seen that the maximum of the temperature field derivative ut increases exponentially as t⟶T∗, while the solution quenches peacefully. The rapid increase of averaging temperature in the combustion chamber is particularly important to engineering applications, since it indicates the success of the quick release of the chemical energy from fuel [3]. A y-semi-logarithmic scale is used to show greater details of the phenomenon, especially for the derivative function max0≤x,y≤1utx,y,T.
In Figure4, three-dimensional surfaces of u and ut, together with their x-u/ut plane projections, immediately before the quench are given. It can be seen that while the peak temperature, which is at the center of the combustion chamber, gradually approaches the fuel ignition, distribution of u tends to occupy the entire space of the chamber. Original influence of the dual sparks disappears. On the other hand, however, rates of the temperature change, ut, increase exponentially to the infinity at the center. The semi-discretized scheme is highly stable and captures the phenomenon satisfactorily in about seventy thousand time steps.Figure 1
Numerical solutionux,y,T at T=2τ≈1.36602691×10−5 (a); T=10,000τ≈0.06830134 (b); T=20,000τ≈0.13660269 (c); and T=30,000τ≈0.20490403 (d). The corresponding maximal values are 0.09996733, 0.13774652, 0.20421692, and 0.28487089, respectively.
(a)(b)(c)(d)Figure 2
Temporal derivative function, or rate of temperature change,utx,y,T prior to quenching at T=69,600τ≈0.47537736 (a); T=69,700τ≈0.47606037 (b); T=69,720τ≈0.47619698 (c); and T=69,726τ≈0.47623796 (d). The corresponding maximal values are 22.39748639, 43.58288207, 65.42874536, and 82.82435985, respectively.
(a)(b)(c)(d)Figure 3
Trajectories of the maximal (blue curves), mean values (red dotted curves) ofux,y,T,0≤x,y≤1, (a); and the maximal (blue curves), mean values (red dotted curves) of utx,y,T,0≤x,y≤1, (b) for T∈0,69725τ. Explosive increments of the maximal derivative function values as t approaches the quenching time 68725τ≈0.47623113<T∗≈69727τ≈0.47624479 are again visible. This is consistent with quenching criteria (6) and (7).
(a)(b)Figure 4
Numerical solutionux,y,T (a, b) and its temporal rate of change function utx,t (c, d), together with their respective x−u/ut plane projections, are given at the reference T≈69727τ≈0.47624479 immediately before quenching. We have max0≤x,y≤1ux,y,T≈98952207, max0≤x,y≤1utx,y,T≈87.32452308. Again, plots are over internal mesh points without showing the boundary data (2).
(a)(b)(c)(d)Example 2.
As a highly realistic extension, we consider another favorable natural quenching case withσ=5 and a much simplified four ignition spark initial condition [2, 3]:(60)u0x,y=1401−cos4πx1−cos4πy,0≤x,y≤1.
We continue settingXu≡1. Again, we consider corresponding one-dimensional subproblems via exponential splitting and then adopting Strategies 1 and 2, respectively. The same spatial mesh density controller and Courant-Friedrichs-Lewy constant are used for most of the computations until t is extremely close to the quenching time. Standard exponentially evolving grid adaption is adopted in the final stages of solution calculations [12]. Cubic splines are again used for passing data on each temporal level since spatial grids are moved. We only plot solutions on internal mesh points without showing the homogeneous boundary data (2).
Figure5 is devoted to four key moments of the solution procedures. Distributions of the temperature u and its temporal derivative ut inside the rectangular combustion chamber are simulated. The evidence of initial multiple sparks is obvious at T≈0.00683013. However, the influence apparently becomes fuzzy at T≈0.01366026 and fuzzier at T≈0.48121712, while the maximal temperature, max0≤x,y≤1ux,y,T, increases monotonically from 0.13965273 to 0.20441452 and then 0.97828162. We also notice that the maximal value locations move from the initial four sparks to the central of the spatial domain. The phenomenon can be seen more clearly in the last figure of u, for which T≈0.48142203. This is further supported by the corresponding set of ut surfaces in Figure 5. The facts indicate a unique quenching/blow-up position in the combustion chamber which is physically correct [4, 5, 11].
In the final phase atT≈0.48142203, values of u≈1, which is the singular point of the source function fu, in almost the entire combustion domain. On the other hand, we also observe that, while the distribution profile of ut is similar to that of u, the amplitude of rate function ut jumps to approximately 1.5×105. This is an indication that, during the final ignition, the fuel combustion speed must be tremendously high, while the dimensionless temperature stays under the unity.
We show both the maximal and mean values of the numerical solutionu and its temporal derivative ut,x,y∈D¯, in Figure 6. It is interesting that, while both mean values increase monotonically, the maximal derivative function decreases slightly in the beginning but quickly picks up the momentum and increases rapidly and becomes unbounded as t⟶T∗. However, both maximal values remain positive throughout computations. The simulation features well agree with the monotone solution property predicted in the theory of thermal combustion [2, 3].
Finally, in Figure7, profiles of the solution u and its temporal derivative ut are given at y=0.5,T=70483τ≈0.48140837 immediately before quenching. It can be seen that while the temperature inside the chamber increases uniformly towards the unity, the rate-of-change function ut, in other words, the velocity of the temperature field evolution, grows rapidly to the infinite at the quenching location x∗,y∗=0.5,0.5∈D.Figure 5
Numerical solutionux,y,T (a, c, e, g) and its temporal rate of change function utx,t (b, d, f, h) at references T=1000τ≈0.00683013, 2000τ≈0.01366026,70455τ≈0.48121712, and 70485τ≈0.48142203, respectively. For the last pair of figures, we have max0≤x,y≤1ux,y,T≈0.99014166, max0≤x,y≤1utx,y,T≈1.44966641×105. Again, boundary data are omitted for better clarity.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 6
Trajectories of the maximal (blue curves), mean values (red dotted curves) ofux,y,T,0≤x,y≤1, (a); and the maximal (blue curves), mean values (red dotted curves) of utx,y,T,0≤x,y≤1, (b) for T∈0,0.48140837. They are similar to those in Figure 3. Explosive increments of the derivative function values as t approaches the quenching time 70483τ≈0.48140837<T∗≈70485τ≈0.48142203 are again visible. The phenomena agree with (6) and (7) satisfactorily.
(a)(b)Figure 7
Profiles of the numerical solutionux,y,T (a) and utx,y,T (b), at y=0.5,T=0.48141520 before quenching. The profiles are highly symmetric with peaks at the center of the scaled dimensionless combustion domain. Boundary values are not included.
(a)(b)


## 5. Conclusions and Future Endeavours

In this article, we have built and investigated two interconnected second-order finite difference schemes for the numerical solution of the stochastic quenching-combustion partial differential equation (1) equipped with Dirichlet boundary conditions and initial condition (2) and (3). This is for the first time, to the authors’ best knowledge, to use a second-order finite difference scheme for effective approximations of quenching solutions on arbitrary spatial meshes. Numerical results obtained are not only consistent with existing results, but also offering more structure related details. Parameter determination procedures are discussed and analyzed.Computer simulation experiments are focused on two examples with solid combustion-reaction backgrounds [1, 5]. Evenly distributed multiple sparks are deployed in combustion chambers. It is found that both the temperature and its temporal derivative distributions increase monotonically inside the chamber once the ignition starts. The evolution of the temperature field remains smooth and bounded during computations. However, on the other hand, the derivative distribution in the chamber increases rapidly and becomes unbounded as t reaches the neighborhood of the quenching time T∗. This is an indication of a physical combustion reaction of the fuel utilized.The new second-order methods developed on arbitrary spatial grids are extremely straightforward and easy to use in programming implementations. They capture successfully the unique quenching location of the solution of (1)–(3). The numerical results acquired correctly reflect the fact that the final quenching-combustion location is unique and in the center of the camber domain, rather than the multiple initial ignition spark positions suggested by initial functions.Future endeavors in this direction involve continuing numerical analysis and the implementation of numerical algorithms for more rigorous stochastic structure coalitions and conditions. The authors will also pay special attention to the degeneracies, since these are indications of engine metal fatigue [12, 14]. They are extremely important in applications. Higher-order continuous and discontinuous Galerkin methods [15, 16] may also be explored for solving singular quenching problems including (1)–(3). It is also the authors’ intention to inspire more collaborations in this very meaningful and promising territory of advanced theory and numerical methods for the natural world.



---

*Source: 5530744-2021-05-06.xml*