This tutorial paper provides a step -by-step guide to developing a co mprehensive underst anding of the different forms of the deformation gradient used in Abaqus, and outlines a number of key issues that must be considered when developing an Abaqus user defined material subroutine (UMAT) in which the Cauchy stress is compute d from the deformat ion gradient. Firstly , we examine the “classical” forms of global and local deformation gradients. We then show that Abaqus/ Standard does not use the classical form of the local deformation gradient when continuum elements are used, and w e highlight the im portant implications for UMAT development. We outline the key steps that must be implemented in developing an anisotropic fibre -reinforced hyperelastic UMAT for use with continuum elements and local orientation systems. We also demonstrat e that a classical local deformation gradient is provided by Abaqus/ Standard if structural (shell and membrane) elements are used, and by Abaqus /Explicit for all element types. We emphasise, however, that the majority of biomechanical simulations rely on the use of continuu m elements with a local coordinate system in Abaqus /Standard, and therefore the development of a hyperelastic UMAT requires an in -depth and precise understanding of the form of the non -classical deformation gradient p rovided as input by Abaqus. Several wo rked examples and case studies are provided for each section , so that the details and implications of the form of the deformation gradient can be fully understood. For each worked example in this tutorial paper the so urce files and code (Abaqus input files , UMATs, and Matlab script files) are provided, allowing the reader to efficiently explore the implications of the form of the deformation gradient in the development of a UMAT . Keywords : finite element analysis; a baqus; anisotropic; fibers * Correspond ing Author s: dnolan4@tcd.ie , Parsons Building, Trinity College Dublin, Dublin 2, Ireland ; **patrick.mcgarry@nuigalway.ie , Depa rtment of Biomedical Engineering, Alice Perry Engineering Building, NUI, Galway, Galway Ireland lallyca@tcd.ie , Parsons Building, Trinity College Dublin, Dublin 2, Ireland 2 Nomenclature 𝑮𝑖 global basis vector in the 𝑖 direction. 𝑬𝑖 material bas is vector in the reference configuration in the 𝑖 direction. 𝒆𝑖 material basis vector in the current configuration in the 𝑖 direction. 𝑰 identi ty matrix . 𝝌 deformation motion . 𝑿, 𝒙 position vector in the reference and current configuration respe ctively . 𝑭 deformation gradient tensor. 𝐹𝑖𝑗𝑬𝑬 components of the “global” deformation gradient in a global 𝑬𝑖⊗𝑬𝑗 basis system. 𝐹𝑖𝑗𝒆𝑬 components of a classical “local” deformation gradient in a local 𝒆𝑖⊗𝑬𝑗 basis system . 𝐹𝑖𝑗𝑎𝑙 components of an “Abaqus local” deformation gradient when a local coordinate system is specified with continuum elements. 𝐽 determinant of the deformation gradient 𝑹 rotational part of the deformation gradient , in co -rotational system us ed to map 𝑬𝑖 to 𝒆𝑖. 𝑼, 𝑽 right and left stretch tensors, respectively 𝑼𝑎𝑙 Abaqus definition of local right stretch tensor 𝑼̅ isochor ic part of the right stretch tensor 𝑪 right Cauchy -Green deformation tensor 𝑪′ classical definition of the l ocal right Cauchy -Green deformation tensor 𝑪𝑎𝑙 Abaqus definition of the local right Cauchy -Green deformation tensor 𝑩 left Cauchy -Green deform ation tensor 𝑩′ classical definition of the local left Cauchy - Green deformation tensor sdfsdfsdf 𝑩𝑎𝑙 Abaqus definition of the local left Cauchy - Green deformation tensor 𝑩̅ isochoric part of the left Cauchy -Green deformation tensor 𝜽, 𝑸, 𝚽, 𝝆 generic rotation matr ix 𝝈 Cauchy stress tensor 𝑨𝛽, fibre vector in the reference configuration, t he sub script 𝛽 denotes its basis system. 𝒂𝛽, fibre vector in the current configuration, the subscript 𝛽 denotes its basis system. 𝐼𝑓 fibre invariant, square of the stretch of the fibre 3 1 Introduction The implementation of new constitutive laws throug h user defined material su broutines (UMATs) in the finite element software Abaqus represents a powerful tool for researchers in the fields of solid mechanics, biomechanics, structural mechanics and material science. A UMAT is called at each integrati on poi nt with the deformation gr adient matrix provided as input. The user must compute the Cauchy stress tensor as the output from the UMAT (a consistent tangent matrix must also be provided for UMATS (Abaqus/Standard) , but not for VUMATs (Abaqus/Explicit) ). Whi le the logarithmic strain tensor is also provided as input to a UMAT, it should be noted in hyperelastic constitutive laws for large deformation kinematics that the Cauchy stress tensor is usually constructed from the deformation gradient. Examples i nclude constitutive laws for sof t tissue (Holzapfel et al., 2000; Weiss et al., 1996) , elastomers and rubber (Mooney, 1940; Rivlin, 1948) . It is often necessary or convenient to define a local material coordinate system (basis) who se orientation varies spatially across a body. For example, in a model of an aneurysm, it is useful to define a local ba sis system where the basis vectors are tangential to the aneurysm wall so that circumferential and axial stresses may be easily analysed . A lo cal coordinate basis co -rotates along with any rigid body rotation or rotational deformations as the material defo rms from the reference to the current configuration. This change of basis through out the simulation must be accounted for. A finite elem ent mo del in Abaqus is defined by an input file which contains all of the relevant information required to run the simul ation. An input file consists of several keywords which the Abaqus compiler recognises and interprets. A keyword is usually proceeded by numer ical data, for example the position of a node in 3D space, or a switching variable which turns a simulation featur e on or off. When viewing an input file, a keyword is easily identified by the asterisk character (*) which precedes it. A local coordin ate sy stem is specified in a n input file by using the *orientation keyword. If this keyword is not present then a global basis scheme is used throughout the simulation (Abaqus Analysis User’s Guide; 2.2.5 Orientations, 2019; Abaqus Keywords Reference Guide; *orientation, 2019) . In this paper we demonstrate that the form of the deformation gradient depends on the choice of coordinate basis system (i.e. local or global) and on the element type . The fo rm of the deformation gradient provided as input to an Abaqus user defined material subroutine has not previously been systematically outlined, to the best of our knowledge. In particular we demonstrate that if a local coordinate system is used in conjunct ion wi th continuum elements, the form of the deformation gradient provided as input to a UMAT is different to the classi cal definition of a local deformation gradient matrix. In Section 2 we outline the key differences between the class ical form of a local deformation gradient matrix and t hat provided as input to a UMAT. In Section 3 we then systematically outline the key steps that must be taken by the user in order to develop a UMAT 4 where the local Cauch y stre ss must be constructed from the deformation gradient using the e xample of an isotropic hyperelastic material. The use of anisotropic hyperelastic constitutive models based on structural tensors has become widespread in biomechanics in the past two decad es (Holzapfel et al., 200 0; Spe ncer, 1984; Weiss et al., 1996) . This class of constitutive mod el has been used to simulate soft tissues such as: arterial (Creane et al., 2012; Famaey et al., 2013; Ghasemi et al., 2018; Nolan and McGarry, 2015) , tendon (Khayyeri et al., 2016; Shearer, 2015) , cartilage (Nagel and Kell y, 2010) , skin (Annaidh et al., 2012) , myocardium (McEvoy et al., 2018) , annulus fibrosis (Eberlein et al., 2004) to name but a few . These constitutive models ar e based on fibre vectors which indicate the direction in which structurally important fibres are orientated , for example collagen or elastin fibres in arterial tissues (Gaul et al., 2017; Whelan et al., 201 9). This offers an attractive, intuitive method to define the anisotropy of a material. The findings of the current study are also relevant to several other fields of biomechanics modelling, including anisotropic deformation of bone behaviour (Feerick and McGarry, 2012) , stress fibre formation in cells (Reynolds et al., 2014) , and collagen alignment in tissue engineered scaffo lds (Wang et al., 2015) . Fibre vector based constitutive models often take advantage of local coordinate systems as it allows directional mechanica l propertie s determined in straightforward mechanical tests to be deployed in complex geometries. For example, one can calibrate a two -fibre hyperelastic model using biaxial and uniaxial mechanical testing data on arterial tissue in the circumferential and axial dire ction s. If one wish es to simulate an arch shaped arter y whose diameter is decreasing along its axial length, a local coordinate system will allow the use of this directionally calibrated model in this complex geometry. In Section 4 we demonstrate the key steps that must be taken by t he user to correctly account for fibre rotations in the construction of an anisotropic hyperelastic UMAT. The crea tion of advanced UMATs is not straightforward . One must have a strong command of solid mechanics and a comprehensive knowledge of how a cons titutive model must be implemented. Additionally, o ne must be familiar with Fortran , a coding language in which it is easy to make minor but consequential mistakes and is difficult to debug by modern standards. This tutorial paper has been created for thre e purposes: a) to introduce th e reader to the important differences in the continuum mechanics of a body where a global basis system or a local basis system is used, b) to highlight some issues we have found where a local basis system is used in conjunctio n with continuum elements in Abaqus, c) to provide materials which will act as a starting point to reader s who wish to create their own UMAT s. 1.1 Document structure This tutorial paper is intended as a step -by-step guide to developing a comprehensive understa nding of the construction of Abaqus UMATs in which the Cauchy stress is co mputed from the deformation gradient. 5 In Section 2 we examine the “Classical” forms of global and local deformation gradients. We then show that Abaqus do es not use the classical form of the local deformation gradient when continuum elements are used. We highlight the important implications this has for w riting user -defined material subroutines (UMATs) for isotropic hyperelastic material laws (Section 3), and anisotropic fibre - reinforced hyperelastic material laws (Secti on 4). Further case studies are presented in Section 5 highlighting the errors wh ich will occur if the local deformation gradient in Abaqus is not correctly managed . Several worked examples are provided for each section. The user is provided with all required code (Abaqus input files, UMATs, and Matlab script files) to systematically w ork through each example. 1.2 Acces s to Codes The codes created for this tutorial paper are stored and will be maintained at the GitHub site https://github.com/LallyLabTCD/localBa sisAbaqus/releases/tag/1.3 . Additionally, the codes have been archived at zenodo.org and have the DOI : https://doi.org/10.5281/zenodo.1968356 6 2 Lesson A: The deformation gradient in global and local coordinate systems Figure 1: 2-D schematic demonstrating the coordinate bases used in the reference and current configurations in Abaqus. The basis 𝑮𝑖 is the global basis, 𝑬𝑖 is the basis in the reference configuration at the material point defined by position vector 𝑿, and 𝒆𝑖 is the basis in the current configuration of the same material point which has undergone a motion 𝝌; this poin t is defined by position vector 𝒙. A) In the case where a local coordinate basis is used, this basis need not be orthogonal to the global basis in the reference configuration; i.e. 𝑬𝑖 need not be orthogonal to 𝑮𝑖. Importantly, under a motion 𝝌 the local basis co -rotates with any rigid bo dy rotations. B) The default setting in Abaqus where the global coordinate basis is used throughout the simulation. Both 𝑬𝑖 and 𝒆𝑖 are colinear to the global basis 𝑮𝑖. Figure 1 shows the typical kinematic setup in large deformation continuum mechanics where we consider a generic body in two configurations. The first is the reference configuration which is the starti ng position of the body. In this paper we assume that in the reference configuration the body is undeformed and holds no stress or strain (while this convention is widely adapted, it should be noted that it is not necessary that the reference configuration is a stress/strain free configuration) . Next, as shown in Figure 1, starting from the reference configuration a deformation 𝝌 is applied to the body , after which we consider the deformed body to be in the current configuration. In this paper we refer to three distinct coordinate systems ( orthogonal basis vector systems) which are shown in Figure 1. Firstly, there is the set of global basis vectors 𝑮𝑖 where 𝑖 is the number of spatial dimensions; this is typically a Cartesian coordinate system. Second ly, we describe a set of orthogo nal basis vectors at a material point in the reference configuration which we denote as 𝑬𝑖. Finally, we define a set of orthogona l basis vectors at a material point in the current configuration which we denote as 𝒆𝑖. Figure 1 illustrates the key difference between global and local coordinate systems in terms of the relationship between the basis vector systems and the material def ormation. When a local scheme is 7 used (Figure 1A), one may define a unique ly orientated local orthogonal material basis system 𝑬𝑖 at any point in the body in the reference configuration. This local basis system maps to the local basis system in the current configuration 𝒆𝑖 which may have its own unique orientation . In a co -rotational system, 𝒆𝑖 is determined by a rot ational transformation of 𝑬𝑖 based on the rotational deformation of the material . In the following sections we attempt to provide an in -depth explanation and understanding of the implications of such transformation of the basis vectors . In Abaqus, a local orthogonal material basis system 𝑬𝑖 is defined using the *orientation keyword . In the case of a global sche me (Figure 1B) the relationship is straightforward; material basis vectors in the reference and current configurations 𝑬𝑖 and 𝒆𝑖 are colinear with the global basis vectors 𝑮𝑖. Essentially, 𝑬𝑖=𝒆𝑖=𝑮𝑖 regardless of deformations. 2.1 Basic kinematics The deformation gradient is the fundamental unit of deformation in finite strain and hyperelastic continuum mechanics. It is defined as, 𝑭=𝜕𝒙 𝜕𝑿 , (1) where 𝑿 is a position vector in the reference configuration, and 𝒙 is a position vector in the current configuration (Abaqus Theory Guide; 1.4.1 Deformation, 2019) . The deformation gradi ent may be multiplicatively decomposed into two parts, one representing the r otational part of the deformation and the other representing the stretching part. This is known as a polar decomposition and is given as, 𝑭=𝑹∙𝑼=𝑽∙𝑹 (2) where the rotational part of the deformation gradient is the rotation matrix 𝑹, and the stretching part of the deformation with respect to (w.r.t) the reference configuration is the right stretch tensor 𝑼 and the stretching part of the deformation w.r.t. the current configuration is the left stretch tensor 𝑽.1 The right Cauchy -Green deformation tensor , commonly used in finite deformation continuum mechanics and hyperelasticity , is given as 𝑪=𝑭𝑇∙𝑭 (3) Using (2), this can be written in terms of the right stretch tensor , 1 Both 𝑼 and 𝑽 are symmetric tensors, and any rotation matrix 𝑹 has the property 𝑹−1=𝑹𝑇. 8 𝑪=𝑼𝑇∙𝑹𝑇∙𝑹∙𝑼=𝑼∙𝑼. (4) Similarly, the left Cauchy -Green deformation tensor, is given as 𝑩=𝑭∙𝑭𝑇 (5) Using (2), this can be written in terms of the right stretch tensor, 𝑩=𝑹∙𝑼∙𝑼𝑇∙𝑹𝑇=𝑹∙𝑪∙𝑹𝑇. (6) 2.2 Local and global basis vectors and defo rmation gradient components Consider a local scheme , as shown in Figure 1 A, whereby the basis vectors 𝒆𝑖 used to describe deformed material points in the current configuration are not the same as the basis vectors 𝑬𝑖 used to describe the corresp onding material points in the undeformed reference configuration (i.e. 𝒆𝑖≠ 𝑬𝑖). In a local scheme the basis vectors 𝑬𝑖 are fixed, while 𝒆𝑖 can change as the material deforms/rotates . Assuming that 𝑬𝑖=𝑮𝑖≠𝒆𝑖, the deformatio n gradient can be computed as follows (Bonet and Wood, 2008)2: 𝑭=∑𝜕𝑥𝑖𝒆 𝜕𝑋𝑗𝑬(𝒆𝑖⊗𝑬𝑗)3 𝑖,𝑗=1=∑ 𝐹𝑖𝑗𝒆𝑬(𝒆𝑖⊗𝑬𝑗)3 𝑖,𝑗=1; (7) where the superscripts are introduced to emphasise that the coordinates 𝑥𝑖𝒆 of the deformed position vector 𝒙 are w.r.t. the basis vectors 𝒆𝑖 and the coordinates 𝑋𝑖𝑬 of the undeformed position vector 𝑿 are w.r.t. the basis vectors 𝑬𝑖, whereby 𝒙=𝑥𝑖𝒆𝑒𝑖 and 𝑿=𝑋𝑖𝑬𝐸𝑖. It should be noted that the deformation gradient is a two -point tensor, given that it refers to both the reference configuration (with vectors 𝑬𝑖) and the deformed configuration (with vectors 𝒆𝑖). In a global scheme, a s shown in Figure 1B, the same basis system is used for the reference and current configurations, i.e. 𝑬𝑖=𝒆𝑖=𝑮𝑖, such that 𝒙=𝑥𝑖𝐸𝐸𝑖 and 𝑿=𝑋𝑖𝐸𝐸𝑖. In a global scheme 𝒆𝑖 (=𝑬𝑖) is fixed and does not rotate with the material , in contrast to a local schem e. The deformation gradient can be computed as , 𝑭=∑𝜕𝑥𝑖𝑬 𝜕𝑋𝑗𝑬(𝑬𝑖⊗𝑬𝑗)3 𝑖,𝑗=1=∑ 𝐹𝑖𝑗𝑬𝑬(𝑬𝑖⊗𝑬𝑗)3 𝑖,𝑗=1; (8) Noting that the dyadic product 𝑬𝑖⊗𝑬𝑗 has a value of one at the ( 𝑖, 𝑗) component and zero everywhere else, the components of 𝑭 in a global scheme are given as 𝜕𝑥𝑖𝑬 𝜕𝑋𝑗𝑬 . 2Unless otherwise stated we assume th at 𝑬𝑖=𝑮𝑖. In Section 2.6 we consider the special case where 𝑬𝑖≠𝑮𝑖. 9 Returning again to the general case of a local system where 𝒆𝑖≠𝑬𝑖 (as illustrated in Figure 1B ), the basis vector s of the current and referenc e configuration can be related through a standard transformation such that 𝒆𝑖=𝜽∙𝑬𝑖 or in index notation, 𝒆𝑖=∑𝜃𝑘𝑖3 𝑘=1𝑬𝑘 (9) Substituting into (8) we can relate the deformation gradient components in the local 𝒆𝑖⊗𝑬𝑗 basis system and the global 𝑬𝑖⊗𝑬𝑗 system , 𝑭=∑𝜕𝑥𝑖𝒆 𝜕𝑋𝑗𝑬(𝒆𝑖⊗𝑬𝑗)3 𝑖,𝑗=1=∑𝜕𝑥𝑖𝒆 𝜕𝑋𝑗𝑬𝜃𝑘𝑖𝑬𝑘⊗𝑬𝑗=∑𝜕𝑥𝑖𝑬 𝜕𝑋𝑗𝑬(𝑬𝑖⊗𝑬𝑗)3 𝑖,𝑗=13 𝑖,𝑗,𝑘=1 (10) Expanding the third and fourth expressions in equation (10) above it is straightforward to show (see Appendix A. 1) that, ∑𝜃𝑖𝑘3 𝑘=1𝜕𝑥𝑘𝒆 𝜕𝑋𝑗𝑬=𝜕𝑥𝑖𝑬 𝜕𝑋𝑗𝑬 (11) 𝐹𝑖𝑗𝑬𝑬=𝜃𝑖𝑘𝐹𝑘𝑗𝒆𝑬 (12) If we drop the 𝑖,𝑗,𝑘 index notation we can write (12) using tensor notation 𝑭𝑬𝑬=𝜽∙𝑭𝒆𝑬 (13) Or the reverse relationship , 𝑭𝒆𝑬=𝜽𝑇∙𝑭𝑬𝑬 (14) And it should also be noted that , 𝑭𝑬𝑬=𝑭 (15) A key point to note from equ ation (13) is that , because the deformation gradient is a two -point tensor , forme d from the reference configuration (with basis vector s 𝑬𝑖) and the current configuration (with basis vectors 𝒆𝑖), only one operation is requ ired to transform 𝑭𝑬𝑬 to 𝑭𝒆𝑬. In other words, it is clearly not correct to state that transformation of the deformation gradient 10 requires two operations, i.e. 𝑭𝒆𝑬≠𝜽𝑇∙𝑭𝑬𝑬∙𝜽, in contrast to a conventional second order tensor 𝑺 (e.g. Cauchy stress) where the standard transformation is given as 𝑺𝒆𝒆=𝜽𝑇∙𝑺𝑬𝑬∙𝜽. We shall return to this important point in Section 2.3. 2.2.1 Co-rotational local system: In a co-rotational system , the rotation of the basis vectors is given by 𝑹 from the polar decom position 𝑭=𝑹∙𝑼, the current basis vector is given by, 𝒆𝑖=𝑹∙𝑬𝑖 (16) Following from (9) and (13) above, we ca n write , 𝑭𝑬𝑬=𝑹∙𝑭𝒆𝑬 (17) Noting that the rotation matrix 𝑹 is orthogonal, we can write t he co-rotational local deformation gradient 𝑭𝒆𝑬 as 𝑭𝒆𝑬=𝑹𝑇∙𝑭𝑬𝑬=𝑹𝑇∙𝑹∙𝑼=𝑼 (18) This demonstrates that in a co -rotational system the local deformation gradient is the Cauchy -Green stretch tensor. This result is expected, given that a co-rotational system is not affected by rigid body rotations. The right and left Cauchy -Green deformation tensors can be written in terms of the local basis system using the definitions in (3), (5), (18) and the polar decomposition in (2), 𝑩𝒆𝒆=𝑭𝒆𝑬∙(𝑭𝒆𝑬)𝑇 (19) 𝑩𝒆𝒆=𝑹𝑇∙𝑭𝑬𝑬∙(𝑭𝑬𝑬)𝑇∙𝑹 𝑩𝒆𝒆=𝑹𝑇∙𝑹∙𝑼∙𝑼𝑇∙𝑹𝑇∙𝑹=𝑼∙𝑼 where once again the superscripts indicate the basis vectors of the tensor . Similarly, 𝑪𝑬𝑬=(𝑭𝒆𝑬)𝑇∙𝑭𝒆𝑬=𝑼∙𝑼 (20) Therefore , 𝑪𝑬𝑬=𝑩𝒆𝒆 (21) 11 2.3 Abaqus/ Standard form of the local deformation gradient for continuum elements in a local basis syste m When continuum elements are used in Abaqus /Standard (Abaqus Analysis User’s Guide; 28.1 General Purpose Continuum Elements, 2019) in conjunction with a user defined material subroutine (UMAT), the deformat ion gradient and other tensors are provided as input to the UMAT are defined w.r.t. the global basis system , provided that a local basis system is not defined through a *orientation statement3. Essentially, at all material points 𝑬𝑖=𝒆𝑖=𝑮𝑖, as shown in Figure 1A. Importantly, the deformation gradient components that are input to the UMAT are correctly given as 𝜕𝑥𝑖𝑬 𝜕𝑋𝑗𝑬 and are fully consistent with t he expression for 𝑭 given in equation (8). In Abaqus one may define a local coordinate frame where the basis vectors in the reference configuration , 𝑬𝑖, are directly specified by the user. This feature is invoked when the *orientation keyword is use d. In such a scheme , the basis vectors in the current configuration 𝒆𝑖 co-rotate with the material, as described in equation (18). We denote 𝑭𝑎𝑙 as the deformation gradient that is passed to a UMAT when such a local coordinate system is used in conjunction with continuum elements in Abaqus/Standard. The superscript “ al” denotes “abaqus local”. As we will describe in the following, 𝑭𝑎𝑙 does not correspond to the classical form of a global (𝑭𝑬𝑬, equation (8)) or local (𝑭𝒆𝑬, equ ation (11) ) deformation gradient. Instead, 𝑭𝑎𝑙 is related to the global deformation gradient 𝑭 (in the 𝑬𝑖⊗𝑬𝑗) basis system through the following operation: 𝑭𝑎𝑙=𝑹𝑇∙𝑭∙𝑹 (22) As discussed in Section 2.2 above, such a transformation is only appropriate for a single basis tensor (Cauchy stress is an example of such a tensor ) whose basis is to be changed from 𝑬𝑖⊗𝑬𝑗 to 𝒆𝑖⊗ 𝒆𝑗, rather than a two -point tensor where the basis is given as 𝒆𝑖⊗𝑬𝑗. Recalling equation (14), we once again emphasise that the correct transformation of the deformation gradient from a global to a local co -rotational scheme is given as 𝑭′=𝑹𝑇∙𝑭=𝑼. The unconventional f orm of 𝑭𝑎𝑙 has important consequences for the calculation of deformation tensors and for fibre vector kinematics. Consider, for example, the Cauchy -Green deformation tensors, 𝑪𝑎𝑙=(𝑹𝑇∙𝑭∙𝑹)𝑇∙ (𝑹𝑇∙𝑭∙𝑹) (23) 𝑪𝑎𝑙=𝑹𝑇∙𝑭𝑇∙𝑭∙𝑹=𝑹𝑇∙𝑪∙𝑹=𝑹𝑇∙𝑼∙𝑼∙𝑹 and 3 In an Abaqus/Standard UMAT the deformation gradient components are stored under the variable names DFGRAD1 (current increment) and DFGRAD0 (previous increment). 12 𝑩𝑎𝑙=(𝑹𝑇∙𝑭∙𝑹)∙(𝑹𝑇∙𝑭∙𝑹)𝑻 (24) 𝑩𝑎𝑙=𝑹𝑇∙𝑭∙𝑭𝑇∙𝑹=𝑹𝑇∙𝑩∙𝑹 so that 𝑩𝑎𝑙=𝑼∙𝑼𝑇=𝑼∙𝑼 (25) Therefore, in contrast with the classical definition of the local deformation tensors in (21), in Abaqus , 𝑪=𝑩𝑎𝑙≠𝑪𝑎𝑙 (26) This result has important consequences for the i mplementation of hyperelastic constitutive models and it is crucial that mechan icians writing user -defined materials in Abaqus are aware of it. The difference between 𝑭′ and 𝑭𝑎𝑙 also has important consequences for the use of the Abaqu s utility subroutine sprind, as discussed in Appendix A.2. 2.3.1 Relating the classical local deformation gradient 𝑭𝒆𝑬 to the abaqus local deformation gradient 𝑭𝑎𝑙 Let us probe a little further into (22), using (2), 𝑭𝑎𝑙=𝑹𝑇∙𝑭∙𝑹=𝑹𝑇∙𝑹∙𝑼∙𝑹=𝑼∙𝑹 (27) Contrast this with the classical expression for the local deformation gradient in (18). Though (27) is similar to the polar decomposition in outlined in (2), the order ing of 𝑼 and 𝑹 is incorrect. A relationship for the local right stretch tensor can also be established using (27), 𝑭𝑎𝑙=𝑹∙𝑼𝑎𝑙=𝑼∙𝑹 (28) 𝑼𝑎𝑙=𝑹𝑇∙𝑼∙𝑹 (Note that no such relationship exists in the classical description of local basis system , whereby 𝑭𝒆𝑬=𝑼) One straightforward approach to writing an isot ropic hyperelastic UMAT in Abaqus/Standard when a local orientation system is used in conjunction with continuum elements is as follows: (i) Perform a polar decomposition of 𝑭𝑎𝑙 (which is provided as input to the UMAT) in order to determine 𝑹. (ii) Compute 𝑭𝒆𝑬=𝑼=𝑭𝑎𝑙∙𝑹𝑇. (iii) From 𝑭𝒆𝑬(=𝑼), compute the Cauchy stress in the deformed local coordinate system. 13 We will now step through some examples which will underscore the theo retical aspects outlined in this section. 2.4 Worked Example # 1: Struct ural Elements versus Continuum Elements We will now go through a worked example to illustrate the implications of the u se of global and local coordinates in Abaqus. Complementary Abaqus an d Matlab files are included in the supplementary data in the folder called “ Lesson A - deformation gradient ” which will aid in understanding this example. We first highlight the differenc e between the treatment of the deformation gradient in structural elements (Abaqus Analysis User’s Guide; 29.6 Shell E lements, 2019) such as a shell element , and in continuum elements. Specifically, we compare a four noded shell element to a four noded plane stress element, given that the out of plane stress is zero for both element types. For direct comparison with a p lane stress element, the shell element is constructed in the x -y plane and is constrained to deform only in thi s plane. A uniaxial stretch with a 45o rotation is imposed on the element, as shown in Figure 2. In the case of the mem brane the following components of the deformatio n gradient is passe d to a UMAT, regardless of whether an orientation is prescribed in the membrane section: 𝑭𝒆𝑬=[1 0 0 01.1 0 0 0 0.9091] (29) Clearly this deformation gradient matrix corresponds to the right stretch tensor 𝑼, and is correct for a local co -rotational 𝒆𝑖⊗𝑬𝑗 basis system . In other words, for membrane elements a default local coordinate system that rotates with the element is always used, so that the resultant local deformation gradient does not contain any information on material rotation. It should be noted that the valu e of 𝐹33=0.9091 is simply a consequence of the default implementation of material incompressibility . 14 Figure 2 Uniaxial stretch of 10% accompanied by a 45o rotation. This deformation field is applied to membrane and shell elemen t, and to and plane stress a nd plane strain elements. In contrast, for a plane stress element, if *orientation is not prescribed and an orientation is not prescribed in the solid section, a global scheme is implemented and the following global deformation gradient matrix is passed t o a UMAT: 𝑭=𝑭𝑬𝑬=[0.7071 0.7778 0 −0.7071 0.7778 0 0 0 0.9091] (30) A polar decomposition yields , 𝑼=[1.0 0 0 1.1 0 0 0 0.9091] (31) And 𝑹=[0.7071 0.7071 0 −0.7071 0.7071 0 0 0 1] (32) However, if a local orientation is specified in the solid section state ment, such that the *orientation statements prescribes a local coordinate system that originally coincides with the global system, the fol lowing deforma tion gradient is passed to a UMAT: 𝑭𝑎𝑙=[0.7071 0.7071 0 −0.7778 0.7778 0 0 0 0.9091] (33) As outlined previously, this form of the deformation gradient is related to the global deformation gradient through the following operation: 𝑭𝑎𝑙=𝑹𝑇∙𝑭∙𝑹 (34) It should be noted that all of the results for plane stress elements outlined above can also be obtained for plane strain elements, with the exception of the 𝐹33 component, which is equal to 1 due to the imposition of zero out -of-plane strain. • The Matlab script Worked_example_ 1.m demonstrates the calculation of the Classical and Abaqus definition of the local deformation gradients outlined in Sections 2.2 and 2.3. These correspond to the local deformation gradient when either a structural or continuum element is used respectively. • The deformation gradient tensor is calculated by a continuum element in Abaqus/Standard with and without the *orientation keyword are given in input files WE1_PStress_Ori .inp and WE1_PStress_NoOri .inp . 15 • The deformation gradient tensor is calculated by a str uctural element (shell) in Abaqus/Standard with and without the *orientation keyword are given in input files WE1_Shell_Ori.inp and WE1_Shell_NoOri.inp . 2.4.1 Large rotations and large deformations: We next consider a uniaxial stretch to a higher nominal strain of 100%, again accompanied by a 45o rigid body rotation, as shown in Figure 3A. Figure 3 (A) Uniaxial stretch of 100% accompanied by a 45∘ rotation. (B) Uniaxial stretch of 100% without rotation. As expected, a membrane or shell element return the correct local deformation gradient of 𝑭𝒆𝑬=[10 0 02 0 000.5] (35) regardless of whether *orientation is specified in the shell or membrane section statement. Similarly, in the case of a plane stress continuum element, it a local orientation is not specified in the solid section statement, Abaqus p asses the fol lowing global deformation to a UMAT: 𝑭=𝑭𝑬𝑬=[0.7071 1.414 2 0 −0.7071 1.414 2 0 0 0 0.5] (36) 16 Clearly a polar decomposition such that 𝑭=𝑹∙𝑼 yields 𝑼=[10 0 02 0 000.5] (37) And 𝑹=[0.7071 0.7071 0 −0.7071 0.7071 0 0 0 1] (38) However, if a local orientation is specified, the following deformation gradient is passed to a UMAT 𝑭𝑎𝑙=[0.7312 0.6846 0 −1.4367 1.3902 0 0 0 0.5] (39) In this case, some errors are introd uced such that 𝑭𝑎𝑙≠𝑹𝑇∙𝑭∙𝑹. The magnitude of the error can be determined from the relationship 𝑭𝑎𝑙=𝑼𝑎𝑙∙𝑹, so that the inaccurate stretch matrix is given as: 𝑼𝑎𝑙=[1.0011 −0.0329 0 −0.0329 1.9989 0 0 0 0.5] (40) The following additional points should be noted: (i) Such inaccurate values for 𝑭𝑎𝑙 (and, as a result, for 𝑼𝑎𝑙) are provided by Abaqus for plane stress, plane strain and 3D continuum elements when a local orientation is specifi ed and large deformations are accompanied by large rotations . Such inaccuracies are not observed for large deformations without rotations (e.g. for the deform ation shown in Figure 3B, 𝑭𝑎𝑙=[1 0;0 2] as expected). (ii) Despit e the inaccuracy in the values of 𝑭𝑎𝑙, accurate values of log strain and nominal strain are provi ded by Abaqus. In the example above, accurate values of maximum principal log strain (0.6931) and maximum principal nominal strain (1.0) are obtained, re gardless of the choice of element or global/local basis system. For the accom panying files review the deformation given in (35): • The Matlab script Worked_exam ple_1_largeDef .m demonstrates the calculation of the Classical a nd Abaqus definition of the local deformation gradients outlined in Sections 2.2 and 2.3. These correspond to the local deformation gradient when either a structural or continuum element is used respe ctively. • The deformation gradient tensor is calculated by a continuu m element in Abaqus/Standard with the *orientation keyword is given in input file WE1_PStress_Ori _largeDef .inp . 17 • The deformation gradient tensor is calculated by a structural (shell) elemen t in Abaqu s/Standard with the *orientation keyword is given in input file WE1_Shell_Ori_largeDef .inp . 18 Figure 4 Schematic outlining the reference and current configuration of a unit cube highlighting the relevant basis systems f or Worked E xample #2. The global basis vectors are 𝑮𝑖, and in this particular scenario the material basis vectors in the reference config uration are colinear with the global basis, 𝑬𝑖=𝑮𝑖. The basis vectors in the current configuration 𝒆𝑖 are shown at each of the integration points in the unit cube. 2.5 Worked Example # 2: Abaqus local deformation gradient in a 3 -D continuum element Consi der a unit cube which is subjected to a motion 𝝌 which entails a rotation component, and is described in the global coordinate system by the deformation gradient, 𝑭𝑬𝑬=[1.10 0.10 0.00 0.05 0.90 0.15 0.20 0.00 1.20] (41) We now wish to calculate the deformation gradient when a local basis is used. Assuming the local basis i s aligned with th e global basis in the reference configuration , using (18) the following is obtained4, 𝑹=[0.995 0.035 −0.089 −0.029 0.997 0.077 0.091 −0.074 0.993] (42) 𝑭𝑎𝑙=[1.114 0.105 0.011 0.053 0.896 0.123 0.212 −0.025 1.190] (43) To verify that this is indeed the deformation gradient that Abaqus uses when a local coordinate system is invoked, we created a simple FE model of a single 3D unit cube element with full 4 All presented values have been rounded to 3 decimal places. 19 integration being deformed according to deformation gradient in (41). The displacement at each node is calculated using 𝒖=𝑭∙𝑿−𝑿, where 𝒖 is the displacement of a node whose position in the reference configuration is giv en by 𝑿. Figure 4 shows the global basis 𝑮𝑖 which in this case is equal to the local material basis in the reference configuration 𝑬𝑖 in the unit cube. Also shown is the rotated local material basis 𝒆𝑖 in the current configuration . Figure 5 shows the values of t he deformation gradient output to the Abaqus .dat file (see Appendices for how to output the deformation gradient). When a global coordinate system is used we obtain the global defo rmation gradient given in (41), and when a local coordinate system is used we get the result (43). • The Matlab script Worked_example_ 2.m demonstrates the calculation of the global and local deformation gra dients outlined a bove. • The global deformation gradient tensor is calculated by the Abaqus/Standard input file WE2_global.inp . • The local deformation gradient tensor is calculated by the Abaqus/Standard input file WE2_local.inp . • The post -processing file plotLocalCoordinates. py plots the material basis vectors in resulting output database (.odb) file from WE2_local.inp . Figure 5 Output in the .dat file of the deformation gradient from the Abaqus FE simulation of a cube under deformat ion using A) a global coordinate system, and B) a local coordinate system. The deformation gradient output in ( A) correspon ds to (41) The deformation gradient output in ( B) corresponds to (43). 20 2.6 Worked Example # 3: Abaqus local deformation gradient in a 3 -D continuum element where 𝑬𝑖≠𝑮𝑖 Figure 6 Schematic for Worked Example # 3 showing a unit cube in ; A) the reference configuration with basis vectors 𝑬𝑖 rotated 20∘ about the global 𝑮3 basis vector, and B) the current configuration showi ng the co - rotated basis vectors 𝒆𝑖. Next , we consider the same deformation as Worked Example # 2, but in this case the reference configuration local coordinate basis is not the same as the global coordinate basis, i.e. 𝑬𝑖≠𝑮𝑖, instead 𝑬1 and 𝑬2 are rotated clockwise by 20∘ about the 𝑮3 axis. Starting with the global deformation gradient 𝑭 given in (41) and using a standard rotation matrix 𝑸, 𝑸=[cos(20∘)−sin(20∘)0 sin(20∘) cos(20∘) 0 0 0 1] (44) we perform a change of basis of (41) from 𝑮𝑖 to the new 𝑬𝑖, 𝑭20𝑬𝑬=𝑸𝑇∙𝑭𝑮𝑮∙𝑸 (45) where the subscript acts as a remin der that the reference basis vectors are rotated by 20∘. The reader w ill recognise that the transformation in equation (42) above is similar to that for a single point tensor. This is the case because we have merely transformed th e global basis vector 𝑮𝑖⊗𝑮𝑗 21 to 𝑬𝑖⊗𝑬𝑗. By pe rforming a polar decomposition o n this deformation gradient, 𝑭20=𝑹20𝑼20, we can calculate the local deformation gradient in this new basis system5, 𝑭20𝑎𝑙=(𝑹20)𝑇∙𝑭20𝑬𝑬∙𝑹20 (46) The local deformation gradient, 𝑭20𝑎𝑙, is passed into the UMAT subroutine by Abaqus when a local material basis is used. Abaqus does not provide 𝑭𝑬𝑬 or 𝑭20𝑬𝑬. • The Matlab script Worked_example_ 3.m demonstrates the calculation of the local deformation gradient in the rotated reference material bas is outlined above. • The c orresponding Abaqus simulation is performed by running the file WE3_local.inp . On inspection of the resulting .dat file from this simulation, one will observe that the deformation gradie nt is identical to (46). In the two worked examples above , we start with a global deformation gradient and map to a local deformation gradient. It is important to n ote that the inverse of this process is achievable, one can start with a local deformation gradient and map back to a global deformation gradient . 5 On running the Matlab code ( Worked_example_3.m ) you will notice that 𝑹20=𝑹; this is because the rotation between 𝑬𝑖 and 𝒆𝑖 is the same in both of these examples. 22 3 Lesson B: The calculation of stress in an isotropic hyperelastic material Hyperelastic materia l models are i deal for use in applications where there are large deformations. Unlike linear elastic or hypoelastic materials, the stress state in a hyperelastic material is a unique function of its deformation, i.e. for a given 𝑭 there is a unique 𝝈. Consequently, for a closed -loop deformation, the total work is zero and the starting stress and energy states are recovered. This thermodynamic consistency is evident in the objectivity of stress in hyperelastic materials. This is where the stress state in material is independent of the chosen coordinate system, and no erroneous stresses or stress rates will be computed when the material undergoes rigid body rotation. Practically speaking, one of the easiest ways to ensure that stre ss objectivity has been preserved when validating a UMAT is to deform a body, then check that the stress invariants6 remain constant under a subsequent rigid body rotation. When a user -defined material (UMAT) is used in Abaqus /Standard , the stress which i s returned to the FE solver must be given w.r.t. the coordinate basis of that solution increment . As we have seen in the previous lesson , simulations can be posed using a fixed global coordinate basis or a co - rotat ional local coordinate basis. Regardless o f the choice of basis system, for the application of a displacement field 𝝌 must result in the calculation of a unique stress state (in terms of stress invariants) and a unique strain energy density . In this lesson we will look at the calculation of stress in an isotropic neo-Hookean hyperelastic material . Here, the stress in the global basis is given as, 𝝈NH=2 𝐽𝐶10(𝑩̅−tr(𝑩̅) 3𝑰)+2 𝐷1(𝐽−1)𝑰 (47) where 𝐽 is the determinant of the deformation gradient 𝑭, the left Cauchy -Green tensor 𝑩 is given in (5) and its isochoric part 𝑩̅=𝐽−23⁄ 𝑩. The material parameters 𝐶10 and 𝐷1 control the isochoric and volumetric response of the stress respectively. The operator tr(∙) is the trace of a tensor , which in this context is the sum of the diagona l elements. Let’s now consider the construction of the global 𝑩 tensor in terms of the Abaqus local deformation gradient 𝑭𝑎𝑙 that is input to a UMAT when a local orientation is specified for a continuum element. Recalling the definition of 𝑭𝑎𝑙 from Equation (22), the Abaqus local form of the left Cauchy -Green tensor, 𝑩𝑎𝑙, is related to the global 𝑩 tensor as follows: 𝑩=𝑭∙𝑭𝑇= (𝑹∙𝑭𝑎𝑙∙𝑹𝑇)∙(𝑹∙𝑭𝑎𝑙∙𝑹𝑇)𝑇 (48) 6 For example the pri ncipal stresses, von Mises stress, and pressure stress are all stress invariants 23 = 𝑹∙𝑭𝑎𝑙∙𝑹𝑇∙𝑹∙𝑭𝑎𝑙𝑇∙𝑹𝑇 = 𝑹∙𝑭𝑎𝑙∙𝑭𝑎𝑙𝑇∙𝑹𝑇 = 𝑹∙𝑩𝑎𝑙∙𝑹𝑇 So, the left Cauchy -Green tensor 𝑩 can be mapped from a global to a local coordinate basis in a similar manner as the Cauchy stress is. The isochoric left Cauchy -Green tensor 𝑩̅ can also be mapped from a global to a local coordinate basis in this manner , 𝑩̅=𝑹∙𝑩̅𝑎𝑙∙𝑹𝑇. (49) If we consider the expression for neo -Hookean stress in (47) in the global basis and substitute for (49) for 𝑩̅ using (34) we obtain the global stress tensor in terms of the Abaqus local deformation tensor 𝑩̅𝑎𝑙, 𝝈NH= 2 𝐽𝐶10(𝑹∙𝑩̅𝑎𝑙∙𝑹𝑇−𝑹∙tr(𝑩̅𝑎𝑙) 3∙𝑹𝑇)+2 𝐷1(𝐽−1)𝑹∙𝑰∙𝑹𝑇 (50) = 𝑹∙{2 𝐽𝐶10(𝑩̅𝑎𝑙−tr(𝑩̅𝑎𝑙) 3𝑰)+2 𝐷1(𝐽−1)𝑰}∙𝑹𝑇 = 𝑹∙{𝝈NH𝑎𝑙}∙𝑹𝑇 The term in the brackets is the local stress which simply follows the rules for the change of basis of a second order tensor. In the context of a UMAT, 𝝈NH𝑎𝑙 is to be returned to the solver. Another alternative perspective is to consider the polar decomposition of 𝑭=𝑹∙𝑼, where in much the same way as (48)7, 𝑩̅=𝑹∙𝑼̅2∙𝑹𝑇 where 𝑼̅=𝐽−1/3𝑼. Again, if we consider the expression for neo-Hookean stress in the global coordinate frame, 𝝈NH= 2 𝐽𝐶10(𝑹∙𝑼̅2∙𝑹𝑇−𝑹∙tr(𝑼̅2) 3∙𝑹𝑇∙𝑰)+2 𝐷1(𝐽−1)𝑹∙𝑰∙𝑹𝑇 (51) = 𝑹∙{2 𝐽𝐶10(𝑼̅2−tr(𝑼̅2) 3𝑰)+2 𝐷1(𝐽−1)𝑰}∙𝑹𝑇 = 𝑹∙{𝝈NHcr}∙𝑹𝑇 The term in brackets is called the co -rotational stress 𝝈NHcr. This is the stress th at should be returned to the solver for all elements when using a user -material in Abaqus/Explicit (i.e. via a VUMAT) , or for 7 Note th at as 𝑼 is symmetric, 𝑼∙𝑼𝑇=𝑼∙𝑼=𝑼2 24 membrane and shell elements in Abaqus/Standar d (i.e. via a UMAT) regardless of whether a local basis is defined in the model . All stress and strain quantities are defined w.r.t. the co -rotational basis. The use of a corotational stress is particularly useful in dynamic and rate dependent simulation as it avoids the computation of erroneous stress rates. It is noted that the stress i n the Abaqus local coordinate basis, (50), is the same as the co -rotational stress, (51): i.e. 𝝈NH𝑎𝑙=𝝈NHcr. This is because the local coordi nate basis in (50) is the same as a co - rotational basis. 3.1 Worked Example #4: Calculation of the local stress in a continuum element consisting of an isotropic material We will now go through some of the principles outlined in Lesson B and will calculate the neo - Hookean stress in a m aterial using a global and local basis. We start with the global deformation gradient 𝑭 shown in Worked Example # 2 and given in (41). Equation (47) is used to calculate the global Cauchy stress tensor using 𝑩̅, with material properties 𝐶10=0.2 MPa and 𝐷1=2.0 MPa-1. The resulting stress tensor is 𝝈NH=[0.198 0.044 0.066 0.044 0.082 0.057 0.066 0.057 0.276] (52) Using the Abaqus local Cauchy -Green deformation tensor 𝑩𝑎𝑙 (48) in the expressi on for neo - Hookean stress given in (47) we arrive at the expression for the Abaqus local stress tensor , 𝝈NH𝑎𝑙=[0.207 0.047 0.075 0.047 0.077 0.041 0.075 0.041 0.271] (53) Additionally, the re sult in (53) is calculated when 𝑼̅2 is used in place of 𝑩̅𝑎𝑙 in (47). This is the co - rotational stress. It can a lso be shown that 𝝈NH𝑎𝑙 can be directly determined from 𝝈NH using the relationship given in (50). • The Matlab script Worked_example_ 4.m demonstrates the calculation of the global, local, and co -rotational stress tensors outlined above. • The global stress tensor is calculated by the Abaqus /Standard input file WE4_global.inp . • The local stress tensor is calculated by the Abaqus /Standard input file WE4_local.inp . 3.2 Worked Example # 5: Calculation of the local stress in a continuum element consisting of an isotropic material where 𝑬𝑖≠𝑮𝑖 We now consider the calculation of stress using a neo-Hookean model for the scenario in Worked Example # 3. Here the local basis is rotated by 20∘ about the 𝑮3 axis in the reference configuration (refer back to Figure 6), i.e. 𝑬𝑖≠𝑮𝑖. As we are using the same deformation as Worked Example # 4 25 the global stress is the same as that given in (52). However , in the case of a local basis, we use the deformation gradient 𝑭20𝑎𝑙 given in (46) to calculate the stress, resulting in, 𝝈NH𝑎𝑙,20=[0.222 0.006 0.084 0.006 0.062 0.013 0.084 0.013 0.271] (54) The result for the local stress in (54) is also arrived at by using the right stretch tensor 𝑼20, determined from the polar decomposition of the global deformation gradient in the rotated basis 𝑭20𝑬𝑬=𝑹20𝑼20, instead of the deformation gradient for the calculation of stress. This is the co - rotational stress and, as noted previously, upon comparison of (50) and (51) the local stress and co - rotational stress are identical, provided they share the same coordinate basis in the refere nce configuration. To map the stress calculated in this r otated local basis to the global stress calculated in (52), we may use the transformation, 𝝈NH=𝑸∙𝑹20∙𝝈NH𝑎𝑙,20∙𝑹20𝑇∙𝑸𝑇. (55) where was given prev iously in (44). • The Matlab script Worked_example_ 5.m demonstrates the calculation of the local, and co-rotational stress tensors in the rotated reference basis outlined above. • The local stress tensor is calculated by the Abaqus/Standard input file WE5_local.inp . 26 4 Lesson C: The calculation of stress in an anisotropic hyperelastic material Lastly, we examine the fibre kinemati cs and calculation of stress in global and local coordi nate bases for a fibre -based anisotropic hyperelastic material. In this scenario, one must take care that the fibre orientation vector is defined w.r.t. the appropriate coordinate basis system . Figur e 7: Schematic outlining the key concepts of fibre vector kinematics in a global and local coordinate frame. A) Global coord inate frame. As we move from the reference to the current configuration, the basis system does not rotate, 𝒆𝑖=𝑬𝑖. The change of the fibre vector in the reference configuration 𝑨 to the current configuration 𝒂 is calculated from the global defor mation gradient 𝑭. B) Local coordinate frame. In this case the basis system co-rotates as per (16). The local deformation gradient is (partially) defined w.r.t. this rotated coordinate frame and so too should the deformed fibre v ector 𝒂. This calculated using (58). 4.1 Fibre K inematics The orientation of a fibre in the reference configuration is defined by a unit vector 𝑨𝛽, where we use a capital 𝑨 to denote that the fibre vector is in the reference configuration, and the subscript 𝛽 to denote its basis system . As the ma terial is deformed the fibre vector may stretch and rotate in an affine manner such that 𝒂=𝑭∙𝑨, (56) 27 where lower case 𝒂 is the fibre vector in the current configuration. Let us consider the deformat ion of the fibre vector when a global basis system is used , as i llustrated in Figur e 7A. The global basis is used in both the reference and current configuration s; i.e. 𝑬𝑖=𝑮𝑖. Hence the fibre vector in the curren t configuration (56) becomes 𝒂𝑬𝑖=𝑭𝑬𝑬∙𝑨𝑬𝑖. (57) In the case of a global scheme in Abaqus/Standard, t he user typically defines the fibre direction unit vector 𝑨𝑬𝑖 in the global basis , which is typically provided as an input to a UMAT through the user specified material properties . This straightforward approach to the specification of fibre directions is widely used in the literature on fibre based anisotropic hyperel asticity (Gasser et al., 2006; Guo et al., 2006; Holzap fel, 2000; Holzapfel et al., 2000; Weiss et al., 1996) . Figur e 7B illustrates the case of a local basis system. In this case the user typically prescribes an undeformed fibre vector 𝑨𝑬𝑖 in the reference local coordinate system, 𝑬𝑖. The local basis sys tem co- rotates with the material, again with the rotated basis system denoted as 𝒆𝑖. Of course, the affine assumption enforces that the fibre also rotat es with the mater ial, The deformed fibre vector is the n simply obtained using the local deformation gradient such that 𝒂𝒆𝑖=𝑭𝒆𝑬∙𝑨𝑬𝑖. (58) Recall also that 𝑭𝒆𝑬=𝑼. As previously mentioned, in the case of structural elements in Abaqus/Standard, and in the case of a ll element types in Abaqus explicit, the user is provided with the classical form of the local deformation gradient 𝑭𝒆𝑬=𝑼. From this, 𝒂𝒆𝒊 can be directly calculated . However , in Abaqus /Stand ard for the case of a continuum element used with a local basis syste m (through a *orientation statement), a UMAT provides 𝑭𝑎𝑙(=𝑹𝑇∙𝑭∙𝑹=𝑼∙𝑹) as input, so that the following steps must now be taken to compute 𝒂𝒆𝑖: (i) Perform a polar deco mposition of 𝑭𝑎𝑙 (which is provided as input to the UMAT) in or der to determine 𝑹. (ii) Compute 𝑼=𝑭𝑎𝑙∙𝑹𝑇 (=𝑭′). (iii) Compute 𝒂𝒆𝑖=𝑼∙𝑨𝑬𝑖. A key component in anisotropic hyperelasticity is the 2nd order structural tensor forme d by the operation 𝒂⊗𝒂. This operation is written in index notation as, 𝒂⊗𝒂=𝑎𝑖𝑎𝑗, and may be more simply stated in matrix notation as, 𝒂⊗𝒂=𝒂∙𝒂𝑇. 28 4.2 Anisot ropic Structural Tensors , Stress , and Invariant Once again let us first consider the case of the global coordinate system . The user must compute the structural tensor in the global coordinate system, which follows from Equation (57) as: 𝒂𝑮𝑖⊗𝒂𝑮𝑖=(𝑭∙𝑨𝑮𝑖)∙(𝑭∙𝑨𝑮𝑖)𝑇=𝑭∙𝑨𝑮𝑖∙𝑨𝑮𝑖𝑇∙𝑭𝑇 (59) In standard approaches to modelling fibre reinforced soft tissue s (e.g. HGO, Weiss, MA) the stress is given as scalar function of the deformation ( let us name it 𝜓) multiplied by the st ructural tensor. For a UMAT i n a global scheme , the stress is required w.r.t. the global basis system (denoted here as 𝝈𝑮𝑖) such that 𝝈𝑮𝑖=𝜓(𝑭) {𝒂𝑮𝑖⊗𝒂𝑮𝑖}=𝜓(𝑭) {𝑭∙𝑨𝑮𝑖∙𝑨𝑮𝑖𝑇∙𝑭𝑇} (60) In a classical corotational local formulation in which the local basis system 𝒆𝑖 rotates with the material, such that the global defor mation gradient is given as 𝑭=𝑹∙𝑼 and the local deformation gradient is simply 𝑭′=𝑼. In the case of a UMAT in conjunction structural membrane and shell elements in Abaqus/Standard, and for all ele ment types in a VUMAT in Abaqus/Explicit, the user is provided with 𝑭′=𝑼, and the structural tensor directly follows from equation (58) as follows: 𝒂𝒆𝒊⊗𝒂𝒆𝒊=𝑭′∙𝑨𝑬𝒊∙𝑨𝑬𝒊𝑇∙𝑭′𝑇=𝑼∙𝑨𝑬𝒊∙𝑨𝑬𝒊𝑇∙𝑼 (61) In a local scheme the stress is required w.r.t. the deformed local basis system (denoted here as 𝝈𝒆𝑖), and is simply given as 𝝈𝒆𝑖=𝜓(𝑭′) {𝒂𝒆𝑖⊗𝒂𝒆𝑖}=𝜓(𝑼) {𝑼∙𝑨𝑬𝑖∙𝑨𝑬𝑖𝑇∙𝑼𝑇} (62) Finally, i n Abaqus /Standard for the case of a continuum element us ed with a local basis system (through a *orientation statement), a UMAT provides 𝑭𝑎𝑙 as input and 𝝈𝒆𝑖 must be computed as output. Recalling that 𝑭𝑎𝑙=𝑹𝑇∙𝑭∙𝑹=𝑼∙𝑹, such that the deformed fibre vector is given as 𝒂𝒆𝒊=𝑭𝑎𝑙∙𝑹𝑇∙𝑨𝑬𝒊, the structural tensor and stress tensor in the deformed local basis system are given as 𝒂𝒆𝒊⊗𝒂𝒆𝒊=𝑼∙𝑨𝑬𝒊∙𝑨𝑬𝒊𝑇∙𝑼𝑇=(𝑭𝑎𝑙∙𝑹𝑇)∙𝑨𝑬𝒊∙𝑨𝑬𝒊𝑇∙(𝑹∙𝑭𝑎𝑙𝑇) (63) and 𝝈𝒆𝑖=𝜓(𝑭𝑎𝑙) {𝒂𝒆𝑖⊗𝒂𝒆𝑖}=𝜓(𝑭𝑎𝑙){𝑭𝑎𝑙∙𝑹𝑇∙𝑨𝑬𝒊∙𝑨𝑬𝒊𝑇∙𝑹∙𝑭𝑎𝑙𝑇} (64) Finally, we examine the scalar anisotropic invariant associated with the fibre vector and resulting structural tensor. In the literature it is commonly referred to as 𝐼4 and here we shall call it 𝐼𝑓, where 29 𝑓 is the fibre vector it is associated with8. 𝐼𝑓 is defined as the square of th e stretc h of the fibre vector and is be calculated in a global basis as, 𝐼𝑓=𝑨𝑮𝑖𝑇∙𝑪∙𝑨𝑮𝑖=𝑨𝑮𝑖𝑇∙𝑭𝑇∙𝑭∙𝑨𝑮𝑖 (65) For shell and membrane elements in Abaqus/Standard, and for all elements in Abaqus/Explicit the fibre invariant is computed from 𝐼𝑓=𝑨𝑬𝑖𝑇∙𝑼2∙𝑨𝑬𝑖 (66) And finally, fo r the case of continuum elements with a local basis system in Abaqus/Standard , 𝐼𝑓=𝑨𝑬𝑖𝑇∙𝑪𝑎𝑙∙𝑨𝑬𝑖 =𝑨𝑬𝑖𝑇∙𝑭𝑎𝑙 𝑇∙𝑭𝑎𝑙∙𝑨𝑬𝑖 (67) Clearly the same operations are impleme nted to compute 𝐼𝑓, regardless of whether a UMAT/VUMAT provides 𝑭𝑬𝑬, 𝑭𝒆𝑬 or 𝑭𝑎𝑙 as input. The invariant may be calculated alternatively by taking the trace of the structural tensor, 𝐼𝑓=tr(𝒂𝒆𝑖⊗𝒂𝒆𝑖) (68) where in the case of the global basis 𝒂𝒆𝑖=𝒂𝑮𝑖. As one might expect 𝐼𝑓 should remain constant , regardless of the choice of basis system be it global or local. This makes 𝐼𝑓 a useful quantity to check when verifying a UM AT as it should not change if the basis system is changed, or if a body is subject ed to a simple rigid ro tation. 4.3 Specific form of anisotropic hyperelastic constitutive law used in worked examples below: The expression for stress in fibre -reinforced aniso tropic hyperelastic constitutive models can be additively decomposed into isotropi c and anisotropic parts , 𝝈=𝝈iso+𝝈aniso. We will use the modified anisotropic (MA) model proposed by Nolan et. al (Nolan et al., 2014) to demonstrate local and global implementations. The isotropic part of the Cauchy stress is represented by the neo - Hookean co nstitutive model given in (47) in either the local or global basis, i.e. 𝝈iso=𝝈NH. The anisotropic part of the stress is given by the second term in (69) below. The summation term indicates that th ere are two families of fibres ( 𝑓=4,6) present in this constitutive model each with their own unique fibre vector. T he expression for the total Cauchy stress in the MA model is, 8 It is possible to have more than one fibre vector in this class of constitutive model. 30 𝝈MA=𝝈NH+2𝐽−1𝑘1∑ (𝐼𝑓−1)exp [𝑘2(𝐼𝑓−1)2](𝒂𝒆𝑖𝑓⊗𝒂𝒆𝑖𝑓) 𝑓=4,6, (69) where 𝑘1 and 𝑘2 are material parameters, 𝒂𝒆𝑖𝑓 is the fibre vector in the current configuration for fibre family 𝑓; this fibre vector may be in either a local or g lobal basis (in the global case 𝒆𝑖=𝑮𝑖). 31 4.4 Worked Example #6: Fibre kinematics in a continuum element with a local basis system We will again consider the same deformation motion 𝜒 as the previous Worked Examples and calculate the fibre vector and structural tensor in the current configuration for the case s where global and local bases are used. In this example we will impose the condition that the local basis vectors in the reference configuration are aligned with the global basis vectors, 𝑬𝑖=𝑮𝑖. Starting with the scenario of a global coordinate basis, we choose the fibre vector in the reference configuration to be a unit vect or in the 𝑮1−𝑮2 plane, at 30∘ to the 𝑮1 axis, 𝑨𝑮𝑖=[cos(30∘) sin(30∘) 0], (70) The fibre vector in the current configuration 𝒂𝑮𝑖 is then calculated using (41) and (70) in (57). So too is the structural tensor, 𝒂𝑮𝑖⊗𝒂𝑮𝑖, which is calculated using (59). The fibre invariant is calculated using (68), 𝒂𝑮𝑖=[1.003 0.493 0.173], 𝒂𝑮𝑖⊗𝒂𝑮𝑖=[1.005 0.495 0.174 0.495 0.243 0.085 0.174 0.085 0.030], 𝐼𝑓=1.2786 . (71) Next, we examine the local basis case and use the reference configuration fibre vector in (70) as the starting point, setting 𝑨𝑬𝑖=𝑨𝑮𝑖. The fibre vector in the c urrent configuration 𝒂𝒆𝑖 is then calculated using (42), (43) and (70) in (58). So too is the structural tensor, 𝒂𝒆𝑖⊗𝒂𝒆𝑖, which is calculated using (63). The fibre invariant is calculated using (68). 𝒂𝒆𝑖=[1.000 0.514 0.121], 𝒂𝒆𝑖⊗𝒂𝒆𝑖=[0.9995 0.5141 0.174 0.5141 0.2644 0.0623 0.174 0.0623 0.0147], 𝐼𝑓=1.2786 . (72) • The Matlab script Worked_example_ 5.m demonstrates the calculation of the global and local fibre vectors, structural tensors, and invariants in the scenario outlined above. 4.5 Worked Example # 7: Fibre kinematics in a continuum element with a local basis system where 𝑬𝑖≠𝑮𝑖 Next , we examine th e scenario where the local basis system in the reference configuration 𝑬𝑖 is rotated by 20∘ about the 𝑮3 axis, as per Worked Example #3. To remain consistent with Worked Example #6, we establish the fibre in the 𝑬1−𝑬2 plane and to be at 10∘ w.r.t. the 𝑬1 axis. Therefore, in a global sense the fibre is at 30∘ w.r.t. the 𝑮1 axis; the same as in Worked Example #6. The fibre vector in the reference configuration is, 32 𝑨𝑬𝑖=[cos(10∘) sin(10∘) 0]. (73) Once again, we start by using the change of basis in (45) to map the global deformation gradient in 𝑮𝑖 to the global deformation gradient in 𝑬𝑖. The local deformation gradient 𝑭20𝑎𝑙 is then computed using (46). Using 𝑭20𝑎𝑙 and the rotation matrix 𝑹20 (used to map 𝑭20 to 𝑭20𝑎𝑙) in (58) the fibre vector in the current configuration is calculated. The structural tensor is calculated using (63) and the fibre invariant using (68). 𝒂𝒆𝑖=[1.115 0.141 0.121], 𝒂𝒆𝑖⊗𝒂𝒆𝑖=[1.2440 0.1576 0.1351 0.1576 0.0200 0.0171 0.1351 0.0171 0.0147], 𝐼𝑓=1.2786 . (74) Note that the fibre invariant remains constant across (71), (72), and (74). This is a result that we should expect as there is the same deformation motion and fibre orientation in all of these scenarios. • The Matlab script Worked_example_ 7.m demonstrates the calcula tion of the local fibre vectors, s tructural tensors, and invariants in the scenario outlined above. 33 4.6 Worked Example # 8: Calculation of stress in a continuum element consisting of an anisotropic material We will now work through an example of how the Ca uchy stress is calculated using th e MA constitutive model for the deformation motion we have used for all of the Worked Examples previously. We assume that the material parameters 𝑘1=0.05 MPa and 𝑘2=2.0. Starting with the global case, we assume there are two families of fibres whos e vectors are ±20∘ to the 𝑮1axis, in the 𝑮1−𝑮2 plane. The isotropic stress is the same as the neo -Hookean stress in (52), the structural tensor in the current configuration is calculated using (59), and the fibre invariant is calculated using (68). Using these quantities in (69) results in the Cauchy stress , 𝝈MA=[0.2063 0.0456 0.0679 0.0456 0.0824 0.0576 0.0679 0.0576 0.2762] (75) Next, we consider the case of a local coordinate basis where the reference basis 𝑬𝑖=𝑮𝑖. We can use the local deformation gradient given in (43), and the isotropic stress given in (53). The structural tensor in the current configuration is comput ed using (63) and the fibre invariant using (68). Using these quantities in (69) gives the Cauchy str ess, 𝝈MA𝑎𝑙=[0.2168 0.0488 0.0755 0.0488 0.0779 0.0409 0.0755 0.0409 0.2708] (76) We can next compute the co -rotational stress using the right stretch tensor 𝑼 determined from the polar decomposition of the global defo rmation gradient . When 𝑼 is used in place of the global deformation gradient in the constitutive equation in (69) we compute the local stress in (76). Finally, t he stress in the local basis can be mapped back to the global stress using the change of basis operation , 𝝈MA=𝑹∙𝝈MA𝑎𝑙∙𝑹𝑇 (77) • The Matlab script Worked_example_ 8.m demonstrates the calculation of the global, local, and co -rotational stress tensors outl ined above. • The global stress tensor is calculated by the Abaqus/Standard input file WE8_global.inp together with the UMAT file umat_MA_global.for . • The local stress tensor is calculated by the Abaqus/Standard input file WE8_local.inp together with the UMAT file umat_MA_local.for . 34 4.7 Worked Example # 9: Calculation of stress in a continuum element consisting of an anisotropic material Figure 8 A) Deformation of a unit cube according to the deformation gradient in (41). B) Rigid body rotation of the deformation in A) given by (78) and (79). Note that the basis vectors in the current configuration 𝒆𝑖 have also rigidly rotated. C) The stress component 𝜎22 when the local UMAT is used, for the deformation in A). D) The stress componen t 𝜎22 when the local UMAT is used, after the rigid rotation deformation in B). Note that the stress remains constant, as we expect for a rigid body rotation in a co -rotati onal ba sis. E) The stress component 𝜎22 when the global UMAT is used, for the deformation in A). F) The resulting stress component 𝜎22 when the global UMAT is used, after the rigid rotation deformation in B). Note that the stress changes between E) and F) , and that this is incorrect. The happens because the g lobal UMAT fails to take account of rotation of the basis system 𝒆𝑖. Finally, we will demonstrate that the principal stresses, deformation and fibre invariants remain constant under a rigid body r otation deformation. In the specific case where a local co ordinate system is used, we will show that the stress tensor computed before and after rotation is identical. We start with the global and local stresses calculated in Worked Example #8 given in (75) and (76) respectively. We shall name 𝑭 given in (41) to be State 1. This body can be rigidly rotated using a rotation matrix 𝚽 to a State 2 using, 𝑭2=𝚽∙𝑭1 (78) where we use the subscript to indicate which state the deformation gradient exists. In this example we rotate the deformed unit cube −30∘ about the 𝑮3 axis. 𝚽=[cos(−30∘)−sin(−30∘)0 sin(−30∘) cos(−30∘) 0 0 0 1] (79) When 𝑭2 is used to calculate stress i n the MA constitutive model using the same parameters as the global case in Worked Example #8 we calculate the tensor, 35 𝝈MA=[0.2148 −0.0309 0.0876 −0.030 9 0.0739 0.0159 0.0876 0.0159 0.2762]. (80) When we calculate the pri ncipal stresses and pressure stress of global stress tensors in (75) and (80) we find that they are identi cal. This is an important result and demonstrates that there is the same stress state before and a fter the rigid body rotation. Though the individual components of the stress tensor are different in (75) and (80) this is only because rigid body rotation has changed the orientatio n of the stress compon ents w.r.t. the global coordinate basis. To calculate the stress when a local coordinate frame is used we first perform a polar decomposition of 𝑭2 from w hich we determine the rotation matrix 𝑹2. The local deformation gradient in State 2 may then b e calculated using a change of basis, 𝑭2𝑎𝑙=𝑹2𝑇∙𝑭2∙𝑹𝟐 (81) Following the same procedure for the calculation of the local stress as Worked Example #8 we calculate the local stress, 𝝈MA𝑎𝑙=[0.2162 0.0488 0.0792 0.0488 0.0779 0.0409 0.0755 0.0409 0.2708] (82) This stress tensor is identical to (76) demo nstrating that when a local coordinate basis is used, the stress tensor co -rotates with a rigid bo dy rotation. The co -rotational stress in State 2 is calculated by using the right stretch tensor 𝑼2, determined from a polar decomposition of 𝑭2, in p lace of the deformation gradient in (69). This leads to the calculation of the co -rotational stress tensor and is identical to (82). Figure 8 shows the stress component 𝜎22 computed when a local material orientation is implemented and the local UMAT is used ( Figure 8C and D). As you can see the stress does not change with rigid body rotation and 𝜎22 is as per (76) and (82). This is the expected, correct result. On the other hand, when the global UMAT is used, as the body rotates from Figure 8E to F the stress changes. This is incorrect and a clear indication that the UMAT is not frame invariant. • The Matlab script Worked_example_ 9.m demonstrates the calculation of the global, local , and co -rotational stress tensors outlined above. • The global stress tensor is calcul ated by the Abaqus/Standard input file WE9_global.inp togeth er with the UMAT file umat_MA_global.for . • The local stress tensor is calculated by the Abaqus/Standard input file WE9_local.inp together with the UMAT file umat_MA_ local.for . 36 Figure 9 Torque of an anisotropic hyperelastic tube. Anisotropic invariant 𝐼4 when A) the correct, local UMAT is used, and B) when a global UMAT is used incorrectly. Likewise, the circumferential stress 𝜎𝜃𝜃 C) and D), and shear stress 𝜎𝜃𝑍 E) and F), are calculated using the local and global UMATs respectively. 5 Case studies The two case studies presented in this example sho w the consequences of incorrectly implementing a local coordinate system for a fibre reinforced hyperelastic cons titutive model. These examples have been run using Abaqus/St andard 2017 and the global and local UMATs referred to in Section 4.6 and 4.7. The simulation files required to run these examples are inc luded with the accompanying codes in the folder Case Studies . 5.1 Torque and axial stretch of a tube Consider a tube whose inner radius is 2 mm, outer radius is 3 mm, and whose length is 10 mm . One end of the tube is fixed for all degrees of freedom. A to rque is applied to t he tube by rotating the opposite end by 1 rad ian, and an axial stretch by displacing this end by 0.5 mm. The fibre reinforced constitutive model given in (69) is impl emented the rep resent the material of the tube , with parameters 𝐶10=0.2 MPa , 𝐷=2.0 MPa, 𝑘1=0.2 MPa, 𝑘2=2.0. Two families of fibres are used at ±30∘ to the circumferential direction , using a local cylindrical coordinate system to define these directions . This problem is simulated using FE analysis using b oth the global and local UMATs presented in Section 4.6. Results for this simulation are given in Figure 9 when both the local and global UMATs. Figure 9A) and B ) show the fibre i nvariant 𝐼4 when the local and global UM ATs are used respectively. In a similar fashion we show the circumferential stress 𝜎𝜃𝜃 C) and D), and the shear stress 𝜎𝜃𝑍 E) and F). 37 In the cases where there is the local UMAT is used, we see that there is a s ymmetric pattern in each of the co ntours as expected. However, when the global UMAT is used with a local coordinate system, we see that this not the case. The contou rs on the top half of the tube are different to those on the bottom. This is because there has been large rotation that the g lobal UMAT does not account for, as we go axially from bottom to top, the rotation becomes larger, and the error is higher. 5.2 Closin g an artery segment to create residual stress It is widely understood that there is a resid ual stress in arteries in their ze ro-pressure configuration. One method to measure this experimentally is to cut arteries and measure their opening angle (Chuong and Fung, 1986) . Conversely, to recreate the residual stress in a FE simulation the opened geometry is restored to the closed geometr y, see Figure 10A). To define fibre directions, it is useful to d efine a local, concentric basis sy stem based on polar coordinates. Two families of fibres are defined at ±45∘ to the circumferential direction, and the constitutive model given in Section 5.1 is implemented again here. When the artery is closed the stress an d strain in the artery wall should also be concentric. Figure 10B) shows the circumferential stress 𝜎𝜃𝜃 when the local UMAT is used to compute the stress; as expected we calculate a concentric circumferential stress contour. On the other hand, Figure 10C) shows the circumfe rential stress 𝜎𝜃𝜃 when the global UMAT is used to compute the stress. In this case an eccentric stress contour is calculated. This result is incorrect, and it is useful to n ote that the more rotation an element undergoes, the larger the error in the stress. This is because the global UMAT fails to correctly account for the deformation of the fibre w.r.t. the local coordinate system. This example again highlights the importanc e of including rotations when a local coordinate basis is used in conjunctio n with a fibre vector based UMAT. 38 Figure 10 A) Schematic of the closing of an artery segment to create a residual stress in the artery wall. B) Contour plot showing concentric isolines of the resulting circumferential stress 𝜎𝜃𝜃 when the local UMAT is used. C) Contour plot showing isolines of circu mferential stress 𝜎𝜃𝜃 when the global UMAT is used. Note that these are eccentric w.r.t. the geom etry. 39 6 Discussion and Conclusions In this tutorial paper we have discusse d the form of the deformation gradient when either a global or co -rotational local basis system is used to describe the deformation. Key points arising from this paper are summarised as follows : 1. There are important differences in the way in which the local deformation gradient is calculated depending on whether a structural element (e.g . a shell or membrane element) or a continuum element is used in Abaqus. (i) A structural element im plements what we term the ‘Classical’ definition of the local deformation gradient, i.e. where 𝑭′=𝑼. (ii) A continuum element by default uses a global coordinate system with a global deformation gradient 𝑭=𝑹∙𝑼. (iii) If a local orientation is speci fied for a continuum element a different form of the deformation gradient is provided. We refer to this f orm as the Abaqus local deformation gradient , which is given as 𝑭𝑎𝑙=𝑹𝑇∙𝑭∙𝑹. The consequences of the construction of the Abaqus local defo rmation gradient are discussed in Sections 2.2 and 2.3. 2. This important difference has consequences for the interpretation of the local deformation tensors , such as the Cauchy -Green tensors, depending on whether a continuum or structural element is used. In particular, one s hould be aware that the eigenvectors o f these tensors may not those expected (see Appendix A 1). 3. When writing a user -defined material subroutine (UMAT/VUMAT) to compute the stress from a given deformation gradient, one must ensure that the Cauchy stress is returned in the current configuration basis system, be that global or local. A summary is provided in Table 1. 4. The construction of the Abaqus local deformat ion gradient for continuum element s with a local basis system has importan t consequences for the calculation of fibre vectors in the current configuration for anisotropic hyperelastic models . The fibre vector in the reference configuration must have its ba sis changed from 𝑬𝑖 to 𝒆𝑖 leading to Equation (58) where 𝒂𝒆𝑖=𝑼∙𝑨𝑬𝑖=𝑭𝑎𝑙∙𝑹𝑇∙𝑨𝑬𝑖. 5. When very large deformations are accompanied by very large rotations , a slight numerical error in the Abaqus local deforma tion gradient is observed (Section 2.4.1). 𝑭𝑎𝑙=𝑹𝑇∙𝑭∙𝑹. However , no such numerical error is observed in the log strain tensor. Such slight numerical errors are not observed for cases of large deformation with small rotations, or for cases of small deformations with large rotations. An incomplete understanding of the constriction of the global ( 𝑭), classical local ( 𝑭′), and Abaqus local (𝑭𝑎𝑙) deformation gradient when writing a UMAT/VUMAT may lead to insidious error s that will not be detectable through uniaxial or biaxial verification test . Test cases that apply material 40 deformation in conjunction with rigid body rotation must be performed to ensure that the deformation gradient is correctly interpreted. It is ess ential for the progress of the field of computational mechanics that user -defined materials are rigorously verified and validated before their emergent results are published. The authors hope that this paper will provide the reader with a good guide to cor rectly implementin g user defined material subroutines with local orientation and finite deformations. 41 Continuum Elements Structural Elements Abaqus/Standard Abaqus/Explicit# *orientation used? No Yes No Yes No Yes Deformation tensor (input) 𝑭 𝑭𝑎𝑙 𝑼𝑮𝑖 𝑼𝑬𝑖 𝑭′ 𝑭′ Fibre vector in current configuration 𝒂𝑮𝑖=𝑭∙𝑨𝑮𝑖 𝒂𝒆𝑖=𝑭𝑎𝑙∙𝑹T∙𝑨𝑬𝑖 𝒂cr=𝑼𝑮𝑖𝑨𝑮𝑖 𝒂cr=𝑼𝑬𝑖𝑨𝑬𝑖 𝒂𝒆𝑖=𝑭′𝑨𝑬𝑖 𝒂𝒆𝑖=𝑭′𝑨𝑬𝑖 Cauchy stress tensor (output) 𝝈𝑮𝒊 𝝈𝒆𝒊 𝝈cr* 𝝈cr 𝝈𝒆𝒊 𝝈𝒆𝒊 Table 1 Guide to the deformation inputs and expected stress o utputs for a given coordinate basis when a user -defined material is used in either Abaqus/Standard (UMAT) or Abaqus/Explicit (VUMAT). #Both the deformation gradient and the stretch tensor are input s in a VUMAT ; the Abaqus manual recommends the use of the stretch tensor as this will result in the correct corotational stress being computed. *The Fortran subroutine must return the co -rotate d stress, however Abaqus will internally rotate the stress to the global basis 42 Appendices A.1 Relating 𝑭𝐸 to 𝑭′ 𝑭=∑𝜕𝑥𝑖𝒆 𝜕𝑋𝑗𝑬(𝒆𝑖⊗𝑬𝑗)3 𝑖,𝑗=1=∑𝜕𝑥𝑖𝒆 𝜕𝑋𝑗𝑬𝜃𝑘𝑖𝑬𝑘⊗𝑬𝑗=∑𝜕𝑥𝑖𝑬 𝜕𝑋𝑗𝑬(𝑬𝑖⊗𝑬𝑗)3 𝑖,𝑗=13 𝑖,𝑗,𝑘=1 (A1.1) where the components of the local deformation gradient in the 𝒆𝑖⊗𝑬𝑗 basis system are 𝐹𝑖𝑗′=𝜕𝑥𝑖𝒆 𝜕𝑋𝑗𝑬 and the components of the global deformation gradient in the 𝑬𝑖⊗𝑬𝑗 basis system are 𝐹𝑖𝑗𝐸=𝜕𝑥𝑖𝑬 𝜕𝑋𝑗𝑬 Expanding the summat ion over i, j and k for the th ird expression in equation (A1.1) above gives , ∑ 𝐹𝑖𝑗′𝜃𝑘𝑖𝑬𝑘⊗𝑬𝑗3 𝑖,𝑗,𝑘=1= 𝐹11′(𝜃11𝑬1⊗𝑬1+𝜃21𝑬2⊗𝑬1+𝜃31𝑬3⊗𝑬1) +𝐹12′(𝜃11𝑬1⊗𝑬2+𝜃21𝑬2⊗𝑬2+𝜃31𝑬3⊗𝑬2) +𝐹13′(𝜃11𝑬1⊗𝑬3+𝜃21𝑬2⊗𝑬3+𝜃31𝑬3⊗𝑬3) +𝐹21′(𝜃12𝑬1⊗𝑬1+𝜃22𝑬2⊗𝑬1+𝜃32𝑬3⊗𝑬1) +𝐹22′(𝜃12𝑬1⊗𝑬2+𝜃22𝑬2⊗𝑬2+𝜃32𝑬3⊗𝑬2) +𝐹23′(𝜃12𝑬1⊗𝑬3+𝜃22𝑬2⊗𝑬3+𝜃32𝑬3⊗𝑬3) +𝐹31′(𝜃13𝑬1⊗𝑬1+𝜃23𝑬2⊗𝑬1+𝜃33𝑬3⊗𝑬1) +𝐹32′(𝜃13𝑬1⊗𝑬2+𝜃23𝑬2⊗𝑬2+𝜃33𝑬3⊗𝑬2) +𝐹33′(𝜃13𝑬1⊗𝑬3+𝜃23𝑬2⊗𝑬3+𝜃33𝑬3⊗𝑬3) (A1.2) Expanding the summation over i and j for the fourth expression in equation (A1.1) above gives , ∑ 𝐹𝑖𝑗𝐸(𝑬𝑖⊗𝑬𝑗)3 𝑖,𝑗=1= 𝐹11𝐸𝑬1⊗𝑬1+𝐹21𝐸𝑬2⊗𝑬1+𝐹31𝐸𝑬3⊗𝑬1 +𝐹12𝐸𝑬1⊗𝑬2+𝐹22𝐸𝑬2⊗𝑬2+𝐹32𝐸𝑬3⊗𝑬2 +𝐹13𝐸𝑬1⊗𝑬3+𝐹23𝐸𝑬2⊗𝑬3+𝐹33𝐸𝑬3⊗𝑬3 (A1.3) Equating the terms in (A1.3) and (A1. 2) which have the same dyadic product in common, and writing in matrix form we get , 43 [𝐹11𝐸𝐹12𝐸𝐹13𝐸 𝐹21𝐸𝐹22𝐸𝐹23𝐸 𝐹31𝐸𝐹32𝐸𝐹33𝐸]= [𝜃11𝐹11′+𝜃12𝐹21′+𝜃13𝐹31′𝜃11𝐹12′+𝜃12𝐹22′+𝜃13𝐹32′𝜃11𝐹13′+𝜃12𝐹23′+𝜃13𝐹33′ 𝜃21𝐹11′+𝜃22𝐹21′+𝜃23𝐹31′𝜃21𝐹12′+𝜃22𝐹22′+𝜃23𝐹32′𝜃21𝐹13′+𝜃22𝐹23′+𝜃23𝐹33′ 𝜃31𝐹11′+𝜃32𝐹21′+𝜃33𝐹31′𝜃31𝐹12′+𝜃32𝐹22′+𝜃33𝐹32′𝜃31𝐹13′+𝜃32𝐹23′+𝜃33𝐹33′] (A1.4) Clearly, we can write , 𝑭𝐸=𝜽∙𝑭′ (A1.5) A.2 Consequences for use of utility subroutine SPRIND The difference between the Classical and Abaqus definition of the local deform ation gradient has important implications for the use of the inbuilt utility subroutine sprind in a UMAT. This subroutine is use d to determine the principal directions (eigenvectors) of stress, strain, and deformation tensors. If a global basis system is u sed, 𝑪 does not co ntain any information on rigid body rotation of the material (see (4)) and therefore the principal directions are simply the reference basis vectors, i.e. 𝑰. On the other hand, 𝑩 contains information on the materia l rigid body r otation (see (6)) and therefore the use of 𝑩 in the sprind statement returns the principal strain directions w.r.t. the global basis system. This outcome is expected from the Classical continuum mechanics framework. In contrast, when *orientation is used, the construction of 𝑪𝑎𝑙 from the local deformation gradient given in (27), with its reversed order of 𝑹 and 𝑼, results in a quantity that now does contain information on material ro tation, and th e use of 𝑪𝑎𝑙 in sprind computes the principal directions w.r.t. the global coordinate system. The opposite is now the case for the 𝑩𝑎𝑙 matrix which does not contain rotational information (see (25)) and c annot now be u sed to determine principal material orientations w.r.t. the global directions. A.3 Table of Accompanying Codes Filen ame Description Worked_example_1 Calculation of 2D deformation gradients looking in particular at the difference between structu ral elements, which calculate the local deformation gradient using the Classical definition, and continuum elements, which use a different transformation as outlined in Section 2.3 Worked_example_2 Looking at the difference in the deformation gradient be tween using a global and local coordinate system. 44 Worked_ex ample_3 Same as above, except that the local coordinate system is offset by 20 degrees about the 𝑮3-axis in the reference configuration. i.e. 𝑬𝑖≠𝑮𝑖. Worked_example_4 Calculating the neo-Hookean s tress for the global and local coordinate systems, as per Worked_example_1 Worked_example_5 Calculating the neo -Hookean stress for the global and local coordinate systems, as per Worked_example_2 Worked_example_6 Calculatin g the fibre ori entation vecto rs and structural tensor when a local or global coordinate scheme is used Worked_example_7 Same as above, except that the local coordinate system is offset by 20 degrees from the x -axis in the reference configuration. i.e. 𝑬𝑖≠𝑮𝑖. Worked_example _8 Calculate the stress for the MA hyperelastic model in both the local and global coordinate schemes. abaqus j=WE 8_global user=umat_MA_global.for abaqus j=WE 8_local user=umat_MA_ local.for Worked_example_9 Same as 7, except perform second step which is a rigid body rotation. Show that when a local scheme is used that there is no change in the stress upon this rigid rotation. abaqus j=WE 9_global user=umat_MA_global.for abaqus j=WE 9_local user=umat_MA_ local.for Case Study 1 Torque and stretch of a tube. abaqus j=localTorque inp=torque_stretch_cylinder.inp user=umat_MA_local.for abaqus j= globalTorque inp=torque_stretch_cylinder.inp user=umat_MA_ global.for Case Study 2 Creating a residual stress in an artery segment. abaqus j= localResidual inp=residualStress .inp user=umat_MA_local.for abaqus j= globalResidual inp=residualStress .inp user=umat_MA_ global.for A.4 Outputting the deformation gradient Abaqus does not write out the deformation gradient by default, and only outputs it when hyperelastic and user -defined ma terials are used. To output the deformation gradient the keyword options below must be added to the Abaqus .inp file. For example in Wor ked Example #2, file WE1_local.inp , line 118 has: *EL PRINT DG 45 This outputs the deformati on gradient in each element, at each integration point into the resulting .dat file for the simulation. The Fortran code below will print out the matrix tens of size m × n to the .da t file when called from a (V)UMAT. subroutine kprinter(tens, m, n) INCLUDE 'ABA_PARAM. INC' intent(in):: tens, m, n dimension tens(m,n) write(6,*) do i = 1,m do j = 1,n write(6,'(e19.9)' ,advance= 'no'),tens(i,j) end do write(6,*) end do write(6,*) return end subroutine kprinter 46 Acknowledgements This research was funded by Irish Research Council Enterprise Partnership Postdoctoral Fellowship EPSPD/2017/111 and Science Foundation Ireland unde r grant number s: SFI 12/IP/1723 , SFI 18/ERCD/5481 and SFI 13/CDA/2145. The authors wish to thank Prof. Michel Destrade for useful discussion on the Classical definition of the local deformation gradient, and Prof. Nele Famaey and Dr. Eoin McEvoy for discus sions on Abaqus UMATs /VUMATs . 47