Modeling of Heat Distribution in Porous Aluminum Using Fractional Differential Equation

The authors present a model of heat conduction using the Caputo fractional derivative with respect to time. The presented model was used to reconstruct the thermal conductivity coefficient, heat transfer coefficient, initial condition and order of fractional derivative in the fractional heat conduction inverse problem. Additional information for the inverse problem was the temperature measurements obtained from porous aluminum. In this paper, the authors used a finite difference method to solve direct problems and the Real Ant Colony Optimization algorithm to find a minimum of certain functional (solve the inverse problem). Finally, the authors present the temperature values computed from the model and compare them with the measured data from real objects.


Introduction
Fractional calculus is a part of mathematical analysis and has a lot of applications in technical science.One of the most popular books about fractional calculus is Reference [1].References [2][3][4] provide information about fractional calculus, fractional differential equations, approximations of fractional derivatives and numerical methods.There are also a lot of articles about fractional calculus-for example, [5][6][7].
Various phenomena in nature can be modeled using fractional derivatives [8][9][10][11][12][13][14][15][16]-for example, in [9,11], the authors surveyed fractional-order electric circuit models, Reference [12] shows applications of fractional derivatives in control theory, and, in [13,14,16,17], we can find information about application fractional derivatives in heat conduction problems.In [16], the authors present an algorithm to solve the fractional heat conduction equation.In the presented model, the heat transfer coefficient is reconstructed based on measurements of the temperature.The direct problem is solved by using the implicit finite difference method.To minimize the functional defining the error of the approximate solution, the Nelder-Mead algorithm is used.In [18], the authors consider the inverse problem of recovering a time-dependent factor of an unknown source on some sub-boundary for a diffusion equation with time fractional derivative.The authors present two regularizing schemes in order to reconstruct an unknown boundary source.Another paper where authors solved the inverse problem with fractional derivatives is Reference [19].Zhuang et al. considered a time-fractional heat conduction inverse problem with a Caputo derivative in a three-layer composite medium.To solve the direct problem, they used a finite difference method, and, for the inverse problem, the Levenberg-Marquardt method was applied.The results show that the time-fractional heat conduction model provides an effective and accurate simulation of the experimental data.In addition, in [20], it is considered an inverse fractional heat conduction problem.The authors show that the model with fractional derivatives better describes the process of heat conduction in ceramics media.
In this paper, the authors consider the heat conduction inverse problem.A mathematical model describing the heat transfer phenomenon in porous aluminum is given by fractional differential equation with initial-boundary conditions.In this case, we used the Caputo fractional derivative.The algorithm consists of two parts: solution of direct problem and solution of inverse problem by finding the minimum of the functional.Additional information for the inverse problem was the temperature measurements obtained from porous aluminum.The direct problem was solved using a finite difference method and approximations of Caputo derivatives [17,21].In the inverse problem, the heat transfer coefficient, thermal conductivity coefficient, initial condition and order of derivative were sought.In order to do that, we need to minimize the functional describing the error of approximate solution.The functional was minimized by a Real Ant Colony Optimization algorithm [22,23].
More about heat conduction inverse problems can be found in [24][25][26][27][28]. Zielonka et al. in [25] solved the one-phase inverse problem of alloy solidifying within the casting mould.The authors also include their paper shrinkage of the metal phenomenon.
The investigated inverse problem consists of reconstruction of the heat transfer coefficient on the boundary of the casting mould on the basis of measurements of temperature read from the sensor placed in the middle of the mould.In [28], Stefan problems relevant to burning oil-water systems are considered.The author used the heat balance integral method.

Fractional Heat Conduction Equation
In this section, we would like to present the fractional differential equation defined in region D = {(x, t) : x ∈ [0, L], t ∈ [0, T], L, T ∈ R + }, where c is specific heat, is density of material, λ is thermal conductivity coefficient, α is order of derivative and u is the function describing the distribution of temperature.In literature, this equation is called the Time Fractional Diffusion Equation (TFDE).For a more precise description of the model, we still need initial-boundary conditions.
In this case, we used a Neumann boundary condition for x = 0, and a Robin boundary condition for x = L. Below, we present initial-boundary conditions: where f is function describing initial condition, q is the heat flux, h is heat transfer coefficient and u ∞ is ambient temperature.By solving models ( 1)-( 4), we obtain the temperature values at the points of the domain D. To model process of heat conduction in porous media, we used Caputo fractional derivative of order α ∈ (0, 1), which, in our case, is defined by the formula [1]: In the next part of the article, we use the considered model to solve the fractional heat conduction inverse problem.

Formulation of the Problem
In an inverse problem, some information in the considered model is unknown; in our case, we do not know: thermal conductivity coefficient λ, heat transfer coefficient h, initial condition f and order of derivative α.We have additional information, called input data, which is the temperature measurements.We denote it by: where N 1 is number of measurements from thermocouples.
If we solve the direct problems ( 1)-( 4) for fixed values of the sought parameters, then we obtain calculated values of the temperature at certain fixed points of the domain D-in our case, we denote it by U j (λ, h, f , α).Using the calculated values of the temperature U j (λ, h, f , α), input data U j , we create functionals defining the error of approximate solution: By minimizing the functional (7), we find the approximate values of the sought parameters.
The measured temperatures (input data) U j were obtained from a sample of porous aluminum.This sample was formed by pressurizing the powders' aluminum of medium size 0.8 mm in the plate hydraulic press.Powders were pressed at 150 bar pressure.Sample was heated to 300 • C at speed of 1 K/s and then cooled to ambient temperature.During that time, sample temperature measurements were done, which measurements were used as input data for algorithm.

Method of Solution
The solution of the inverse problem can be divided into two parts: first-solution of direct problems (1)-( 4), and second-finding minimum of the functional (7).

Solution of the Direct Problem
In solving the inverse problem, we have to solve the direct problem many times.To solve direct problems (1)-( 4), we used implicit finite difference scheme.In order to do that, we create grid with size (N + 1) × (M + 1) and steps ∆x = L/N, ∆t = T/M.Caputo fractional derivative ( 5) is approximated by the formula [17]: where We also need to approximate boundary conditions of the second and third kinds.The following approximations were used: Using all approximations (8)-( 10) and the differential quotient for the derivative of second order with respect to space we get the following differential equations where is the thermal diffusivity coefficient.Solving the system of equations gives us approximate values of function u in points of grid S.More about stability of the presented method can be found in [17].

Minimum of the Functional
The second part of the presented algorithms is finding the minimum of the functional (7).In this paper, we used the Real Ant Colony Optimization algorithm (Algorithm 1).The algorithm was inspired by the behavior of swarm of ants in nature.Pheromone spots, which are L, are identified with solutions.Firstly, they are distributed randomly in the considered area.Then, we rank them according to their quality-better solution means stronger pheromone spot and greater probability to choose it by ant.In this way, we create the solution archive.In every iteration, the one of M ants constructs one new solution (new pheromone spot) using the probability density function (in this case, Gaussian function).The ant chooses with the probability a pheromone spot (solution) and transforms it by sampling its neighborhood using the Gaussian function.Then, the solutions archive is updated with new solutions and sorted according to the quality; next, M worst solutions are rejected.The described algorithm was adapted for parallel computing.For the description of the algorithm, we will introduce the symbols: Setting input parameters of the algorithm L, M, I, nT, q, ξ.
2. Randomly generate L pheromone spots (solutions) and assign them to set T 0 (starting archive).
3. Calculate values of the minimized function F for each pheromone spot and sort the archive T 0 from best to worst solution.

Iterative process
4. Assigning probabilities to pheromone spots (solutions) according to the following formula: where weights ω l are associated with l-th solution and expressed by the formula 5. Ant chooses a random l-th solution with probability p l .6. Ant transforms the j-th coordinate (j = 1, 2, . . ., n) of l-th solution s l j by sampling proximity with the probability density function (Gaussian function) , Repeat steps 5-6 for each ant.We obtain M new solutions (pheromone spots).8. Divide new solutions on nT groups.Calculate values of minimized function F for each new solution (parallel computing).9. Add to the archive T i new solutions, sort the archive by quality of solutions, remove M worst solution.10.Repeat steps 4-9 I times.

Results
In this section, we present the obtained results.We consider models (1)-( 4) with the following data: In the presented model, we lack the following data: which is reconstructed based on measurements of temperature.We assume that Solving the direct problem, we used two different grids 100 × 1995 (∆x = 0.03825, ∆t = 0.036) and 100 × 3990 (∆x = 0.03825, ∆t = 0.018).In the case of the real ACO algorithm, we used the following parameters: L = 12, M = 16, I = 55, nT = 4, q = 0.9, ξ = 1.
The real ACO algorithm is probabilistic, so we decide to execute them ten times to check the stability.
Table 1 presents values of reconstructed parameters a 1 , a 2 , . . ., a 6 in the case of two different grids.In both grids, parameters have similar values except parameter a 1 -the thermal conductivity coefficient.For 100 × 1995 grid, a 1 is equal to 300 and, for 100 × 3990 grid, the parameter has value 237.91.The initial condition is approximately equal to 569 K and the order of derivative is equal to 0.20 (100 × 1995) or 0.21 (100 × 3990).As we can see in Table 1, if the 100 × 1995 grid is used, then, for the most of sought after parameters, the standard deviation is less than for the 100 × 3990 grid.This means that, for the first grid, the obtained results from the ACO algorithm were slightly less dispersed than in the case of the second grid.In Figure 1, we can see plots of reconstructed function h for two grids.Both plots are similar.The values of the function for the grid 100 × 1995 are slightly larger than for the second grid.Table 2 presents errors of reconstruction temperature in measurement points.In both cases, these errors are similar and temperature is reconstructed very well.Average errors are a little bit smaller in the case of 100 × 3990 grid, but, in the case of maximal errors, it is otherwise.At the end of this section, we would like to present plots of reconstructed temperature and distribution of errors of this reconstruction (Figures 2 and 3).Temperature reconstruction and distribution of errors are very similar.

Conclusions
In summary, the paper presents a fractional heat conduction equation with Caputo derivative.Based on this equation, the inverse problem is solved using temperature measurements from porous aluminum.Next, computed values of temperature were compared with real data.First of all, the results obtained were very good.Average relative error of reconstruction was equal to 1.06% for 100 × 1995 grid, and 1.02% for 100 × 3990 grid.More density grid gives us smaller relative errors, but a little bit higher maximal errors.Using less density, grid results obtained from the ACO algorithm were less dispersed.The sought order of fractional derivative is approximately 0.20.In the future, we would like to take under consideration other models of heat conduction using different fractional derivatives.

FAlgorithm 1 :
(x) minimized function, x = (x 1 , . . ., x n ) ∈ D n dimension (number of variables) nT number of threads M = nT • p number of ants in population I number of iterations L number of pheromone spots q, ξ parameters of the algorithm Parallel Real ACO algorithm Initialization of the algorithm 1.

Figure 2 .
Figure 2. (a) Measurements of the temperature (black line) and reconstructed temperature (blue line); and (b) distribution of errors for this reconstruction (100 × 1995) .