Bayesian spatiotemporal inference of trace gas emissions using an integrated nested Laplacian approximation and Gaussian Markov random fields. Geoscientific

The increasing reliance of aerospace structures on numerical analyses encourages the development of accurate, yet computationally eﬃcient, models. Finite element (FE) beam models have, in particular, become widely-used approximations during preliminary design stages and to investigate novel concepts, e.g. aeroelastic tailoring. Over the last 50 years, developments in hp-FE methods based on elements of variable size (h) and polynomial degree (p) have helped reduce the computational cost of numerical analyses. Concurrently, many structures, including aircraft wings and wind turbine blades, have gradually increased in length, slenderness, and complexity. As a result, modern blades and wings regularly op-erate beyond the range of linear deformation, requiring non-linear analyses for which hp-FE methods are often not readily applicable. The aim of this paper is, therefore, to derive a co-rotational ﬁnite element formulation for enriched 3-, 4- and 5-noded beam elements, suitable for non-linear hp-FE reﬁnement. To this end, we derive the mathematical formulation to incorporate enriched elements within the co-rotational FE beam framework. The proposed formulation is then validated against multiple non-linear benchmark problems and an experimental case study.


I. Introduction
Modern aircraft wings and wind turbine blades, driven by competition for increasing performance, are more slender and flexible than their predecessors. The ability of numerical models to accurately predict the flexibility of aerospace structures has therefore become central to further improve aeroelastic performance, e.g. aeroelastic tailoring [1][2][3]. Furthermore, considering the reliance on computer software and the critical choices that must be made during early design phases, the development of accurate, yet rapid, structural analysis methods is of paramount importance for aerospace applications [4]. Finite element (FE) beam models, which are known to be reliable at predicting global structural behaviours with relatively few degrees of freedom (DoFs), have consequently become widely used approximations for the preliminary design and aeroelastic tailoring of aerospace structures [5,6].
Over recent decades, prompted by successful numerical and experimental applications [7][8][9][10], the aeroelastic tailoring of variable stiffness structures has gained significant academic and industrial interests. Designing structures to purposefully withstand and avail of non-linearities has also been shown to provide a means to further improve the performance of aerospace structures [11,12]. Although both non-linear and variable stiffness designs hold promising potential for the field of aerospace, the development of accurate and rapid analysis design and optimisation methods for non-linear structures with spatially variable properties remains a challenge.
The work on FE beam modelling presented in this paper is aimed at the non-linear analysis of slender, beam-like, variable stiffness structures, with primary focus on the use of refined elements and computational efficiency. In a first attempt to address this challenge, within a linear framework, the authors developed a procedure for the automated generation of 3-, 4-and 5-noded linear beam elements as well as a spanwise integration method [13]. The authors showed that these enriched elements could significantly reduce spurious strains and improve the computational performance of linear FE beam models (through a two-to threefold reduction in DoFs). The aim of the present paper is to extend the co-rotational finite element beam formulation to enable the use 3-, 4-and 5-noded elements for non-linear structural analyses. The expected advantages include, convergence with fewer D.O.Fs, reduced spurious strains, and most importantly the capability to transfer the advantages of many dedicated linear elements to the analyses of non-linear problems.
Although other non-linear beam formulations are available (e.g. total Lagrangian and updated Lagrangian), our objective is to provide designers with the possibility to utilise existing enriched linear elements for nonlinear analyses. The co-rotational method solves non-linear structural problems by splitting the deformation of beam elements into rigid body motions and local deformations. In contrast to linear FE which assumes small rotations, large rotations are captured, in the co-rotational approach, by rigid body rotation matrices. Linear beam elements can, therefore, be used as long as strains remain small within each element. As such, the co-rotational framework and its ability to employ linear beam elements is conveniently aligned with our objectives.
The literature on the topic include many co-rotational beam element frameworks, starting with the pioneering work by Wempner, Oran and Belytschko [14][15][16], dating back to the '60s and '70s. Following their work, in the '90s, Crisfield et al. [17] developed a unified co-rotational framework for solids, shells and beams. More recently, Li et al. [18] proposed a mixed co-rotational 3-noded beam element employing vector rotations as variables, but assuming a bi-symmetric cross-section. Interesting developments of the co-rotational formulation have also been proposed for shells, in particular with regards to computational efficiency [19,20] and laminated structures [21]. Herein, the 2-noded co-rotational beam element framework proposed in 1991 by Nour-Omid and Rankin [22], and further developed by Pacoste and Eriksson [23] and Battini [24], is used as a starting point for the extension to 3-, 4-and 5-noded elements. First, we propose undeformed and deformed element configurations from which we derive the mathematical relationship needed to split global and local displacements, through the evaluation of rigid body motion matrices. Second, we equate global and local work in order to obtain the tangent stiffness matrix. The extended framework is applied and validated against multiple benchmark case studies and experimental data. Note that, due to extensive amount of derivations, the present manuscript focuses primarily on the mathematical formulation required to enable the use of enriched linear beam elements within the co-rotational framework. The authors demonstrated the benefits arising from the use of enriched linear elements for linear analyses in [13]. From first principles of structural analysis and numerics, it is possible to deduce that these benefits will translate to the non-linear regime, possibly with even more pronounced effects.
The rest of this paper is structured as follows. The co-rotational framework is briefly put into the context of non-linear analysis in Section II, followed by the analytical derivations of the 3-, 4-and 5-noded beam elements. The validation of the proposed elements using benchmark case studies is carried out in Section III, while concluding remarks are provided in Section IV. Prolonged mathematical derivations have been moved to the appendices for the sake of brevity, clarity, and completeness.

II. A Co-rotational Framework for 3-, 4-and 5-Noded beam elements
Throughout the rest of this paper, vectors are shown in bold or enclosed in round brackets, e.g. d or (x, y, z) , while matrices are enclosed in square brackets, e.g. [K]. The upper script • denotes the transpose operator.
The co-rotational method is well suited to the analysis of geometrically non-linear elastic problems, which are, under static equilibrium, expressed as where R denotes a residual, F int (x) and F ext are internal and external forces respectively, and x is the vector of structural DoFs. Typically, the non-linear deformation of structures subject to monotonically increasing loads is found by using the Newton-Raphson method, iteratively solving for the change in displacements between the current and the next iteration ∆x i+1 , where [K t (x)] is the tangent stiffness matrix.
The co-rotational beam formulation, based upon the combined use of linear elements and local reference frames, provides a simple means to evaluate [K t (x)] and F int (x). Because the co-rotational formulation heavily relies on rotation and transformation matrices between reference frames, the definition of coordinates system is introduced first, in the following section. Background knowledge, including a brief introduction to rotation and spin vectors, is provided in Appendix A. The keen reader is referred to [25] for a detailed introduction to the co-rotational formulation.

A. Reference Frames and Rigid-Body Motions
Multiple coordinate systems are necessary to describe the rigid motions and structural deformations of beam elements. Without loss of generality, the global inertial reference frame, the triad [T g ], is taken as the identity matrix where E 1 , E 2 and E 3 are the global reference frame triad vectors. For simplicity in the following derivations, the undeformed elements are assumed to be straight with uniformly spaced nodes, as illustrated in Figure 1.
The initial, i.e. undeformed, orientation of an element is described by the orthonormal basis where E 10 points along the element axis. Fig. 1 Undeformed elements, L0 ij is the initial length between node i and node j The non-linear deformation of an element is described by a global rigid body motion, and a linear deformation expressed in a local reference frame that rigidly rotates and translates with the element. The element local reference frame is represented by the triad attached to one of the element nodes, acting as the local origin (0, 0, 0) , and where the subscript • r stands for rigid. The influence of the local origin node on the computational efficiency of the formulation is not investigated in the present work, but the reader is referred to the work of de Veubeke [26] for further details. The normalised vector r 1 is defined by the line connecting the local origin node (0, 0, 0) to the following node as illustrated in Figure 2. The remaining vectors r 2 and r 3 are calculated based on r 1 to obtain an orthonormal basis for [T r ], details of which are provided in Appendix B. Nodal rotations are represented by the structural node triads where k = 1, ..., N refers to the triad node, and i = 1, 2, 3 refers to the local triad axes. Last, the length between two nodes i and j is denoted L ij . The rotation matrices depicted in Figure 3 provide a means to transform vectors between coordinate systems. The first matrix, [R 0 ], rotates the global reference triad [T g ] to the undeformed triad [T 0 ] Next, the [R r ] rotation matrix Finally, the element triad [T r ] is rotated to the nodal triads with the rotation matrices R i , expressed in local reference frame, as

Fig. 3 3-noded beam element reference frames
The proposed notation is also applied to the 4-and 5-noded elements, as illustrated in Figure 4. The global and local nodal displacement vectors of an N -noded beam element are respectively denoted d g andd, and defined as where u ig and θ ig are the nodal translation and rotation vectors of node i, expressed in the global reference frame, and where the • g subscript stands for global. Similarly to their global counterparts,ū i andθ i are the nodal translation and rotation vectors of node i, but expressed in the local element reference frame. Throughout this manuscript, the bar symbol,•, is used to denote values expressed in the local element reference frame. Following this notation, the global and local force vectors are in which F and M are the nodal forces and moments.
A key step during the application of the co-rotational framework involves the extraction of the rigid body motions and local deformations from the global displacement vector d g . To this end, we look for a matrix [B] between variations in global and local DoFs such that Because three dimensional rotations are non-linear, [B] is a local approximation and must be updated as the structure deforms. The internal force vector F g , expressed in the global reference frame, is found by equating local and global works due to these variations and substituting Eq. (15) into (16) to obtain a direct relationship between global and local forces The tangent stiffness matrix is then found by taking the variation of Eq. (17) where [K] is the local element stiffness matrix satisfying δF = [K]δd. By identification, we find the symmetric global tangent stiffness matrix [K g ] to be The tangent stiffness matrix is generally not symmetric due to the first term on the right hand side of Eq. (18). However, if symmetry is recovered at equilibrium then truncating the tangent stiffness matrix to its symmetric part should not significantly impede the convergence rate of non-linear solvers [17]. For completeness, both the symmetric and non-symmetric tangent stiffness matrices are derived.
The following subsections are dedicated to the derivation of the transformation matrix which is at the core of the co-rotational method, i.e. [B]. The procedure used to evaluate [B] is divided into five steps.
(1) In Section C, the local displacement vectord is extracted from the global displacement vector d ḡ We then progressively expand the algebraic expression for δd to express it as a function of δd g and concurrently evaluate the global tangent stiffness matrix [K g ], and the force vector F g .
(2) In Section D, rotation vectors are replaced by their equivalent spins δθ → δw. The spin variables are local approximations of rotations, simplifying intermediate derivations to ultimately find the relationship between local rotations δθ and the global displacement vector δd g .
The intermediate displacement vector δd a is defined as (3) Variations of local translations δū and spins δw are expanded as functions of global translations δu g and spins δw g , in Section E, to find [B g ] such that (4) Global spins are transformed back into rotations δw g → δθ g , completing the procedure, C. From Global to Local Displacements -Step (1) In this section, local deformations are recovered from the global displacement vector.

Nodal Translations
We start with the 3-noded beam element translations of nodes 1, 2, and 3. The origin (0, 0, 0) of the local reference frame is set to correspond to the deformed position of node 1 as shown in Figure 2. Hence, the translation of node 1 in the local reference frame are always nullū The element triad vector r 1 connects nodes 1 and 2 with a straight line. The translation of node 2 expressed in the local element reference frame is thereforē where L 12 and L0 12 refer, respectively, to the deformed and undeformed lengths between node 1 and node 2. In other words, the global deformation of node 2 simplifies to a translation along r 1 in the local element reference frame. The local displacement of node 3,ū 3 , as illustrated in Figure 2, is obtained by expressing the difference between the rigid and deformed positions of node 3 in the local reference framē Following similar derivations for the 4-and 5-noded elements, the nodal translations expressed in the local element reference frame are summarised in Table 1.

Nodal Rotation Vectors
We continue with the recovery of local rotations from the global displacement vector. Recall, with the help of Figure 3, that [R i ] is a local rotation matrix which rotates the local element triad [T r ] onto the node triad T k i . The [R i ] matrix is evaluated as Local rotation vectors are then retrieved from the rotation matrices employing the matrix logarithm The local displacement vector is obtained by combining translations and rotations from Table 1 and Eq. (26), which results ind for the 3-noded beam element.
In the following sections, the algebraic expression for the virtual displacement δd is expanded to be expressed as a function of δd g in order to evaluate the tangent stiffness matrix [K g ].
D. From Local Rotation to Spin -Step (2) Variations of local rotational DoFs in δd are converted into their equivalent spins δw in this section. Using spin conveniently provides a first order approximation for the variation of rotation matrices The matrix [B a ] links variations of local rotations δθ i to their equivalent spins δw i such that with in which local variations of rotations and spins are linked via the T s −1 function, defined in appendix Eq. (A. 15), The non-symmetric tangent stiffness matrix, K a ns such that δF a = K a ns δd a , is found by taking the variation ofF and substituting δF = [K]δd to obtain The symmetric stiffness matrix is identified to be The non-symmetric tangent stiffness matrix (subscript "ns") is with δ[B a ] and δ T s −1 (θ i ) , given in appendix Eq. (A.17).

E. From Local to Global Displacements -Step (3)
In this section, the displacement δd a is transformed into its global counterpart δd a → δd g,spin = (δu 1g , δw 1g , . . . , δu 3g , δw 3g ) , by virtue of the matrix [B g ] such that The derivation of [B g ] involves lengthy calculations which are briefly summarised in the following subsections. Detailed mathematical derivations are provided in Appendix C.

Nodal Translations
We start with local translations and global displacements The variation of nodal translations introduced in Table 1 are summarised in Table 2. Table 2 Variations of nodal translations The individual values of [C i ] are retrieved by substituting δL ij and [δ [R r ] ] in the corresponding rows of Table 2. The variation of δū 2 for the 3-noded element is, for instance, found to be where for clarity, and only in this equation, 0 is a 1 × 3 row vector.

Nodal Spin Vectors
Similarly, the relationship between variations of local spins and global displacements is for which the variation of the rotation matrix with,

Tangent Stiffness Matrix [K g,spin ]
The global tangent stiffness matrix is defined such that δF g,spin = [K g,spin ]δd g,spin , where according to Eqs. (38) and (40) we have Taking the variation, the symmetric tangent stiffness matrix is identified as and the tangent stiffness matrix is For simplicity, we suggest to evaluate ∂[B g ] ∂d g,spin numerically using finite difference.
This step is straightforward and involves, once more, the relationship between spin and angle variables The tangent symmetric stiffness matrix is then given as and the non-symmetric stiffness matrix is The derivation of δ ([T s (θ ig )]) is provided in Appendix (A.16).

G. Procedure Summary -Step (5)
Combining the transformations introduced in previous steps, we obtain Finally, the symmetric tangent stiffness matrix is whilst the non-symmetric tangent stiffness matrix is

III. Applications
In this section, the derived co-rotational formulation for 3-, 4-, and 5-noded beam elements is validated against increasingly difficult benchmarks found in the literature, and one empirical wind turbine static blade test. Timoshenko's beam element formulation, as previously extended by the authors in [13], is employed in this study. The elements shape functions for 2, 3, and 4 nodes are shown in Figure 5 for convenience.

A. Non-linear Bending of a Slender Beam
Our first application is the bending analysis of a slender cantilever beam subject to a tip moment, illustrated in Figure 6. This example is used to verify that the 2-and 3-noded elements have correctly been derived and implemented. The material Young's modulus, Poisson's ratio and shear factor are respectively denoted E, ν and κ. The beam is 100 meters long. Figure 7 presents the results compared with data from Li and Vu-Quoc [18]. It shows that the non-linear displacements obtained with ten uniformly spaced 3-noded beam elements agree well and converge towards the literature benchmark.  The second benchmark example is the clamped 45 • curved rectangular beam subjected to a bi-axial concentrated tip load P proposed by Bathe and Bolourchi [27], reproduced in Figure 8. In this example, we compare the predictions of 4-and 5-noded elements with those of 2-and 3-noded elements. The bend has a radius R of 100 in, a cross-sectional area of 1 in 2 , and lies in the X-Y plane. Nodal tip deflections obtained with 10 uniformly spaced elements are presented in Figure 9 as functions of the non-dimensionalised load parameter k = PR 2 EI . The structural analysis is first carried out with conventional 2-noded elements, results are shown in Figure 9a. As this figure shows, the conventional co-rotational approach employing 2-noded linear Timoshenko beam elements is subject to locking and consequently shows significant discrepancies with the expected results. By contrast, results obtained for the enriched 3-, 4-and 5-noded elements, with higher order shape functions, show close convergence to the results reported in [27]. This benchmark demonstrates that the refined elements are free of locking and can readily be used within the proposed extended co-rotational beam element formulation. The third example, depicted in Figure 10, and taken from Mathisen et al. [28], is a typical benchmark problem for non-linear beam formulations under combined bending, shear and torsion for non-smooth threedimensional geometries. It consists of three straight beams of equal length connected at right angles. The structure is clamped at one end and free at the other, where two point loads are applied. The nodal tip displacements (i.e. U x , U y , U z ) and global rotation vector values (i.e. R x , R y , R z ) predicted with eighteen 3-noded elements, and twelve 4-noded elements are compared against the literature in Figure 11, which shows excellent agreement. Although results for 5-noded elements are not presented for the sake of brevity, they are also found to agree with the reported data. Fig. 10 Three-legged beam, reproduced from Mathisen et al. [28], R x , R y and R z illustrate the convention used for positive rotations.

D. Deployable Ring
The deployable ring example illustrated in Figure 13, introduced by Goto [29] and reused in Battini's work [24], is our fourth benchmark. The ring is clamped at one end and twisted at the other through the application of a moment M at point A. The ring is also constrained at point A such that only translation and rotation with respect to the X axis are allowed. This example is ideal for benchmarking the proposed formulation under large rotations. As the rotation at point A increases, the ring slowly wraps around itself and transforms into a smaller ring with one-third of its original radius. As the rotation continues further, the ring is brought back to its initial configuration. The results are compared against the literature in Figure 15 and the half ring model deformation is show in Figure 14.

E. Right-Angle Frame
The right-angle frame example illustrated in Figure 16 and found in [30] is used as our fifth and last academic benchmark case study. It is sufficient to model half of the frame, employing the symmetry condition about the y-z plane through the apex, to capture the frame non-linear deformations. Translation along the x axis and rotation about the z axis are unconstrained at the support, whilst a point load moment about the z axis is applied. The cross section has a thickness/height ratio of 1/50. The material properties are given as E = 71240N/mm 2 for the young's modulus and ν = 0.31 for Poisson's ratio. As the load increases towards the first bifurcation point, at about 620 N/mm the frame buckles. A small geometrical imperfection, in the form of an offset in the z-direction, is introduced at the support to continue past the otherwise singular tangent stiffness matrix. As the load continues to increase into the post-buckling regime, the frame deforms out-of-plane, rotating about the x-axis until a full revolution brings it back to its initial positions under the opposite buckling load as shown in Figure 17, in blue. At this stage, the solver stops as it encounters a second singularity corresponding to a higher buckling mode, incompatible with the original geometric imperfection, as explained in [30]. Note that the 3-Noded beam element results are solely presented for the sake of brevity, whilst the 2-Noded linear beam element failed to capture the structure non-linear response.

F. Wind Turbine Blades
Last, but not least, the proposed co-rotational formulation is applied to a real-life case study: The static bending test of a modern wind turbine blade, which has been observed to exhibit small yet noticeable geometrical non-linearities under extreme loading conditions. In comparison to the prismatic structures used in the previous academic benchmarks, the experiment showcases the framework's ability to model realistic and complex non-prismatic beam-like structures. The test set-up and results, presented in Figures 18  and 19, have been normalised due to proprietary agreements. The blade is clamped at its root whilst six concentrated loads distributed along the blade span are applied to loading frames through cables linked to hydraulic cylinders, on the floor. Two load cases are evaluated. First, the blade is bent to achieve compression on the flapwise suction side. Second, the blade is rotated 180 • and bent under a different load distribution and magnitude to test the suction side under tension (i.e. positive and negative flapwise root bending moments).
In spite of the model complexity, the numerical deflections, predicted by the co-rotational formulation based on 80 3-Noded elements, show good agreement with the experimental measurements. Note, the number of elements is much higher in this case study, in comparison to previous ones, as it is necessary to accurately capture the span-wise variations of blade properties.
Bla de Ax is Pulleys Fig. 18 Wind turbine blade bending test.

IV. Conclusion and Future Work
The mathematical derivation needed to implement enriched 3-, 4-and 5-Noded linear beam elements into the co-rotational framework for non-linear beam analyses has been presented in this paper. Although the methodology needed to derive these elements was known in theory, the detailed derivation required for implementation had not yet been reported in the literature. Herein, the lengthy derivations are described in full for the sake of reproducibility.
The proposed formulation has been successfully validated against five non-linear benchmark problems and an experimental wind turbine blade static test. The application of enriched elements within non-linear structural analyses is shown to resolve the classical issue of locking experienced by conventional Timoshenko beam elements and agree well with the reference data for all benchmarks. The proposed framework provides a simple means to transfer the advantages of enriched linear beam elements to non-linear analyses. The expected advantages from hp refinements in a non-linear setting include a reduction of the DOFs, and hence the computational cost, needed for convergence-the reason being the gain in accuracy due to the use of higher order shape functions.
The use of enriched linear beam elements in non-linear co-rotational formulations, coupled with hierarchical theories such as the Unified Formulation by Carrera and co-workers [31], could lead to accurate, yet inexpensive, descriptions of 3D strain and stress fields in complex beam structures. Ongoing and future work by these authors include, for instance, the application of the analyses method presented in this paper to the design of geometrically and materially bend-twist-coupled wind turbines blades, focusing on computationally efficient aeroelastic optimisation. [31] E. Carrera and G. Giunta, "Refined beam theories based on a unified formulation," International Journal of Applied Mechanics, vol. 2, no. 01, pp. 117-143, 2010.

Appendix A. Rotation and Spin Vector
This section serves as a short introduction to rotation and spin vectors.

Axis-angle Representation
An axis-angle representation is used to describe 3D rotations as vectors in whichn denotes the normalised rotation axis and θ is the angle of rotation, as illustrated in Figure A.1. Note that θ 1 , θ 2 and θ 3 in Equation. (A.1) are generally not equal to the respective rotations around the E 1 , E 2 and E 3 axes. The rotation matrix corresponding to the θ vector is obtained employing Rodrigues' formula where [θ] is the skew matrix of θ defined as

Spin
Consider the 2D example illustrated in Figure A.2, in which the vectors P i , P 0 , and P n respectively denote a randomly orientated initial vector, its image after being rotated by θ 0 , and its image after by being further rotated by an infinitesimally small rotation δw. One can define P n in three equivalent ways Note the key difference between the non-additive rotation (or spin) δw and the additive rotation δθ. We can further rewrite the right hand side of Eq. (A.4) as where in 2D δP is equal to a vector of length δw||P 0 || pointing in the direction normal to P 0 as shown in Figure A.2 and expressed as δP = δw||P 0 ||n = δw||P 0 || − sin(θ 0 ), cos(θ 0 ), 0 , (A. 6) or, where where T s and its inverse are defined as It should be noted that Eq. (A.14) ceases to be bijective for θ = 2kπ with k = 0, 1, ..., n and that care must be taken during its numerical implementation.
The variation of (T s (θ)) is (A. 16) and the variation of its inverse T s −1 (θ) is (A. 17) In addition, we have and, In this section the calculation required in order to compute the element and node triads, [T r ] and T k i respectively, are introduced for the 3-Noded beam element presented in Figure 2.

Node Triads
The node triads are defined in the global reference frame as which are calculated based on global rotations such that The global rotation matrix [R ig ] connecting the un-deformed and deformed triads of node i is given using the Rodrigues' formula where θ ig is the global rotation vector associated with node i, whilst the angle of rotation and the tilde denotes the skew matrix.

Element Reference Frame
The where P 1g and P 2g are the global position vectors of node 1 and 2, respectively. Furthermore, the global position vectors are and with u 1g and u 2g the global translation vectors of nodes 1 and 2.

C.1. Element Length
Starting with the variation of element lengths δL ij , the 3-Noded element length L 12 defined as the Euclidean distance between the deformed positions of node 1 and 2 is where P 20g and P 10g denote the undeformed position vectors of node 1 and 2, expressed in the global reference frame. Taking the variation of Eq. C.1 we obtain and Employing similar derivations for the 4, and 5-noded elements, the results are summarised in Table 3.
The expressions presented in Table 3 can be substituted in Table 2 to obtain the direct relationship between variations in local translations δū and their global counterparts δu g .

C.2. Element Directional Vector
We continue with the derivation of the 3-noded element vector r 1 starting as expanding into and substituting the value of δL 12 from Table 3 δr where [I] is the 3x3 identity matrix. Similar expressions for the 4, and 5-noded elements are summarised in Table 4.
and isolating the spin variables we find The values of [G] are found at the end of this paper in Table 5, and derived in the following appendices.
Once δw re is found, its equivalent in global reference frame is obtained by transforming it back using [R r ] as where S(•) denote the skew matrix.

C.4. Nodal Rotation and Spin
The variation of spin vectors are detailed in this section. We are looking for an expression of the form of Equation (40) recalled here δw i = [D i ]δd g,spin , (C. 16) which relates the variation between the local spin vector at node i to the global displacement vector. We start by noting that as demonstrated in the optional proof below. Moreover, we already know that δw re = [G] [R r ] diag δd g,spin . (C.18) After employing the chain rule on δw ie we obtain δw ie = ∂w ie ∂d g,spin δd g,spin = ∂w ie ∂d e,spin ∂d e,spin ∂d g,spin δd g,spin , (C. 19) in which we recognise that ∂d e,spin /∂d g,spin corresponds to a change of coordinate between the global system and local  Recall the variation of the vector attached to the element reference frame δr 1 for the 3-noded element from Table 4 to be Finally, the evaluation of δw re1 , slightly more complicated, can be found in [32] and results in δw re1 = −r 2 δ r 1 × q ||r 1 × q|| (C.39)