A finite element method for non-linear hyperelasticity applied for the simulation of octopus ARM motions

147 IV International Conference on Computational Methods for Coupled Problems in Science and Engineering COUPLED PROBLEMS 2011 M. Papadrakakis, E. O˜nate and B. Schrefler (Eds) A FINITE ELEMENT METHOD FOR NON-LINEAR HYPERELASTICITY APPLIED FOR THE SIMULATION OF OCTOPUS ARM MOTIONS VASILEIOS VAVOURAKIS∗,†, ASIMINA KAZAKIDI†, DIMITRIOS P. TSAKIRIS†AND JOHN A. EKATERINARIS‡,§ ∗Department of Civil and Environmental Engineering University of Cyprus, Nicosia 1678, Cyprus e-mail: [email protected] †Institute of Computer Science, Foundation for Research and Technology-Hellas Heraklion, Crete 71110, Greece ‡Department of Mechanical and Aerospace Engineering University of Patras, Rio 26500, Greece §Institute of Applied and Computational Mathematics, Foundation for Research and Technology-Hellas Heraklion, Crete 71110, Greece Key words: FEM; sceletal muscles; muscular hydrostats Abstract. An implicit non-linear finite element (FE) numerical procedure for the simulation of biological muscular tissues is presented. The method has been developed for studying the motion of muscular hydrostats, such as squid and octopus arms and its general framework is applicable to other muscular tissues. The FE framework considered is suitable for the dynamic numerical simulations of three-dimensional non-linear nearly incompressible hyperelastic materials that undergo large displacements and deformations. Human and animal muscles, consisting of fibers and connective tissues, belong to this class of materials. The stress distribution inside the muscular FE model is considered as the superposition of stresses along the muscular fibers and the connective tissues. The stresses along the fibers are modeled as the sum of active and passive stresses, according to the muscular model of Van Leeuwen and Kier (1997) Philos. Trans. R. Soc. London, 352: 551-571. Passive stress distribution is an experimentally-defined function of fibers’ deformation; while active stress distribution is the product of an activation level time function, a force-stretch function and a force-stretch ratio function. The mechanical behavior of the surrounding tissues is determined adopting a Mooney-Rivlin constitutive model. The incompressibility criterion is met by enforcing large bulk modulus and by introducing modified deformation measures. Due to the non-linear nature of the problem, 1 148 V. Vavourakis, A. Kazakidi, D. P. Tsakiris and J. A. Ekaterinaris approximate determination of the Jacobian matrix is performed, in order to utilize the full Newton-Raphson iterative procedure within each time-step. In addition, time discretization is performed via the implicit Newmark method. We developed an open-source finite element code that is capable of simulating large deflection maneuvers of muscular hydrostats. The proposed methodology is validated by comparing the numerical results with existing measurements for the squid arm extension. The efficiency and robustness of the proposed numerical method is demonstrated through a series of octopus arm maneuvers, such as extension, compression and bending. 1 INTRODUCTION Muscular tissues are complex, non-linear, anisotropic, incompressible, viscoelastic materials undergoing large deformations. A wide class of muscular tissues deforms voluntarily and simulation of muscles behavior through the finite element method (FEM) has been the subject of many numerical investigations. Amongst the pioneering works in skeletal muscles’ computational mechanics through the FEM, Kojic et al. [1] proposed a FEM numerical algorithm for the determination of muscle response. In this work, Hill’s threeelement model was used as basis for the mechanical description of fibers and accounted for non-linear force-displacement relation and change of geometrical shape. In the same year, Martins et al. [2] developed a three-dimensional FE methodology for the simulation of skeletal muscles, where the constitutive relation adopted is a generalization of the model proposed by [3], being compatible with the passive and active behavior of skeletal muscles of [4]. Johannsson et al. [5] proposed a mixed total Lagrangian FE formulation for simulating muscular behavior, based on non-linear continuum mechanics, where contractile active and passive properties of skeletal muscles were considered and the stress distribution was assumed to be equal to the superposition of passive and active ones. Yucesoy et al. [6] considered the skeletal muscle in two domains: the intracellular and extracellular matrix domain, represented by two separate meshes linked elastically to account for the trans-sarcolemmal attachments of muscle fibers, which allows force transmission between these domains. Furthermore, Oomens et al. [7] proposed a FE approach where physiological reasoning and ”cross-bridge” kinetics via a two-state Huxley model was adopted, in order to examine the mechanical behavior of a tibialis anterior of a rat. In the recent work of Liang et al. [8], the governing equations for a muscle element, were based upon the approach of [9], and incorporated in a commercial general purpose explicit FE program. It was also, Martins et al. [10], in succession to [2], who introduced a multiplicative split of the fiber stretch into contractile and (series) elastic stretches, and they considered the simultaneous presence of the series elastic element, the dependence of the contractile stress on the strain rate and the activation level function. R¨ohrle and Pullan [11] presented a simulation framework of an anatomically realistic model of the human masseter muscles and associated bones, in order to investigate the dynamics of chewing. In the work of 2 149 V. Vavourakis, A. Kazakidi, D. P. Tsakiris and J. A. Ekaterinaris Tang et al. [12], Hill’s muscle theory coupled with fatigue was proposed to describe the mechanical behavior of skeletal muscles, where the force developed by a fatigued muscle was described by a muscle fatigue formula. Stojanovic et al. [13] proposed an extension of Hill’s three-component model of [1], in order to take into account for different fiber types. They presented a model consisting of different type sarcomeres coupled in parallel with the connective tissues, where each sarcomere was modeled by one non-linear elastic element connected in series with one non-linear contractile element. Most recently, Tang et al. [14] presented a three-dimensional FEM for skeletal muscles that was developed to simulate their mechanical behavior during lengthening or shortening. The constitutive relation of the muscle was determined by using a strain energy approach and active contraction behavior of the muscle fiber was modeled through the Hill’s three-element muscle concept. In addition, Lu et al. [15] developed a visco-hyperelastic model for skeletal muscle, where the constitutive relation was based on the definition of a Helmholtz free-energy function, while their model involves fourteen material parameters. In the present work, motivated by the OCTOPUS project (FP7-231608) that aims to the development of octopus-like robotic arms, an implicit non-linear finite element numerical procedure has been developed for the accurate simulation of biological muscular tissue dynamic motion. In the next section, the constitutive equations adopted and the FE framework developed in this work are presented in detail. The numerical results obtained with the implementation of these constitutive equations are initially validated, and presentation and discussion of results representative of octopus arm motions concludes the third section. The conclusions of this work are summarized in the last section. 2 FINITE ELEMENT METHODOLOGY Assume a three-dimensional non-linear, homogeneous elastic continuous medium and denote by V 0the volume and by S 0the surface in its unloaded state that undergoes large deformations, as shown in Fig. 1. The volume and bounding surface of the body in its current (deformed) state is V tand S t, respectively. The equilibrium equations for a solid, subject to finite deformations, are identical to those for small deformation analysis, i.e. ∂σij ∂x tj =ρ t¨ui,(1) where body forces have been neglected and ρ tis the material density. Prescribed displacements are assumed on S tUand prescribed tractions are assumed on the portion S tT of the bounding surface to obtain a well-posed boundary value problem. The stress equilibrium Eq. 1 is replaced by an equivalent principle of virtual work. Application of the variational theorem for finite elasticity in the current configuration yields [16] 3 150 V. Vavourakis, A. Kazakidi, D. P. Tsakiris and J. A. Ekaterinaris V 0,S 0 V t,S t Figure 1: Deformed (red) and undeformed (blue) octopus arm configurations. δΠ=V t δviρ t¨uidV +V t δLij σij dV −S tT δvi¯ tidS =0,(2) where ¯ tiis the prescribed traction vector on S tT,δvibeing an admissible velocity variation that satisfies the condition: δvi= 0 on S tU, and δLij =∂δvi /∂x tjthe virtual velocity gradient matrix. The integrals in the virtual work relation are considered with respect to the current configuration. However, proper transformations can be applied in order to transform the integrals into the reference configuration, since the undeformed elastic body configuration is known. These transformations are the following: ρ 0=Jρ t,dx ti= Fij dx 0j,dS t=JγdS 0and dV t=Jd V 0, where Fij =∂x ti /∂x 0jthe deformation gradient tensor (expressed in terms of the initial and current coordinates of a material point x 0i, x ti, respectively), Jthe deformation gradient determinant, Bij =FinFjn the left Cauchy- Green deformation tensor, and Jγ=JˆniB−1 ij ˆnj. In the present work, muscles are assumed as a composite material, comprising of muscle fibers and connective tissues. Therefore, one can safely assume that the stress distribution inside the muscle is the superposition of the stress distribution in the connective tissues and the fibers, respectively, i.e. σij =σ(ct) ij +σ(f) ij . Inserting the aforementioned transformations into Eq. 2, the virtual work variational balance equation can be expressed into the reference configuration δΠ=V 0 δviρ 0¨uiJ dV +V 0 δLij σij J dV −S 0T δvi¯ tiJγdS =0,(3) where ρ 0the initial elastic body material density. Finite element discretization of Eq. 3 using Lagrange polynomial shape functions δvi= αN(α)δv(α) i,∂δvi /∂x tj=α(∂N(α)/∂x tj)δv(α) iyields: 4 151 V. Vavourakis, A. Kazakidi, D. P. Tsakiris and J. A. Ekaterinaris V 0 N(α)ρ 0¨uiJ dV +V 0 ∂N(α) ∂x tj σij J dV −S 0T N(α)¯ tiJγdSδv(α) i=0⇔ V 0 N(α)ρ 0¨uiJ dV +V 0 ∂N(α) ∂x 0k F−1 kj σij J dV −S 0T N(α)¯ tiJγdSδv(α) i=0.(4) This discrete, non-linear virtual work equation is solved using Newton-Raphson iteration assuming a corrected updated solution u(α) i+∆u(α) i, where u(α) iis the solution at the end of the previous time increment or the previous iteration solution. After linearization of the integral Eq. 4 and some algebra it can be obtained V 0Cijkl ∂N(α) ∂x 0r F−1 rj ∂N(β) ∂x 0l F−1 lk −σij ∂N(α) ∂x 0m F−1 mk ∂N(β) ∂x 0n F−1 nj J dV ∆u(β) k =S 0T N(α)¯ tiJγdS −V 0 σij ∂N(α) ∂x 0p F−1 pj J dV , (5) where Cijkl is the tangent stiffness tensor given by Cijkl =∂σij ∂Fks Fls +σij δkl =C(ct) ijkl +C(f) ijkl ,(6) is the sum of the corresponding stiffnesses of the connective tissues and fibers contribution, respectively. Numerical evaluation of the domain and boundary integrals of Eq. 5 yields a system of linear equations. In the present work, the Newmark implicit method is utilized as time-integration procedure. Given the displacement, velocity and acceleration from the previous time-step, the displacement increment is calculated through Eq. 5 and the updated displacement, velocity and acceleration components are evaluated accordingly [17]. At this point, it is necessary to address the constitutive description of the muscular finite element model. As stated above, the muscle consists of connective tissues with biofluids and muscular fibers. Their material description is provided in detail in the following subsections. 2.1 Muscle fibers material description The fiber is the component that takes active role in the fiber and muscle contraction, resulting into muscle deformation. The fiber comprises of parallel bundles of myofibrils, which in turn are divided longitudinally by the Z-discs into sarcomeres. Sarcomeres are the basic contractile units of a muscle and they are also responsible for the passive properties of the muscle [18]. 5 152 V. Vavourakis, A. Kazakidi, D. P. Tsakiris and J. A. Ekaterinaris Assume a direction vector along the fiber, denoted with ˆniin the undeformed configuration, which is explicitly defined, and denote with ˆmi=1 /λ(Fij ˆnj) the direction vector in the deformed (current) configuration, where λ=ˆniFkiFkj ˆnjthe fiber stretch ratio. The nominal fiber strain εm 0is defined by the change of length divided by the reference length of the fiber, i.e. εm 0=λ−1. Therefore, the corresponding volume preserving fiber strain tensor can be written as ε(f) ij =1 2εm(3 ˆmiˆmj+δij). The corresponding Cauchy stress tensor has the form σ(f) ij =σmˆmiˆmj,(7) where the nominal axial stress σm 0in the fiber is defined in terms of the Cauchy true fiber stress as: σm=σm 0(εm 0+ 1), given that the volume of the fiber is preserved. According to the approach of [19], the fiber nominal axial stress in muscles is defined as the accumulation of passive σ(pass)and active axial stress. The latter part is considered as the product of the maximum isometric stress at fiber optimum length σ(max), a normalized active state function f(a), a force-length function f(l)and a force-velocity dependent function f(v) σm 0=σ(pass)+σ(max)f(a)f(l)f(v).(8) The activation state f(a)describes the activation pattern and is a function of time. The filamentary overlap function f(l)describes the dependence of active stress on the nominal fiber strain εm 0, and the force-velocity function f(v)is a rate-dependent function that relates the active muscle stress and the nominal fiber strain rate ˙εm 0. Inserting the above stress tensor relations in Eq. 6, one can analytically evaluate the tangent stiffness tensor as follows C(f) ijkl =λ∂σm 0 ∂εm 0 (1 + εm 0)+λσm 0−2σm 0(1 + εm 0)ˆmiˆmjˆmkˆml +σm 0(1 + εm 0)[δik ˆmjˆml+ˆmiδjk ˆml+ˆmiˆmjδkl].(9) 2.2 Connective tissues material description In the present work the connective tissues are described through a hyperelastic Mooney- Rivlin constitutive relation, where a modified stored energy function ¯ Wis used in order to evaluate stresses, provided by the generalized relation: ¯ W=c1(¯ I1−3) + c2(¯ I2−3) + K/2(J−1)2, where c1,c2,Kare material constants. The modified invariants of the left Cauchy-Green deformation tensor are introduced in ¯ W:¯ I1=I1J−2/3and ¯ I2=I2J−4/3 based on the modified deformation gradient tensor ¯ Fij =FijJ−1/3, due to the nearly or fully incompressible behavior of biological tissues the constraint J≈1 must be satisfied. The calculation of the Cauchy stress tensor can be obtained by differentiation of the stored energy density function with respect to deformation as 6 153 V. Vavourakis, A. Kazakidi, D. P. Tsakiris and J. A. Ekaterinaris σ(ct) ij =2 J5/3∂¯ W ∂¯ I1 +¯ I1 ∂¯ W ∂¯ I2Bij −2 J7/3 ∂¯ W ∂¯ I2 BimBmj −2 3J¯ I1 ∂¯ W ∂¯ I1 +2¯ I2 ∂¯ W ∂¯ I2δij +∂¯ W ∂J δij ,(10) with the Cauchy stress tensor consisting of the sum of a purely isochoric contribution and a purely volumetric one. Substituting the modified generalized Mooney-Rivlin constitutive material relation of ¯ Win Eq. 10, the Cauchy stress tensor for the connective tissues is σ(ct) ij = 2(c1+c2¯ I1)J−5/3Bij −2c2J−7/3BimBmj −2 /3J(c1¯ I1+2c2¯ I2)+K(1 −J)δij ,(11) where for the present analysis the material constant Krepresents the bulk modulus of elasticity, while c1is equal to the shear modulus and c2= 0. The corresponding tangent stiffness tensor of the connective tissues can be evaluated analytically through Eq. 6 and is provided below C(ct) ijkl =4 /9J(c1¯ I1+4c2¯ I2)+K(2J−1)δijδkl +2c2J−7/3(2BijBkl −BilBjk −BikBjl) −4 /3(c1+2c2¯ I1)J−5/3(Bijδkl +δijBkl)+8 /3c2J−7/3B2 ijδkl +δijB2 kl + 2(c1+c2¯ I1)J−5/3(δikBjl +Bilδjk)−2c2J−7/3δikB2 jl +B2 ilδjk,(12) where B2 ij =BikBkj. 3 NUMERICAL EXAMPLES In order to validate the finite element methodology, the squid arm extension during the strike to catch prey is simulated. The squid arm is modeled as a simplified cylindrical geometry, consisting of an active stalk and a passive club. The stalk consists of transverse muscles inside the arm and an outer layer consists of longitudinal muscles. Detailed description of the squid musculature is provided by Van Leeuwen and Kier [9], while the material properties of the muscle fibers and surrounding tissues are identical to the ones used by previous FE approaches [5, 8]. Due to symmetry, only one quadrant of the cylindrical arm is modeled. Symmetry boundary conditions are applied on the symmetry planes and the outer surface is considered traction-free. The discretized quadrant consists of 246 four-node hexahedral and 41 six-node triangular base prismatic elements, as seen in Fig. 3(a). The applied activation signal for the current simulation is a step function [9] having a 40 msec activation time for maximum activation level, while the total time duration of the simulation is 100 msec. In Fig. 3(a) a comparison of the squid arm-length growth in time is presented. The FE numerical results, obtained by the proposed methodology (diamonds), are compared 7 154 V. Vavourakis, A. Kazakidi, D. P. Tsakiris and J. A. Ekaterinaris (a) (b) Figure 2: (a) Undeformed and (b) final deformed squid arm. (a) (b) Figure 3: Comparison of experimental and numerical results for the squid arm extension: (a) tentacle length and (b) arm-tip velocity. with the experimental data (solid line) of the squid arm extension and the corresponding simulations (circles) obtained by Van Leeuwen and Kier [9]. It is observed that the overall agreement, both qualitative and quantitative, of the present FEM numerical results with the experimental measurements is very good. In Fig. 3(b), it can be noticed that a relatively lower arm-tip velocity is predicated through the proposed FE analysis. However, similar observations were made by other investigators [5, 8, 14], who used the same values of the material parameters. Next, a conical geometry resembling an octopus arm is considered. The arm is 10 cm long, extending along the zaxis, and has 1 cm root diameter, as seen in Fig. 4(b). The arrangement of muscles in the octopus muscular hydrostat is very different to that of the squid, and is depicted in Fig. 4(a). The musculature consists mainly of four groups of longitudinal muscles that extend along the arm and the transversal muscles. In addition, oblique muscles are present, which have helically aligned fibers around the arm, and the central axis of the arm is occupied by the axial nerve cord. Due to the lack of experimental data for the octopus muscular hydrostat, the same material parameters utilized for the squid arm simulation are taken for the octopus arm as well. The finite element mesh of the octopus arm consists of 420 four-node hexahedral and 8 155 V. Vavourakis, A. Kazakidi, D. P. Tsakiris and J. A. Ekaterinaris (a) (b) Figure 4: (a) Octopus muscular hydrostat structure: longitudinal muscles (L), transverse muscles (T), oblique muscles (O), axial nerve cord (N), and (b) finite element discretized model of the octopus arm 420 six-node triangular base prismatic elements (see Fig. 4(b)). The core of the truncated conical geometry contains the transverse muscles and the rest of the domain contains the longitudinal muscles. Oblique muscles are neglected in the present analysis because they have minor contribution to bending motion of octopus arms and they are hard to incorporate in robotic arm models. The root of the arm is allowed to move on the x-y plane and is fixed at the origin point (0,0,0), while the rest of the boundary is taken traction-free. The applied activation signal for the current simulation is a step function, as follows f(a)=         0,t⩽ti �1 2(1 + sin (πt /ta−π/2))�3.5,t⩽ti+ta 1,t⩽ti+td 0, t>t i+td ,(13) In Eq. 13 the activation time is set equal to ta=0.5 sec, while the total simulation duration is one second. Furthermore, tiand tdare the initialization and duration time of the activation function. As seen in the previous example, (Fig. 3(a)) the squid hydrostat can perform an extension maneuver if all transverse muscles are activated simultaneously. In order for the octopus arm to perform a bending or/and reaching move, primarily longitudinal muscles have to be activated. Initially, it is assumed that one group of longitudinal muscles is activated uniformly (ti= 0); then it is assumed that the same muscle is activated nonuniformly (ti=¯z/0.2, ti=¯z/0.6 and ti=¯z), given the normalized axial position ¯zof a material point within the muscle. The time duration of the activation level is equal to td= 1 sec. In Figs. 5(a) it is shown how the octopus arm deforms when one longitudinal mus- 9