Variational Method for Solving the Viscoelastic Deformation
Problem in Biomaterials with Fractal Structure
Volodymyr Shymanskyia and Yaroslav Sokolovskyyb
a
Department of Artificial Intelligence, Lviv Polytechnic National University, Lviv, Ukraine
b
Department of Computer-Aided Design Systems, Lviv Polytechnic National University, Lviv, Ukraine
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.
Keywords 1
Fractal structure, viscoelasticity, variational formulation, finite element method, Galerkin’s
method.
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].
Information Technology and Implementation (IT&I-2021), December 01–03, 2021, Kyiv, Ukraine
EMAIL: vshymanskiy@gmail.com (V. Shymanskyi); sokolowskyyyar@yahoo.com (Ya. Sokolovskyy)
ORCID: 0000-0002-7100-3263 (V. Shymanskyi); 0000-0003-4866-2575 (Ya. Sokolovskyy)
©️ 2022 Copyright for this paper by its authors.
Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
CEUR Workshop Proceedings (CEUR-WS.org)
125
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 integro-
differentiation 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]:
1 f ( , y, z ) d
x
1
Dx f , (1)
(1 ) a 1
( x )
1 f ( , y, z ) d
x
1
I x f
( ) a 1
( x )
, (2)
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]:
Dx Dx f Dx f , (3)
I I f I
x
x
x f. (4)
Dx f g Dx f Dx g ,
(5)
I x f g I x f I x g.
(6)
Provided that f ( x, y, z) xa 0 the following property is valid [7, 15, 16]:
Dx I x f I x Dx f f . (7)
The parts fractional integration can be written as [7, 15, 16]
f I g dx g I f dx.
b b
x x (8)
a a
It is possible to generalize the derivative of the whole order to fractional derivatives (Riemann-
Liouville derivatives and Caputo derivatives). Fractional derivatives appear in new physical, technical
and chemical problems arising in research activities.
126
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 in the corresponding directions. And also surface forces FV X V , YV , ZV
T 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 - stress, ε , , , , , -
T
x y z xy xz yz
T
deformation and displacement u u , , , which are satisfying the equilibrium equation in
T
elementary volume [9]:
D xj ij Fi 0 (9)
and the equilibrium conditions on the surface [17, 18, 9]:
FV ij cos(n, x j ) (10)
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
2
1
D x j 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 Rt 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]:
σt εt I t R t rel εt , (12)
εt σ t I t Πt rel σt .
(13)
Using the following equations makes it possible to express stress due to deformation:
x x I t R11 t rel x y
I t R12 t rel y z I t R13 t rel z ,
(14)
y x I t R21 t rel x y
I t R22 t rel y z I t R23 t rel z ,
(15)
z x I t R31 t rel x y
I t R32 t rel y z I t R33 t rel z ,
(16)
I R t ,
xy xy t 44 rel xy (17)
xz xz I t R55 t rel xz ,
(18)
yz yz It R66 t rel yz .
(19)
127
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]:
σ ε dV F u dV F u dS.
V
ij ij
V
i i
S
Vi i (20)
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, , Au, , u, , dV u, , dS , (21)
V S
u, , Xu Y Z , (22)
u,, XV u 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]:
Au, , div u 2
D u * I R t D u
x
t 11 rel
x
R12 t rel D y R13 t rel D z
D y * I t R21 t rel D x u
R22 t rel D y R23 t rel D z
D z * I t R31 t rel D x u
R32 t rel D y R33 t rel D z (24)
D u D D u D
y
x
2
y
x
2
D D
z
y
2
D u D * I R t D u D
y
x t
44 rel
y
x
D u D * I R t D u D
z
x
t 55 rel
z
x
D D * I R t D D
z
y
t 66 rel
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
t
y
I R12 t rel D I R13 t rel D z cosn, x
t
D u D I R t D u D cosn, y
y
x
t 44 rel
y
x
(25)
D u D I R t D u D cosn, z
z
x t
55 rel
z
x
128
YV div u I t R21 t rel D x u
t
I R22 t rel D I R23 t rel D z cosn, y
y
t
D u D I R t D u D cosn, x
y
x
t 44 rel
y
x
(26)
D D I R t D D cosn, z
z
y
t 55 rel
z
y
Z div u I R t D u
V
t 31 rel
x
I R t D I R t D cosn, z
t
32 rel
y
t 33 rel
z
D u D I R t D u D cosn, x
z
x
t 44 rel
z
x
(27)
D D I R t D D cosn, y
z
y
t 55 rel
z
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:
Au, d X V u YV d (28)
Here A u, will look like:
Au, div u 2
D u * I R t D u R t D
x
t 11 rel
x 12 rel
y
D * I R t D u R t D
y
t 21 rel
x 22 rel
y (29)
D u D D u D
y
x
2
y
x
2
D u D * I R t D u D
y
x
t 33 rel
y
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 cosn, x
(30)
D u D I R t D u D cosn, y
y
x
t 33 rel
y
x
Y div u I R t D u
V
t 21 rel
x
I R t D cosn, y
t 22 rel
y (31)
D u D I R t D u D cosn, x
y
x
t 33 rel
y
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, 8000T .
Considering that the force FV acts in the direction opposite to the axis conditions (30)-(31) takes
the form:
129
D u D I R t D u D F
y
x
t 33 rel
y
x V 1 (32)
div u I R t D u I R t D F
t 21 rel
x t
22 rel
y V 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 u 0 x, y a i u i x, y (34)
i 1
n
x, y 0 x, y bii x, y (35)
i 1
The u 0 , 0 must satisfy the boundary conditions (32)-(33). u i , i are linearly independent basis
functions. We use Robotnov's fractional-exponential operators [3] to approximate the functions u 0 , 0 .
u0 x, y a01 x a02 a x a03 y a04 b y (36)
0 x, y b01 x b02 a x b03 y b04 b y (37)
x k ( 1)
x 1k 1
k 1
(38)
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 min (39) a0 j , b0 j
f a01 , a02 , a03 , a04 , b01 , b02 , b03 , b04
D u D
y 0
x 0
I R t D u D
F 2
t 33 rel y 0 x 0 x , y 2 V 1 (40)
div u I R t D u
F
0 t 21 rel x 0
I R t D
2
t 22 rel y 0 x , y 2 V 2
and satisfy the condition
u0 x, y x , y 0 (41)
1
0 x, y x , y 0 (42)
1
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 above-
described 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.
130
Figure 1: Distribution of stress component y depending on space with taking into account the
fractal structure and without
Figure 2: Distribution of stress component x depending on space with taking into account the fractal
structure and without
Let’s 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 of the axis y , i.e. F 0, FV 0, 8000 . The dynamics of
T
the stress components depending on the spatial coordinate y are shown in Figure 3 and Figure 4.
131
Figure 3: Distribution of stress component y depending on space with taking into account the fractal
structure and without
Figure 4: Distribution of stress component x depending on space with taking into account the fractal
structure and without
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%.
132
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
[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
133
[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 advection-
diffusion-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
134