Numerical Investigation of a Combustible Polymer in a Rectangular Stockpile: A Spectral Approach

: Despite the wide application of combustion in reactive materials, one of the challenges faced globally is the auto-ignition of such materials, resulting in ﬁre and explosion hazards. To avoid this unfortunate occurrence, a mathematical model describing the thermal decomposition of combustible polymer material in a rectangular stockpile is formulated. A nonlinear momentum equation is provided with the assumption that the combustible polymer follows a Carreau constitutive relation. The chemical reaction of the polymer material is assumed to be exothermic; therefore, Arrhenius’s kinetic theory is considered in the energy balance equation. The bivariate spectral local linearization scheme (BSLLS) is utilized to provide a numerical solution for the dimensionless equations governing the problem. The obtained results are validated by the collocation weighted residual method (CWRM), and a good agreement is achieved. Dimensionless velocity, temperature, and thermal stability results are presented and explained comprehensively with suitable applications. Some of the obtained results show that thermal criticality increases with increasing power law index ( n ) and radiation ( Ra ) values and decreases with increasing variable viscosity ( β 1 ) and material parameter ( W e ) values.


Introduction
In engineering and industries, burning combustible materials in a slab is important for storing cellulosic materials, solid combustion, refuse cremation, heavy oil recovery, and other processes [1][2][3]. Polymers can undergo combustion and release a large amount of energy, which can then be used for transportation, generating electric power, and providing heat for various domestic and industrial applications. Compared with other energy sources, like solar cells, wind generators, and turbines, polymers are quite inexpensive. However, there are some challenges with polymer burning, such as fire ignitions due to human negligence and the physical characteristics of hydrocarbon polymers. These have caused extensive property destruction and claimed the lives of an unknown number of people. From this perspective, several authors have been prompted to study the causes of fire ignition in the combustion process and how it can be controlled.
Drysdale [4] described ignition as the process of initiating a fast exothermic reaction, which then propagates and causes changes in the materials involved and as well generates temperatures in excess. Also, ref. [4] distinguished between two kinds of ignition: (1) piloted ignition, in which a flammable mixture is ignited by a pilot such as an electrical spark or an autonomous flame, and (2) auto-ignition, in which flame develops spontaneously within the mixture. Shi and Chew [5] studied polymers' responses to fire under auto-ignition conditions within a cone-shaped calorimeter and concluded that the CO and CO 2 production process for polymers is a two-step reaction. In the first-step reaction, CO Mathematics 2023, 11,3510 2 of 17 and other flammable substances are generated following Arrhenius's law. The second stage is the oxidation reaction of CO in the presence of air. The whole process is expressed as follows: Step1 : Polymer + O 2 → H 2 O + CO + others product Step2 : CO + 0.5O 2 → CO 2 .
It was also discovered in [5] that the CO and CO 2 emissions from flaming combustion are greater than those from non-flaming combustion.
Moreover, the combustion of polyethylene and polyvinyl chloride was examined in [6][7][8]. These investigations used experimental settings, allowing researchers to acquire more information and factors for pyrolysis, ignition, and combustion. Geschwindner et al. [9] incorporated a mix of high-speed planar laser-induced fluorescence of the HO radical (OH-PLIF) and a thermal decomposition analysis to examine the combustion of micrometer-sized polypropylene (PP) particles. They found that the highest density of flame-retardant polymer particles decreases during ignition and the early stage of burning. Lohrer et al. [10,11] investigated the effects of physical factors such as material wetness, atmospheric humidity, and concentrations of oxygen on the combustion of reactive materials and discovered that water in the reactive material increases auto-ignition.
Some mathematical models that are less expensive and faster than experimental approaches have been implemented in the literature to explain the auto-ignition of combustible materials in a stockpile. For instance, a one-step combustion process of heat transfer in a spherical channel was investigated in [12,13]. According to their reports, the system maintains stability as heat escapes into the atmosphere. In addition, enhancement in the chemical reaction rate leads to increased heat generation in the stockpile, resulting in quick auto-ignition. Lebelo et al. [14][15][16] examined the two-step thermal decomposition of combustible materials in a sphere. They identified that elevation in the two-step kinetic parameter diminishes the rate of heat loss on the sphere's surface, which, in turn, speeds up auto-ignition.
In the literature above reviewed, the authors did not consider flow behavior in their studies. However, an increase or decrease in flow speed contributes to the heat transfer performance of combustible materials. It was discovered in the literature that the Carreau fluid model well describes the flow behavior of polymeric solutions because of its shear rate properties [17]. Several studies on Carreau fluid constitutive relations have been documented. For instance, Siska et al. [18] examined the terminal velocity of non-spherical particles falling through a Carreau fluid and concluded that the Carreau fluid model can well characterize the rheology of various polymeric solutions, such as 1% methylcellulose tylose in glycerol solutions and 3% hydroxyethyl-cellulose Natrosol HHX in glycerol solutions. Ref. [19] reported an intriguing study on the entropy production of Carreau fluid in the presence of infinite shear rate viscosity. Also, the behavior of a Carreau fluid flow past a stretching sheet was extensively analyzed in [20,21]. In addition, the peristaltic movement of a Carreau fluid was also extensively studied in [22][23][24]. For more on the Carreau fluid model with different configurations, see [25][26][27].
Motivated by the reviewed literature in [12][13][14]22,28], this study focuses on an investigation of the thermal decomposition of Carreau fluid in a rectangular stockpile with variable thermophysical properties. It is believed that this present study has not been reported in the literature. However, the outcome of this study could be useful for engineers dealing with the combustion of polymers by determining the conditions necessary for explosions and how to control them. The rest of the article is structured as follows: a mathematical model for the unsteady, fully developed flow and temperature of the polymer is presented in Section 2; Section 3 deals with the application of BSLLM to the dimensionless initial-boundary value problem; in Section 4, an extensive discussion of the obtained findings is provided; and concluding remarks are provided in Section 5.

Mathematical Analysis
A transient laminar flow of a reactive incompressible Carreau fluid material in a combustible stockpile positioned at a distance of 2h apart is considered (see Figure 1). x − axis is parallel to the flow direction, andŷ − axis traversed to it. Initially, the fluid is assumed to be fully developed in the stockpile of the temperature, T 0 , and the material's viscosity and thermal conductivity, denoted as µ = µ 0 and κ = κ 0 , respectively, are assumed to be constant. At timet > 0, the combustion process begins, and the material properties become temperature-dependent: µ = µ(T) and κ = κ(T). The chemistry involved in this problem assumes two-step Arrhenius kinetics. The reactant consumption is assumed negligible in this model. We also assume that the means of heat loss into the environment of the temperature, T a , is mainly via radiation and convection. The influence of density variation with temperature is approximated following Boussinesque approximation. The equations governing the momentum and energy balance under the assumptions above are as follows [14][15][16]20,26]: with the initial-boundary conditionŝ Mathematics 2023, 11, x FOR PEER REVIEW 10 of 21  u¯axial velocity,P¯modified pressure, T¯absolute temperature, T a¯a mbient temperature, T 0¯i nitial temperature of the stockpile, µ 0¯m aterial's dynamic viscosity at temperature T 0 , κ 0¯m aterial's thermal conductivity at T 0 , c p¯s pecific heat at constant pressure, ρ¯material's density, Q 1 , Q 2¯fi rst-and second-step heat of reaction, A 1 , A 2¯fi rst-and secondstep rate constants, C 1 , C 2¯fi rst-and second-step reactant's concentration, E 1 , E 2¯fi rst and second step activation energies, ε¯stockpile's emissivity (0 < ε < 1), σ¯Stefan-Boltzmann constant, β¯volumetric coefficient, g¯gravitational acceleration, h t¯c oefficient of heat transfer, Γ¯time constant, n¯dimensionless power law index. n < 1 represents shear-thinning fluids, n = 0 represents Newtonian fluids, and n > 1 represents shear-thickening fluids The variable viscosity and thermal conductivity, µ(T) and κ(T), are expressed, respectively, as where b 1 and b 2 are dynamic viscosity and thermal conductivity variation parameters. We then introduce the below dimensionless parameters to Equations (1)-(3): The following dimensionless equations are then obtained: where β 1 , β 2 , W e , Pr, λ, 1 , 2 , A, Gr, V d , Ra, ω, Bi, θ a are, respectively, the variable viscosity parameter, the variable thermal conductivity parameter, the material parameter, the Prandtl number, the Frank-Kamenetskii parameter, the activation energy parameter, the activation energy ratio parameter, the pressure gradient, the Buoyancy parameter, the viscous heating parameter, the radiation parameter, the two-step exothermic reaction parameter, the Biot number, and the ambient temperature parameter. We also considered the Nusselt number, defined as follows:

Solution Method
In this section, the BSLLS is implemented to provide a numerical solution for Equations (6)-(8), as outlined in [29,30]. For further studies on the convergence analysis of the BSLLS, see [31]. To adopt the BSLLS, the nonlinear equations, (6) and (7), and the nonlinear convective boundary conditions (8) are, respectively, represented by F, Θ, and Bc.
The iteration technique (quasi-linearization method) is applied independently to Equations (9)- (11) to arrive at where The next step is to transform the physical domains, t ∈ [0, T] and y ∈ [0, 1], respectively, into domains τ ∈ [−1, 1] and x ∈ [−1, 1] using transformations t = T(τ+1) 2 and y = (x+1) 2 , with collocation points It is assumed that solutions to u(τ, x) and θ(τ, x), in the form of bivariate Lagrange's interpolating polynomials, are defined as: where function L p (x) represents the Lagrange cardinal polynomial of the Chebyshev-Gaus-Lobatto grid points, The L q (τ) function is defined similarly. The derivative values at the Chebyshev-Gaus-Lobatto points (x i , τ j ) are computed as follows: where r is the order of the derivative, and D r = 2 r D r i,p andd j,q = 2 T d j,q are Chebyshev differentiation matrices (N x + 1) × (N x + 1) and (N τ + 1) × (N τ + 1) respectively. u j and θ j are defined as Superscript T denotes a transpose. Substituting Equations (17) and (18) into (12) and (13) yields Mathematics 2023, 11, 3510 and I is an identity matrix. Applying spectral collocation to the boundary conditions (8) and the convective boundary condition (14), we have Imposing boundary conditions on Equation (19) for j = 0, 1, . . . , N τ − 1, we obtain the following N τ (N x + 1) × N τ (N x + 1) system of matrices: with The vectors u r+1,N τ and θ r+1,N τ correspond to the initial condition given in Equation (8). The matrices (20) are solved iteratively until suitable results are obtained.

Convergence Analysis
The convergence of the BSLLS is evaluated by considering the error norms between two successive iterations. Error norms are defined as E u and E θ decrease swiftly as the number of iterations increases (see Figure 2a. This shows that the BSLLS converges within a few iterations. Also, residual error norms are computed to show the accuracy of the BSLLS. Residual error norms are expressed as where ∂ u and ∂ θ are corresponding nonlinear partial differential equations: Equations (9) and (10), respectively. Figure 2b depicts residual errors R u and R θ against the number of iterations. Residual errors are found to decrease rapidly with an increasing number of iterations.

Results and Discussion
In this section, we employ the parameter values, as default values, unless otherwise stated in graphs and tables. The results obtained by using the BSLLS are validated with the ones obtained using the collocation weighted residual method (see Table 1), and a good agreement is observed.

Results and Discussion
In this section, we employ the parameter values, n = 0.5, W e = 0.5, A = 1, β 1 = 0.1, β 2 = 0.1, λ = 0.1, 1 = 0.1, 2 = 0.1, ω = 1, Ra = 0.1, Pr = 10, Gr = 1, V d = 0.5, θ a = 0.1, A = 1, T = 120, as default values, unless otherwise stated in graphs and tables. The results obtained by using the BSLLS are validated with the ones obtained using the collocation weighted residual method (see Table 1), and a good agreement is observed.   (Figure 3b) profiles rise until they reach steady-state maximum values. Furthermore, the velocity profile reaches a steady state faster than the temperature profile. This is expected since velocity acts as a source of heat for the combustion process and, hence, increases the temperature profile. Temperature   Figures 3 and 4 show the time development of the velocity and temperature profiles. As time passes, the velocity ( Figure 3a) and temperature (Figure 3b) profiles rise until they reach steady-state maximum values. Furthermore, the velocity profile reaches a steady state faster than the temperature profile. This is expected since velocity acts as a source of heat for the combustion process and, hence, increases the temperature profile.  Figures 3 and 4 show the time development of the velocity and temperature profiles. As time passes, the velocity ( Figure 3a) and temperature (Figure 3b) profiles rise until they reach steady-state maximum values. Furthermore, the velocity profile reaches a steady state faster than the temperature profile. This is expected since velocity acts as a source of heat for the combustion process and, hence, increases the temperature profile.

Solution Blow-Up Profile
The bifurcation plot of the maximum temperature, (θ(0)), versus the Frank-Kamenetskii parameter, λ, is provided to examine the thermal criticality condition of the system. The critical value, λ c , is computed at a steady state (when combustion is independent of time) to explain auto-ignition during the combustion process. Figure 5 displays a bifurcation diagram explaining the thermal behavior of the system as λ increases. The bifurcation diagram illustrates that the solution to Equation (7), at a steady state, is finite for λ between interval 0 and λ c . Auto-ignition is observed at the upper limit of the interval, λ c , and a real solution does not exist when λ > λ c .
To prevent or control spontaneous ignition, the impact of the thermophysical parameters on the thermal criticality is examined (see Table 2). Thermal criticality increases with increasing values of the power law index (n), radiation (Ra), and the Biot number (Bi) and decreases with increasing values of variable viscosity (β 1 ), the variable thermal conductivity parameter (β 2 ), and the material parameter (W e ). This implies that these factors are important in minimizing auto-ignition. In other words, for thermal stability to be maintained during combustion processes, these parameters should be made significant. ent of time) to explain auto-ignition during the combustion process. Figure 5 displays a bifurcation diagram explaining the thermal behavior of the system as λ increases. The bifurcation diagram illustrates that the solution to Equation (7), at a steady state, is finite for λ between interval 0 and c λ . Auto-ignition is observed at the upper limit of the interval, c λ , and a real solution does not exist when c λ λ > . To prevent or control spontaneous ignition, the impact of the thermophysical parameters on the thermal criticality is examined (see Table 2). Thermal criticality increases with increasing values of the power law index ( ) n , radiation ( ) Ra , and the Biot number ( ) Bi and decreases with increasing values of variable viscosity ( 1 β ), the variable thermal conductivity parameter ( 2 β ), and the material parameter ( e W ). This implies that these factors are important in minimizing auto-ignition. In other words, for thermal stability to be maintained during combustion processes, these parameters should be made significant.   Table 2.

Dependence of Velocity and Temperature Profiles on Flow Parameters
Figures 6 and 7 depict the impacts of W e on the velocity and temperature profiles, respectively. Both profiles are elevated as the W e values increase. This is ascribed to the fact that an increase in W e makes the fluid thinner, and the resistance force to the flow decreases; hence, the velocity profile increases. Furthermore, internal heat generation, as a result of the viscous term, is high when W e increases. This leads to an enhancement of the temperature profile. Figures 8 and 9 present the effect of the power law index, n, on the velocity and temperature distributions, respectively. As such, the velocity distribution decreases with increasing n values. When n > 1, the resistance force to the flow for shear-thickening fluids becomes maximal, and the velocity profile decreases (see Figure 8). A decrease in fluid speed results in a decrease in viscous heating, and the temperature profile reduces significantly (see Figure 9). decreases; hence, the velocity profile increases. Furthermore, internal heat generation, as a result of the viscous term, is high when e W increases. This leads to an enhancement of the temperature profile. Figures 8 and 9 present the effect of the power law index, n , on the velocity and temperature distributions, respectively. As such, the velocity distribution decreases with increasing n values. When 1 n > , the resistance force to the flow for shear-thickening fluids becomes maximal, and the velocity profile decreases (see Figure  8). A decrease in fluid speed results in a decrease in viscous heating, and the temperature profile reduces significantly (see Figure 9).        β , as observed in Figure 10. The reason for this is that fluid viscosity reduces as 1 β increases, and the kinetic energy of the fluid molecules increases, leading to velocity distribution enhancement. An increase in the velocity profile naturally enhances the heat source term in the energy equation, resulting in an elevation in the temperature field (see Figure 11). The impact of variable thermal conductivity, 2 β , on the velocity and temperature distributions is provided in Figures 12 and 13, respectively. An increase in the values of 2 β increases the temperature profile (see Figure 13). This is attributed to the fact that a rise in 2 β leads to a reduction in the thermal conductivity term, 2 e β θ − . This leads to the slow, random movement of fluid molecules and thus facilitates heat transfer through the fluid, which consequently enhances the temperature profile. A significant rise in temperature as 2 β increases results in a reduction in the fluid viscosity and thus enhances the flow speed (see Figure 12). Figure 14 depicts the impact of the radiation parameter, Ra , on the fluid temperature. It can be observed that the temperature profile reduces as the Ra values rise. This indicates that more heat exits the stockpile through radiation and, thus, reduces the fluid temperature profile. The same scenario is seen in Figure 15 as the temperature profile reduces with an increasing Biot number, Bi , because more heat escapes the stockpile through the walls. Enhancement in the temperature profile is observed in Figure 16 as the two-step parameter, ω , values increase. This is due to extra heating created by higher values of ω . Figures 17 and 18 illustrate the variation in the Nusselt number, Nu (heat transfer rate at the wall) for different combinations of thermophysical parameters. Nu increases with higher values of the Frank-Kamenetskii parameter, λ , and the two-step reaction parameter, ω , and decreases with higher values of the radiation parameter, Ra , and the Biot number, Bi . These results agree with the work in [32].   Figure 10. The reason for this is that fluid viscosity reduces as β 1 increases, and the kinetic energy of the fluid molecules increases, leading to velocity distribution enhancement. An increase in the velocity profile naturally enhances the heat source term in the energy equation, resulting in an elevation in the temperature field (see Figure 11). The impact of variable thermal conductivity, β 2 , on the velocity and temperature distributions is provided in Figures 12 and 13, respectively. An increase in the values of β 2 increases the temperature profile (see Figure 13). This is attributed to the fact that a rise in β 2 leads to a reduction in the thermal conductivity term, e −β 2 θ . This leads to the slow, random movement of fluid molecules and thus facilitates heat transfer through the fluid, which consequently enhances the temperature profile. A significant rise in temperature as β 2 increases results in a reduction in the fluid viscosity and thus enhances the flow speed (see Figure 12). Figure 14 depicts the impact of the radiation parameter, Ra, on the fluid temperature. It can be observed that the temperature profile reduces as the Ra values rise. This indicates that more heat exits the stockpile through radiation and, thus, reduces the fluid temperature profile. The same scenario is seen in Figure 15 as the temperature profile reduces with an increasing Biot number, Bi, because more heat escapes the stockpile through the walls. Enhancement in the temperature profile is observed in Figure 16 as the two-step parameter, ω, values increase. This is due to extra heating created by higher values of ω. Figures 17 and 18 illustrate the variation in the Nusselt number, Nu (heat transfer rate at the wall) for different combinations of thermophysical parameters. Nu increases with higher values of the Frank-Kamenetskii parameter, λ, and the two-step reaction parameter, ω, and decreases with higher values of the radiation parameter, Ra, and the Biot number, Bi. These results agree with the work in [32].

Concluding Remarks
This article considered the combustion of polymer material in a rectangular stockpile. The rheology of the polymer was assumed to follow a Carreau fluid constitutive relation. The spectral local linearization method was implemented to provide a numerical solution to the problem. The impacts of various thermokinetics factors on the flow and thermal behaviors were examined. From the obtained results, it was found that some parameters (

Concluding Remarks
This article considered the combustion of polymer material in a rectangular stockpile. The rheology of the polymer was assumed to follow a Carreau fluid constitutive relation. The spectral local linearization method was implemented to provide a numerical solution to the problem. The impacts of various thermokinetics factors on the flow and thermal behaviors were examined. From the obtained results, it was found that some parameters (W e , β 2 , and ω) improve the combustion process since the temperature profiles increase as the values of these parameters increase. This may speed up thermal ignition and, thus, lead to an explosion. The opposite case can be observed for an increase in the values of n, Ra, and Bi. This slows down the chemical reaction and, thus, minimizes the combustion process.