General one-dimensional model of the time-fractional diffusion-wave equation in various geometries

This paper deals with the analysis of the time-fractional diffusion-wave equation as one-dimensional problem in a large plane wall, long cylinder, and sphere. The result of the analysis is the proposal of one general mathematical model that describes various geometries and different processes. Finite difference method for solving the time-fractional diffusion-wave equation using Grünwald-Letnikov definition for homogeneous or inhomogeneous material and for homogeneous or inhomogeneous boundary conditions is described. Dirichlet, Neumann and Robin boundary conditions are considered. Implementation of numerical methods for explicit, implicit, and Crank-Nicolson scheme were realised in MATLAB. Finally, illustrative examples of simulations using the developed toolbox are presented.


Introduction
The fundamental phenomena, such as diffusion and wave propagation, are described by the following partial differential equations [4], i.e.
Generalising the above partial differential equations by replacing the integer-order time derivatives with derivatives of fractional order α ∈ (0; 2] the following timefractional diffusion-wave equation is obtained: The standard (integer-order) diffusion or wave equation used for the description of many materials such as visco-elastic materials, granular and porous materials, composite materials, etc., is often not sufficiently adequate. In such cases, the description requires the development of more suitable models using fractional-order derivatives [13,22,28,29,44]. The causes are mainly roughness or porosity of the material [21][22][23], fractality and chaotic behavior of systems [3,16,24,42], and memory of the systems and ongoing processes [25,26,35,36,41].
The time-fractional diffusion-wave equation as an adequate model requires the creation of new methods and tools for its solution. This topic is currently relevant and many important authors have contributed to its development, such as Fourier, Abel, Leibniz, Grünwald, Letnikov, Liouville, and Riemann [18,22,34].
Therefore, the present work aims at the design, numerical solution and MAT-LAB implementation of the general one-dimensional model of the time-fractional diffusion-wave equation applied to the time-fractional heat conduction problem in various geometries. This is presented in Fig. 1, which can be considered as an outline and a graphical abstract.

Time-fractional heat conduction models
The one-dimensional time-fractional heat conduction model in different geometries can be expressed as follows: -in a large plane wall -and in a sphere The examination of the one-dimensional time-fractional heat conduction model in large plane wall, long cylinder, and sphere reveals, that all three equations can be expressed in a compact form as a general model In the case of a constant value of the thermal conductivity λ, the general model has the form where γ is a geometric factor (large plane wall γ = 0, long cylinder γ = 1, sphere γ = 2), k 2 = λ/( c p ) is the thermal diffusivity, is density, and c p is specific heat capacity.

Initial and boundary conditions
The general time-fractional heat conduction model (2.2) requires two boundary conditions, as well as one initial condition. The initial condition specifies the temperature distribution in the material at the begin of the time, that is, The boundary conditions specify the temperature or the heat flux at the boundaries of the region. The following boundary conditions ( where h a , h b are heat transfer coefficients (W·m −2 ·K −1 ); T sa , T sb are the surrounding medium temperatures (K) at the boundary surface.

Numerical solution
The fractional derivative according to time in the general time-fractional heat conduction model (2.2) is discretised by the backward Euler method and the Grünwald-Letnikov definition Δt α by using the principle of "short memory" (Fig. 3) [26], where t L is the memory length, Δt is the time step and the value of N f shall be determined by the following relation For the calculation of the binomial coefficients c j the following relation can be used The derivative according to the spatial coordinate ( Fig. 4) is discretised as follows ∂ T (r , t) ∂r The values ξ i or ζ i represent the coefficients for calculation of the area A i+ 1 2 or A i− 1 2 for the selected geometry as follows Similarly, we can express the relation between area and volume as

Steady-state
The general time-fractional heat conduction model (2.1) for inhomogeneous material and under steady-state conditions (∂ T /∂t = 0) takes the form and the matrix form is as follows where the coefficients are The specified temperatures of the boundary conditions are incorporated by simple assigning the given surface temperatures to the boundary nodes and therefore for the Dirichlet boundary conditions the coefficients are When other boundary conditions such as the heat flux or convection are specified at the boundary, the finite difference equation for the node at that boundary is obtained by writing an energy balance of the elementary volume at that boundary (Fig. 5). The energy balance is expressed as and therefore For the Neumann boundary conditions the coefficients are where R a , R b are internal thermal resistances (m 2 ·K·W −1 ).
In the case of Robin boundary conditions the energy balance is expressed as For the Robin boundary conditions the coefficients are where β a , β b are Biot numbers or ratios of internal and external thermal resistances at the boundary surface.

Unsteady-state
According to the type of differential expression finite difference methods can be divided into explicit, implicit and Crank-Nicolson schemes.

Explicit scheme
The explicit scheme for the general time-fractional heat conduction model (2.2) has the form where module τ i is determined by the relation while the value of τ i with respect to the stability of the solution must have values less or equal to 0.5. The matrix form of equation (4.1) is as follows where the coefficients are For the Dirichlet boundary conditions the coefficients are In the case of Neumann boundary conditions the energy balance is expressed as and therefore For the Neumann boundary conditions the coefficients are In the case of Robin boundary conditions the energy balance is expressed as For the Robin boundary conditions the coefficients are where R a , R b are internal thermal resistances (m 2 ·K·W −1 ); β a , β b are Biot numbers or ratios of internal and external thermal resistances at the boundary surface.

Implicit scheme
The implicit scheme for the general time-fractional heat conduction model (2.2) has the form

In matrix form this is
where the coefficients are For the Dirichlet boundary conditions the coefficients are For the Neumann boundary conditions the coefficients are For the Robin boundary conditions the coefficients are

Crank-Nicolson scheme
The Crank-Nicolson scheme for the general time-fractional heat conduction model (2.2) has the form

In matrix form this is
where the coefficients are For the Dirichlet boundary conditions the coefficients are For the Neumann boundary conditions the coefficients are For the Robin boundary conditions the coefficients are

Implementation
The implementation was realised in the programming environment MATLAB  TFDWEg_CN (Uin,alpha,Nf,k,dr,a,b,dt,ts,TBC,BC,g) where the function inputs are: Uin is the vector of the initial conditions, alpha is the order of the time derivative, Nf is the memory length, k is the process coefficients vector, dr is the spatial step, a and b are the left, right side coordinates of the body, respectively, dt is the time step, ts is the time of simulation, TBC is the vector of the boundary conditions type, BC are the values of the boundary condition, and g is the geometry type. The function outputs are: R is the vector of coordinates, T is the vector of time, and MU is the output values matrix. The toolbox is published at MathWorks, Inc., MATLAB Central File Exchange [40].

Simulations
Time-Fractional Diffusion-Wave Equation simulations in various geometries for homogeneous or inhomogeneous material and for homogeneous or inhomogeneous boundary conditions are realised using the script TFDWEg_test.m which is a part of the toolbox TFDWEg [40]. In the following examples the behaviour of temperature in space and time for various geometries is illustrated.
Let us propose an example, where the initial conditions are been defined in the form of the constant temperature (T (r , 0) = 0 • C) in the whole cross-section, considering copper as the material (with density = 8900 kg/m 3 , thermal conductivity λ = 400 W/m/K, specific heat capacity c p = 380 J/kg/K), left (a = 0.01 m) and right (b = 0.03 m) side coordinate. The material was divided into twenty parts in the space, with time step (0.001 s), total time simulation (2 s), and Dirichlet boundary conditions with temperatures at the edges (T a = 0 • C and T b = 15 • C). The simulation results using Crank-Nicolson scheme are shown in Fig. 6, where the left coloumn shows the behaviour of temperature in time (where the parameter is the position in space), the middle coloumn shows the behaviour of temperature in space (where the parameter is time), and the right coloumn displays 3D plots of temperatures in space and time for hollow body (Fig. 2a), i.e. large plane wall (Fig. 6a), long cylinder (Fig. 6b), and sphere (Fig. 6c). In the case of setting the left coordinate to zero (a = 0.0 m), it is possible to perform the same simulation but for a full symmetrical body (Fig. 2b), as shown in figure Fig. 7.
The graphs in Fig. 8 show the evolution results solving the time-fractional heat conduction using the fractional-order derivative α = 0.5, 1.0, and 1.5, respectively. Comparing the graphs, one can observe that the derivative of order 0.5 exhibits fast temperature rise in the beginning and slow temperature rise later. Moreover, it is evident that the temperatures propagate and diffuse with time, which means that the temperature continuously depends on the fractional derivative, therefore, when α = 1.5, both diffusion and wave response can be observed.
Let us use the previous settings to perform simulations demonstrating the use of different homogeneous boundary conditions for homogeneous material. The behaviour of the temperature for the Dirichlet, Neumann, and Robin boundary conditions is shown in Fig. 9. In the case of the Dirichlet boundary conditions, the surface temperatures are T a = T b = 15 • C; in the case of the Neumann boundary conditions, the thermal flux and the internal thermal resistance are i a · R a = i b · R b = 0.75 • C, and in the case of the Robin boundary conditions the surrounding media temperatures and the ratio of internal and external thermal resistances are T sa = T sb = 15 • C, β a = β b = 1. The graphs in Fig. 10 show inhomogeneous boundary conditions for selected combinations, i.e. Neumann-Dirichlet, Robin-Dirichlet, and Neumann-Robin, respectively.
The script TFDWEg_test.m also allows to simulate homogeneous as well as inhomogeneous material. For example, in the case of thermal diffusivity for twenty layers and two materials (brass and copper) of the same thickness a twenty-component thermal diffusivity vector is needed (ten components for brass and ten components for copper). Results of the simulations are shown in Fig. 11a for homogeneous and in Fig. 11b for inhomogeneous material, where in the first half of the specimen length there is a lower   temperature rise than in the other half. This is due to the fact that in the first half the material there is copper whose thermal diffusivity is about three times larger than that of brass.

Conclusion
The general one-dimensional model of the time-fractional diffusion-wave equation in various geometries is designed and described in this contribution. The numerical solution using the Grünwald-Letnikov definition of fractional derivative according to time for homogeneous and inhomogeneous boundary conditions and for homogeneous and inhomogeneous material is presented. The numerical schemes were implemented in MATLAB in the form of a library of functions. The possibilities of using this library are illustrated by examples and simulations. Additionally, besides modeling heat conduction processes, the library allows to model other processes such as diffusion, wave distribution, and so on. The library is a suitable tool for the implementation of simulations of a wide class of processes and also for the creation of complex models in the MATLAB environment. Benefits of this work are the design, numerical solution and implementation of one general mathematical model that can be used for modeling various different processes, in various geometries, materials and boundary conditions.