Modelling heat transfer in heterogeneous media using fractional calculus

This paper presents the results of modelling the heat transfer process in heterogeneous media with the assumption that part of the heat flux is dispersed in the air around the beam. The heat transfer process in a solid material (beam) can be described by an integer order partial differential equation. However, in heterogeneous media, it can be described by a sub- or hyperdiffusion equation which results in a fractional order partial differential equation. Taking into consideration that part of the heat flux is dispersed into the neighbouring environment we additionally modify the main relation between heat flux and the temperature, and we obtain in this case the heat transfer equation in a new form. This leads to the transfer function that describes the dependency between the heat flux at the beginning of the beam and the temperature at a given distance. This article also presents the experimental results of modelling real plant in the frequency domain based on the obtained transfer function.


Introduction
In many thermal problems, the temperature of the body is related to the heat flux. One example of such a problem is the temperature of electrical windings, which strongly depends on the heat flux generated by the currents  [1]. In such a case, the dynamical relation between the temperature and the heat flux is described by a partial differential equation (heat diffusion equation) [2]. In this paper, following the approach from Poinot et al. [1], a model based on fractional calculus of the heat diffusion process is given and its properties are discussed. The fractional (non-integer) calculus is an extension of the traditional differential calculus for a case when derivatives and integrals are in noninteger orders. This extension allows us, for example, to obtain more accurate results in modelling systems with complex structure and systems based on diffusion processes. The fractional order calculus was especially useful for modelling thermal processes [3][4][5][6]. In the study by Gabano & Poinot [5], the modelling and estimation algorithms, which were based on the state-space realization of fractional order integrator and numerical simulations, were presented. In the studies by Dzieliński and colleagues [6,7], experimental results of modelling the heat transfer process using fractional order calculus were presented. Another example of the systems that can be efficiency modelled using fractional calculus are ultracapacitors. The results of ultracapacitor modelling were presented by Dzieliński et al. [7]. A comparison between integer and fractional order modelling of ultracapacitors was presented by Dzieliński et al. [8].
While introducing the mathematical model of non-homogeneous beam heating, the effect of anomalous diffusion was taken into account (see [9] for details of anomalous diffusion, and [10] for subdiffusion parameter measurement). On the other hand, the control of such thermal plants has been discussed previously [11][12][13], in which a fractional order proportional-integralderivative-type controller has been employed to solve the temperature control problem. In this paper, the use of a fractional order calculus to model the heating process using a thermoelectric module also known as the Peltier module is presented. This paper is organized as follows. In §2, the basic facts on fractional calculus are described. Section 3 gives a mathematical description of the heating process. In §4, experimental verification of the heating process model is described. Section 5 concludes this article with some comments and ideas for further work.

Fractional calculus
A fractional order differential calculus is a generalization of the integer order integral and derivative to real or even complex order. This idea first emerged at the end of the seventeenth century and has been developed in the area of mathematics throughout the eighteenth and nineteenth centuries in the works of, for example, Liouville, Riemann, Cauchy, Abel, Grünwald and many others. More recently, by the end of the twentieth century, it turned out that some physical phenomena are modelled more accurately when fractional calculus is used. There exist two (in fact three) main definitions of the fractional order integrals, derivatives and differences: Riemann-Liouville, Caputo and Grünwald-Letnikov [3,14]. Some others are also present in the literature but are less commonly used in these applications. To be precise, the Riemann-Liouville and Caputo definitions concern both fractional derivatives and integrals.
To define the fractional order differ-integral, the definition of the Γ (x) function is needed. The Γ (x) function is given in the following way [3]: where (x) > 0.
In this paper, we will consider Caputo's definition, which can be written as [3] 0 for n − 1 < α < n. The initial conditions for the fractional order differential equations with the Caputo derivatives are in the same form as for the integer order differential equations. For the Caputo partial fractional derivative of order α of a function f (t, λ) with respect to variable t, we will use notation of the form ∂ α f (t, λ)/∂t α , which is often used in the related literature.

Mathematical description of the heating process
The ideal heating process, without energy loss, of a semi-infinite beam can be described by the following partial differential equation [2]: and with the following boundary conditions: where T(t, λ) is the temperature of the beam at time instant t and space coordinate (distance) λ, and 1/a 2 is a beam material conductivity.
Using the Laplace transform with respect to time and equation (3.1), we obtain which, using the boundary conditions, can be transformed to the following form: The solution of this equation has the form T(s, λ) = C 1 (s) e aλs 0.5 + C 2 (s) e −aλs 0.5 , (3.5) where C 1 (s) and C 2 (s) are the constants of the equation. The boundary conditions imply that C 1 = 0, because the solution must be bounded at ∞. This means that The heat flux H(t, λ) at time t and at distance coordinate λ is defined as follows: which using the Laplace transform and equation (3.7) gives the following relation: As a result, we achieve the following relation between the heat flux and the temperature at the desired point: (3.10) Using the above-mentioned steps, we can also try to describe the heating process of a semiinfinite beam as a subdiffusion process with the following partial differential equation: and with the following boundary conditions: where T(t, λ) is the temperature of the beam at time instant t and space coordinate (distance) λ, and a α is a parameter depending on beam parameters, such as heat conductivity and density; however, in fractional form, α is an anomalous diffusion order with relation α < 1, and H(t, λ) is the heat flux at time t and at distance coordinate λ. Again applying the Laplace transform with respect to time to equation (3.11), the following equation is obtained: For zero initial conditions T(0, λ) = 0, equation (3.12) can be rewritten as The solution of equation (3.13) is where C 1 (s) and C 2 (s) are the constants of the equation. The boundary conditions imply that C 1 = 0, because the solution must be bounded at infinity. This means that The heat flux H(t, λ) at time t and at distance coordinate λ is defined as follows: which using the Laplace transform and equation (3.16) gives the following relation: As a result, we achieve the following relation between the heat flux and the temperature for anomalous diffusion:

(a) Heat transfer equation for the heat flux loss case
In laboratory experiments of beam heating, it was found that the ideal equation does not describe the heat transfer process with good enough accuracy. This can be caused by non-ideal thermal insulation of the system from the environment, in which part of the heat flux is dispersed to the air around the beam. Additionally, when we assume that the beam is built from non-homogeneous material, the heating process is accompanied by an effect called anomalous diffusion. In such a case, the heat flux is described by the following equation: where α = 1 for anomalous diffusion (α > 1 for hyperdiffusion; α < 1 for subdiffusion, which is our case); the last part of the right-hand side of this equation represents the dispersed heat flux and b = a −1 . Hence, by first-order differentiation with respect to λ of both sides of equation (3.20), and by using equation (3.8), the following fractional order partial differential equation is achieved: By applying the right-hand side of equation (3.20), we obtain a fractional order partial differential equation of the following form: To obtain the solution of this equation in the transfer function domain, let us apply the Laplace transform with respect to time to equation (3.21), The solution of this equation (for H(0, λ) = 0) is given as follows: By applying the Laplace transform to equation (3.20), we obtain the following relation: For ease of analysis in practical implementation of e −Ts α , the following well-known power series expansion is frequently used: The power series is fast converging, thus the number of terms does not have to be too big. In fact, it can be reduced to just a few first terms. This reduction affects the accuracy of the implementation but makes the implementation and analysis easier and faster. We consider two examples of the approximation of the exponential transfer function. First, taking into account only three terms of the power series expansion, we achieved e −Ts α = 1 e Ts α ≈ 1 (T 2 /2)s 2α + Ts α + 1 , (3.29) and, second, taking into consideration four terms we achieved e −Ts α ≈ 1 (T 3 /6)s 3α + (T 2 /2)s 2α + Ts α + 1 .

Experimental verification of the heating process model
Heat distribution process modelling by fractional order partial differential equations and their respective counterparts in the frequency domain has been verified by experiments with a real physical thermal system. The results obtained from the mathematical model proposed have been compared with those obtained from the experiment.
The construction of the experimental set-up allows us to put a non-homogeneous material into the pipe and test materials with different parameters. The results obtained in the presented experiments used two types of materials, (i) buckshot (small metal balls of approx. 1 mm diameter) and (ii) a 50 : 50 mixture of buckshot and couscous groats (approx. diameter 2 mm).
The Peltier module was used as the source of the heat flux H(t, 0) at the beginning of the pipe (beam). The power amplifier was used as a voltage to current converter in order to easily control the current in the Peltier module. All sensors and control signals for the power amplifier were connected to the dSPACE control card DS1104.

(b) Thermoelectric module: Peltier module
The Peltier module is a semiconductor device that allows us to control not only the value of the heat flux but also the direction of that flux. In the experiments, the TM-127-1.0-3.9-MS Peltier module, produced by SCTB NORD, was used. The main parameters of the module are as follows:  A Peltier module can be described by the following equation [15]: where H is the value of heat flux on the heat side of the module, α, R and K p are the Peltier module parameters, I is the current applied to the module, and T h and T c are the temperatures of both sides of the module (hot and cold sides). A simplified diagram of the Peltier module used is given in figure 2.
Based on the module characteristics provided by the producer the following parameters were obtained: R = 4.8879 Ω, c = 0.0565, K p = 0.4978. From equation (4.1) the following relation for obtaining the module current for the desired value of heat flux is then achieved: As can be seen in figures 3 and 4, the accuracy of the modelling is very high for both cases. Only at the end of the magnitude plots can some significant values of error be observed. These can be caused by measurement errors because measured gains for these frequencies are very low (approx. −30 dB).

(d) Model validation in the time domain
In this section, the validation of the obtained frequency domain models in the time domain is presented. For numerical realization of the fractional order derivatives, the recursive Oustaloup's algorithm [14] was used (in the implementation included in [16]). The parameters of the Oustaloup's algorithm that were used were ω l = 0.001, ω h = 1000 and N = 10. Additionally,  ≈ k e −λ 1 1/k as α/2 + 1

Conclusion
In this paper, the results of modelling the heat transfer process in heterogeneous materials are presented. The mathematical model in the form of fractional partial differential equations for a non-homogeneous beam heating process was established. The experimental set-up for beam heating with a Peltier module as the heat flux source was built. This set-up allowed us to study beams containing different materials. Two of them (buckshot and a mixture of buckshot and couscous groats) were used in experiments presented in this paper. For both cases, the parameters of the mathematical models were found. The results obtained from measurements were compared with those resulting from the mathematical model. The comparison shows, in general, a good match between the mathematical model and the physical experiment.
In further work, we will verify the results obtained in the frequency domain with the time domain experimental measurements by using a well-known method for the numerical simulations [3,17].