UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A physics-based model for wrinkling skin Harrison, Darcy 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2015_november_harrison_darcy.pdf [ 7.5MB ]
JSON: 24-1.0166733.json
JSON-LD: 24-1.0166733-ld.json
RDF/XML (Pretty): 24-1.0166733-rdf.xml
RDF/JSON: 24-1.0166733-rdf.json
Turtle: 24-1.0166733-turtle.txt
N-Triples: 24-1.0166733-rdf-ntriples.txt
Original Record: 24-1.0166733-source.json
Full Text

Full Text

A Physics-Based Model For Wrinkling SkinbyDarcy HarrisonB.Sc., University of Manitoba, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Computer Science)The University of British Columbia(Vancouver)September 2015© Darcy Harrison, 2015AbstractWrinkling of human skin significantly affects the realism of computer generatedcharacters. Wrinkles convey emotion and expression, provide clues of age andhealth, and indicate interaction between the skin and external objects. Wrinklingis caused by compression: an elastic material buckles out-of-plane in order to pre-serve length and volume. Human skin buckles in a distinctive pattern, character-ized by sharp valleys with rounded peaks. Many techniques used in visual effectsrequire artists to directly produce wrinkles through sculpting or painted displace-ment maps, while automated techniques are generally designed for adding detailto coarse, cloth-like simulations which are usually not consistent with human skin.The layered structure of skin, and the properties of each layer are critical to pro-ducing the buckling patterns observed in real life.In this work a simulation of wrinkling skin is developed that is physicallybased, while also simple enough for use in computer graphics. A novel constitutivemodel suitable for large compressive strain is derived and applied to a three-layeredmodel of skin, with a thin shell outermost layer (stratum corneum), and volumet-ric dermis and hypodermis layers. Finally, we present a modified Newton schemeand linear finite elements for simulating equilibrium configurations of skin undercompressive strain.iiPrefaceThis dissertation is original, unpublished, independent work by the author, DarcyHarrison, under the guidance and supervision of Dr. Dinesh K. Pai.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.1 Skin Wrinkles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2 Physics of Buckling . . . . . . . . . . . . . . . . . . . . . . . . . 82.3 Physical Skin Models . . . . . . . . . . . . . . . . . . . . . . . . 92.4 Computer Graphics Models . . . . . . . . . . . . . . . . . . . . . 113 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133.1 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133.2 Stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17iv3.3 Quasistatics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183.4 Potential Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . 193.5 Constitutive Models . . . . . . . . . . . . . . . . . . . . . . . . . 204 Skin Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 244.1 Structural Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 244.1.1 Modeling the Stratum Corneum . . . . . . . . . . . . . . 254.1.2 Modeling the Dermis and Hypodermis . . . . . . . . . . . 264.2 Constitutive Model . . . . . . . . . . . . . . . . . . . . . . . . . 274.2.1 A Novel Stretching Model . . . . . . . . . . . . . . . . . 274.2.2 Bending Model . . . . . . . . . . . . . . . . . . . . . . . 284.3 Prestrain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 295 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315.1 Linear FEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325.1.1 Stratum Corneum . . . . . . . . . . . . . . . . . . . . . . 325.1.2 Dermis & Hypodermis . . . . . . . . . . . . . . . . . . . 365.2 Optimization Framework . . . . . . . . . . . . . . . . . . . . . . 385.3 Elastic Strain Energy Derivatives . . . . . . . . . . . . . . . . . . 405.4 Indefiniteness . . . . . . . . . . . . . . . . . . . . . . . . . . . . 436 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 486.1 2D results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 486.1.1 Experiment 1 . . . . . . . . . . . . . . . . . . . . . . . . 496.1.2 Experiment 2 . . . . . . . . . . . . . . . . . . . . . . . . 516.1.3 Experiment 3 . . . . . . . . . . . . . . . . . . . . . . . . 516.1.4 Experiment 4 . . . . . . . . . . . . . . . . . . . . . . . . 556.1.5 Experiment 5 . . . . . . . . . . . . . . . . . . . . . . . . 556.2 3D Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 577 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 627.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 627.2 Limitations and Future Work . . . . . . . . . . . . . . . . . . . . 657.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67vBibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68A Supporting Materials . . . . . . . . . . . . . . . . . . . . . . . . . . 74A.1 Membrane & Volume Energy Derivatives . . . . . . . . . . . . . 74A.1.1 Gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . 74A.1.2 Hessian . . . . . . . . . . . . . . . . . . . . . . . . . . . 75A.2 Bending Energy Derivatives . . . . . . . . . . . . . . . . . . . . 78A.2.1 Gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . 78A.2.2 Hessian . . . . . . . . . . . . . . . . . . . . . . . . . . . 81viList of TablesTable 6.1 Base parameters used in 2D experiments. . . . . . . . . . . . . 49Table 6.2 Parameters specific to Experiment 3. . . . . . . . . . . . . . . 51Table 6.3 Parameters specific to Experiment 4 . . . . . . . . . . . . . . . 55Table 6.4 Parameters used in our 3D results . . . . . . . . . . . . . . . . 60viiList of FiguresFigure 1.1 Structure of the skin, from Fro¨hlich et al. [17] . . . . . . . . . 3Figure 2.1 The wrinkle profile given in Hatzis [24] as one might see in thefurrows caused by contracting the muscles of the forehead. . . 8Figure 3.1 Elastic potential energy as a function of deformation gradi-ent, for a 1D strand utilizing Saint Venant-Kirchoff (ST. V-K),Corotational linear, and Hencky models. . . . . . . . . . . . . 21Figure 3.2 Force response as a function of deformation gradient, for a 1Dstrand utilizing ST. V-K, Corotational linear, and Hencky models. 22Figure 4.1 Comparison of constitutive models linear in Hencky strain vs.Bazˇant’s approximation. . . . . . . . . . . . . . . . . . . . . 28Figure 5.1 A triangular element of the stratum corneum . . . . . . . . . 33Figure 5.2 A bending element of the stratum corneum . . . . . . . . . . 35Figure 5.3 Bending energy density as a function of dihedral angle for lin-ear, 2sin(θ2)and 2tan(θ2)models discussed by Garg et al. [18]. 35Figure 5.4 A tetrahedral element of the dermis or hypodermis . . . . . . 37Figure 5.5 Plot of Equation 5.46 as a function of endpoint position x . . . 44Figure 5.6 The osculating paraboloid of Equation 5.46 about the pointx = (0.11,−0.68). In (a) the shape is a hyperbolic paraboloid,indicating an indefinite Hessian. In (b) the modified Hessian isan elliptic paraboloid because it has no negative eigenvalues. 45viiiFigure 5.7 Plot of Equation 5.46 with two strand segments, forced to buckle.The red ‘x’-mark indicates a saddle point, confirming the non-convexity of the system. . . . . . . . . . . . . . . . . . . . . 46Figure 6.1 A 3cm long, 2D cross-section of our 3-layered skin model un-der uniform compression. . . . . . . . . . . . . . . . . . . . 50Figure 6.2 A 3cm long, cross-section of our 3-layered skin model using acorotational linear constitutive mode for the stratum corneum(SC) membrane, dermis and hypodermis. . . . . . . . . . . . 52Figure 6.3 A 3cm long, cross-section of our 3-layered skin model using alinear bending mode for the SC. . . . . . . . . . . . . . . . . 53Figure 6.4 A 3cm long, cross-section of a 2-layered skin model of the SCand underlying substrate. . . . . . . . . . . . . . . . . . . . . 54Figure 6.5 A 3cm long, cross-section of a 3-layered skin model of the SC,living epidermis and dermis. . . . . . . . . . . . . . . . . . . 56Figure 6.6 A 3cm long, cross-section of our 3-layered skin model withvaried SC thickness, compressed to 86% its original size. . . . 58Figure 6.6 A 3cm long, cross-section of our 3-layered skin model withvaried SC thickness, compressed to 86% its original size. . . . 59Figure 6.6 A 3cm long, cross-section of our 3-layered skin model withvaried SC thickness, compressed to 86% its original size. . . . 60ixGlossarySC stratum corneumST. V-K Saint Venant-KirchoffPDE partial differential equationSVD singular value decompositionxAcknowledgmentsMy deepest thanks and appreciation to Dr. Dinesh K. Pai for his excellent supervi-sion. Thanks also to Dr. Uri Ascher and Ye Fan for their feedback.xiDedicationTo my love, Christine Schultz.xiiChapter 1IntroductionSkin is the largest and most visible organ of the human body, so it is fair to con-sider a virtual character as being defined by the shape and motion of its skin. Oneof the primary goals of the visual effects industry is to create believable, virtualhumans for film, broadcast, visualization, advertising, etc. An artist equipped witha realistic virtual human can create visuals that defy reality, explore design choicesthat relate to human clothing or equipment, or even place their virtual character insituations where it is impossible or unsafe to use real actors. All of these scenariosrequire an accurate model for not only the visual aspects of skin (i.e., how to syn-thesize an appropriate texture and lighting behavior), but also how skin reacts tothe motion of the underlying body or interactions with external objects. Aestheticsresearchers describe the uncanny valley, a psychological phenomenon wherein au-dience comfort drops significantly as a virtual character approaches, but does notquite reach a healthy, natural likeness. The subtle motions and wrinkling of skinmay be an important, remaining barrier towards passing through this valley.Wrinkles in human skin cover a wide range of phenomena grouped together byfact that they are all buckling of the skin due to compression. This includes expres-sive wrinkles due to the compression caused by the contraction of muscles, as wellas interaction wrinkles caused by some external compressive force applied to theskin. Permanent wrinkles develop as the result of physical changes to the skin’scomposition due to age, as well as permanent expressive wrinkles that develop dueto repeated, long-term buckling in a consistent pattern.1To understand why and how skin buckles, we need an understanding of thestructure of skin. Human skin is a non-homogeneous, elastic material, generallystructured into layers. A wrinkle forms as the thinner, stiffer outer layers bucklerelative to the thicker, softer underlying layers with a characteristic wavelength. Inhumans, we can readily observe wrinkles on the scale of ∼ 2mm and higher. Thefact that there are more than two layers leads to a distinct, hierarchical bucklingpattern which differs significantly from a simple sinusoidal profile.For animation purposes, muscles and bones are often considered most impor-tant to model, with the motion of skin derived procedurally. Much effort is spentmanually tweaking weights, poses, and sculpts to fix the artifacts produced bytreating skin as a function of the bones and muscle. Additionally, the fine scale ofwrinkles is often a computational burden and is usually handled in a non-physicalmanner. Artists paint wrinkle textures, or sculpt wrinkle poses for repetitive mo-tions and facial expressions. Extremely fine scale wrinkling, or wrinkling due toexternal interactions, are often handled with cloth simulation technology produc-ing wrinkles unlike those found in real skin. A physically-based simulation ofhuman skin is necessary to increase realism. Realistic models are studied in thebiomechanical and cosmetics literature; however, there is no clear consensus onwhat aspects of the structure and behavior of skin are truly critical for producinglife-like, believable wrinkles. Finding the minimal set of mechanical layers, a rel-atively simple constitutive model, and an efficient discretization that lead to theexpected buckling profile will offer significant improvement over the skin modelscurrently used in computer animation.In this thesis, we present a minimal physically based skin model that bucklesunder compressive load in a manner that is consistent with the shape of real skinwrinkles. We aim to identify the subset of structural phenomena without whichthe buckling patterns are inconsistent, then validate these choices by developing asimulation of compressive load applied to small patches of skin. By developing askin model that buckles in a manner visually consistent with skin we can increaserealism in artificial characters, especially those experiencing highly variable exter-nal interactions and non-repetitive motions. There are many phenomena related toskin that would also benefit from a physically based skin model. For example, sen-sation could be simulated since many nerve endings in the skin respond to pressure2Figure 1.1: Structure of the skin, from Fro¨hlich et al. [17]and other applied stresses.With the goal of generating plausible wrinkle shapes in mind, we must furtherconsider the structure and properties of skin. Maurel et al. [38] describe skin as anon-homogeneous, anisotropic, non-linear, viscoelastic, multi-component materialranging in thickness, by location, from 0.2mm to 6mm. Skin covers most of thebody, in some areas tightly coupled to the underlying muscles and loosely drapingbones in others (usually at joints like the elbow or knuckles). The skin is usuallyconsidered as consisting of three major layers: the outermost epidermis, the dermis,and the underlying hypodermis. This is illustrated in Figure 1.1.The epidermis is further split into a thin outermost layer of dead cells called thestratum corneum (SC) that make a stiff, nearly inextensible sheet of approximately20µm [14]. The remainder of the epidermis is the viable epidermis consisting ofliving basal cells which undergo mitosis and migrate to the stratum corneum [35].Additionally there are pores, glands, etc.The dermis is a dense network of collagen and elastin fibers embedded in agel-like medium. The majority (75%) is collagen with much less elastin (4%) andthe remaining portion consisting of the gel of proteins and water, blood and lymphvessels, nerve endings, hair follicles and glands. The collagen fibers are composedinto fibrils of diameter 100 nm, which are further arranged into bundles of diam-eter from 1 to 40µm. Flynn [13] writes that at rest, collagen fibers are randomlyoriented and unextended and do not significantly contribute to the behavior of theskin. As skin is stretched however, the collagen bundles uncrimp and align withthe direction of loading until eventually dominating the elastic behavior as they3become fully extended. The elastin fibers do not form fibrils or bundles, and in-stead form a scattered network between the collagen. Daly [10] notes that elastin(and the ground substance it is embedded in) dominates the behavior of skin undercompression and minor tension. There are also a variety of other structures suchas blood vessels, hair follicles and glands. The dermis is generally subdivided intothe thinner, less fiber-dense papillary dermis, and the reticular dermis.The hypodermis is a subcutaneous fatty tissue providing the loose, flexibleconnection with to the underlying muscle and soft tissues [38]. Some authors em-phasize the viscoelastic effects of the hypodermis [e.g., 14]. There is relativelylittle concrete information about the hypodermis in the biomechanics literature.Classical theory predicts that a thin, elastic sheet bound to a thicker substratewill buckle into a sinusoidal profile of specific wavelength, with amplitude de-pendent on compressive strain. This is relevant to skin, since the SC layer of theepidermis is very thin and stiff. The layers under the SC are much thicker, andgenerally monotonically less stiff according to most authors [e.g., 14, 29]. The factthat there are several distinct layers underneath lead to skin buckling in a uniquemanner. When human skin buckles, the shape it assumes has broad, rounded peakswith narrow troughs, unlike the sinusoidal profile predicted in simple, two-layeredmaterial.Having established that many techniques currently in use in computer graphicsdo not capture the buckling profile seen in real skin, the question remains as towhich set of skin properties are necessary in order to improve simulation realism.Many authors disagree on which layers of the skin must be modeled, lumpingall, or several of the layers into some sort of average. Certainly the literature agreesthat a single or two layer model is insufficient for achieving multi-scale bucklingeffects. Assuming one chooses a three layer model, there is disagreement on whichlayers the three should be. The SC appears to be a common choice, while the twoor more layers beneath it are chosen inconsistently. The number of layers, and therelative thickness of the chosen layers will determine which methods are appropri-ate for discretizing the skin. These choices additionally impact the computationalefficiency of a simulation.Once a layer combination is chosen, the constitutive models used to capturethe material properties of each layer must be considered. While linear elasticity4is definitely not applicable, it would be wise to find a relatively simple non-linearmodel that behaves similarly to the compressive response of skin. The constitutivemodel will determine if each layer is a homogenous material, or if it contains adistribution of fibers (i.e., anisotropy) or other sub-components. In fact, there aremany state of the art constitutive models that explicitly model the orientation andextension of a distribution of fibers within the dermis. While it seems that thesefibers are mostly related to the behavior of skin under tension, they may contributeto the buckling behavior as well. Strain rate dependent effects (i.e., viscoelasticity)are also the responsibility of the constitutive model.Plasticity, the permanent deformation of a material, is also a relevant phe-nomenon for skin, especially as it relates to the development of permanent wrin-kles. It is desirable to include a mechanism that approximates plastic deformation.Li et al. [33] describe a reduced coordinates, membrane model of skin, gener-ating a realistic elastic response to deformation that causes skin to slide and stretchbut discarding any motion out of plane (such as wrinkling). Their model is partic-ularly effective for computer graphics as it utilizes an efficient linear finite elementformulation, and implicitly keeps the skin constrained to lay on the underlyingbody. We extend their technique with a physical, layered model of the skin. Givenan initial solution to the coarse, membrane-only motion of the skin, we remove anyresidual compressive strain by solving for buckling out-of-plane.A three layered model of skin will be proposed in Chapter 4, consisting of thestratum corneum, dermis and hypodermis. The living epidermis is eliminated (orincluded into the dermis layer), while the papillary and reticular dermis layers arecombined into one. Owing to the small relative thickness of the SC, it is modeledusing the theory of thin shells which reduces the volume into a flat, 2D surface withmembrane (i.e., stretching) and bending components. The dermis and hypodermisare modeled volumetrically since shear strain, and strain through the thickness isimportant to retain. A novel constitutive model, quadratic in an approximation tothe classic logarithmic strain tensor, is derived for use in the membrane componentof the SC as well as the dermis and hypodermis. This constitutive model is wellequipped to handle large compressive strain compared to many similarly simplemodels. Anisotropy and viscoelastic effects are omitted. The concept of prestrainis included which allows for the skin model to include the natural state of tension5that is found in real skin. Additionally, prestrain enables a simple approximationfor plasticity and permanent wrinkles.In Chapter 5 the proposed skin model is implemented via linear finite elements,and included in a quasistatic simulation that utilizes a modified Newton minimiza-tion algorithm. Linear finite elements partition the skin surface (i.e., the SC) intotriangles, with the dermis and hypodermis partitioned into tetrahedra. The lay-ers are fully coupled, meaning the separate elements utilize shared vertices at theboundaries. The constitutive model is implemented for the linear finite elementsusing a simple 1 point quadrature rule, while the bending forces are implementedusing non-linear hinge elements. A bending energy appropriate for large bend-ing strain is associated to the hinge elements. A quasistatic simulation is chosensince the inertial effects in skin are minimal, being driven primarily by boundaryconditions. An energy minimization problem results, which is addressed with amodified Newton method. Unfortunately the problem is not inherently convex, soconvergence of a Newton method is not guaranteed. The fixes for local indefi-niteness that are commonly used in computer graphics are considered, but shownto be insufficient when applied naively so an alternative modification strategy isprovided.The remainder of this document is organized into chapters which progress to-wards a minimal, physically based model of skin wrinkling. Chapter 2 will surveythe related works of previous authors, including the cause and categorization ofwrinkles, the theory of buckling, and current state of the art skin models from thebiomechanics literature. Chapter 3 will cover the principles of elasticity, and sur-vey some non-linear constitutive models. The chosen skin model is described inChapter 4 while the implementation within a quasistatic simulation is described inChapter 5. Several experiments will be outlined in Chapter 6, whose results willillustrate the impact of the decisions made in the skin model. Finally, we finish upwith Chapter 7 wherein the applicability of the skin model is discussed as well aslimitations and future work.6Chapter 2Related Work2.1 Skin WrinklesThere is significant difficulty even describing and classifying the variety of wrin-kles present in human skin. Pie´rard et al. [44] notethe literature appears quite confusing because there is no consensusabout the definition of such terms as wrinkle, rhytide, crinkle, crease,groove, line and furrow.They are primarily interested in the wrinkle phenomena that emerge due to the ag-ing process as the elastin and collagen in the skin is chemically modified. Skinmicrorelief, the typical patterns of intersecting lines and regularly arranged hairfollicle and duct openings [44], are an important aspect of the skin’s appearance.Over time this microrelief will develop into age-related wrinkles. In contrast wewould like to study the buckling behavior of skin due to contraction of the under-lying muscles, orientation of the joints, or external interactions. We are interestedthough, in the permanent, expressive wrinkle caused by repeated muscle contrac-tions leading to stiffening of hypodermis adjacent to the muscle.Magnenat-Thalmann et al. [36] separate wrinkles into the categories of expres-sive wrinkles (buckling or folding of the skin surface due to skin deformation,which disappear after removal of the deformation) and permanent, and aging re-lated wrinkles. They note that repeated expressive wrinkling at the same location7Figure 2.1: The wrinkle profile given in Hatzis [24] as one might see in thefurrows caused by contracting the muscles of the forehead.will lead to the development of permanent wrinkles in time (i.e., Type 3 wrinklesof Pie´rard et al. [44]).It is important to note that expressive wrinkles do not have a sinusoidal pro-file (unlike the results of Le´veˆque and Audoly [32], and the two layer model inMagnenat-Thalmann et al. [36]). The cross-sectional curve of real wrinkles haspeaks similar to a sine wave, but much sharper valleys [e.g., 24, 36]. Figure 2.1shows a wrinkle profile that is expected in the forehead region, or can be observedwhen significantly compressing the skin of the forearm. This is obviously over-simplified however, since Efimenko et al. [12] show that buckling skin would haveat least five nested buckling wavelengths present.2.2 Physics of BucklingCerda and Mahadevan [8] analyze the asymptotic buckling behavior of a thin,stretched, elastic sheet under uniaxial compression. They find the critical ingredi-ents to be an inextensible thin layer with bending stiffness B, an elastic foundationwith stiffness K, and an imposed compressive strain ε . The balance between bend-ing and foundational energies produces wrinkles with wavelength λ and amplitudeA proportional toλ ∝(BK) 14(2.1)A ∝ (ε)12 λ (2.2)They conclude that the geometric constraints (i.e., inextensibility) lead to the for-mation of wrinkles, the bending resistance of the sheet penalizes short wavelengths,and the foundation layer penalizes longer wavelengths.8Efimenko et al. [12] experiment with a sandwich of a thin (5nm) stiff skin withuniaxial prestrain of 30, 50, and 70% atop a substrate (0.5mm) of reduced stiff-ness. When the sample is removed from the specified tension conditions that leavethe surface skin at zero strain, the material buckles to maintain approximate inex-tension of the much stiffer skin. They observe competition between the bendingresistance of the skin, and the stretching of the substrate which suggests a fixedwavelength of buckling and leave the amplitude free. As the tension is further re-moved, the nonlinear effects of the substrate causes the amplitude to saturate andforce the material to buckle at a larger wavelength. They observe five distinct gen-erations of wrinkles from this two layer material which they argue is similar tohuman skin.Kuwazuru et al. [29] proposed a five-layered skin model for analysis of multi-stage buckling. They divide the skin into three main layers: the epidermis, dermis,and hypodermis. They further split the epidermis into the thin, extremely stiff stra-tum corneum and the remaining part is referred to as the viable epidermis. Thedermis is split into a softer papillary dermis on top, and a less flexible reticularlayer. They chose a simple, homogeneous, isotropic model for each layer sincethe purpose is to analyze the general buckling properties of a layered medium,as opposed to matching experimentally observed data. As shown by Cerda andMahadevan [8] a thin, stiff layer atop a softer foundation is the recipe for buck-ling hence in this model the stratum corneum buckles relative to the foundationof the viable epidermis. Under further compression, the viable epidermis (alongwith the stratum corneum) buckles relative to the papillary dermis followed by afinal buckling stage with the hypodermis as foundation. Their model confirms thehypothesis that each buckling stage occurs separately with increasing strain, andproduces wrinkles with increasing wavelength.2.3 Physical Skin ModelsMaurel et al. [38] describe the critical features that a physically based skin modelmust consider: non-homogeneity, anisotropy, non-linearity viscoelasticity, multi-ple layers/components.Many state-of-the-art constitutive models are reviewed in [35]. Some mod-9els are purely phenomenological, considering the skin as a homogeneous materialdescribed by a single relationship between strain and stress. Others consider thestructure and geometric configuration of the micro-structural elements of the skin.For example, they often include a distribution of fibers (i.e., collagen) within theskin and use a phenomenological model for these individual fibers. The commonfeature of the presented models is the emphasis given to the biphasic nature ofcollagen within the skin. The Gasser-Ogden-Holzapfel model, for example, com-bines a low-strain linear elastic portion attributed to elastin with an exponentialterm associated to the collagen [35]. It also includes a two direction anisotropyrelated to fiber orientation distributions. The other models presented have similarfeatures, using different techniques to model the collagen response, anisotropy, vis-coelasticity, plasticity, etc. The focus appears to be mostly on the behavior of skinunder tension, and does not mention any discretization of skin into layers and thusshould not be expected to be useful in a buckling simulation (without combinationor modification, that is).Flynn [13] provides a detailed survey of state of the art structural models specif-ically for the dermis. They specifically treat the collagen and elastin fiber distribu-tions, as well as the response of the gel matrix these fibers are embedded within.Many of the presented models have been used within simulations of wrinkling skin.Magnenat-Thalmann et al. [36] developed a 2D model featuring the epidermisand the dermis only that did not produce wrinkling with the expected shape. Theyproposed an additional three layer model of the stratum corneum, living epider-mis, and deep dermis which better matched experimental data. The correct shapeof wrinkle was only produced when they included a wavy interface between theepidermis and dermis.Le´veˆque and Audoly [32] investigate a similar model: a thin, stiff SC with softliving epidermis, and somewhat stiffer deep dermis. They were primarily interestedin a parameter study of how the elastic properties of the SC affected the wrinklingperiod and amplitude. Their results indicate that their model produced only a singlebuckling wavelength, so it is would not be effective for the a large compressivestrain simulation.Flynn and McCormack [14–16] developed a three-layer finite element modelfeaturing the SC, dermis and hypodermis modeled as neo-Hookean, orthotropic10visco-hyperelastic, and viscoelastic-Yeoh materials respectively. Their simulationproved consistent with the literature, and the predictive properties were shown to beuseful for matching experimental data from compression of the forearm in severaltest subjects. Their results show multiple buckling wavelengths and qualitativelymatch our desired wrinkle shape.2.4 Computer Graphics ModelsIn visual effects, skin is commonly simulated by wrapping a muscle and bonesystem with a tightly fitted, thin cloth layer. The shape change of the muscles andfat underneath cause the skin to stretch and slide, or buckle in regions experiencingcompression. The results are often unsatisfactory, however, and tend to look likecloth.A related technique, due to Li and Kry [34], Re´millard and Kry [46], involvescoupling one or more skin layers (as thin shells) to an underlying body using area-averaged constraints. These constraints are chosen such that they allow the thinshell to buckle independently of the body at wavelengths consistent with humanskin, as predicted by Cerda and Mahadevan [8], while also following the coarsemotion of the body. Their results are visually plausible, with Li and Kry [34] es-pecially producing multi-wavelength buckling by utilizing multiple, layered shells.The underlying body, and the constraints between the multiple thin shell surfacesare not derived from the actual structure of skin however.Li et al. [33] describe a technique for treating the skin as a membrane (i.e., ashell which experiences in-plane strain but does not resist buckling). Their methodimplicitly constrains the skin to lie on the surface (by way of a parametrization) ofthe underlying muscle and fat. The membrane will slide and stretch realistically,but does not produce any buckling in regions subject to compression.Other techniques include modeling a variety of wrinkled poses and using themas blendshapes [e.g., 26], or using strain maps to reveal (painted, or synthesized)wrinkle textures as bump or displacement maps [e.g., 23, 53, 54]. For arbitrarydeformation a single, static wrinkle texture is insufficient so often several wrinkleimages could be combined based on measurement of the strain tensor [e.g., 28].Many techniques are described in the literature for simulating nearly inextensi-11ble cloth [e.g., 1, 4, 7, 20]. Several authors describe techniques for adding detailedwrinkles to coarse cloth simulations. Several of these work by generating explicitwrinkle curves (automatically or from an artist) that are revealed based on the par-ticular strain configuration [e.g., 9, 31, 45, 48]. Others upsample a coarse simu-lation using high frequency wrinkles from another dataset [e.g., 27] or a looselyconstrained, higher resolution shell [e.g., 40]. In general, techniques designed foradding wrinkles to clothing do not produce wrinkle shapes or patterns consistentwith human skin.12Chapter 3BackgroundIn this section we provide a brief introduction and discussion of non-linear elastic-ity. For a more thorough treatment, refer to Bonet and Wood [6].An elastic medium is defined by a relationship between deformation and theforces which seek to eliminate it. Deformation is measured by a strain tensor whichencapsulates changes in length and angles between vectors in material space andthe equivalent vectors in the deformed physical space. The elastic material on oneside of an oriented, differential area exerts a force on the other side (and vice versa).In general, the orientation can be freely chosen so that the forces are characterizedby a stress tensor which converts an orientation into a traction (force per unit area).The traction at the material boundary is uniquely determined since its orientationis provided by the surface normal. The surface traction of the entire medium canbe related to the forces in the volume via the divergence theorem.3.1 KinematicsKinematics is the name associated to the tools which measure and describe motionwithout an consideration as to the cause. For elasticity we are interested in mea-suring deformation (i.e., change in shape) which is best illustrated with a 1D strandoriginally of length L that has been stretched (or compressed) to length l. The most13straightforward way to measure the deformation is the relative length changeε =lL−1. (3.1)An alternative which has nice properties when moving to higher dimensions isgiven byE =l2−L22L2=12((lL)2−1), (3.2)where the constant factor of one half is chosen so that the Taylor series expansionabout l = L matches Equation 3.1 to first order.Another useful strain measure which has the symmetric property that a defor-mation and its inverse have equivalent strain (up to sign) is given byH =∫ lLdll= log(lL). (3.3)This is the natural extension of the small strain, linear model to finite elasticity.The equivalent three dimensional measure is known as the Hencky strain tensor.The fact that these measures of deformation have different values for differentstretches illustrates that there is no definitive way to measure a deformation. Thismeasure of deformation is called a strain tensor.Consider a continuous object parametrized by position X in some reference (ormaterial) pose. The mapping from the reference pose to the current (or physical)pose is given by the map x(X). The pointwise quantity which is the analogue ofthe stretch lL in the above description, is the deformation gradientF =dxdX. (3.4)A small displacement dX in the material pose is transformed to physical spacedx via the deformation gradient asdx = FdX =dxdXdX . (3.5)A general measure of strain can be derived from observing two material space14perturbations dX1 and dX2 and the inner product of their corresponding physicalspace perturbations. The inner product will measure the lengths and angle betweenthese vectors in physical space.〈dx1|dx2〉= dxT1 dx2 (3.6)= (FdX1)T FdX2 (3.7)= dXT1 FT FdX2 (3.8)= dXT1 C dX2. (3.9)The change in length and angle is then given by〈dx1|dx2〉−〈dX1|dX2〉= dXT1 C dX2−dXT1 dX2 (3.10)= dXT1 (C− I)dX2 (3.11)and the relative change in length and angle by〈dx1|dx2〉−〈dX1|dX2〉〈dX1|dX2〉 =dXT1 (C− I)dX2dXT1 dX2. (3.12)Since the perturbations dX1 and dX2 are arbitrary, the strain measure E (calledthe Green or Langrange strain tensor) is completely determined byE =12(C− I) , (3.13)with Cauchy strain tensorC = FT F. (3.14)The constant factor of one half, as in Equation 3.2, is to match the first order be-havior of a strain measure which is linear in stretch.Analogous to Equation 3.1, a strain measure which is linear in stretch (calledthe Biot strain tensor) is given asB =√C− I. (3.15)15A critical, yet intuitive, property of any strain tensor is invariance under affinetransformation. This means that an arbitrary rotation and/or translation of the entireobject must not have any impact on the strain tensor; a rigid motion does not actu-ally lead to a change of shape, just a change of frame. This can be made explicitby decomposing the deformation gradient into the product of a rotation (possiblyincluding a reflection) R and a symmetric, positive-definite stretch tensor U via apolar decompositionF = RU (3.16)Since C = FT F = URT RU = U2, clearly the Green strain from Equation 3.13and Biot strain from Equation 3.15 can be written asE =12(C− I) = 12(U2− I) , (3.17)B =√C− I =U− I. (3.18)illustrating their invariance under rotations. The invariance under translations is as-sured by the nature of the deformation gradient which discards any constant trans-lation in the map x(X).There is an entire class of strain tensors, the Seth-Hill (also called Doyle-Ericksen) family, which includes the Green and Biot tensorsE(m) = 1m (Um− I) if m 6= 0lnU if m = 0. (3.19)The tensor E(0), analogous to Equation 3.3, is also known as Hencky, logarith-mic, or true strainH = E(0) = lnU =12lnC. (3.20)A convenient and efficient approximation to the Hencky strain tensor due toBazˇant [5] is given byHˇ =12(U−U−1) . (3.21)This approximation will be used extensively in Chapter 4.163.2 StressThe stress in a deformed body is measured point-wise as the force exerted by oneside of an oriented, differential area, dA, on the other (and vice versa). If theresulting force on the area is ∆ f then the traction vector, t, is given byt(n) = limda→0∆ fda. (3.22)This leads to the definition of the Cauchy stress tensor, σ , which produces a tractionfrom the outward facing normal vector, n, in physical spacet(n) = σn. (3.23)The first Piola-Kirchoff stress tensor P, related to the Cauchy stress tensor byP = det(F)σF−T , (3.24)produces the same physical space traction but in terms of the outward facing, ma-terial space normal, N, ast(N) = PN. (3.25)Consider a body with volume Ω and surface ∂Ω experiencing an external bodyforce fext per unit volume. Due to elastic stress, the body experiences a traction tat the boundary. Translational equilibrium requires the total experienced forces bezero. ∫∂Ωt dA+∫Ωfext dV =∫∂ΩPN dA+∫Ωfext dV = 0. (3.26)By the divergence theorem, the traction term is converted to a volume integral∫∂ΩPN dA+∫Ωfext dV =∫Ω(∇ ·P+ fext) dV = 0. (3.27)The internal force fint per unit volume is defined asfint(X) = ∇ ·P(X). (3.28)The principle of virtual work dictates that static equilibrium will be obtained17when the variation of work (w) due to any virtual displacement (δu) is zero.δw = 0 =∫Ωδu · (∇ ·P+ fext) dV (3.29)=∫Ωδu · fext dV +∫Ωδu ·∇ ·P dV (3.30)=∫Ωδu · fext dV −∫ΩP :ddXδu dV +∫Ω∇ · (PTδu) dV (3.31)=∫Ωδu · fext dV −∫ΩP : δdudXdV +∫∂Ωδu ·PN dA (3.32)=∫Ωδu · fext dV −∫ΩP : δF dV +∫∂Ωδu · t dA (3.33)∴∫ΩP : δF dV =∫Ωδu · fext dV +∫∂Ωδu · t dA. (3.34)Note that the tensor contraction operator (A : B = tr(AT B)) has been used. Equa-tion 3.34 is useful as the weak form of the partial differential equation (PDE) thatwill solved via the finite element method in Chapter 5. It will also be expanded onin Section QuasistaticsIn a simulation including dynamics, the principle of virtual work would be aug-mented (via the principle of d’Alembert) to include the force of inertia. In skinhowever, the inertial effects are somewhat negligible due to low mass and highstiffness. In reality, the motion of the skin is largely determined by boundary condi-tions and the timescale of the strain waves which might propagate are significantlybelow the time discretization that we are interested in. This justifies a quasi-staticsimulation wherein at each time instant a static equilibrium is sought, as opposedto integrating the positions and velocities in time.183.4 Potential EnergyFor certain classes of elastic materials (i.e., hyperelastic) it turns out that the virtualwork due to internal forces is the negative variation of some function Ψ such thatδwint =∫ΩP : δF dV =−∫ΩδΨ dV. (3.35)Since δΨ = dΨdF : δF this identifies the first Piola-Kirchoff stress as the negativederivative of the elastic potential energy density functionP =−dΨdF. (3.36)In this sense, the definition of the elastic potential energy density function Ψ com-pletely defines the behavior of a hyperelastic material since it provides the relation-ship between strain and stress. Thus, we refer to such a function as the constitutivemodel for a hyperelastic material.The implications of this relationship to elastic potential energy density is thatthe stress is path-independent, i.e., a function of only the initial and current con-figurations. Furthermore, there are restrictions on the potential energy function: itmust be invariant to physical space rotations. An isotropic material (i.e., one whichhas a response unbiased by material space orientation) will have the further re-striction that the potential energy function must also be invariant to material spacerotations. Thus, for arbitrary rotation matrices R and QΨ(RFQT ) =Ψ(F). (3.37)A consequence of this is that given the singular value decomposition (SVD) of F asF =UΣV T it is required thatΨ(F) =Ψ(UT FV ) =Ψ(Σ). (3.38)Clearly such a constitutive model is actually a function of the principle strains, i.e.,the singular values of F . Many constitutive models are defined directly in termsof the principal strains, as opposed to utilizing other invariants of the strain tensorsvia the trace and determinant.193.5 Constitutive ModelsSeveral common constitutive models in the literature are simply quadratic functionsof some strain tensor A, in the formΨ(F) = µ‖A(F)‖2F +λ2tr2(A(F)). (3.39)Models defined this way yield stresses which are linear in the chosen strain tensor.The quantities µ and λ are the Lame´ coefficients, derived from Young’s modulusY (a measure of stretch resistance) and Poisson’s ratio ν (a measure of compress-ibility).µ =12Y(1+ν), (3.40)λ =ν(1−2ν)E(1+ν). (3.41)Young’s modulus and Poisson’s ratio for many materials are provided in the engi-neering literature.Models of this form include that of Saint Venant-Kirchoff (ST. V-K) (using theGreen strain tensor E from Equation 3.13), co-rotational linear elasticity (using thestretch tensor U from the Biot strain Equation 3.15), or the Hencky strain tensor Hfrom Equation 3.20.St. Venant-Kirchoff Ψ(F) = µ‖E‖2F +λ2tr2(E) (3.42)= µ∥∥12 (C− I)∥∥2F +λ2tr2(12 (C− I)).Co-rotational Linear Ψ(F) = µ‖B‖2F +λ2tr2(B) (3.43)= µ‖U− I‖2F +λ2tr2(U− I).Hencky Ψ(F) = µ‖H‖2F +λ2tr2(H) (3.44)= µ‖ ln(U)‖2F +λ2tr2(ln(U)).The ST. V-K model Equation 3.42 is computationally efficient, but suffers a200 U0W (U)1Ψ(U) = 14 (C− I)2Ψ(U) = (U− I)2Ψ(U) = ln2 (U)Figure 3.1: Elastic potential energy as a function of deformation gradient, fora 1D strand utilizing ST. V-K, Corotational linear, and Hencky models.fatal flaw in that under increasing compression it lacks a monotonic increase inassociated stress, violating Drucker’s first stability criterion[11]. In Figure 3.1 thisis indicated by lack of convexity as the material is infinitely compressed (i.e., goesto zero volume). In Figure 3.2 the same flaw is illustrated as a non-monotonicstress response as compression is increased. Since the behavior under compressionis of utmost importance in a buckling scenario, ST. V-K is not useful for wrinklingskin.The co-rotational linear elastic model Equation 3.43 is quite common in com-puter graphics, likely due to intuitive linear relationship between strain and stress.Computationally it requires a polar decomposition (or more generally an SVD)an operation with significant expense. The stress derivative dPdF can also be trickydue to the derivatives of the polar decomposition. Formulae are available in theliterature, however[39]. The main flaw of this model is that it has a finite stressresponse when infinitely compressed which, numerically speaking, will allow a210 U0ddU W (U)1Ψ(U) = 14(U2− I)2Ψ(U) = (U− I)2Ψ(U) = ln2 (U)Figure 3.2: Force response as a function of deformation gradient, for a 1Dstrand utilizing ST. V-K, Corotational linear, and Hencky models.volume element to invert (i.e., have negative volume, a decidedly non-physical be-havior). This is especially relevant for a model of buckling skin where we expectsignificant compression.The constitutive model of Equation 3.44 has many desirable features, such assymmetric responses to equivalent tension and compression, as well an infinitelystrong response as the material is subjected to an infinite compression. The latterproperty may be especially important since the buckling of skin implies the pres-ence of significant compressive strain. Furthermore, the second term of this model(i.e., λ2 tr2 (lnU)) is equivalent to λ2 ln2 (det(F)) which measures the relative changein volume of the material, enabling its use with nearly incompressible materials.The main downside of this model is that it is computationally expensive (evalu-ation of logU implemented via spectral decomposition or Taylor series expansion)and has complicated derivatives (however formulae are published in [43])There is another family of constitutive models that are commonly used to22model isotropic rubber, polymers or biological tissues. The polynomial hypere-lastic model due to Rivlin and Saunders [47] is a double infinite series of the formΨ(F) =∞∑i=0, j=0Ci j tr(C− I)i tr(C2− I) j (3.45)parametrized by constants Ci j, with C00 = 0. For specific truncations of these se-ries, and choices of constants we obtain a variety of common constitutive modelsdue to Mooney, Rivlin, Ogen, Yeoh, etc. As specified above, this model is de-signed for purely incompressible materials. It assumes that the pointwise, relativevolume (measured by J = det(F) =√det(C)) is unchanged during deformation.In order to extend this to the compressible regime we add another infinite serieswith coefficients Dk such thatΨ(F) =∞∑i=0, j=0Ci j tr(J−23 C− I)itr(J−43 C2− I) j+∞∑k=1Dk(J−1)2k. (3.46)This family is popular in the literature because there are several standardizedexperiments for determining the material constants of rubber-like materials. Par-ticular loading conditions are applied such that the relation between stress and de-formation becomes straightforward leading to discovery of best-fit coefficients fora specific truncation of the polynomial hyperelastic model. Ogden [42] describestypical deformations to incompressible rubbers, measuring uniaxial tension/com-pression, equibiaxial tension/compression, and pure shear.In practice a number of models are fairly simple truncations of the polynomialhyperelastic model:Neo-Hookean Ψ(F) =C10 tr(C− I), (3.47)Mooney-Rivlin Ψ(F) =C10 tr(C− I)+C01 tr(C2− I). (3.48)While the Neo-Hookean model is similar to the ones discussed previously, thehigher order truncations produce consitutive models with the freedom to bettermatch the multiphasic tensile response of rubber-like materials.23Chapter 4Skin ModelIn this chapter we choose a specific combination of skin layers to model, as wellas an appropriate consitutive model.4.1 Structural AnalysisGiven the constitutive options available for describing elastic materials, the ques-tion remains: which model has the breadth to capture the most prominent mechan-ical behaviors of buckling skin? Furthermore, is the geometric arrangement of skinmaterials important? In fact, the material properties of different skin layers are sig-nificantly different. The SC is critical to the buckling behavior of skin, but is oftenmodeled as a stiff, isotropic, Neo-Hookean solid [e.g., 14–16, 36], while the der-mis has received many treatments (surveyed in [35] and [13]) that emphasize theanisotropy, fiber distributions, and structural properties. The hypodermis, when in-cluded in the model [e.g., 14, 16], may have significant viscoelastic behaviors. Inearlier models, skin was treated as a single elastic membrane with no real thickness[e.g., 36] but most modern computational models agree that skin is best modeledas a collection of heterogeneous layers with differing constitutive properties [e.g.,14, 36].While some works utilize a two layer model [e.g., 8, 36], most recent publica-tions have identified three layers required for capturing the buckling characteristicsof skin. The SC, dermis, and hypodermis are the usual layers [e.g., 14–16, 32].24Some models ignore the hypodermis in favor of splitting the dermis into the pap-illary dermis and reticular dermis [e.g., 36]. The two layers discussed in [34] areconsidered along with their connection to an underlying body so it should be con-sidered a three layer model as well. At a minimum, the relationship of wavelengthto the relative stiffness of skin and substrate layers given by Cerda and Mahadevan[8] implies at least that there must be three layers to have the multi-scale bucklingobserved in human skin. A two layer model would only generate a single bucklingwavelength with a sinusoidal deflection profile.Following the lead of Flynn and McCormack [14] we consider a three layermodel of the SC, the dermis, and the hypodermis. While their model accounts forthe anisotropy of the dermis as well as the viscoelasticity of the hypodermis, weignore these effects and choose an isotropic hyperelastic model for both layers. Forcomputer graphics purposes the rate dependent effects of the hypodermis should benegligible. The anisotropy of the dermis is mostly due to collagen fiber networksand their behavior under tension so may therefore be ignored when consideringcompressive buckling models.4.1.1 Modeling the Stratum CorneumThe thickness of the SC (≈ 20µm) relative to the other two dimensions (which,depending on the patch, are on the order of centimeters to meters) suggest thata thin shell model is applicable. The high relative stiffness further supports useof a thin shell which is based on the assumption of small strain and isometricdeformations[2]. The isotropic Neo-Hookean model used in [14] and [32] is alsoconsistent with thin shell assumptions in this scenario. With consideration to thediscretization of the various layers, a full volumetric model would require ex-tremely small elements; a common rule of thumb is to require at least two elementsthrough the thickness. Thus volumetric elements in the SC would be much smallerthan those used in dermis and hypodermis which is known to negatively effect thequality of the finite element method. In computer graphics, thin shell models arequite commonly used for cloth [e.g., 7, 21], thin walled structures [e.g., 21], andhas recently been used in wrinkling models [e.g., 34, 46].An important implication of a thin shell model, is that the elastic potential25energy of the shell is decomposed into the sum of a membrane potential energy(i.e., a stretching term, derived from the 2D limit of the full volumetric equationsgoverning elasticity) and a bending potential energy.The SC is generally modeled as nearly incompressible, with Poisson ratio 0.495[19, 25]. The constitutive model chosen for this layer must have a componentwhich measures volume change.4.1.2 Modeling the Dermis and HypodermisThe dermis and hypodermis are significantly thicker than the SC (≈ 50×) and so afull volumetric treatment is required in order to capture the significant shearing andcompression through the thickness of these layers. The dermis layer is ≈ 1.2mmthick while the hypodermis layer is ≈ 1.5mm thick, as per Flynn and McCormack[14], though the numbers vary in the literature and by the location on the body.The structure of the dermis suggests an anisotropic material, and many physi-cally based models are parametrized on the directionality of collagen fibers (via asingle, pair, or distribution of fiber orientations). For a model of buckling however,we expect the compressive behavior of the dermis to dominate and therefore thecollagen orientation should not be as critical. Thus we choose an isotropic consti-tutive model for the dermis layer. We expect the dermis to have a Young’s modulus≈ 1000× less than the SC as described in Cerda and Mahadevan [8], but there areconflicting ratios in the literature (e.g., 10× [36] and 100× [29]).In anticipation of the time independence of our proposed quasistatic simulation,we discard the viscoelastic effects (i.e., dependence on strain rate) of the hypoder-mis from Flynn and McCormack [14]. The Yeoh model they propose is similar to apolynomial expansion of a log-strain model with parameters as the coefficients ofeach term. In the next section we will use a direct approximation of the log-straininstead.The SC is tightly coupled to the dermis, and the dermis to the hypodermis; noseparation of the layers is allowed. Fascia connects the hypodermis to the underly-ing muscle with varying rigidity (and it doesn’t at all in some areas where the skindrapes the bone). This coupling to the underlying muscle is accomplished in ourmodel via linear springs with varying stiffness, up to an infinite stiffness for a rigid26attachment.4.2 Constitutive Model4.2.1 A Novel Stretching ModelAs described in Section 3.5, the Hencky strain tensor from Equation 3.20 is wellsuited to situations where significant compression is expected. A simple and in-tuitive constitutive model is one that is linear in this strain tensor (i.e., Equa-tion 3.44). The computational challenges of the Hencky strain tensor (i.e., thematrix logarithm is expensive, numerically unstable with a naive implementation,and the derivatives, especially the Hessian, are difficult to compute) are a limita-tion. Bazˇant [5] provides an efficient approximation of the Hencky strain tensorHˇ =12(U−U−1) , (4.1)that maintains its critical properties: an infinite magnitude when compressed tozero volume, and a symmetry between compression and tension. Since the traceof the Hencky strain is merely the determinant of the deformation gradient (i.e., itmeasures the relative volume change) we should maintain that term exactly whileusing Bazˇant’s approximation for the part measuring length and angle change (i.e.,the µ part). Luckily the derivatives of the volumetric part of the linear Henckystress model are easy to compute so the efficiency gains of the approximation aremaintained.Ψ(F) = µ‖Hˇ‖2F +λ2tr2(H) (4.2)= µ‖Hˇ‖2F +λ2ln2(detU) (4.3)=µ4‖U−U−1‖2F +λ2ln2(detU) (4.4)=µ4tr(C+C−1−2I)+ λ8ln2(detC). (4.5)Here the Cauchy strain tensor C = FT F is trivial to compute, so the identity C =U2is used to avoid computing the polar decomposition of F = RU . Bazˇant suggests270 U0W (U)1W (U) = ln2 (U)W (U) = 14(U−U−1)2Figure 4.1: Comparison of constitutive models linear in Hencky strain vs.Bazˇant’s approximation.this approximation will be very accurate for 2× compression or stretch, which isalso supported by Figure 4.1. Outside of this range the behavior is similar enoughfor graphics purposes. In fact, this approximation fixes one of the numeric issues ofthe log-stress model wherein it is only convex for stretches less than e. Equation 4.5is strictly convex in this scenario.We use this constitutive model to completely describe the elastic behavior ofthe dermis and hypodermis, as well as the membrane component of the SC.4.2.2 Bending ModelThe bending model of the SC is provided by the bending potential energy describedin Audoly and Pomeau [2]. While the Fo¨ppl-von Karman equations are derivedassuming a linear relationship between strain and stress, our log-based relationship(i.e., Equation 4.5) is consistent with a linear model under small strain (i.e., the firstorder Taylor expansion is the same). Audoly and Pomeau [2] show that the bending28energy must only be accurate for nearly inextensible deformations of the shell’scenter surface, thus implying a small membrane strain and ensuring consistencywith our non-linear membrane model. When the membrane strain is not small,the bending terms are dominated and therefore of little consequence. Our bendingenergy density is a function of the local principal curvatures of the surface and isnot dependent on the membrane strain.The thin shell model eliminates explicit strain through the shell thickness, re-placing it with an elastic response to bending. Only the so-called “center surface”of a thin shell is modeled, lumping the material thickness. When the shell is bent,the material on either side of the center experiences stretching and compressionrespectively. The amount of stretching/compression experienced is proportional tothe curvature of the surface[2]. Thus, the energy associated to bending is a func-tion of the mean curvature H of the surface in physical space, relative to the meancurvature H¯ in material space.Ψb =∫∂ΩD(H− H¯) dA (4.6)The bending stiffness D is derived from the shell thickness h and usual materialproperties given in the literature: Young’s modulus Y and Poisson’s ratio νD =Y h312(1−ν2) (4.7)The nature of the bending energy, and a modification for large bending strainwill be given in Chapter 5.4.3 PrestrainFor a given physical space deformation, not all strains correspond to a valid config-uration in material space. It may be convenient to have the strain-free configurationbe different from the material space configuration. This can be accomplished withthe idea of a prestrain which encodes some residual strain which would exist evenif the physical space were to match the material space configuration. The bend-ing term of the SC is particularly easy, merely substituting a modified curvatureH0 instead of the material space curvature H¯. For the stretching terms, we right-29multiply the deformation gradient by a stored prestrain deformation gradient F0.The deformation gradient becomesF =dxdXF0. (4.8)It may be useful to consider F0 as the mapping from a new, strain-free space X0, tothe material space. As such we can writeF0 =dXdX0, (4.9)though the quantity F0 may not correspond to any realizable configuration in thestrain-free space. Interpreting it as a derivative is useful conceptually but not re-quired.In vivo the skin is never in a strain-free configuration, instead it exists in a nat-ural state of tension [38]. Langer’s lines [30] indicate the directions and presenceof prestrained collagen and elastin fibers in the skin. The exact strains are not in theliterature according to Flynn and McCormack [14], but the effect is that the criticalcompression which leads to buckling is delayed relative to a simulation that doesnot account for prestrain.Prestrain, especially of the SC curvature, is an effective way to model perma-nent wrinkles. An artist could locate a buckled configuration that they wish to con-vert into a permanent wrinkle, then using painted weights, capture some portion ofthe current configuration as the stored prestrain values. A similar application is apseudo-plasticity effect (i.e., developing permanent wrinkles from repeated, char-acteristic motions). For some characteristic animation wherever the quasistaticsimulation finds a solution where the final strain exceeds some critical value, theprestrain can be increased by some fraction.30Chapter 5ImplementationInitially the COMSOL®multi-physics package was used to implement the qua-sistatic simulation and skin model proposed in Chapter 4. Unfortunately, the ex-treme stiffness of the nearly-incompressible stratum corneum prevented the com-mercial package from producing usable solutions regardless of the discretization(linear or quadratic triangular as well as bilinear hexahedral elements were tried).A similar approach using their dynamics solvers was also found unsuitable. There-fore a custom buckling simulation was implemented as a Matlab®program using acustom Inexact Newton optimization algorithm. The skin is discretized using lin-ear finite elements, a thin shell for the stratum corneum coupled to the volumetricelements of the dermis and hypodermis. The Hencky-like constitutive model fromEquation 4.5 was used for modeling the elastic response of the skin, and a non-linear, discrete-hinge bending model [7, 21, 50] was used for the shell elements.Formulae for the first and second derivatives of the potential energy function wereinitially produced by the symbolic differentiation toolbox of Matlab®, but theywere found to have significant numerical difficulties. A derivation by hand was per-formed and tested against finite difference techniques. This symbolic expressionwas helpful for obtaining insights into removing indefiniteness from the Hessianmatrix. It is also expected to perform significantly faster than computer generatedversions [50].315.1 Linear FEMAssuming no external forces or tractions (and discarding inertial effects), Equa-tion 3.34 yields the weak form of the equation we will use the finite element methodto solve ∫ΩP : δF dV = 0. (5.1)As mentioned in Section 3.4, a hyperelastic material has an elastic strain (i.e., po-tential) energy function W (F) from which the first Piola-Kirchoff stress tensorP = dΨdF is derived. For such a material Equation 5.1 indicates we are locating acritical point of the total elastic strain energy functionW=∫ΩΨ(F) dV. (5.2)In fact, this is called the “minimum total potential energy principle”. In the sectionsto follow, we discuss the discretization of the skin which allows us to evaluate thisintegral and solve for the displacements which minimize Equation Stratum CorneumDiscretizing the thin shell model described in Section 4.1 using linear, triangularelements is appropriate for computer graphics and has been well studied in theliterature (mostly in relation to cloth). The point-wise quantity of interest for aLagrangian elasticity simulation is the physical space position corresponding toeach fixed material point. For a finite element implementation, these values arecomputed at the nodes (i.e., triangle mesh vertices) and interpolated within theelements by a shape function. A linear triangular element makes use of a linearshape function (i.e., linear or barycentric interpolation). The element is illustratedin Figure 5.1. For all points ξ ∈ {R2+ | ξ1 + ξ2 ≤ 1} within an element T withnodal positions x0,x1,x2 ∈ R3x(ξ ) =x0 x1 x2N(ξ ), (5.3)32x0x1x2Figure 5.1: A triangular element of the stratum corneumwith linear shape functionN(ξ ) =1−ξ1−ξ2ξ1ξ2 . (5.4)StretchFor the stretching term of the stratum corneum, we require the calculation of thedeformation gradient within the elements. The linear shape function ensures thephysical position derivative Dx is constantDx =∂x(ξ )∂ξ=x0 x1 x2−1 −11 00 1 . (5.5)The material position is interpolated similarly, producing the material positionderivative DXDX =dX(ξ )dξ=X0 X1 X2−1 −11 00 1 . (5.6)Therefore, the deformation gradient F for the element isF =∂x∂X=∂x∂ξdξdX= Dx(DX)†, (5.7)33where † denotes the pseudo-inverse.For isotropic elasticity the material space rotation is irrelevant, so Li et al. [33]suggest that the pseudo-inverse of DX should be computed via QR decomposition(i.e., DX† = R−1QT ) and then modified to remove the rotation part Q. This yieldsa 3 × 2 reduced deformation gradientF3×2 = DxR−1. (5.8)As a consequence of the reduced dimensionality of F , the Cauchy strain tensorC = FT F is a 2×2 tensor, significantly reducing the numerical complexity of themembrane energy calculations for the shell.As mentioned in Section 4.2.1 the elastic strain energy function Ψm(F) forthe SC membrane is defined by Equation 4.5. Since the deformation gradient Fis constant within an element, the total energy for an element T (with area A inmaterial space) of the surface ∂Ω isWm,T = A[µ4tr(C+C−1−2I)+ λ8ln2(detC)]︸ ︷︷ ︸Ψm(F)(5.9)and the total energy for the entire surface is found by summing over all elementsWm = ∑T∈∂ΩWm,T . (5.10)BendingAs discussed in Section 4.2.2 the bending energy of the SC is derived from the meancurvature [18, 21]. Obviously, a triangular element has a flat interior so the only lo-cations with non-zero curvature are the edges and vertices. The hinge element[21],popular in computer graphics, defines bending elements at the shared edge be-tween two adjacent triangles. The dihedral angle (i.e., the angle between trianglenormals) encodes the integrated mean curvature across the shared edge. This isconverted to a pointwise quantity and re-integrated over the area of support (i.e.,one third of each adjacent triangle) [22]. The bending element for two adjacent34x0x1x2 x3e0e1 e2n1 n2(a) 3D viewx0/x1x2 x3n1 n2θ(b) Side view (i.e., along e0)Figure 5.2: A bending element of the stratum corneumθW¯ (θ)Ψb(θ) = θ 2Ψb(θ) = 4sin2(θ2)Ψb(θ) = 4tan2(θ2)(a) Bending energy densityθW¯ (θ)Ψb(θ) = θ 2Ψb(θ) = 4sin2(θ2)Ψb(θ) = 4tan2(θ2)(b) Bending stressFigure 5.3: Bending energy density as a function of dihedral angle for linear,2sin(θ2)and 2tan(θ2)models discussed by Garg et al. [18].triangles T1 and T2 is illustrated in Figure 5.2. The shared edge e0, identifies thebending element. The discrete mean curvature is derived from the dihedral angle θbetween normals n1 and n2 across the edge. The quantity θ¯ refers to the dihedralangle in the reference pose, and the bending stiffness D is given by Equation 4.7.Several formulae for discrete bending energy are provided in the literature,though they are all compatible within first order to an energy proposed in Grinspunet al. [21]Ψb =D2(θ − θ¯)2 , (5.11)while Garg et al. [18] utilizeΨb =D2(2sinθ2−2sin θ¯2)2(5.12)35and Tamstorf and Grinspun [50] proposeΨb =D2(2tanθ2−2tan θ¯2)2. (5.13)We choose Equation 5.13 because it produces infinite strain energy as the bend-ing angle approaches 180◦, while Equation 5.11 is contrasted in Chapter 6. The2sin(θ2)based bending energy(in the same manner as the ST. V-K stretch model)does not yield monotonically increasing stress for increasing bending strain andis therefore inappropriate. The different discrete bending energies are plotted inFigure 5.3.The total elastic bending energy for an element ei is given byWb,ei =3‖e¯i‖2A¯i[D2(2tanθ2−2tan θ¯2)2]︸ ︷︷ ︸Ψb, (5.14)where A¯i is the sum of the areas of the two adjacent triangles, and ‖e¯i‖ is the lengthof the edge, both in material space. The total energy for the entire surface is foundby summing over the set of all edge elements E = {ei} in the triangle meshWb = ∑ei∈EWb,ei . (5.15)5.1.2 Dermis & HypodermisSimilar to the shell model of the stratum corneum, the dermis and hypodermisare discretized using linear, tetrahedral elements. The physical space positioncorresponding to each fixed material point is computed at the nodes (i.e., tetra-hedral mesh vertices) and interpolated within the elements by a similar, linear,shape function. The element is illustrated in Figure 5.4. For all points ξ ∈ {R3+ |36x0x1x2x3Figure 5.4: A tetrahedral element of the dermis or hypodermisξ1 +ξ2 +ξ3 ≤ 1} within an element T with nodal positions x0,x1,x2,x3 ∈ R3x(ξ ) =x0 x1 x2 x3N(ξ ), (5.16)with linear shape functionN(ξ ) =1−ξ1−ξ2−ξ3ξ1ξ2ξ3 . (5.17)As in the stratum corneum, we require the calculation of the deformation gra-dient within the elements. The linear shape function ensures the physical positionderivative Dx is constantDx =∂x(ξ )∂ξ=x0 x1 x2 x3−1 −1 −11 0 00 1 00 0 1 . (5.18)The material position is interpolated similarly, producing the material position37derivative DXDX =dX(ξ )dξ=X0 X1 X2 X3−1 −1 −11 0 00 1 00 0 1 . (5.19)Therefore, the deformation gradient F for the element is the 3×3 matrixF =∂x∂X=∂x∂ξdξdX= Dx(DX)−1. (5.20)As in the SC, the deformation gradient F is constant within an element so thetotal energy for an element T (with volume V in material space) of the body Ω isWv,T =V[µ4tr(C+C−1−2I)+ λ8ln2(detC)]︸ ︷︷ ︸Ψv(5.21)and the total energy for the entire volume is found by summing over all tetrahedrain the skin volume ΩWv = ∑T∈ΩWv,T . (5.22)5.2 Optimization FrameworkThe total energy of the skin is found by summing the total membrane and bendingenergies of the stratum corneum, with the total elastic energy of the dermis andhypodermisW= ∑T∈∂ΩWm,T + ∑ei∈EWb,ei + ∑V∈ΩWv,T . (5.23)The nonlinear minimization problem with objective function described in Equa-tion 5.23 can be solved in a variety of ways such as steepest descent, via a Quasi-Newton scheme (e.g., BFGS) or the full Newton method. The expectation of supe-rior convergence (quadratic within some region of the solution) associated with theNewton method is offset by the difficulty that is associated to calculating the matrixof second derivatives of the energy function. This Hessian matrix may be difficult38to derive by hand, or numerically unstable if derived via automated means (auto-matic differentiation or symbolic computations) [50]. Furthermore, the Hessianmatrix may be indefinite when not near a solution thus leading to search directionsthat may not consistently make progress towards a local minimum. In fact, for thepurposes of simulating buckling skin we can be assured that the Hessian matrix islocally indefinite due to the bifurcating nature of the problem [50]. This impliesthat we should expect an indefinite Hessian quite often. Results indicate that pro-jecting the system Hessian matrix onto a nearby positive definite version(as donein [39, 51]) will improve the conditioning of the system in addition to ensuringconsistent progress towards a minimum.Note that in the following sections, xk, pk, ∇W(xk), etc.will refer to a singlevector in which the per-vertex values have been stacked (i.e., a 3n vector, with nbeing the number of vertices in the discretization).A general optimization framework requires an initial guess for the solution x0and generates a sequence of iteratesxk+1 = xk +αk pk, (5.24)converging to the configuration with minimum elastic strain energy. At each itera-tion, the search direction pk and the step size αk can be calculated by a variety ofstrategies. Effectively, the search direction must be a descent direction (i.e., mustlead to a decrease in energy for some small, unknown step) which is assured by anegative inner product with the objective function gradient〈pk|∇W(xk)〉< 0 (5.25)and the step size should be chosen to approximately minimize the objective func-tion along the search directionW(xk +αk pk)<W(xk). (5.26)In general there are further conditions (i.e., the Wolfe conditions) [41] on the stepsize that must be met in order to assure convergence. For arbitrary constants k1 and39k2 such that 0 < k1 < k2 < 1, there is the condition of sufficient decreaseW(xk +αk pk)≤W(xk)+ k1αk∇W(xk)T pk, (5.27)and curvature condition∇W(xk +αk pk)T pk ≥ k2∇W(xk)T pk. (5.28)For a Newton’s method based minimization algorithm, the search direction pkis found by solving∇2W(xk)pk =−∇W(xk), (5.29)requiring the solution of a set of linear equations obtained from the Hessian matrixof second derivatives ∇2W(xk). Once pk is found, a suitable step size αk is chosenby a simple, backtracking line search from Nocedal and Wright [41]. The curva-ture condition of Equation 5.28 is implicitly satisfied by starting from an initialassumption of αk = 1 and scaling it back towards 0 by a constant multiplier, untilthe sufficient decrease condition of Equation 5.27 is satisfied. After meeting theseconditions, Equation 5.24 is used to produce the next solution estimate. The itera-tions are terminated and xk accepted as the solution once the norm of the gradient‖∇W(xk)‖ is sufficiently close to zero.5.3 Elastic Strain Energy DerivativesThe strain energy density gradient and Hessian are best computed with help of thevec operator [37] which converts a matrix to a vector by stacking the columns ofthe matrix.vec(A) = vec( a11 a12a21 a22 ) =( a11a21a12a22), (5.30)and the Kronecker product (A⊗B).The gradient of the stretch energy density is derived in Appendix AdΨdvec(x)=dΨdvec(F)dvec(F)dvec(x), (5.31)40with vectorized first Piola-Kirchoff stressdΨdvec(F)= vec(FC−1(µ12(C−C−1)+λ ln(J)I))T , (5.32)and the vectorized derivative of the deformation gradient (recall that R is from theQR decomposition of the material position derivative DX)dvec(F)dvec(x)=(R−T(−1 1 0−1 0 1)⊗ I3) . (5.33)The Hessian requires use of the symmetrizing tensor Sym which maps a vector-ized matrix onto its symmetric part. That is,Sym vec(A) = vec(A+AT2). (5.34)The Hessian also has a derivation in Appendix A, provided here for reference.d2Ψdvec(x)T dvec(x)=dvec(F)dvec(x)T d2Ψdvec(F)T dvec(F)dvec(F)dvec(x), (5.35)with first Piola-Kirchoff stress derivatived2Ψdvec(F)T dvec(F)=µ2((I−C−2)⊗ I3)+µ(I⊗F)Sym((C−2⊗C−1)+(C−1⊗C−2))(I⊗FT )+λ ln(J)(C−1⊗ I3)−2λ ln(J)(I⊗F)Sym(C−1⊗C−1)(I⊗FT )+λ (I⊗F)vec(C−1)vec(C−1)T (I⊗FT ). (5.36)The bending energy gradient is defined in terms of the edges ei, normals n j, anddihedral angle θ of the bending element described in Section 5.1.1 and illustratedin Figure 5.2. A full derivation is provided in Appendix A with the final result41copied here.dΨbenddvec(x)=dΨbenddθdθdvec(e)dvec(e)dvec(x), (5.37)with energy functional derivativedΨbenddθ= D(2tan(θ2)−2tan(θ¯2))sec2(θ2), (5.38)and dihedral angle derivativedθdvec(e)=dθde0Tdθde1Tdθde2TT=‖e1‖‖n1‖(eˆ0 · eˆ1)nˆ1 +‖e2‖‖n2‖(eˆ0 · eˆ2)nˆ2− ‖e0‖‖n1‖ nˆ1− ‖e0‖‖n2‖ nˆ2T, (5.39)and edge vector derivativedvec(e)dvec(x)=((−1 1 0 0−1 0 1 0−1 0 0 1)⊗ I3). (5.40)The bending Hessian is rather large and unwieldy. At a high level it isd2Ψbenddvec(x)T dvec(x)=dvec(e)dvec(x)T d2Ψbenddvec(e)T dvec(e)dvec(e)dvec(x)(5.41)d2Ψbenddvec(e)T dvec(e)=d2Ψbenddθ 2dθdvec(e)T dθdvec(e)+dΨbenddθd2θdvec(e)T dvec(e), (5.42)42withd2Ψbenddθ 2=D(2− cos(θ)− sin(θ) tan(θ¯))sec4(θ2), (5.43)andd2θdvec(e)T dvec(e)=d2θdeT0 de0d2θdeT0 de1d2θdeT0 de2d2θdeT0 de1T d2θdeT1 de10d2θdeT0 de2T0 d2θdeT2 de2 , (5.44)where the formulae for d2θdeTi de jare presented in Appendix A.5.4 IndefinitenessThe k-th iteration search direction pk defined by∇2W(xk)pk =−∇W(xk) (5.45)may not be a descent direction since the Hessian matrix∇2W(xk) may be indefinite.In fact, the Hessian is only positive definite in some vicinity of a minimum since thetotal potential energy function W is generally not convex. This can be illustratedwith a simple 1D example.Consider a segment of a 1D elastic strand embedded in a 2D space. Fix oneendpoint at the origin with the other having position x ∈ R2. Assume the materialspace length of the segment to be 1 and the current length of the segment is ‖x‖.Using a constitutive model quadratic in log-strain, the total elastic strain energy(illustrated in Figure 5.5) is given byW =12log2 (‖x‖) . (5.46)Taking derivatives, the elastic force isdWdx=log(‖x‖)‖x‖2 xT , (5.47)43−2−1012−2−10120123456Plot of W = 0.5 log(|x|)2Figure 5.5: Plot of Equation 5.46 as a function of endpoint position xwith Hessian matrixd2Wdx2=log(‖x‖)‖x‖2 I +1−2log(‖x‖2)‖x‖4 xxT . (5.48)Let xˆ = x‖x‖ and nˆ be a unit vector perpendicular to xˆ. We can substitute I =nˆnˆT + xˆxˆT which leads tod2Wdx2=log(‖x‖)‖x‖2 nˆnˆT +1− log(‖x‖)‖x‖2 xˆxˆT . (5.49)Thus we have constructed a spectral decomposition of the Hessian which has neg-ative eigenvalues under two conditions. The lesser condition corresponds to thesecond term being negative when the current length ‖x‖ exceeds e. The othercondition is more troubling: under any compression the first term will become44(a) Using the indefinite Hessian matrix (b) Using the nearest positive semi-definitematrixFigure 5.6: The osculating paraboloid of Equation 5.46 about the point x =(0.11,−0.68). In (a) the shape is a hyperbolic paraboloid, indicating anindefinite Hessian. In (b) the modified Hessian is an elliptic paraboloidbecause it has no negative eigenvalues.negative. In this configuration, this corresponds to a rotation of the segment aboutthe fixed endpoint. Clearly the energy is defined in a manner that is independent ofrotation, as is required of any constitutive model, but the quadratic approximationthat is solved at each Newton step is unable to preserve this property. The quadraticapproximation at an iterate is illustrated in Figure 5.6.The previous example is illustrative, but not necessarily a strong motivationsince applying Newton’s method will converge to the same solution as the projectedversion. The saddle point of the true Hessian corresponds to the minimum normsolution to the projected version. By adding another strand segment (with materialspace length 1 as well), connected to the free vertex x and a new fixed vertex at(√22 ,−√22 ) we construct a buckling scenario. The free vertex connecting the twostrands must buckle out-of-plane in order to return to the reference lengths. Theenergy of this system is plotted in Figure 5.7. Without modification of the Hessian,Newton’s method may converge to the indicated saddle point, or one of the trueminima to each side depending on the initial guess.Nocedal and Wright [41] describe a variety of positive definite modifications tothe Hessian. By considering the spectral decomposition of the Hessian matrix, they45Figure 5.7: Plot of Equation 5.46 with two strand segments, forced to buckle.The red ‘x’-mark indicates a saddle point, confirming the non-convexityof the system.identify the negative eigenvalues and corresponding eigenvectors as the problem-atic directions. One strategy is to take the absolute value of these eigenvalues, orto clamp them to some small positive value. This technique is rarely used directly,since the spectral decomposition of a large matrix is prohibitively expensive. Inthe case of a finite element formulation, the Hessian is a sum of “smaller”, elementHessian matrices Ki such that∇2W(xk) =∑iKi. (5.50)These element Hessians are smaller in the sense that they are low rank, with veryfew non-zero elements due to the compact support of the element shape functions.46Positive semi-definiteness requiresyT(∇2W(xk))y = yT(∑iKi)y =∑iyT Kiy≥ 0. (5.51)Teran et al. [51] enforce positive semi-definiteness on the element Hessians(i.e., Ki) directly, by clamping their eigenvalues to be non-negative. McAdamset al. [39] enforce positive semi-definiteness on the matrix-vector product in thesame manner. Since the element Hessians are of low rank compared to the fullHessian, this technique is significantly more efficient. In general, this modificationwill reduce the expected convergence to merely super-linear instead of quadratic.McAdams et al. [39] find that this method is convergent and in practice effectivelyretains the quadratic convergence of the true Newton method, even for large non-smooth deformations.Unfortunately this technique is ineffective when applied to separately to themembrane, bending, and volumetric elements. By projecting each of these elementHessians onto the nearest positive semi-definite matrix before summing, the result-ing matrix is quite different from summing the elements first and then projecting.Even very near a critical point (i.e., gradient norm approaching zero) the individ-ual element Hessians are not positive definite, only their sums are. Summing themembrane and bending element Hessians before projecting is difficult because theelements are defined on different components of the simulation’s triangle mesh. Inpractice projecting the individual element Hessians appears to be effective, onlywhen far from a solution we expect the unprojected system Hessian to becomepositive definite when near the solution. The technique used here is to begin thesimulation by projecting the individual element Hessians, disabling the projectionas the gradient norm becomes small. If the unprojected system Hessian gener-ates a non-descent direction (indicated by a non-negative directional derivative),we re-enable the projection until the gradient norm decreases further. Eventually,quadratic convergence is expected.Curiously, partitioning the membrane Hessian into three pieces and addingthem to the attached bending element before projection did not improve the be-havior near the solution.47Chapter 6ResultsAs discussed in the introduction and Chapter 2, the shape of buckling skin is dis-tinctive, with a hierarchical collection of wavelengths and an overall shape that hasbroad, rounded peaks and sharp valleys. To illustrate the results of our skin model,we have implemented the quasistatic simulation from Chapter 5 in both 2D and 3D,for a patch of skin undergoing uniaxial compression. We show that our results pro-duce the expect wrinkling profile, while alternative models do not generally meetour wrinkle shape requirements.6.1 2D resultsWe construct a two dimensional patch of skin using our three layer model, withprescribed displacements for the bottom, left, and right boundary vertices. The in-terior and top boundary vertices (i.e., the SC interior) remained unconstrained. Thepatch is slowly compressed by displacing all vertices (included the fixed boundary)towards the center of the patch, over 40 equal increments. After each incrementaldisplacement the quasistatic simulation is run until equilibrium of internal elasticforces is found. Since the boundary vertices with prescribed physical space po-sitions are changing over time, this is sufficient for our quasistatic simulation toproduce an animation of the patch compression, similar to pinching the skin of theforearm with the fingers.Using this setup, the following experiments investigate the impact of several48Parameter ValueMaterial Length 2.94 cmPrestrain 1.02×SC Young’s Modulus (Ysc) 100000 PaSC Poisson’s ratio (νsc) 0.495SC number of elements 250SC thickness 50 µmDermis Young’s Modulus (Yd) 100 PaDermis Poisson’s ratio (νd) 0.42Dermis number elements (thick) 3Dermis thickness 1 mmHypodermis Young’s Modulus (Yh) 10 PaHypodermis Poisson’s ratio (νh) 0.42Hypodermis number elements (thick) 3Hypodermis thickness 1 mmTable 6.1: Base parameters used in 2D experiments.design choices of our skin model. The baseline parameters are listed in Table 6.1.The relative stiffnesses are from Cerda and Mahadevan [8], the thicknesses aremodified from Flynn and McCormack [14], and the SC is nearly incompressible asin Geerligs et al. [19], Hendriks et al. [25]. Other parameters are chosen for effectwhen not available in the literature.6.1.1 Experiment 1Experiment 1 compares our large compression, non-linear constitutive model to acorotational linear model. We expect the linear model to allow elements of thedermis and hypodermis to invert (i.e., to take on negative volume) due to the finitestrain response to strong compression. The simulation results using our nonlinearconstitutive model are shown in Figure 6.1, while the corotational linear constitu-tive model is shown in Figure 6.2. Smooth peaks and sharp valleys are present inboth models, but the nonlinear strain response in our constitutive model appearsto be crucial for capturing multiple buckling wavelengths. Additionally, the tensilestrain in the dermis adjacent to wrinkle valleys is interesting to observe. Figure 6.2bshows that under further compression the linear model does not prevent inversion;49(a) Skin patch compressed to 86% its original size.(b) Skin patch compressed to 80.5% its original size.Figure 6.1: A 3cm long, 2D cross-section of our 3-layered skin model underuniform compression.50Parameter ValueSubstrate Young’s Modulus (Yd) 100 PaSubstrate Poisson’s ratio (νd) 0.42Substrate number elements (thick) 6Substrate thickness 2 mmTable 6.2: Parameters specific to Experiment 3.elements in the vicinity of the valleys have negative area.6.1.2 Experiment 2Experiment 2, similar to Experiment 1, compares the 2tan(θ2)SC bending modelto one which is linear in dihedral angle. Refer to Section 5.1.1 for a comparisonof these bending energy density functions. Figure 6.3 confirms that the bendingenergy model has a fairly subtle impact on the simulation result. We observe thatthe large wrinkle furrows (i.e., where the dermis buckles relative to the hypodermis)are much sharper using the linear bending model. In Figure 6.1 the same furrowsare slightly smoothed out.6.1.3 Experiment 3This experiment compares our three layered skin model, to a two layered modelconsisting of the stratum corneum and a substrate modeling an averaged livingepidermis, dermis and hypodermis. Essentially the dermis of our skin model isextended to replace the hypodermis. In Figure 6.4 we observe a single bucklingwavelength, even under large compression. This is consistent with Cerda and Ma-hadevan [8] which predicts a fixed wavelength proportional to the relative stiffnessof the SC and substrate. The parameters of the substrate used in this experimentare listed in Table 6.2, while the SC parameters are the same as the previous exper-iments.51(a) Skin patch compressed to 86% its original size.(b) Skin patch compressed to 80.5% its original size.Figure 6.2: A 3cm long, cross-section of our 3-layered skin model using acorotational linear constitutive mode for the SC membrane, dermis andhypodermis.52(a) Skin patch compressed to 86% its original size.(b) Skin patch compressed to 80.5% its original size.Figure 6.3: A 3cm long, cross-section of our 3-layered skin model using alinear bending mode for the SC.53(a) Skin patch compressed to 86% its original size.(b) Skin patch compressed to 80.5% its original size.Figure 6.4: A 3cm long, cross-section of a 2-layered skin model of the SCand underlying substrate.54Parameter ValueLiving Epidermis Young’s Modulus (Ye) 10 PaLiving Epidermis Poisson’s ratio (νe) 0.42Living Epidermis number elements (thick) 2Living Epidermis thickness 0.4mmDermis Young’s Modulus (Yd) 100 PaDermis Poisson’s ratio (νd) 0.42Dermis number elements (thick) 8Dermis thickness 1.6 mmTable 6.3: Parameters specific to Experiment 46.1.4 Experiment 4In Experiment 4, we consider an alternative three layered model proposed byLe´veˆque and Audoly [32], consisting of the SC, living epidermis and deep dermis.As explained previously, the living epidermis is comparatively less stiff than the SCabove and the deep dermis below it; this contrasts with our model where stiffnessdecreases for each successive layer. Since this intermediate layer is less stiff thanthe one below it, the theory of Cerda and Mahadevan [8] would not suggest theoccurrence of any buckling. Kuwazuru et al. [29] include a similar configuration intheir analysis, but find that it produces buckling wavelengths that are basically thesame as that of the SC buckling relative to its foundation. This suggests that thismodel will not be appreciably different from the two layer model from Experiment3. Figure 6.5 supports this prediction, since only a single buckling wavelength ispresent under various compressions. The parameters for the living epidermis anddeep dermis are listed in Table 6.3, while the SC parameters are the same as theprevious experiments.6.1.5 Experiment 5Wrinkling occurs when the bending effects of the SC compete with the stretch andshear effects of the dermis and hypodermis underneath [12]. The SC penalizesshort wavelength buckling while the stretching of the foundation layers penalizeslong wavelengths. Experiment 5 investigates these predictions by comparing our55(a) Skin patch compressed to 86% its original size.(b) Skin patch compressed to 80.5% its original size.Figure 6.5: A 3cm long, cross-section of a 3-layered skin model of the SC,living epidermis and dermis.56simulations while varying the thickness of the SC layer. Our results in Figure 6.6show an intuitive relationship between the specified thickness and the resultingsize and smoothness of the wrinkles, for a fixed compression. We observe pleasingresults with thickness between 50µm and 25µm. As the thickness approaches20µm [14] the buckling wavelength becomes difficult to resolve.6.2 3D ResultsAnalogous to the previous section, we construct a three dimensional patch of skinusing our three layer model, with the same boundary conditions and incrementallyapplied compressive strain. As in two dimensions, we observe multiple, hierarchi-cal wavelengths in Figure 6.7a and the sharp valleys and rounder peaks that wedesire. In Figure 6.7b we randomize material properties of the SC to break thesymmetry inherit in uniaxial compression to achieve wrinkles which vary throughthe width of the patch.When utilizing the fully 3D model, we find that the parameters suggested byFlynn and McCormack [14] produce excellent results unlike our 2D model whichneeded modification. The parameters used when producing our 3D results arelisted in Table 6.4. For the most part, these are consistent with Flynn and McCor-mack [14] with relative stiffness of layers from Cerda and Mahadevan [8].57(a) SC thickness 100 µm(b) SC thickness 50 µmFigure 6.6: A 3cm long, cross-section of our 3-layered skin model with var-ied SC thickness, compressed to 86% its original size.58(c) SC thickness 37.5 µm(d) SC thickness 33.3 µmFigure 6.6: A 3cm long, cross-section of our 3-layered skin model with var-ied SC thickness, compressed to 86% its original size.59(e) SC thickness 25 µmFigure 6.6: A 3cm long, cross-section of our 3-layered skin model with var-ied SC thickness, compressed to 86% its original size.Parameter ValueMaterial Size 3cm × 2cmSC Young’s Modulus (Ysc) 1 MPaSC Poisson’s ratio (νsc) 0.495SC number of elements 258 × 172SC thickness 22.5 µmDermis Young’s Modulus (Yd) 1 kPaDermis Poisson’s ratio (νd) 0.4Dermis min. elements (thick) 3Dermis thickness 1.2 mmHypodermis Young’s Modulus (Yh) 13 kPaHypodermis Poisson’s ratio (νh) 0.4Hypodermis min. elements (thick) 3Hypodermis thickness 1.5 mmTable 6.4: Parameters used in our 3D results60(a) A 3cm by 2cm patch of skin under large uniaxial compression. Multiple bucklingwavelengths are present.(b) A 3cm by 1cm patch of skin under small uniaxial compression. By randomizing mate-rial properties, the symmetry of the buckling is broken.61Chapter 7Discussion7.1 DiscussionAt the heart of this thesis, we pose the question of which properties of skin arenecessary and sufficient for producing plausible wrinkles when a simulated patchof skin has compressive strain applied. As mentioned in Chapter 1, skin is a non-homogeneous, anisotropic, non-linear, viscoelastic, multi-layered material. All ofthese properties are present in many state-of-the-art models used in the biome-chanics and dermatology literature. On the other hand, few of these properties arepresent in the models typically found in computer graphics. Unsurprisingly, thishas left artists using techniques that are often intended for cloth simulation, pro-ducing wrinkles which are unsatisfying or otherwise nonphysical; they generallyresort to manual sculpting methods for producing pleasing wrinkles. Our goal isto develop an efficient model for wrinkling skin that is physically plausible whileeliminating phenomena that are not necessary.As is evident when pinching the skin of the forearm, human skin buckles withat least two distinct wavelengths. The profile of these buckles is also decidedly notin the shape of sine waves, but has rounded peaks and sharp valleys in general. Ourmodel presented in Chapter 4 captures this phenomenon (see Figure 6.7a), whilesignificantly simplifying the constitutive models and properties compared to thatof Flynn and McCormack [14] and Magnenat-Thalmann et al. [36] which inspiredit.62Flynn [13] and Limbert [35] describe several constitutive models used in biome-chanics and cosmetics. They emphasize the role of the collagen fibers in the der-mis, which are known to be the dominant structural component of skin under ten-sion. The distribution and orientation of these fibers are a significant source ofanisotropy as well. Our skin model has elided these properties in exchange for amuch simpler formulation and we have found that the overall behavior under com-pression is still sufficiently captured. This is likely because the collagen fibers donot contribute much to the elastic behavior of skin when compressed, since they re-main crimped, unextended, and randomly oriented [13]. The collagen and groundsubstance of the dermis provide the elastic and volume preserving response whichare critical during compression [10].The viscoelastic effects that are present in the skin are ignored by our choiceof a quasistatic formulation. The explicit time-dependence is removed, making theconcept of rate dependence superfluous. Our results show that a multi-wavelengthbuckling pattern of the correct overall shape does not require this phenomenon tobe modeled explicitly. On the other hand there are artifacts associated to the qua-sistatic simulation, rapid configuration changes (i.e., snap through), that might beimpacted by a strain rate dependence. Testing the impact of viscosity is some-what beyond the scope of our work, so no definitive conclusions should be drawn.Nonetheless, it seems likely that plausible wrinkles can be generated without ratedependent effects.Real human skin is composed of several layers with distinct mechanical prop-erties. The SC is extremely thin and stiff, while the living epidermis is much lessso. The dermis overall continues the trend of decreasing stiffness through the skinthickness though the papillary dermis is somewhat less stiff than the reticular der-mis owing to it being less dense in elastic fibers. The hypodermis is the least stiffof all the layers. The work of Cerda and Mahadevan [8] supports a choice of atleast three layers in order to capture dual wavelength wrinkling patterns, thoughthe literature has been conflicted over the choices for these layers. Section 6.1.3confirms this prediction, where a simulated two layer model produces only sinu-soidal buckling with a fixed wavelength.Our three layer model, based closely on that of Flynn and McCormack [14], hassuccessfully reproduced the desired wrinkle profile by modeling the SC, dermis and63hypodermis. Another three layered configuration is considered in Section 6.1.4,based on the work of Le´veˆque and Audoly [32] and Magnenat-Thalmann et al. [36],but we find that it is not sufficient for producing multiple buckling wavelengths.This is not surprising when viewed through the establish buckling theory whichprovides a fixed wavelength from the ratio of stiffness of a thin sheet atop a lessstiff substrate. With a low stiffness layer between the SC and dermis we wouldnot expect any meaningful buckling wavelength to emerge, beyond that of the SCbuckling relative to the layers below it. In the work of Magnenat-Thalmann et al.,they require a wavy interface between the living epidermis and dermis, which wehave not included. An extension of this idea to three dimensions is unclear.It is also illuminating to compare our results to that of Li and Kry [34]. In theirwork, they use two thin shells representing the epidermis and dermis connectedto an underlying body. The two shell layers are connected to themselves and thebody by specialized constraints that are designed to allow the outer shell to bucklerelative to the inner shell with the wavelengths predicted by Cerda and Mahade-van and skin properties derived from the literature. Their constraints eliminate thevolumetric behavior that we have explicitly modeled. While their published resultsshow multiple buckling wavelengths, the overall shape is sinusoidal so we con-clude that the volumetric properties of the dermis and hypodermis are necessaryfor simulating wrinkles with the correct shape.The novel constitutive model proposed in Chapter 4 is designed to react stronglyto large compression, while also being computationally efficient. Most notably itignores the complex treatment of collagen fiber distributions that are often presentin modern skin models. As previously discussed, this simplification appears to bereasonable since it does not affect the overall buckling profile. Our results sup-port the necessity of a strongly non-linear model however; in Section 6.1.1 a linearstrain constitutive model is shown to be insufficient for maintaining multiple buck-ling wavelengths. The linear model is also shown to have a threshold of compres-sion, after which it fails by allowing some elements to take on negative volume. Asimilar finding for the bending model is also shown in Section Limitations and Future WorkIt has been previously mentioned that the literature is uncertain on the values forthe stiffness and other properties of the various skin layers. Practically every paperhas a different idea for the relative stiffness of the SC to the dermis, etc. No doubtthis stems from the lack of consensus on which layers are used in skin models, aswell as the variation of skin’s properties by individual, location on the body, rel-ative humidity, and other factors. Additionally, the tools used to measure skin invivo are likely unable to measure the properties of a single layer without confound-ing the properties of the others as well. With this in mind, we should note that wedid not intend to compare our results to measurements in the literature or validatethem by experiment. Our work can certainly be improved by using realistic mea-surements for skin thickness and the material properties of our chosen layers. Asimple experiment, not unlike Flynn and McCormack [14], would also be effectivefor testing our skin model. With data captured from such an experiment we couldfurther validate our model, fit parameters, or justify more complicated constitutivemodels.The SC is generally modeled as a nearly incompressible material, with Poissonratio of 0.495 [e.g., 19, 25]. Our variational approach to incompressibility requiresa large stiffness with respect to the measure of volume change (i.e., that associatedto the Lame´ coefficient λ ). As the name suggests, this leads to a numerically stiffformulation indicated by a Hessian with large condition number (in the hundredsof thousands or higher) which affects the efficiency of our Newton based solver.In terms of optimization techniques, this is equivalent to the penalty method forenforcing constraints, which is known to be generally inefficient. Our methodcould be significantly improved by handling this constraint explicitly, perhaps viaLagrange multipliers. Bonet and Wood [6] suggest that an efficient technique fornearly incompressible elastic materials is to use a mixed finite element formulationdue to Simo et al. [49].When we purposefully chose to discard the effects of inertia and rate depen-dence, we enabled the use of a quasistatic formulation where at each time instantan equilibrium configuration was sought instead of integrating the system in timefrom a previous state. This allowed us the elegant option of interpreting our skin65motion as a minimization problem, and applying the machinery of optimization. Abenefit of this approach is that we no longer have any stability limitation on timestep (since time is not explicitly tracked). An unpleasant side effect is that qua-sistatic solutions adjacent in time (assuming time-varying boundary conditions andmaterial properties) may have significantly different configurations, known as pop-ping artifacts in animation. This is especially prevalent as the skin instantaneouslybuckles into a new shape after a critical compression is surpassed. Thomaszewskiet al. [52] describe a strain limiting framework, which could be applied to enforce amaximum strain change relative to the previous solution. Alternatively, Bonet andWood [6] discuss arc length methods that constrain the iterative solution to followa certain route toward the equilibrium path.The excessive number of tetrahedral elements required in the dermis and hy-podermis are a limiting factor for the application of our skin model. Since theSC is involved in high frequency wrinkling, it needs to be finely tessellated. Theunderlying dermis and hypodermis, on the other hand, are generally deformed ina much smoother manner except for a small region directly adjacent to the SC.This suggests an adaptive scheme for filling the volume may be effective. On theother hand, there are well defined layers (i.e., the dermis and hypodermis) wherethe material properties change discontinuously so there are restrictions on the dis-cretization. Namely, the interface between the dermis and hypodermis must bespecified. Constructing a scheme for adaptively discretizing the volume would bea useful addition to this model. Alternatively, it may make sense to consider socalled ”long elements” from Balaniuk and Salisbury [3]. In this finite element for-mulation, the number of elements is proportional to the square of the length of aside rather than its cube as in a standard meshing using cubes or tetrahedra. Whilethe simple one point integration scheme would be lost, the savings in total elementcount could be substantial.The indefiniteness of the Hessian matrix for our problem is discussed in Sec-tion 5.4 where we determine that the standard technique used in computer graphicsdoes not apply to our problem. Specifically, the partially overlapping nature ofthe bending and membrane elements does not admit an indefiniteness fixing mod-ification on a per-element basis. Our solution applies the standard, non-optimalfix for a few initial iterations then proceeds to use the unmodified Hessian when-66ever it does not generate a descent direction. While this technically preserves thequadratic convergence behavior of Newton’s method when sufficiently close to asolution, we are unaware of an efficient method for determining when we are in factsufficiently close. The weakness of this approach is that when using the unmod-ified, indefinite Hessian we may generate a direction that is technically a descentdirection, but unnecessarily increases energy in some elements where the negativeeigenvalues are present. A future work might investigate efficient mechanisms forassuring positive definiteness of the Hessian, while also minimizing the differencewith the true Hessian so that quadratic convergence is retained.7.3 ConclusionOur results clearly produce the overall buckling shape illustrated in Figure 2.1. Fur-thermore, we observe multiple wrinkling wavelengths as evidenced by Figure 6.1and Figure 6.7a. We conclude that modeling the SC, dermis, and hypodermis withour log-strain based constitutive model is sufficient for achieving a physically plau-sible wrinkling model. Our experiments with fewer layers (i.e., Section 6.1.3) andalternative layer choices (i.e., Section 6.1.4) confirm that three layers and our rel-ative organization of those layers is necessary for correct wrinkling behavior. Ad-ditionally a nonlinear constitutive model that responds strongly to compression isalso necessary. We believe that we have succeeded in producing a minimal, plau-sible model of wrinkling skin that is motivated by the actual structure and physicsof skin.67Bibliography[1] U. M. Ascher and E. Boxerman. On the modified conjugate gradient methodin cloth simulation. The Visual Computer, 19(7-8):526–531, 2003. → pages12[2] B. Audoly and Y. Pomeau. Elasticity and geometry: from hair curls to thenon-linear response of shells. Oxford University Press, 2010. → pages 25,28, 29[3] R. Balaniuk and K. Salisbury. Dynamic simulation of deformable objectsusing the long elements method. In Haptic Interfaces for VirtualEnvironment and Teleoperator Systems, 2002. HAPTICS 2002. Proceedings.10th Symposium on, pages 58–65. IEEE, 2002. → pages 66[4] D. Baraff and A. Witkin. Large steps in cloth simulation. In Proceedings ofthe 25th Annual Conference on Computer Graphics and InteractiveTechniques, SIGGRAPH ’98, pages 43–54, New York, NY, USA, 1998.ACM. ISBN 0-89791-999-8. doi:10.1145/280814.280821. URLhttp://doi.acm.org/10.1145/280814.280821. → pages 12[5] Z. P. Bazˇant. Easy-to-compute tensors with symmetric inverseapproximating hencky finite strain and its rate. Journal of engineeringmaterials and technology, 120(2):131–136, 1998. → pages viii, 16, 27, 28[6] J. Bonet and R. D. Wood. Nonlinear continuum mechanics for finite elementanalysis. Cambridge University Press, second edition, 2008. → pages 13,65, 66[7] R. Bridson, S. Marino, and R. Fedkiw. Simulation of clothing with folds andwrinkles. In Proceedings of the 2003 ACM SIGGRAPH/EurographicsSymposium on Computer Animation, SCA ’03, pages 28–36, Aire-la-Ville,Switzerland, Switzerland, 2003. Eurographics Association. ISBN1-58113-659-5. URL http://dl.acm.org/citation.cfm?id=846276.846281. →pages 12, 25, 3168[8] E. Cerda and L. Mahadevan. Geometry and physics of wrinkling. Physicalreview letters, 90(7):074302, 2003. → pages 8, 9, 11, 24, 25, 26, 49, 51, 55,57, 63, 64[9] L. D. Cutler, R. Gershbein, X. C. Wang, C. Curtis, E. Maigret, L. Prasso, andP. Farson. An art-directed wrinkle system for cg character clothing. InProceedings of the 2005 ACM SIGGRAPH/Eurographics Symposium onComputer Animation, SCA ’05, pages 117–125, New York, NY, USA, 2005.ACM. ISBN 1-59593-198-8. doi:10.1145/1073368.1073384. URLhttp://doi.acm.org/10.1145/1073368.1073384. → pages 12[10] C. H. Daly. Biomechanical properties of dermis. Journal of InvestigativeDermatology, 79:17–20, 1982. → pages 4, 63[11] D. C. Drucker. A definition of stable inelastic material. Technical report,DTIC Document, 1957. → pages 21[12] K. Efimenko, M. Rackaitis, E. Manias, A. Vaziri, L. Mahadevan, andJ. Genzer. Nested self-similar wrinkling patterns in skins. Nature materials,4(4):293–297, 2005. → pages 8, 9, 55[13] C. Flynn. Fiber matrix models of the dermis. In Computational Biophysicsof the Skin, pages 133–160. CRC Press, 2014. → pages 3, 10, 24, 63[14] C. Flynn and B. A. McCormack. Finite element modelling of forearm skinwrinkling. Skin research and technology, 14(3):261–269, 2008. → pages 3,4, 10, 24, 25, 26, 30, 49, 57, 62, 63, 65[15] C. Flynn and B. A. McCormack. A three-layer model of skin and itsapplication in simulating wrinkling. Computer methods in biomechanics andbiomedical engineering, 12(2):125–134, 2009. → pages[16] C. Flynn and B. A. McCormack. Simulating the wrinkling and aging of skinwith a multi-layer finite element model. Journal of biomechanics, 43(3):442–448, 2010. → pages 10, 24[17] J. Fro¨hlich, S. Huclova, C. Beyer, and D. Ernib. Accurate multiscale skinmodel suitable for determining the sensitivity and specificity of changes ofskin components. pages 353–393. CRC Press, 2014. → pages viii, 3[18] A. Garg, E. Grinspun, M. Wardetzky, and D. Zorin. Cubic shells. InProceedings of the 2007 ACM SIGGRAPH/Eurographics Symposium onComputer Animation, SCA ’07, pages 91–98, Aire-la-Ville, Switzerland,69Switzerland, 2007. Eurographics Association. ISBN 978-1-59593-624-0.URL http://dl.acm.org/citation.cfm?id=1272690.1272703. → pages viii, 34,35[19] M. Geerligs, L. Van Breemen, G. Peters, P. Ackermans, F. Baaijens, andC. Oomens. In vitro indentation to determine the mechanical properties ofepidermis. Journal of biomechanics, 44(6):1176–1181, 2011. → pages 26,49, 65[20] R. Goldenthal, D. Harmon, R. Fattal, M. Bercovier, and E. Grinspun.Efficient simulation of inextensible cloth. In ACM SIGGRAPH 2007 Papers,SIGGRAPH ’07, New York, NY, USA, 2007. ACM.doi:10.1145/1275808.1276438. URLhttp://doi.acm.org/10.1145/1275808.1276438. → pages 12[21] E. Grinspun, A. N. Hirani, M. Desbrun, and P. Schro¨der. Discrete shells. InProceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium onComputer Animation, SCA ’03, pages 62–67, Aire-la-Ville, Switzerland,Switzerland, 2003. Eurographics Association. ISBN 1-58113-659-5. URLhttp://dl.acm.org/citation.cfm?id=846276.846284. → pages 25, 31, 34, 35[22] E. Grinspun, M. Desbrun, K. Polthier, P. Schro¨der, and A. Stern. Discretedifferential geometry: an applied introduction. ACM SIGGRAPH Course, 7,2006. → pages 34[23] S. Hadap, E. Bangerter, P. Volino, and N. Magnenat-Thalmann. Animatingwrinkles on clothes. In Proceedings of the conference on Visualization’99:celebrating ten years, pages 175–182. IEEE Computer Society Press, 1999.→ pages 11[24] J. Hatzis. The wrinkle and its measurement: A skin surface profilometricmethod. Micron, 35(3):201–219, 2004. → pages viii, 8[25] F. Hendriks, D. Brokken, C. Oomens, D. Bader, and F. Baaijens. Therelative contributions of different skin layers to the mechanical behavior ofhuman skin in vivo using suction experiments. Medical engineering &physics, 28(3):259–266, 2006. → pages 26, 49, 65[26] J. Jimenez, J. I. Echevarria, C. Oat, and D. Gutierrez. GPU Pro 2, chapterPractical and Realistic Facial Wrinkles Animation. AK Peters Ltd., 2011. →pages 1170[27] L. Kavan, D. Gerszewski, A. W. Bargteil, and P.-P. Sloan. Physics-inspiredupsampling for cloth simulation in games. In ACM SIGGRAPH 2011Papers, SIGGRAPH ’11, pages 93:1–93:10, New York, NY, USA, 2011.ACM. ISBN 978-1-4503-0943-1. doi:10.1145/1964921.1964988. URLhttp://doi.acm.org/10.1145/1964921.1964988. → pages 12[28] S. Kimmerle, M. Wacker, and C. Holzer. Multilayered wrinkle textures fromstrain. In VMV, pages 225–232, 2004. → pages 11[29] O. Kuwazuru, J. Saothong, and N. Yoshikawa. Mechanical approach toaging and wrinkling of human facial skin based on the multistage bucklingtheory. Medical engineering & physics, 30(4):516–522, 2008. → pages 4, 9,26, 55[30] K. Langer. On the anatomy and physiology of the skin: I. the cleavability ofthe cutis. British journal of plastic surgery, 31(1):3–8, 1978. → pages 30[31] C. Larboulette and M.-P. Cani. Real-time dynamic wrinkles. In ComputerGraphics International, 2004. Proceedings, pages 522–525. IEEE, 2004. →pages 12[32] J. L. Le´veˆque and B. Audoly. Influence of stratum corneum on the entireskin mechanical properties, as predicted by a computational skin model. SkinResearch and Technology, 19(1):42–46, 2013. → pages 8, 10, 24, 25, 55, 64[33] D. Li, S. Sueda, D. R. Neog, and D. K. Pai. Thin skin elastodynamics. ACMTrans. Graph., 32(4):49:1–49:10, July 2013. ISSN 0730-0301.doi:10.1145/2461912.2462008. URLhttp://doi.acm.org/10.1145/2461912.2462008. → pages 5, 11, 34[34] P. Li and P. G. Kry. Multi-layer skin simulation with adaptive constraints. InProceedings of the Seventh International Conference on Motion in Games,MIG ’14, pages 171–176, New York, NY, USA, 2014. ACM. ISBN978-1-4503-2623-0. doi:10.1145/2668084.2668089. URLhttp://doi.acm.org/10.1145/2668084.2668089. → pages 11, 25, 64[35] G. Limbert. State-of-the-art constitutive models of skin biomechanics. InComputational Biophysics of the Skin, pages 95–131. CRC Press, 2014. →pages 3, 9, 10, 24, 63[36] N. Magnenat-Thalmann, P. Kalra, J. L. Le´veˆque, R. Bazin, D. Batisse, andB. Querleux. A computational skin model: fold and wrinkle formation.IEEE Transactions on Information Technology in Biomedicine, 6(4):317–323, 2002. → pages 7, 8, 10, 24, 25, 26, 62, 6471[37] J. R. Magnus and H. Neudecker. Matrix differential calculus withapplications in statistics and econometrics. John Wiley & Sons, thirdedition, 2007. → pages 40, 75, 76, 77[38] W. Maurel, Y. Wu, N. M. Thalmann, and D. Thalmann. Biomechanicalmodels for soft tissue simulation. Springer, 1998. → pages 3, 4, 9, 30[39] A. McAdams, Y. Zhu, A. Selle, M. Empey, R. Tamstorf, J. Teran, andE. Sifakis. Efficient elasticity for character skinning with contact andcollisions. In ACM SIGGRAPH 2011 Papers, SIGGRAPH ’11, pages37:1–37:12, New York, NY, USA, 2011. ACM. ISBN 978-1-4503-0943-1.doi:10.1145/1964921.1964932. URLhttp://doi.acm.org/10.1145/1964921.1964932. → pages 21, 39, 47[40] M. Mu¨ller and N. Chentanez. Wrinkle meshes. In Proceedings of the 2010ACM SIGGRAPH/Eurographics Symposium on Computer Animation, SCA’10, pages 85–92, Aire-la-Ville, Switzerland, Switzerland, 2010.Eurographics Association. URLhttp://dl.acm.org/citation.cfm?id=1921427.1921441. → pages 12[41] J. Nocedal and S. Wright. Numerical optimization. Springer Science &Business Media, 2006. → pages 39, 40, 45[42] R. Ogden. Large deformation isotropic elasticity-on the correlation of theoryand experiment for incompressible rubberlike solids. In Proceedings of theRoyal Society of London A: Mathematical, Physical and EngineeringSciences, volume 326, pages 565–584. The Royal Society, 1972. → pages 23[43] M. Ortiz, R. Radovitzky, and E. Repetto. The computation of theexponential and logarithmic mappings and their first and secondlinearizations. International Journal for Numerical Methods in Engineering,52(12):1431–1441, 2001. → pages 22[44] G. E. Pie´rard, I. Uhoda, and C. Pie´rard-Franchimont. From skin microreliefto wrinkles. an area ripe for investigation. Journal of cosmetic dermatology,2(1):21–28, 2003. → pages 7, 8[45] T. Popa, Q. Zhou, D. Bradley, V. Kraevoy, H. Fu, A. Sheffer, andW. Heidrich. Wrinkling captured garments using space-time data-drivendeformation. In Computer Graphics Forum, volume 28, pages 427–435.Wiley Online Library, 2009. → pages 1272[46] O. Re´millard and P. G. Kry. Embedded thin shells for wrinkle simulation.ACM Trans. Graph., 32(4):50:1–50:8, July 2013. ISSN 0730-0301.doi:10.1145/2461912.2462018. URLhttp://doi.acm.org/10.1145/2461912.2462018. → pages 11, 25[47] R. S. Rivlin and D. Saunders. Large elastic deformations of isotropicmaterials. vii. experiments on the deformation of rubber. PhilosophicalTransactions of the Royal Society of London A: Mathematical, Physical andEngineering Sciences, 243(865):251–288, 1951. → pages 23[48] D. Rohmer, T. Popa, M.-P. Cani, S. Hahmann, and A. Sheffer. Animationwrinkling: Augmenting coarse cloth simulations with realistic-lookingwrinkles. In ACM SIGGRAPH Asia 2010 Papers, SIGGRAPH ASIA ’10,pages 157:1–157:8, New York, NY, USA, 2010. ACM. ISBN978-1-4503-0439-9. doi:10.1145/1866158.1866183. URLhttp://doi.acm.org/10.1145/1866158.1866183. → pages 12[49] J. Simo, R. Taylor, and K. Pister. Variational and projection methods for thevolume constraint in finite deformation elasto-plasticity. Computer Methodsin Applied Mechanics and Engineering, 51(1):177–208, 1985. → pages 65[50] R. Tamstorf and E. Grinspun. Discrete bending forces and their jacobians.Graphical Models, 75(6):362–370, 2013. → pages 31, 36, 39[51] J. Teran, E. Sifakis, G. Irving, and R. Fedkiw. Robust quasistatic finiteelements and flesh simulation. In Proceedings of the 2005 ACMSIGGRAPH/Eurographics Symposium on Computer Animation, SCA ’05,pages 181–190, New York, NY, USA, 2005. ACM. ISBN 1-59593-198-8.doi:10.1145/1073368.1073394. URLhttp://doi.acm.org/10.1145/1073368.1073394. → pages 39, 47[52] B. Thomaszewski, S. Pabst, and W. Strasser. Continuum-based strainlimiting. In Computer Graphics Forum, volume 28, pages 569–576. WileyOnline Library, 2009. → pages 66[53] P. Volino and N. Magnenat-Thalmann. Fast geometrical wrinkles onanimated surfaces. In Seventh International Conference in Central Europeon Computer Graphics and Visualization (Winter School on ComputerGraphics), 1999. → pages 11[54] Y. Zhang and T. Sim. Realistic and efficient wrinkle simulation using ananatomy-based face model with adaptive refinement. In Computer GraphicsInternational 2005, pages 3–10. IEEE, 2005. → pages 1173Appendix ASupporting MaterialsA.1 Membrane & Volume Energy DerivativesThe vector of first derivatives, and matrix of second derivatives of the various elas-tic energies are required to implement the Newton method for minimization. Thederivation of these quantities for the chosen stretching and bending energy densi-ties is described in the following sections.From Section 4.2.1 we have the elastic energy potential density used in the SCmembrane and the dermis/hypodermis volume elementsΨ=µ4tr(C+C−1−2I)+ λ8ln2(detC). (A.1)A.1.1 GradientRecalling the definition of the differentialdΨ=dΨdx: dx, (A.2)74and making use of differential rules from [37]d tr(C) = tr(dC) (A.3)dC = FT dF +dFT F (A.4)dC−1 =−C−1(dC)C−1 (A.5)d ln(detC) = tr(C−1dC), (A.6)we can derive the derivative of the energy density Ψ by observingdW =µ4dtr(C+C−1−2I)+ λ8dln2(detC) (A.7)=µ4tr(dC−C−1(dC)C−1)+ λ4ln(detC) tr(C−1dC)(A.8)=µ4tr(dC−C−2dC)+ λ2ln(J) tr(C−1dC)(A.9)=µ4tr((I−C−2)dC)+ λ2ln(J) tr(C−1dC)(A.10)=µ2tr((I−C−2)FT dF)+λ ln(J) tr(C−1FT dF) (A.11)= tr((µ12(I−C−2)+λ ln(J)C−1)FT dF) (A.12)= tr((µ12(C−C−1)+λ ln(J)I)C−1FT dF) (A.13)=[FC−1(µ12(C−C−1)+λ ln(J)I)] : dF (A.14)=dWdF: dF∴ dWdF= FC−1(µ12(C−C−1)+λ ln(J)I) . (A.15)A.1.2 HessianFor computing the Hessian of this energy function we rely on the second differen-tiald2A = d(dA) = dx :d2Adx2: dx. (A.16)75To alleviate the complexity of derivatives of second order tensors with respect tovectors or other second order tensors, we vectorize our tensor expressions with thevec operator [37] which converts a matrix to a vector by stacking the columns ofthe matrix.vec(A) = vec( a11 a12a21 a22 ) =( a11a21a12a22). (A.17)Making use of various properties[37] of vec and the Kronecker product (A⊗B)vec(AXB) =(BT ⊗A)vec(X), (A.18)tr(AT B)= vec(A)T vec(B), (A.19)Sym vec(A) = vec(A+AT2). (A.20)For example, the differential of the Cauchy stress tensor dC becomesvec(dC) = vec(FT dF +dFT F)(A.21)= 2Sym vec(FT dF)(A.22)= 2Sym(I⊗FT )vec(dF). (A.23)The per-element Hessian matrix is much easier to work with by first vectorizingthe element’s degrees of freedom (i.e., vertex coordinates x) and then computingthe second derivative of energy density with respect to vec(x). By the chain rule:d2Ψdvec(x)T dvec(x)=dvec(F)dvec(x)T d2Wdvec(F)T dvec(F)dvec(F)dvec(x). (A.24)By choosing linear elements, the derivative of the deformation gradient is con-76stant within the element (i.e., dN/dξ = const) giving:dvec(F)dvec(x)=ddvec(x)vec(DxDX†)(A.25)=ddvec(x)vec(xdNdξDX†)(A.26)=ddvec(x)((dNdξDX†)T⊗ I3)vec(x) (A.27)=((dNdξDX†)T⊗ I3), (A.28)where dN/dξ =(−1 −11 00 1)for triangles and dN/dξ =(−1 −1 −11 0 00 1 00 0 1)for tetrahedra.While quite cumbersome, we derive the Hessian matrix by massaging the sec-ond differential of energy density into a canonical form[37]d2Ψ=µ2tr((I−C−2)dFT dF)+µ2tr((C−2(dC)C−1 +C−1(dC)C−2)FT dF)+λ ln(J) tr(C−1dFT dF)−λ ln(J) tr(C−1(dC)C−1FT dF)+λ tr2(C−1FT dF) (A.29)=µ2vec(dF)T ((I−C−2)⊗ I3)vec(dF)+µ vec(dF)T (I⊗F)Sym((C−2⊗C−1)+(C−1⊗C−2))(I⊗FT )vec(dF)+λ ln(J)vec(dF)T (C−1⊗ I3)vec(dF)−2λ ln(J)vec(dF)T (I⊗F)Sym(C−1⊗C−1)(I⊗FT )vec(dF)+λ vec(dF)T (I⊗F)vec(C−1)vec(C−1)T (I⊗FT )vec(dF) (A.30)=vec(dF)Td2Ψdvec(F)T dvec(F)vec(dF),77leading tod2Ψdvec(F)T dvec(F)=µ2((I−C−2)⊗ I3)+µ(I⊗F)Sym((C−2⊗C−1)+(C−1⊗C−2))(I⊗FT )+λ ln(J)(C−1⊗ I3)−2λ ln(J)(I⊗F)Sym(C−1⊗C−1)(I⊗FT )+λ (I⊗F)vec(C−1)vec(C−1)T (I⊗FT ). (A.31)A.2 Bending Energy DerivativesThe bending strain energy density from Section 5.1.1 was chosen asΨb =D2(2tan(θ2)−2tan(θ¯2))2(A.32)A.2.1 GradientGiven the normals n1 = e0× e1 and n2 = e2× e0 of T1 and T2 respectively, thesigned dihedral angle θ is given byθ = atan((n1×n2) · eˆ0n1 ·n2), (A.33)which has derivativedθdei=datan(tanθ)dtanθ(dtanθdn1dn1dei+dtanθdn2dn2dei). (A.34)Equation A.34 shows the intermediate quantities needed to calculate the dihe-dral angle derivative. We first show formulae for the derivatives of the angle withrespect to the two triangle normals. Then, the chain rule is used to extend theseto the three edge vectors of the element, and finally to the vertex positions them-selves. Firstly, the angle is calculated from the inverse tangent function which has78derivativedatan(tanθ)dtanθ=11+ tan2 θ(A.35)= cos2 θ (A.36)=(n1 ·n2‖n1‖‖n2‖)2. (A.37)The derivative of the tangent is computed with respect to the two triangle normals,first n1dtanθdn1=(1n1 ·n2dn1 · (n2× eˆ0)dn1+d(n1 ·n2)−1dn1n1 · (n2× eˆ0))(A.38)=1(n1 ·n2)2((n1 ·n2)(n2× eˆ0)T −n1 · (n2× eˆ0)n2T)(A.39)=n1T(n1 ·n2)2(n2(n2× eˆ0)T − (n2× eˆ0)nT2)(A.40)=n1T(n1 ·n2)2 (E((n2× eˆ0)×n2)) (A.41)=n1T(n1 ·n2)2 (E(n2× (eˆ0×n2))) (A.42)=n1T (Eeˆ0)(n1 ·n2)2 ‖n2‖2. (A.43)Combining Equation A.37 and Equation A.43 via the chain ruledθdn1=n1T‖n1‖2(Eeˆ0) =(nˆ1× eˆ0)T‖n1‖ . (A.44)Similar calculations for n2 showdθdn2=− n2T‖n2‖2(Eeˆ0) =(eˆ0× nˆ2)T‖n2‖ . (A.45)Finally, the derivatives with respect to the edges can be calculated from the normal79derivativesdθde1=‖e0‖‖n1‖ nˆT1 (Eeˆ0)(Eeˆ0) (A.46)=−‖e0‖‖n1‖ (eˆ0× (nˆ1× eˆ0))T (A.47)=−‖e0‖‖n1‖ nˆT1 , (A.48)dθde2=‖e0‖‖n2‖ nˆT2 (Eeˆ0)T (Eeˆ0)T (A.49)=‖e0‖‖n2‖ (eˆ0× (eˆ0× nˆ2))T (A.50)=−‖e0‖‖n2‖ nˆT2 , (A.51)anddθde0=dθdn1dn1de0+dθdn2dn2de0(A.52)=‖e1‖‖n1‖ nˆT1 (Eeˆ0)(Eeˆ1)T +‖e2‖‖n2‖ nˆT2 (Eeˆ0)(Eeˆ2)T (A.53)=‖e1‖‖n1‖ (eˆ1× (nˆ1× eˆ0))T +‖e2‖‖n2‖ (eˆ2× (nˆ2× eˆ0))T (A.54)=‖e1‖‖n1‖(eˆ0 · eˆ1)nˆT1 +‖e2‖‖n2‖(eˆ0 · eˆ2)nˆT2 . (A.55)By combining these derivatives into the 1× 9 row matrix dθdvec(e) , we get anexpression for the derivative in terms of vertex positionsdθdvec(x)=dθdvec(e)dvec(e)dvec(x)(A.56)=(dθde0dθde1dθde2)((−1 1 0 0−1 0 1 0−1 0 0 1)⊗ I3). (A.57)80A.2.2 HessianThe Hessian matrix for the bending energy has the formd2Ψbdvec(x)T dvec(x)=d2Ψbdθ 2dθdvec(x)T dθdvec(x)+dΨbdθd2θdvec(x)T dvec(x)(A.58)The first term (containing d2Ψbdθ 2 ) is an outer product (sincedθdvec(x) is a 1× 12 rowvector) of the gradient discussed previously. The second term containsd2θdvec(x)T dvec(x)=dvec(e)dvec(x)T d2θdvec(e)T dvec(e)dvec(e)dvec(x). (A.59)Since dΨbdθ anddvec(e)dvec(x) are given already, the remaining quantity to derive is the 9×9matrix of second derivatives of θ with respect to the element edge vectors. Thismatrix is derived by observing symmetry and then computing the upper triangle of3×3 blocks d2θdeTi de jd2θdeT1 de1=−‖e0‖ ddn1(n1‖n1‖2)dn1de1(A.60)=−‖e0‖(1‖n1‖2I +n1d‖n1‖−2de1)dn1de1(A.61)=−‖e0‖2‖n1‖2(I−2nˆ1nˆT1)(Eeˆ0) (A.62)=−‖e0‖2‖n1‖2((Eeˆ0) nˆ1nˆT1 + nˆ1nˆT1 (Eeˆ0)T)(A.63)=−‖e0‖2‖n1‖2(mˆ01nˆT1 + nˆ1mˆT01), (A.64)81where the edge normals for e0 in triangles T1 and T2 are respectively mˆ01 = eˆ0× nˆ1and mˆ02 = nˆ2× eˆ0. A similar expression can be found for e2d2θdeT2 de2=−‖e0‖2‖n2‖2(mˆ02nˆT2 + nˆ2mˆT02). (A.65)Since e0 is contained in both triangles, it leads to a more complicated secondderivatived2θdeT0 de0=(eˆ0 · e1) ddn1(n1‖n1‖2)dn1de0+nˆ1‖n1‖dde0(eˆ0 · e1)+(eˆ0 · e2) ddn2(n2‖n2‖2)dn2de0+nˆ2‖n2‖dde0(eˆ0 · e2) (A.66)=‖e1‖2‖n1‖2(eˆ0 · eˆ1)(mˆ1nˆT1 + nˆ1mˆT1)− nˆ1mˆT01‖e0‖2+‖e2‖2‖n2‖2(eˆ0 · eˆ2)(mˆ2nˆT2 + nˆ2mˆT2)− nˆ2mˆT02‖e0‖2 , (A.67)where the remaining edge normals are mˆ1 = nˆ1× eˆ1 and mˆ2 = eˆ2× nˆ2. The re-maining, off-diagonal terms ared2θdeT0 de1= (eˆ0 · e1) ddn1(n1‖n1‖2)dn1de1+nˆ1‖n1‖dde1(eˆ0 · e1) (A.68)=‖e0‖‖e1‖‖n1‖2(eˆ0 · eˆ1)(mˆ01nˆT1 + nˆ1mˆT01)+nˆ1eˆT0‖n1‖ , (A.69)andd2θdeT0 de2= (eˆ0 · e2) ddn2(n2‖n2‖2)dn2de2+nˆ2‖n2‖dde2(eˆ0 · e2) (A.70)=‖e0‖‖e2‖‖n2‖2(eˆ0 · eˆ2)(mˆ02nˆT2 + nˆ2mˆT02)+nˆ2eˆT0‖n2‖ . (A.71)Note that d2θdeT1 de2and its transpose are both zero because of the intermediate normalquantities each only depend on one of e1 or e2.82The Hessian matrix is constructed from these componentsd2θdvec(e)T dvec(e)=d2θdeT0 de0d2θdeT0 de1d2θdeT0 de2d2θdeT0 de1T d2θdeT1 de10d2θdeT0 de2T0 d2θdeT2 de2 . (A.72)83


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items