Structural response of existing spatial truss roof construction based on Cosserat rod theory

Paper presents the application of the Cosserat rod theory and newly developed associated finite elements code as the tools that support in the expert-designing engineering practice. Mechanical principles of the 3D spatially curved rods, dynamics (statics) laws, principle of virtual work are discussed. Corresponding FEM approach with interpolation and accumulation techniques of state variables are shown that enable the formulation of the C0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^{0}$$\end{document} Lagrangian rod elements with 6-degrees of freedom per node. Two test examples are shown proving the correctness and suitability of the proposed formulation. Next, the developed FEM code is applied to assess the structural response of the spatial truss roof of the “Olivia” Sports Arena Gdansk, Poland. The numerical results are compared with load test results. It is shown that the proposed FEM approach yields correct results.


Introduction
The paper presents the use of nonlinear theory of Cosserat rods and the resulting finite element computer code as tools for analysis of statics and dynamics of engineering structures, as the main example the spatial truss being the roof of "Olivia" Sports Arena in Gdansk, Poland, is analyzed. The formulation presented here has been used as a part of an advanced structural health monitoring system (SHM) [1,2], yet this paper focuses only on the theoretical formulation and the comparison of the numerical results with test load measurements.
Structural health monitoring systems are nowadays commonly used to confirm theoretical assumptions made in design phase [3][4][5][6][7]. The main purpose of the SHM system here was to ensure the safety of visitors by initiation of emergency procedures in case of dangerous near-failure states that might lead to the total collapse. Because of this it was necessary to implement solution that would be able to detect of type of near-failure modes of the roof. The idea of the SHM system is shown in Fig. 1. Briefly, the designed SHM system relied on own Fortran FEM code that was constantly performing nonlinear calculations evaluating the state variables of the roof under the current snow weight. The use of the own code (called B6 [2]), in comparison with existing commercial FEM packages, was motivated by its open structure and full integration with the expert monitoring system. Unquestionably, such solution was cheaper in comparison with the use of commercial FEM system and did not require costly renewal of the license. Other relevant approaches can be found in [8,9]. As indicated previously, for brevity the details of SHM system are omitted in this paper. These details will be presented in the forthcoming papers. Here the theoretical aspect of formulation, numerical implementation and validation of the FEM results is presented.
As the basic structural element of the truss systems, the Cosserat curve model is used. This model is of the general model of rods and beams presented in the literature, see [2,[10][11][12][13][14][15][16][17][18][19][20][21]. Developed from the Cosserat rod theory, B6 code naturally accommodates nonlinear analysis where translations and rotations are unlimited from the theoretical viewpoint. The underlying Cosserat theoretical foundation of the code presents modern approach to the rod theory. Apart from necessary mathematical limitations (smoothness, differentiability, continuity), the theory is free from redundant kinematical assumptions that are typical for classical theories. The approach is exact as far as kinematical and dynamical resultants for the rod-like body are concerned. The only stage of formulation which has approximate character is the constitutive relation. The truss structure under study could be treated also as a continuous body after a proper homogenization procedure as, for example, presented in [22][23][24][25][26][27][28][29]. These kinds of procedures may produce higher gradient or non-local models with uncommon mechanical responses which should be taken into account because some localization effects could not be negligible at macroscale for the size of the length scale at micro-level (see, e.g., [30][31][32][33][34][35]).
The structure of the paper is as follows: • In Sect. 2, relevant mechanical principles of the 3D spatially curved rods, dynamics (statics) laws and the principle of virtual work (PVW) are presented. The latter one reveals the kinematics of the rod which is equivalent to that of the Cosserat curve. • In Sects. 3-5 are discussed details of linearization of the principle of virtual displacements, interpolation and accumulation techniques of state variables, and formulation of beam and truss (cable) Lagrange finite elements, in the context pertaining to the FEM approximation.
• Test examples presented in Sects. 6 and 7 confirm the correctness and suitability of the proposed formulation.
The developed B6 code is applied to the analysis of the Olivia Sports Arena roof. The obtained numerical results presented in this section are compared with in situ measurements of static deflection and natural frequencies of the analyzed structure. The comparison of the results shows the excellent performance of the B6 code in the engineering practice.

Kinematics of rod, weak form
The idea of rod as the Cosserat directed curve receives major attention in the literature, see [2,[10][11][12][13][14][15][16][17][18][19][20][21]. In comparison with work [14], which was developed within research team of Gdansk University of Technology and creates the background of the present work, the study is extended to accommodate the mass matrix. However, the dynamic analysis is limited to natural frequency extraction and time integration algorithms on SO(3) group are not discussed. Some interesting references in this area are to be found, for instance, in [36,37].
To fix the ideas of the paper, let us define a rod-like body B (Fig. 2) as a (directed) curve l imbedded in the 3D Euclidean space E 3 , (o, {e i }), i = 1, 2, 3. Each point s ∈ l is endowed with non-coplanar triad of directors {t i (s, t)}, assumed here as orthonormal, or equivalently, with a proper orthogonal structure tensor (3), it can be obtained as a proper orthogonal transformation of the fixed frame e i , i.e., t i (s, t) = T (s, t)e i . Each point s ∈ l is associated with the cross-sectional area A > 0. The motion of l is understood as a sequence of configurations l(t), t ∈ T , where t denotes the time. As the reference configuration, the position of l at t = t 0 is taken. Let s be an arc length along l 0 . Then denotes the position vector of an arbitrary point of B. Here x(s) denotes the initial placement of the axis and ζ 0 (s) = ξ α t 0 α (s), ξ ≡ {ξ α }, α = 1, 2 denotes the position of arbitrary point of the cross section A(s) ≡ A 0 (s). For the purposes of this paper, it is assumed that the tangent vector to the initial axis of the rod l ≡ l 0 is defined by To establish dynamic equilibrium conditions, the idea from [38,39] was followed, which was later used in [40] to derive a complete set of field equations for the 6-parameter nonlinear shell theory, or in [14] for the 6-parameter rod theory. For the clarity of the exposition, the reader should refer to the cited references for the details pertaining to the equations of motion.
To formulate the weak form necessary for FEM implementation, it was proceed in the manner discussed in [14] (compare also [40,41] for shells). This not only yields the weak statement but also reveals the formulae for virtual strain measures. Some details about the general structure of the strong formulation may be found, for example, in [42][43][44][45]. Ultimately, the following functional of the principle of virtual work is used where v and w are kinematically admissible virtual displacements and virtual rotations, respectively. The term G i (l, t; v, w) = l (n · (v − w × y) + m · w)ds being the internal virtual work defines vectors of virtual strain measures (stretching strains and change of curvatures) in the spatial representation which are energetically conjugated with internal stress n and couple vectors m, respectively. Here y denotes the current placement of l(t). It can be shown (cf. for instance [14] for rods and [41,42] for shells) that corresponding Here ad −1 denotes the axial vector of skew-symmetric tensor (. . .). Theoretical studies [38][39][40][41]43] reveal that the kinematical model described by (5) is identical with the 1-D Cosserat continuum, where the motion is described by two sets of equations (see Fig. 3): Differentiating the first equation of (7) 2 with respect to time and assuming the spatial representation yields the formula for skew-symmetric tensor of angular velocity Ω (cf. for instance [46]) and its axial vector ω The skew-symmetric tensor of angular acceleration in spatial representation A and its axial vector have the following formula (cf. for instance [46]) Let us recapitulate the present developments by setting energy-conjugated objects in matrix notation: generalized displacements and their virtual counterparts, velocity and accelerations, respectively External point loads, external distributed loads, boundary terms and inertia forces are set in the following vectors The strain vector, the virtual strain vector and the internal (generalized) stress vector are defined Furthermore, the following operators are introduced: that enable the following definitions of virtual relations: With this notation, the principle of virtual displacements (4) can be rewritten as

Linearization
Linearization of the internal virtual work G i [u, t; w] and the external virtual work G e [u, t; w] in (16) has been extensively described in [14], and it also holds in the present formulation. Without losing generality, attention is confined to the dead load for which δG z = 0. Thus, the attention should be focused to the virtual work of inertia forces. It is necessary to calculate the variation of This is a complex calculation. Therefore, here are only the main steps required to arrive at the linearized form presented. In the first step, ω and a are treated as independent variables and The linearized tensor of inertia in view of implication is written as Hence, the term in bracket in (21) after some manipulation reads The term in [.] bracket from (21) follows from linearization of angular velocity and angular acceleration and may be written as As a consequence of the first step, the following expression is obtained: where the following matrices are introduced In the second step, ω and a are expressed as the functions of velocity ẇ and acceleration ẅ, (29) After some manipulations, and by making use of Lie bracket

Constitutive relations
In this paper, the linearly isotropic elastic rods are discussed. In such a case, the constitutive relations are given by the formulae in material representation (n, m, ε, κ; In (32), the constitutive matrices C 1 and C 2 do not depend on curvature of the rod axis. The component form of (32) follows from decomposition of the stress and couple resultants Here n 3 ≡ n is the axial force and n 1 , n 2 are the shear forces. Similarly, m 3 ≡ m is the torsional moment, while m 1 and m 2 are bending moments. The explicit form of (32) is Here μ, λ are the classical Lame constants and A i = k i A 0 , where k i is the shear correction factor and I i are the moments of inertia about axes i = 1, 2, respectively. Kinetic constitutive relations are assumed as where J 1 (x) is averaged time-independent mass tensor, while J 2 is the inertia tensor in spatial representation. In writing (36), it is assumed that the center of mass coincides with the axis of the rod.

Shape functions
Assuming that approximation functions are regular and smooth enough to satisfy necessary C 0 interelement continuity conditions, the spatial approximation of (16) may be carried out over typical finite elements In the reminder therefore the attention is confined to the single finite element. The In C 0 interpolation of vector variables, the following relations are employed: where L a (ξ ) are the Lagrange polynomials satisfying usual conditions n w a=1 L a (ξ ) = 1, L a (ξ b ) = δ ab , ξ ∈ π (e) , a = 1, 2, . . . , n w . In the standard FEM approach, nodal and element shape function matrices are defined byL with 1 as the second-order identity matrix. Since typical element Π (e) ⊂ E 3 is parameterized by ξ ∈ π (e) ≡ [− 1, + 1] ⊂ R, to calculate the strain field, it is necessary to find a transformation rule between derivatives ∂(.)/∂s ≡ (.) and ∂(.)/∂ξ . Through the chain rule, it is to find The notations a = ∂ x/∂ξ , t 0 3 ≡ t 0 = T 0 (ξ )e 3 are used. The transformation of the elementary length is given by dl ≡ ds = α(ξ )dξ .

Interpolation in SO(3) group, virtual rotations
From the theoretical viewpoint, the main difficulty of analyzed theory stems from the fact that definition of the configuration space that emerges in the nonlinear theory of rods contains Cartesian product of the Euclidean space and the rotation group SO(3), so it does not possess the character of linear space. This necessitates the use of nonstandard techniques of approximation and interpolation of variables.
The crucial aspect of C 0 interpolation is associated with objects from the rotation group, e.g., T (ξ ), T 0 (ξ ), Q(ξ ) ∈ SO (3). Typical interpolation (38) is an outer operation in SO(3) group. Hence, twostage interpolation originally proposed in [47] (compare with [48]) is employ indirect. The scheme is described in detail in [49] and references given there.
The vector of virtual rotations w ∈ T Q SO(3) ⊂ E 3 can be calculated in the fixed base e j by making use of the transformation rule t i (ξ, t) = T (ξ, t)e i . Formally therefore the following relation holds or equivalently in the matrix notation typical for FEM The treatment of the latter ones is more demanding as it depends on the employed parameterization SO(3) → R 3 of the rotation group. In this paper, the canonical interpolation is employed where W = ad w = wad(e), ||e|| = 1. As a consequence, at each node a = 1, 2, . . . , n w the resulting finite element has six engineering degrees of freedom (see Fig. 4), which variations form a nodal vector Considering the n-node element, one can write Similarly, vectors of velocity and acceleration are introduced Vectors of increments q and q are formulated by analogy to (44) and (45) that is q a ⇒ q (e) and q a ⇒ q (e) . In writing interpolation rules, it should be remembered that the object such as w(ξ, t), ω(ξ, t), a(ξ, t) belongs to different tangent spaces, i.e., different placements along the element length, and is being parameterized by time t ∈ [0, +∞): u(t) = (u(t), Q(t)). Therefore, the general formulae relating the vectors from different tangent spaces are given by Here (ξ ) depends on parametrization of rotations and can be defined in various ways, see, for instance, [50,51]. In this paper, (ξ ) = T(ξ) is chosen and thus Therefore, the interpolation formulae read In the case of generalized accelerations a, it is necessary to take into account that (ξ ) = T(ξ) also depends on time t through Q(ξ, t). Therefore, since a =ω = ẅ +˙ ẇ is obtained and the interpolation formula takes the form

Element matrices, numerical integration
Standard FEM argumentation leads to the following relations between expression for variation of the internal virtual work and its discrete counterpart on the element level where material and geometric stiffness matrices are given, respectively, by Similarly, for inertia forces (30) there is where the element mass matrix, the extended gyroscopic matrix and the centrifugal matrix are defined by The internal virtual work is calculated as In the same fashion, the element inertia force vector is obtained Eventually, on the element level the following equation is formulated . All the matrices, including the mass matrix, and vectors are computed numerically using Gauss quadrature. In order to minimize the locking effect in low-order elements and to avoid spurious zero-energy modes that accompany reduced integration (RI) 4-node elements with full integration (FI) are used. Standard aggregation procedure is then applied to obtain global equations for the whole structure.
To conclude the description of the formulation, it is mentioned in passing that the ultimate structure of the developed Fortran code of B6 program includes also the cable (truss) elements. These are formulated as 2-node Lagrange C 0 elements based on the classical Green-Lagrange strain tensor, cf. for instance [53].

Validation of the program
In [2], a number of tests were carried out to verify the correctness of the proposed approach and its FEM implementation. For the sake of brevity, only two of them in this paper are reported.

Right-angle cantilever
Following [52], stability analysis of the right-angle cantilever shown in Fig. 5 is performed. The following values are used: L = 240 mm, b = 30 mm, h = 0.6 mm, E = 71,240 N/mm 2 , ν = 0.31, P ref = ± 0.01 N . This is the very popular benchmark problem for both rod and shell formulations, see, for instance [12,52,[54][55][56]. To evaluate the critical load, the nonlinear stability analysis with imperfection introduced as the perturbation load P impf = P/1000 is used. The load deflection curves are presented in Fig. 6, whereas Table 1 shows the   [40,47,49]. Presented results, labeled as B6, are in good agreement with the reference solutions.

45-Degree bend cantilever
This example has been studied by various researches, e.g., [12,14,[57][58][59][60]. The radius of cantilever is R = 100, α = 45 • , see Fig. 7. The unit square cross section is assumed A = b × h = 1 × 1 = 1 made of material with properties E = 10 7 , ν = 0. In the FEM simulation, full integration (FI) or uniform reduced integration (URI) is used. Four 4-node B6 elements are employed. The analysis revealed that there is no significant difference between FI and URI results within own formulation B6, see Fig. 8. Again, good agreement between own and the reference solution is obtained, see Table 2.

Validation of B6 code, the roof of "Olivia" Sports Arena
The analyzed structure is created on the basis of reinforced concrete stands and a rod-cable roof structure with slender truss girders of lenticular shape. The main load bearing components of the roof are formed in a bold, lightweight and attractive spatial structure, built of eleven flat, inclined to each other, compressed lenticular truss girders, see Figs. 9, 10 and 11.
The span of the girders in plane is 78.54 m, and the structural height at the chord centerlines changes according to the parabola, in theory, from 0 m at the girders to 4 m in the middle of the span. With the compression of the lower chords, low structural height of 1/20 of the span was achieved, thus reducing the building space so important for air-conditioning of ice rinks. The upper chords are formed of a system of I-beams, and the lower chords are made of variably spaced channel bars in a parabolic arrangement, from      Consequently in 2008, technical survey of the structural health of the entire building was commenced, with particular focus on the roof structure [1]. The findings of the survey as well as detailed numerical analysis indicated the risk of "domino-type" collapse.
Since it was not possible to repair the construction immediately and the authorities could not close the Arena due to public pressure, it was decided to keep the Arena operational provided that a SHM system "safeguarding" the structure is enforced.
As mentioned in Introduction, the SHM system was enforced [1] where B6 was the core ingredient. For that purpose, prior to SHM enforcement, technical validation of the FEM formulation was carried out.

Single lenticular girder and FEM model
Firstly, a single lenticular girder was analyzed since there existed the results of the in situ tests of a single girder performed in 1969 on the building site. To account for complex spatial geometry, for the pre-and postprocessing phase SOFiSTiK [61,62] a general-purpose FEM code was used. This eased the development of the complicated spatial structure of the girder, see Fig. 12. In the calculations with B6 code, the following finite elements were employed: • 4-node spatial rod elements with full integration accounting for offset, • 2-node pre-stressed spatial truss-cable elements with no-tension feature, • non-dimensional contact element (between rod and cable).
The main analysis of the truss girder was preceded by convergence analysis which revealed that the mesh consisting of 2684 nodes, 800 rod elements and 632 cable elements yields appropriate convergence. For comparison purposes, SOFiSTiK [61] was also employed with the mesh built from 2 684 nodes and the following finite elements: • 2400 2-node spatial beam elements. Since [61] manual does not provide information, a detailed study of element properties was undertaken revealing that these are Timoshenko-type beam elements accounting for shear effect and axis eccentricity, • 632 2-node C 0 spatial truss-cable elements with no-tension feature.
In all calculations, the value of Young's modulus was assumed as 205 GPa and Poisson's ratio 0.3. The following load cases were studied: • gravity load, • pre-tension force 1030 kN in cables of a single-channel bar in the lower cord, • pre-tension force 1.96 kN of X bracing, • concentrated forces 5 kN on the upper cord at every node connecting the chord and vertical post along the total chord length (Fig. 13, case A), • as above, however along half of the total chord length (Fig. 13, case B).
All calculations were carried out in geometrically nonlinear range. The obtained results are presented in Fig. 14, where comparison is given to the original in situ measurements taken in 1969. The curves show a good qualitative agreement. The quantitative assessment revealed the 9% discrepancy between in situ and numerical results which is also fine. Table 3 presents some representative values to support this observation.
In the next step, the natural frequency extraction was performed. Here different numerical integration rules to obtain element mass matrix (55) 1 were tested. The obtained values are depicted in Fig. 15.   Fig. 16.
In load test, the static load was created by pairs of hanging bags filled with stones, and the weight of each bag was 9 kN. Their placement is shown in Fig. 17 along with the actual photograph in Fig. 18. The following load protocol was applied • 2 readings of displacements at every 15 min prior to load, • 1 reading directly after the placement of the next row of the hanging bags, • 1 reading directly after the placement of the whole load, • series of successive readings at each 15 min under the whole load, the load was kept as long as the increment of displacement in two consecutive readings did not differ more than 2% from the total load, • 1 reading after removal of the row of load, • 1 reading directly once the last portion of load was removed, • series of successive readings at each 15 min until the stabilization of the structure.   After convergence analysis, the overall B6 model consisted of 17156 nodes, 5250 4-nodes beam elements and 3708 2-nodes string elements. Some representative time histories of the displacements are compared in Fig. 19. Table 4 presents detailed numerical values.
Finally, natural frequencies of the whole roof were assessed. To this end, impulse load was generated by cutting off a single bag (9 kN). The representative acceleration time history is shown in Fig. 20 where FFT transformations is also given. The extracted frequencies along with associated modes are shown in Fig. 21 where comparison with numerical results is presented.  The above results confirmed the correctness of proposed formulation, the resulting FEM code and assumed discretization of the structure of the roof. Thus, B6 code was placed into SHM system as shown in Fig. 4.
The purpose of the installed SHM system [1,2] was to control the basic parameters of the roof structure and signaling its state in automatic way to appropriate authorities. It was novel technological solution as, in comparison with existing commercial offers, it included sophisticated computational module of structural analysis. The module was based on detailed theoretical model of the roof along with integrated expert system. Even though the number of sampling points was small, due to the performed real-time computations by B6, it was possible to supervise online the current state of the whole roof at arbitrary point. The advantage of own FEM code within SHM system was its low price (no renewal of license), robustness (unlimited translations and rotations) and small size of the executable file (below 1MB). The similar applications of the presented idea of structural health monitoring system may be found in [63][64][65].

Conclusions
Here is discussed the direct approach to the modeling of rod behavior which are free of any simplifying kinematical assumptions. Except of necessary mathematical requirements-continuity, differentiability and appropriate smoothness of fields entering the formulation-no other assumptions typical for classical formulations (such as Bernoulli hypothesis) are incorporated. The presented formulation is kinematically unique in the sense that the strain measures follow uniquely from the principle of virtual work and are not postulated. The only approximation on the theoretical level appears when material law is presented. The derived inertia terms are suitable for analysis of small oscillations superposed on finite deformations with unlimited translations and rotations. One of the crucial aspects of the formulation is the lack of linear structure of the configuration space C(E 3 × SO(3)) that is not a vector space. This necessitates the use of special interpolation techniques. Theory and the spatial approximation are consistent in regard to all simple states, i.e., tension/compression, bending, shear and torsion. Thus, connecting of non-planar elements is natural in the aspect of kinematics and dynamics and does not require artificial techniques.
Based on the Cosserat rods theory, the new nonlinear B6 FEM program presented in this article after calibration was finally used in SHM system that kept Sports Arena "Olivia" operational from 2008 till the construction works in 2010. The most important element of the system was the analysis and expert modules based on the updated nonlinear FEM models. The obtained detailed FEM model of the roof structure was obtained due to methodical derivation of the general geometrically nonlinear FEM model. The model was updated based on in situ tests as well as on numerous studies on the work done by temperature changes and snow load distribution. The program was also validated with few benchmark problems. Application of the code to the large engineering structure such as Sports Arena "Olivia", which is one of the largest rod-like structures, together experimental studies demonstrates the efficiency of the code.