Automatic differentiation for solid mechanics
AAutomatic differentiation for solid mechanics
Andrea Vigliotti ∗ and Ferdinando Auricchio Innovative Materials Laboratory, Italian Aerospace Research Center, 81043 Capua, Italy Department of Civil Engineering and Architecture, University of Pavia, 27100 Pavia, Italy
January 22, 2020
Abstract
Automatic differentiation (AD) is an ensemble of techniques that allow to evaluate accuratenumerical derivatives of a mathematical function expressed in a computer programming language.In this paper we use AD for stating and solving solid mechanics problems. Given a finite elementdiscretization of the domain, we evaluate the free energy of the solid as the integral of its strainenergy density, and we make use of AD for directly obtaining the residual force vector and the tan-gent stiffness matrix of the problem, as the gradient and the Hessian of the free energy respectively.The result is a remarkable simplification in the statement and the solution of complex problemsinvolving non trivial constraints systems and both geometrical and material non linearities. To-gether with the continuum mechanics theoretical basis, and with a description of the specific ADtechnique adopted, the paper illustrates the solution of a number of solid mechanics problems,with the aim of presenting a convenient numerical implementation approach, made easily availableby recent programming languages, to the solid mechanics community.
The general problem of solid mechanics hinges on the following familiar pointwise Cauchy’s equilibriumequations [Asaro and Lubarda, 2006] σ ij,j + b i = 0 in V (1a) σ ij n j − t i = 0 on S, (1b)which describes the equilibrium of a body subjected to a system of external forces, where i ∈ { , , } are the coordinate directions of a Cartesian frame of reference, V is the current volume of the body, S ≡ ∂V is the boundary of V , σ ij are the components of the Cauchy stress tensor, with the lettersafter the comma in the subscript denoting derivatives along spatial directions and the repeated in-dex denoting summation; furthermore t i denotes the components of the tractions on S , and b i thecomponents of the body forces in V . Together with stress-strain relationships, and with the kine-matic compatibility relations, equations (1) define a boundary value problem that allows to find thedeformed configuration of a body, given its initial configuration, its material’s constitutive laws and aset of boundary conditions.For the simple cases it is possible to find analytical solutions for the above equilibrium equations [Tim-oshenko and Goodier, 1987]. Nevertheless, in the general practice, solutions are usually sought throughnumerical methods, such as the Finite Element (FE) method [Zienkiewicz et al., 2005, Bathe, 2014].FE methods are based on a twofold discretization of the problem: (i) the domain and its boundary arediscretized in elements connected at nodes; and (ii), the solution is approximated as the weighted sum ∗ [email protected], [email protected] a r X i v : . [ m a t h . NA ] J a n f a finite set of shape functions, associated with the elements. The weighting coefficients that controlthe solution define the Degrees of Freedom (DoFs) of the problem. In this framework, the problemcan be stated and solved, using weighted residual methods such as Galerkin [Zienkiewicz and Taylor,2005], that are based on selecting the values of the DoFs that minimize a given norm of the error overthe domain of interest. Thus, the problem reduces to finding the solutions of a system of non-linearequations in terms of the DoFs. The error on each equation defines the residual force vector, and aparticular solution can be found bringing the residual to zero, by means of iterative techniques, suchas Newton-Raphson, that make use of the tangent stiffness matrix, which coincides with the Jacobianof the residual vector, to update the trial solution on the basis of its residual.Despite the procedural and algorithmic nature of the general approach, the computer implementa-tion of FE methods for the solution of complex boundary value problems includes many challengingaspects. With particular reference to solid mechanics (but similar considerations apply to any contin-uum mechanics problem), the actual statement of equilibrium, (1), involves the modelling of complexmaterial behaviours and requires expressing several vector and tensor quantities, which are best de-fined in specific, and distinct, reference systems. The implementation of these steps can in principle beautomated, and a number of methodologies have been proposed to facilitate and, at various extents,automate the generation of computer programs capable of efficiently state and solve FE problems. Tothis end, the FEniCS project [Logg, 2007, Logg et al., 2012] set the goal of automating Computa-tional Mathematics Modelling (CMM) problems in general, including the FE methods. The projectaimed at the mechanization of the essential discretization steps of any CMM problem, by means of asuite of general purpose, high level, C++ and python libraries that allow to deal with the numericalimplementation of general physical models in a quite abstract, yet efficient, manner, and provide aninterface for the definition of the variational problem, its boundary conditions, and its solution.One essential undertaking in the FE implementation of complex solid mechanics problems is trans-lating the mathematical models of the physical processes into a form that can be incorporated into theFE formalism. This step generally involves analytical manipulations of various sorts, that are normallydone by hand and can be source of errors. In order to address this matter, Korelc and Wriggers [2016]proposed an integrated methodology based on the use of a symbolic engine for the automation of codegeneration starting from an abstract mathematical statement of the physics under consideration. Inthe approach proposed by Korelc and Wriggers the handling of the symbolic expressions is carriedout through the use of AceGen , a package within the Mathematica software suite, whose end-productis the source code implementing the models, in different programming languages, for the use in FEprograms. The entire environment includes different components that are capable to produce efficientsource code for the generation of the FE residual vector and stiffness matrix.The approach presented here makes use of
Automatic Differentiation (AD) for the numerical eval-uation of the FE residual force vector and stiffness matrix. In fact, if the material of the body is a
Green elastic type, for which the deformation work is an exact differential, it is possible to write theexpression for the free energy of the solid, and the gradient and the Hessian of the free energy, withrespect to the DoFs of the problem, take the meaning of the residual force vector and of the tangentstiffness matrix. According to the method described in the paper, the gradient and the Hessian ofthe free energy are not explicitly calculated, but are automatically obtained, through AD, from thefunction that evaluates the free energy. The resulting formulation is particularly streamlined and in-sightful, with the surprising consequence that, with the approach described here, it is possible to writea finite element program without introducing the concept of a stress tensor.AD is an ensemble of techniques that allows for the numerical evaluation of the derivatives of afunction with the same accuracy of the function itself. AD differs from finite differences, because itdoes not approximate the continuous derivatives with discrete differences, thus it does not suffer fromtruncation error, and the only source of error is the inevitable round-off error, due to the finite precision2epresentation of real numbers. The differentiation techniques based on AD rely on the assumptionthat the numerical evaluation of the value of a function, regardless of the complexity of the functionitself, is always decomposed into in a sequence of elementary sub-expressions by computers. Therefore,if the analytical derivatives of the sub-expressions are know, it is possible to evaluate the derivativesof the entire function, with respect to the given independent variables, by operating on the partial re-sults. In this way, AD allows the evaluation of the gradient, along with higher order derivatives, of anycomputable differentiable function, without having to explicitly write computer code for the derivatives.The paper is structured as follows. First, we recall the equivalency between the pointwise, strongform statement of equilibrium equations and free energy stationarity. Subsequently, we discuss theproperties of the dual numbers systems, which is the numerical method chosen to implement automaticdifferentiation in the present paper. Lastly we discuss a number of examples of solid mechanicsproblems that have been solved using AD. The approach presented in here has been implemented inthe Julia programming language, an open source programming language of recent introduction thatcombines high level of abstraction, high expressiveness and fast execution time [Bezanson et al., 2017,Perkel, 2019]. All the scripts developed for producing the examples presented in this paper are availableunder a collaborative licence [Vigliotti, 2019].
In conservative systems the deformation work is an exact differential, and it is possible to use the con-cept of free energy of the system for finding its equilibrium configurations. In fact, thermodynamicsguarantees that all, and only, the configurations that make the free energy stationary are equilibrated[Callen, 1985]. Therefore, it is possible to directly write the equilibrium statement by setting thegradient of the system’s free energy to nought. The advantage in using free energy minimization isthat free energy is always a scalar quantity, independent of the frame of reference, which is generallywell defined and can be calculated using a FE discretization. In the following we will illustrate theequivalence between free energy minima and equilibria as stated by equation (1), which will allow usto introduce all the relevant solid mechanics quantities as well as to expose the connections betweenthe stress tensor and the gradient of free energy density.Let B be a deformable body, occupying a region of an Euclidean space. Let’s assume that B is subjected to some external actions from the surrounding environment, in the form of body andsurface forces, and in the form of mechanical constraints that restrict its motion. In its interactionswith the environment, the body deforms and can take different configurations. Let us define the reference configuration of B as the configuration taken by the body when all of external actions areremoved. Given the reference configuration, the Cartesian coordinates of the points of B in anydeformed configuration are given by x i = X i + u i , (2)where X i are the coordinates of the points of B in the reference configuration and u i are the componentsof a displacement field mapping the position of the points of B from the reference configuration to the current configuration . In association with u i , it is possible to introduce the deformation tensor F , asa pointwise measure of the deformation in B , whose components are given as F ij = x i,j = δ ij + u i,j (3)where δ ij is the Kronecker symbol and the letter after the comma in the subscript denotes differentiationalong direction i . If we assume that B is made of a Green elastic material, it is possible to define astrain energy density function, φ , such that the total deformation energy can be expressed asΦ = (cid:90) V φ d V , (4)3here V is the volume occupied by B in the reference configuration and φ is a function of F ij in V ,with units of energy per reference unit volume. Let us also assume that B is subjected to a systemof conservative body and surface forces, whose potentials, per unit reference volume and per unitreference surface, are b and t respectively. In accordance with the above assumptions the total freeenergy of the body is given as Ψ = (cid:90) V ( φ − b ) d V − (cid:90) S t d S , (5)with S ≡ ∂V . Thermodynamics minimum free energy principle ensures that equilibrium configu-rations coincide with the configurations that make the free energy stationary. Therefore, u i is anequilibrium configuration if and only if δ Ψ = 0 ∀ δu i , (6)where δu i is an arbitrary variation in the space of the configurations compatible with the boundaryconditions. We now recall that, through a mere change of variables, the integral (5), and its variations,can be evaluated in any arbitrary configuration, provided that a mapping exists between the points ofthe reference configuration and the given configuration. Therefore the following holds δ Ψ = (cid:90) V (cid:20) ∂φ∂F ij δF ij − ∂b ∂u i δu i (cid:21) d V − (cid:90) S ∂t ∂u i δu i d S = (7a)= (cid:90) V (cid:20) J − ∂φ∂F ij δF ij − ∂b∂u i δu i (cid:21) d V − (cid:90) S ∂t∂u i δu i d S, (7b)with δF ij = δu i,j (8a) J = det ( F ) (8b) b = b J − (8c) t = t n i F ik F jk n j J − (8d)where V and S denote the current configuration, as in (1), and n i are the components of the localnormal on S . We observe that with the substitutions (8c) and (8d), b and t take the meaning ofthe potential of the external forces per unit current volume and unit current surface, respectively.We also observe that through mathematical manipulations, equation (7b) can be written in terms ofdisplacement variations taken with respect to the current configuration, δ ˜ u i , as follows δ Ψ = (cid:90) V (cid:20) J − ∂φ∂F ij F kj F − hk δF ih − ∂b∂u i F − ik F kj δu j (cid:21) d V − (cid:90) S ∂t∂u i F − ik F kj δu j d S == (cid:90) V [ σ ij δ(cid:15) ij − b i δ ˜ u i ] d V − (cid:90) S t i δ ˜ u i d S , (9)where F − ij are the components of the inverse of F , thus F kj F − hk = δ jh , with the following substitutions:4 ij = J − ∂φ∂F ih F jh (10a) δ ˜ u i = F ij δu j (10b) ∂ · ∂ ˜ u i = ∂ · ∂u k F − ki (10c) δ(cid:15) ij = δ ˜ u i,j = F − hj δF ih (10d) b i = ∂b∂ ˜ u i (10e) t i = ∂t∂ ˜ u i , (10f)After observing that the following identity holds( σ ij δ ˜ u i ) , j = σ ij,j δ ˜ u i + σ ij δ ˜ u i,j , (11)we can express the first variation of the free energy of B , making use of the divergence theorem, as δ Ψ = − (cid:90) V ( σ ij,j + b i ) δ ˜ u i d V + (cid:90) S ( σ ij n j − t i ) δ ˜ u i d S = 0 ∀ δ ˜ u i . (12)Standing the arbitrariness of δ ˜ u i , it follows that each of the integrands in equation (12) have to beseparately equal to nought everywhere in d V and on d S . Therefore, equation (12) is equivalent toequation (1). We also observe that the same procedure, starting from equation (7a), also leads to theequilibrium equation in terms of the nominal stress tensor, or first Piola-Kirchoff tensor, as follows P ij,j + b i = 0 in V (13a) P ij N j − t i = 0 on S , (13b)where N j are the component of the local normal to the surface on the undeformed configuration, with P ij = ∂φ∂F ik (14a) b i = ∂b ∂u i (14b) t i = ∂t ∂u i . (14c)We have thus obtained the equivalence between mechanical equilibrium, in the Newton’s laws sense,and the equilibrium in the thermodynamic sense, as the minima of free energy. In the following we discuss various means for numerically evaluate the derivatives of a function andillustrate the general aspects of automatic differentiation. In particular, we will discuss with greaterdetail the dual number system, which is the frame used to implement AD in the present study.5 .1 Discrete derivatives approximations
Finding accurate estimates of the derivatives of multivariate functions at low computational costs isessential in many fields of applied sciences and engineering. The simplest numerical estimate for aderivative is finite difference . Finite difference is based on the truncated Taylor expansion formula foranalytical functions, and produces the following approximation for a first order derivative ∂f∂x i = f ( x + ∆ x i ı i ) − f ( x )∆ x i + O ( | ∆ x i | ) , (15)where f is a scalar function of the vector x = x i i i while i i are the independent directions of thespace of x and ∆ x i is a finite increment for the i -th component. The estimates obtained throughequation (15) require one additional function evaluation per each independent variable, and sufferfrom a truncation error of order O ( | ∆ x i | ). We remark that the term O ( · ) does not represents anactual quantification of the approximation of the formula, but it rather represents the convergencerate to the exact value as ∆ x i collapses. In addition to truncation error, another important source ofinaccuracy is represented by the inevitable round-off, due to the finite precision of the floating pointrepresentation of real numbers. The effects of round-off error are particularly significant in expressionsof the type of equation (15), which involves small differences of finite quantities on the numerator,and on the ratio of two small numbers. As a consequence, it is not possible to arbitrarily enhancethe accuracy in the estimate of the derivatives by simply reducing the step length ∆ x i . One way toimprove the accuracy in the estimation of first derivatives for a fixed ∆ x i is through central differencescheme as follows ∂f∂x i = f ( x + ∆ x i ı i ) − f ( x − ∆ x i ı i )2 ∆ x i + O (cid:16) | ∆ x i | (cid:17) , (16)at the cost of two additional function evaluations per variable. In a similar way, other formulas canbe devised that offer improved estimates of the derivative at the cost of a larger number of functionevaluation. We also observe that although the above formulas can be applied recursively for the cal-culation of higher order derivatives, the accuracy of such estimates rapidly deteriorates since they arebased on already noisy estimates of lower order derivatives. Complex step is an alternative technique that allows to mitigate round-off error from finite difference[Lyness and Moler, 1967, Lyness, 1968]. The complex step method allows to estimate the derivativesof the analytical functions that can be evaluated on the complex plane as follows. Given the followingTaylor series expansion of the function along the imaginary axes of the i -th component, f ( x + ı ∆ x i ı i ) = f ( x ) + ∂f∂x i ı ∆ x i − ∂ f∂x i ∂x j ∆ x i ∆ x j
2! + O (cid:0) (cid:107) ∆ x (cid:107) (cid:1) , (17)where ı is the imaginary unit, an estimate of the first partial derivative with respect to x i can beobtained from the imaginary part of the above series as follows ∂f∂x i = Im { f ( x + ı ∆ x i ı i ) } ∆ x i + O (cid:16) | ∆ x i | (cid:17) . (18)As we can observe, equation (18) provides an approximation O (cid:16) | ∆ x i | (cid:17) of the first derivatives in asingle complex valued evaluation of f ( x ), which does not suffer from the round-off error due to thedifference on the numerator. Nonetheless, formulas based on (18) still need to be evaluated once pereach independent variable, at the higher cost of the evaluation on the complex field. As opposite to finite difference, symbolic differentiation treats mathematical expressions as strings ofsymbols, and applies the analytical differentiation and simplification rules to produce mathematical6ymbolic expressions that can be evaluated in any programming language. The availability of robustsymbolic differentiation engines has prompted the development of computational approaches tendingto the automation of FE code generation [Korelc and Wriggers, 2016]. However, the expressionsof the derivatives obtained by symbolic differentiation are often far from optimal with respect tocomputation and memory allocation cost, as they might include redundancy and repetitions thatsymbolic simplification steps are not capable of removing. In addition, symbolic differentiation is notdirectly suitable for dealing with non mathematical functions, like algorithms or computer programsthat might include for loops or if-then-else constructs that are common practice in the computerimplementation of numerical problems.
An alternative approach for calculating the numerical values of the derivatives of a function is through
Automatic Differentiation (AD). AD is based on the idea that any mathematical expression is evalu-ated by computers as a sequence of elementary algebraic operations, or call to mathematical functions,with the result being accumulated as the sub-expressions are evaluated [Margossian, 2018]. In contrastto finite difference that is based on the discretization of the derivative operator, AD techniques do notry to approximate the differentials but numerically evaluate the derivatives of the sub-expressions,alongside the value of the function itself, using the analytic rules of calculus.Therefore, if the derivatives of all the functions and operators used in the main expression are known, itis possible to evaluate the derivatives of the results with respect to the operands alongside to the valueof the function itself. As a result, AD is not affected by round off errors and allows for the numericalevaluation of derivatives within the same accuracy of the function itself, with the error only beinglimited by the machine’s representation of floating point numbers. In addition, AD can be appliedrecursively to evaluate higher order derivatives with no accuracy detriment, or error build up, becauseof re-using noisy estimates of lower order derivatives.Two general approaches for the implementation of AD are mostly employed, forward mode and reversemode . In forward mode AD the derivatives of the function, with respect to each of its independentvariables, are evaluated along with the main function, in the same order they are encountered, for eachsub-expression. Accordingly, the cost for derivatives evaluation is roughly the cost of one function callper independent variable. In reverse mode AD the evaluation of the function and of its derivativestakes place in two separate steps. First, the main function is parsed into sub-expressions that are indi-vidually evaluated, whose result is stored, alongside with their derivatives with respect to the argumentof each individual sub-expression. Secondly, the derivatives of the function value with respect to theindependent variables are reconstructed using the intermediate derivative of the sub-expressions thathave been evaluated, and stored, in the first step. Thus, at the end of the second step, all the partialderivatives are available in one full run. Because of its structure, the cost of the reverse mode AD isa few times the cost of the function evaluation alone, depending on how the main expression tree isstructured and on how the partial sub-expressions are interconnected, and it is essentially indifferentto the number of the partial derivatives required. In typical situations we can expect roughly up toa few tens ot times the cost of the evaluation of the function value alone, for a number of partialderivative that can be in the order of more than a few hundred thousands. Nonetheless, reverse modeAD suffers from the necessity of allocating, and keeping available to the CPU, all of the intermediateresults until the entire function is evaluated. Griewank and Walther [2008], Hogan [2014] and Elliott[2018] provide detailed and thorough description of the forward and the reverse mode AD algorithms.Both forward and reverse mode AD have been traditionally implemented as algorithmic differenti-ation techniques [Bischof et al., 1996, Corliss et al., 2002, Naumann, 2012, Forth et al., 2012], whichwould take a function’s source code as input, and produce the source code for the derivative of thefunction as output. While reverse mode AD necessarily requires to operate on the function’s sourcecode, forward mode AD also allows a different type of implementation. Since forward mode AD isbased on a single passage, and it only requires to access the value and the derivatives of the operands at7ach operation singularly, it can also be implemented by purposely defining a data type that is capableto store both the value and the derivatives of a variable. Therefore, in the programming languagesthat allow to extend the ordinary maths operators and function to deal with user defined data types,through a functionality known as operator overloading , once the arithmetic of the extended numericaltypes is defined, it is possible to call the same code that evaluates a function on the ordinary floatingpoint numeric types, with the newly defined data types, and obtain both the value and the derivativesof the result. Appendix A illustrates one application of such technique in the Julia programminglanguage [Bezanson et al., 2017, Perkel, 2019].In this paper we implement forward mode AD through operator overloading. Therefore, in the caseof conservative systems, it is only necessary to write the code for evaluating the free energy at the in-tegration points of the elements, and the result will also include its gradient and Hessian, that coincidewith the element nodal forces and stiffness matrix, respectively. The choice of operator overloadingforces us to use forward mode AD, which is less efficient then reverse mode AD when the number ofindependent variables increases. Nonetheless, since we operate on the model element-wise, we onlydeal with a reduced number of DoFs each time, for instance a QUAD element in 2D involves only 8DoFs, while a HEXA element in 3D involves 24 DoFs. At such number of DoFs forward mode is stillmore efficient than reverse mode because of the reduced no overhead cost needed for preprocessing,and the absence of additional costs for memory allocation and garbage collection, due to the storageof the intermediate results until the end of the function.
The particular implementation of the forward mode AD adopted here is based on the recurs to dualnumbers, an enriched number system, whose elements have multiple, higher dimensional parts thatcan deal with the derivative information up to a desired order. Dual numbers, together with therelated algebra, extend real numbers in a way similar to complex numbers and quaternions. As thecomplex field allows dealing with expressions that include the square root of negative numbers, andthe arithmetic of the quaternions allows to simplify the treatment of rotation in three dimensions, thearithmetic of dual numbers allows the simultaneous calculation of the value of a given expression andits derivatives up to an arbitrary order.In introducing the dual numbers and their properties, we follow the general treatment of higherdimensional number systems as given by Shenitzer et al. [2011]; a similar treatment of the subject isgiven by Fike and Alonso [2011]. Howbeit, the treatment of the cited authors did not cover numberswith multiple, separate, higher-order components, with different dimensionality, while the treatmentpresented in this paper makes use of such structures to deal with derivatives with arbitrary differenti-ation order. For the sake of simplicity, in this sections we refer to dual numbers of the second order,and we leave to Appendix B the generalization to dual numbers of arbitrary order.In the present study, dual numbers are the structures chosen to store, and operate, both on thevalue of a given parameter, x , and on its derivatives with respect to the independent variables of theproblem. We define dual numbers of the second order and dimension N the quantities of the kind x ≡ x + x i ı i + x ij ı ij with i ∈ . . . Nj ∈ i . . . N , (19)where x is the value taken by x , x i are the values taken by the first derivatives of x with respect tothe i -th independent variable, and x ij are the second derivatives of x with respect to the i -th and the j -th independent variable; furthermore, the symbols ı i are the elements of the canonical base of R N ,while ı ij are defined as ı ij ≡ ı i ⊗ ı j + ı j ⊗ ı i , (20)8ith ı i ⊗ ı j being the dyadic product of ı i over ı j , thus the following statement holds ı ij ≡ ı ji , (21)which translates the symmetry of the Hessian and, more in general, the independence of higher deriva-tives from the order of differentiation. In the following we refer to the first summand in equation (19), x , as the real part of x , and to the second summand, x i ı i , as its first order dual part and to x ij ı ij asthe second order dual part or x . As we can observe, the first order dual part of x is the vector spaceof real numbers of dimension N , while x ij ı ij is the vector space of real symmetric square matrices ofdimension N.We will now show that over the set of dual numbers, it is possible to define the operations ofaddition as well as subtraction, and the operation of multiplication and division. The neutral elementfor the sum is the dual zero, i.e. the dual number whose components are all nought, while the neutralelement for the product is the dual unity, i.e. the dual number whose value is one, and the derivativescomponents are all nought. Therefore, the set of dual numbers as defined above is a field, and the dualzero and unity coincide with the zero and the unity of the real field.By definition, two dual numbers are equal if all of their components are equal, therefore, thefollowing equality statement for dual numbers holds x = y ⇐⇒ y = x y i = x i y ij = x ij ∀ i ∈ . . . N ∀ j ∈ i . . . N . (22)For the sake of brevity of notation we will omit to specify the limits value for i and j in the following.We define the sum/the difference of two dual numbers as the sum/the difference of their components,component by component, as follows z = x ± y ⇐⇒ z = x ± y z i = x i ± y i z ij = x ij ± y ij . (23)As we can observe, the above definition of sum also induces the definition of the opposite of a dualnumber as the symmetric with respect to zero for the sum operation. In a similar way, we define theproduct of two dual numbers as the sum of all the mixed products of their components, where thefollowing product rules apply for the symbols ı i and ı ij ı i ı j ≡ ı ij (24a) ı ij ı k ≡ , (24b)with ı ij defined in equation (20). With the above positions, the product of two dual numbers followsas z = xy ⇐⇒ z = x y z i = x i y + x y i z ij = x ij y + x i y j + x j y i + x y ij , (25)we observe that the product of two dual numbers is commutative and associative, thus the followingholds xy = yx (26) x ( yz ) = ( xy ) z . (27)9s we can observe, given the above definition of multiplication we can define the quotient of two dualnumbers z = x / y as the solution to the following equation zy = x , therefore the following holds z = xy ⇐⇒ yz = x ⇐⇒⇐⇒ z = x y z i = x i y − x y i y z ij = x ij y − x i y j + x j y i y + 2 x y i y j y − x y y ij , (28)which also defines the inverse of a dual number, as the symmetric to unity with respect to multipli-cation, obtained by replacing x with 1 in the above equation. We observe that, with the definitionsgiven above, the dual numbers field, similarly to the complex number field, is an associative algebra.We now observe that, by recursively applying the identity (25), it is possible to extend the operationof raising to integer power over the field of second order dual numbers as y = x n ⇐⇒ y = x n y i = n x n − x i y ij = n ( n − x n − x i x j + n x n − x ij . (29)More in general, it is possible to extend any continuous, twice differentiable, function f ( x ) over thesecond order dual number field by making use of the chain rule for the derivatives as follows f ( x ) = f ( x ) + ∂f∂x x i ı i + (cid:18) ∂ f∂x x i x j + ∂f∂x x ij (cid:19) ı ij , (30)where all the derivatives are evaluated in x = x .As an example let’s assume that we are interested in evaluating the expression y ( x , x , x ) = x x + x (31)over the second order dual number field, i.e. by treating x i as second order independent dual quantities,we will show that the result will be a dual quantity itself, retaining the value of the function, and itsderivatives with respect to the x i , up to the second order. Assuming that x , x and x are theindependent variables, by definition their first derivative with respect to themselves is one, and anyother derivatives is zero, therefore their dual representation is the following x = x + ı x = x + ı x = x + ı , (32)and the expression (31), evaluated as a dual quantity, takes the following value y = ( x + ı ) ( x + ı ) + ( x + ı ) == ( x + 3 x ı + 3 x ı )( x + 2 x ı + ı ) + x + 2 x ı + ı == x x + x + 3 x x ı + 2 x x ı + 2 x ı + 3 x x ı + 6 x x ı + x ı + ı == y + y i ı i + y ij ı ij , (33)10ith: y = x x + x ,y i ı i = 3 x x ı + 2 x x ı + 2 x ı ≡ x x x x x ,y ij ı ij = 3 x x ı + 6 x x ı + x ı + ı ≡ x x x x x x x
00 0 2 , as we can observe y i and y ij coincide with the gradient of y with respect to x i and with its Hessian,respectively. Appendix A shows an implementation of the above example in the Julia programminglanguage, carried out with numerical values for x , x and x .The dual numbers system can be readily implemented in the programming languages that allowthe users to define data types and overload of existing arithmetic operators over the newly definedtypes. The dual numbers type should include data members to hold the real value of the number,and as many higher dimensional arrays up to the desired order of differentiation, with the number ofelements in each dimension equal to the number of independent parameters the real part of the numberdepend on. Once the data type and the operators have been implemented, scripts with mathematicaloperations carried out on dual numbers take the same form of the scripts operating on real numbers,and virtually no change is necessary. The source code of all the script developed for producing theresults presented in this paper are available from the web repositories indicated in [Vigliotti, 2020]. In this section we will implement AD for the solution of solid mechanics problems. In first placewe will recap how a continuum mechanics problem is generally stated and solved within the FEframework. We will then use the same FE discretization for evaluating only the free energy of a solidin a given configuration, and we will recognize how the Jacobian and the Hessian of the system’s freeenergy coincide with the residual force vector and the tangent stiffness matrix, respectively. Thus, wewill show how the complexity of the direct calculation of the residual force vector and the tangentstiffness matrix in FE contrasts with the simplicity of the calculation of the free energy alone. We willalso discuss the implementation, within the AD-assisted framework, of some important components inclassical solid mechanics problem, such as non trivial boundary conditions and the hyperelastic materialmodels, while the treatment of geometric non linearity is essentially built-in the AD framework.
The FE method is generally presented starting from the statement of the principle of virtual work(PVW) [Bathe, 2014, Zienkiewicz et al., 2005, Bonet and Wood, 2008], which is equivalent to equations(7) and (9), in the reference and the current configuration, respectively. In particular, with referenceto (7a), the PVW can be expressed, in terms of quantities defined on the reference configuration, as (cid:90) V P ij δF ij d V = (cid:90) V ∂b k ∂u k δu k d V + (cid:90) S ∂t k ∂u k δu k d S ∀ δu k with δF ij = ∂F ij ∂u k δu k , (34)where δu k is a virtual displacement field compliant with the boundary conditions, δF ij is the corre-sponding virtual deformation gradient, and it has been made use of the positions (14) for the remainingsymbols. The FE approach consists in approximating the evaluation of the integrals in equation (34)by discretizing both the domain of integration and the functional space over which the solution is11ought. The domain of integration is partitioned into elements connected in nodes, and the solutionis expressed as the weighted sum of a finite set of shape functions. The shape functions are definedover the elements, and are selected in a way that guarantees a number of requisites, such as adequatedifferentiability, continuity across element boundaries and convergence to the analytical solution asthe size of the elements collapses. To this end, the classical choice in FE approaches are multivariatepolynomials, which are used to interpolate the nodal values of the unknown function on the interiorof each element. Nonetheless other options are possible, such as the Isogeometric Analysis approach[Hughes et al., 2005] where the geometry is described by Non-Uniform Rational B-Splines, and thesame rational functions used for the geometry take the role of the shape functions that interpolate thesolution over the domain, while the control weights coincide with the coordinates of the vertices of thecontrol polygon for the unknown field.Once the domain has been discretized in elements, and a suitable set of shape functions has beenselected, the integrands in equation (34) are a function of a discrete number of degrees of freedom only.In this framework the components of the displacement field, and of the displacement gradient at anypoint of a given element of the domain can be written as u i = N i · u (35a) F ij = N i,j · u + δ ij (35b)where u is the array of the DoFs, N i is the array of the shape functions for the i -th componentof the displacement field, · is the dot product, and N i,j is the j -th component of the gradient of N i , as usual. We remark that in this section vector and matrices are marked in bold face, and theproduct to two vector quantities should be interpreted as the dyadic product, whose result is a matrix.Therefore, for a given u , given N i,j , it is possible to evaluate u i , F ij and any quantity depending onthe displacement field, at any point of any element of the FE model of the domain. In the same waythe virtual displacement field can be obtained by means of the same interpolation as δu i = N i · δ u (36a) δF ij = ∂F ij ∂ u · δ u = N i,j · δ u , (36b)where the arbitrariness of δu k over the functional space of the displacement fields that are compatiblewith the boundary conditions of the problem translates into the arbitrariness of the components of δ u .Thus, by means of equations (35) and (36) it is possible to rewrite equation (34) as (cid:90) V P ij ∂F ij ∂ u d V = (cid:90) V ∂b ∂ u d V + (cid:90) S ∂t ∂ u d S , (37)which is a system of non linear equations, in the unknown unconstrained components of u . Therefore,the differential problem of the equilibrium, as stated in equation (1), is translated into a system ofnon-linear equations, where the unknowns are represented by the unconstrained DoFs. The integralsin equation (37) can be numerically evaluated by means of quadrature rules, and the solution of theFE problem can be found as the zero of the residual vector given by the following r = N BE (cid:88) m =1 N mBW (cid:88) l =1 w ml (cid:20) P ij ∂F ij ∂ u − ∂b ∂ u (cid:21) r ml − N SE (cid:88) m =1 N mSW (cid:88) l =1 v ml (cid:20) ∂t ∂ u (cid:21) r ml = , (38)where N BE is the number of volume elements, N mBW is the number of integration points of the m -thvolume element, w ml is the l -th volume integration weight of the m -th element; while N SE is thenumber of surface elements, v ml is the l -th surface integration weight of the m -th surface element, andthe subscripts of the square bracket indicate that the quantities enclosed are evaluated at the point r ml , the position of the w ml integration weight. 12n the finite element formulation equation (38) can be solved through Newton-Raphson iterativeschemes, after an expression for the Jacobian of r , or the tangent stiffness matrix, has been obtainedby differentiating equation (38) with respect to u , as follows ∂ r ∂ u = N BE (cid:88) m =1 N mBW (cid:88) i =1 w ml (cid:20) ∂P ij ∂F hk ∂F hk ∂ u ∂F ij ∂ u − ∂ b ∂ u ∂ u (cid:21) r ml − N SE (cid:88) m =1 N mSW (cid:88) i =1 v ml (cid:20) ∂ t ∂ u ∂ u (cid:21) r ml , (39)where we made us of the fact that, since F ij is linear in u , the following holds ∂ F ij ∂ u ∂ u = . (40)The calculation of the summands in equation (38) and (39) is the core of the FE methods and representthe most challenging aspect of the computer implementation of the method. In particular the calcu-lation the components of the stress tensor, P ij , and their derivatives ∂P ij /∂F hk , is in general a quitesophisticated task, since it requires dealing with second order and fourth order tensors respectively. We now turn our attention to the use of the AD for the solution of the equilibrium problem. We beginby assuming that the material of the solid is a Green elastic material, because under this assumptionthe resulting formulation is particularly simple and insightful. For Green elastic solid the deformationwork is an exact differential, and the free energy function for the system is given by equation (5).In order to numerically evaluate the integrals in equation (5) we can make use of the same twofolddiscretization used for FE, obtaining the following expressionΨ( u ) = N BE (cid:88) m =1 N mBW (cid:88) i =1 w ml [ φ + b ] r ml + N SE (cid:88) m =1 N mSW (cid:88) i =1 v ml [ t ] r ml , (41)where it has been highlighted that the free energy, within the FE discretization, is a function of thearray of DoFs, u . As discussed in section 2, equilibrium configurations are those that satisfy equation(6), which in the FE discretization can be written as δ Ψ = ∂ Ψ ∂ u · δ u = 0 , ∀ δ u ⇐⇒ ∂ Ψ ∂ u = . (42)The expression above is a system of non-linear equations in the unknown u , whose residual andJacobian, are given as the gradient and Hessian of Ψ, respectively, as r = ∂ Ψ ∂ u (43a) ∂ r ∂ u = ∂ Ψ ∂ u ∂ u . (43b)We remark that the residual on equation (38), which derives from the PVW statement given in (34),and the residual on equation (42), which is obtained as the first the deformation work in the equilibriumconfiguration given in (41), are both work-conjugated through the same virtual nodal displacements, δ u , hence they must coincide. Therefore the gradient of the free energy coincides with the residualforce vector of the finite element problem. At the same time, the tangent stiffness matrix coincideswith the Hessian of the free energy, being both the derivative of r with respect to u .We now remark that the gradient and the Hessian of Ψ, in equations (43a) and (43b), respectively,can be both immediately calculated by the same computer program that evaluates equation (41),13hrough automatic differentiation, if u is treated as an array of dual numbers, and the dual numberalgebra has been implemented in the programming language. Therefore, AD allows to evaluate theresidual force vector and the tangent stiffness matrix by simply calculating the numerical integral ofthe free energy density over the domain.We finally remark that the same approach can be used in the cases when a functional relation existsbetween P ij and F ij , but no elastic potential can be defined. For these materials equation (43b) isreplaced by (38), through equations (35) and (36), with the components of u being independent dualquantities, while equation (43b) still holds and it is obtained as the first order dual components of r .In the sections that follow we will illustrate how some of the fundamental elements in a solidmechanics problem, such as non linear constitutive laws, or complex boundary conditions, can be easilyincluded in the problem formulation with recurs to automatic differentiation for their implementation. Non trivial boundary conditions can be applied with the use of Lagrange multipliers using automaticdifferentiation technique for the direct evaluation of the gradient and of the Hessian of the Lagrangefunction. We recall that, following the Lagrange multipliers technique, the minimization of a function,in the presence of constraints can be achieved by weighting the residuals of the constraint equationsthrough unknown factors, the Lagrange multipliers, and adding them to the function to be minimized,as follows L ( u , λ ) = Ψ ( u ) − λ · g ( u ) , (44)where Ψ is the function to be minimized in the first place, which in our case is the free energy of thesolid, u are the degrees of freedom of the problem, g ( x i ) is the array of the constraint equations and λ is the array of the Lagrange Multipliers.Since the problem of minimizing L with respect to u and λ , is essentially identical to the problemof minimizing Ψ with respect to u only, it can be treated in the same way. Nonetheless, since L ( u , λ )is linear in λ , it is not necessary to treat λ as dual quantity, but suffices to evaluate Ψ ( u i ) and g ( u i )over the dual number field of u , while the components of the augmented gradient and Hessian of L can be obtained from the expression of the first and second variation of L , which is given as δL = (cid:18) ∂ Ψ ∂ u − λ · ∂ g ∂ u (cid:19) · δ u − g · δ λ (45a) δ L = (cid:18) ∂ Ψ ∂ u ∂ u − λ · ∂ g ∂ u ∂ u (cid:19) : δ u δ u − ∂ g ∂ u : δ u δ λ + − δ λ · ∂ g ∂ u · δ u , (45b)which yields the following, in block matrix notation, ∇ L = ∂ Ψ ∂ u − λ · ∂ g ∂ u − g (46a) ∇ L = ∂ Ψ ∂ u ∂ u − λ · ∂ g ∂ u ∂ u − ∂ g ∂ u T − ∂ g ∂ u . (46b)Therefore the problem reduces to solving the following ∇ L = 0 , (47)where equation (46b) takes the meaning of the tangent stiffness matrix of the problem.14 .4 Hyperelastic material models We now turn our attention to the most common expressions for the strain energy density functionsof materials. In very simple cases, such as for the small deformations of strut or beam elements, thestress is uniaxial, the deformation state of the solid is adequately described by a single scalar quantity,and the material behaviour can be treated as linear elastic, with a strain energy density function ofthe type Φ H = 12 E s (1 + (cid:15) n ) , (48)where E s is the Young modulus of the material, and (cid:15) n the component of the nominal strain tensorconjugated to only non zero stress tensor component . In these cases the calculation of the deformationwork is particularly simple; however, for a general solid mechanics problem, the state of deformation hasarbitrary principal directions and distinct principal stretches, therefore, more sophisticated expressionsfor the deformation energy density are used. Green elastic materials are a quite general class of materialmodels for which the strain energy density function is assumed as a local function of the componentsof the right Cauchy-Green deformation tensor, C = F T F [Ogden, 2013], as φ G = φ G ( C ij ) , (49)where C ij are the components of C . Among Green elastic materials, one class of material modelsthat are of great interest for engineering applications are the isotropic hyperelastic materials, whosebehaviour is invariant to rigid rotations of the applied strain. For these material models, the strainenergy density function can be expressed as a function the invariants of C only, defined as φ iso = φ iso ( I , I , I ) with: I = C + C + C I = C C + C C + C C − C − C − C I = C C C + 2 C C C − C C − C C − C C (50)Incompressible isotropic materials are subject to the internal isochoric constraint, I = 1, therefore φ iso depend on I and I only, and its general expression is of the type [Ogden, 2013] φ inc = ∞ (cid:88) p,q =0 c pq ( I − p ( I − q with I = 1 , (51)where p and q are non negative integers and c pq are non negative real parameters. Expression (51) isonly valid if the isochoric constraint is explicitly enforced, however, since such a constraint often yieldsto convergence problems in the finite element formulation, when using an incompressible materialmodel, the following decomposition of the C can be assumed, C = CC v where C v = J / I C = J − / C , (52)where J = det ( F ) and I is the identity tensor. It is straightforward to verify that det (cid:0) C (cid:1) = 1,therefore equation (52) decompose the total deformation into a isochoric deformation, representedby C , and a purely volumetric deformation, given by C v . Following the decomposition (52), Φ inc isapproximated as Φ = ∞ (cid:88) p,q =0 c pq (cid:0) I − (cid:1) p (cid:0) I − (cid:1) q + f ( J ) , (53)where I are the invariants of C , and f ( J ) is a positive function of J that effectively penalizes volumevariations. In the example section of this paper we will use, in particular, the Mooney-Rivlin and the15eo-Hokkean models, whose strain energy density expression is given by Neo-Hookean φ NH = c ( I −
3) + G ( J − (54a) Mooney-Rivlin φ MR = c ( I −
3) + c ( I −
3) + G ( J − (54b)where c , c and G are constant, non negative, parameters that define the material behaviour. We remark that in the approach presented here, since the residual force vector and the tangent stiffnessmatrix are automatically obtained from the free energy, we never explicitly calculate the componentsof the stress tensor. However, the value of the entries of the stress tensor are still important quantities,since resistance criterion, such as Von Mises, are based on it. Nonetheless, they can always be evaluated,as a post processing step, from the equalities (10a) and (14a), by means of the applicable expressionfor the strain energy density, by treating the components of F as independent dual quantities, whosevalue is obtained from the displacement field of the equilibrium configuration. In this section we present a selection of solid mechanics problems whose solution has been found withthe recurs to the automatic differentiation techniques described in the paper. The examples presentedinclude structural elements, such as rods and beams (section 5.1 and 5.2), continuous plane stresselements (section 5.3), a problem with cylindrical symmetry (section 5.4), and a full three-dimensionalproblem (section 5.5).All the problems included the effects of geometric non-linearities, non trivial boundary conditions, andthe hyperelastic material models described in section 4.4. The non trivial boundary conditions wereintroduced using the Lagrange multipliers technique, as described in section 4.3.All problems presented here were solved using the Julia programming language [Bezanson et al., 2017,Perkel, 2019], and the script files used for the solution have been made available to the reader [Vigliotti,2020].
In this section we consider the equilibrium of a tridimensional structure made of prismatic elements,connected at their endpoints to form a truss. We also assume that the material of the struts is linearelastic, with Young modulus E s , and that cross section deformations are negligible with respect tothe axial deformation of the elements. Under these assumptions, the elements can only store elasticenergy by variations of their length, and the deformation energy of a single element, is given as φ rod = A l Φ H , (55)where A is the cross section area, l is the reference length, Φ H is defined by equation (48), with (cid:15) n = l/l −
1, and l is the length of the element in the current configuration. The total deformationenergy of the truss can then be readily obtained as the sum of the deformation energy of all of itselements, and it is given as Φ truss = (cid:88) i φ rod i , (56)16 r u u i i i (a) i i i (b) i i i (c) Figure 1: (a) Rod element, r and r are the nodes positions in the reference configuration, u and u are the displacement vectors. (b) Sketch of the regular octet unit cell, all struts have the same length L .(c) Sketch of the entire truss structurewhere φ rod i is the strain energy of the i − th element. With reference to figure 1.a, let r and r be thepositions of the end nodes of a rod in the reference configuration, and u and u be the displacementvectors of the nodes, l and l can be easily obtained for a given element as l = (cid:107) r − r (cid:107) l = (cid:107) r + u − ( r + u ) (cid:107) . (57)Therefore, since a truss can be idealized as network of rods connecting in nodes, given the topologyof the connections, and the cross section and material properties of the struts, the deformation en-ergy of a truss can be easily computed as function of the components of the displacements of the nodes.Here we consider the equilibrium of a structure obtained by replicating the regular octet unit cell,shown in figure 1.b along the directions ı , ı and ı , without duplicating the coincident rods. Theregular octet is a well known structure, which is characterized for its lightness and strength [Fuller,1966, Deshpande et al., 2001a]. Since the regular octet topology is both statically and kinematicallydetermined, no mechanisms arise from its deformation, and it can withstand any external load bystretching of its elements only [Deshpande et al., 2001b]. In particular we consider a structure madeof 2 unit cells in the directions ı , N = 2, three units in direction ı , N = 3, and 10 cells in direction ı , N = 10, as shown in figure 1.c.The selected boundary conditions produce the bending of the domain around an axis parallelto ı by constraining the nodes on both the the ends of the domain to remain on planes that aresymmetrically rotated around an axes parallel to ı , for a prescribed angle θ . Under such boundaryconditions the rigid body translation along direction ı is still available to the structure, and it isremoved by constraining any motion of the center of gravity in direction ı . In summary, the boundaryconditions are implemented through Lagrange multipliers and are expressed as( X + u ) cos ( θ ) + ( X + u ) sin ( θ ) − L = 0 (58a)( X + u ) cos ( θ ) − ( X + u ) sin ( θ ) + 12 L = 0 (58b) (cid:88) u = 0 , (58c)17igure 2: Shape of the lattice structure for different rotation angle of the top and bottom planes. Therotation on the top and the bottom planes have equal amplitude and opposite sign.where the condition (58a) applies to the nodes at X = L /
2, condition (58b) applies to the nodes at Z = − L /
2, while the summation in (58c) extends to all the nodes of the model and serves the purposeof removing the residual rigid body degrees of freedom. The model comprised 368 nodes and 2160rods, boundary conditions were applied in π/ π/ π/
8. In all cases convergencewas achieved in 5 iterations, except in steps 6 and 7 where convergence was achieved in 6 and 7 stepsrespectively, which confirmed that the Hessian of the free energy was accurate. Simulation results fordifferent values of θ between 0 and 7 π/ In a similar way it is possible to analyse the response of structures made of Euler beams. Given aprismatic bar, the Euler beam model assumes that each cross section rigidly rotates around an axisorthogonal to the beam axis, passing through the centre of gravity of the section, neglecting any shearcontribution to deformation energy and load bearing. With reference to figure 3, in a reference framewith ı aligned with the beam axis, and ı aligned with the axis of rotation of the cross section, underthe assumption of small local cross section rotation, the following displacement model for the pointsof the beam holds u = u − X u , u = u , (59)where u i are the components of the displacement of the points on the beam axis, X i are the coordinateof the points of the beam in the reference configuration, and the subscript after the comma denotesdifferentiation with respect to the coordinate X i . Since the shear contributions to deformation energyare neglected, under finite displacement assumption, the elastic energy is only a function of the firstcomponent of the Green Lagrange deformation tensor, which, in accordance to equation (59), is givenas E = 12 (cid:104) (1 + u , ) + u , − (cid:105) == 12 (cid:104) (1 + u , − X u , ) + u , − (cid:105) . (60)Under the assumption that a linear elastic model is adequate for representing the material behaviour,the deformation energy of an individual beam and of a structure made of beam elements are given18 eformed shapereference shape i i Figure 3: Euler beam displacement modelrespectively as φ i = 12 E s (cid:90) V i E d V (61)Φ = (cid:88) i φ i . (62)Therefore, once a parametric representation for the displacement of the points of the axis of the beamsis given, the deformation energy of the structure can be expressed as a function of the chosen parame-ter, and the equilibrium configuration can be found by minimizing the free energy of the structure. Acommon parametric representation for u assumes the nodal displacement and rotations as the degreesof freedom of the elements, and takes the axial component of the displacement as a linear function of X , and the transverse components as cubic functions of X [Bathe, 2014].Figure 4 shows the equilibrium configuration of a bidimensional hexagonal lattice under tensionobtained using equation (62) for evaluating the deformation energy and the automatic differentiationapproach described in this paper to find the stationary energy configurations. Displacement boundaryconditions have been applied to the nodes on the top and bottom of the model, by preventing horizontaldisplacement and rotation around the axis orthogonal the plane of the model, and prescribing thevertical displacement of the top nodes. All elements have the same length, L , and the same squarecross section with side t = L /
10. Figure 4a shows the deformed configuration for a prescribed totaldisplacement in the vertical direction of ∆ L = 0 . L , where L is the initial length of the model inthe direction ı , while figure 4b shows the total reaction force, obtained as the sum of the residualsconjugated to vertical displacement for the top nodes of the model, normalized by the Young Modulusof the material, E s , and the cross section area, A = t , as a function of the applied displacement. Thesolution was obtained in 20 steps, with every step taking between 6 and 9 Newton-Raphson iterations. In this section we consider the equilibrium of an hyperelastic plate subjected to in-plane boundaryconditions. In particular we consider the domain, and the boundary conditions shown in figure 5, wherea vertical displacement ∆ u is applied to the points of the top boundary, while the displacements ofthe points of the bottom boundary are fixed. The domain features two types of internal boundaries.The internal boundaries with radii R are empty cavities, where standard traction free boundaryconditions applied. The internal boundaries with radii R I are rigid, circular, frictionless inclusions,whose points are constrained to remain at a fixed distance from the centre of the inclusion, which isfree to move in both directions. The material of the domain is assumed to be a Mooney-Rivlin type themodulus c = 10 and G = 10 . The presence of the inclusions has been introduced through Lagrangemultipliers, by constraining the nodes lying on each inclusion boundary to remain on a circle having19 i (a) (b) Figure 4: Simulation results for a bidimensional hexagonal lattice made of Euler beam elements,the beam elements have square cross section with side t = L /
10. (a) Deformed and undeformedconfiguration of the lattice for total applied. (b) Normalized total reaction force vs applied displacementcurve, R is the sum of the residual conjugated to the vertical component of the displacement of thetop nodes, E s is the Young Modulus of the material, A = ( L / is the cross section area of theelements. i i L = 20 L = 20 Inclusions, R I = 1.25Openings, R O = 1.50 Figure 5: Plane stress problem, undeformed domain with boundary conditions, ∆ u is the applieddisplacement, the radius of the openings is R O = 1 . R I = 1 . R I , whose centre’s coordinates were an unknowns, introduced through additional boundaryconditions, determined at the equilibrium. Figure 6.a shows the plot of the normalized constraintreaction as a function of the displacement, while figure 6.a shows the equilibrium configuration. Asit can be observed, while the cavities dramatically changed their shape, both in compression and intension, the inclusions maintained their circular shape. The solution was achieved in 150 incrementsfor the compressive branch and 100 increments for the tensile branch, with each increment convergingin 4 or 5 iterations. In this section we discuss the solution of problems with cylindrical symmetry. In the cases wherethe geometry of the domain, the material, and boundary conditions all have cylindrical symmetry, theproblem can be significantly simplified incorporating the symmetry conditions within the displacementmodel. Assuming that the axis of symmetry coincides with the ı axis, in the absence of torsion, the20 a) (b) (c) Figure 6: Plane stress problem, simulation results. (a) Normalized force-displacement plot, ∆ u is thevertical displacement applied to the top side of the boundary, R is the total reaction force conjugated to∆ u , L is the side of the domain, t is the thickness. Deformed configuration at maximum displacementin compression (b) and tension (c), c is the modulus of the Neo-Hookean material model, the colormapis based on the first invariant of the deformation tensor.deformation gradient takes the following form F cyl = u , u , u , u ,
00 0 1 + u X , (63)and it is invariant to rotations around the symmetry axis.In this example we consider the domain shown in figure 7a, obtained by completely a completerotation of the section shown in figure 7b around the axis ı . We assume that the material followsa Neo-Hookean model with modulus c = 10 .
0, as per equation (54a). In addition to the boundaryconditions illustrated in figure 7b, we consider two cases for the behaviour of the internal cavity. Inone case we assume that the cavity is filled with an incompressible fluids. For this case no particularshape is enforced, but only the value of the cavity’s volume is kept constant during the solution. Inthe second case we consider the cavity as an empty volume that can take any shape, with no otherconstraint than the external boundary conditions.For any given configuration of the body, the volume of the internal cavity is given as V c = 2 π (cid:90) Σ c x dΣ = π (cid:90) Γ c x x d x (64)where Σ c is the intersection of the internal cavity volume with the plane ı − ı in the initial configu-ration, and Γ c its boundary, as shown in figure 7b. We recall that the second equality in equation (64)holds thanks to the Gauss-Green theorem, and allows replacing the area integral with a curvilinearintegral, without the need of discretizing the interior of the cavity. The constraint on the volume of theinner cavity can then be introduced by means of Lagrange multipliers by requiring that the followingholds V c ( u k ) = V c , (65)where V c ( u k ) is the current volume of the cavity and V c is the initial volume. Therefore, the expressionof the Lagrange functional to be minimized in order to solve the problem is the following L cyl ( u k , λ c ) = Φ cyl ( u k ) − λ c [ V c ( u k ) − V c ] , (66)where Φ cyl is the deformation energy of the body, given asΦ cyl = 2 π (cid:90) Σ φ (cid:0) F cyl (cid:1) X dΣ , (67)21here Σ is the section of the domain on the plane ı − ı , X , the first coordinate, is the distance fromthe rotation axis, and φ is the strain energy density of the material. Figure 8 shows the results of the i i i (a) i i R =50 u D =20 D =50 (b) Figure 7: Axi-symmetric problem, domain geometry, dimensions and boundary conditions, R is thedistance of the centre of the section from the symmetry axes, the fillet radius are 2.0 and 5.0. Thethicker line in (b) marks the portion of the boundary where displacement constraints were applied.simulations for the axisymmetric problem for a Neo-Hookean material model with modulus c = 10.As expected, the presence of an incompressible fluid within the cavity, introduced through the con-straint (65), results in a general macroscopic stiffening of the solid, which enforces a lager widespreadof the deformation across the domain, and lager average value of the deformation. In both cases thesolution was achieved in 200 steps with each steps taking 5 iterations to converge in the case with thecavity volume constraint, and 6 iterations in the cases without the constraint. (a) (b) (c) Figure 8: Axi-symmetric problem, simulation results. (a) Normalized force displacement plot, ∆ u is the displacement applied in direction 2, L is the initial height of the domain, R is the totalreaction force, A is the area over which the boundary condition is applied, c is the modulus of theNeo-Hookean material. Line 1. is the response with internal volume constraint (i.v.c.), line 2. is theresponse without i.v.c. . (b) Deformed cross section with i.v.c., and (c) without i.v.c. at maximumdeformation.We remark that the cylindrical symmetry, in the above example, was introduced by simply incor-porating it in the displacement model and in the definition of F , in equation (63), while at no point, inthe statement of the elastic problem, it was necessary to express the equilibrium equation in cylindrical22oordinates. In this section we analyse the response of a three dimensional hyperelastic solid undergoing largedisplacements. Figure 9a shows the domain geometry, which consists in a right-handed helicoidal solidwith circular cross section. The helix has radius R e = 20, pitch p = 20 and height h = 40. In thisexample the boundary conditions have been applied constraining the displacement of the centre ofmass of the two ends of the helix to move along direction ı , increasing the height of the helix, hasshown in figure 9a. Therefore the following set of equation was imposed on the nodes of the end crosssections (cid:90) A btm u d A = 0 , (cid:90) A top u d A = 0 (68a) (cid:90) A btm u d A = 0 , (cid:90) A top u d A = 0 (68b) (cid:90) A btm u d A = − ∆ h , (cid:90) A top u d A = ∆ h A top and A btm are the top and the bottom cross section, respectively. The boundary conditionshave been applied, similarly to the previous example, through Lagrange multipliers, adding a constraintequation for each of the equations 68. We observe that the value of Lagrange multipliers conjugatedto equations (68c), at equilibrium correspond to the total constraint reactions in direction ı on thebottom and the top faces, respectively, which are the active forces, producing the deformation of thespring. Because of the symmetry of the domain and of the boundary conditions, at equilibrium boththe active reaction forces have the same value, R . Figure 9b shows the normalized plot of R as afunction of ∆ h , and figure 9c shows the deformed configurations of the helix corresponding to the pointsmarked in figure 9b. As we can observe the plot of the reaction force shows the expected hardeningbehaviour due to the alignment of the helix with the applied force. The solution was obtained in 300increments, with each increment taking 3 iterations to converge. i i i A btm h A top (a) (b) i i abcdef (c) Figure 9: Three-dimensional solid. (a) Domain’s geometry, A top and A btm are the top and the bot-tom end cross sections of the helix, respectively; (b) Reaction force, R , A is the cross section area,versus normalized applied displacement, c is the modulus of the Neo-Hookean model; (c) Deformedconfigurations at the stages of the simulation marked with dots in (b).23 Concluding remarks
Automatic differentiation (AD) techniques allow for the accurate and efficient numerical evaluation ofthe derivatives of a multivariate function. In this paper, AD has been used for stating and solvingnon-linear finite element solid mechanics problems. The approach presented here focuses in particularon Green elastic materials, for which the deformation work is an exact differential, and the solid can betreated as a proper conservative thermodynamic system. In these cases, the residual force vector andthe tangent stiffness matrix of the model coincide, respectively, with the gradient and the Hessian of thesystem’s free energy, which can be numerically evaluated through AD. Therefore, with the approachpresented here no explicit calculation of the stress tensor, nor of the elasticity tensor is required, norit is necessary to implement the complex kinematics that link the degrees of freedom of the model tothe internal forces and their derivatives. The same framework can also be applied with arbitrary, nonconservative, material models, although, here the explicit calculation of the components of the stresstensor is required, while the calculation of the elasticity tensor and of the tangent stiffness matrix canstill be automated. In the same way, sophisticated constraints equations and boundary conditions,can be introduced by means of Lagrange multipliers, and treated through AD. The method has beenpresented with a number of examples that illustrate the solution of selected non-linear problems,featuring hyperelastic material models, and complex constraints, along with the computer programsdeveloped for producing the results included in this article.
References
Robert Asaro and Vlado Lubarda.
Mechanics of Solids and Materials . Cambridge University Press,2006. doi: 10.1017/CBO9780511755514.S. Timoshenko and J.N. Goodier.
Theory of Elasticity . McGraw-Hill Book Company, 1987. ISBN978-0070701229.O.C. Zienkiewicz, R.L. Taylor, and J.Z. Zhu.
The Finite Element Method: Its Basis and Fundamentals .Elsevier Science, 2005. ISBN 0750663200.K.J. Bathe.
Finite Element Procedures . Klaus-J¨urgen Bathe, 2014. ISBN 9780979004957.O.C. Zienkiewicz and R.L. Taylor.
The Finite Element Method for Solid and Structural Mechanics(Sixth Edition) . Elsevier Butterworth-Heinemann, sixth edition edition, 2005. ISBN 0-7506-6321-9.Anders Logg. Automating the finite element method.
Archives of Computational Methods in Engi-neering , 14(2):93–138, Jun 2007. ISSN 1886-1784. doi: 10.1007/s11831-007-9003-9.Anders Logg, Kent-Andre Mardal, and Garth N. Wells, editors.
Automated Solution of DifferentialEquations by the Finite Element Method . Springer, 2012. ISBN 978-3-642-23098-1. doi: 10.1007/978-3-642-23099-8.J. Korelc and P. Wriggers.
Automation of finite element methods . Springer, 2016. doi: 10.1007/978-3-319-39005-5.J. Bezanson, A. Edelman, S. Karpinski, and V. Shah. Julia: A fresh approach to numerical computing.
SIAM Review , 59(1):65–98, 2017. doi: 10.1137/141000671.J.M. Perkel. Julia: come for the syntax, stay for the speed.
Nature , 572(7767):141–142, 2019. doi:10.1038/d41586-019-02310-3.Andrea Vigliotti. Automatic differentiation for solid mechanics in julia. https://github.com/avigliotti/AD4SM , 2019. 24.B Callen.
Thermodynamics and an Introduction to Thermostatistics . Wiley, 1985. ISBN9780471862567.J. Lyness and C. Moler. Numerical differentiation of analytic functions.
SIAM Journal on NumericalAnalysis , 4(2):202–210, 1967. doi: 10.1137/0704019. URL https://doi.org/10.1137/0704019 .J.N. Lyness. Differentiation formulas for analytic functions.
Mathematics of Computation , 22(102):352–362, 1968. doi: 10.1090/S0025-5718-1968-0230468-5.Charles C. Margossian. A review of automatic differentiation and its efficient implementation.
CoRR ,abs/1811.05031, 2018. URL http://arxiv.org/abs/1811.05031 .A. Griewank and A. Walther.
Evaluating Derivatives: Principles and Techniques of Algorithmic Dif-ferentiation . Other Titles in Applied Mathematics. Society for Industrial and Applied Mathemat-ics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104), second edition, 2008. ISBN9780898717761.R.J. Hogan. Fast reverse-mode automatic differentiation using expression templates in c++.
ACMTransactions on Mathematical Software , 40(4), 2014. doi: 10.1145/2560359.Conal Elliott. The simple essence of automatic differentiation.
Proc. ACM Program. Lang. , 2(ICFP):70:1–70:29, July 2018. ISSN 2475-1421. doi: 10.1145/3236765. URL http://doi.acm.org/10.1145/3236765 .Christian Bischof, Peyvand Khademi, Andrew Mauer, and Alan Carle. Adifor 2.0: Automatic dif-ferentiation of fortran 77 programs.
IEEE Comput. Sci. Eng. , 3(3):18–32, September 1996. ISSN1070-9924. doi: 10.1109/99.537089.G. Corliss, C. Faure, A. Griewank, L. Hascoet, and U. Naumann.
Automatic Differentiation of Algo-rithms: From Simulation to Optimization . Springer New York, 2002. ISBN 9780387953052.U. Naumann.
The Art of Differentiating Computer Programs: An Introduction to Algorithmic Dif-ferentiation . Software, Environments, and Tools. Society for Industrial and Applied Mathematics,2012. ISBN 9781611972061.Shaun Forth, Paul Hovland, Eric Phipps, Jean Utke, and Andrea Walther, editors.
Recent Advancesin Algorithmic Differentiation . Springer Science & Business Media, 2012.A. Shenitzer, I.L. Kantor, and A.S. Solodovnikov.
Hypercomplex Numbers: An Elementary Introductionto Algebras . Springer New York, 2011. ISBN 9781461281917.J. A. Fike and J. J. Alonso. The development of hyper-dual numbers for exact second-derivativecalculations. In
AIAA paper 2011-886, 49th AIAA Aerospace Sciences Meeting , page n.a., 2011.Andrea Vigliotti. Automatic differentiation for solid mechanics, example scripts, mendeley data, v2. http://dx.doi.org/10.17632/ybbsszpbss.2 , 2020.Javier Bonet and Richard D. Wood.
Nonlinear Continuum Mechanics for Finite Element Analysis .Cambridge University Press, 2 edition, 2008. doi: 10.1017/CBO9780511755446.T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs. Isogeometric analysis: Cad, finite elements, nurbs,exact geometry and mesh refinement.
Computer Methods in Applied Mechanics and Engineering ,194(39):4135 – 4195, 2005. ISSN 0045-7825. doi: https://doi.org/10.1016/j.cma.2004.10.008.R.W. Ogden.
Non-Linear Elastic Deformations . Dover Civil and Mechanical Engineering. DoverPublications, 2013. ISBN 9780486318714.R. B. Fuller. U.s. patent serial no. 2, 986, 241, 1966.25. S. Deshpande, N. A. Fleck, and M. F. Ashby. Effective properties of the octet-truss lattice material.
Journal of the Mechanics and Physics of Solids , 49(8):1747–1769, 2001a.V. S. Deshpande, M. F. Ashby, and N. A. Fleck. Foam topology: Bending versus stretching dominatedarchitectures.
Acta Materialia , 49(6):1035–1040, 2001b.J. Revels, M. Lubin, and T. Papamarkou. Forward-mode automatic differentiation in julia. arXiv:1607.07892 [cs.MS] , 2016. URL https://arxiv.org/abs/1607.07892 .26 ppendicesA Implementation of dual number systems in the Julia pro-gramming language
In this section we illustrate a possible implementation of the dual number system in the Julia pro-gramming language. Julia is a dynamically typed scientific programming language, whose semanticis particularly suitable for the description of physical problems [Bezanson et al., 2017, Perkel, 2019].Aside to user-defined data types, Julia permits to overload existing operators or functions to be eval-uated on the new types. Therefore the same script that evaluates a numeric function on floating pointnumbers, can be used to operate on dual numbers, once their arithmetic has been implemented, andproduce dual number as a result.The script blocks reported below show a possible implementation of dual numbers in Julia. Thedual number type is called D2 and it is defined in Listing 1 as having a scalar component v , that storesthe current value of the variable, a one dimensional array, d1 , that stores all of the first derivatives of v , and a two dimensional array, d2 , that stores all of the second derivatives of v . s t r u c t D2 < : Number v : : F l o a t 6 4 d1 : : Array { Float64 , 1 } d2 : : Array { Float64 , 2 } end Listing 1: Definition of dual numbers in JuliaThe code in Listing 2 extends some ordinary maths operators to function with the D2 type. The firstline in the script block informs the language that the scope of the mentioned operators, defined in theBase module, will be extended, while the remaining lines implement the arithmetic of dual numbersas defined in section 3.3.1, where each component of a dual number value is accessed through the dotsyntax ( . ), and the single quote ( ' ) denotes array transposition. import Base : + , − , ∗ ,/ ,ˆ +(x : : D2 , y : : D2) = D2( x . v+y . v , x . d1+y . d1 , x . d2+y . d2 ) − (x : : D2 , y : : D2) = D2( x . v − y . v , x . d1 − y . d1 , x . d2 − y . d2 ) ∗ ( x : : D2 , y : : D2) = D2( x . v ∗ y . v , x . d1 ∗ y . v+y . d1 ∗ x . v , x . d2 ∗ y . v+x . v ∗ y . d2+ x . d1 ∗ y . d1 ' +y . d1 ∗ x . d1 ' ) / ( x : : D2 , y : : D2) = D2( x . v/y . v , x . d1/y . v − (x . v/y . v ˆ 2 ) ∗ y . d1 , x . d2/y . v − (x . d1 ∗ y . d1 ' +y . d1 ∗ x . d1 ' ) / y . vˆ2+2x . v ∗ ( y . d1 ∗ y . d1 ' ) / y . vˆ3 − ( x . v/y . v ˆ 2 ) ∗ y . d2 ) ˆ ( x : : D2 , n : : I n t 6 4 ) = D2( x . vˆn , ( n ∗ x . v ˆ ( n − ∗ x . d1 , ( n ∗ ( n − ∗ x . v ˆ ( n − ∗ ( x . d1 ∗ x . d1 ' ) + ( n ∗ x . v ˆ ( n − ∗ x . d2 )Listing 2: Operators overloading27n Listing 3 the function given by equation (31) is defined in the first line, in the following line it isevaluated for the double precision floating point values x1=x2=x3=1.0 and the result is printed. y ( x1 , x2 , x3 ) = x1 ˆ3 ∗ x2 ˆ2 + x3 ˆ2 p r i n t l n ( ” \ n y0 : ” , y ( 1 . , 1 . , 1 . ) ) y0 : 2 . 0 Listing 3: Numerical example with Float64 argumentsIn Listing 4, the variables x1 , x2 and x3 are defined as dual quantities, of the D2 type, where thefirst argument of the call to the D2 constructor is the value of the variable, the second argument is thegradient of each variable, and the third argument is the Hessian. We remark that independent variablesare defined by properly specifying the components of their gradient and Hessian. In fact, independentvariables are such that the derivative with respect to themselves is one, while all other derivatives arenought. Thus x1 is the independent variable that occupies the first position of the gradient, x2 thesecond and x3 the third. As a result, the derivative with respect to x1 of any operation involving x1 , x2 and x3 , with will be stored in the first component of the d1 array of the result, and similarly forthe derivatives with respect to x2 and x3 , and for higher order derivatives. Following in Listing 4, thesame function y, defined in Listing 3, is called with the dual quantities just defined and the result,which is a dual quantity itself, is printed. As we can observe, the result returned by the function thistime includes both the value of the function y(x1,x2,x3) , and its gradient and Hessian. x1 = D2 ( 1 . , [ 1 . , 0 , 0 ] , z e r o s ( 3 , 3 ) ) x2 = D2 ( 1 . , [ 0 , 1 . , 0 ] , z e r o s ( 3 , 3 ) ) x3 = D2 ( 1 . , [ 0 , 0 , 1 . ] , z e r o s ( 3 , 3 ) ) p r i n t l n ( ” \ n yd : ” , y ( x1 , x2 , x3 ) ) yd : D2 ( 2 . 0 , [ 3 . 0 , 2 . 0 , 2 . 0 ] , [ 6 . 0 6 . 0 0 . 0 ; D2 argumentsWe observe that having overloaded the operators involved in the definition of y(x1,x2,x3) allowed usto call the same function with both data type without making any modification or having to add anyspecification to the function itself.We remark that the implementation of dual numbers in Julia as presented in this section is anattempt to provide a brief and clear illustration of a possible computer implementation of AD, throughoperators overloading, nonetheless in this form it does not exploit any of the powerful features offeredby the Julia programming language, like parametric types, and macros [Perkel, 2019, Bezanson et al.,2017]. The implementation developed for the solution of the example presented in the paper, availablethrough [Vigliotti, 2020], which is based on Revels et al. [2016], makes a better use of Julia’s featuresand functionalities and ensures better performances than the example presented in this section.28 Arbitrary order dual number systems
In this section we briefly generalize the definition of dual numbers to an arbitrary order of differentia-tion. Let x be a dual number of dimension N and order K x ≡ x + x i ı i + x i i ı i i + x i i i ı i i i + · · · + x i ...i K ı i ...i K with i ∈ . . . Ni ∈ i . . . Ni ∈ i . . . N ... i K ∈ i K − . . . N (69)with ı j the canonical base of R N , with j ∈ . . . N , and ı i , ı ij , ı ijk . . . , ı i ...i K are symbols defined as ı i i ≡ ı i ⊗ ı i + ı i ⊗ ı i ı i i i ≡ ı i ⊗ ı i ⊗ ı i + ı i ⊗ ı i ⊗ ı i + ı i ⊗ ı i ⊗ ı i + ı i ⊗ ı i ⊗ ı i + ı i ⊗ ı i ⊗ ı i + ı i ⊗ ı i ⊗ ı i ... ... ı i ...i K ≡ (cid:88) I K ∈ Π( K ) ı I K ⊗ · · · ⊗ ı I KK (70)where I K is a permutation of the indices 1 . . . K , I Ki are its elements, and Π( K ) is the set of all thepermutations of the indices 1 . . . K . With respect to equations (70) we observe that the following holds ı ij = ı ji ı ijk = ı ikj = ı jik = ı jki = ı kij = ı kji ... ... ı I K = ı J K ∀ I K , J K ∈ Π( K ) . (71)The quantities x , x i , x ij , x ijk . . . , x i ...i K are real scalars and are the components of x , x is the realpart of x , the remaining are dual parts of order 1 , , . . . K . Two dual numbers of dimension N andorder K are identical if and only if all of their components are identical, as follows x = y ⇐⇒ x = y x i = y i ... ... x i ...i K = y i ...i K (72)The sum of two dual numbers is defined as the dual number whose components are the sum of thecomponents, as follows z = x + y ⇐⇒ z = x + y z i = x i + y i ... ... z i ...i K = x i ...i K + y i ...i K (73)29he product of two dual numbers is a dual number obtained as the sum of the mixed products of theircomponents, where the following rules applies for the product of the symbols ı i , ı ij , ı ijk . . . , ı i ...i K ı i ı j ≡ ı ij ı i ı j ı k = ı i ı jk ≡ ı ijk ... ... ı . . . ı K = ı ı ...K ≡ ı ...K , (74) ı i ı ...K ≡ , (75)where equations (74) produce the contribution to higher terms in the product as results of the productsof lower order terms in the factors, and equation (75) ensures that no component with order higherthan K appears in the result. The components of the product are given as z = xy ⇐⇒ z = x y z i = x i y + x y i z ij = x ij y + x i y j + x j y i + x y ij z ijk = x ijk y + x ij y k + x i y jk + x y ijk ... ... z i ...i K = x i ...i K y + x i ...i K − y i K + · · · + x y i ...i K , (76)With reference to the quotient of two dual numbers, we observe that this operation is equivalent tothe product of the first time the inverse of the second, where the inverse of a dual number is obtainedby solving the following 1 x = y ⇐⇒ yx = 1 , (77)from which it results1 x = y ⇐⇒ y x = 1 x i y + x y i = 0 x ij y + x i y j + x j y i + x y ij = 0 x ijk y + x ij y k + x i y jk + x y ijk = 0... ... x i ...i K y + x i ...i K − y i K + · · · + x y i ...i K = 0 , (78)where we observe that the right hand side of the definition (78) is an lower diagonal linear system inthe unknowns y ······