Evaluation of contact definitions in a Finite Element model of the human cervical musculature

(1)

Evaluation of contact

definitions in a Finite

Element model of the

human cervical

musculature

Victor S. Alvarez

Degree project in

Solid Mechanics

Second level, 30.0 HEC

Supervisors:

Peter Halldin

Svein Kleiven

(2)

Evaluation of contact definitions in a Finite

Element model of the human cervical

musculature

Victor S. Alvarez

Degree project in Solid Mechanics Second level, 30.0 HEC Stockholm, Sweden 2012

Abstract

The human neck is especially vulnerable to severe injuries why it is of interest to gain a better understanding of the injury mechanisms involved. A 3D Finite Element (FE) model including the geometry of the individual neck muscles has been developed at Kungliga Tekniska Högskolan (KTH) and gives the possibility to study the strains inside the muscles. However, in this model the muscles are only hindered to interpenetrate and can glide relative each other with a roughly estimated coefficient of friction. The focus of this thesis has been to modify the FE‐model in order to model the connectivity on a physiological basis using methods available in the used FE‐software. The main reason for this was to get a more realistic muscle displacement without separations that are present in the current model.

The connective tissue surrounding and connecting the muscles was modeled with two principal approaches. The first was based on a combination of FE‐contacts where the material properties were coupled to the stiffness of the FE‐contacts. The other approach was based on a new element type called Smoothed Particle Hydrodynamics (SPH) elements in combination with FE‐contacts. Both approaches showed that the structural stiffness did not increase significantly, but the strain levels where somewhat elevated. Stability issues arouse with deleted elements and negative element volumes causing the FE‐contact based approach to fail prematurely. The SPH‐based approach had fewest deleted elements and completed the simulation but increased the calculation time with approximately 50 %.

It was concluded that the implementation of a connection between the muscles had a relatively low influence on the strains and kinematics of the neck and could be used to avoid muscle separations.

The influence on stability of the model was however more evident and the most stable result increased the calculation time significantly.

(3)

2

Utvärdering av kontaktdefinitioner i en Finita

Element model av den mänskliga

nackmuskulaturen

Victor S. Alvarez

Examensarbete i Hållfasthetslära Avancerad nivå, 30 hp Stockholm, Sverige 2012

Sammanfattning

Den mänskliga nacken är en väldigt sårbar del av kroppen och det är därför av stort intresse att lära sig mer om de olika skademekanismerna inblandade. På Kungliga Tekniska Högskolan (KTH) har en 3D Finita Element (FE) modell utvecklats som inkluderar musklernas geometri och ger möjlighet att studera töjningarna inuti musklerna. Men i denna modell hindras musklerna endast från penetration och kan glida relativt varandra med en grovt uppskattad friktionskoefficient. I detta examensarbete har fokus legat på att modifiera FE‐modellen för att modellera sammanbindningen mellan musklerna på en fysiologisk basis med hjälp av tillgängliga metoder i den använda FE‐mjukvaran. Huvudmålet med detta var att uppnå mer realistiska muskelrörelser utan separationer som uppstår i den befintliga modellen.

Den bindvävnad som omger och binder samman musklerna modellerades med två grundstrategier.

Den första var baserad på en kombination av FE‐kontakter där materialegenskaperna kopplades till styvheten i FE‐kontakterna. Den andra strategin baserades på en ny elementtyp kallad Smoothed Particle Hydrodynamics (SPH) i kombination med FE‐kontakter. Båda strategierna visade att den strukturella styvheten inte ökade med några signifikanta nivåer, men töjningsnivåerna ökade något mer. En del stabilitets problem uppstod som gav raderade element och negativa elementvolymer vilket ledde till att den FE‐kontakt baserade strategin kraschade innan simuleringen slutförts. Den SPH‐baserade strategin resulterade i minst antal raderade element och slutförde simuleringen men ökade beräkningstiden med c.a. 50 %.

Det fastställdes att implementeringen av de nya metoderna hade en relativt liten påverkan på töjningar och kinematiken hos nacken och kan användas för att undvika muskelseparationer.

Inverkan på stabiliteten hos modellen var dock mer märkbart och det stabilaste resultatet ökade beräkningstiden betydligt.

(4)

Introduction

The human neck is responsible for the mobility and flexibility of the head and the general kinematics of the head and neck is a product of all components that comprises it, such as muscles, vertebrae and vertebral discs. But this flexible construction also makes it vulnerable to trauma and since the neck connects the brain with the rest of the body this can lead to serious injuries. If it affects the central nervous system this can lead to partial or total paralysis and even death, as well as more diffuse damages to the spine or soft tissues that still can have grave impact on a person’s life. It is therefore important to study and learn more about the injury mechanisms of the head and neck which amongst other thing could lead to better protection systems. Traditionally the main sources of information aimed to understand injury mechanisms and develop protection systems have been experiments with volunteers, cadavers or detailed dummies. The problem with these methods is that they have many deficiencies and are often time consuming and expensive (Hedenstierna, 2008). A powerful supplement to the conventional methods is numerical computer models where loading conditions and system variables are easily controlled. They can also be used in cases where it is impossible to use voluntary experiments or cadaver tests.

The main problem with numerical modeling is the complex nature of human tissues with a non‐

homogenous material composition and very non‐linear material behavior. And in the muscles there is not only the passive elasticity to take into account as in conventional engineering materials, but the active part that contracts the muscle depending on signals being sent throughout the body. Several mathematical models of the human body have been constructed with different approaches and levels of detail. Amongst the ones that include the muscles of the neck, the predominant approach has been to use discrete spring elements with non‐linear constitutive relationships, where later, the active and viscoelastic properties have been integrated (de Jager, 2000; Brolin, 2002). With increasing computational power it has become possible to increase the complexity and detail of the models.

One of the latest models of the head and neck, is a 3D Finite Element (FE) model developed at Kungliga Tekniska Högskolan (KTH) with continuum material properties that has been first to include the actual geometry of the cervical musculature with the use of MR images (Hedenstierna, 2008).

This approach has many advantages such as improved kinematic compliance and possibility to study tissue injuries. It can also be used in educational purposes where different injury scenarios are simulated so that for example medical students can get a better understanding of the injury mechanisms. There is currently a project being developed at KTH partly for this purpose (STH, department of Neuronics, 2012). The model has had a good level of success in validation with volunteer testing (Hedenstierna, 2008) but some aspects can still be improved. For example, the activation pattern of the muscles still remains to be modeled on a more physiological basis. Another aspect is that the material model currently used is a combination of two existing program defined models available in the used software. A more complex material model might increase the certainty in the material behavior and thus the certainty in the response on a global and local scale.

Another aspect that has not yet been modeled in detail is how the muscle groups are connected to each other and other surrounding tissue. In the current KTH model there is only a simple interaction in the form of FE‐contacts that avoids interpenetrations and implements friction but with a roughly estimated coefficient of friction. This contact leads to the possibility of separations between the

(7)

6

muscle parts and is observed in collision simulations where the muscles buckle and separate on large scales. This is not considered physiological and makes it hard to use for educational purposes. An improvement on this aspect has the possibility of producing a more realistic muscle deformation and making the model more usable as a visualization tool. It is also possible that it changes overall passive stiffness of the neck and the strain distribution within the muscles.

Objective

The main focus of this thesis has been to find methods that can model the connectivity of the muscles on a more physiological basis in the 3D FE‐model created by (Hedenstierna, 2008) using the FE‐software LS‐DYNA. The main goal of this is to improve the displacement pattern of the muscles and hence reduce the muscle separation. A sub‐objective has also been to investigate what effects these changes can have on general kinematics and muscle strains. To confine the problem the restrictions have been to use contact definitions, elements and materials available in LS‐DYNA.

Dynamic Finite Element Method

The Dynamic FE‐method, similar to the conventional FE‐method, is a numerical approximate method for solving Partial Differential Equations (PDE). In solid mechanics this is used to solve the stresses, stains and deformation of a body. In the following section this is briefly explained and for a more detailed description the reader is referred to (Kleiven et al., 2001) or other literature on the subject.

The concept of stresses and strains can most easily be explained in one‐dimension, but can be generalized to three‐dimension with little effort.

Figure 1. A one‐dimensional body exerted to the body force q(x) and external forces at boundaries.

The stress σ(x) in a point in a body in figure 1 is defined as the force F(x) in a cross section divided by its area A(X) as

(1)

and the strain ε(x) is defined as

. (2)

where u(x) is the displacement. The stress and strain is then coupled with a constitutive relation as

L

q(x) x, u(x)

f f

dx

F(x) F(x+dx)

dx A(x)

dx+du(x) After deformation

q(x)

(8)

(3)

where E(x) is the material description. If the relationship between σ and ε is linear then E is called the Young’s modulus.

Combining equations (1)‐(3) with the equilibrium equation leads to a second order differential equation in the displacement u called the strong formulation of the governing equation (Kleiven et al., 2001). In the principal of virtual work this equation is multiplied with an arbitrary scalar function called a weight function v(x) and integrated over the entire domain resulting in the weak formulation and is for the one‐dimension case

0 0 . (4)

This equation could be solved with a displacement assumption, but for most boundary value problems it is not possible to find a displacement assumption that satisfies the governing differential equation and all the boundary conditions. Therefore the domain can instead be divided into subdomains, or finite elements, for which simpler assumptions can be made that also fulfills continuity over subdomain boundaries. For the one‐dimensional example the approximation for the entire domain can be written as

(5)

where u_i is a discrete value at point i, or nodal value, of u and N_i is a so called element shape function.

Applying this to equation (4) yields

(6)

where

(7)

and

0 0 (8)

Equation (6) can also be written on matrix form as

(9)

and if the shape function is chosen with certain properties the integration can be made over each element instead of the entire structure and the components are summed up in so called stiffness matrix assembly (Kleiven et al., 2001). The resulting equation however is still as equation (9) but with larger matrices. Similar derivation can be made for two and three dimensions.

(9)

8

A general three‐dimensional body as in Figure 2 can be approximated as an assemblage of discrete finite elements interconnected with nodal points.

Figure 2. General three dimensional body with an 8‐node three‐dimensional element (Kleiven et al., 2001).

The displacement in the local coordinate system x, y, z can then be written as

, , , , (10)

where N^(m) is the displacement interpolation matrix for element m and Û is the vector of the three global displacement components at all nodal points. Then the element strains become

x, y, z ∂ , , ∂ , , , , (11)

and the stresses

(12)

where C^(m) is the elasticity matrix of element m and σ^I(m) are the given element initial stresses. Now the principal of virtual of virtual displacement can be applied which states that for a body to be in equilibrium it is required that for any compatible virtual displacement, the total internal virtual work be equal to the total external virtual work (Kleiven et al., 2001). By including the element inertia and by approximating the element velocities and accelerations in the same way as the element displacement in equation (10), the principal of virtual displacement leads to the equilibrium equations as

(13)

where M is the mass matrix given by

(14)

where V^(m) is the volume of element m and K is the stiffness matrix given by

(15)

(10)

and F is the load vector as

,…

(16)

where f ^B(m) is the body forces for element m and f ^S(m) is the surface forces for element m. The equations in (13) can then be integrated in a step‐by‐step procedure where the time derivatives are approximated by taking the differences of displacements at various instants of time with a method of choice.

For large deformations, finite strains F_ij in multiple dimensions can be described by the deformation gradient as

(17)

where x_i is the deformed coordinate and X_j is the reference coordinate. A strain tensor commonly used for solid elements, where large deformations occurs (Hallquist, 2006), is the Green‐St. Venant strain tensor described by

1

2 (18)

where I_ij is the identity matrix.

FEContact

When describing the contacts in FE‐Analysis terms frequently used are slave and master. This is a way of distinguishing between two entities in contact and can be in terms of nodes, surfaces or segment (segments can be a general area comprised of a number of nodes). The reason to call it slave and master is because the nodes on the slave side are the ones that are in some way forced to follow the master side. This is however dependent on the formulation used and not always a necessity. The constraints used to force the bodies to interact are applied on the nodes, as all other external forces.

In LS‐DYNA there are three different methods for handling contacts between bodies. These are, the kinematic constraint method, the distributed parameter method and the penalty method (Hallquist, 2006).

The kinematic constraint method imposes constraints on the global equations by a transformation of the slave nodes displacements components along the contact interface. It is a very stiff formulation because the nodes are constrained to move with the surface. No way has been found to control this stiffness why this method is not considered to be of use in this thesis.

(11)

10

The distributed parameter mass method is also based on constraints, but uses a method of imposing constraints on the accelerations and velocities of the slave nodes to insure their movement along the master surface. This method is disregarded for the same reasons as the kinematic constraint method.

The last method will be the one of interest for the applications in this thesis since it allows for the stiffness to be controlled and is briefly described below.

Penalty Method

The penalty method checks each slave node for penetration through the master surface, see Figure 3.

Figure 3. Schematic figure of a contact surface (Kleiven et al., 2001).

If a slave node is found to penetrate, an interface force is applied between the contact point and the penetrating node. The magnitude of the force fs on the slave node is proportional to the amount of penetration l as

(19)

where n_i is the normal vector of the master segment i and k_i is the stiffness factor.

Hence this can be considered as placing interface springs between the penetrating nodes and the contact surface. The stiffness factor k_i of these springs can in LS‐DYNA be calculated in three different ways. The default setting is to express it in terms of the bulk modulus Ki, the volume Vi and the face area Ai of the element that contains the master segment si that is being penetrated as

(20)

which is for brick elements and fsi is the scale factor for the interface stiffness that can be changed but is normally defaulted to 0.1. For shell elements the stiffness factor is given by

max . (21)

There are also some options for setting the penalty stiffness value and they are:

‐ Using the minimum of the master segment and slave node stiffness. (default)

‐ Using master segment stiffness.

‐ Using slave node value.

‐ Using slave node value, area or mass weighted.

Master surface Material A

Slave node Material B

n_i

(12)

‐ Using slave node value, inversely proportional to the shell thickness.

The second approach is called soft constraint penalty formulation and is based on not only calculating slave and master stiffness according to equation (20) or (21), but also an additional stiffness called the stability contact stiffness k_cs given by

0.5 · · 1

Δ (22)

where SOFSCL is the scale factor for the soft constraint penalty formulation, m^* is a function of the mass of the slave node and the master nodes. ∆t_c is set to the initial solution time step but if the time step grows ∆t_c is reset to the current time step in order to prevent unstable behavior. This stiffness is then compared with the stiffness calculated with the first method and in general the maximum of the two is chosen as the used stiffness.

The third approach is called segment‐based penalty formulation and uses a stiffness as

0.5 · · 1

Δ (23)

where SLSFAC is the same scale factor as fsi and SFS and SFM are a scale factors for slave and master segment respectively defaulted to 1.0. The masses m here are segment masses rather than nodal masses and is for solid element segments equal to half the element mass. Here ∆t_cis also set to the initial solution time step, but it is only updated if the solution time step grows more than 5%.

LSDYNA Contact Definitions

There are many contact definitions in LS‐DYNA that uses different methods and approaches. The ones of interest in this thesis are described in this section. All contacts described use the penalty method.

The first contact definition used is denoted “CONTACT_AUTOMATIC_SINGLE_SURFACE” and will be referred to as singe surface for the remainder of the thesis. This contact only uses a slave side which means that it searches for all slave nodes that penetrate through any surface on the same slave side and applies an interface force on it but does not generate any force at separation. It is possible to apply a coefficient of friction and to use any of the penalty stiffness formulations above. The theory for the friction is not mentioned in detail since it is not a central part of this thesis but it is based on a Coulomb formulation and adds an interaction force to the nodes in contact. For more details the reader is referred to (Hallquist, 2006).

The second contact is denoted “CONTACT_TIEBREAK_SURFACE_TO_SURFACE” and will be referred to as tiebreak for the remainder of this thesis. This contact is originally developed to simulate failure which is not necessary for this application but can also be used to tie segments together using the penalty method unlike other tied contacts. The reaction force for this contact is proportional to the relative displacement between master segments and slave nodes and therefore it applies a penalty force both in separation, penetration and relative sliding. The distance between the nodes in contact can also be large and arbitrary. In this contact master and slave sides have to be defined and it is only possible to use the standard penalty stiffness formulation given by equation (20).

(13)

12

The last contact is denoted “CONTACT _AUTOMATIC_SURFACE_TO_SURFACE” and will be referred to as surface to surface for the remainder of this thesis. It works in a similar way as the single surface contact with the exception that a slave and master segment have to be defined.

Smoothed Particle Hydrodynamics

Smoothed Particle Hydrodynamics (SPH) is a so called mesh‐free computational method and has the advantages over conventional mesh‐based methods in that it provides a simple and accurate solution at large deformations and an easier discretization of complex geometries, amongst others. The idea is that the continuum is described by individual particles that move in space and carry all computed information, and thus the computational frame for solving the partial differential equations describing the continuum (Vesenjak and Ren, 2007). This is achieved by the fact that the exact value of any function f(x) of a three‐dimensional position vector x can theoretically be determined by taking the convolution of the function with the Dirac delta function δ(x‐xi) as

| Ω (24)

where

∞,

0, (25)

But since the Dirac function is theoretical and hard to use in practice a so called smoothing function, or Kernel function, W(x,h) is introduced in place of the Dirac function as a first approximation and is non‐zero in a small domain and zero elsewhere and can for example be defined as

, 1

(26)

where d is the number of space dimension and h is the smoothing length determining the influence domain of the smoothing function. When h goes to zero the kernel function can be approximated by the Dirac function. A kernel function commonly used is the cubic B‐spline (Hallquist, 2006) which is gained by choosing Θ as

1 3

2 | | 1

1

4 2 1 | | 2

0 2 | |

(27)

where C is a constant given by normalization and depends on the number of space dimensions.

Applying the Kernel function and further approximating the function f(x) to be solved numerically by discretization of the continuous domain Ω into small patches as

∆Ω (28)

where m_j is the mass and ρ_j is the density of particle j. The particle i has the position x_i(t) moving in the vector field v. This yields the so called particle approximation of the i‐th particle as

(14)

, (29)

where N is the number of neighboring particles and the approximation of gradients of a function is

(30)

Equations (29) and (30) are applied to the equations of conservation and give the local continuum properties as the sums of particles contributions. For example the momentum conservation equation is given by

1

(31)

where v is the velocity vector and α, β are the space indices, and the particle approximation becomes

, ,

(32)

where

1 || ||

. (33)

The energy conservation equation is

(34)

where P is the pressure and the particle approximation becomes

(35)

The SPH method has some numerical difficulties like particle inconsistency, inaccuracy at domain boundaries and instabilities at tensile stress state (Vesenjak and Ren, 2007). The stability and accuracy also depend on the particle density within the influence domain and the time step of time integration scheme. Another problem at large deformation is the change of number of particles in the influence domain. If they are compressed the number of particles increase, leading to increased calculation times. And at tension they decrease, leading to lowered accuracy and stability problems.

Therefore h varies in time and space in order to maintain approximately the same number of particles in the influence domain (Vesenjak and Ren, 2007). The influence domain is defined as 2h (Lacome, 2001) and the equation that governs the smoothing length is based on keeping the total mass M of n particles within a sphere with volume V constant. The mass in the sphere is given by

32

3 (36)

(15)

14

and too keep this constant over time equation (36) is differentiated in time as

32

3

32

3 (37)

where the left hand side is zero because of mass conservation and with some simplifications this leads to

1

3 (38)

where div(v) is the divergence of the flow. The initial smoothing length h₀is calculated by LS‐DYNA in the initialization phase and is computed by, for each particle, searching the shortest distance to another particle. The largest values of these distances is then denoted L and is then scaled with a factor defaulted to 1.2, giving the initial smoothing length. The initial smoothing length can also be set to a user defined value. The mesh of SPH‐element should be as regular as possible and not contain large discrepancies of the mass between the particles (Hallquist, 2006).

Cervical Region

The human backbone, or vertebral column, is divided into five regions (Seeley et al., 2005) and consists of bones called cervical vertebrae separated by intervertebral disks. The region closest to the head is called the cervical region and the region below is called the thoracic region. The main functions of the vertebral column are to support the head and body, protect the spinal cord, provide attachment sites for the muscles and permit movement of the head and body.

Figure 4. The upper part of the vertebral column from C1 to T2 viewed from the left side.

The vertebrae in the cervical region have smaller bodies than the rest of the column and the cervical region is also more flexible then the other regions making it more vulnerable to dislocations and fractures. The vertebrae in the cervical region are named C1 to C7 with C1 as the first cervical vertebra closest to the head and in the thoracic region they are named from T1 to T12 as seen in Figure 4.

(16)

Cervic

The cerv basic ph of individ The mus layer of muscle f contract tendons bones. S

Figure 5. S

Conne

Surround consists al., 2005 non‐fibr connecti short co return to The con constitut structura enclosin

cal Muscu

vical muscula ysiology as t dual muscle scle fibers ar connective fibers and a tile part of t that consist Surrounding t

Structure of the

ective Tis

ding all the mostly of ex 5). This mate

ous proteins ive tissue, co ollagen fibers

o their origin nective tissu ting compon al support in g and separ

ulature

ature consis the rest of th

fibers and a e long, mult tissue called are comprise the muscle t largely of c the entire m

e skeletal musc

ssue

specialized xtracellular m erial has thre

s, other mole ollagen fiber s that branc nal shape aft ue can be fou

nents depen n the form

ating. As m

ts of the mu he skeletal m re surrounde inuclear cell d endomysiu ed of two lo called the s collagen fibe muscle is a co

cle. Figure mod

organs in th material tha ee mayor com

ecules and fl rs which are ch and form er being stre und in many nding on loca

of bone, co entioned ab

uscles conne muscles. Each ed by a thin s consisting um, see Figu ong and ove sarcomeres ers and trans

nnective tiss

dified from (wik

he body is t t separates mponents: p

uid. Three d e somewhat a supportin etched.

y different pa ation and fu

nnecting tiss bove, it surro

ected to the h muscle con layer of con largely of my ure 5. Myofib rlapping pro (Seeley et a sfer the forc sue called th

kibooks, 2012).

he connecti the cells wit protein fibers ifferent type flexible, ret ng network a

arts of the b unction. Thes

sues to one ounds the sk

Muscl

cervical spin nsists of fasci

nective tissu yofibrils and brils run alo otein chains al., 2005). T ces from the

e epimysium

.

ve tissue. Th thin from on s, ground su es of fibers a icular fibers and finally e

body with va se functions another, cu keletal musc

e fiber

ne and has t icles that are ue called per d are surroun ong the lengt . These buil he fibers en

muscle fibe m, see Figure

he Connecti ne another (S

bstance con are generally which are v elastic fibers

rying amoun s can for exa ushioning, in les in three

the same e bundles imysium.

nded by a th of the d up the nd in the ers to the

5.

ve tissue Seeley et sisting of y found in very fine, that can

nts of the ample be nsulating, different

(17)

16

levels and in some sense creates compartments that separate these muscle components. This can be seen in Figure 6 where the cross‐section of a muscle without its muscle proteins is shown.

Figure 6. Image from scanning electron micrographs after removal of skeletal muscle protein. Top left shows the epimysium (EP) and in the bottom the perimysium (P) and the endomysium (E). The top right shows an enlargement of the endomysium surrounding an individual muscle fiber (Kjaer, 2004).

According to a recent study (Guimberteau et al., 2010) the connective tissue is present as a histological continuum throughout the body with no clear separation between the skin, the hypodermis, the vessels and the muscles. Although structures that allow sliding with very little influence on surrounding tissue are present everywhere, as for example seen on the skin when contracting a muscle. They propose that this continuum consists of a network of microvolumes and microfibrils in a chaotic pattern that adapts when exposed to forces to preserve energy and fill space in an optimal manner, see Figure 7. But since the function and composition of the connective tissue can vary extensively throughout the body, so can the allowed sliding. It is widely accepted that the main contributor for the transmission of the force generated in the muscle to the bones at its binding sites is the myotendinous junction. These are sites that connect the muscle with the extracellular connective tissue (Järvinen et al., 1991). The muscle fibers have also been shown to be connected to the endomysium and the perimysium creating pathways that can transmit forces within the muscles (Purslow, 2002) and therefore the muscle fibers cannot be seen as independently working units.

Figure 7. The multimicrovacuolar system of the connective tissue (Guimberteau et al., 2010).

(18)

It has also been shown that forces can be transmitted to adjacent muscles via the epimysium and the extracellular matrix of a muscle to surrounding non‐muscular elements of so called compartments containing a group of muscles (Maas et al., 2001) and even antagonistic muscles (Huijing, 2009;

Yucesoy et al., 2010). A consequence is that length–force characteristics may differ for an isolated muscle in situ and the same muscles in vivo. That is to say it is dependent on not only its own material properties but also on the surrounding muscles. It is also dependent on the properties, configurations and connective tissue structure in its direct environment. The resulting mechanical behavior of the muscles is therefore determined by the properties of the fibers, the extracellular matrix and the interactions between these (Gao et al., 2008b). In order to study the extramuscular force transmissions a FE‐model was created by (Yucesoy et al., 2003). The muscle extensor digitorum longus from the left anterior crural compartment of a rat was modeled with three times six, solid eight node elements. These were connected with 2‐node spring elements with uniaxial and linear stiffness characteristics to a set of fixed points representing the adjacent muscle. The simulations were performed in the same manner as the experiments and the value of the spring stiffness k was determined so that the simulations provided a good agreement with the experimental results, giving a value of k=0.067 unit force/mm.

In another recent study (Gao et al., 2008a) the epimysium is modeled as a material consisting of ground substance containing pairs of wavy collagen fibers. These are distributed in different angles to the muscle fiber direction and allow straightening and reorientation of the collagen fibers. The model predicts experimental tests fairly well as seen in Figure 8.

Figure 8. Comparison of the stress‐strain relationship of the epimysium between the experimental data (Gao et al., 2008a) and the prediction.

The experimental studies were performed on young and old rats (Gao et al., 2008a). Small sections of the epimysium were dissected from tibialis anterior muscles of rats and each end of the samples were attached to plastic holders with glue and stretched at a velocity of 10 mm/s. It was shown that for low strains the main contributors to the tensile strength are the elastin and the ground substance. At high strains the collagen fibers are straightened and become the main contributor to the tensile strength. These experiments show the tensile strength in plane of the connective tissue, see Figure 9, and though it is not specified in what direction it is assume that it can be approximated as transversely isotropic. However the tissues were also analyzed with electron microscope and thee different layers were identified where the outermost (facing the epimysium of the adjacent muscle) consists of bundles of collagen fibers aligned in a dominant direction, see Figure 10.

(19)

18

Figure 9. The stress‐strain relationship of the epimysium (Gao et al., 2008a).

Figure 10. Schematic representation of the ultrastructure of the epimysium in rat TBA muscle. Inner is facing the muscle, outer is facing epimysium of adjacent muscle (Gao et al., 2008a).

The thickness of the epimysium was also measured and was almost the same for old and young rats, 29.7 ± 8.6 and 34.6 ± 14.6 µm respectively. From this study it is possible to extract an approximate Young’s modulus of 0.7 Mpa, see Figure 11.

(20)

Figure 11.

2008a).

In a stud (see Figu cadavers different the heal with a si displace

Figure 12.

These te standard

. Young’s modu

dy performe ure 12) was s of humans t points and thy specime ize of 5 by 10

ments at a ra

Figure of the h

ests showed d deviation o

ulus E extracte

d by (Osamu extracted fro s without th

the average ens was 0.48 0 mm was m ate of 1 mm/

human carpal t

that in heal of 1.8 kPa.

ed from the slo

ura et al., 20 om human p his condition e was consid mm with a mounted on

/s

tunnel. Figure m

lthy humans

ope of the curv

007) a piece patients with n. The thickn dered as the standard de plastic plate

modified from

s the shear m

ve from the str

of the SubS h so called ca ness of the s thickness of eviation of 0.

s and stretch

(thefreediction

modulus had Car

ress strain curv

Synovial Con arpal tunnel

samples wa f the sample 12 mm. The hed until fail

nary, 2012).

d a mean va pal tunnel

ve plotted by (

nective Tissu syndrome, a s measured e. The mean e faces of SSC

lure under c

lue of 2.7 kP

(Gao et al.,

ue (SSCT) and from at three value for CT pieces ontrolled

Pa with a

(21)

Figure 13.

healthy ca

As ment with ver the finge properti that it is moveme indicativ other via

Dissec

To get a nature o were de The con 2010), th

Figure 14.

. Representativ adaver specime

tioned befor ry little effec ers should be

es of the co s unclear if ent in which ve. But if all,

a the epimys

ction

a better und of the connec

livered with nective tissu hat covers al

Picture showin

ve shear stress en (Osamura et

re it is proba t on surroun e able to mo onnective tis the connect h case the s or some of sium, this val

erstanding a ctive tissue,

head and th ue truly seem

l muscles, te

ng the connect

against shear t al., 2007).

able that for nding tissue, ove independ sue in the n tive tissues shear modu f the neck m

lue should p

and a hands the dissectio he outer skin ms to be a co

endons etc. a

tive tissue on th

20

strain curves.

r tendons th which beco dently. Altho neck, they co in the neck lus taken fr muscles trans robably be h

s‐on view of on of roe dea

removed.

ontinuous ne as a very thin

he neck of a ro

a) Patient spec

he connective mes especia ough none of ould be used have the fu rom (Osamu smit a large higher.

how the m ar necks was

etwork, as su n layer (see F

e deer.

cimens with ca

e tissue sho ally evident in

f these studi d as indicatio

unction of a ura et al., 2 amount of f

uscles are ti s performed

uggested by Figure 14).

arpal tunnel syn

uld allow m n the last stu ies show the on. The main allowing inde 007) might forces betwe

ied together . The dissect

(Guimberte

ndrome, b)

ovement udy since e material n issue is ependent be fairly een each

r and the ted necks

eau et al.,

(22)

It covers each individual muscle but also groups of muscle and finally the whole neck. Although to say that it covers may be misleading since it is not individual capsules, but rather as if the space between the different tissues is filled with this material. The connective tissue can also be said to be a continuation of a certain tissue that gradually changes composition as suggested by (Guimberteau et al., 2010). For example, in the center of a muscle the material composition is as described for muscles but towards the outer part of it the composition changes, the active fiber components disappears and the relative amount of collagen fibers, elastin fibers and background material changes and the material becomes connective tissue, but no clear line can be found. Continuing towards another type of tissue, say a tendon, the composition continues to change and transforms into the tendon tissue, which basically just has a different relative amount of the components that make up the connective tissue (Seeley et al., 2005). How fast this transition is probably differs between tissues, anatomical position and ultimately physiological function but is quite hard to determine. This appears to be the same for all other tissues as blood veins and fat which fills up some of the space in‐between the muscles making all off the tissues one continuous entity.

The resulting mechanical behavior of the interactions between the muscles is, because of the nature of the connective tissue, quite complicated. But during examination of the roe deer neck one observation made was that the muscles appear to be able to slide in relation to one another with little or no impact on the neighboring muscle when the relative displacement is not to large. This can be seen in Figure 15 where a muscle is in rest in the first picture and deformed in the next.

Figure 15. Figure showing the effect on surrounding tissue when deforming a muscle piece. No interference in A and a muscle compressed in B.

At the same time the muscles are packed relatively tightly and even if it is possible to drag them apart it is hard to determine how strongly they are tied. It is also hard to determine how much of the strength lie in the thin layer of connective tissue separating the pieces and how much is a result of all the tissues and complex intertwining of them. This makes it very hard to determine the stiffness normal to the muscle tissue. As an analogy of the response to interaction, one can imagine a thin plastic bag filled with muscles and some moist which is emptied of air. The pieces inside can be moved relative to each other but cannot be truly separated since there is no air in the system. To be able to see how the muscles are bonded, some of the fabric must be destroyed, as seen in Figure 16.

Compressed muscle tissue No deformation

A B

(23)

22

Figure 16. Picture showing the layer between to muscle pieces.

If a small piece of connective tissue is removed from the muscle it reveals that it in fact is very stiff when stretched and also brittle, but the brittleness is probably because it is so thin.

Combining the ideas of these different articles and dissection, it could be assumed that the tensile stiffness in the plane comes from the outermost parts of the connective tissue. But since the amount of fibers decreases towards the middle this is probably what is responsible for the weakness in shearing. The force transmission mentioned by (Maas et al., 2001) is probably varying in different parts of the body but much of it could be assumed to come from the shear stiffness. No information has been found about the stiffness of the connective tissue perpendicular to its surface, though it should probably be lower than the stiffness in the plane of the surface since the shear stiffness is much lower. Nevertheless the resistance against separation, i.e. the apparent stiffness between the muscles is suggestively much higher than the stiffness perpendicular to the plane. The reason is that the resistance would be a result of the entire connective tissue network and not only of a small sheet of tissue between the muscles, since a large portion of the resistance can come from tensile stretching in the plane of the epimysium.

Methods

As for many other soft biological tissue the material properties of the connective tissue is complex, non‐linear, non‐isotropic etc. To implement this material in the current neck model with conventional FE‐methods would be difficult since it requires a very thin sheet of elements in‐

between the muscles if not to remodel them completely. To avoid extremely small brick elements this would have to be done using shell elements but would then not correctly model the stiffness of the connective tissue normal to the muscle surfaces since they are two dimensional. But from the information above a few ideas arose concerning how to mimic the connectivity amongst the tissues in the neck with other approaches which are presented below.

Tiebreak Contact

One approach is to use the predefined contact definitions available in LS‐DYNA and manipulate different options in order to get appropriate properties. The stiffness can then be represented through the penalty stiffness and will be a linear stiffness which the connective tissue is evidently not.

But since the total stiffness of the neck (with vertebrae, muscles etc.) is several times larger this

(24)

would probably not be very noticeable. The stiffness of the connective tissue is also fairly linear in a range of strain that the tissue could be assumed be exposed to during these simulations.

The simplest way to implement this idea is to use a tiebreak contact between the surface nodes of the solid elements and manipulate the penalty stiffness. But since there is only one penalty stiffness factor this could lead to interpenetration amongst the muscles if the stiffness is set to low. From the literature studied it is not clear if the connective tissue has the same stiffness in normal and tangential direction, but from the dissection it is assumed that the apparent resulting stiffness is not equal in shearing and in separation.

The first issue to avoid a linear stiffness cannot be countered without making a new contact definition, and ports outside of this thesis. Interpenetrations on the other hand can be countered by adding an extra single surface or surface to surface contact that prohibits the penetration but does not affect the stiffness in separation and sliding. To easily implement the single surface contact the solid elements can be covered with shell elements of very low thickness and a so called null material that has very little influence, something that is already used in the current model to implement the single surface contact. Another way to implement these extra contacts is to create segments on the same nodes as these shell elements but removing the actual shell elements afterwards. With the surface to surface contact a master and slave surface have to be defined.

One practical difficulty in implementing the tiebreak contact is that every individual muscle part needs to be set as slave or master in contact with surrounding muscles. And since several muscles are in contact with multiple muscles this leads not only a tedious work but also to the fact that many muscles will have multiple contact definitions on single nodes. Moreover, since the tiebreak contact ties nodes with large separation, the nodes on the side of the muscle not facing the adjacent muscle can get tied to the contact surface. To counter these problems a small MATLAB program was created that finds the nodes within a certain distance and outputs the master segments for each part and the slave segments within the given distance, see Appendix A.

SPHElements

The other approach is to make use of the SPH‐elements and the fact that they describe a continuum material without actually being tied together. This makes it possible to model the very narrow space confined between the muscles as a 3D material. An easy way to achieve this is to define the nodes on the surface of the muscles as SPH‐elements and set the volume of influence so that the neighboring nodes on the same muscle and on the adjacent muscle fall inside this volume. The SPH‐elements can then be given a material property of choice to describe the connective tissue. However, not all materials were found to be compatible with the SPH‐elements but two materials that are suitable for the application and were found to work are an orthotropic linear‐elastic and a simple isotropic linear‐

elastic. With this modeling approach there are also some risks of interpenetration if using the isotropic material for low values on the Young’s modulus that can be countered in the same way as with the tiebreak‐based approaches above using additional contact definitions.

One issue with this approach is that it was found that, when using the SPH‐elements, the simulation time increased quite significantly. In an effort to counter this, another MATLAB program was created, see Appendix B, that finds only the nodes on the neighboring muscles that lie within a certain

(25)

24

distance, so to reduce to total amount of SPH‐elements used and thereby the simulation time.

Another difficulty with these elements is to implement the orthotropic material. This is because no good way was found to define the material axes since, in LS‐DYNA, this is based on the nodes in the corners of brick elements if they are to be local. The only solution to this was to use a global coordinate system roughly aligned with the desired material axes. The particle distribution in the initial SPH‐mesh is a problem hard to counter since the nodes on the surface of the muscles are not uniformly distributed. Finally, the last issue is that no data was found on the relative amount or density of connective tissue surrounding the muscles and therefore it was assumed to be the same as the muscles, and the mass is then set to the density times the average volume between the muscles.

Since the added mass will be fairly low and the strength of the connectivity is not affected by this, it is not considered a big issue.

Modeling Approaches

Based on these ideas a few approaches where created to model the connective tissue in order to compare their performance and find advantages and disadvantages.

Single surface

The first approach simply uses what is in the current neck model which is a single surface contact applied on shell elements comprised by the surface nodes of the solids muscle elements, as depicted in Figure 17. A coefficient of friction of 0.3 is used and the approach is denoted

“single_surface_original”. This approach is also used without friction and is then denoted

“single_surface_nofric”. The idea with these contacts is to have a reference that can show how the new approaches affect the behavior of the overall model. Both using friction and not using friction can show how much the friction in the current model is actually affecting the behavior.

Figur 1:

Figure 17. Schematic figure showing the nodes on the surface of solid elements constituting also shell elements.

Tiebreak

The second approach uses a tiebreak contact implemented on surface nodes of the solid element and is denoted “tiebreak_E_p”, where E_p is substituted with the equivalent value of the Young’s modulus produced by the penalty stiffness. As for all approaches using the tiebreak contact, the nodes used in this contact are found with the MATLAB program previously mentioned.

Solid element

Shell element Node i

(26)

Tiebreak and single surface

The third approach is to use the same tiebreak contact on the solid elements but also the single surface contact without friction on the shell elements from the original model. This approach is denoted “tiebreak_singlesurf_E_p” where the E_p is the same as above. The approach is also implemented by replacing the shell elements with segments of the same nodes, for which it is denoted “tiebreak_singlesurf_segments_E_p”. An extra approach was also created using the same shell elements only without the single surface contact applied on them and it is denoted

“tiebreak_and_shells_Ep”. This approach does not try to model the connective tissue in any new way, but is rather used in order to determine the effect of using multiple element and contact definitions on the same nodes. The extra contact used in these approaches could also be implemented using the surface to surface contact, but it is basically the same contact as the single surface. Although it could give a slightly different result, because of searching algorithms and other unknown aspects, it would not be possible to determine if it was a better way of modeling the connective tissue. It is also less convenient to implement because a master and slave need to be defined which would be tedious work in the full neck model and is for these reasons omitted.

SPH and single surface

The fifth approach uses the SPH‐elements with the isotropic linear‐elastic material and has the same single surface contact on shell elements as the first approach. It is denoted “SPH_isotropic_singleurf_

Ep” where the Ep is the Young’s modulus used in the material definition for the SPH‐elements. This approach is, as the third approach, also implemented using segments instead shell elements for the single surface contact and is then denoted “SPH_isotropic_singleurf_segments_E_p”. The reason to use an isotropic material was to see if there are any advantages with this approach in comparison to the contact based ones.

SPH

The sixth approach is also based on SPH‐elements but with the orthotropic material and is denoted

“SPH_orthotropic_Eab_Ec”, where Eab are the Young’s moduli for the two global material axes roughly in the plane of the muscle surface facing the adjacent muscle surface. Ec is the Young’s modulus in the direction normal to the same surface.

The two last approaches are divided into two groups, one with SPH‐elements on all the nodes covering the whole muscles surface indicated by a “full” at the end of the model approach description but before the Young’s modulus, as “SPH_***_full_ Ep” and one with only the nodes in‐

between indicated in the same manner but with “red”. This is mostly to verify a relationship between the calculation time and the number of SPH‐elements.

The reason that the orthotropic SPH approach does not have a version with contact is that the idea is to have the stiffness in normal direction high enough to prohibit interpenetration on its own. This approach is included mostly to show on possible use of the SPH‐elements as it is not likely to function very well being that the orthotropy will not follow the muscle movement.

(27)

26

Material Parameters

Basis for the choice of parameters was the shear modulus G minus the standard deviation resulting in G=0.9 kPa found by (Osamura, et al., 2007), the Young’s modulus of E=0.7 MPa extracted from (Gao et al., 2008a) and the spring stiffness of 67 N/m determined by (Yucesoy et al., 2003). To transform the shear modulus to a Young’s modulus, and vice versa, a homogenous and isotropic material is assumed, since the used methods are based on a linear elastic approach, and the relationship is given by

2 1 ν (39)

where ν is the Poisson’s ratio and is set to 0.49 that is the value used for the muscles. This is since no value was found for the connective tissue but is also probable since it is intimately connected to the muscles. From the shear modulus of 0.9 kPa the equivalent Young’s modulus is approximately 2.7 kPa and from the Young’s modulus of 0.7 Mpa the equivalent shear modulus is 0.2 MPa. The values of the Young’s modulus can then be coupled to the penalty stiffness in equation (20) with the following scheme. If the space between two 8‐node cubic solid elements is modeled with springs the total force f between the solids can be described as

4 Δ (40)

where ∆x is the displacement in the direction of the force, the factor 4 is because there are four nodes on the surface and k is the spring stiffness which is equivalent to the penalty stiffness k_i. Equation (40) is divided with the distance t between the solids and the area A of the solid element as

4 Δ

(41)

where

(42)

and

Δ

. (43)

The distance t should be the thickness of the samples used in the material tests performed by (Osamura et al., 2007) and (Gao et al., 2008a). But since both values will be used in a range it is chosen as an educated guess to 0.2 mm which lies closer to the study by (Gao et al., 2008a) since it is for the epimysium of a rat muscle and the other study is from the human hand but not the epimysium.

Combining equations (41), (42) and (43) with Hooke’s law given by

(44)

and some minor manipulation the relation between the penalty stiffness and the Young’s modulus Ep becomes

Evaluation of contact definitions in a Finite Element model of the human cervical musculature

Evaluation of contact

definitions in a Finite

Element model of the

human cervical

musculature

Victor S. Alvarez

Degree project in

Solid Mechanics

Second level, 30.0 HEC

Supervisors:

Peter Halldin

Svein Kleiven

Evaluation of contact definitions in a Finite

Element model of the human cervical

musculature

Victor S. Alvarez

Abstract

Utvärdering av kontaktdefinitioner i en Finita

Element model av den mänskliga

nackmuskulaturen

Victor S. Alvarez

Sammanfattning

Contents

Introduction

Objective

Dynamic Finite Element Method

FE­Contact

Penalty Method

LS­DYNA Contact Definitions

Smoothed Particle Hydrodynamics

Cervical Region

Cervic

Conne

cal Muscu

ective Tis

ulature

ssue

Dissec

ction

Methods

Tiebreak Contact

SPH­Elements

Modeling Approaches

Material Parameters

FEContact

LSDYNA Contact Definitions

SPHElements