Numerical Modeling and Analysis of Physical Properties in Biomaterials with Fractal Structure

Abstract

The synthesized two-dimensional mathematical models of non-isothermal humidity transfer in media with fractal structure taking into account the effects of memory and spatial correlation. Built explicit and implicit difference schemes for equations related heat-mass transfer in two-dimensional domain with boundary conditions of the third kind. These algorithmic aspects for the realization of the obtained difference equations using method predictor-corrector and analyzed the conditions of stability of difference schemes.

Keywords

• heat-mass transfer
• fractal structures
• difference schemes
• derivatives of fractional order
• stability
• approximation fractional derivatives

-

Actuality of the topic. The development of the theory and methods of mathematical or computer modeling of processes and systems in various fields of human activity has always been based on the use of new ideas, approaches from the field of analysis, applied and computational mathematics. One of the most important tasks that arises during modeling is the adequacy of a mathematical model for an object or phenomenon from the real world.

Dynamic systems, as an object of modeling, have traditionally been studied using integro-differential equations in integer and fractional derivatives. Classical analysis assumes that integrals and derivatives have integer orders. Nevertheless, it has already been discovered by observing that the behavior of a number of objects and processes does not fully correspond to its mathematical models (with integrals and derivatives of integer orders), that is, in some models the question of adequacy is raised. Accordingly, models and their solutions for the dynamical systems with integrals and derivatives of real orders began to be developed. The notion of an integral and derivative of noninteger orders lie at the basis of the integral and derivative of fractional orders.

At the moment, the fractional calculus is at a stage of great development, in the theoretical and practical application. One of the main questions faced by scientists is the interpretation and application of integration and differentiation of fractional orders for different models. For today it's hard to say what an integral or derivative of a fractional order is in terms of geometric and physical interpretations. Nevertheless, this section of mathematical analysis has become a tool for mathematical modeling of complex dynamic processes, in ordinary and fractal environments, and allows us solve various problems.

Here is an incomplete list of tasks in which the fractional derivatives are effective: automatic control; signal processing; physics and electronics; biology and medicine; economy and finance; classical mechanics; hydrodynamics (movement of the body in a viscous liquid); thermal conductivity (dynamics of heat flows); diffusion processes; dynamics of turbulent environment; visco-elasticity (rheology of polymers); static optics; radio-physics; radio engineering; dynamic chaos; and other…

However, despite the practical application of fractals and their analysis, there are still many unresolved problems and problems associated with the geometric, physical and probabilistic interpretations of the fractional derivatives apparatus and integrals. 2

Analysis of mathematical apparatus for differentiation of fractional order Fractal calculus has involved a lot of famous scientists after Lopital and Leibniz. Fourier, Euler, Laplace - the most famous among those who were engaged in this mathematical apparatus. Many of them introduced their own notation and methodology that would match the concept of integration and differentiation of fractional order.

Typically, for the description of non-stationary processes, the integration and differentiation operators are used, which determine the overlay of certain conditions on processes and generalize their properties. Today, in many branches of science, there are new structures for which the use of ordinary differential equations is insufficient. Instead, they could be adequately described with the aid of the mathematical apparatus of integration and differentiation of the fractional order. A fractal calculus is called the domain of mathematical analysis, where the operators of differentiation and integration of any real order are studied [ 3, 6, 8, 9, 10, 12, 19 ]. In the last few decades there was an urgent need to use this apparatus in various fields of science, such as: classical and quantum physics, field theory, solid state physics, fluid dynamics, aerodynamics, stochastic analysis, image processing [ 1, 4, 5, 11, 13, 20 ].

Equations containing integrating or differentiating operators of fractional order are widely used to describe the behavior or state of a real physical environment or process. There are many phenomena and processes that have a characteristic fractal or memory effect. The mathematical apparatus of integration and differentiation of fractional order is the best method for constructing models of such systems [18, 21]. The memory mechanism may be different depending on the type of process, by the way the phenomenological description of many processes with this property may have one basis. A fractional calculation in the theory of such systems becomes irreplaceable, which could be compared with the classical analysis in continuum mechanics [ 2, 7 ].

Recently, operators of integration and differentiation of fractional order in the theory of visco-elasticity are widely used [ 8, 9, 14, 16 ]. The use of operator data to describe the relations between stresses and deformations made it possible to take into account the existence of irreversible phenomena due to the rheological properties of the material [ 8, 16, 19 ]. Research of the stress-strain state in visco-elastic bodies play an important role in estimation of their strength and reliability during technological processing.

An analysis of scientific sources suggests that the definition of derivatives of fractional order is based mainly on three approaches. The first is based on the generalization of the well-known Cauchy formula, which allows us to construct a multiple integral of an integer order to a single [ 8, 9, 19 ]. The second approach is developed in the works [ 6 ] and generalized in [ 7, 20 ] about the definition of a fractional derivative by the boundary of a finite-difference relation. There are also known a number of generalizations and modifications of such approaches. The main difference of fractional derivatives from integers is their non-locality, that is, the dependence of the results of differentiation on the values of functions at all points of a certain segment or numerical line, and not on the values of functions at points from the small circle of a given point - as in the case of ordinary differentiation. Also known studies on the generalization of fractional differentiation operators, in particular [ 8, 21 ], the fractional order is described by the function of time, and in [21] by a random variable.

There are various options for introducing integration and differentiation operations of fractional order, in particular Riemann-Liouville, Kaputo, Grunwald-Letnikov's approaches, and their various modifications [ 8, 9, 19 ]. On certain classes of functions, these operations lead to identical results. As an example, we can give a fractional integration and differentiation of a completely integrable function on a finite segment of Riemann-Liouville, which coincides with the corresponding operations for GrunwaldLetnikov [20]. Let's consider more detailed modifications of these operators.

Studies are devoted to the construction of mathematical models and software for physical and mechanical fields in capillary-porous materials with fractal structure. Such fractional order models describe the evolution of physical systems with residual memory and the very similarity of a fractal structure that occupy an intermediate position between Markov systems and systems that are characterized by complete memory. In particular, the fractional indicator indicates the share of system states that are stored throughout the process of its operation 3

Mathematical model for transferring heat and humidity in environments with fractal structure The mathematical model for heat and humidity transfer in an environment with a fractal structure is described by a system of differential equations in partial derivatives with a fractional order over time t and a spatial coordinate x. with initial conditions: c  T  , x 

  U  , x   

   a  T  , x 

x  U  , x  x   0r U  , x 

   a  T  , x  x

, T U  0  0  T  x  ,

0  U 0  x  , and boundary conditions of the third type:   T x x0,l a  T x   0 1     U x0,l  U p     T

x0,l  tc  ,  a  U x x0,l x0,l    U p  U x0,l T;U - unknown functions, T - temperature, U - humidity, c - specific heat capacity, ρ density, λ - coefficient of thermal conductivity, ε - coefficient of phase transition, ρ0 basic function, r - specific heat of steam generation, a - coefficient of conductivity, δ thermo gradient coefficient, Tc - temperature of the environment, Up - relative humidity of the external environment, σ - humidity transfer coefficient, ω - heat transfer coefficient, 0    1 - fractional order of derivative over the time, 1    2 , 0    1 fractional order of derivative over the spatial coordinates. 4

Numerical algorithm Let's make space-time discretization in the domain D:   ,h  { k , x(n)  : x(n)   n 1 h, k  k  , n  1,..., N ; h  k  0,1,..., K ;   t K

}.

Using the Riemann-Liouville formula: l N 1 ; (7) (1) (2) (3) (4) (5) (6)  f      k 

1  1      k 1  k  f  k 

 k 1  k  k|1    f  

 d  ,   We can write difference approximation of fractional derivative of order α 0   1 In the case when   1 , we obtain an explicit finite-difference scheme, and when   0 - an implicit scheme with discrete initial conditions:

Tn0  T0  x(n)  , U n0  U 0  x(n)  , where h  xn1  xn ,  h -integer part of a number x.

Then the difference approximation of the fractional derivative β for the coordinate x will look like:  u x x(n)

1  h n  q jun j1 , j0 where q0  1, q j   1 j    1...  j  1 .

j!

Taking into account (9), (11) we will obtain an explicit difference scheme for the numerical implementation of the system of differential equations (1),(2): с

Tnk 1  Tnk    2      h n  q jTnkj1  0r j0

U k 1 U k

n n ,   2     , on line segment [ k , k 1 ] :  u    k 

u k 1 u k 2    ,    k 1  k , where   - Gamma function.

To determine the fractional derivative of order from one to two 1    2 , we can use the Grunwald-Letnikov formula:  f   x  lhim0 h

  1 h  1 j j0

   1   j  1    j  1 f  x  j  1, (8) (9) (10) (11) (12) (13)

  0 1     U1k 1  U р     T1k 1  tс  , a a

T k 1  T1k 1 2   2   h  a

U k 1  U k 1

2 1   2   h    U р U k 1  , 1  TNk21  TNhk11   0 1     U Nk 1  U p     TNk 1  tc  ,

TNk 1  TNk11  a U Nk1  U Nk11    U   2   h   2   h p U k 1 .

N

To find the numerical solutions of the resulting system of difference equations, use the predictor-correction method. In the role of the predictor, we will use an implicit difference scheme (12), (13) with the step ht/2, and in the role of the corrector will be an explicit difference scheme (12), (13).

с

T k 12  Tnk

n   2      2    n   q jTnkj121  0r h j0

U k 12 U k

n n   2     

 2 , U k 12 U k

n n     2      2 a n   q jUnkj121  h j0 a n k 12 . h j0 q jTn j1 Equation (17) will be rewritten in the form for convenience:

 с U k 12   A1q0Tn1

k 12   n  0r  A1q1 Tnk 12  A1 q jTnkj121  Unk  0r Tnk  , (19)

n  с  j2   where A  1  2    2 

 0rh

U1k 12    0 1   

  B1 T1k 12  B1T2k 12  U р 

 tс  0 1   

Similarly, we write in the matrix form the equation (18) and the boundary conditions (15), (16), which correspond to it: C  cij , B  bij , i, j  1, N .

BTnk  12  CU nk  12   2 U k  0 , n where B  1 0 1   2   h T Unk  0,U2k ,...,U Nk 1, 0T , Tnk  0,T2k ,...,TNk1, 0T , 1  U p  0 1tc   ,0,...,0,U p  0 1tc   .  The components aij , i, j  1, N of the matrix A are determined by expressions:       aij    c     0r

0, j  i  2; 0, i  N ,1  j  N  2;

  A1q1  , i  j  1  N;

      0 1  

  B1  , i  j  1;       0 1  

  B1  , i  j  N;

 B1, i  1, j  2; B1 , i  N , j  N 1;  A1qi j1, інше.

(22) (23)  0, j  i  2;  0, i  N ,1  j  N  2;   Z1q1 1 , i  j  1  N ;    a     2   h  , i  j  1;  cij     a     2   h  , i  j  N ;  a, i  1, j  2;  a , i  N , j  N 1;   

Z1qi j1, інше. a  2     2    0, j  i  2;  0, i  N ,1  j  N  2;  Zq1 , i  j  1  N ; bij  a , i  j  1; i  N , j  N 1;  a , i  j  N ; i  1, j  2;   Zqi j1 , інше.

Z1 

h

Substitute (22) into (23) and obtain a system of equations that are solved relative to the function T:

,  2     2   h U p , 0,..., 0,   2   h U p T  B  CATnk  12   c

CTnk   C  U nk  1  2  0 .

Finding from (24) the set of solutions - T1k 12 ,T2k 12 ,...,TNk112 ,TNk 12 , (k  0,1,..., K 1) , looking for from (22) the set of solutions U1k 12 ,U2k 12 ,...,U Nk112 ,U Nk 12 , (k  0,1,..., K 1)

For solutions on the whole interval  , use concealer, which is implemented in an explicit finite difference scheme: (24) (25) (26) с

T k 1  T k

n n   2       h n  q jTnkj121   0r j0 n k  12  a j0 q jU n j1 h n  q jTnkj121 . j0 Thus, from (26) find the great number of decision - U k 1 : k  0, K 1; n  1, N , n and from (25) will get the great number of decision - T k 1 : k  0, K 1; n  1, N . n conditions (14) will find the value of functions Tn0 ,U n0 .

2. For a condition 0  k  K carry out next operations: 2.1. REALIZATION OF THE PREDICTOR METHOD. 2.1.1. For n  1,..., N ; carry out cycles for k  k  1 ; 2 2.1.2. Determined from the matrix equation (24) the value Tnk 12 ; 2.1.3. From the matrix equation (22) the value U k 12 ; n 2.2. REALIZATION OF METHOD CORRECTOR. 2.2.1. For k  k 1; n  2,..., N 1 in cycles we find the value T k 1 and U k 1 of n n equations (25) and (26); 2.2.2. Wanted values T1k1 , TNk1 , U k1 ; U k1 we find from the boundary conditions 1 N (15), (16).

3. Increasing the time step to k  k 1 and for the condition 0  k  K carry out realization of sub-items 2.1 - 2.2, that is, the method of predictor-corrector. In the opposite case, that is, if the condition is not fulfilled 0  k  K , complete the implementation of the calculations.

The conditions of stability. To determine the stability conditions of the obtained difference equations of the connected heat-and-mass transfer, the method of conditional assignment of some known functions of the system is used, according to which the following relation is found: t  С1   h1  С2   1 С3 h2     2      2    (27) where С1  1, a1; C2   2 , a2 ; C3  c   0 r , 1   1 .

Supposing that fractal parameters  ,  take integer values, an analysis and comparison have been made, according to which the obtained stability condition (21) coincides with the condition of stability for classical equations of thermal conductivity. 5

Software for implementation of mathematical models For numerical solving of discretized model it was created a programming software, by using Python programming language. When the program loads you can see interface Fig. 1, where you can set your data which is needed by model, and run by pressing calculate button. It automatically sets all data inside the code, calculate the results of temperature and humidity which are unknown and draw graphics. After that, you can take cutting of the graphics over t or x variable and see 2D cutted graphics from 3D.

Software development for implementation of mathematical models and analysis of simulation results The coefficients included in the mathematical models for non-isothermal humidity transfer in materials with a fractal structure are considered to be temperature- and humidity-dependent. They include factors such as thermal conductivity, heat transfer, water exchange, conductivity of humidity, thermo-gravity coefficient, equilibrium humidity and density.

For the numerical experiment, pine wood has been chosen with the following values of the physical parameters of the material and the parameters of the environment.  basic density ρ0 = 530 kg/m3  the initial value of humidity content u0 = 0.4 kg/kg  initial temperature T0 = 20 oC  temperature of the environment Tc = 70 0C  relative humidity φ = 60%

We will show the dynamics of temperature and humidity of the numerical experiment of a one-dimensional mathematical model of heat and mass transfer processes taking into account the fractal structure (Fig. 2-3) for a biophysical material with a conditional density ρ0 = 360 kg/m3 for boundary conditions of the third kind and heterogeneous initial conditions, namely ( T 0, x  4  x2  20x  4.7 ).

The temperature of environment tc changes in relation to time on such law tc  80 1  exp  t   4 . The fractal parameters of the mathematical model were chosen as follows:   0.3,  1.9,  0,1. Equilibrium moisture content, depending on time t, is described by the formula: U p  0.0003t3  0.0063t 2  0.0193t  0.6 , the initial moisture content is equal to the moisture content at the border of satiation of cellular walls of material. Taking into account of fractal structure gives an opportunity to extend the great number of realization of mathematical model of unisothermal moisture transfer by the choice of fractional indexes, and also to take into account the evolution of temperature and humidity during all time of flowing of physical processes. Synthesized mathematical models of transferring heat and humidity in environments with fractal structure, taking into account non-localities in time and spatial correlation.

Finite-difference approximations of the system of differential equations of fractional order with boundary conditions of the third kind are obtained.

The explicit and implicit difference schemes for the realization of a mathematical model of heat and humidity transferring with derivatives of fractional order are developed.

The algorithmic aspects of their implementation are based on the use of the predictor-corrector method.

The synthesized one-dimensional mathematical models of non-isothermal humidity transfer in media with fractal structure taking into account the effects of memory and spatial correlation. Builted explicit and implicit difference schemes for equations related heat-mass transfer with boundary conditions of the third kind. These algorithmic aspects for the realization of the obtained difference equations was using method predictor-corrector. 11. Povstenko. Y.Z.: Two-dimensional axisymmetric stresses exerted by instantaneous pulses and sources of diffusion in an infinite space in a case of time-fractional diffusion equation.

International Journal of Solids and Structures, 44, pp. 2324-2348. (2007) 12. Pyanylo Ya. D. Lopuh N. B.: Numerical analysis of models with fractional derivatives for gas filtration in porous media. Coupled Syst. Multiscale Dyn. 2, 15-19. (2014) 13. Sokolovskyy Ya., Levkovych M., Mokrytska O., Atamanyuk V.: Mathematical modeling of anisotropic visco-elastic environments with memory based on integro-differentiation of fractional order. 14th International Conference on Advanced Trends in Radioelectronics, Telecommunications and Computer Engineering (TCSET), pp. 816-820. Lviv (2018) 14. Sokolowskyi, Ya., Shymanskyi, V., Levkovych, M., Yarkun, V.: Mathematical and software providing of research of deformation and relaxation processes in environments with fractal structure, Proceedings of the 12th International Scientific and Technical Conference on Computer Sciences and Information Technologies (CSIT), pp. 24-27. Lviv (2017) 15. Sokolovskyy, Ya., Shymanskyi V., Levkovych M.: Mathematical Modeling of Non-Isothermal Moisture Transfer and Visco-Elastic Deformation in the Materials with Fractal Structure, Computer Sciences and Information Technologies - Proceedings of the 11th International Scientific and Technical Conference (CSIT), p. 91-95. Lviv (2016) 16. Sokolovskyy, Ya., Shymanskyi, V., Levkovych, M., Yarkun, V.: Mathematical Modeling of Heat and Moisture Transfer and Reological behavior in Materials with Fractal Structure using the parallelization of Predictor-Corrector Numerical Method, 1-st International Conference Data Stream Mining Processing (DSMP), pp.108-111. Lviv (2016) 17. Shymanskyi, V., Protsyk Yu.: Simulation of the Heat Conduction Process in the ClayditeBlock Construction with Taking Into Account the Fractal Structure of the Material, Computer Sciences and Information Technologies - Proceedings of the 13th International Scientific and Technical Conference (CSIT), p. 151-154. Lviv (2018) 18. Tong Liu.: Creep of wood under a large span of loads in constant and varying environments.

Pt.1: Experimental observations and analysis, Holz als Rohund Werkstoff 51, pp. 400-405. (1993) 19. Uchajkin V.: Method of fractional derivatives Ulyanovsk: Publishing house "Artishok".

(2008) 20. Vasilyev V.V., Simak L.A.: Fractional calculus and approximation methods in the modeling of dynamic systems. Scientific publication Kiev, National Academy of Sciences of Ukraine, p. 256. (2008) 21. Zavada P.: Operator of fractional derivative in the complex plane, Communications in Mathematical Physics, Vol.192, pp. 261-285. (1998)

References

1. Datsko B.Y. , Gafiychuk V.V. : Different types of instabilities and complex dynamics in reaction-diffusion systems with fractional derivatives . Computational and Nonlinear Dynamics . ( 2012 ) DOI No: CND-09-1119

2. Diethelm K. , Ford N. , Freed A .: predictor-corrector approach for the numerical solution of fractional differential equations . Nonlinear Dynamics , 29 . pp. 3 - 22 . ( 2002 ).

3. Gafiychuk V. , Datsko B. : Mathematical modeling of different types of instabilities in time fractional reaction-diffusion systems , Computers and Mathematics with Applications , 59 . pp. 1101 - 1107 ( 2010 )

4. Kirichenko L.O. , Chalaia L.A. : An integrated approach to the study of fractal time series , International Journal "Information Technologies & Knowledge" , Vol. 8 , No 1, pp. 22 - 28 . ( 2014 )

5. Kostrobij P. P. , Grygorchak I. I. : Mathematical modeling of subdiffusion impedance in multilayer nanostructures , Math. Model. Comp . Vol. 2 , No. 2 , pp. 154 - 159 . ( 2015 )

6. Mainardi F. , Rusev I. , Kiryakova V.: Considerations on fractional calculus: interpretations and applications . Transform Methods and Special Function, Varna' 96 , Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, ISBN 954-8986-05-1 . p. 594 - 597 . Varna ( 1998 )

7. Nigmatullin R. R. : Fractional integral and its physical interpretation . Theoretical and Mathematical Physics , Vol. 90 , No. 3 , pp. 354 - 368 . ( 1992 )

8. Podlubny I. : Fractional Differential Equations , of Mathematics in Science and Engineering , Academic Press, San Diego, Calif., USA. ( 1999 )

9. Podlubny I. : Geometric and physical interpretation of fractional integration and fractional differentiation . Fractional Calculus and Applied Analysis , Vol. 5 , p. 367 - 396 . ( 2002 )

10. Povstenko Y.Z. : Fractional heat conduction equation and associated thermal stress . Journal of Thermal Stresses , 28 , pp. 83 - 102 . ( 2005 )