A non-field analytical method for heat transfer problems through a moving boundary

This paper presents an extension of the non-field analytical method—known as the method of Kulish—to solving heat transfer problems in domains with a moving boundary. This is an important type of problems with various applications in different areas of science. Among these are heat transfer due to chemical reactions, ignition and explosions, combustion, and many others. The general form of the non-field solution has been obtained for the case of an arbitrarily moving boundary. After that some particular cases of the solution are considered. Among them are such cases as the boundary speed changing linearly, parabolically, exponentially, and polynomially. Whenever possible, the solutions thus obtained have been compared with known solutions. The final part of the paper is devoted to determination of the front propagation law in Stefan-type problems at large times. Asymptotic solutions have been found for several important cases of the front propagation.

It is worth noting here that the non-field method, discussed in this study, yields the solution of Eq. (1) given in the form of Eq. (2) for any boundary conditions imposed on Eq. (1). Hence, there is no need to specify a set of boundary conditions for the said equation. Moreover, as can be easily seen from Eq. (4), the method can be used to transform the Direchlet boundary conditions into the Neumann ones and vice versa 12,13 .
Equations (2) and (4) can be obtained because the original transport Eq. (1) is reducible to the corresponding fractional partial differential equation where the operator is given as a series with respect to fractional derivatives.
The latter case is modelled by the heat equation, in which the coefficients α, β, and γ depend on time only. In such a case, the fractional operator can be expressed in an alternative form, which enables improved convergence of final results. Section "An alternative representation of fractional operators: improved convergence" discusses this alternative form of the fractional operator.
Section "Heat transfer through a moving boundary" is devoted to establishing the non-field solution of the heat equation, which models the process of heat transfer into a semi-infinite domain through the arbitrarily moving boundary. Some particular cases of the solution are considered in the end of the section. Among them are such cases as the boundary speed changing linearly, parabolically, exponentially, and polynomially.
In section "Determination of the front propagation law in Stefan-type problems at large times", some Stefantype problems are considered. Most of this section is devoted to determination of the front propagation law at large times. Asymptotic solutions have been found for several important cases of the front propagation.

An alternative representation of fractional operators: improved convergence
First, consider the case, when all the coefficients in the original Eq. (1) are constant, that is, α = const , β = const , and γ = const . In this case, the solution, obtained by the non-field method can be compared with the exact solution found by means of Laplace transforms. Indeed, Eqs. (4) and (7)  Notice that, in Eqs. (8) and (9), T s = 1 for the sake of simplicity but without lose of generality. Taking into account that L t [t ν ] = Ŵ(ν + 1)/p ν+1 , Eqs. (8) and (9) are identical. Upon invoking Eq. (3), the series in the right side of Eq. (8) can be summed up, that is, In comparison with Eq. (8), the latter expression is much more convenient, especially for large values of t. The latter result suggests that, in the case of α = const , β = const , and γ = const , the operator in the fractional representation of the transport equation, Eq. (2), can be written as For Eq. (2) with the arbitrary coefficients, the corresponding fractional form is where the coefficients a * n differ from those in Eqs. (4) and (6) but can be determined by the same methods as the coefficients a n . In the following sections, these coefficients are determined separately for each case considered.
As can be easily seen from Eq. (12), the corresponding non-field solution of the original transport Eq.
(1) at the surface becomes www.nature.com/scientificreports/ It is worth noting here that, for practical purposes, convergence of the series in Eqs. (12) and (13) for large values of t is better than in Eq. (6)-in particular, when α, β, andγ are slowly varying functions.
To conclude this section, it is worth noting that the result, presented above, is important from the methodological point of view. This is so, because the alternative form of the fractional operator given by Eq. (12) is established here for the first time. In comparison with the generalisation of the method reported reported earlier 18 , the alternative form warrants improved convergence of the method. The latter is explored in the second part of sections "Heat transfer through a moving boundary" and "Determination of the front propagation law in Stefan-type problems at large times".

Heat transfer through a moving boundary
Consider the process of heat transfer into a semi-infinite domain through the boundary, which moves arbitrarily as l = x 0 − t 0 β(ζ )dζ . Then, in the coordinate system related to the boundary, the heat transfer problem becomes identical to Eq. (1) with α = const , β = β(t) , γ = 0 . The latter problem is of great importance while modelling processes of propagation of the phase transition fronts, chemical reactions, or combustion.
If the coefficients in the original Eq. (1) depend on time only, as is the case for the problem in question, Eq. (5) reduces to Consequently, the recurrent expressions given by Eq. (7) simplify into Then, from Eqs. (4) and (15), the expression for the temperature gradient at the surface, q s = (∂T/∂x) x=0 , is in the form where dots denote derivatives with respect to time.
As an example, consider the case T s = 1 and β = C/ √ t with C = const. Taking into account that fractional derivatives of power functions have a very simple expression, namely 19 : from (16), it follows that As can be easily seen, in this case, all the terms in series (18) have the same power 1/ √ t. To obtain solutions in more complicated cases, it is much more useful to look for these solutions in the form with improved convergence.
To find the corresponding alternative form of the fractional operator (12) (20), follow the recurrent expressions for the coefficients a * n : Hence, the solution is given by Eqs. (13) and (22). To illustrate how the method works in practice, consider several particular examples. First, consider the case β = bt + d . Assume α = 1 in the original heat equation (this is easily achieved by rescaling of the spatial variable as x → x √ α). In this case, the fractional operator given by Eq. (19) satisfies Eq. (5) with where Hence, the solution is given by Eqs. (13), (23), and (24).
Next, consider the case, when β(t) is a decreasing function. If, for instance, the decreasing function β(t) is given in the form of a series with respect to exponential functions then it is possible to find the operator D from Eq. (5) in the form different from Eq. (6), which is the main advantage of the method leading to the form given by Eq. (19).
For the sake of illustration, consider a particular case Substituting D = F + β(t)/2 into Eq. (5) yields the expression for F The operator F is sought in the form ...  β n e −nct , c, β n = const, where the kernels K n can be found by Laplace transforms in the course of the solution procedure.
For the case given by Eq. (26), it is possible to find the first two terms in the series of Eq. (28) in a simpler way-directly from Eq. (16)-after noticing that, for c → ∞ , it is necessary to keep only those of them, which contain e −ct (except the first one). Thus, Treating the operator ∂ as a constant, Eq. (30) can be written as

the solution becomes
The equivalence between Eqs. (32) and (30) can be checked by using the product rule for fractional operators 12 , that is, To conclude this section, consider a very important case, when β → εβ(t) , where ε is a parameter. The equation for the fractional operator D in Eq. (14) becomes In this case, the operator D can be found in the form of a power series with respect to the parameter ε: Upon substituting Eq. (35) into (34) and equating expressions at the equal powers of ε , follow recurrent expressions for the operators D n : By direct substitution, it is easy to see that the expression for D 1 is given by a simple formula The solution of a more general operational equation can be found in the same way, rendering  www.nature.com/scientificreports/ The latter formula allows one to find finite expressions for all D n ( n > 1 ), if β(t) is a polynomial in t.

Determination of the front propagation law in Stefan-type problems at large times
It is known that problems of finding the front propagation speed can be modelled by a functional equation, which relates the dependent variable (e.g., temperature, mass concentration), its gradient at the moving boundary, and an auxiliary physical condition characterising the process under consideration 3 . For small and medium values of t , such a relation is given by Eqs. (16) and (19). An approximate solution of the corresponding inverse problem-determining β(t)-can be found out by using finite intervals of series in the mentioned equations. This section is devoted to determination of the front propagation law l = t 0 β(ζ )dζ for large values of t . As will be seen from the following analysis, the asymptotic equations relating T s , q s , and l are quite simple. The latter provides an efficient way to determine the front propagation speed for processes, which are characterised by a nonlinear boundary condition on the front of a phase or chemical change.
Consider Eq. (16) assuming that T s (t) and l (t) are monotonically increasing functions-T s (t) = const is also possible-and l (t) has all derivatives with respect to t . Assume also that, for t → ∞ , the surface temperature increases not faster than t µ , 0 ≤ µ < ∞ ; the derivatives of T s (t) increase not faster than a power function as well.
Assume now that, for t → ∞ , l (t)/ √ t → 0 (the front propagates slower than √ t ). From substituting l = t ν into Eq. (16) and taking into account Eq. (17), it follows that, for t → ∞ and 0 ≤ ν < 1/2 , the first term in the series dominates and, hence, The same conclusion can be drawn upon assuming l = √ t/ln ε t , ε > 0. The latter equation does not contain the front propagation law and is identical with the steady case (the boundary is at rest). This can be easily understood after recalling that, within a medium at rest, heat propagates as √ t ; hence, if the surface temperature increases quite slow, the slowly moving front has no influence upon the heat transfer process at large t.
In the case when l (t)/ √ t → ∞ for t → ∞ (the front propagates faster than √ t ), the values of the terms in the series of Eq. (16) increase with n . However, one of the terms in each multipliers at each fractional derivative ∂ (1−n)/2 dominates. Consequently, neglecting all non-dominant terms yields The latter equation can be written in the finite form as where the operators K 1 and K 2 are The symbol () * means that, while calculating ∂ 1/2 by Eq. (3), β 2 /(4α) is to be treated as a constant independent of t.
Because, for t → ∞ , e −ct ∂ 1/2 e ct T s = √ ∂ + cT s ≈ T s √ c if c = const and T s (t) is a power function, it follows that K 1 → β 2 √ α T s and, using the mean value theorem, In an intermediate case, when, for instance, the boundary moves by a parabolic law l (t)/ √ t → const for t → ∞ , all terms in the series of Eq. (16) are of the same order of magnitude. In this case, the relationship between T s and q s can be written in the form similar to Eqs. (40) and (43) The latter expression, in contrast with Eqs. (40) and (43) contains constant coefficients C 1 and C 2 , which can be found by the methods developed for the case β = · l ∼ 1/ √ t 3 . Now assume that, for t → ∞ , T s (t) increases exponentially as e αt or t µ e αt , 0 ≤ µ < ∞. First, consider the general case l (t)/t → 0 for t → ∞ (the front propagates slower than t ). Upon substituting l = t ν ( 0 ≤ ν < 1 ) into Eq. (16) and taking into account that ∂ ν t µ e αt → α ν t µ e αt , it follows that the first term dominates the entire series. Hence, the asymptotic solution is given by Eq. (40).
In the case when l (t)/t → ∞ for t → ∞ (the front propagates faster than t ), analysis of Eq. (16) renders the asymptotic solution given by Eq. (43). www.nature.com/scientificreports/ In the intermediate case of l (t)/t → β = const for t → ∞ (the front propagates with a constant speed), the problem has the straightforward solution The latter expression can be written in one of the forms given by Eq. (44), too.
Thus, the asymptotic expressions for the front propagation speed, for which the law relating T s , q s and β changes its form, have been obtained. The regimes change for various values of ν depending on how the surface temperature T s (t) increases.
To conclude this section, consider an illustrative example T s = At with A = const and the auxiliary condition on the boundary in the form where Q > 0 and k > 0 are both constants.
It is necessary to determine the front propagation law for t → ∞. Consider two assumptions, l / √ t → 0 and l / √ t → ∞. In the first case, combining Eqs. (46) and (40) 20 . However, the latter requires solving multiple intermediate and auxiliary differential equations. And, which is even more important, the knowledge of the entire temperature field is needed. On the contrary, the way, proposed here, yields an identical result from a single formula. In this lies a great advantage of the method employed in this study.
It is worth noting, however, that the regimes logarithmically close to the critical ones, for instance, l ∼ √ tln ε t or l ∼ tln ε t for T s ∼ t µ and T s ∼ t µ e αt , respectively, are not yet studied.

General discussion
The non-field solution to the problem of heat transfer through a boundary moving with a constant speed has been obtained by the method of Kulish earlier 14 . In a later study, the behaviour of this solution has been investigated to establish a set of conditions, under which chaotic solutions become possible 15 . However, the case of the boundary moving arbitrarily remained outside the study, mostly due to the fact that the series form of the fractional operator was not known until recently 18 .
Some attempts to apply the methods of fractional calculus to problems with phase change on the boundary 21 and combustion 22 have also been made. But again, the cases, considered in those studies, are rather special and lack the generality presented in this work. In particular, the solutions, reported in the mentioned works, are approximate, because, to find them, non-compact fractional operators with slow convergence have been used.
The extension of the method of Kulish, presented in this work, allows not only to establish a compact form of the fractional operator, which ensures much faster convergence of results, but also obtain non-field solutions, which relate local values of the temperature and its gradient, in the case of the boundary moving arbitrarily. In most of the cases, these solutions either cannot be obtained by other analytical methods, or the procedures leading to these solutions are too laborious (e.g., involve a large volume of calculations), which makes them impractical.
Among the cases investigated in this study are such cases as the boundary speed changing linearly, parabolically, exponentially, and polynomially. Whenever possible, the solutions thus obtained have been compared with known solutions. In addition, section "Determination of the front propagation law in Stefan-type problems at large times" is fully devoted to determination of the front propagation law in Stefan-type problems at large times.
Last but not least, it has been shown that, for Stefan-type problems, where applications of the non-field method are practically impossible due to a large volume of necessary computations, it is still possible to analyse the solution behaviour at large times and, in many cases, find asymptotic solutions. Obviously, nowadays all necessary computations can be fulfilled by built-in algorithms on public software (e.g., Matlab, Maple, Mathematica, Derive, etc.); however, in many practical engineering applications, the knowledge gained from asymptotic analysis is sufficient to provide desired estimates.