Mathematical analysis of a combustible viscoelastic material in a cylindrical channel taking into account induced electric ﬁ eld: A spectral approach

: In spite of the enormous applications of heating combustible materials due to exothermic chemical reactions, scientists and engineers still face a problem with these materials ’ tendency to self-ignition, which can result in signi ﬁ cant property damage if serious precautions are not taken. Therefore, the thermal decomposition of combus-tible viscoelastic material in a cylindrical channel is investigated in this study. With a third-order constitutive model considered for viscoelastic ﬂ ow behavior, momentum and energy-balanced governing equations are provided. The chemical reaction of the material is assumed exothermic and thus follows Arrhenius ’ s kinetics. A numerical solution is provided for the boundary value problem via the bivariate spectral collocation method, and the impact of kinetics parameters on the combustible viscoelastic material is investigated. Our ﬁ ndings show that self-ignition is delayed with higher values of activation energy parameter ( ε ) and third-order parameter ( γ ), as well as lower values of magnetic ﬁ eld parameter ( M ), current density parameter ( δ ), and electrical conductivity exponent ( m ).


Introduction
The relevance of non-Newtonian fluids in technical and industrial processes has prompted an increase in non-Newtonian fluids research during the last few decades.When it comes to non-Newtonian fluid behavior, many constitutive equations have been presented in the scientific literature.Among them, differential constitutive equations have gotten a lot of attention.Differential constitutive models such as the second-order model, the third-order model, the Maxwell model, the Oldroyd-B model, and the Giesekus model are often used to explain polymer behavior [1].The second-and third-order models are subclasses of the Rivlin and Ericksen models.The second-order model, which is the simplest subclass, has received a lot of attention.However, this model cannot account for shear thinning and thickening in a unidirectional steady flow problem.Hence, the third-order model provides a better understanding of viscoelastic materials behavior than the second-order model [2].Extensive work on third-order model can be found [3][4][5][6][7].
Heating of combustible materials as a result of exothermic reaction is found useful in engineering and industries for storage of cellulose, heavy oil recovery, pyrolysis of biomass and coal, incineration of waste product, combustion of solids, etc.Despite the usefulness of the heating of combustible materials, one issue the world faces is the self-ignition of such materials, which may cause significant property damage and countless lives if necessary precautions are not followed.In view of this, several authors have reported in literature the causes of this self-ignition and how they can be controlled.For instance, Zhdanova et al. [8] looked into the impact of compartment fire behavior at ignition development stages on the working condition of fire detector.They discovered that the fuel sample mass had no effect on the flame detector reaction time during steady combustion.In addition, as the bulk of the fuel sample increased, the smoke detector activation time reduced.Lebelo et al. [9] examined spontaneous combustion due to exothermic chemical reaction in a sphere by considering some kinetic parameters controlling the system.The research findings demonstrate that variables like activation energy and heat loss parameters limit exothermic reactions, which lowers oxygen consumption in the system and delays the system's self-ignition.In a combination of air and nitrous oxide, Mitu et al. [10] examined the ignition temperature of a flammable liquid.They came to the conclusion that when the concentration of nitrous oxide increases, so does the ignition temperature.Using a heated plate and dried sewage particles, Polka et al. [11] studied spontaneous combustion and ignition.They also noted that the amount of stored material affects the selfignition temperature value.Korolchenko et al. [12] analyzed ignition of a gas explosion in a premise.They concluded that combustible materials in the corners of the premise are less affected by the flame and combustion products of the gas.Furthermore, combustible material ignites when external air enters after explosion.The thermal decomposition of combustible materials due to exothermic reaction in a sphere, cylinder, and slab with two-step combustible materials are examined in [13][14][15][16][17].
Much attention has been paid to a better understanding of the electrically conducting fluid due to their wide range of applications in electrolyte, polymer, hydroponic, aquaculture, molten metals, and so on.Based on these various applications, lots of work on how to enhance magnetohydrodynamics materials for manufacturing and engineering purposes have been reported in the literature.For instance, Salawu and Fatunmbi [19] numerically explored the second-order fluid with thermal stability and varying electrical conductivity in the Poiseuille device.Makinde and Onyejekwe [20] investigated thoroughly heat transfer of hydromagnetic generalized Couette flow with varying electrical conductivity and viscosity.Entropy minimization in transient hydromagnetic Couette flow with varying electrical conductivity was studied numerically by Chinyoka and Makinde [21].An unsteady hydromagnetic Powel-Eyring fluid with varying electrical conductivity was studied by Salawu et al. [22] who utilized finite deference approach to determine thermal stability and entropy generation of transient reactive magnetohydrodynamic (MHD) Powell-Eyring fluid with variable properties.Thermal stability of second-order and Oldroyd 8-constant MHD fluid flow in a channel taking into account varying electrical conductivity was studied by Salawu et al. [18,23], respectively.
Inspired by the aforementioned literature, and with the assumption that combustible viscoelastic material follows third-order constitutive relation, this study considers thermal decomposition of unsteady third-order fluid with induced electric field in a horizontal cylindrical channel.Bivariate spectral collocation approach (BSCA) is implemented to provide numerical solution for the time-dependent differential equations governing the system.The advantage of this method over other methods such as finite difference method and finite element method is that it requires less discretization to obtain suitable solution.

Problem formulation
We consider transient flow and heat transfer of an incompressible reactive third-order fluid driven by a pressure gradient in a cylindrical channel (Figure 1).The fluid is considered electrically conducting and also assumed flowing in the presence of a magnetic field with induced electric field taken into consideration.Magnetic Reynold number is considered to be minor; hence, the induced magnetic field is negligible.At time = t ¯0, the flow is considered fully developed, and the combustion process begins at > t ¯0.Heat loss to the environment through the wall channel is assumed to follow Newton law of cooling defined as: ), where h t is the heat transfer coefficient, k is the fluid thermal con- ductivity, T is the absolute temperature, and T a is the ambient temperature.

Flow analysis
Following the above assumption, the equation governing third-order fluid in a cylindrical channel is given as follows [24][25][26][27]: Here, V ¯is the velocity, t ¯is the time, τ connotes stress tensor, ρ is the density, and × J B (where B is magnetic field and J , is the current density) is the Lorenz force.J is assumed to be a function of variable electrical conductivity, i.e., where σ T ( ) is the temperature-dependent electrical con- ductivity, E f is the electric field, and = B B 0 is the magnetic field.The Cauchy stress tensor for a thermodynamic compatible third-order fluid is given as follows [25]: where S, the stress tensor, is written as follows: and Here, P stands for pressure, I denotes tensor identity, μ α α , , 1 1 2 , and β 3 are fluid constants, Tr represents trans- pose, A 1 and A 2 are Rivlin-Ericksen tensors, and For an incompressible fully developed flow, continuity and momentum equations become where Continuity Eq. ( 7) implies that = V w r t ¯¯¯, ( ).In view of Eqs ( 11) and ( 12), Eqs ( 8) and ( 9), respectively, yield subjected to the following initial and boundary conditions: Once the flow field is determined from Eq. ( 14), the actual pressure field can be obtained from Eq. ( 13).

Heat transfer analysis
By neglecting reactant consumption, the heat transport equation is given as follows [28,29]: where c p is the specific heat, T stands for the absolute temperature, and k is the thermal conductivity.We also have Q, the heat generated in an exothermic reaction pro- cess, denoted as follows: where Q 0 is the reaction heat, A represents the rate constant, C stands for the initial concentration, E is the activation energy, and R is the universal gas constant.
In view of Eqs ( 3), (12), and ( 17), Eq. ( 16) results into the following: with initial and boundary conditions as follows: We express variable electrical conductivity as follows [18,20]: where σ 0 is the initial electrical conductivity of the fluid and m is the electrical conductivity variation index.The following dimensionless variables are utilized in Eqs ( 14), ( 15), (18), and (19).
The non-dimensional momentum and energy equations are, therefore, obtained as follows: with dimensionless initial and boundary conditions as follows: where A, Pr, V d , α, γ, δ, ε, λ, Bi, and θ a are, respectively, pressure gradient, Prandtl number, viscous heating, second-order material, third-order material, loading electric field, activation energy, Frank-Kamenetski and magnetic flux, Biot number, and ambient temperature.Another quantity of interest considered in this work is a non-dimensional hydromagnetic flow current density C z ( ) defined as follows: 3 Numerical solution Spectral collocation method has been adopted by various researchers in handling ordinary differential equation and partial differential equation problems [30][31][32]34,35] due to its simplicity and rapid convergence, compared to other methods in the literature such as Keller box finite deference method, finite element method, and homotopy analysis method.In this study, spectral collocation method implemented by Uddin et al. [30] is modified into bivariate to solve partial differential Eqs ( 22)-( 24).The method is explained and applied as follows: Consider a function v x τ , ( ) with collocation points for x and τ defined as follows: The derivative of v x τ , ( ) at collocation points x i and τ j is given in matrix form of Chebyshelve differentiation matrices "D" and "d" as follows [30,33]: , with entries of matrices "D" and "d" given as: , and where x τ [ ( ) ( ) ( ) ( ) ( )] and ′ stands for transpose.The general derivative of order(s) ) is represented as follows: In order to apply BSCA on PDEs ( 22)-( 24 With respect to Eqs ( 26)- (30), Eqs ( 22)-( 24) are collocated at the grid points )and = τ j N 0, 1, …, j τ ( ) to obtain the following discretized equations: with discretized initial and boundary conditions as follows: where . Eqs (31)-(33) constitute a system of nonlinear algebraic equations ) with unknown vector functions, w w w 0, 0 , 0, 1 , 1, 0 ).These equations are solved using Newton method by calling FindRoot in mathematical symbolic package MATHEMATICA 11.3.For a given set of parameter values, the reactive third-order material develops over time until a steady state solution is reached.To assess the method's accuracy, we compute absolute residual errors R R , (these parameter values are kept fixed throughout the computation unless otherwise stated in the graphs or tables).Figures 2 and 3 show the plot of respective absolute residual errors R w and R θ defined as follows: 4 Discussion of results In this section, numerical computation of the obtained results are provided.We chose the fluid Prandtl number to be 10, which represents liquid with momentum diffusivity dominant over thermal diffusivity.Other parametric parameter values are picked at random to investigate their  influence on the flow and heat transfer behavior.We also compute the critical values of Frank-Kamenetski λ c ( ) at the steady sate (when the combustion process is time-independent).Critical values are circumstances that should not be surpassed during an exothermic reaction process in order to prevent self-ignition.Figure 4 displays bifurcation curve of = θ θ 0 max ( ) against λ.The graph illustrates that the system is thermally stable when ≤ < λ λ 0 c .Auto-ignition occurs at the upper limit of λ c , and the system becomes thermally unstable.Real solution does not exist when > λ λ c .
The findings in Table 1 indicate the impacts of various parameters on thermal criticality levels.The thermal criticality magnitude increases with increasing values of ε γ , and reduces with increasing values of δ m , , and M , as shown in the table.This clearly illustrates that self-ignition is delayed with an increase in activation energy and viscoelastic characteristics but hastened with an increase in electric field loading, electrical conductivity, and magnetic field intensity.

Transient velocity and temperature profiles
Figures 5-8 depict the time development of the velocity and temperature components in the channel.As time passes, the fluid's velocity and temperature rise until they reach a steady state maximum value, as seen in Figures 5 and 6.Furthermore, velocity profiles achieve a steady state sooner than temperature profiles.This is to be expected given that velocity acts as a source of heat for temperature.The flow transports thermal energy, resulting in an improved thermal profile.Likewise, the energy Eq. ( 23) clearly shows that the fluid velocity and velocity gradient, in terms of joule heating and viscous dissipation, give an extra heat source to the system.

Parameter dependence of velocity and temperature profiles
Figures 9 and 10 reveal the impact of the magnetic field (M ) on the magnitude of velocity and temperature, respectively.
An increase in M correlates to an enhancement in the Lorentz force resisting flow motion and thereby decreases velocity profiles (Figure 9).Increasing values of M improve the heat source strength in the energy equation and thus increases the temperature profiles (Figure 10).Figures 11  and 12 depict an intriguing phenomena when the fluid temperature rises in proportion to the increase in V d and λ.This may be understood by the fact that when V d increases, the influence of both viscous and joule heating increases, causing the temperature profiles to rise.Likewise, when λ improves, so does the internal heat produced inside the flow system as a result of the exothermic reaction, resulting in a rise in fluid temperature.Mathematical analysis of a combustible viscoelastic material  7 The behavior of velocity and temperature distributions with increasing values of third-order parameter γ ( ) is seen in Figures 13 and 14, respectively.The fluid's viscoelastic characteristics rise with an increase in γ.This makes the material to be more thicker, causing an increase in flow resistance and, as a consequence, a decrease in fluid velocity (Figure 13).Reduced velocity minimizes the magnitude of viscous heat source in the temperature equation.This results in a drop in temperature profiles (Figure 14).
Figure 15 depicts the influence of electrical conductivity exponent (m) on velocity and temperature distributions.As m increases, the current density decreases.Hence, the Lorentz force known as retarding force diminishes.The net outcome is that the flow intensifies as m increases (Figure 15).Furthermore, when m grows, the term θ m reduces, and the intensity of the heat source component in the temperature equation declines, leading to lower temperature profiles (Figure 16).
The circuit design may also impact the behavior of an MHD flow.A short circuit configuration occurs when the external load is zero = δ 0 ( ); else, the circuit is classified as open.A short circuit design = δ 0 ( ) reduces the current density magnitude, and there is less resistance in the flow direction and more flow in the channel.This accounts for higher velocity profile when = δ 0 in Figure 17.In   addition, viscous and joule heating dissipation are less in magnitude in a short circuit condition, and hence less temperature loss is observed (Figure 18).Higher values of Bi mean more heat is lost to the environ- ment by convection, leading to temperature reduction.

Current density profiles
The impact of δ m γ , , , and M on current density profiles is depicted in Figures 21-24.Current density improves as δ grow but declines as γ, M and m increase.
This shows clearly that higher values of electric field loading enhance current density, while higher values of   Mathematical analysis of a combustible viscoelastic material  9 electrical conductivity, magnetic field intensity, and non-Newtonian parameter reduce current density.

Conclusion
The numerical investigation of an unsteady hydromagnetic reactive third-order material with variable electrical conductivity was carried out using BSCA.The impacts of different thermophysical factors on velocity, temperature, current density and conditions for self-ignition are correctly determined.Significant results are summarized as follows: • The velocity distributions reach steady state faster than the temperature distributions.• The material temperature rises when parameter values of M , λ, δ, and V d increase but drop as parameter values of m and γ grow.• The viscoelastic material velocity increases with increasing values of m but decreases as the values of M , δ, and γ increase.• The magnitude of current density increase with higher values of δ but reduces with an increase in the values of M , m, and γ • Self-ignition is delayed with increase in ε and γ values but hastened with higher values of M , δ, and m.

Figure 11 :
Figure 11: V d variation for temperature distributions.

Figure 19
reveals the impact of ε on temperature profiles.Temperature profiles reduces with increasing ε values.This is attributed to the fact that increasing ε's values lead to a reduction in heat source term the energy equation and thus reduced the temperature profiles.The influence of Bi on temperature distributions is presented in Figure20.Temperature profile seems depreciated with increasing values of Bi.

Figure 16 :
Figure 16: m variation for temperature distributions.