Annual Technical Report No. VSR-04.15

Biomechanical analysis of skeletal muscles is an important task in the development of a digital human system. Standard finite element method (FEM) can be used for muscle modeling, however, in an interactive digital human environment an FEM muscle can be unnecessarily complicated and impractical in real-time interaction. In this project, we developed a new modeling method based on the combination of Non-Uniform Rational B-spline (NURBS) geometric representation and the Galerkin methods, which we refer as NURBS-FEM. NURBS primitives are the de facto industrial standard for geometric modeling, graphics and graphical information exchange. The basic idea of NURBS-FEM is to utilize the NURBS primitives that come with a typical geometric model of muscle to construct the kinetic equations of the discretized system. In other words, we derive finite element equations directly on the basis of existing NURBS geometric description, without additional meshing. Salient features of the method include also (1) it provides a better (smoother) geometric description, and (2) it involves less degree of freedoms. The method by itself is sufficiently general, applicable to general deformable bodies and has the potential of real-time simulation. In addition, the formulation retains most the desirable properties of the finite element method.

We have developed a fully nonlinear formulation of NURBS-FEM suitable for the analysis of large-strain motion of general deformable body. The formulation is validated quantitatively or qualitatively against motions of which either the analytical solution or some features of the motion are known.

As an indispensable step towards biomechanical analysis of muscles in an interactive digital human environment, we have developed isolated models for representative human muscle in the upper limb. The muscle surface is extracted from the Visible Human Data Set [16] and represented using NURBS primitives. Once the surface representation is obtained, the muscle body is parameterized as a NURBS solid. Muscle is active, anisotropic material. For now in this project, the muscle’s stress response is described by a hyperelastic, anisotropic constitutive equation that includes an active component. Using the isolated model, we have simulated the motion of the muscle under different load condition and activation level. Two types of characteristic contractions of skeletal muscle, namely isometric and isotonic contraction, are tested. The simulation results agree qualitative with the characteristic features of the contractive motion.

The most important future work is to integrate the muscle model into Santos and interacts with other components. Outstanding technical developments on the method itself include contact treatment, the development of a full scale model for the upper limb that include all major muscles in mutual contact, and better constitutive description that takes into account of strain rate effect and other biological features.

IV.1 Introduction

Muscle modeling and simulation is an important task in biomechanical analysis of digital humans. Standard finite element stress analysis can be used for muscle modeling, however, in an interactive digital human system an FEM muscle can be unnecessarily complicated and impractical in real-time interaction. In this project, we developed a new modeling method based on the combination of NURBS geometric representation and the Galerkin methods, which we refer as NURBS-FEM. The method can adequately predict muscle stress and deformation, and at the same time, keeping the model size and complexity at a tractable level. The method by itself is sufficiently general, applicable to general deformable bodies bounded by smooth surfaces. Being essentially a finite element method, the formulation retains most the desirable properties of the latter. In particular, different constitutive laws (e.g. materials) can be implemented easily, allowing us to experiment different biomechanics models for skeletal muscle.

Over the past decades, physically based deformable modeling has been an interest for the computer graphics community, in medical applications and other technical fields. Several major methods have been proposed [12]: mass-spring models, finite element methods, finite volume models and other low degree approximated continuum models. As an example of deformable objects, muscle is modeled with various approaches. Porcher-Nedel and Thalmann used an action line to represent the force produced by the muscle, and a surface mesh of muscle as a mass-spring network [9] which can be deformed under applied force. Ng-Thow-Hing used the trivariate B-spline solid to model individual muscle in animals and humans [17]. The muscular deformation is also achieved by the embedding a mass-spring-damper network defined in the B-spline solid. The mass-spring model is simple and fast, but less accurate. A recent trend is to use FEM to simulate muscle behavior based on the realistic muscle shape and pointwise stress-strain relations. Chen and Zeltzer [3] developed a finite element model of skeletal muscle to simulate muscle forces with Zajac’s dimensionless biomechanical model of muscle [4,5]. The low DOF prismatic bounding box of muscle shape is used as the mesh for finite element simulation. The resultant embedding muscle deformations were visualized by free-form deformation [13] defined by the mesh. J. Teran and S. Blemker etc. [8] used finite volume method to perform rigorous large deformation analysis. Tetrahedral meshes of the biceps and triceps are generated form Visible Human Data Set [16]. As alluded before, FEM method is suitable for stand-alone analysis, but can be prohibitively expensive in an interactive environment because of the resulting large number of DOFS. Some researchers have explored ways to do physically based models with low DOF by adding physical behavior to traditional parametric surface patches. One of the excellent works is Dynamic NURBS (D-NURBS) developed by Terzopoulos and Qin [2]. D-NURBS is a physically based dynamics model where the NURBS control points used as generalized coordinates, and the dynamical equations are derived from the Hamiltonian principle.

The basic idea of NURBS-FEM is to utilize the NURBS primitives that come with a typical geometric model of muscle to construct finite element discretization. In other words, we sought to derive finite element equations directly on the basis of existing NURBS geometric description, without additional meshing. NURBS is widely used to represent free-form curve and surfaces. It has become the de facto industry standard for the representation, design and data exchange of geometric information processed by computers.

			Control points

(a)

(b)

Figure. 1 NURBS curve and surface (a) NURBS curve (b) NURBS Surface

As illustrated in Figure 1, the NURBS geometry is defined by the positions of its control points and parametric mappings known as basis or shape functions. The NURBS shape can be changed by adjusting the positions of control points or the shape functions. For a geometrical standpoint, the deformation of a NURBS body can be described by displacing the control points. In NURBS-FEM, the displacement field constructed from NURBS interpolations is directly used to facilitate mechanical analysis. The governing equations for the motions of the control points are derived from mechanical balance laws using the Galerkin method. The formulation is essentially an isoparametric finite element with the NURBS function employed to define both the geometric mapping and deformation (namely the deformed geometry). The formulation bears some similarities with the meshfree method [14] in that the shape functions do not have the d-property and hence, the “nodal” parameters, namely the control positions, are not the nodal values of the dependent variable. On the other hand, unlike meshless method, the nodal points do not just span a meshless interpolation function, they still define elements. Indeed, the domain is still discretized into a patch of elements. In contrast to the standard FEM, the NURBS Control points, including boundary and inner control points, are generalized nodes not necessarily lying on the domain. However, the “nodes” still carry a clear geometric connotation since they control the shape of the body.

Muscle is the motivation for us to develop the NURBS-FEM, and on the other hand it’s also a good application for NURBS-FEM due to muscle’s geometric physics complexity. Muscles are anisotropic solids that would be thought as a nonlinear elastics material embedded with active contractive fibers. The surface shape of major muscles in human limb can be readily obtained from the Visible Human Data. In this project, we use trivariate NURBS solid extended from the NURBS surface of muscle boundary to represent the muscle geometry. An anisotropic, hyperelastic constitutive equation with an active component is used for muscle stress analysis. The directions of the contractile fibers can be represented as the tangent of appropriate NURBS curves in the muscle solid.

This report is organized as follows. Section 2 contains a brief review of trivariate NURBS solids and how could it be extended from NURBS surface. In section 3, the NURBS-FEM formulation is discussed in the context of linear and nonlinear elasticity. Since the procedure for 2D case has been discussed in the half year report [19], we will focus on the 3D case and delineate the details of both NURBS interpolation and the derivation of finite element equations. Section 4 contains numerical examples showing the application of NURBS-FEM in 3D, large deformation analysis. Section 5 gives a detailed account of the muscle model used in the project. In section 6, numerical simulations of muscle contraction are presented. Some concluding remarks and a brief plan for future works are given in Section 7.

IV.2 Trivariate NURBS Solid

In computer graphics, 3D solid objects are typically represented by their bounding surfaces. However, surface boundaries do not contain information about the interior of 3D objects. The problem of defining a solid from its surface has long been identified as an important research topic. One of the approaches is the trivariate parametric solid [6].

IV.2.1 Trivariate NURBS Solids

A NURBS solid representation is the generalization of NURBS representation of curves and surfaces. It defines not only the surface boundary of an object but also its interior. A position of a generic point in the solid is defined by

(2.1)

where are the control points, are the weights associated with the control points, and is the shape functions. The first index relates to the number of the control point, and the second index refers to the degree of the shape function (see later). Here are three parametric variables defined on three non-decreasing sets of knots

The valid parametric domain is and .

The shape functions for parameter are defined recursively as

where

The same applies to shape functions for parameterand .

Example

Some examples of B-spline and NURBS shape functions are shown in figure 2. Here, “rational” means weights are not all equal to 1, and “uniform” means all the interval of knot vector have equal length.

(a)

(b)

(c)

(d)

Figure. 2 NURBS Shape functions in one parameter (u/v/w)

(a) Degree 1 uniform B-spline shape functions (all weights equal to 1) with knot vector .

(b) Degree 3 uniform B-spline shape functions with knot vector ;

(d) Degree 3 non-uniform rational B-spline shape functions with knot vector and weights

It is clear from Equation (2.1) that the shape function in 3D interpolation are obtained by the trivariate tensor-product of B-spline shape functions defined by the knot sequences. Introducing the piecewise rational shape functions

the NURBS solid representation (2.1) can also be written as

(2.2)

Properties of NURBS Shape functions

Some important properties of the functions, which are the cornerstone for the NURBS-FEM, are summarized below.

(Partition of unity): The shape functions satisfy

for all and .

(Local support): if
In any given domain , there are at mostshape functions which are nonzero. In particular, the are nonzero for and .
Non-rational B-spline, Bezier solids are special cases of NURBS solids. If for all , then the trivariate NURBS solid becomes a trivariate B-spline solid.

IV.2.2 Extension from NURBS Surface

In most applications the solid shape is given in terms its bounding surface, given for example with free-form NURBS surface

Adding another parametric domain and its corresponding new control points, a NURBS solid represented by equation (2.2) can be obtained. From the perspective of computer aided geometric design, an extended NURBS representation means that all the exiting algorithms can be applied in the domain of . The methods that are commonly used in surface construction, such as ruling, skinning, sweeping and swing, can be applied to solid model construction [6]. In addition, a method called shrinking [1] has been used frequently when dealing with a close or periodic NURBS surface. A NURBS solid can be constructed by shrinking a surface to a point or a curve. For example, a spheroid is derived by shrinking a sphere to its center point, and a cylinder solid is obtained by shrinking a cylinder surface to its centerline.

IV.3 NURBS Based Finite Element Method

IV.3.1 Isoparametric Mapping

The basic idea of NURBS-FEM is to use the NURBS representation as an isoparametric mapping in finite element formulation. In equation (2.2), is the geometric shape function for a NURBS solid, which maps the parametric domain

{: |,}

to the real geometry. For any point, there is a spatial material point in the NURBS solid corresponding to it. In light of this, is called the parametric coordinates of the spatial material point. Under certain restrictions on the position of control point, the mapping is one-to-one, meaning that for a given parametric point, there is one and only one corresponding material point.

In Figure 3(a), the cube is the entire parametric domain, denoted by. The whole cube is partitioned into a patch of small cubics, which are called parametric elements or mesh. The partition of the whole parametric domain into elements is constructed by the knot spans of the knot vectors in and. In other words, any small domain in the form of the trivariate tensor product , denoted by , is an element in the parametric domain. In order to obtain non-degenerated element, some restrictions have to be placed on the knots values. For example, and . The details will not be discussed here.

Remark 3.1

In fact, each parametric element is mapped into a real geometric element, such that these geometric elements are the real mesh of the solid, as shown in figure 3(b). It’s because that every parametric coordinate is to be mapped into one and only one real point in the solid, which ensures that no intersection or overlapping between the real elements. In mesh generation community, this mapping has been used to generate finite element mesh. However, In NURBS-FEM formulations, the real spatial mesh will not be used; instead, the governing equations will be developed directly for the control variables. The numerical computation will be performed in the parametric domain.

From the property 3 (in Section 2.1), it is easy to show that in any parametric element , there are at most control points (,) having influence on it, namely the corresponding shape functions of are nonzero at this element. Similar to standard finite element method, here and associated with the parametric element are treated as ‘nodes’ and ‘shape functions’ of this element. Also from the local support property of shape functions, any control point has limited influence domain of , which may at most include parametric elements. In other words, the elements in a sub-domain will share the control point as their common node.

The key idea of NURBS-FEM is to use the interpolation (2.2) as an isoparametric mapping. We represent the coordinates of a generic material point in the reference (i.e., undeformed) configuration as

where is the coordinates of control point . AT the same time, the coordinate of the same material point in the current (i.e., deformed) configuration is constructed by

where is the coordinates of after deformation. Thus, the displacement field, defined as the difference of the current to the reference position, is given by

whereis the displacement of the control point , and is the displacement of a spatial point whose parametric coordinates is .

For the convenience of notation and the convention of FEM, in the following text the subscripts are changed to one simple, global index given by the formula . In such a way, the geometric mapping can be written as

(3.1)

Where is the coordinates of control point . Similarly we can write the displacement interpolation

(3.2)

Likewise,

(3.3)

With the kinematic mapping in hand, the basic finite element formulations can be applied on the NURBS-solid by following the procedure of the traditional isoparametric finite element [10].

IV.3.2 Small Strain Formulation

The balance equation for a static boundary value problem is expressed as

in domain

, on essential (i.e. displacement) boundary

, on natural (i.e. force) boundary

where is stress tensor, b is the body force, is the boundary traction, and n is the normal vector of natural boundary. As a standard procedure in FEM, the governing equation is cast into the weak form

where is any admissible variation of deformation, and stands for the symmetric part of the gradient of . For small deformation, the difference between reference configuration and current configuration is neglected, and measure of the strain and other quantities with respect to current and reference configuration is treated identical. If the material is linear elastic, the stress vector is represented as

(3.4)

where D is the material or constitutive matrix ,is the strain vector defined by . Using the interpolated displacement, the strain inside an element is computed as

(3.5)

In equation (3.5), ,, is the displacement in x, y, z direction, respectively. With the displacement interpolation given in (3.2), the strain-displacement matrix is written as, where

where. The comma followed by a coordinate denotes the derivative with respect to that coordinate. The derivative relative the spatial coordinate is computed with the aid of the chain rule. Starting from

where u,v,w is the parameters of the isoparametric mapping, and multiplied by , we get

The Matrix J is known as the Jacobian of the geometric mapping. In terms of the shape function defining the coordinate transformation, we have

Finally, the stiffness matrix for the linear small deformation formulation is

where are the weights for Gauss quadracture and is the Gauss point in the parametric domain, and is the order of gauss quadratures. is the determinant of the Jacobian matrix of the reference geometric mapping.

(3.6)

The external forces vector on a control points is determined as

whereis the global index of a control point.

Finally, after assembling all the elements equations into the entail domain, we obtained a system of linear equations in matrix form as

where f_g is the assembled global external force vector applied on all of control points with dimension of , d_g is the global vector of displacements of all control points with dimension of , and K_g is global stiffness matrix with dimension by .

After solve these linear equations, the displacements of control points, d_g, is obtained.

By updating the position of control points, we get the final shape after deformation. The stress and strain can be computed using equations (3.4) and (3.5).

Remark 3.2

We are at the point to discuss the similarity and difference of NURBS-FEM with the standard FEM. Like the finite element method, the displacement interpolation is constructed with the aid of a parametric domain. An underlying coordinate patch exists, with which the body is discretized into a collection of elements. However, the shape functions do not possess the Kronnecker delta property, therefore, the control points are generalized coordinates, not directly the nodal displacement as in FEM. Nevertheless, the NURBS-Fem shape functions have the partition of unity property, as the standard finite element methods and meshfree methods. This property ensures that any linear displacement can be interpolated exactly. Essentially, it guarantees the satisfaction of the patch test, which requires that any constant strain field be recovered exactly. An example of patch test is shown in Figure 4, where a constant strain field produced by applying a linear displacement on the boundary control points is computed exactly

Figure 4. Patch test: A constant strain field obtained by prescribing boundary control points according to a linear displacement field.

IV.3.3 Large Deformation Formulation

The small strain theory is inadequate for muscle mechanics because muscles in its normal function undergo large deformation, and because muscle’s stress-strain relation is highly nonlinear. In this project, a fully nonlinear formulation of NURBS-FEM is used for muscle analysis. The formulations are derived based on the large deformation elasticity theory. By background, the fundamental measure of deformation in large deformation is the deformation gradient

The strain of the body at a material point can be described by the right Cauchy-Green deformation tensor , or alternatively by the Lagrangian strain . The constitutive equation for a nonlinear elastic solid is typically given by an equation , where is the 2^nd Piola Kirchhoff stress which is related to the true (Cauchy) stress by , where is the determinant of the deformation gradient (). If the material is hyperelastic, the constitutive equation is specified with a strain energy density , such that .

(3.7)

The balance equation, i.e. the strong form of the boundary value problem, is the same as that of the small deformation. The weak form of the balance equations written relative to the current configuration is given as

Whereis the current mass density of the material. An equivalent form, defined in the reference configuration, is expressed as

(3.8)

where is the reference domain, and is the reference density, .

Using the geometric representations (3.1) and (3.3), the deformation gradient is computed according to

where is the referential gradient of the shape function, computed again using the chain rule.

The finite element equation can be established, equivalently, using either (3.7) or (3.8). Typically the equation is derived from (3.7). In the case, the finite element equation retains the form

The stain-displacement matrix B is identical to that in small deformation formulation. However, the derivative is taken relative to the deformed coordinated, which is unknown. For this reason, the equation is typically solved using iteration solver such as the Newton-Raphson method, in which a series of linearized equations are solved. The displacement increment at i-th iteration is computed from

The derivative in the bracket gives the so-called tangent stiffness matrix . In nonlinear elasticity the tangent stiffness consists of two terms,

where the first term is the material tangent , in which is the matrix form of the material elasticity tensor tangent in the current configuration [15]. The second term, , defines a tangent term arising from the non-linear strain-displacement relation and is often called the geometric stiffness matrix. Thus, we obtain a form of the finite deformation problem which is identical to that of the small deformation problem, with a geometric stiffness term added to the stiffness matrix. Here, we assume the external force f is not a function of displacement d, such that its derivative with respect to d vanishes. More detailed description of the stress, strain, and tangent moduli will be given in the section 5 with a constitutive model of muscle.

IV.4 Examples of Deformable Body Modeling

To demonstrate the qualitative behavior of NURBS-FEM, some basic NURBS primitives are deformed under various boundary conditions, as shown in Figures 5-7. In all the examples, we use only first order (linear) shape function for extended parameter (). This is equivalent to shrinking the boundary surface to its centerline. In general, the order of shape function and number of extended control points can be increased so that the interior interpolation is smoother. Numerically, this will increase the accuracy of finite element computation due to increases “mesh” resolution.

Figure 5 shows a sphere under applied stretch and pinching at its poles. The surface is described by 26 control points originally. By adding only one center point, a solid sphere obtained with totally 27 control points, and 81 DOF.

Figure 6 shows a torus under out-of-plane shear. Initially, there are 100 non-repeated control points for describing the torus surface. Ten additional control points are separated, which is initially repeated, from the origin control points for detaching the two end surfaces. And another 11 inside control points are added by extending the torus surface to solid. Thus, there are totally 121 control points for the torus solid and 363 DOF.

Figure 7 shows a cylinder under stretch. The cylinder surface is described by 40 control points. By shrinking the cylinder to the center line, another 5 control points are added for the solid representation. Thus, totally 45 control points describe the cylinder solid, with 135 DOF.

IV.5 Skeletal Muscle

Muscle tissue has a highly complex material behavior: it is active, nonlinear, incompressible, anisotropic, and hyperelastic [8, 18]. The microstructure of muscle is determined by the arrangement of the contractile elements (muscle fibers) and passive background tissue within a muscle [20].

IV.5.1 Constitutive Model

In this project muscle is modeled as a hyperelastic solid. The strain energy function is assumed to be the sum of two parts

where the first part is the strain energy associated with the passive ground substance, also known as the matrix material. This part is typically modeled as isotropic. The second part represents the active and passive muscle fiber strain energy. The ground substance consists of connective tissue, water, etc. In this report, it is modeled as a neo-Hookean material with an energy form

whereis the first principal invariant of C, are the material constants, and is the determinant of the deformation gradient. The constant may be best understood as a penalty parameter for incompressibility: nearly incompressibility can be approximately modeled with a large value of .

The muscle fiber strain energy is assumed a form

where is the active strain energy of the muscle fiber and is the passive strain energy of the muscle fiber due to stretch. The active strain energy is a function of the muscle fiber stretch and the muscle activation level , where and N is the fiber direction in the reference configuration. The first derivative of with respect to is defined as

and

With these terms the well known nonlinear stress-strain relationship of muscle fibers [4,5] (figure 8) can be modeled in the continuum level. and are the coefficients for the passive property of muscle fiber.

The second Piola Kirchhoff stress in a hyperelastic material is related to the strain by

where

And is the second order identity tensor, and

The Cauchy (true) stress related to S through takes the form

The material elasticity tensor is defined as. For hyperelastic material specified by a strain energy function W,

. In component form,

Specializing to defined muscle material model, we have

Where

The contribution from fiber energy gives

where

The spatial form of the material tensor, needed for computing the material stiffness matrix, is related to the elasticity tensor through the relation

, in component form,

which can be readily computed. The spatial tangent matrix is obtain by transforming the forth order tensor to a second order matrix using Voigt notation. The details are omitted.

IV.5.2 Fiber Representation

The fiber arrangement of skeletal muscles can be classified into several categories [11] : parallel-fibred muscles(muscle fibers oriented in parallel to the muscle line-of-action), pennate-fibred muscles (all muscle fibers oriented at the same angel relative to the muscle line-of-action), fusiform muscles (long, fusiform like muscle with fibers attach to the two ends), triangular muscles (fibers radiate form a narrow attachment at one end to a broad attachment at the other end) and others.

In this report, we use constant direction vector N in the undeformed configuration to model parallel-fibred muscles and unipennate-fibred muscles which have only one distinct fiber direction. In the work of Ng-Thow-Hing [17], B-spline solid models were built to capture muscle architectural details of internal fiber arrangements by carefully fitting the isocurves (two parameter are held constant while the third varies) of B-spline solid to the real fiber directions in actual muscle specimens. This idea is also used in our work. For non-uniform fiber distribution, the fiber direction at any point is acquired by computing the tangent of an appropriate isocurve through that point. For example, if parameter u is chosen to vary while the other two parameter v and w are held constant for an isocurve, there is an analytical expression

Taking the derivative of the above equation with respect to u, the normalized tangent vector of this isocurve is expressed as

Since the isocurve is used to model a muscle fiber, is used as the muscle orientation. For the gauss integration of our finite element formulation, which is carried out in each parametric element, the fiber orientation at a gauss integration point

Annual Technical Report No. VSR-04.15

Virtual Soldier Research (VSR) Program

Center for Computer-Aided Design

The University of Iowa

116 Engineering Research Facility

Iowa Cit, IA 52242

Project

Digital Human Modeling and Virtual Reality for FCS

CONTRACT/PR NO. DAAE07-03-D-L003/0001

Chapter IX

Table of Contents