CC BY
12
УДК 531.31:519.6
Discrete a Non-linear Hamiltonian Dynamics Models of Hyper Elastic Deformable Media
Vladimir A. Petushkov*
Blagonravov' Mechanical Engeneering Research Institute RAS, S. Kharitonievsky lane, 4, Moscow, 101990
Russia
Received 27.10.2012, received in revised form 03.11.2012, accepted 03.12.2012 A method of mathematical modeling of the dynamics of three-dimensional nonlinear deformable hypere-lastic media is developed in this paper. The method is based on the Hamiltonian description of discrete classical mechanics and on symplectic integration method for the solution at each instant of time. Comparative results of the solution of a model problem are presented. The results of solution of the topical problem of dynamic behavior of an artificial aortic artery are also presented.
Keywords: deformable media, finite deformation, Hamiltonian description, symplectic integrator, point approximation, mathematical simulation.
Introduction
Discrete classical mechanics with the Hamiltonian description of the motion of a system of bodies, a set of particles or local vortexes has been successfully used for long time in celestial mechanics, mechanics of machines and mechanisms, plasma physics, molecular dynamics and hydrodynamics, etc. [1-3]. The latest results of the use of corpuscular approach in the study of polycrystalline bodies confirm that the methods of discrete classical mechanics can be also applied to micromechanics [4]. The advantage of such an approach is that it allows one to describe adequately the objects of complex geometry and local features of solutions. It also allows the use of symplectic integrators to obtain stable in time solutions of Hamiltonian system. This approach also greatly simplifies the equations integration. For simulation of deformable media, this approach is rarely used. Usually finite and boundary element methods, finite difference methods are used in their traditional form. Meshfree methods with the point representation of computational domain and with a local approximation of the solution in weak (variational) form are also used [5,6]. In this work a new and rather universal method for modeling the dynamics of continuous media in which the dissipative processes are of minor importance is presented. The method is based on discrete Hamiltonian mechanics. It combines the advantages of meshfree Galerkin' method [5] and the discrete representation of the gradients of field functions which follows from the special form of the method of least squares [7].
A domain occupied by a deformable medium is approximated by a randomly distributed set of points. Approximate solutions are constructed by means of the local groups of material points. Every point has volume and weight and is not topologically connected with other points. Therefore, there is no need to introduce a grid. This approach is especially useful for the study of nonlinear dynamics of biological objects using their real anatomical geometric models (images)
*[email protected] © Siberian Federal University. All rights reserved
obtained with the use of tomography [8]. The set of spatial pixels that forms the picture of an object image can be directly used as a point (discrete) model of the object. The effectiveness of the proposed approach is confirmed by the solution of some model problems. The approach allows the construction of the dynamic models of aorta and heart with synthetic or native aortic valves. These models are useful in understanding the complex dynamics of such objects, especially in the presence of abnormalities (pathologies), and in developing the future generations of artificial heart valves. Some of the results of solution of the problem of artificial aortic artery dynamics are presented in this paper.
1. Formulation of the problem. Description of the mathematical model
Let a non-linear deformable elastic body of volume Vand arbitrary shape occupy at the initial moment of time t0 an open domain Q c R3 bounded by the surface S = Sa U Su, Sa n Su = 0, where Sa and Su are parts of the surface S where forces and displacements are specified. The initial (primary) configuration of the domain is denoted by k0(Q) and it is related to the Cartesian coordinate system Xi. With each point of the deformable medium we associate a convective curvilinear coordinate system n k in such a way that Xi = Xi(n k, t), and introduce the Rieman-nian metric aij(n k,t) = aiaj, where ai are tangent to the coordinate curve n k vectors. In what follows the conventional in tensor calculus usage of indexes and operations is adopted.
The movement (deformation) of the medium with respect to k0(Q) at any arbitrary point in time t > t0 is defined by the following C1 mapping:
xi = ¿(Xi, t),
Xi e k0(Q), k0 : Q ^ R3, (1)
t e Dt = (t0, T),
for any actual configuration kt, t > t0, ai = (X,t), aij (n k ,t) = aiá¿.
where x1 = xl(rjk, t) are the spatial coordinates of the point in the deformed medium. Therefore,
d^
d'rf
The time change of the initial geometry of the body is determined by displacement vector u1 = ui(nk ,t)
ui = xi - Xi (2)
movement of the deformable medium is determined by velocity field
Vi
= x i(X k ,t) = iii(X k ,t). (3)
Deformation measure of the medium is the gradient F, which may be represented as follows
F = (X,t) = aiai, FT = ai j = F, J = det(F) > 0. (4)
dX
Therefore, tensor of strains of the medium can be written in the following form
2E = (aij - aij)ai ® aj = (FTF - a), (5)
where a = aijai <g> aj is the metric tensor. Using (4) and (2) we get the Green deformation tensor
2 Eij = (ViUj + Vj ui + Vi ur Vj ur), (5a)
here
VjMj = dj — rpjiUp = u < V (6)
i da da ' da ■ ■ \
is the displacement gradient and j = 1/2ap s i sj + —j — —1 is the Christoffel's symbol.
y ox oxj ox J
The objective finite deformation rate tensor is defined as
Eij = FkiDkiFij, (7)
1 dv
where D is the stretch rate tensor with components Dij = 2(Lj + Lji) and Lij = -X = FikF—1
is the gradient of strain rate defined by (3).
Local stresses result when deformation of the medium occurs. They are defined by the second symmetric Piola-Kirchhoff tensor
Tij = JF-k1aklFjl 1 = Frk1rklFjl 1, (8)
where aij is the true Cauchy stress, rkl = Jakl is the Kirchhoff stress tensor. It should be noted that tensor T is referred to the initial configuration k0 of the medium but tensors a and t are referred to the current configuration kl.
Material derivative of the true stress is not objective or invariant with respect to the displacement of the medium as a rigid body. Therefore, the stress rate is generally understood as the Jaumann-Noll rate [9]
0
Tjk = Tjk — Tjr Wrk — Tkr Wrj, (9)
where Tjk, is the material derivative of stresses Wij = 2(Lij — j) are the components of the vorticity tensor.
The equations describing the non-stationary nonlinear deformation of the medium follow from the laws of conservation of mass, momentum and energy. These equations may be represented in Q, x Dt in the form:
p(xi,t)J = po(Xi),
(TklFjl ),k +Pobj = PoVj, Tkl = Tlk, (10)
poev9 = —Vq + T : E + poh, where p is the medium density b is the vector of mass forces 9 is the temperature cv is the specific heat, h is the density of the internal heat sources, q is the vector of heat flux. Equations (10) must be supplemented with an equation of state and with appropriate boundary conditions.
In the case of nonlinear hyperelastic behavior of the medium with large displacements and finite strains, thermal effects can be ignored. The existence of a potential energy function e that depends only on the strain gradient (4) or its corresponding measure (for example (5)) is assumed, i.e. e = e(F) or e = e(E).
Stress (8) and the equation of energy conservation (103) can be expressed in term of this function:
T = pof
poe = T : E. (IO3, a)
Constitutive relation for an isotropic medium in the form of a generalized Hooke's law follows from the expression for the strain energy Win = p0e
T = d-WEr = C : E, (11)
2
where C = 2pI + (K - 3p)I < I is the elasticity tensor, I is the unit tensor, p and K are the isotropic elastic constants. For elastic polymers relation (11) is usually used in the neo-Hookean form of Mooney-Rivlin type [8]
T = -pI + pB, (12)
where p is the local pressure of the medium, B = F • FT is the Finger deformation tensor. The boundary conditions for equations (10-(12) are written in the classical form
ui(xj ,t) = u* (xj ,t), xj e Su, t e Dt,
fi(xj ,t) = Tik (xj ,t) • nk, nk , xj e Sa, t e Dt, (13)
ui(xj,0) = ui(xj); vi(xj, 0) = vi(xj), xj e
2. Discrete Hamiltonian dynamics model of a deformable body
We replace the motion of the nonlinear deformable medium with the motion of a finite number of material particles. A domain occupied by the medium is presented as an arbitrary set of points. Let QN besomesetof N points xi that fill the domain QN = {xi, x2, ...,xi, ..., xN}, xi e
All these points are considered as material particles qi, i= 1, 2, ... N. For example, they have spherical shape of radius hi, that is qi := {y e R3 : \\xi - y\\R3 < hi}^. A set of all particles qi determines an open covering QN := {qi }N=1 of the domain and parameter h = min h
i=1, N
characterizes the spatial solvability of the model.
The size of parameter h or the cloud density of the non-uniformly scattered points is determined by the required accuracy of domain approximation and available computational resources. More particles per unit volume are needed in areas with the gross shape variation and/or considerable field variable gradient.
N
Every particle qi has volume vi and mass mi provided that vi ~ V and the sum of particle
i=i
masses is equal to total mass. Particle position in R3 at time t e Dt determined by the vector xi and linear momentum p-i = m^xi, where xi = dxi/dt is the particle velocity. Then each new position of the particle defined by the mapping (1) relative to the initial position characterizes the motion of the deformable medium and defines the discrete displacement vector (2) which is the main variable of the physical field.
To calculate derivatives of the functions (1), (2), or derivatives of deformation gradients with respect to discrete values obtained at each point xi e QN, we use the "moving" form of the least squares method [7]. In this case we consider not only the given point i but also the set of n neighboring points or the cloud c QN.
The distance between particle j and particle i in the initial configuration is r0j = (x0 - x0) and in the current configuration is rij = (xj - xi). Let us find the vector Fh that approximates at the point i the strain gradient F so that the residual
Ji = E \Fi4 - rijfSij (r) (14)
jeOn
tFor particles of other form, hi is the radius of circumsphere.
is minimal with a weight function Sj(r) defined on a compact support [r G |r| < r*j,
Sij (r) = S (| rj | ; r*) and Sj (r) = 0 for all |rj |
where r* is a suitable radius, Sij(r) = S(|rjI ; r*) and Sij(r) = 0 for all IrjI > r*, otherwise
Sij (r) = l - \rj \ /r\.
The selected weight function satisfies all the necessary conditions [7] that is the conditions of continuity, positivity: Sij (r) > 0 in the domain , compactness Sij (r) = 0 outside this region. The size of the compact support is normally chosen from the condition r* « (2 ^ 6)h.
dJi
Taking into account (4) and condition 0 = —f = 2 ^ (F^j-rij ) ® rjSij (r) we obtain
dFi je On
Fh = (VV)i = ^ rij ® r°jSij(r)G-1, (15)
where Gi = rij- ® rijSij (r), and from (5) we obtain the discrete Green deformation tensor.
jeOn
The kinetic and potential energy of the discrete system that moves in space
R3 are
N 2 N
K =E|m; =£ ^(f*), (16)
where ei is the density of strain energy of the particle qi. Strain energy depends on the discrete deformation gradient Fh. In this case elastic potential Win is the function of particle position x only.
In the absence of external forces the Hamiltonian is expressed in the following form (16)^:
H (x,p) = K (p) + Win(x), (17)
where H(x,p) is in general a piecewise smooth function. Hence we turn from the study of motion of the deformable medium in space R3 to the study of motion of discrete system in phase space R2N.
N
If domain Q in R2N has a canonical symplectic structure u = dxi A dpi then basic second
i=i
order equation (102) can be changed for the Hamiltonian first order system of ordinary differential equations of the form:
dH (x,p) dH (x,p) x = p =--(18)
Taking into account (15) we obtain the following equation of particle motion
d H
Pox = - ^ = £ (Tj Fj G-1r% + TiFiG7% )Sij (r). (19)
ij
This equation guaranties the energy conservation given by (103) and (103 , a). Symplectic integrators are usually used for numerical solution of Hamiltonian systems [2,3]. These integrators are conservative schemes, and they provide stable solutions in temporal domain Dt. We choose 1st order explicit symplectic scheme [10]:
xn+1 = xn + xnAt, xn+1 = xn + (dx/dt)n+i At. (20)
The stability condition for scheme (20): At = h/c, where c is the speed of sound in the deformable medium.
tin the presence of external forces their potential is added to H.
However, internal instability due to high-frequency non physical oscillations of the particles may occur in the system. The occurrence of instability depends on the number of points (size of the particles) and on the choice of the weight function and the radius of a compact support (14). That leads to a loss of accuracy of the solution. In this case, an additional potential of external pseudo dissipative force (see, for example, [10,11]) is included in the right-hand side of equation (19).
3. Numerical results and discussion
Example 1. Dynamics of a cantilever subjected to end instantaneous displacement. Dimensions of the beam and the loading characteristics are shown in Fig. 1. The cantilever material is aluminum with E = 2.0 • 105Pa and Poisson's ratio v = 0.33. The value of the initial displacement at t = 0 is u3= 0.1 mm. A various density of particles with h= (1^4)-10-3m. is used for the discrete approximation of the beam. Displacements of the beam point P with respect to the equilibrium position as a function of time are shown in Fig. 1.
Fig. 1
The conservatism of the numerical scheme is confirmed by the constancy of the total energy (see Fig. 2a, where (1) is the change in time of the kinetic energy, (2) is the change in time of the potential energy and (3) is the change in time of the total energy of the beam). Fig. 2b shows the movement of point P on the phase plane. These results confirm the stability of the symplectic integrator within the specified interval of time. The results are in good agreement with the FEM solution. With the same number of nodes the difference between the obtained results and the FEM results in this case does not exceed 2-6%.
Example 2. Pulsation of the aorta. Let us consider the problem of nonlinear dynamics of a shell of complex shape, namely, dynamics of aortic artery under the action of pulsating pressure (Fig. 3). The real 3D geometry of the artery (Fig. 3a) is obtained from computer tomography images with the help of pattern recognition methods [12,13]. As the picture of the cross-section A indicates the geometry data is the set of pixels with different color depth.
Fig. 2
Fig. 3
Number of computer model points out of the total number of pixels is determined mainly by the required accuracy of the solution and the available computer resources. The choice of a point is controlled by a given depth of the computer tomography image color tone. Each of these points represents a small material volume, piece of the aorta. The cloud of the points is the discrete representation of the aorta. Fig. 3b or Fig. 4a shows the discrete representation of some cross-section of the aorta.
The distance between particles characterizes the spatial resolution of computational domain. In this case the distance is 0.2 mm that corresponds to about 1000 particles in the cross-section shown in Fig. 4a.
o
Fig. 4
The experimental equation-of-state for hyper elastic material is taken in the form [14]:
Win
0.5m(Ii - 3) + 0.5K(log(FT))2,
where ^=1.17 106 Pa, k= 1.2 108 Pa, Ii is the first invariant of the deformation tensor (5a).
Let us compare the proposed computational model with the FEM model. Fig. 5 shows the time response of the aorta wall at point A (Fig. 4). The results demonstrate a good accuracy of the proposed discrete Hamiltonian model. The value of the wall displacement is somewhat overestimated (line 2 in Fig. 5) in comparison with the FEM results (line 1 in Fig. 5). This is the characteristic feature of meshfree methods.
12 10 8 6 4 2 0
■ 1Q4, MM
|2 .
1 \
I t,c I
0
10
Fig. 5
In order to choose a proper material for artificial aorta it is important to know the magnitude and distribution of stresses in the wall of the aortic artery under the influence of pulsating pressure. The distribution of dynamic stresses in the wall of the non-damaged aorta at some instant of time is shown in Fig. 4b.
Conclusions
3-D boundary value problem of nonlinear deformation of hyperelastic media of arbitrary shape under large displacements and finite strains is formulated. Meshfree particle computational model
for computer vision problems is presented. It is based on the point representation of a media and does not require mesh to represent domain occupied by the media. This approach makes it easier to handle large deformation and domain discontinuities with the desired numerical accuracy. The efficiency of the proposed computational model for study of the dynamic deformation of solids and objects of biomechanics is also demonstrated.
The work is supported by the Russian fund for fundamental research (grant 13-01-00977-a).
References
[1] V.I.Arnold, V.V.Kozlov, A.I.Neustadt, Mathematical aspects of classical and celestial mechanics. Dynamic Systems - 3. Results of science and tech., Ser. Modern Probl. Math. Fund. Direction. VINITI, 1985 (in Russian).
[2] B.Leimkuhler, S.Reich, Simulating Hamiltonian Dynamics, Cambridge, Cambridge Univ. Press, 2004.
[3] Integration Algorithms and Classical Mechanics, Ed. Marsden J.E., Patrick G.W. Shadwick W.F., AMS, Fields Institute Communications, 1996.
[4] A.M.Krivtsov, Deformation and fracture of solids with microstructure, M., PhysMatGis, 2007.
[5] I.Babuska, U.Banerjee, J.E.Osborn, Survey of meshless and generalized finite element methods: A unifed approach, Acta Numerica, 12(2003), 1-125.
[6] S.Li, W.K.Liu, Meshfree and particle methods and their applications, Appl Mech Rev., 55(2002), 1-34.
[7] P.Lancaster, K.Salkauskas, Surfaces generated by moving least squares methods, Mathematics of Computation, 37(1981), 141-158.
[8] L.A.Taber, Nonlinear theory of elasticity — applications in biomechanics, London, World Scientific, 2004.
[9] C. A.Truesdell, A First Course in Rational Continuum Mechanics, 2 ed., London, Acad. Press, 1991.
[10] J.C.Simo, N.Tarnow, K.K.Wong, Exact energy-momentum conserving algorithms and sym-plectic schemes for nonlinear dynamics, Comput. Methods Appl. Mech. Eng., 100(1992), 63-116.
[11] M.Kondo, S.Koshizuka, Y.Suzuki, Application of symplectic scheme to three-dimensional elastic analysis using MPS method, Trans. of the Japan Society of Mech. Eng., Part A, 72(2006), 425-431.
[12] G.Xu, M.Li, A.Gopinath, C.Bajaj, Computational Inversion of Electron Tomography Images Using L2-Gradient Flows, ICES REPORT 11-02, The Institute for Computational Engineering and Sciences, The University of Texas at Austin, 2011.
[13] A.I.Nadareyshvili, K.B.Chandran, J.Lu, Imaged-based explicit Dynamic Simulations in biomechanics, 10th International Conference on Dynamical Systems. Theory and Application, Lodz, Poland, 2009.
[14] G.A.Holzapfel, T.C.Gasser, R.W.Ogden, A new constitutive framework for arterial wall mechanics and a comparative study of material models, J. Elasticity, 61(2001),1-48.
Дискретные модели нелинейной гамильтоновой динамики гиперупругих деформируемых сред
Владимир А. Петушков
Разработан метод математического моделирования динамики трехмерных нелинейно деформируемых гиперупругих сред, основанный на Гамильтоновом описании дискретной классической механики и симплектическом методе интегрирования решения на временном слое. Представлены сравнительные результаты решения модельной задачи и решение актуальной задачи биомеханики о динамике аортальной артерии.
Ключевые слова: деформируемые среды, конечные деформации, точечная аппроксимация, Гамиль-тоново описание, симплектический интегратор, математическое моделирование.
