Variational Method for Solving the Viscoelastic Deformation Problem in Biomaterials with Fractal Structure

Abstract

The mathematical model of the stress-strain state for biomaterials with a fractal structure is considered. The dependences between the stresses, strains and displacements components in the biomaterial with taking into account the existing effects of self-organization, deterministic chaos and spatial nonlocality are obtained. The apparatus of integro-differentiation of the fractional order to take into account the above properties of the material during construction of mathematical models was used. A method for obtaining an approximate solution of the considered problem has been developed which uses a constructed variational formulation of the problem of the stress-strain state of biomaterials with taking into account their fractal structure. The distribution of stresses, strains and displacements components in the material for partial cases of the problem is analyzed. Fractal structure, viscoelasticity, variational formulation, finite element method, Galerkin's

1. Introduction

The study of the stress-strain state of bodies in the framework of the three-dimensional elasticity theory began a long time ago. The various methods for studying these problems are known now and many important results in this area have been obtained. The method for solving plane boundary problems of the elasticity theory for simply connected and multiply connected domains, which later became classical, proved effective. Another method for studying three-dimensional problems of the elasticity theory is the method of integral equations that serve as the basis for the development of algorithms for the numerical solution of the problem [ 1-3 ].

Using analytical methods for studying the stress-strain state of complex configuration bodies is associated with very significant mathematical difficulties [ 4 ]. Therefore, taking into account the expansion of the capabilities of computer technology, recently various numerical methods (finite elements, finite differences, variational-difference, etc.) have begun to be widely used [ 5-8 ].

At present, in medicine, the studying of the stress-strain state of biomaterials is an urgent scientific task, particularly in cardiovascular surgery. The person's properties of blood vessels change with age including mechanical (Young's modulus, Poisson's ratio). As a result, the walls of the vessels lose their elasticity. The blood flow leads to a qualitative redistribution of stresses inside the vessel walls. An experimental study of this problem is very difficult. Therefore, the numerical study of the influence of the vascular clamp on the stress-strain state of the wall of blood vessels is an urgent scientific problem. Such models can be used in medicine instead of laborious and difficult experimental studies in the future. A practically important result is to obtain the mechanical characteristics of the material of the arteries at various pressures [ 9, 2 ].

2022 Copyright for this paper by its authors.

Also to the material that is constantly exposed to stress is the bones of the skeleton. They are subjected to different stresses every day. A bone fracture may occur depending on the condition of the bone tissue and the nature of the stress. The possibilities for conducting direct experiments in this area are significantly limited. There is a need to use the capabilities of mathematical modeling. Bones are biological structures with complex structural organization and complex geometric shapes. The bone material is heterogeneous, with significantly pronounced anisotropy and nonlinear mechanical characteristics [ 10, 11 ].

A significant number of real processes do not fit into the concepts of continuum mechanics and require the involvement of ideas about the fractality of the environment in which these processes occur. Such processes, for example, are the rheological behavior of biomaterials. The appropriately modified Hooke's law was used to describe them, which requires the use of the mathematical apparatus of fractional integro-differential calculus. Yu.N. Rabotnov introduced a generalization of the rheological equation to describe the behavior of hereditary media in the case of fractional derivatives [ 12, 3 ].

Thus, since biomaterials are characterized by a complex structure, existing effects of spatial nonlocality and deterministic chaos it is advisable to use a mathematical apparatus of integrodifferentiation of fractional order to build mathematical models of their rheological behavior during loading which will allow to take into account these properties.

2. Intego-Differentiation of Fractional Order Tools

One of the problems with fractional derivatives is that there is no single definition. The numerical methods used to obtain an approximate solution of the problems described by equations with derivatives of fractional order are closely related to the type of derivative [ 13, 14 ].

Let us consider the fractional order integro-differentiation operators integral of the function f ( x, y, z) over the variable x in Caputo's understanding in more detail [ 15, 16, 7, 10 ]: Dx f 

I x f 

1 (1   ) 

a 1 ( )  a x   1 f ( , y, z)

  1 x 1  f ( , y, z)  1 

d (x  )  ,

d (x  )  ,  where      ,   N, 0    1, ( )   x 1ex dx - gamma function. 0

Relation (1) determines the differentiation operator of fractal order  in the Caputo sense, and (2) integration.

It also needs to note the following properties of these operators [ 7, 15, 16 ]:

Dx Dx f   D  f ,

x

I x I x f   I x  f .

Dx  f  g   Dx f  Dx g,

I x  f  g   I x f  I x g.

Provided that f (x, y, z) xa  0 the following property is valid [ 7, 15, 16 ]:

Dx I x f   I x Dx f   f .

The parts fractional integration can be written as [ 7, 15, 16 ] b b  f I x g dx   g I x f dx.

a a

It is possible to generalize the derivative of the whole order to fractional derivatives (RiemannLiouville derivatives and Caputo derivatives). Fractional derivatives appear in new physical, technical and chemical problems arising in research activities. (1) (2) (3) (4) (5) (6) (7) (8)

3. Production of a Problem

Let’s consider the viscoelastic deformation problem in biomaterial with taking into account its fractal structure. Suppose that a body that is in equilibrium is affected by mass forces F  X , Y , Z T in the corresponding directions. And also surface forces FV  X V ,YV , ZV T with corresponding projections on the axis x, y, z . Let’s find the components of the stress-strain state of the body, namely vectors σ   x , y , z , xy , xz , yz T - stress, ε   x , y , z , xy , xz , yz T deformation and displacement u  u, , T , which are satisfying the equilibrium equation in elementary volume [ 9 ]: and the equilibrium conditions on the surface [ 17, 18, 9 ]:

Dj ij  Fi  0

x

FV   ij cos(n, x j ) where n - outer normal to the surface of the body S .

The relationship between displacements and deformations will be written as follows [ 1, 18 ]:  ij  12 Dxj u i  Dxi u j . (11)

Currently, the development of the fractal concept, using the mathematical apparatus of fractional integro-differentiation has caused a tendency to revise the basic provisions of the mechanics of biomaterials. This helps to model complex structures with a reasonable degree of adequacy. Using these tools it is possible to take into account the complex nature of nonlinear phenomena: the memory effects described above and the correlations of spatial dependences. In this case, in the passage to the limit, it is naturally possible to obtain previously known solutions, which makes it possible to formulate nontrivial generalized laws [ 19-21, 3 ].

Thus, using the updated Boltzmann-Voltaire theory for describing hereditary viscoelasticity the relationship between stresses and strains can be described through integral equations involving creep Πt  rel  and relaxation Rt  rel  kernels. Thus, the use of the apparatus of integro-differentiation of fractional order makes it possible to take into account the fractal properties of the considered biomaterial [ 22, 12 ]:

Using the following equations makes it possible to express stress due to deformation: σt   εt   I t Rt  rel εt , εt   σ t   I t Πt  rel σt .  x   x  It R11t  rel  x    y 

It R12 t  rel  y   z  It R13t  rel  z ,  y   x  It R21t  rel  x    y 

It R22 t  rel  y   z  It R23t  rel  z ,  z   x  It R31t  rel  x   y 

It R32 t  rel  y   z  It R33t  rel  z ,  xy   xy  It R44t  rel  xy ,

  xz   xz  It R55t  rel  xz ,  yz   yz  It R66t  rel  yz .

 (9) (10) (12) (13) (14) (15) (16) (17) (18) (19)

Thus, we obtain the formulation of the viscoelastic deformation problem in biomaterials with a fractal structure using the above approach.

4. Variation Formation of the Problem

Today, many finite-difference numerical methods allow to obtain an approximate solution of the problem at the points of area discretization. However, they are conditionally convergent to the exact solution. Moreover, large computing resources to ensure proper accuracy are required.

In turn, the use of the variational formulation of the problem will make it possible to obtain a continuous approximate solution. This will avoid area discretization [ 5, 6, 11, 18 ].

All general theorems for small deformations are based on the equation of virtual works [ 23 ]:  σijεijdV   FiuidV   FV iuidS. (20)

V V S where V - body volume, S - body surface.

Thus, among all permissible displacements u, , that satisfy the boundary conditions provide a minimum of full potential energy functional  : u, ,    Au, ,   u, , dV  u, , dS,

V  u, ,   Xu  Y  Z ,

S  u, ,   XVu  YV  ZV. (23)

Taking into account the dependencies between stresses and strains components (12) and using relation (11) the energy function of the deformation components can be written as [ 24, 23 ]: Au, ,   div u2  D u * I t R11t  rel D u 

x x D * I t R21t  rel D u  y x

R12 t  rel D  R13 t  rel D 

y z R22 t  rel D  R23 t  rel D 

y z D * I t R31t  rel D u 

z x D u  D 2  D u  D 2  y x y x

R32 t  rel D  R33 t  rel D 

y z D  D 2 

z y D u  D * I t R44 t  rel D u  D 

y x y x D u  D * I t R55 t  rel D u  D 

z x z x D  D * I t R66 t  rel D  D 

z y z y The boundary conditions of equilibrium on the surface can be written as:

X V  div u  I t R11t  rel Dx u

I t R12 t  rel D  I t R13 t  rel D cosn, x 

y z D u  D  I t R44 t  rel D u  D cosn, y 

y x y x D u  D  I t R55 t  rel D u  D cosn, z  z x z x (21) (22) (24) (25) Dz  Dy  I t R55 t  rel Dz  Dy cosn, y 

Thus, substituting relation (24) into the functional (21), we obtain a variation formulation of the viscoelastic deformation problem of biomaterials with taking into account their fractal structure.

5. Obtained Results

5.1. Method of Constructing an Approximate Solution

Using the obtained variational formulation of the constructed mathematical model will be studied the distribution of stresses, strains and displacements components in biomaterial with a fractal structure. Consider the cross-section of a bone with such spatial dimensions   x; y 0; a0; b 0;0.010;0.15..

Functionality (21) takes the form: Here Au,  will look like: 

    Au, d   X V u  YV d Au,   div u2  D u * I t R11 t  rel D u R12 t  rel D 

x x y Dy * I t R21 t  rel Dx u R22 t  rel Dy   D u  D 2  D u  D 2 

y x y x Dy u  Dx * I t R33 t  rel Dy u  Dx  The boundary conditions (25)-(27) will be written as follows:

X V  div u  I t R11 t  rel Dx u

I t R12 t  rel Dy cosn, x  YV  div u  I t R21 t  rel Dx u

I t R22 t  rel Dy cosn, y  Dy u  Dx  I t R33 t  rel Dy u  Dx cosn, y Dy u  Dx  I t R33 t  rel Dy u  Dx cosn, x

Let's fix the biomaterial on the boundary 1  x; y 0; ay  b. Then we will calculate stress in the sample after applying the force 8 KN / m 2 to the boundary 2  x; y 0; ay  0 in the direction opposite to the axis y, i.e. F  0, FV  0,  8000T .

Considering that the force FV acts in the direction opposite to the axis conditions (30)-(31) takes the form: (26) (27) (28) (29) (30) (31) div u  It R21t  rel Dx u It R22t  rel Dy   FV 2 (33)

Let’s find an approximate solution of the minimum of full potential energy functional (28) in the next form: n ux, y  u0 x, y   ai ui x, y (34) i1 n  x, y 0 x, y   bii x, y (35)

i1

The u0 , 0 must satisfy the boundary conditions (32)-(33). ui , i are linearly independent basis functions. We use Robotnov's fractional-exponential operators [ 3 ] to approximate the functions u0 , 0 .

u0 x, y  a01  x a02  a  x a03  y a04  b  y 0 x, y  b01  x  b02  a  x  b03  y  b04  b  y

 xk( 1)  x   (38)

k 1  1k 1

We formulate an optimization problem for finding unknown coefficients for the functions (36)-(37) that would satisfy the boundary conditions (32)-(33) in such way:

f a01, a02, a03, a04,b01,b02,b03,b04 a0 ,b0 j j min (39) f a01, a02 , a03 , a04 , b01, b02 , b03 , b04   Dy u0  Dx  0 

 I t R33 t  rel Dy u0  D    x 0 div u0  I t R21t  rel Dx u0  I t R22 t  rel D  y 0 x,y2  F  2

V 2 x,y2  F  2 

V 1 and satisfy the condition u0 x, y x,y1  0  0 x, y x,y1  0

Robotnov's piecewise fractional-exponential functions were chosen as the basis for constructing the approximate solution (34)-(35). Thus, the problem is to find such coefficients ai ,bi for which functions (34)-(35) give a minimum to the functional (28). 5.2.

Obtained Numerical Results

The femur bone of a 40 and 80 years old person was chosen as a model for numerical experiments. The cross-section was considered.

The distributions of the stresses, strains and displacements components of the described sample subjected to external loads are calculated. Developed software was used for this purpose. [ 25, 26 ].

Let's analyze the obtained values of stresses in the biomaterial which was subjected to the abovedescribed load. The dynamics of the stress components  y and  x depending on the spatial coordinate y are shown in Figure 1 and Figure 2 accordingly. We can conclude that the stresses at the boundary y  0 with taking into account the fractal structure and without differs by less than 2% but at the boundary y  0.15 at the   0.9834 this value is 3.2% and at   0.9368 this value is 9.7%. This is due to the existing effects of spatial nonlocality and the ability of the environment to remember the stress state. Figure 1 shows the transverse compressive stresses due to longitudinal tension. (32)

Let’s calculate stress in the sample after applying the force 8 KN / m2 to the boundary 2  x; y 0; ay  0 in the direction of the axis y, i.e. F  0, FV  0, 8000T . The dynamics of the stress components depending on the spatial coordinate y are shown in Figure 3 and Figure 4.

Analyzing those stresses dependences we can note that at the boundary y  0 obtained  y with taking into account the fractal structure and without is differ less than 1% but at the boundary y  0.15 at the   0.9834 this value is 1.7% and at   0.9368 this value is 6.9%.

Analyzing the graphical dependences in Figure 2 and Figure 4, we can see that the maximum by module values of the normal stress  x are an order of magnitude smaller than the maximum by module values of the normal stress  y . However, analyzing the influence of the fractal structure on the obtained results, we see that the maximum difference between the stress component  x with taking into account the fractal structure and without is equal to 18.7%.

The fractality degree parameter  was determined by approximating the experimental data on the creep of biomaterials using the Mittag-Leffler functions. In particular, it was determined that for the bones of a 40-year-old person the average   0.9834 , and for an 80-year-old -   0.9368 .

We conclude that the compressive stresses also accumulate due to the nonlocal structure of the medium.

The obtained results are also consistent with the results of other scientists conducting research in the direction of modeling the rheological unit of bones [ 27 ].

6. Conclusion

In this paper, we solve the actuality problem of constructing mathematical models of the viscoelastic deformation problem for biomaterials with taking into account their fractal structure. The main relations of rheological behavior of biomaterials under the action of external loads are obtained, which take into account the available effects of memory, self-organization and deterministic chaos in the material. Variational formulation of the problem will make it possible to obtain a continuous approximate solution and avoid diskretting the area. For partial cases, obtained results showed their adequacy and compliance with the simulated process. The dependence of the stress-strain state components on the fractality parameter is analyzed. The obtained results indicate that for materials with a higher degree of fractality the absolute values of the stress-strain state components increase.

7. References

References

[1] Carpinteri A. , Cornetti P. , Sapora A . Nonlocal elasticity: an approach based on fractional calculus . Meccanica 49 ( 11 ), 2014 , pp. 2551 - 2569 .

[2] Li L. , Hu Y . Buckling analysis of size-dependent nonlinear beams based on a non-local strain gradient theory . Int. J. Eng. Sci. 97 , 2015 , pp. 84 - 94 .

[3] Palaniappan D. A general solution of equations of equilibrium in linear elasticity . Applied Mathematical Modeling , 35 , Elservier , 2011 , pp. 5494 - 5499

[4] Sweilam N. H. , Khader M. M. , Al-Bar R. F . Numerical studies for a multi-order fractional differential equation . Phys. Lett, A 371 , 2007 , pp. 26 - 33

[5] Boffi D. Finite Element Methods and Applications. Springer Series in Computational Mathematics, 2013 , 575 p.

[6] Chyr I.A. , Shynkarenko H.A. Well-Posedness of the Green-Lindsay Variational Problem of Thermoelasticity . Journal of Mathematical Sciences , Vol. 226 , No 1, 2017 , pp. 11 - 27

[7] Hilfer R. , Luchko Y. , Tomovski Z. Operational method for the solution of fractional differential equations with generalized Riemann-Liouville fractional derivatives . In: Fract. Calc. Appl . Anal, № 12 , 2009 , pp. 299 - 318

[8] Sokolovskyy Y. , Boretska I. , Gayvas B. , Shymanskyi V. , Gregus M . Mathematical modeling of heat transfer in anisotropic biophysical materials, taking into account the phase transition boundary . CEUR Workshop Proceedings 2488 , 2019 , pp. 121 - 132

[9] Shymanskyi V. , Sokolovskyy Y. Variational Formulation of Viscoelastic Problem in Biomaterials with Fractal Structure . 3rd International Conference on Informatics & Data-Driven Medicine , IDDM-2020, Växjö, Sweden, November 19 - 21 , 2020 , pp. 360 - 369 .

[10] Cattani C. , etal. Fractional Dynamics. Emerging Science Publishers, Berlin, 2015

[11] Faraji Oskouie M. , Ansari R. , Rouhi H . Bending analysis of functionally graded nanobeams based on the fractional non-local continuum theory by the variational legendre spectral collocation method . Meccanica 53 ( 4 ), 2018 , pp. 1115 - 1130

[12] Rabotnov Yu.N. Equilibrium of an elastic medium with after-effect . Fractional Calculus and Applied Analysis , vol. 17 , no. 3 , 2014 , pp. 684 - 696 . DOI: 10 .2478/s13540-014-0193

[13] Tarasov , V.E. No nonlocality. No fractional derivative . Commun. Nonlinear Sci. Numer . Simul. 2018 , 62 , pp. 157 - 163 [doi:10.1016/j.cnsns. 2018 . 02 .019].

[14] Uchaikin V.V. Fractional Derivatives for Physicists and Engineers . Springer-Verlag: Berlin,Germany ( 2013 )

[15] Kilbas A. , Srivastava H. M. , Trujillo J. J. Theory and Applications of Fractional Differential Equations . Elsevier, 2006

[16] Podlubny I. Fractional Differential Equations . Academic Press, SanDiego, 1999

[17] Ogorodnikov E.N. , Radchenko V.P. , Ungarova L.G. Mathematical modeling of heredi-tary elastically deformable body on the basis of structural models and fractional integro-differentiation Riemann-Liouville apparatus . Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz. Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.] ,, vol. 20 , no. 1 , 2016 , pp. 167 - 194 . DOI: 10 .14498/vsgtu1456

[18] Savula Y. , Turchyn Y. Replacements in finite element method for the problem of advectiondiffusion-reaction . Journal of computing and applied mathematics. Series "Computational Mathematics" , Iss . 1 ( 130 ), 2019 , pp. 84 - 98

[19] Lewandowski R. , Chorążyczewski B . Identification of the parameters of the Kelvin Voigt and the Maxwell fractional models, used to modeling of viscoelastic dampers . Computers and Structures , vol. 88 , no. 1-2 , 2009 , pp. 1 - 17 . DOI: 10 .1016/j.compstruc. 2009 . 09 .001

[20] Mainardi F. Fractional calculus and waves in linear viscoelasticity. An introduction to mathematical models . Hackensack, NJ, World Scientific, 2010 , 347 p. DOI: 10.1142/9781848163300

[21] Sasso М. , Palmieri G. , Amodio D. Application of fractional derivative models in linear viscoelastic problems . Mech. Time-Depend Mater , Vol. 15 , 2011 , pp. 367 - 387

[22] Boltzman L. Zur Theorie der elastischen Nachwirkung / Wissenschaftliche Abhandlungen . vol. 2 / Cambridge Library Collection; ed. Friedrich Hasenӧhrl . Cambridge: Cambridge University Press, 2012 , pp. 318 - 320 . DOI: 10 .1017/CBO9781139381437.015

[23] Washizu K. Variational Methods in Elasticity & Plasticity, Pergamon Press, New York, 3rd edition ( 1982 )

[24] Di Pietro D. A. , Nicaise S. A locking-free discontinuous Galerkin method for linear elasticity in locally nearly incompressible heterogeneous media , J. Applied Numerical Mathematics , 63 , 2013 , pp. 105 - 116

[25] Shymanskyi V. , Sokolovskyy Ya ., Boretska I. , Sokolovskyy I. , Markelov O. , Storozhuk O . Application of FEM with Piecewise Mittag-Leffler Functions Basis for the Linear Elasticity Problem in Materials With Fractal Structure . IEEE XVIIth International Conference on the Perspective Technologies and Methods in MEMS Design (MEMSTECH 2021 ), Lviv, 2021 , pp. 16 - 19 . DOI: 10 .1109/memstech53091. 2021 . 9468027 .

[26] Shymanskyi V. , Sokolovskyy Ya . Variational Formulation Of The Stress-Strain Problem In Capillary-Porous Materials With Fractal Structure . XV International Scientific and Technical Conference on Computer Science and Information Technologies- CSIT2020 , Vol. 1 , 2020 , pp. 200 - 203 .

[27] Nekhozhin A. V. Dvukhslojnaya matematicheskaya model` shejki bedra cheloveka dlya issledovaniya napryazhyonnogo sostoyaniya pri armirovanii implantatami razlichnoj konstrukczii , Vestn. Sam. gos. tekhn. un-ta. Ser. Fiz.-mat. nauki , 3 ( 32 ), 2013 , pp. 129 - 135 . DOI: 10 .14498/vsgtu1254