Geometrically exact elastoplastic rods: determination of yield surface in terms of stress resultants

This work addresses the determination of yield surfaces for geometrically exact elastoplastic rods. Use is made of a formulation where the rod is subjected to an uniform strain field along its arc length, thereby reducing the elastoplastic problem of the full rod to just its cross-section. By integrating the plastic work and the stresses over the rod’s cross-section, one then obtains discrete points of the yield surface in terms of stress resultants. Eventually, Lamé curves in their most general form are fitted to the discrete points by an appropriate optimisation method. The resulting continuous yield surfaces are examined for their scalability with respect to cross-section dimensions and also compared with existing analytical forms of yield surfaces.


Introduction
The theory of geometrically exact rods and its implementation have been applied in a wide range of research fields to model slender structures whose dimension in one direction is dominant compared to the other two dimensions. In contrast to classical beam theories such as the Euler-Bernoulli and Timoshenko beam theories, the theory of geometrically exact rods is not restricted to the geometrically linear case, but is able to describe large deformation as well as large rotation of the rod's cross-section. This is the reason that rods are often used to model long nano wires and tubes [1,2], cables [3], hair and other biological fibres [4]. The concept of modelling large rotation in slender structures goes back to the Cosserat brothers [5]. They consider every material point of the rod to be connected to the centerline and a triad of orthonormal vectors characterizing the orientation of the rod's cross-section [6].
The modelling of three-dimensional elastoplasticity, both for the geometrically linear and nonlinear case, has been widely explored in the literature. As a key reference, the B Ludwig Herrnböck ludwig.herrnboeck@fau.de 1 groundbreaking work of Simo and Hughes [7] can be mentioned. Attempts to reduce these theories to slender structures such as beams and rods have been made. For the case of linear beam theory, elastoplastic deformations have been taken into account in [8]. For the case of geometrically nonlinear beams, we refer to [9][10][11][12]. Attempts towards a direct approach to model viscoplastic and thermoelastoplastic rods are, e.g., [13] and [14,15], respectively. There, the assumption is made that the strain prescriptors (the rod's translational strains and curvatures) split additively into elastic and plastic parts. All attempts to model elastoplastic rods require profound knowledge of an appropriate yield criterion defined by a yield surface. Obtaining the yield surface in terms of stress resultants (the rod's internal (contact) forces and internal (contact) moments) for rods of arbitrary cross-sectional shape and with arbitrary constitutive response still remains a challenge. Yield surfaces have been derived for different cross-sections in [10,[16][17][18][19], applying the upper and lower bound techniques of limit analysis. There, the rod's crosssection at yield is assumed to be fully plastified at the state of yield limit. These yield surfaces are often used to investigate inelastic properties and limit loads of steel frames [20,21]. The objective of this work is to present a framework enabling to determine a yield surface in terms of stress resultants for an arbitrary cross-section and deformation state. Especially, the rod's cross-section need not be fully plastified at yield: the yield limit is defined in terms of plastic work done. This distinguishes the yield surface obtained here from the above Table 1 In the sequel, strain measures and their energetically conjugates will be denoted differently for geometrically exact rods and for threedimensional continuum mechanics, respectively Geometrically exact rods Three-dimensional continuum mechanics

Strain measures Strain prescriptors Strains
Energetically conjugates

Stress resultants Stresses
Strain prescriptors collectively denote the rod's translational strain and curvature, whereas stress resultants collectively denote the rod's internal (contact) forces and internal (contact) moments mentioned ones and ensures versatility in the definition of the yield limit. In practice, the determination of such yield surface is made a priori for a specific cross-section and can then be incorporated into an elastoplastic rod model. The paper is structured as follows. In Sect. 2, a brief review of the theory of geometrically exact rods with rigid cross-sections is given and extended to the inelastic case. For uniformly deformed rods, the deformation of a cross-section can be evaluated, exclusively depending on strain prescriptors. An appropriate model is outlined in Sect. 3 for both the elastic and the inelastic case. Forces and moments as well as the rod's stiffnesses are evaluated from the deformed (warped) cross-section and the elastoplastic response of geometrically exact rods is investigated. The model presented in this section is used to evaluate the yield surfaces later on.
In Sect. 4, a framework is presented aiming to obtain yield surfaces in terms of stress resultants. Therefore, discrete points on yield surfaces resulting from numerical simulations are computed. Generalised Lamé curves are fitted to the discrete points on the yield surfaces to obtain a continuous yield surface in terms of stress resultants. Finally, the resulting yield surfaces are comprehensively discussed, examined regarding scalability of the cross-section dimension and compared with existing ones from the literature. Table 1 summarises the different designation of strain measures and their energetically conjugates for geometrically exact rods and for three-dimensional continuum mechanics, respectively. This nomenclature will be used in the sequel.
The origin of the terms strain prescriptors and stress resultants is explained in chapter 3. Throughout the paper, the arc length of the rod in its undeformed configuration B 0 is denoted by s. The derivative of an arbitrary quantity (•) with respect to s is defined by All numerical simulations carried out to achive the results shown in this contribution base on the open source finiteelement library deal.II [22]. Illustration of a geometrically exact rod deforming from its straight material configuration B 0 into its deformed spatial configuration B t . Any point in the material configuration is denoted by X. The position of the centerline of the deformed rod is denoted by r. Material points in the spatial configuration are denoted by x

Geometrically exact rods with rigid cross-sections
In this section the theory of geometrically exact rods with rigid cross-sections is briefly reviewed [23,24]. First we present the kinematics, then we introduce stress resultants both for the case of elastic rods and elastoplastic rods. Figure 1 shows a typical rod deforming from its undeformed straight material configuration B 0 into its deformed spatial configuration B t .

Kinematics
Here, X = X i E i with X 3 = s denotes the position of any material point in the undeformed material configuration. In the sequel, Latin indices run from 1 to 3, whereas Greek indices run from 1 to 2. The coordinates of the undeformed cross-section 0 are denoted by X α . The position of the deformed centerline is denoted by r(s). The orthogonal directors d 1 (s), d 2 (s) span the plane of the deformed cross-section t , assumed planar. Together with d 3 (s) = d 1 (s) × d 2 (s), an orthogonal basis is constructed. It is not a requirement that d 3 (s) is tangential to the deformed centerline r(s). Taken together, the directors d i describe the orientation of the deformed cross-section. In the material configuration the orientation of the cross-section is described by the directors D i (s), where D i (s) = E i . Usually the directors D i (s) are oriented along the principal directions of the undeformed cross-section 0 . The directors in the spatial configuration are related to the directors in the material configuration via Here R(s) denotes a three-dimensional rotation tensor and is a member of the special orthogonal group SO (3). The quantities r(s) and R(s) are the kinematic variables of this theory. Assuming a rigid cross-section, the position x of any material point in the spatial configuration can be described as witĥ For the sake of readability, the dependencies on s and X α are omitted from now on. In this restricted setting, the crosssection of the deformed rod remains rigid. A framework to model the cross-sectional deformation depending on strain measures only is outlined in Sect. 3, based on [25,26]. Two strain measures for the rod are introduced: the translational strain v, and the rotational strain k (also called curvature): Their rotational pull backs are defined as v 0 and k 0 in the following way: Here, K and K 0 are skew-symmetric tensors whose axial vectors are k and k 0 , respectively. We collectively denote v 0 and k 0 as strain prescriptors. Incorporating the strain prescriptors, the deformation gradient of map (2) can be written as

Elastic rods
Internal (contact) force and internal (contact) moment are denoted by n and m, respectively: Their rotational pull backs are which are also called stress resultants. For the elastic case, assuming the existence of a scalar valued stored energy density function ψ rod (v 0 , k 0 ), the stress resultants become The second derivative of ψ rod (v 0 , k 0 ) with respect to the strain prescriptors returns the 6 × 6 tangent stiffness matrix C 0 which relates the change in the stress resultants to the change in the strain prescriptors as

Elastoplastic rods
In the sequel, the key components to model geometrically exact elastoplastic rods are briefly outlined (see [14,15] for details). The strain prescriptors are assumed to decompose into an elastic and a plastic part as follows: Assuming that the stored energy density function ψ rod is purely dependent on the elastic strain prescriptors, the free energy density rod is introduced as where H rod (ξ rod ) describes hardening depending on internal variables ξ rod . For the special case of ideal elastoplasticity H rod (ξ rod ) = 0 and the dependency on ξ rod can be omitted. From now on, for elastoplastic rods, we restrict to the ideal elastoplastic case. The dissipation power D rod is defined as the difference between the stress resultant power and the rate of the free energy density, i.e., where a superposed dot denotes the material time derivativė (14) holds for any admissible process, which renders the constitutive relation and the reduced form of the dissipation inequality, i.e., Let us denote by rod (n 0 , m 0 ) the yield surface for rods. Maximising Eq. (16) under the constraint rod ≤ 0, one obtains the following evolution equations for the plastic strain prescriptors: Here,γ represents the Lagrange multiplier, ensuring that the stress resultants are admissible and satisfy the following Karush-Kuhn-Tucker (KKT) conditions The 6 × 6 elastoplastic tangent stiffness matrix C ep 0 relates the change in the stress resultants to the change in the total strain prescriptors as We note that a detailled knowledge of a suitable yield surface in terms of stress resultants is crucial to successfully implement geometrically exact elastoplastic rods. This contribution is dedicated to the determination of such yield surfaces.

Incorporating cross-sectional warping
An approach to incorporate cross-sectional deformations and nonlinear elastic constitutive relations into the formulation of geometrically exact rods has been presented in [25] and further developed in [26]. It has, however, not yet been generalised to the case of inelasticity. Here, the approach is extended to elastoplasticity.

Outline of cross-sectional warping
The theory presented in [25] and [26] allows for the crosssection 0 of a rod to deform due to strain prescriptors, i.e., the cross-section may warp. In general, a rod deforms such that the strain prescriptors vary along the rod's arc length s. Thus, the deformation of the cross-section at s depends not only on the strain prescriptors at s but also on the strain prescriptors in the neighborhood of s. If, however, the strain prescriptors vary slowly along s, the deformation can be assumed to depend only on the strain prescriptors at s. Eventually, a uniformly strained rod is constructed, i.e., v 0 (s) = v 0 and k 0 (s) = k 0 : the strain prescriptors are no more dependent on the arc length. A uniformly strained rod takes the shape of a helix [25]. The deformation map of such a rod is where the expressions of r(s) and R(s) for a uniformly strained rod are derived and presented in detail by Kumar et al. [25]. Here,X(X α ) denotes the deformed material configuration of the cross-section. Due to the assumption of uniform deformation,X is not dependent on s and can be rewritten aŝ Here, U i (X α ) denotes the cross-sectional displacement (warping) of the material configuration. According to the map (20), the deformed cross-section at s is obtained by rotatingX by R(s).
To evaluateX for the purely elastic case, an energy approach is applied: the stored energy density ψ F in the rod's cross-section is minimised with respect toX, thus leading to The deformation gradient F of map (20) is If cross-sectional displacement (warping) is suppressed, i.e. U i = 0, one gets back to the deformation gradient introduced in Eq. (6). Equation (23) highlights why v 0 and k 0 are collectively denoted as strain prescriptors: they prescribe the deformation and thus the local strains in the cross-section.
To assure uniqueness, the minimisation problem (22) has to be solved under two sets of constraints. Furthermore, the constraints ensure that the imposed strain prescriptors do not change during minimisation (for more details, the reader is referred to [25]). The first constraint fixes the mass center of the cross-section: The second constraint ensures that the average rotation about E 3 vanishes and that the principal axes of the cross-section align with E 1 and E 2 . where Here, ρ denotes the mass density. Consequently, to solve the constrained minimisation problem, a constrained energy functional is defined: with the Lagrange multipliers λ and μ enforcing the constraints (24) and (25). The first variation of the functional (27) leads to an Euler Lagrange equation describing the cross-sectional balance of linear momentum with traction free boundary: Here, P denotes the Piola stress andN denotes the outwards pointing normal to ∂ 0 . Div α denotes the divergence in the plane of the undeformed cross-section. Furthermore we introduce the following abbreviation Restricting to the case of isotropy in this paper, the stored energy density ψ depends on the left Cauchy-Green strain b = FF T only. The Piola stress then becomes P = τ F −T with the Kirchhoff stress following from τ = 2b ∂ψ ∂b (b).
For inelasticity, the energy approach (27) can not be applied. However, one can make use of the cross-sectional balance of linear momentum (28) which holds always. Together with the constraints (24) and (25) one can then solve for the elastoplastic cross-sectional deformation. The constitutive equation resulting from the dissipation inequality renders where denotes the free energy density. We note that for inelasticity, τ and thus P is no more a function of b only, but a function of the elastic left Cauchy-Green strain b e and an internal variable ξ . An associated plasticity model, incorporating the evolution equations forḃ e andξ is comprehensively discussed in "Appendix A".
In the sequel the problem to solve for the elastoplastic cross-sectional deformation is denoted as Q(v 0 , k 0 ).

Stress resultants and stiffnesses
In three-dimensional continuum mechanics, the normal projection of the Piola stress P on to the undeformed boundary returns the traction T. The cross-section in its material configuration has a normal aligned with the E 3 direction, thus here T = PE 3 . Integrating T over the cross-section returns the (contact) force Integrating the cross product of T with the position vectorX over the cross-section results in the (contact) moment Here, it is obvious why n 0 and m 0 are collectively denoted as stress resultants. They result from integration of stress measures over the undeformed cross-section and exclusively depend on the strain prescriptors. As demonstrated in [26], the rod's stiffness is obtained by taking the derivative of the stress resultants with respect to the total strain prescriptors. Here, the elastoplastic stiffness is obtained. Let p ∈ v 0 1 , v 0 2 , v 0 3 , k 0 1 , k 0 2 , k 0 3 denote any individual component of the total strain prescriptors. Elastoplastic stiffnesses are then obtained with and The elastoplastic tangent A ep = ∂P ∂F in the above expression is not straightforward to obtain due to the dependence of P on internal variables (see "Appendix B.3" for details). Similarly, obtaining ∂X ∂ p and ∂F ∂ p is not trivial, the reader is referred to [26]. Eventually, the 6 × 6 elastoplastic tangent stiffness tensor C ep 0 is arranged as Note that C ep 0 is symmetric for associated plasticity.

Elastoplastic response of rods
Now the elastoplastic response of rods is studied. Therefore, stress resultant versus strain prescriptor plots for cyclic loading are discussed and the progression of plastification in the cross-section is analysed. Later, the elastoplastic longitudinal and bending stiffnesses are examined. The continuum inelastic constitutive response of the material is characterised by the J2-flow theory [27] which is summarised in "Appendix B.2". The free energy density takes the form with the stored energy density ψ(b e ) being a function of the elastic left Cauchy-Green strain b e = F e F e T . Here, the stored energy density ψ(b e ) is defined in terms of the logarithmic principal stretches as shown in "Appendix B.2", Eq. (B.11). H(ξ ) denotes a hardening contribution, depending on the internal variable ξ and a constant hardening parameter H . Here, The continuum yield criterion, describing the region of admissible stresses, is defined as the von-Mises yield criterion Cyclic loading is performed to obtain the well-known hystereses loop in terms of stress resultants and strain prescriptors. Simulations are firstly performed with H = 0, i.e., for ideal plasticity, and subsequently with increasing H > 0. Figure 2 shows results for different values of the hardening parameter. For ideal plasticity, it is obvious that the curves for renewed loading coincide with the curves describing previous loading. The slope of the curves is horizontal in the plastic region. When incorporating hardening, however, the slope rises. This trend is clearly discernible in all four plots. The major difference in the four plots is between (a) and the remaining (b) to (d). In (a), the longitudinal straining yields uniform deformation in the cross-section, which leads to simultaneous plastification of the whole cross-section. In the other load cases (b) to (d), twist, bending and shearing do not lead to uniform deformation of the cross-section. Thus, different regions of the cross-section plastify earlier or later, leading to a gradual transition from purely elastic to plastic behaviour, when expressed in stress resultants.
The plastification of the different regions in the crosssection can be visualised by the distribution of the plastic work density ψ p . The plastic work density is obtained with The rate of plastic work density is defined in the appendix by Eq. (A.10). Figure 3 displays the distribution of ψ p for a cross-section undergoing longitudinal strain, shear, twist and bending respectively. Here, H = 0.1μ. The discrete figures show the distribution for t = 0.2, t = 0.8 and t = 1. The prescribed strain prescriptors are set to v 0 3 = 0.01t, v 0 1 = 0.01t, k 0 3 = 0.01t and k 0 1 = 0.01t. The color plots nicely visualise that the distribution of ψ p is not uniform in the whole cross-section, except for the case of longitudinal straining. The nonuniform distribution results in the effect displayed by the hystereses curves in Fig. 2: the transition from elastic to plastic deformation is gradual for the case of shearing, bending and twisting, but instantaneous for the case of longitudinal straining.
In the sequel, the stress resultants and the corresponding stiffnesses are analysed in detail as a function of strain Increasing hardening leads to increasing stiffness in the plastic region.
In contrast to the remaining curves, the curves in (a) show a prompt transition between the elastic and the plastic deformation state. There, the whole cross-section plastifies simultaneously, whereas for shearing, bending and twisting, the plastification is not simultaneous throughout the whole cross-section. Thus the transition between elastic and plastic region is gradual  Fig. 4 Axial force and longitudinal stiffness as a function of longitudinal strain. A dramatic decrease in stiffness is observed when plastic deformation is attained. The higher the hardening parameter H , the lower the decrease in stiffness. It is remarkable that the stiffness gets not equal zero but negative for ideal plasticity. This phenomenon is explained with the presence of cross-sectional shrinking, a geometrical softening effect. Due to tension, the cross-section of the rod reduces in area prescriptors. First, the resulting axial force and the axial stiffness for longitudinal straining is displayed in Fig. 4. In (a) one can observe that the axial force increases linearly in the elastic region. Linear behaviour is also observed in the plastic region. In (b), the axial stiffness is displayed. It is equivalent to the slope of the curves shown in (a). It is obvious that the stiffness decreases drastically, when the material starts yielding. The decrease is less pronounced when the hardening parameter is high. It is notable that for the case of ideal plasticity, the stiffness gets slightly negative. Ideal plasticity would suggest that the stiffness reduces to zero. However, due to geometric nonlinearity, the cross-section of the strained rod shrinks. Thus, the cross-sectional area is reduced and geometrical softening is observed. Figure 5 displays the bending moment as well as its corresponding stiffness as a function of bending. In the elastic region, a linear relation between moment and curvature gets visible (a). Then, the transition from purely elastic to plastic behaviour is gradual, due to nonuniform plastification of the cross-section. In contrast to the longitudinal stiffness (Fig. 4), the bending stiffness decreases gradually and remains positive, even for the case of ideal plasticity (b). This is not the case for the longitudinal stiffness, where it has been shown that the stiffness gets slightly negative. For completeness, stress resultants and their corresponding stiffnesses for the case of shear and twist are displayed in "Appendix C".
Highly pronounced negative longitudinal stiffness can be obtained by increasing the longitudinal strain and raising the yield stress to σ y = 45. Even if this material does not represent a physical material, it visualises nicely, the effect that can occur. Figure 6a depicts the axial force as a function of longitudinal strain for different hardening parameters. The high strain results in a pronounced shrinking of the cross-  Fig. 4 is intensified. Nonlinearities due to geometrical softening are observed in the elastic as well as in the plastic deformation state. Negative stiffnesses occur in the plastic region section. Geometric softening is observed both for elastic as well as for plastic deformation. In (b) it is observed that the stiffness decreases from the beginning and reaches negative values when the rod is deforming plastically. In contrast to the results shown in Fig. 4, negative stiffness is present for all hardening parameters at hand. Negative stiffness can cause difficulties in load prescribed simulations, since the load can not be related uniquely to the deformation state.

Yield surfaces in terms of stress resultants
In Sect. 2, we showed that a suitable yield surface in terms of stress resultants rod = rod (n 0 , m 0 ) is required to properly model geometrically exact elastoplastic rods. The yield surface describes a five-dimensional surface in the sixdimensional stress resultants' space. Eventually, the yield criterion describes the region of admissible stress resultants as rod (n 0 , m 0 ) ≤ 0. In the following, a framework to determine yield surfaces in terms of stress resultants is outlined and the resulting yield surfaces are comprehensively discussed.

Determining yield surfaces in terms of stress resultants
Aiming to obtain a smooth yield surface in terms of stress resultants, an analytic function is fitted to discrete points on the yield surface, which are characterised by combinations of stress resultants n 0 y and m 0 y at the yield limit.
Here, the stress resultants at yield result from discrete simulations of the cross-sectional warping problem Q(v 0 , k 0 ). As the entire cross-section does not plastify simultaneously, the yield limit is defined to be the point when the plastic work W p in the cross-section reaches a prescribed limit W if W p < W p limit then increase load factor l; else set current stress resultants to stress resultants at yield n 0y m 0y T = ld;

break; end end
Here, W p limit is set to be the plastic work corresponding to v p 0 3 = 0.002 for longitudinal loading, which corresponds to the yield strength. Other definitions of the yield limit are also conceivable. For example, instead of defining the limit in terms of the plastic work, one could define it in terms of the plastification ratio of the cross-section.
In the following, discrete points on a yield surface are obtained for rods with different cross-sections. Material parameters are chosen as summarised in Sect. 3.3 and the hardening parameter is set to H = 0.1μ. With non-zero hardening, geometric softening as shown in Figs. 4 and 6 is avoided. Furthermore, to determine the discrete points on the yield surface, the direction in the stress resultants' space d is chosen such that d lies in the x − y plane of the six-dimensional stress resultants' space, with x, y ∈ {n 0 1 , n 0 2 , n 0 3 , m 0 1 , m 0 2 , m 0 3 }. In every plane 48 discrete points on the yield surface are evaluated. Together they define the discrete yield surface.
Eventually, the yield surface in terms of stress resultants is described by a smooth function. Here, a generalised fivedimensional Lamé curve with twelve unknown parameters {a, b, c, d, e, f , α, β, γ, δ, , ζ } is fitted to the discrete data. Later in this contribution, it will get visible that above function fits the discrete yield surface in a satisfactory way and stands out for its simple structure.
In the sequel, the index • 0 is generally omitted for the sake of better readability. Figure 7 shows the discrete yield surface data resulting from the cross-sectional warping problem Q(v, k) for different combinations of stress resultants (blue dots). Here, the crosssection is circular with the radius r = 1. Observe that the shapes of the discrete yield surfaces reflect the symmetry of the cross-section.

Fitting a yield surface to the discrete data resulting from Q(v, k) for a circular cross-section
Here, symmetry reduces Eq. (41) to a Lamé curve with only eight parameters such as With a least-squares optimisation method, one can now fit Eq. (42) to the numerically obtained discrete yield surface.
To ensure that the yield surface rod is convex, the following condition has to be fulfilled Two-dimensional projections of the continuous yield surface are shown in Fig. 7 in red lines. The projections are in good agreement with the discrete yield surface resulting from Q(v, k). Discrepancies can be observed in the projections which deviate from a smooth elliptic shape (shearingbending and tension-bending). However, these deviations are considered as reasonably small. In the following, Eq. (42) is interpreted from a geometrical and mechanical point of view. Geometrically, the values of the variables a, c, d and f represent the length of the half axes of the six-dimensional Lamé curve. From a mechanical point of view, they represent the magnitude of the stress resultants at yield for the load cases where all stress resultants except for one equal zero. The exponents define the shape of the Lamé curve and can not be interpreted in mechanical terms. The generalised form of the yield surface eventually becomes where the stress resultants with subscript • y denote the values at yield. This equation shows that it is sufficient to determine the yield limit for simple loading, where only one stress resultant is not-equal to zero, and fit the exponents to get the analytic yield surface.

Remark
We note that this contribution focusses on a yield criterion for geometrically exact ideal elastoplastic rods, i.e., without considering any hardening at the rod level. With hardening, the structure of the yield surface (45) would have to be extended by including the thermodynamic forces conjugated to internal hardening variables One such possibility is where f (q 0 ) describes isotropic hardening and the remaining force conjugates describe kinematic hardening. To extend the yield surface by the thermodynamic force conjugates, as shown in (47), H rod ξ rod has to be known. Such a function can be evaluated by comparing the yield surface resulting from different values of the limit plastic work W p limit . This will be a topic of subsequent publications.
Having derived a yield surface for a rod with circular cross-section with r = 1, following question arises: how does the yield surface change, if the cross-section dimension is scaled by a scalar factor ϑ? In general it is expected that the stress resultants at the yield limit increase for ϑ > 1 and decreases for ϑ < 1. To compare the resulting discrete yield surfaces, one has to normalise the stress resultants with respect to ϑ. Therefore, Eqs. (31) and (32) are recalled, where it is obvious that the forces scale with the square of ϑ and the moments scale with its cube. The discrete yield surfaces resulting from the cross-sectional warping problem Q(v, k) for different ϑ are plotted in Fig. 8, where the forces are normalised by the square, and the moments by the cubes of ϑ (for the sake of visibility, only the connections between the One can observe that the normalised yield surfaces do not show a major dependency on ϑ discrete points are plotted). It appears that the normalised curves fit one another almost perfectly. The reason for the good agreement of the curves is probably in the choice of the material parameters at hand. Here, yielding already appears at small deformations. Thus, the deformation state is almost geometrically linear and the relation between stresses and deformation scales linearly.
The fact that the normalised discrete yield surfaces are virtually identical leads to the conclusion that the scaling shall nearly not affect the values of the exponents in the yield surfaces. However, comparing in detail the exponents resulting It is obvious that the forces scale with ϑ 2 and the moments with ϑ 3 . The exponents determining the shape of the yield surface vary slightly from the fit of the general yield surface (Eq. 45) to the discrete yield surface, one can not fully confirm this assumption (see Table 2). While the half axes (n 1 y , n 2 y , n 3 y , m 1 y , m 2 y , and m 3 y ) scale with ϑ 2 and ϑ 3 , respectively, the exponents (α, γ , δ and ζ ) show slight deviations. Nevertheless, the range of the exponents remains similar. This observation allows for the following statement: If the yield surface of a rod with given circular cross-section is known, one can easily obtain the yield surface for a circular cross-section of different size without big loss in accuracy by simply scaling the yield surface. Later, it can be seen that the scaling also applies for squared cross-sections (see Fig. 9).

Fitting a yield surface to the discrete data resulting from Q(v, k) for a squared cross-section
Now, discrete points of the yield surface for a rod with squared cross-section are determined. The symmetry of the cross-section will be reflected in the discrete yield surface. Eventually, the parameters of the continuous yield surface (45) are fitted to the discrete data. Figure 9 shows the normalised discrete yield surfaces obtained from the cross-sectional warping problem Q(v, k) for different crosssectional sizes (again, for the sake of better visibility, only the connection lines between the discrete points are plotted). The edge length of the squared cross-section is d = 2ϑ. Again it appears that ϑ does not show any considerable impact on the normalised yield surfaces. The parameters of the general yield surface (Eq. 45) are fitted to the data visualised in Fig. 9. The resulting parameters are summarised in table 3. In analogy to table 2, it is obvious that the forces and moments at yield scale with ϑ 2 and ϑ 3 , respectively. The exponents do not remain equal but vary slightly with ϑ.

Comparing the fitted yield surface for a squared cross-section to existing formulations of yield surfaces
Let us now compare the fitted yield surface of a rod with a squared cross-section to existing formulations of yield surfaces. Especially, we compare to a formulation used by Simo et al. [9] and another used by Duan et al. [16]. The crosssection has the edge length d = 2.
Simo et al. [9] suggests a yield surface for a rod deforming in the E 1 − E 3 plane. The yield surface is obtained by applying upper and lower bound techniques of limit analysis, which does not require to solve the cross-sectional problem. Thus, the yield surface there results from a different approach than the approach introduced in this contribution. The yield surface in Simo et al. reads: where the stress resultants subscripted with • p denote fully plastic stress resultants, i.e. the stress resultants leading to a fully plastified cross-section. However, in the present framework the definition of yield does not necessarily lead to a fully plastified cross-section (see Sect One then can look for similarities between Eqs. (49) and (45). In both cases, the stress resultants are normalised with their corresponding yield limit and raised to the power of a variable exponent. The obvious difference in both expressions is in the mixed term n 3 n 3 y 2 n 1 n 1 y 2 . This term does not appear in Eq. (45). In Fig. 10, projections of the yield surface (49) are compared to projections of the yield surface obtained from the fitting of Eq. (45). The discrete yield surface resulting from the cross-sectional warping problem Q(v, k) serves as the reference. In the n 1 − n 3 and n 1 − m 2 projections, both definitions of the yield surface fit the discrete yield surface equally well. We also note that the modified yield surface from Simo et al. matches the discrete yield surface slightly better in the n 3 − m 2 projection. However, the difference is not large.  Fig. 8, the shape of the curves differ slightly. In analogy, however, it can again be observed that the normalised yield surfaces do not show a considerable dependency on ϑ It gets obvious that the forces scale with ϑ 2 and the moments with ϑ 3 . The exponents determining the shape of the yield surface vary slightly

Fig. 10
Comparison between the yield surface used in Simo et al. [9] and the yield surface presented in this contribution with the discrete yield surfaces resulting from the cross-sectional warping problem Q(v, k).
The squared cross-section has edge length d = 2. The plots show that Simo et al. approximates the yield surface slightly better than the above presented yield surface Duan et al. [16] proposes a yield surface describing biaxial bending and stretching as follows In analogy to Eq. (48), the yield surface is described in terms of fully plastic stress resultants The parameters α 1 , α 2 , β 1 and β 2 are introduced as α 1 = α 2 = 1.7 + 1.3 n 3 n 3 p and β 1 = β 2 = 2 for squared cross-sections [16]. Again, we naively replace the fully plastic stress resultants in expression (50) with the stress resultants at yield evaluated from Q(v, k). The resulting yield surface is rewritten as Comparing Eq. (51) to the formulation presented in this framework (45), a significantly more complex structure stands out. In Fig. 11 projections of the yield surface (51) are compared to projections of the yield surface obtained from the fitting of Eq. (45). The discrete yield surface serves as a reference. It gets obvious that the yield surface proposed by Duan et al. fits the discrete points slightly better in the n 3 −m 1 and n 3 − m 2 projection than the yield surface introduced in this contribution. Again, the difference is not large. When comparing different formulation of yield surfaces one has to keep in mind that they possibly emerge from different approaches. In Simo et al. [9] and Duan et al. [16], the yield surface is described in terms of fully plastic stress resultants. In the present framework the yield surface is described in terms of stress resultants at yield defined by the plastic work. It is all the more astonishing that the yield surfaces used in Simo et al. and in Duan et al. fit the discrete yield surface, when naively replacing the fully plastic stress resultants with the stress resultants at yield. The yield surfaces described by Simo et al. and by Duan et al. just define a twodimensional surface in three-dimensional stress resultants' space. In contrast, the yield surface derived in this contribution shows a five-dimensional surface in six-dimensional stress resultants' space. Thus, it is more general in the stress state of the rod. Furthermore, the yield surface presented here is smooth and of very simple structure. From the point of view of implementing geometrically exact elastoplastic rods, this is convenient, since derivatives with respect to the stress resultants are required. With the presented Lamé curve, the remaining projections also match the numerical data. This can not be investigated with the yield surface used in Simo et al. or in Duan et al.. Finally, here, the cross-sectional problem is solved. This ensures that the deformation of the crosssection is taken into account, when deriving yield surfaces.
As a final note, let us remark that the shape of the discrete yield surface may change depending on the chosen yield limit. In fact, its shape depends on the geometry of the crosssection, too. At the same time, it is also not guaranteed that the discrete yield function can be fitted by a six dimesional convex Lamé curve for all possible cross-sections. In that case, one has to use a different function that replicates the discrete yield surface best. The process to obtain the discrete yield surface, however, would remain unchanged.

Summary
A framework to obtain a yield surface in terms of stress resultants has been presented with the objective to model geometrically exact elastoplastic rods. Based on the crosssectional deformation of a strained rod, stress resultants are computed. Here, the cross-section obeys a J2-elastoplastic constitutive model. The discrete yield surface results from performing sets of load controlled simulations of the crosssection. Eventually, smooth analytic functions are fitted to the discrete yield surface. The agreement between the fitted functions and the discrete yield surface obtained from simulations is remarkable. The fitted yield surfaces have been obtained both for circular and squared cross-sections. We have shown that in both cases, the yield surfaces are scalable with the size of the cross-section. For an exemplary squared cross-section, the resulting yield surface has been compared with the existing description used by Simo et al. and Duan et al., showing excellent agreement. The yield surface obtained from this contribution can now be used in a finite element formulation for geometrically exact elastoplastic rods such as the one presented in [15] for a more accurate simulation. The extension of the proposed scheme to yield surfaces with hardening will be investigated in a forthcoming contribution. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indi-cate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix A. Three-dimensional finite deformation elastoplasticity theory
This appendix reviews elastoplasticity at finite strains. The elastic domain is specified by means of a yield criterion, which is expressed in stresses defined in the deformed spatial configuration. For the case of small strains, the strain tensor ε can be additively decomposed as ε = ε e + ε p . The resulting stress is then only dependent on the elastic strains ε e and internal variables ξ . For the case of large deformations, an additive split of the strain measures is no more permitted. Instead, a multiplicative decomposition of the deformation gradient is used as introduced in [28] and well-rooted in the modelling of crystal-plasticity e.g. [29].

Appendix A.1. Kinematical background of multiplicative plasticity
To characterise the kinematics of finite deformations, the deformation map ϕ(X, t) is introduced. It maps any point X from the undeformed material configuration B 0 to its new position x = ϕ(X, t) in the deformed spatial configuration It is assumed that the stored energy density ψ depends on the elastic deformation. Restricting to isotropy, ψ is a function of the elastic left Cauchy-Green strain b e only, i.e.
A possible hardening contribution H depends only on the internal variable ξ . Then, the free energy density is introduced as follows: The dissipation power D is defined as the difference between the stress power and the rate of the free energy density and is greater or equal zero Here, τ denotes the (symmetric) Kirchhoff stress. Equation with L υ b e := FĊ p −1 F T denoting the Lie derivative of b e . Equation (A.7) has to hold for any admissible process. This standard argument renders the constitutive relation and the X ϕ(X, t)  The derivative in the second sum is the pull back of its spatial counter part derived in [27] ∂ M tr The derivation of a ep is given in [27] and is summarised as Implicit to the above equation is: Obviously, for the case of elasticity γ = 0 and a ep = a. The transition from purely elastic to plastic behaviour is gradual. For ideal plasticity as well as for plasticity with hardening, the decrease in stiffness is smooth tion of shear and twist, respectively. Their course is similar to the case of bending discussed in Sect. 3.3.