An Implicit Finite Element Method For Elastic Solids In .
An Implicit Finite Element Method for Elastic Solids in ContactGentaro Hirota, Susan Fisher, Andrei State, Chris Lee*, Henry FuchsDepartment of Computer Science, University of North Carolina at Chapel Hill{hirota sfisher andrei fuchs}@cs.unc.edu*University of Colorado Health Sciences Centerchris@chs.uchsc.eduAbstractThis work focuses on the simulation of mechanicalcontact between nonlinearly elastic objects such as thecomponents of the human body. The computation of thereaction forces that act on the contact surfaces (contactforces) is the key for designing a reliable contacthandling algorithm. In traditional methods, contactforces are often defined as discontinuous functions ofdeformation, which leads to poor convergencecharacteristics. This problem becomes especially seriousin areas with complicated self-contact such as skin folds.We introduce a novel penalty finite elementformulation based on the concept of material depth, thedistance between a particle inside an object and theobject’s boundary. By linearly interpolating precomputed material depths at node points, contact forcescan be analytically integrated over contact surfaceswithout raising computational cost. The continuityachieved by this formulation supports an efficient andreliable solution of the nonlinear system.This algorithm is implemented as part of our implicitfinite element program for static, quasistatic anddynamic analysis of nonlinear viscoelastic solids. Wedemonstrate its effectiveness on an animation showingrealistic effects such as folding skin and sliding contactsof tissues involved in knee flexion. The finite elementmodel of the leg and its internal structures was derivedfrom the Visible Human dataset.1. IntroductionWhen animal and human bodies move, they deformdue to mechanical contact between components such asskin, muscles or bones. As body posture changes, organspush and slide against each other, changing the shape ofthe body. As a joint bends, the skin surface around it maystretch and fold, creating complicated geometry.Simulation of such phenomena would provideautomatic methods to generate the deformation. Inanimation, this capability could help create believablelooking details of organic bodies, eliminating timeconsuming manual intervention [7]. It could also benefittraining physicians and applications such as medicalimage registration, for example between pre-operativelyacquired CT or MRI datasets and intra-operativeultrasound or X-ray imagery. Surgical simulation oftenrequires procedure-specific postures, which can bederived by deforming generic models such as the “VisibleHuman” dataset [36].A short version of this paper was presented as [18].2. Previous workWhen simulating deformations of elastic objects, oneof the major challenges is avoiding penetration ofdeformable structures. In the following, we describe twomajor approaches in dealing with this problem.2.1. Kinematic approachesMost deformation techniques employed in computeranimation use kinematic approaches. Their majoradvantage is interactive performance due to the relativelysmall computational cost. Two examples are free-formdeformation (FFD) [34] and “skinning” or skeletonsubspace deformation (SSD) [23], which is the smoothblending of multiple rigid transformations. FFD and SSDbelong to a group of algorithms that employ “spacedeformation” [5], which can be viewed as a 3Dtransformation. One can also deform objects by directlymoving the control points of surfaces [11].In these methods, the impenetrability constraint issatisfied by heuristic techniques, often requiringextensive user interaction to produce the desired effects.
Space deformation can be applied to the human bodywithout causing penetration between the organs [7], butsliding effects between organs cannot be obtained withthis method.Many commercial software packages allow animatorsto embed formulae to express specific needs fordeformations [1]. These “deformers” can be written insuch a way that penetration between objects isminimized, but only for specific, limited scenarios.By using pose space deformation [23], one canpartially automate the process. For example, a user can“teach” the system to avoid penetration of the skinaround an elbow by directly adjusting the skin geometry.Henceforth,thesystemautomatically reducespenetration. The drawback of the method is that it ishighly inflexible since the user must instruct the systemhow to handle every new contact scenario.2.2. Simulation of physical lawsTo overcome the drawbacks of kinematic methods,many techniques employ the notion of force and energy[14, 21]. Wilhelms et. al. simulate a sliding skin layer byrelaxation of a spring mesh [39]. The relaxation schemedoes not account for buckling, hence realistic foldingdoes not occur.Accurate physical simulation, which has been studiedin computational mechanics, can provide a powerful toolfor automatically generating realistic deformations. Theyhave been traditionally very expensive in terms ofcomputation time but are becoming increasinglyaffordable with the continually growing performance ofcomputer hardware. The basic approach for contactproblems is to detect penetration between objects andcompute an appropriate response that eliminates,minimizes or reduces penetration. (This is alsocharacteristic of our approach.)Graphics researchers have demonstrated animationsof elastic bodies in contact [2,29,37,40]. As mentioned,the contact problem (i.e., avoiding penetration) has beenextensively studied in the engineering community [19].Solution methods for the contact problem can becategorized by their approach to satisfying theimpenetrability constraint. We now discuss the two majorapproaches.2.2.1. Hard constraint. Most algorithms utilizing hardconstraints use the concept of slave nodes and mastersurfaces to define constraints. Consider two collidingobjects. Nodal points on the surface mesh of one objectare designated as slave nodes, whereas master surfacesare defined from the surface mesh of the second object. Ifa slave node penetrates a master surface, a constraint tokeep the node on the surface is created. A set ofconstraints is formed by all of the penetrating slavenodes. With each iteration, new surface geometry isobtained by solving this constrained optimizationproblem, and the set of constraints is updatedaccordingly. The process repeats until the constraint setbecomes stable.This method suffers from two major drawbacks, thefirst being that it requires frequent constraint updates.Frequent constraint updates occur when highlytessellated surfaces are in sliding contact. As soon as aslave node travels across one master surface to another,the existing constraint becomes invalid, and a new onemust be created. Thus the lifetime of each constraint inthe optimization process is very short.The second drawback is known as the lockingproblem: the surface becomes artificially stiff due to anexcessive number of constraints. Hallquist et al. handlethis problem by partitioning the surface of each objectinto master and slave regions [15]. However, in thepresence of evolving self-contact, it becomes extremelydifficult to maintain the distinction between these twotypes of surface elements.Another method uses heuristic rules to find masterslave pairs based on the history of movements, andcomputes the exact time of collision, [8,16]. The twocolliding objects are treated symmetrically. In thismethod, even if the collision times are computed, somepenetrations are tolerated for several time steps beforethey are eliminated. To completely prevent penetrations,the exact birth and death times of impenetrabilityconstraints for slave-master pairs must be computed atevery step. Such computation is prohibitively expensive.Belytschko et al. [6] proposed the “splitting pinball”method, which uses repulsive forces between hierarchicalbounding spheres around individual surface elements.Although interference checks between spheres are veryefficient, the algorithm’s applicability to complex contactis not clear. The fine details of a surface such as foldingskin would require that the bounding spheres tional cost.Ideally, if two objects are in contact, they share asingle contact surface. However, because it is virtuallyimpossible for two independently discretized surfaces tohave common surface geometry, the methods discussedabove cannot prevent small penetrations or gaps. Thefrequency of constraint updates is also directly related tothe resolution of surface discretization, and can thereforelead to very high computational cost.
2.2.2. Penalty methods. By definition, penalty methodsallow small amounts of penetration to occur. To resolvepenetration, penalty forces proportional to the depth ofpenetration are calculated [3,9,30]. Unlike constraintmethods, a slave node is allowed to travel withoutkeeping track of all master surfaces on its path. Thus, theeffect of the resolution of surface discretization oncomputational cost is not as dramatic as with constraintmethods.2.2.3. Gap function computation. Most conventionalmethods (including constraint methods) seek projectionsto evaluate the gap function (i.e. the negative depth ofpenetration) and its derivative [12]. Because of thegeometric complexity of a projection search, thesemethods are not appropriate for handling complicatedboundary surfaces.Most methods (except for the pinball method)examine penetration at slave nodes only. Because of their“point sampling” nature, the contact force applied on aslave node becomes a discontinuous function ofdeformation (i.e. of node movements), which oftencauses a convergence problem.In section 5, we introduce a robust method to computea continuous gap function for even very complex surfacegeometry. We also show that, based on our gap function,contact penalty forces can be analytically integrated ascontinuous functions.3. Static analysisDue to their low mass density and fairly low viscosity,biological tissues tend to rapidly converge to a final statewhen subjected to external forces. For example, it is hardto flex a finger or change facial expression quicklyenough to be able to observe viscosity or inertia effectssuch as creep and oscillations. In fact, for manyapplications (including animation), the user is onlyinterested in the (static) equilibrium shapes of flexibletissues for a given posture or other specified constraints.In animation applications, dynamic postures andconstraints can be used to steer the animation. Staticanalysis is adequate for these applications, and, hencemethods optimized for the static problem must bedeveloped.In static analysis, the geometry of an elastic objectdepends solely on the forces applied to the object. Therelationship between geometry and forces is described bya differential equation defined on the continuous domainof the elastic object.We assume the resting (undeformed) shape of theobject is known, and elect it as a reference configuration.The current (deformed) shape is referred to as the currentconfiguration.At the equilibrium state, an elastic object in contactsatisfies the following equation:(div T f ) δu dV p n δu da 0 δu(1)ΩΓIn this equation, Ω is the interior of the object, Γ is itsboundary, T is the first Piola-Kirchhoff stress tensor,div is the divergence operator, f is the density of bodyforces (such as gravity), u is the displacement ofparticles, δu is an arbitrary variation of u , dV is thedifferential volume in the reference configuration, p isthe pressure on the contact surface, n is the normal of thecontact surface, da is the differential area in the currentconfiguration, and p n is the surface traction force on thecontact. Since we assume frictionless contact, the tractionforce is normal to the surface. p n is integrated over theboundary. (Actually, its value is non-zero only on thecontact area.)3.1. DiscretizationBecause of its complex boundary conditions andnonlinearity, Eqn. (1) cannot be solved analytically.Instead, it is discretized using a finite element method,and an approximate solution is sought [22].We use tetrahedral elements for the interior andtriangular elements for the boundary of objects. Thetriangular elements are chosen to be a subset of the sidesof the tetrahedral elements.The displacements of particles (internal materialpoints) are obtained by linearly interpolatingdisplacements at nodes. (The interpolation functions arecalled shape functions.) Elastic forces at nodes arecomputed by substituting the virtual displacement uwith the corresponding shape functions.The derivatives for all forces must also be computed toconstruct a stiffness matrix, which is crucial for Newtoniteration (described in subsection 3.2).The details of the contact force computation areexplained in section 5.As a result of discretization, we obtain a nonlinearequation of the form:R (u ) 0(2)Here, u represents the displacement vector(u1,1, u1,2 , u1,3 , . , un,1 , un,2 , un,3 ) , where (ui,1 , ui ,2 , ui,3 )denotes the displacement of the ith node and n is thenumber of nodes. We can impose boundary conditions fordisplacement by assigning fixed values to thecomponents of u.
3.2. Solution of the nonlinear systemBy solving Eqn. (2), we can obtain a new shape of theelastic object as the displacement vector u . There arethree factors that contribute to the nonlinearity of thissystem: Finite deformation: To handle finite deformation (asopposed to infinitesimal deformation), the forcedisplacement relationship must be described bynonlinear equations (geometric nonlinearity). Nonlinear material: Realistic elastic materials are allnonlinear. The choice of materials is explained insubsection 3.3. Collision and contact of objects: Collision andcontact are events that introduce additionalnonlinearity into the system.On top of the nonlinearity, the large size of the systemincreases the complexity; a typical finite element analysisof our interest produces a system with tens of thousandsof variables. We have developed a robust and efficientalgorithm to solve large nonlinear systems.3.2.1. Selecting a search direction. Nearly all solutionmethods start with an initial guess and proceed in a givensearch direction in a step-by-step manner. Finding theproper direction in the high-dimensional search space iscritical for the algorithm’s efficiency. Several differentmethods exist to determine the best search direction.The maximum gradient descent method chooses thedirection of the forces (i.e. the negative of the residualR (u ) ) for each step. This method turns out to be veryslow because it takes many steps for a local force topropagate through the entire mesh. Furthermore it isimpossible for this method to predict rapid force changescaused by the deformation of objects.On the other hand, the Newton method uses thederivative of the forces (i.e. the stiffness matrix), whichprovides information about how the forces vary as afunction of deformation. Each Newton step consists ofcomputation of the residual, computation of the stiffnessmatrix, and solution of a linear system. The processcontinues until the residual drops below a giventolerance. We chose this method due to both its speedand reliability.3.2.2. Avoiding illegal steps. Two questions that remainto be answered are: how to obtain the initial guess andhow far to proceed in a selected search direction.If there are no displacement boundary conditions, theobject is deformed only by external forces. In this case,an obvious choice of the initial guess is the initial shape,i.e. u 0.But if displacement boundary conditions are specified,some components of u are externally given. In suchcases u initially contains zeros for free components andfinal values for constrained components. Thecorresponding mesh may contain a tetrahedral elementwhose orientation is reversed, resulting in an illegalconfiguration. For this reason, constrained componentsmust be gradually incremented, and as soon as an illegalconfiguration is detected, the step size must be reduced(adaptive incremental loading). Each incremental stepperforms a Newton iteration and the solution is cascadedinto the next step. After the second step, the initial guessis computed by extrapolating the previous two solutions(two-point predictor).Given an initial guess, we must determine how far toproceed in a given search direction. If R (u ) is smooth,and the initial guess is close to the solution, a full steptowards the solution of the linear system can safely betaken. However, the nonlinearity of the system oftencauses a full step to lead to divergence. A full Newtonstep can also bring the mesh to an illegal configuration.Our method checks if the full Newton step is legal and ifthe residual decreases. If both of these are true, the step istaken. If not, the step size is halved until the twoconditions are satisfied. The scaling value of the step sizeis called a damping factor. This line search strategygreatly improves the robustness of our method.The resulting algorithm can be summarized as threenested loops:LOOP1: Adaptive Incremental Loading2-Point PredictionLOOP2: Newton IterationLinear System ConstructionLinear System SolutionLOOP3: Line Search3.2.3. Linear system solution. Newton iteration dependson the stable solution of linear systems. The accuracy ofthe solution is traded for speed. Since the degree ofnonlinearity of the equation is high, the residual is notgreatly reduced by a single Newton step. Consequently,computing an exact solution for each linear system doesnot improve the rate of convergence. For this reason, aniterative method is better suited for the linear systemsolution than are direct methods.The linear systems constructed in our finite elementmethod are symmetric, sparse, and usually positivedefinite. Around a bifurcation point, however, the matrixof the system is sometimes not only indefinite but alsonearly singular. The (bi)conjugate gradient method tendsto behave in an erratic manner. We found that animplementation of the Generalized Minimum Residualmethod (GMRES) with diagonal preconditioning [33] isboth efficient and stable. The iteration is terminated
when either the residual reduces to one tenth or apredetermined iteration limit is reached. This strategykeeps the computation time for solving the linear systemsrelatively low without slowing down the convergence ofthe Newton iteration.We use another technique to prevent failure in thelinear system solver. Each node is assigned a “viscosity”proportional to the sizes of the surrounding elements. Awell-chosen value for viscosity gives the algorithm thetendency to pick an energy-minimizing solution at abifurcation point, i.e. to prefer a stable equilibrium point.Also, this method allows starting with an obviously illconditioned initial configuration such as an object placedin mid-air, which will fall until it hits the ground.3.3. Material modelsThe property of an elastic material is defined as therelationship between stress and strain (Constitutive Law).In Eqn. (1), the relationship provides a formula thatassociates the stress tensor T with the displacement u .A few quantities must be introduced to describe therelationship: the deformation gradient F u Id ( Iddenotes the identity matrix), the right Cauchy-Greenstrain tensor C F T F , the invariants of C defined asI1 tr (C ), I 2 12 {tr (C ) 2 tr (C 2 )}, I 3 det( C ), and finallythe stored energy function ψ (F ) . The stress tensor T isgiven as T ψ ( F ) F , thus ψ (F ) actually determinesthe material property.The properties of organic tissues are being activelystudied in biomechanics and several models have beenproposed based on stress-strain data obtained from invivo and in vitro experiments. However, due to thelimitations of measurement technology, those modelshave not been rigorously validated [27].We have implemented compressible variations ofthree different material models. The first model is theMooney-Rivlin material [10,28]ψ (
Instead, it is discretized using a finite element method, and an approximate solution is sought [22]. We use tetrahedral elements for the interior and triangular elements for the boundary of objects. The triangular elements are chosen to be a subset of the sides of the tetrahedral elements.
Finite Element Method Partial Differential Equations arise in the mathematical modelling of many engineering problems Analytical solution or exact solution is very complicated Alternative: Numerical Solution – Finite element method, finite difference method, finite volume method, boundary element method, discrete element method, etc. 9
Finite element analysis DNV GL AS 1.7 Finite element types All calculation methods described in this class guideline are based on linear finite element analysis of three dimensional structural models. The general types of finite elements to be used in the finite element analysis are given in Table 2. Table 2 Types of finite element Type of .
1 Overview of Finite Element Method 3 1.1 Basic Concept 3 1.2 Historical Background 3 1.3 General Applicability of the Method 7 1.4 Engineering Applications of the Finite Element Method 10 1.5 General Description of the Finite Element Method 10 1.6 Comparison of Finite Element Method with Other Methods of Analysis
The Finite Element Method: Linear Static and Dynamic Finite Element Analysis by T. J. R. Hughes, Dover Publications, 2000 The Finite Element Method Vol. 2 Solid Mechanics by O.C. Zienkiewicz and R.L. Taylor, Oxford : Butterworth Heinemann, 2000 Institute of Structural Engineering Method of Finite Elements II 2
The Generalized Finite Element Method (GFEM) presented in this paper combines and extends the best features of the finite element method with the help of meshless formulations based on the Partition of Unity Method. Although an input finite element mesh is used by the pro- . Probl
describe earlier published methods of creating interpolating implicit surfaces. 2.1 Implicit Surfaces An implicit surface is defined by an implicit function, a continuous scalar-valued function over the domainR3. The implicit surface of such a function is the locus of points at which the function takes on the value zero.
4.3. Finite element method: formulation The finite element method is a Ritz method in that it approximates the weak formulation of the PDE in a finite-dimensional trial and test (Galerkin) space of the form V h:" ' h V0, W h:" V0, (4.6) where ' h is a ane o set satisfying the essential BC of (4.1) and V0 h is a finite-dimensional .
nite element method for elliptic boundary value problems in the displacement formulation, and refer the readers to The p-version of the Finite Element Method and Mixed Finite Element Methods for the theory of the p-version of the nite element method and the theory of mixed nite element methods. This chapter is organized as follows.
