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.
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 an -dimensional heat engine, . Let be the temperature distribution inside of its combustion chamber and be the initial temperature distribution due to the sparks. Assume that combustion occurs at the unit temperature. Further, let , 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:where , is its boundary, is the -dimensional gradient vector, the coefficient is positive, , and is the stochastic indicator for [5, 6].
The nonlinear source function is strictly increasing with respect to and
A typical example of such reaction term iswhere is the internal combustion index utilized [3, 5–7].
The solution of (1)–(4) is said to quench if there exists a finite such that
Such a value 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
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 and , there exists a critical value 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 , a numerical method must be capable of evaluating simultaneously correct critical value , quenching time , quenching location , and extremely sensitive and possibly unbounded derivative function [5, 6, 11]. This can be difficult since values of remain bounded and well-behaved until time reaches a certain neighborhood of , if it exists. Furthermore, we may also observe from properties (6) and (7) that, while grows exponentially, or faster than exponentially, in the aforementioned neighborhoods, the solution 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 , or by using fixed mesh steps (cf. [4, 7, 8, 12–14] and references therein). In these procedures, critical values, quenching times, solution , and the rate-of-change function 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 to , where is the number of dimensions and 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 and . In the circumstance, problems (1)–(3) can be simplified to
We consider an arbitrary partition of the spatial domain , that is, gridsfor which . Denote . Apparently, such sets are -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 . At any spatial location , we adopt usual notations or without confusion.
2. Semi-Discretization Approximations
At the time level , we setfor any available index . Without loss of generality, we drop the superindex for simplicity. We have expansions of the spatial derivatives
Strategy 1. From the above, we must haveNote that due to the homogeneous condition (9). Substituting the above into (8)–(10), we acquire thatwhereThe band matrix is large in size since 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):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 matrix A 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 , we also have values 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 approximationwhereDropping the truncation error term from the differential system, we obtain another second-order semi-discretization scheme for solving (8)–(10):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 values can be computed through some direct expansions. To this end, we have
We observe readily that
For the simplicity in calculations, we may set . Thus,
To use the above two formulations for a quadratic order approximation, we need to balance coefficients related to and at internal spatial mesh points. Recall (13) and (15) and set
Recall (25). To find the coefficients , we let
It follows therefore
It turns out that
From the last two equations,
Hence,
Consequently,
Similarly, to find coefficients , we consider
Balancing the coefficients, we have
It follows subsequently that
From the last two equations,
Hence,and furthermore,
Now, recall the expansion (25). To determine coefficients in (29), we have
From the above, we obtain the following linear system:
The system can be conveniently solved digitally via any existing solution package. By the same token, for evaluating coefficients , in (31), we find that
The above expansion leads to the following non-homogeneous linear system:
Now, recall (14). To determine coefficients , we have
Therefore, to guarantee overall second-order approximations of the required function values for , we must ask thatwhich can be compressed into a matrix formwhere
Furthermore, for (13), we observe that
Therefore, the result is as follows:
It can be rewritten aswhere
Now, for (15), we have
This gives us the linear systemwhere
Based on the above investigations, we acquire the following result.
Theorem 1. If matrices , and are nonsingular based on the selection of , then the coefficient in (23) is uniquely determined.
We note that properties, such as the spectrums, of matrices and still remain to be studied based on arbitrary sets of at each temporal level of . Bear in mind that nonuniform mesh steps in are determined through particular adaptive procedure employed, such as the arc-length monitoring functions investigated in [7, 13, 17, 19].
4. Simulation Experiments
Let and . 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 andWe further set 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 be , while the temporal discretization parameter is fixed, where is the Courant-Friedrichs-Lewy constant, and 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 distribution 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, [7, 21]. The evolution of accelerates rapidly as time is approaching the quenching time , which is approximately . To see more details, we show the rate of change function in four temporal levels: , and 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 of and , we plot evolutional trajectories of the maximal values and mean values of them for in Figure 3. It can be seen that the maximum of the temperature field derivative increases exponentially as , 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 -semi-logarithmic scale is used to show greater details of the phenomenon, especially for the derivative function .
In Figure 4, three-dimensional surfaces of and , together with their - 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 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, , 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(a)
(b)
(c)
(d)
Example 2. As a highly realistic extension, we consider another favorable natural quenching case with and a much simplified four ignition spark initial condition [2, 3]:We continue setting . 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 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).
Figure 5 is devoted to four key moments of the solution procedures. Distributions of the temperature and its temporal derivative inside the rectangular combustion chamber are simulated. The evidence of initial multiple sparks is obvious at . However, the influence apparently becomes fuzzy at and fuzzier at , while the maximal temperature, , 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 , for which . This is further supported by the corresponding set of 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 at , values of , which is the singular point of the source function , in almost the entire combustion domain. On the other hand, we also observe that, while the distribution profile of is similar to that of , the amplitude of rate function jumps to approximately . 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 solution and its temporal derivative , 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 . 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 Figure 7, profiles of the solution and its temporal derivative are given at immediately before quenching. It can be seen that while the temperature inside the chamber increases uniformly towards the unity, the rate-of-change function , in other words, the velocity of the temperature field evolution, grows rapidly to the infinite at the quenching location .
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
(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 reaches the neighborhood of the quenching time . 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.
Data Availability
The data are available from the corresponding author upon any reasonable request.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
The first, second, and last authors acknowledge the constant support from Baylor University during the realization of this work. The third author acknowledges the financial support from the National Council for Science and Technology of Mexico (CONACYT) through grant A1-S-45928.