Exact nonclassical symmetry solutions of Arrhenius reaction-diffusion
EExact nonclassical symmetry solutions of Arrheniusreaction-diffusion.
P. Broadbridge , B.H. Bradshaw-Hajek and D. Triadis , .1. (author to whom correspondence should be sent), School of InformationTechnology and Math. Sciences, University of South [email protected]. Dept. of Mathematics and Statistics, La Trobe University, Victoria,Australia.3. Institute of Mathematics for Industry, Kyushu University, Japan. Abstract
Exact solutions for nonlinear Arrhenius reaction-diffusion are con-structed in n dimensions. A single relationship between nonlinear dif-fusivity and the nonlinear reaction term leads to a nonclassical Lie sym-metry whose invariant solutions have a heat flux that is exponential intime (either growth or decay), and satisfying a linear Helmholtz equa-tion in space. This construction extends also to heterogeneous diffusionwherein the nonlinear diffusivity factorises to the product of a function oftemperature and a function of position. Example solutions are given withapplications to heat conduction in conjunction with either exothermic orendothermic reactions, and to soil-water flow in conjunction with waterextraction by a web of plant roots. Keywords:
Arrhenius, reaction-diffusion, plant root extraction, heat conduction,exact solutions, nonclassical symmetries
Assuming that volumetric heat capacity is constant and neglecting reagent con-sumption, the temperature of a reactive mixture satisfies a nonlinear reaction-diffusion equation ∂θ ( r , t ) ∂t = ∇ . [ D ( θ ) ∇ θ ] + R ( θ ) , r ∈ Ω ⊂ R , t ∈ [0 , t ) , (1)where Ω is compact and connected, ∇ is the usual gradient operator, and t > .R ( θ ) is real-valued at any value of absolute temperature θ ∈ [0 , ∞ ). In orderthat heat conduction contributes to increasing entropy, D must be positive (e.g.[1]). One of the most important forms for the reaction function is the Arrheniusreaction term R = (cid:26) R e − B/θ θ > , θ = 0 , (2)with B a positive constant. This follows from the Boltzmann-Gibbs equilibriumcanonical distribution governing the probability of a particle overcoming the ac-tivation energy barrier E of a reaction. The parameter B may then be identifiedwith E/k B where k B is Boltzmann’s constant. The rate factor ρCR , ρ being1 a r X i v : . [ m a t h . A P ] A ug ensity and C being specific heat at constant volume, is the amount of heatenergy released per unit volume per unit time at very high temperature. This isusually described by a weak power-law dependence on temperature, R = S θ m .For example, kinetic theory of a hard-sphere gas predicts m = 0 . m = 0. This will be assumed here, but a straightforwardgeneralisation of the following analysis can cover the case m (cid:54) = 0 with ρC alsodepending on θ .For first-order reactions, E is the energy of bond dissociation into reactivemolecular components such as free radicals. For example, for the exothermicdecomposition of diethyl peroxide, measured reaction rates agree well with theArrhenius law over a broad range of temperatures [3]. In particular, as a com-bustible mixture is controlled by extracting heat through the boundaries of thecontainer, the Arrhenius reaction form is expected to remain appropriate.A full Lie point symmetry classification of (1) was made by Dorodnitsyn et al. [4]. Various combinations of power-laws and exponential functions for D ( θ ) and R ( θ ) lead to an expansion of the Lie group of point symmetry transformations of(1) beyond the common Euclidean isometries in space, and translations in time.These in turn open possibilities of invariant solutions that may be obtained byreduction of variables. The invariant solutions indicate a wide range of possibledynamical behaviours, including stable similarity forms for the temperaturewith multi-peaks, extinguishing, and blow-up in finite or infinite time that canoccur even in solutions with compact support [4].The non-analytic expression exp( − B/θ ) is usually approximated by an ex-ponential function, in the Frank-Kamenetskii approximation [5], before it canlead to useful exact solutions of (1). Although unbounded reaction terms lead toinsight on critical parameters for ignition, there is a significant difference in be-haviour when the reaction term is bounded, as in the standard Arrhenius model.The relationship between the full Arrhenius model and the Frank-Kamenetskiimodel can be more conveniently analysed by using the identity e − E/k B θ = e − E/k B θ exp (cid:18) Θ1 + (cid:15) Θ (cid:19) ; (cid:15) = k B θ E , where Θ = ( θ − θ ) Ek B θ is a rescaled temperature with the zero point shifted to some local value θ (e.g.[6]). Wake and Bazley [7] extended the lower bound e − E/k B θ ≥ e − E/k B θ exp (cid:18) Θ1 + (cid:15) Θ max (cid:19) , beyond the maximum temperature to obtain close approximate critical valuesof parameters at ignition. Gustafson [8] constructed steady linear and radialsolutions of the full nonlinear equation by “a shooting method combined witha Newtonian-Raphson technique and certain boundary value expansions”.2lthough the classical Lie point symmetry classification of (1), with D ( θ )and R ( θ ) arbitrary, selects only constant or unbounded reaction functions, thecomplete nonclassical symmetry classification of reaction-diffusion equations,admits a broader range of possibilities [9, 10, 11, 12], even admitting the Arrhe-nius reaction term [12]. Nonclassical symmetries, in the sense of Bluman andCole [13], leave invariant a system consisting of the original governing equation(1) plus the invariant surface condition that restricts the solution set to only in-variant solutions. The concept of nonclassical symmetry extends more generallyto that of compatibility with an invariance condition, that may be of higher or-der than a point transformation or a contact transformation [14]. Nonclassicalsymmetry analysis reveals that with a suitable nonlinear diffusivity function,after a change of variables, the Arrhenius reaction-diffusion equation admitsseparation of variables, resulting in a linear system. In the following Section 2,time-dependent radial solutions of this system are constructed in two and threedimensions. In the first example, the complicated nonlinear diffusivity function,as well as the temperature, are given exactly. With the same strongly increas-ing diffusivity, the exponentially heating solution is mathematically equivalentto that involving conduction through a finite domain with an endothermic re-action and exponential cooling. The strongly increasing diffusivity resemblessoil-water diffusivity, with the sink term representing plant root absorption. Inthe case of ideal maximal cooling at the boundary, the nonlinear diffusivity isbounded and it is the fixed point of a rapidly converging contraction map. Thesimilarity form of the Kirchhoff variable is given exactly and it is asymptoticallyof the same form as the temperature. The diffusivity varies so little that theexact solution for the Kirchhoff variable closely approximates the temperatureat all times. Particular attention is paid to the case of Newton cooling at theboundary, with non-zero Biot number.In Section 3, it is shown that the same type of solution construction applies toa class of nonlinear reaction-diffusion equations with spatially varying diffusivity.The results of Sections 2 are extended to allow for spatially variable diffusivityand a nonlinear sink term that applies to water transport in unsaturated soilwith extraction by plant roots. A full nonclassical symmetry classification of nonlinear diffusion-reaction equa-tions in two spatial dimensions, was given in [12]. Firstly, (1) is expressed interms of the Kirchhoff variable (e.g. [15]), u = u + (cid:90) θ D (¯ θ ) d ¯ θ . (3)3f u = − (cid:82) θ D (¯ θ ) d ¯ θ for some θ ≥
0, then u = (cid:90) θθ D (¯ θ ) d ¯ θ . In that case, a boundary condition u = 0 corresponds to θ = θ .In terms of the Kirchhoff variable, the reaction-diffusion equation is F ( u ) ∂u∂t = ∇ u + Q ( u ) , (4)where Q ( u ) = R ( θ ) and F ( u ) = 1 /D ( θ ). The starting point of the reductionto the Helmholtz equation is the observation that Equation (4) has a simplenonclassical symmetry ¯ u = e Aε u ; ¯ t = t + ε ; ¯ x j = x j (5)with ε ∈ R , whenever F and Q are related by Q ( u ) = AuF ( u ) + κu, (6)for some A, κ ∈ R . This is not a classical symmetry because it does not ingeneral leave the equation (4) invariant, except when one also assumes theinvariant surface condition u t = Au . Then the symmetry reduction leads to ∇ Φ + κ Φ = 0 with u = e At Φ( x ) . (7)In terms of the original temperature variable θ , the relation (6) is R ( θ ) = (cid:20) κ + AD ( θ ) (cid:21) (cid:34) u + (cid:90) θ D (¯ θ ) d ¯ θ (cid:35) . (8)This gives an explicit construction of R ( θ ) from D ( θ ). Some basic combinations( D ( θ ) , R ( θ )), with R (0) = 0, are given in Table 1. Note that the R ( θ ) functionin case (d) of Table 1 agrees asymptotically with the Arrhenius reaction term R e − B/θ as inverse temperature approaches ∞ .However, in order to construct D ( θ ) from R ( θ ), it is necessary to solve afirst-order nonlinear differential equation. From (6), D ( θ ) = u (cid:48) ( θ ) = AuR ( θ ) − κu . (9)By making u ( θ ) the subject, then differentiating, this implies AD (cid:48) ( θ ) = D κ R ( θ ) + D κ A − R (cid:48) ( θ ) R + DA A − R (cid:48) ( θ ) R . (10)The equation (7) for Φ is simply the linear Helmholtz equation if κ = K > κ = 0. If κ = − K < R >
0, from (9) A must be positive. Then the equation for Φ is the modified Helmholtz equationfor which the radial solutions are unbounded either at the origin r = 0 or atinfinity. 4 ( θ ) R ( θ )(a) θ m ( m > − κm + 1 θ m +1 + Am + 1 θ (b) e θ κ [ e θ − − A [ e − θ − θ ) κ sinh( θ ) − A tanh( θ )(d) R κB (cid:18) Bθ (cid:19) e − B/θ − Aκ R θ ( B + θ ) e − B/θ ( R e − B/θ − AB ) R B ( B + θ ) e − B/θ − AB θ Table 1: Reaction-diffusion combinations that admit nonclassical scaling sym-metry (5). κ = 0 In the case κ = 0, equation (9) is linear homogeneous, allowing direct integrationto obtain u ( θ ) = c A (cid:34) exp (cid:32)(cid:90) θ AR (¯ θ ) d ¯ θ (cid:33) − (cid:35) , (11)with D ( θ ) = c R ( θ ) exp (cid:32)(cid:90) θ AR (¯ θ ) d ¯ θ (cid:33) , (12)where c is an arbitrary constant that is positive (negative) if R is a strictlypositive (negative) source (sink) term. However the spatial dependence of u ,given in (7), is then governed by Laplace’s equation. For example, the isotropicradial solutions u ( x , t ) = u ( r, t ) in two or three dimensions can only be a con-stant added to the singular point-source(sink) solutions with the source(sink)strength varying exponentially in time. The isotherms are specified exactly bythe mapping that follows from Equation(11):( t, θ ) (cid:55)→ u (cid:55)→ Φ = e − At u (cid:55)→ r. In two and three dimensions, the radial solutions take the form r = exp( Φ − c c )and r = c Φ − c respectively.For the case of the Arrhenius reaction (2), the solution (7) is valid when thediffusivity is exactly D = c R exp( B/θ ) exp (cid:18) AR θ exp( B/θ ) − ABR E i ( B/θ ) (cid:19) . where E i is the exponential integral. When θ is large, we can use the powerseries for the exponential integral E i ( x ) = γ + ln( x ) + ∞ (cid:88) n =1 x n nn ! , γ = Euler’s constant, to find D ∼ c R e AB (1 − γ ) /R B − AB/R θ AB/R e Aθ/R . When θ is small, we can use the asymptotic expression E i ( x ) ∼ e x x (cid:34) ∞ (cid:88) n =1 n ! x n (cid:35) , to find that D → θ → . This nonlinear diffusivity is plotted in Figure 1b alongside the Arrheniusreaction function, plotted in Figure 1a.Figure 1: (a) Exact temperature-dependent diffusivity allowing separation ofvariables of Arrhenius reaction-diffusion. (b) Arrhenius reaction rate, as a func-tion of temperature.
Example 1. A heat conduction problem with Arrhenius endothermicreaction term
Even though the enthalpy of reaction is negative, an endothermic reactionwill proceed naturally if the reaction products result in a lowering of the Gibbsfree energy. An endothermic reaction may have a single activation energy, sothat the rate of heat absorption is described by an Arrhenius function. Considera spherical vessel containing the reagents of an endothermic reaction. R isnegative. A small spherical decaying radioactive heat source at the centre,supplies a quantity of heat Q. Then, since u > D > θ > c < A <
0. A solution u = e −| A | t [ c − c r ]can be made to satisfy the boundary conditions u = 0 at r = r , − πr u r = | A | Qe −| A | t at r = r . Note that the latter is simply a heat flux condition that can always be expressedas a linear condition in terms of the Kirchhoff variable u .6 xample 2. Porous media flow with plant-root sink term The combination of exponentially growing nonlinear diffusivity, with boundedsink term, applies to water transport through a web of plant roots [16]. Thedomain of the problem is considered to be a cylindrical pile of soil of radius r = r with a vertical cylindrical injection well in the centre ( r = r < r ). Thedependent variable θ now designates the water content above the plants’ wiltingpoint, where the sink term approaches zero as roots fail to draw water. In thisapplication, R is negative. Then, since u > D > θ >
0, (11)–(12)now imply c < A <
0. A solution u = e −| A | t [ c − c ln r ]can be made to satisfy the boundary conditions u = 0 at r = r , − πru r = | A | Qe −| A | t at r = r . (13)These boundary conditions allow for injection of total water volume Q per unitlength of injection well into a large cylindrical soil mound that is exposed to thesoil-controlled second stage of evaporation (e.g. [17]) at its outer boundary.This solution applies analogously to a vertical current-carrying wire alongthe axis of a cylindrical region, supplying energy for an endothermic reaction. κ > Example 3. Heat conduction with exothermic reaction in a compactregion
With κ = K > u = Φ( r ) exp( At ) satisfies the linear Helmholtz equation.A non-negative isotropic solution can satisfy boundary conditions u r = 0 , r = 0 ,u = 0 , r = r , by choosing Φ = j ( Kr ) in three dimensions, Φ = J ( Kr ) in two dimensionsand Φ = cos( Kr ) in one spatial dimension, with K = λ /r , λ being the firstzero of the spherical Bessel function j ( λ = π ), the standard Bessel function J ( λ ≈ . λ = π/
2) respectively. Further,from the separation of variables in (7), that same solution satisfies other linearhomogeneous boundary conditions such as − u r = Bi u, r = r < r , with Bi constant. This resembles the linear condition for Newton cooling, inwhich Bi is the Biot number, after rescaling and non-dimensionalising vari-ables. The left hand side is heat flux but the right hand side is a multipleof the Kirchhoff variable, rather than temperature. Below, it is demonstratedthat the nonlinear diffusivity has bounded variation and is close to being con-stant. This means that the above boundary condition is very close to that of7ewtonian cooling to a very low-temperature environment. If the physical pa-rameters r , D (0) and B i are prescribed, then the solution parameters are givenby A = − K D (0), where K is the unique solution of − Φ (cid:48) ( Kr ) / Φ = B i /K . K must lie between 0, where − Φ (cid:48) ( Kr ) / Φ = 0 and λ /r , where − Φ (cid:48) ( Kr ) / Φ ap-proaches infinity. With these homogeneous boundary conditions, there alwaysexists an extinguishing solution in which the temperature approaches zero uni-formly everywhere. A demonstration of the stability of this similarity solutionis given in the Appendix A.In the case κ >
0, (9) is equivalent to the canonical form for an Abel equationof the second kind, via w = κu − R ( θ ) , z = − Aθ − R ( θ ) ,w dwdz = w + Φ( z ); Φ( z ) = AR ( θ ) A + R (cid:48) ( θ ) . (14)The standard list of known integrable forms of the Abel equation is given in[18]. In principle, this list may be used to furnish R ( θ ) functions that producesolvable forms of equation (14). If g ( z ) a particular solution for an integrableform of (14), and g − ( w ) is the corresponding inverse function, an R ( θ ) functionthat produces (14) for w ( z ) is given implicitly according to θ = R + g − ( − R ) − A . (15)It is not clear if the segmented solution method of [19] is practicable for thegeneral case.A good approximate analytic reconstruction of D ( θ ) may be obtained byapplying a contraction map towards a fixed point. With u = 0, the solution to(9) satisfying initial value D ( θ ) → , θ → , must be a fixed point of the map D n +1 ( θ ) = M D n ( θ ) = − A (cid:82) θ D n ( s ) dsK (cid:82) θ D n ( s ) ds − R ( θ ) . When considering the separated solution (7), | A | and K may be set to 1 bychoosing dimensionless length and time coordinates Kr and | A | t . D n +1 ( θ ) = ¯ D n ( θ )¯ D n ( θ ) − R ( θ ) /θ , (16)where ¯ D n ( θ ) is the running mean value,¯ D n ( θ ) = 1 θ (cid:90) θ D n ( s ) ds. (17)In the case of the Arrhenius reaction term, we may set R ( θ ) to be R e − /θ by using a dimensionless temperature variable θ/B . Then the maximum value8f R ( θ ) /θ is R /e. From (16), using ¯ D (0) = D (0), it follows that D (0) = 1. From (16) it alsofollows that D ( θ ) > θ >
0. If D ( θ ) is bounded for θ ∈ (0 , ∞ ) and mono-tonic for θ sufficiently large, D ( θ ) must have a limit D ( ∞ ) = ¯ D ( ∞ ). In thatcase it follows from (16) that D ( ∞ ) = 1.Consider M as a mapping on the set S of bounded continuous functions f on (0 , ∞ ) such that f − S = { f ∈ C (0 , ∞ ) : ∃ b ∀ x ∈ (0 , ∞ ) 0 ≤ f ( x ) − ≤ b } . Suppose that D , E ∈ {S} , D n = M n D , E n = M n E , where E n +1 ( θ ) and¯ E n ( θ ) are defined in a analogous way to (16) and (17) above. Now | D n +1 ( θ ) − E n +1 ( θ ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:82) θ D n ( x ) dx (cid:82) θ D n ( x ) dx − R ( θ ) − (cid:82) θ E n ( x ) dx (cid:82) θ E n ( x ) dx − R ( θ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) θ ¯ D n θ ¯ D n − R ( θ ) − θ ¯ E n θ ¯ E n − R ( θ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( ¯ D n − ¯ E n ) R ( θ ) /θ [ ¯ D n − R ( θ ) /θ ][ ¯ E n − R ( θ ) /θ ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . If R < e , then inf { D ( θ ) } = 1 > sup { R ( θ ) /θ } = R /e , consequently (cid:107) E n +1 − D n +1 (cid:107) ∞ (cid:107) E n − D n (cid:107) ∞ ≤ sup { R ( θ ) /θ } [inf { D n ( θ ) } − sup { R ( θ ) /θ } ][inf { E n ( θ ) } − sup { R ( θ ) /θ } ] . ⇒ (cid:107) E n +1 − D n +1 (cid:107) ∞ (cid:107) E n − D n (cid:107) ∞ ≤ e/R − R /e . This is a contraction map, provided R < (3 − √ e ≈ . . For example, consider the case R = 1 for which the above estimate of thecontraction factor is 1 / . R ( θ ) = exp( − /θ ), since D (0) = D ( ∞ ) = 1,it is natural to choose D ( θ ) = 1. Using the standard notation β = 1 /θ , theiterative estimates are D = 1 ,D = 11 − β e − β = ∞ (cid:88) n =0 β n e − nβ (noting βe − β ≤ e − < ,D = II − e − β ; I = (cid:90) ∞ β D β dβ = β − − E i ( − β ) + ∞ (cid:88) n =2 ( n − n n − e − nβ n − (cid:88) k =0 n k β k k !9igure 2: D ( θ ) constructed by series about θ = 0 , . ,
10, as well as approxima-tions D ( θ ) , D ( θ )Each term of the I summation is bounded above by ( n − /n n − , which is oforder e − n / √ n as n → ∞ . The n = 5, e − β terms have a combined value of lessthan 0.01.Note that the sum of the e − β terms has a maximum value of less than 0.01.In Appendix B, D ( θ ) is calculated accurately after obtaining exact series formsfor u ( θ ). The simpler approximate diffusivity functions D j ( θ ) are shown inFigure 2 to agree well with the exact solution.The diffusivity function has a single local maximum where its value is nohigher than 45% above its mean value of 1. From (10), a stationary point of D ( θ ) may occur only where D ( θ ) = 1+ R (cid:48) ( θ ) = 1+ R θ − e − /θ . This expressionhas a maximum value, implying D m ≤ R /e ≈ . . (18)In effect, since at small- t , local disturbances extend by diffusion to a depthproportional to √ Dt and √ D varies by only 20%, the diffusivity is effectivelyconstant for some practical purposes. This is seen in the solution, plotted inFigure 3, wherein the temperature θ ( r, t ) asymptotically approaches the Kirch-hoff variable u ( r, t ) which is identical to the temperature when D = 1.The construction of the nonlinear diffusivity function by a contraction maprelied on the fact that with Arrhenius reaction, βQ ( β ) (= R ( θ ) /θ ) has an upperbound less than 1. The same construction will apply to other reaction laws withthis property, for example Q ( β ) = β m exp( − β ) with m > −
1. This includes thecase m = − / | A | t =-1.5,0,1.5 and 2.5 (solid) approach-ing the Bessel function profile of u(r,t) (x crosses) from below, as time increases. κ < . If κ = − K <
0, from (9) with R ( θ ) > , A must be positive. The diffusivity, D ( θ ), has a singularity where at leading order D (cid:48) ( θ ) ∼ ( κ /A ) D and D ∼ ( θ a − θ ) − / at some freely chosen positive value θ = θ a . In the modellingof unsaturated soil-water flow, the location of the singularity is chosen to beslightly greater than the volumetric water content at saturation [21].Also the equation for Φ is the modified Helmholtz equation for which theradial solutions are unbounded either at the origin r = 0 or at infinity. Thereis an exact solution for exothermic reaction-diffusion on the domain r > r , u = e At K ( Kr ) , (19)where K is the modified Bessel function of the second kind. Although u isunbounded, the corresponding value of θ remains less than the singular value θ a , due to the singularity in D ( θ ) within the integrand in the Kirchhoff trans-formation θ → u . This singular nonlinear diffusivity does not occur in heatconduction but it may occur in models of unsaturated soil-water flow. However,in that application, positive distributed water sources are not common. It ismore common to consider distributed sinks, due to plant roots. Example 4. Porous media flow with plant-root sink term: κ < R <
0, due to plantroots. Because of plant root extraction,
A < r = r . From (8), D > , A < κ < R <
0. However in this case, from the solution(19), θ does not reach zero at any finite value of r but instead approaches zeroat infinite distance. It is to be noted that the modified Helmholtz equation haspreviously been applied to steady state solutions of unsaturated soil-water flow,by analogy to external problems of Helmholtz acoustic scattering [22]. Just asin the latter, we may construct analogous exterior solutions for boundaries ofvarious shapes. However, at large values of r it is well known that the solu-tions agree asymptotically with the isotropic ones, so the spherical or circularboundaries are canonical. We now consider the case in which the medium (or substrate) is spatially hetero-geneous. The appropriate governing reaction-diffusion equation can be written θ t = ∇ . [ f ( x ) D ( θ ) ∇ θ ] + R ( θ ) , (20)where x ∈ Ω ⊂ R n . Here, the heterogeneity is represented by a positive differen-tiable function f that is the amplification factor of the diffusivity. By rewritingthis equation in terms of the Kirchhoff variable (3), we obtain F ( u ) u t = ∇ . [ f ( x ) ∇ u ] + Q ( u ) , (the heterogeneous equation corresponding to (4)) where Q ( u ) = R ( θ ) and F ( u ) = 1 /D ( θ ), as before. Equation (20) admits the same simple nonclassicalsymmetry (5), whenever F ( u ) and Q ( u ) are related by Equation (6).This meansthat R ( θ ) can be constructed from D ( θ ) using (8), or D ( θ ) could be constructedfrom R ( θ ) by solving (9). By again setting u ( x , t ) = e At Φ( x ) we derive thesecond order linear equation for Φ( x ), f ( x ) ∇ Φ( x ) + ∇ f ( x ) . ∇ Φ( x ) + κ Φ( x ) = 0 , (21)which is valid in n dimensions. In the case where A <
0, the solution for θ will be asymptotically similar to u . The relationship between D ( θ ) and R ( θ ) isthe same as that described in the preceeding section, and the solutions for D ( θ )found about also hold in the hetergeneous examples described below. κ = 0 Example 5. Porous media flow with a heterogeneous substrate: κ = 0As an example, we now reconsider the porous medium plant-root extractionproblem described in Section 2.1 above, with the inclusion of spatial heterogene-ity representing a changing scale of soil pores. In Miller self-similar soils, f isthe geometric scale factor of the soil pores ([23, 24]), with f ( r ) = 1 at somechosen reference point. 12n this example, the problem is written in (one-dimensional) radially sym-metric coordinates and the heterogeneity is described by f ( x ) = f ( r ). Thedifferential equation for Φ( r ) replacing Φ( x ) in (21), is fr ddr (cid:18) r d Φ dr (cid:19) + dfdr d Φ dr = 0 , which can be directly integrated for arbitrary f ( r ) to give u ( r, t ) = e At Φ( r ) = ce At (cid:90) rf ( r ) dr with c constant. From (13a), we deduce that A <
0. Having freely chosen f ( r )to be 1, the solution, in terms of the Kirchhoff variable, satisfying boundaryconditions (13), is u ( r, t ) = | A | π e At (cid:90) r r sf ( s ) ds. κ > Example 6. Heat conduction in a heterogeneous medium
Let κ = K . The differential equation for Φ( r ) replacing Φ( x ) in (21) is fr ddr (cid:18) r d Φ dr (cid:19) + dfdr d Φ dr + κ Φ = 0 . (22)With f ( r ) any power-law of the cylindrical radius, exact solutions are readilyavailable (e.g. [18]).With heterogeneity described by the factor f ( r ) = r /r , (22) reduces toAiry’s equation, with general solutionΦ = c A i ( − K / r − / r ) + c B i ( − K / r − / r ) . When f ( r ) is of the form f ( r ) = (cid:18) rr (cid:19) , (23)the resulting ODE is merely a homogeneous Euler equation, with general solu-tion Φ( r ) = c r − √ − ( Kr ) + c r − − √ − ( Kr ) . This solution can be made to satisfy boundary condition Φ( r ) = 0 by setting c = − c r √ − ( Kr ) .For the case with Kr >
1, the general real solution isΦ( r ) = 1 r [ c cos( ω log r ) + c sin( ω log r )] ; ω = (cid:112) ( Kr ) − . For the above cases, solution parameters may be chosen so that u ( r ) = 0.However, r = 0 is a singular point where thermal conductivity is either zero orinfinite. The solution may be regarded as exterior to the surface of a hot wireat r = r from where a finite quantity of heat is supplied.13 .3 The case κ < Example 7. Porous media flow with a heterogeneous substrate: κ < κ = − K <
0, the solution with f ( r ) = ( r/r ) is simplyΦ( r ) = c r − √ Kr ) + c r − − √ Kr ) and with f ( r ) = r /r ,Φ = c A i ( K / r − / r ) + c B i ( K / r − / r ) . The parameters c i and K may be chosen so that the exterior solution satisfies thesame boundary conditions on [ r , ∞ ] as in the case of a homogeneous medium. Provided the nonlinear diffusivity and the nonlinear reaction term satisfy asingle relationship, the Kirchhoff variable u , which is the integral of nonlineardiffusivity, admits solutions that are obtainable by separation of variables to alinear system, whose solution is an exponential in time, multiplying an arbitrarysolution of the Helmholtz, modified Helmholtz or Laplace equation in space. Theheat flux is merely −∇ u , which is given explicitly everywhere in these solutions.If the nonlinear diffusivity function is specified, then the compatible reactionfunction can be constructed directly by integration. If the nonlinear reactionfunction is specified, then the diffusivity function is the solution of a differen-tial equation that is equivalent to an Abel equation if u satisfies a Helmholtzequation, or to a separable equation if u satisfies Laplace’s equation. To thebest of our knowledge, this provides the only known exact closed-form solutionsin two and three dimensions, of nonlinear reaction-diffusion equations with theclassical Arrhenius reaction term. Even when D ( θ ) satisfies an Abel equationthat is not of known solvable type, in some circumstances it may be specifiedto arbitrary precision using exact series expansions.This construction is valid in any natural number n of dimensions, and itgeneralises also to heterogeneous extensions of the Helmholtz factor equation.Applications are given for radial solutions of exothermic reactions with nonlinearheat conduction, endothermic reactions with nonlinear heat conduction, andwater flow from a supply well into a cylindrical soil mound with soil-limitedevaporation at the outer boundary. The logistic shape of the negative Arrheniusreaction term resembles the behaviour of distributed plant roots that have amaximum and minimum value of water uptake rate near saturation and wiltingpoint respectively.As is well known from acoustic scattering theory, exterior solutions of theHelmholtz equation typically asymptotically approach radial symmetric solu-tions at large distances from the scattering surface. Hence, the radial solutionsillustrated here are in a sense, canonical. The solution method used here, in-volves a free function of the Helmholtz equation so it would not be difficult to14se known non-radial solutions. The special solutions presented here could atleast be used as bench tests for more broadly applicable approximate solutionmethods.This approach will lead to ongoing investigations of other nonlinear partialdifferential equations with more than one free function, that may admit specialnonclassical symmetry reductions. Conflicting interests:
The authors have no conflicting interests.
Authors’ Contribution
PB conceived the general approach, provided thesolution in the case K=0, provided recursive approximations for diffusivity anddirected the project. BH provided the solutions for heterogeneous models. DTprovided series constructions for nonlinear diffusivity. All authors contributedto the writing, and all agree on the current version.
Acknowlegements.
We thank Dr Joanna Goard, Prof. Joel Moitshekiand Dr Edoardo Daly for useful discussions on the symmetry reductions andplant-root absorption.
References [1] Broadbridge P. 2008 Entropy Diagnostics for Fourth Order PartialDifferential Equations in Conservation Form.
Entropy , pp. 295–311.(doi:10.3390/e10030365)[2] Glassman I. 1996 Combustion , Chap. 2, 3rd edn. San Diego: Academic Press.[3] Fine D. H., Gray P. & MacKinven R. 1970 Thermal effects accompanyingspontaneous ignitions in gasses. II. The slow exothermic decomposition ofdiethyl peroxide.
Proc. Roy. Soc. A (1525), pp. 241–254.[4] Galaktionov V. A, Dorodnitsyn V. A, Elenin G. G, Kurdyumov S. P. &Samarskii A. A. 1988 A quasilinear heat equation with a source: peaking, lo-calization, symmetry exact solutions, asymptotics, structures.
J Sov Math ,pp. 1222–1292.[5] Frank-Kamenetskii D. A. 1955 Diffusion and Heat Exchange in ChemicalKinetics . Princeton: Princeton University Press.[6] Fowler A. C. 1997
Mathematical Models in the Applied Sciences , Chap. 12.Cambridge: Cambridge University Press.[7] Wake G. C. & Bazley N. W. 1981 Criticality in a model for thermal ignitionin three or more dimensions.
J. Appl. Maths. Phys. (ZAMP) , pp. 594–602.[8] Gustafson K. E. & Eaton B. E. 1982 Exact solutions and ignition parametersin the Arrhenius conduction theory of gaseous thermal explosion. J. Appl.Math. Physics (ZAMP) , pp. 392–405.159] Arrigo D. J., Hill J. M. & Broadbridge P. 1994 Nonclassical Symmetry Re-ductions of the Linear Diffusion Equation with a Nonlinear Source. I.M.A. J.Appl. Math , pp. 1–24.[10] Clarkson P. A., Mansfield E. L. 1994 Symmetry reductions and exact solu-tions of a class of nonlinear heat equations. Physica D , 250–288.[11] Arrigo D. J., Hill J. M. 1995 Nonclassical symmetries for nonlinear diffusionand absorption. Stud. Appl. Math. , pp. 21–39.[12] Goard J. M. & Broadbridge P. Nonclassical Symmetry Analysis of Nonlin-ear Reaction- Diffusion Equations in Two Spatial Dimensions. 1996 Nonlin.Analysis : Theory, Methods Applications , pp. 735–754.[13] Bluman G. W. & Cole J. D. 1969 General similarity solution of the heatequation. J. Math. Mech. , pp. 1025–1042.[14] Olver P. J. 1994 Direct reduction and differential constraints. Proc. Roy.Soc. London A , 509–523.[15] Philip J. R. 1969 Theory of infiltration.
Advances in Hydroscience , 215–296.[16] Feddes R. A., Kowalik P. J. & Zaradny H. 1978 Simulation of field wateruse and crop yield , Simulation Monograph Series. Wageningen: Pudoc.[17] Broadbridge P. & Stewart J. M. 1999 Calculation of humidity during evap-oration from soil.
Adv. Water Resour. , pp 495–505.[18] Polyanin A. D. & Zaitsev V. F. 2003 Handbook of Exact Solutions forOrdinary Differential Equations , Sec. 1.3, 2nd edn. Boca Raton: Chapman &Hall/CRC Press.[19] Panayotounakos D. E. & Zarmpoutis T. I. 2011 Construction of ExactParametric or Closed Form Solutions of Some Unsolvable Classes of Nonlin-ear ODEs (Abel’s Nonlinear ODEs of the First Kind and Relative Degener-ate Equations).
Internat. J. Mathematics and Mathematical Sciences ,doi:10.1155/2011/387429.[20] Mathworks 2013
MATLAB2013b user’s guide .[21] White I. & Broadbridge P. 1988 Constant Rate Rainfall Infiltration: AVersatile Non-linear Model. 2. Applications of Solutions.
Water Resour. Res. , pp.155–162.[22] Waechter R. T. & Philip J. R. 1985 Steady Two- and Three-DimensionalFlows in Unsaturated Soil: The Scattering Analog. Water Resour. Res. ,pp. 1875–1887.[23] Miller E. E. & Miller R. D. 1956 Physical theory for capillary flow phenom-ena. J. Appl. Physics (4), pp. 324–332.1624] Broadbridge, P. 1988 Integrable forms of the one-dimensional flow equationfor unsaturated heterogeneous porous media. J. Math. Phys. , pp. 622–27.[25] Weigao G. & Weber R. O. 1995 Comparison theorems of reaction-diffusionequations and their applications. Nonlinear Analysis: Theory, Methods Ap-plics. (5), pp. 655–665.[26] Ozisik, M. N. 1968 Boundary Value Problems of Heat Conduction . Scran-ton: International Textbook Co.
Appendix A. Stability of similarity solution.
It is shown here that the extinguishing radial similarity solution of the Arrhe-nius reaction-diffusion equation given in Section 2.2, is exponentially stable tosmall perturbations. Change variables to a set of canonical variables of thenonclassical symmetry, namely r, t and v = ue − At . In terms of these variables,the reaction-diffusion equation may be expressed F ( u ) v t = ∇ v + K v, with boundary condition v ( r ) = 0 . The similarity solution is a pseudo-steadystate, v s = Φ( r ) = A J ( λ r/r ), corresponding to exponential decrease of theKirchhoff variable u s = Φ( r ) e −| A | t and satisfying Φ (cid:48) (0) = 0 and Φ( r ) = 0 . Nowconsider a perturbed solution, in plane polar coordinates v = Φ( r ) + w ( r, φ, t ) , with twice-differentiable initial condition v = Φ( r ) + w ( r, φ ); | w | < (cid:15) << . Then order- (cid:15) perturbation w satisfies the same homogeneous boundary condi-tions, plus 2 π -periodicity with respect to φ , plus F ( u s ( r, t )) w t = ∇ w + K w. Unlike the original equation for temperature θ , this is a linear equation with nosquared derivative terms, allowing recourse to comparison theory of reaction-diffusion equations (e.g. [25]). Note that 1 /F ( u ) = D ( θ ) and that from Section2.2, D m ≥ D ≥ D (0) = − A/K . Hence by comparison, | w ( r, t ) | must decay atleast as fast as the solution to the linear initial-boundary problem with smallerdiffusivity and larger reaction term, Q t = D (0) ∇ Q + D m K Q ; Q ( r, φ,
0) = w ( r, φ ); Q ( r , φ, t ) = 0 . Let τ = | A | t/κ ; and q = Qe − [ D m /D (0)] κ τ . Then we have the standard linear heat equation q τ = ∇ q with q satisfyingthe same initial and boundary conditions as Q . The solution q ( r, φ, t ) may17e expanded as a standard Fourier-Bessel series (e.g. [26]). Without loss ofgenerality, we neglect the first component in the series expansion of w , which isa multiple of J ( λ r/r ), and which takes the form of the assumed unperturbedsimilarity solution for v . Each of the other terms in the series for Q is of theform J n ( λ n,m r/r ) [ A n,m cos( nφ ) + B n,m sin( nφ )] e ([ D m /D (0)] λ , − λ m,n ) τ/r , (24)where n ∈ N , λ n,m is the m ’th zero of Bessel function J n , while A n,m and B n,m are arbitrary real coefficients. Note that beyond λ , , the next smallest root is λ , . From the estimate (18), each of the terms (24) is decreasing exponentiallyin time, provided R < e [( λ , /λ , ) − / ≈ . . (25)In terms of the original physical parameters of the unscaled boundary valueproblem, this criterion is R r Bλ , D (0) < . , (26)which is satisfied in practical cases. Appendix B. Evaluating D ( θ ) using asymptotic series. The reaction term R ( θ ) = R exp( − B/θ ) is of special interest. With κ = K ,the parameters K , | A | and B in equation (9) may be set to 1 by adoptingappropriate dimensionless variables. By successive isolation of leading-orderterms, the asymptotic behaviour of u ( θ ) near θ = 0 can be shown to be u ( θ ) ∼ − sign( A ) θ + R θ exp( − /θ ) ∞ (cid:88) r =0 θ − r exp( − r/θ ) ∞ (cid:88) m =0 θ m q r,m . (27)This series structure assumes u = 0 but extension to u > q r, − = 0, substitution into (9) shows that q ,m = ( − m m ! , for m ≥
0; (28) q r, = ( − sign( A ) R ) r r + 1 , for r ≥
0; (29) q r +1 ,m +1 = sign( A ) R r + 2 (cid:26) sign( A )( r − m ) R q r +1 ,m + ( r − j − q r,m − ( r + 1) q r,m +1 + r (cid:88) s =0 m (cid:88) l =0 q r − s,m − l (cid:2) ( s + 1) q s,l − ( s − l ) q s,l − (cid:3)(cid:27) , for r, m ≥ . (30)With the set of coefficients { q r,m } determined iteratively as above, the series(27) is divergent, but can still be used to give accurate specifications of u ( θ ) for18 sufficiently small. Partial sums of (27) may be produced by evaluating both r and m sums up to some maximum integer. Performance of the partial sumsis improved by substituting the known asymptotic form − exp(1 /x ) x Ei( − /x ) ∼ ∞ (cid:88) m =0 m !( − x ) m . (31)A Taylor series expansion for u ( θ ) centered on a non-singular point 0 < θ < ∞ is comparatively easy to derive. We adopt u ( y ) = ∞ (cid:88) n =0 λ n y n , y = θ − θθ ; (32)and by substitution, can show that λ i +1 = 1( i + 1)( λ − R exp( − /u )) (cid:26) sign( A ) θ λ i − i (cid:88) n =1 ( i − n + 1) λ i − n +1 (cid:104) λ n + R θ F (cid:0) n + 1 , − /θ (cid:1)(cid:105)(cid:27) , (33)for i ≥ . Here λ = u ( θ ) is assumed known, and the above sum evaluates to zero for i = 0. The asymptotic series form of u ( θ ) as θ → ∞ can also be ascertainedand appears to have a non-zero radius of convergence.We can use the above series to evaluate u for 0 < θ < θ max to a prescribedaccuracy by matching expansions about different points. To demonstrate, take R = 1 and A <
0. Using (27) with r max = m max = 10, we can ascertain u (1 /
10) = 0 . θ = 1 /
2, with u (1 /
2) =0 . . <θ < .
9, and matches the accuracy of the known value at θ = 1 /
10. Introducinga second expansion about θ = 10, with u (10) = 11 . θ = 19 .
5, and is compatible with the known value u (9 /
10) = 1 . D ( θ ) = du/dθ . To obtain a solution of greater accuracy, we can begin with evaluationat a point θ < /
10 according to series (27). A greater number of linked Taylorseries may then be needed to cover 0 ≤ θ ∗ ≤ ..