Assessment of the E ﬀ ect of Wind Load on the Load Capacity of a Single-Layer Bar Dome

: The main purpose of the paper was the assessment of the e ﬀ ect of wind load on the load capacity of a single-layer bar dome. Additionally, which numerical method is appropriate for low-rise single-layer bar domes was checked. In order to explain the e ﬀ ect of the height-to-span ratio on the selection of the appropriate calculation model and method of analysis of the bar dome, an example of the known von Mises truss was proposed. Two cases of von Mises truss di ﬀ ering in the height-to-span ratio were considered. For the shallow structure, a signiﬁcant change in the value of the sti ﬀ ness matrix determinant and the current sti ﬀ ness parameter was observed. A similar tendency in the behavior of the structure can be observed on fragments of larger structures, including shallow single-layer steel domes. These problems are described on the basis of the dome, which is located on top of the building housing the restaurant. This structure is subjected to large displacement gradients and the actual conﬁguration is taken into account in analysis. The analysis showed that there is a change in sti ﬀ ness for these structures, and, therefore, that such structures should be designed according to geometric nonlinear analysis (GNA).


Introduction
For proper assessment of structural safety, it is necessary to select the appropriate computational model and methods of structural analysis. It is essential for the methods to most accurately represent the behavior of the structure. This concerns both the adoption of a static diagram, loading, and also a model of the behavior of the structure, bars, supports, and nodes under load.
The structural analysis means and scope are also highly affected by the behavior of the element cross-section in compression and bending. The process of structural design was streamlined due to the introduction of the cross-section classification. The class of the cross-section depends on the shape, slenderness of the element walls (width-to-thickness ratio c/t), the member yield strength, and the distribution of compressive stress. When compressive stress does not occur, it is unnecessary to specify the class of cross-section [1]. The cross-section class indicates the degree of the element resistance to the local stability failure in the elastic or plastic state. Consequently, it is of key importance for the selection of the computational model and the criterion of the cross-section capacity [2,3].
Global instability of the steel bars is equivalent to their failure. Global instability of compression element, called buckling, can have one of three modes: flexural (the bar is bent in the plane having the lowest stiffness), torsional (the bar twists around the longitudinal axis), and flexural-torsional (the bar bending is at the same time accompanied by its torsion). Global instability of the element in bending, called lateral torsional buckling, occurs when the bending moment reaches the critical value. Beam lateral torsional buckling is initiated by the compression flange buckling out of the bending plane, and that is immediately transformed into beam torsion. The possibility of stability failure depends on that are disregarded. The linear analysis does not account for the occurrence of overall stability failure of bars or the whole system, or the local stability failure of the elements of bar cross-sections [18,19]. This is equivalent to the assumption that the structure has ideal (not showing up any imperfections) initial geometry that does not change (the strain induced by load is not considered). The method takes into account elastic strain. In this method, displacements are proportional to the loads applied. The principle of superposition is followed in the assessment of internal forces in the structure and strength utilization in the cross-sections. In this model, the elastic capacity of the structure is restricted by the elasticity limit of the material (in practice, the yield point). The linear analysis in terms of the finite element method FEM boils down to the solution of a linear system of algebraic Equation (1): where: K L -linear stiffness matrix of the structure; q-vector of nodal displacements; P-vector of nodal load. The FEM software is additionally equipped with linear bifurcation analysis (LBA). The goal of LBA is to solve the eigenvalue problem: where: K G -geometrical stiffness matrix. The geometric stiffness matrix K G is identified as the initial stress stiffness matrix. Eigenvalues µ are critical load multipliers. Eigenvector q illustrates the form of structure deformation. Running the analysis allows the user to quickly verify how close to instability the analyzed load condition is. However, analysis is still conducted in the field of linear physical and geometrical relationships.
Only the geometrically nonlinear analysis (GNA) is capable of a full description of the possible forms of instability. This type of analysis can provide the equilibrium path and information about the post-critical behavior of the structure. For nonlinear discrete systems which are formulated on the basis of FEM, the mathematical model corresponds to a set of nonlinear algebraic equations. The set can be formulated in a total or in an incremental form. In the first case, equations have the form: where: K s -secant stiffness matrix of the structure. In the second case, equations have the form: where: K T (q)-tangent stiffness matrix of the structure; R = P − F-vector of residual forces; F-vector of internal forces. In the equilibrium state R = 0, while in the iterative process, the norm of R defines the distance from the equilibrium state. The iterative process converges if R→0.
The tangent stiffness matrix K T of the structure arises as a result of the assembly of the stiffness matrices of the elements: (K e L + K e G + K e u1 + K e u2 ) (5) where: K e T is the tangent stiffness matrix of the element composed of linear stiffness matrix K e L , geometric stiffness matrix K e G , and nonlinear stiffness matrices: K e u1 and K e u2 . In GNA, it is extremely important to determine changes in structure stiffness in the load process. An often used measure of structural rigidity is the matrix determinant det (K T ). A more effective proposition to describe these changes was proposed by Bergan and Soreide in [20]. They characterize the behavior of a multidimensional nonlinear system by using one scalar quantity, the current stiffness parameter (CSP), representing the ratio of two quadratic forms formulated for the tangent stiffness matrix at the initial K 0 T and current times K i T : This parameter can be used as a measure for the change in the tangent stiffness matrix K T during the motion in the N-dimensional, displacement solution space. Figure 1a shows a typical snap-through problem (load parameters versus some norm of displacement vector ||q||). The associated curve for CSP as a function of ||q|| is traced in Figure 1b. It is noticeable that at the extreme points of the load-displacement curve, CSP has the value zero. In this situation, the incremental stiffness matrix K T is singular. CSP is positive for the stable branches of the load-displacement curve. The instable configurations are characterized by negative values of CSP.
This parameter can be used as a measure for the change in the tangent stiffness matrix K T during the motion in the N-dimensional, displacement solution space. Figure 1a shows a typical snapthrough problem (load parameters versus some norm of displacement vector ||q||). The associated curve for CSP as a function of ||q|| is traced in Figure 1b. It is noticeable that at the extreme points of the load-displacement curve, CSP has the value zero. In this situation, the incremental stiffness matrix K T is singular. CSP is positive for the stable branches of the load-displacement curve. The instable configurations are characterized by negative values of CSP. Analysis of stability of bar structures with the finite element method (FEM) involves solving large systems of nonlinear equations. The relation between displacements and loads in a (N + 1)dimensional space is described by the load-displacement path. Elastic stability is closely related to singularities along this path. These singular points are considered to be critical ones. They include a limit point known as snapping and a bifurcation point known as buckling. The set of equations can be completed by a constraint equation to formulate an extended set of equations. An advantage of such a formulation is equivalent treatment of displacements q(η) and load parameter μ(η) where η is the so-called control parameter. The two most simple ways to control the incremental process are: µ = η and qi = η with the load parameter μ and chosen displacement qi being the controlling parameters. With the first of the two methods, subsequent points are the points of intersection between µ = ηi and the equilibrium path, and with the second method-as the points of intersection between qk = ηi with the same solution curve ( Figure 2). Both procedures are inefficient at the vicinity of local extremum of the controlling parameter. Analysis of stability of bar structures with the finite element method (FEM) involves solving large systems of nonlinear equations. The relation between displacements and loads in a (N + 1)-dimensional space is described by the load-displacement path. Elastic stability is closely related to singularities along this path. These singular points are considered to be critical ones. They include a limit point known as snapping and a bifurcation point known as buckling. The set of equations can be completed by a constraint equation to formulate an extended set of equations. An advantage of such a formulation is equivalent treatment of displacements q(η) and load parameter µ(η) where η is the so-called control parameter. The two most simple ways to control the incremental process are: µ = η and q i = η with the load parameter µ and chosen displacement q i being the controlling parameters. With the first of the two methods, subsequent points are the points of intersection between µ = η i and the equilibrium path, and with the second method-as the points of intersection between q k = η i with the same solution curve ( Figure 2). Both procedures are inefficient at the vicinity of local extremum of the controlling parameter. Riks [21,22] notes that at the limit points, the equilibrium path is tangent to the corresponding hyperplane: the angle θ between the tangent to the equilibrium path and the normal to the hyperplane is π/2. The intersection of the equilibrium path with this hyperplane is thus well defined when the angle θ is equal to or close to zero and poorly defined if it is close or equal to π/2. In this sense, the perfect family of surfaces will intersect the equilibrium path in perpendicular direction so that for every η the angle θ is zero ( Figure 3). An equation describing a "perfect" family of surfaces cannot be constructed because it is based on the equilibrium path which is unknown. However, an approximation can be used so that the angle θ is close to zero. Riks [21,22] proposed the following constraint equation: where: q-generalized coordinate vector; η − ηα-parameter approximating the arc length; the dots appearing above the symbols indicate derivatives with respect to the arc length. Equation (7) defines a hyperplane that is perpendicular to the vector (qα, µα) and distant from the point (qα, µα) by (η − ηα). It will intersect the equilibrium path almost perpendicular when the distance η − ηα is sufficiently small.
In this paper, the incremental-iterative method of the constant arc-length developed by Riks was applied to solve nonlinear equations. To compute the critical points, it is necessary to update the global tangent stiffness matrix at every iteration step. Four years after the arc-length method was proposed by Riks, Crisfield [23] came up with an alternative to that method. The latter was employed in many commercial FEM programs. Another variant of the constant arc-length method, similar to that developed by Crisfield, was put forward by Ramm [24].
In the review of arc-length methods, it is worth mentioning the pioneering work by Wempner [25]. The incremental procedures for tracing the equilibrium paths are still being developed in the works of Ritto-Correa et al. [26], Leon et al. [27], and Rezaiee et al. [28] etc. Riks [21,22] notes that at the limit points, the equilibrium path is tangent to the corresponding hyperplane: the angle θ between the tangent to the equilibrium path and the normal to the hyperplane is π/2. The intersection of the equilibrium path with this hyperplane is thus well defined when the angle θ is equal to or close to zero and poorly defined if it is close or equal to π/2. In this sense, the perfect family of surfaces will intersect the equilibrium path in perpendicular direction so that for every η the angle θ is zero ( Figure 3). Riks [21,22] notes that at the limit points, the equilibrium path is tangent to the corresponding hyperplane: the angle θ between the tangent to the equilibrium path and the normal to the hyperplane is π/2. The intersection of the equilibrium path with this hyperplane is thus well defined when the angle θ is equal to or close to zero and poorly defined if it is close or equal to π/2. In this sense, the perfect family of surfaces will intersect the equilibrium path in perpendicular direction so that for every η the angle θ is zero ( Figure 3). An equation describing a "perfect" family of surfaces cannot be constructed because it is based on the equilibrium path which is unknown. However, an approximation can be used so that the angle θ is close to zero. Riks [21,22] proposed the following constraint equation: where: q-generalized coordinate vector; η − ηα-parameter approximating the arc length; the dots appearing above the symbols indicate derivatives with respect to the arc length. Equation (7) defines a hyperplane that is perpendicular to the vector (qα, µα) and distant from the point (qα, µα) by (η − ηα). It will intersect the equilibrium path almost perpendicular when the distance η − ηα is sufficiently small.
In this paper, the incremental-iterative method of the constant arc-length developed by Riks was applied to solve nonlinear equations. To compute the critical points, it is necessary to update the global tangent stiffness matrix at every iteration step. Four years after the arc-length method was proposed by Riks, Crisfield [23] came up with an alternative to that method. The latter was employed in many commercial FEM programs. Another variant of the constant arc-length method, similar to that developed by Crisfield, was put forward by Ramm [24].
In the review of arc-length methods, it is worth mentioning the pioneering work by Wempner [25]. The incremental procedures for tracing the equilibrium paths are still being developed in the works of Ritto-Correa et al. [26], Leon et al. [27], and Rezaiee et al. [28] etc. An equation describing a "perfect" family of surfaces cannot be constructed because it is based on the equilibrium path which is unknown. However, an approximation can be used so that the angle θ is close to zero. Riks [21,22] proposed the following constraint equation: where: q-generalized coordinate vector; η − η α -parameter approximating the arc length; the dots appearing above the symbols indicate derivatives with respect to the arc length. Equation (7) defines a hyperplane that is perpendicular to the vector (q α , µ α ) and distant from the point (q α , µ α ) by (η − η α ). It will intersect the equilibrium path almost perpendicular when the distance η − η α is sufficiently small.
In this paper, the incremental-iterative method of the constant arc-length developed by Riks was applied to solve nonlinear equations. To compute the critical points, it is necessary to update the global tangent stiffness matrix at every iteration step. Four years after the arc-length method was proposed by Riks, Crisfield [23] came up with an alternative to that method. The latter was employed in many commercial FEM programs. Another variant of the constant arc-length method, similar to that developed by Crisfield, was put forward by Ramm [24].
In the review of arc-length methods, it is worth mentioning the pioneering work by Wempner [25]. The incremental procedures for tracing the equilibrium paths are still being developed in the works of Ritto-Correa et al. [26], Leon et al. [27], and Rezaiee et al. [28] etc.
The theory of stability deals with critical load and deformation of structures which are associated with sudden, quantitative changes of the structure state. The actual configurations have to be taken into account in the theory of stability of structure. This means that stability equations are actually nonlinear. They can be linearized in the case of linear bucking analysis.
In the paper, the dome is subjected to large displacement gradients and susceptible to stability loss from the condition of node snapping. Nonlinear geometrical relations are defined in the Lagrangian description. The constitutive relations are linear. The geometrical nonlinear analysis is performed using the "Autodesk Robot Structural Analysis" software (FEM code) (Autodesk Revit 2020, San Rafael, CA, USA).

Truss Element Description
In order to explain the effect of the height-to-span ratio on the selection of the appropriate calculation model and method of analysis of the bar dome, an example of the known von Mises truss was proposed. These issues are discussed in the pioneering work of Pecknold et al. [34].Additionally, notice that the former analytical solution has been extended to unsymmetric von Mises trusses by Ligarò and Valvo [35], to three-dimensional pyramidal trusses by Ligarò and Valvo [36], and to more general planar and space trusses by Rezaiee-Pajand et al. [37].
The truss is analyzed in the undeformed (initial) configuration 0 C and in the deformed configuration t C and t+∆t C ( Figure 4). The following material and geometric constants were adopted: E-Young's modulus; l 0 -the initial length; A 0 -the initial cross-sectional area; l-the length in the deformed configuration; A-the cross-sectional area in the deformed configuration. The deformation field in the deformed configuration t C is described by four components, q = {u 1 , v 1 , u 2 , v 2 }, and the displacement increment between time t and t+∆t is described by vector ∆q = {∆u 1 , ∆v 1 , ∆u 2 , ∆v 2 }. These displacement vectors correspond to nodal force vectors: Q = {U 1 , V 1 , U 2 , V 2 }, and ∆Q = {∆U 1 , ∆V 1 , ∆U 2 , ∆V 2 }. Displacement, strain, and stress fields in the deformed configuration t+∆t C are as follows: The displacement field u =[u v] T is interpolated by linear shape functions: The strain tensor in this case reduces to one nonzero quantity of the Green-Lagrange tensor that characterizes the elongation of a bar: where: L = 1 l 0 d dξ 0 -matrix of differential operators; t g = u v T -vector of displacement gradients, B 0 = L·N Equation (9) accounts for the effect of geometric nonlinearity, which means that large displacements were taken into account [17,29]. The Green-Lagrange strain tensor is coupled to the Piola-Kirchhoff second stress tensor. We use the initial reference configuration; thus, axial force S defined as S = EAε is not the real force in the bar. It corresponds to component S 11 = σ of the 2nd Piola-Kirchhoff symmetric stress tensor. The real axial force S' = S·l/l 0 . Displacement, strain, and stress fields in the deformed configuration t+∆t C are as follows: where: t q, t ε, t σ-displacements, strain, and stress at time t; ∆q, ∆ε, ∆σ-increments between time t and t + ∆t. Strain increment can be written as: where: Stress increments depend on the physical model adopted. For a linear-elastic material, Hooke's law-i.e., a linear relationship between stress and strain-is valid: Equilibrium equations of the element have the form: F e is the vector of internal forces in the element: Buildings 2020, 10, 179 For a plane truss element, these matrices take the form:

Frame Element Description
In the second example, the Schwedler dome consisting of 224 frame elements and 81 nodes was analyzed. The following is a description of the plane frame element analogous to the description of the truss element presented in chapter 2.2. Extension to space frame elements is possible [29] The frame element description was analyzed in the undeformed (initial) configuration 0 C and in the deformed configurations t C and t+∆t C ( Figure 5). The following material and geometric constants were adopted: E-Young's modulus; l 0 -the initial length; A 0 -the initial cross-sectional area; I 0 -the initial moment of inertia; l-the length in the deformed configuration; A-the cross-sectional area in the deformed configuration; I-moment of inertia in the deformed configuration. The vector of generalized stresses and modular matrix are: where: S-axial force; M-bending moment Stress increments depend on the physical model adopted. For a linear-elastic material, Hooke's law-i.e., a linear relationship between stress and strain-is valid: The increment of the strain vector with the dependencies Δu and Δg taken into account takes the form: The deformation field in the deformed configuration t C is described by six components, t q = {u 1 , v 1 , ϕ 1 , u 2 , v 2 , ϕ 2 }, and the displacement increment between time t and t + ∆t is described by vector ∆q = {∆u 1 , ∆v 1 , ∆ ϕ 1 , ∆u 2 , ∆v 2 , ∆ ϕ 2 }. These displacement vectors correspond to nodal force vectors: The displacement field is interpolated using the matrix of the shape function and the nodal displacement vector: It is visible that the approximation of displacement u, v by means of polynomials of the first and third orders, respectively, leads to Hermitian shape functions.
We assume that the Bernoulli-Euler hypothesis about the undeformable, perpendicular cross-section is valid. The strain at the point of cross-section is determined by the linear function of variable y: where: The vector of generalized stresses and modular matrix are: where: S-axial force; M-bending moment Stress increments depend on the physical model adopted. For a linear-elastic material, Hooke's law-i.e., a linear relationship between stress and strain-is valid: The increment of the strain vector with the dependencies ∆u and ∆g taken into account takes the form: where: The tangent stiffness matrix used in the incremental equilibrium equations of the frame element is presented as follows.

Von Mises Truss Analysis
Two cases of von Mises truss differing in the height-to-span ratio were considered. Figure 6 shows the structure geometry for both cases.  (Table 1).   (Table 1). The value of the determinant of the tangent det (K T ) stiffness matrix of the structure in the case of geometrically nonlinear analysis is 0.2851 × 10 9 , while in the linear analysis, 0.2872 × 10 9 . The value of the current stiffness parameter (CSP) practically did not change. In the linear analysis, this value was 1.0, while in the geometrically nonlinear analysis, 0.9929. The critical load multiplier obtained from the LBA analysis is µ cr = 423.638. According to the recommendations of the PN-EN-1993-1-1 [3] standard, for a critical multiplier value µ > 10, a linear analysis is sufficient. Figure 7a shows the equilibrium path for the von Mises truss case under consideration, while Figure 7b shows the CSP-q dependency graph. The values of the axial force and the vertical displacement of node 2, determined on the basis of LA and GNA, differ significantly (Table 2).  The values of the axial force and the vertical displacement of node 2, determined on the basis of LA and GNA, differ significantly (Table 2). The value of the determinant of the tangent det (K T ) stiffness matrix of the structure for geometrically nonlinear analysis is 0.4977 × 10 8 , while for linear analysis, 0.7999 × 10 8 . The value of the current stiffness parameter decreased significantly. In the linear analysis, this value was 1.0, while in the geometrically nonlinear analysis, 0.6223. The critical load multiplier obtained from the LBA is µ cr = 8.944. According to the recommendations of the PN-EN-1993-1-1 [3] standard, for a critical multiplier value µ < 10, a more accurate geometrically nonlinear structure analysis is recommended. The results obtained clearly confirm the entries in the Eurocode. In the second example, we will use them consistently. Figure 8a shows the equilibrium path for the von Mises truss case under consideration, while Figure 8b shows the CSP-q dependency graph. Based on the analyses carried out, a significant impact of structure geometry on the values of axial forces and displacements of the structure node was observed. In both cases, examples of the equilibrium path of von Mises trusses and the corresponding CSP-q path were presented.
In the first case, the structure with a height-to-span ratio of 0.25 was not susceptible to geometric nonlinearities. The difference in axial forces and displacements between linear and geometrically nonlinear analysis was below 1%. The change value of the stiffness matrix determinant also proved to be insignificant.
For shallow structures, changes in axial forces and displacements turned out to be significant. The axial force value increased from 100.125 kN obtained in a linear analysis to 116.6 kN in a geometrically nonlinear analysis.
Similarly, the values of the displacement of the second node changed from 2.236 to 2.795 cm. There was also a significant change in the value of the stiffness matrix determinant and the current Based on the analyses carried out, a significant impact of structure geometry on the values of axial forces and displacements of the structure node was observed. In both cases, examples of the equilibrium path of von Mises trusses and the corresponding CSP-q path were presented.
In the first case, the structure with a height-to-span ratio of 0.25 was not susceptible to geometric nonlinearities. The difference in axial forces and displacements between linear and geometrically nonlinear analysis was below 1%. The change value of the stiffness matrix determinant also proved to be insignificant.
For shallow structures, changes in axial forces and displacements turned out to be significant. The axial force value increased from 100.125 kN obtained in a linear analysis to 116.6 kN in a geometrically nonlinear analysis.
Similarly, the values of the displacement of the second node changed from 2.236 to 2.795 cm. There was also a significant change in the value of the stiffness matrix determinant and the current stiffness parameter.
A similar tendency in the behavior of the structure can be observed on fragments of larger structures, including shallow single-layer steel domes (Figure 9). Example 2 describes these problems.

Geometry and Loads
The dome with the diameter of 25 m and height of 1 m will top the building housing a restaurant. The location, a village of Święta Katarzyna, lies in the Świętokrzyskie Mountains.
The structure is based on 16 reinforced concrete columns, 5 m high. The whole structure is closed with a reinforced concrete wreath forming a rigid structure, which simulates grillages work.
A mesh of the dome consisting of 81 nodes and 224 elements is illustrated in Figure 10. The moststressed elements of the groups, and the nodes with the most displacement in GNA, are marked on Figure 10, case 1 by green and case 2 by red color. The structure of concern is Schwedler type bar dome. The length of each meridian is 2.511 m. Parallels become shorter as the dome tapers up; the successive lengths, starting at the base, are: 4.877, 3.908, 2.934, 1.958, and 0.979 m. The covering will comprise tubular steel bars and glass panels.
For the dome of concern, the following loads were assumed: permanent load-G (self-weight of the bar structure and the covering); snow load-S for zone 3 (acc. PN-EN 1991-1-3 [38]); wind load-W (acc. PN-EN 1991-1-4) [39]. The combination of loads applied to the structure was in compliance with PN-EN 1990 [40]. Evenly distributed loads were reduced to concentrated forces applied in all nodes. The load is concentrated to nodes according to the load distribution regions (envelope distribution).
To facilitate the description of analysis concerning a structure with a large number of nodes, the nodes located on subsequent rings were grouped. Table 3 shows the grouping of nodes into clusters R1-R5. Furthermore, loads applied to a group of nodes are given, which means each node of a given cluster is under the same load. The keystone (node no 1) does not belong to any of the groups. Table 3. Groups of nodes and the corresponding node numbers.

Group Name
Numbers of Nodes

Geometry and Loads
The dome with the diameter of 25 m and height of 1 m will top the building housing a restaurant. The location, a village ofŚwięta Katarzyna, lies in theŚwiętokrzyskie Mountains.
The structure is based on 16 reinforced concrete columns, 5 m high. The whole structure is closed with a reinforced concrete wreath forming a rigid structure, which simulates grillages work.
A mesh of the dome consisting of 81 nodes and 224 elements is illustrated in Figure 10. The most-stressed elements of the groups, and the nodes with the most displacement in GNA, are marked on Figure 10, case 1 by green and case 2 by red color. The structure of concern is Schwedler type bar dome. The length of each meridian is 2.511 m. Parallels become shorter as the dome tapers up; the successive lengths, starting at the base, are: 4.877, 3.908, 2.934, 1.958, and 0.979 m. The covering will comprise tubular steel bars and glass panels. The average characteristic value of the dome's own weight was assumed to be 1.45 kN/m 2 . The permanent load values for individual groups of nodes are summarized in Table 4. The average characteristic value of the snow load was determined based on the standard [38]. The structure of concern is located in the Świętokrzyskie Mountains at the height of A = 530 m amsl, which is the snow load in zone 3. For the considered structure, the characteristic value of roof snow load is: s = 1.651 kN/m 2 . The snow load values for individual groups of nodes are summarized in Table 5. For low-rise structures (flat roofs) for which the f/d ratio < 0.05, wind load is usually neglected. This is due to the common belief that consideration of wind has a positive effect on the values of cross-sectional forces of the bars. In the considered work, it was decided to present an example scheme for collecting wind load on a shallow bar dome and to verify the impact of uneven wind suction loading on structural behavior. For the dome of concern, the following loads were assumed: permanent load-G (self-weight of the bar structure and the covering); snow load-S for zone 3 (acc. PN-EN 1991-1-3 [38]); wind load-W (acc. PN-EN 1991-1-4) [39]. The combination of loads applied to the structure was in compliance with PN-EN 1990 [40]. Evenly distributed loads were reduced to concentrated forces applied in all nodes. The load is concentrated to nodes according to the load distribution regions (envelope distribution).
To facilitate the description of analysis concerning a structure with a large number of nodes, the nodes located on subsequent rings were grouped. Table 3 shows the grouping of nodes into clusters R1-R5. Furthermore, loads applied to a group of nodes are given, which means each node of a given cluster is under the same load. The keystone (node no 1) does not belong to any of the groups. The average characteristic value of the dome's own weight was assumed to be 1.45 kN/m 2 . The permanent load values for individual groups of nodes are summarized in Table 4. The average characteristic value of the snow load was determined based on the standard [38]. The structure of concern is located in theŚwiętokrzyskie Mountains at the height of A = 530 m amsl, which is the snow load in zone 3. For the considered structure, the characteristic value of roof snow load is: s = 1.651 kN/m 2 . The snow load values for individual groups of nodes are summarized in Table 5. For low-rise structures (flat roofs) for which the f/d ratio < 0.05, wind load is usually neglected. This is due to the common belief that consideration of wind has a positive effect on the values of cross-sectional forces of the bars. In the considered work, it was decided to present an example scheme for collecting wind load on a shallow bar dome and to verify the impact of uneven wind suction loading on structural behavior.
Unlike forces of dead load and snow load, the wind force acts perpendicular to the roof surface. As the structure of concern is a double-curved bar dome, it is necessary to determine angles of inclination of individual meridians-α n -in order to specify the wind force deviation from the vertical axis "z"-δ n . The values of angles α n and δ n are shown in Figure 11. The dome is divided by meridians into 16 equal parts. The angle γ by which the wind force in the XY plane should be rotated is γ = 360/16 = 22.5 • .
Unlike forces of dead load and snow load, the wind force acts perpendicular to the roof surface. As the structure of concern is a double-curved bar dome, it is necessary to determine angles of inclination of individual meridians-αn-in order to specify the wind force deviation from the vertical axis "z"-δn. The values of angles αn and δn are shown in Figure 11. The dome is divided by meridians into 16 equal parts. The angle γ by which the wind force in the XY plane should be rotated is γ = 360/16 = 22.5°. When calculating wind pressure coefficient, the standards applied were in accordance with [39]. The wind load distribution recommended by the standard is based on the assumption that on the dome arch determined by cutting the dome with a vertical plane perpendicular to the direction of the wind, the pressure is constant (broken lines in Figure 12). To determine the external pressure coefficients for the dome of concern, it is necessary to determine the following Equations (25) and (26): Figure 11. The values of angles α n and δ n.
When calculating wind pressure coefficient, the standards applied were in accordance with [39]. The wind load distribution recommended by the standard is based on the assumption that on the dome arch determined by cutting the dome with a vertical plane perpendicular to the direction of the wind, the pressure is constant (broken lines in Figure 12).
To determine the external pressure coefficients for the dome of concern, it is necessary to determine the following Equations (25) and (26): For the h/d Equation (25), the C pe,A , C pe,B , and C pe,C coefficients can be determined by linear interpolation of the data value in the graph in Figure 13 as a function of the f/d Equation (26). When calculating wind pressure coefficient, the standards applied were in accordance with [39]. The wind load distribution recommended by the standard is based on the assumption that on the dome arch determined by cutting the dome with a vertical plane perpendicular to the direction of the wind, the pressure is constant (broken lines in Figure 12). To determine the external pressure coefficients for the dome of concern, it is necessary to determine the following Equations (25) and (26) For the h/d Equation (25), the Cpe,A, Cpe,B, and Cpe,C coefficients can be determined by linear interpolation of the data value in the graph in Figure 13 as a function of the f/d Equation (26).  Based on Figure 13, the relationship (27) necessary for interpolation of Cpe,i indices was determined: y n = y 2 − y 1 For values of: h/d = 0.20 and f/d = 0.04, the following were obtained: C pe,A = − 1.14; C pe,B = − 0.56; C pe,C = − 0.20.
The considered dome was divided into the wind load zones as shown in Figure 14. Based on Figure 13, the relationship (27) necessary for interpolation of C pe,i indices was determined: For values of: h/d = 0.20 and f/d = 0.04, the following were obtained: C pe,A = −1.14; C pe,B = −0.56; C pe,C = −0.20.
The considered dome was divided into the wind load zones as shown in Figure 14.
On the basis of Figure 15, the relationship 28, necessary for interpolation of the wind pressure coefficient at the designated points (Table 6), was determined.
The considered dome was divided into the wind load zones as shown in Figure 14.  On the basis of Figure 15, the relationship 28, necessary for interpolation of the wind pressure coefficient at the designated points (Table 6), was determined.  After determining the pressure coefficients, the pressure of the wind on the dome was determined acc. standard [39]. For domes, it is necessary to determine the reference height that is significant when determining the wind load. The reference height of the dome in question is: z � e =h+ The values of wind load in individual nodes are obtained from multiplication of the peak value of wind speed pressure, wind pressure coefficient, meridian length values, and length values of individual parallels. The results for the sequentially numbered nodes are summarized in Table 7.   After determining the pressure coefficients, the pressure of the wind on the dome was determined acc. standard [39]. For domes, it is necessary to determine the reference height that is significant when determining the wind load. The reference height of the dome in question is: z e = h + f 2 = 5 m + 1 2 m = 5.5 m the base wind speed is v b,0 = 25.04 m/s; while the exposure factor is: c e (z e ) = 1.99. The base value of the wind speed pressure was determined on the level q b = 0.392 kN/m 2 . The peak value of pressure wind speed is equal to q p (z e ) = 0.392·1.99 = 0.780 kN/m 2 .
The values of wind load in individual nodes are obtained from multiplication of the peak value of wind speed pressure, wind pressure coefficient, meridian length values, and length values of individual parallels. The results for the sequentially numbered nodes are summarized in Table 7. In the first case, the situation of a standard approach to the dimensioning of shallow structures (flat roofs) was considered.
For the dome of concern, the following loads were assumed: permanent load-G (self-weight of the bar structure and the covering); snow load-S for zone 3 (acc. [38]).
The static-strength analysis, intended to produce the dimensioning of the steel elements, was performed acc. PN-EN-1993-1-1 [3] using Autodesk Robot Structure Professional 2019 software for spatial frame structure. To increase the accuracy of calculations, all the elements were additionally divided into 10 parts. The elements of the structure were assumed to be made of steel tubes with yield point f y = 235 MPa and Young's modulus E = 210 GPa; Poisson's ratio v = 0.3. In order to dimension the designed structure, three groups of rods were modeled: meridians, parallels, and diagonals.
The values of cross-sectional forces, load capacity for the most-stressed elements of these groups, and the displacement limits of node 79 are summarized in Table 8. The cross-sections taken for different groups of rods are summarized in Table 9. Figure 16 shows the deformation of the structure.
In the second step, the considered structure was subjected to linear buckling analysis. The lowermost critical load multiplier is equal to µ cr = 2.11199. Figure 17 shows the buckling modes for four successive eigenvalues. According to PN EN 1993-1-1/5.2.2 (5)B [3], when the lowermost critical load multiplier µ cr < 3.0, it is necessary to carry out a more accurate second order (GNA) analysis for the structure.
In the last step, the geometrically nonlinear static analysis was carried out. Nonlinear analysis produced the results of displacements and internal forces while taking into account second-order effects. In GNA, the lowermost critical load multiplier is equal to µ cr = 1.562. Exceeding this multiplier results in a peculiar stiffness matrix. GNA is a mathematical tool with which it is possible to generate complete equilibrium paths. Tracking the equilibrium path is inherently related to the analysis of post-critical states. The problem presented in example 2 concerns the design of a real shallow steel lattice structure. These types of structures are susceptible to large displacement gradients, which can only be considered by using geometrically nonlinear relationships. GNA in this case is used to determine the internal forces, displacements, and critical load multiplier. In the case of the analyzed shallow lattice dome, the critical load is related to the phenomenon of the node snap-through. After the critical load is exceeded, only unstable equilibrium states occur. From an engineering point of view, the part of the equilibrium path after reaching the critical load multiplier does not cause significant changes to the structure design process. Figure 18 shows the equilibrium path of the considered bar dome for case 1. In the first case, the situation of a standard approach to the dimensioning of shallow structures (flat roofs) was considered.
For the dome of concern, the following loads were assumed: permanent load-G (self-weight of the bar structure and the covering); snow load-S for zone 3 (acc. [38]).
The static-strength analysis, intended to produce the dimensioning of the steel elements, was performed acc. PN-EN-1993-1-1 [3] using Autodesk Robot Structure Professional 2019 software for spatial frame structure. To increase the accuracy of calculations, all the elements were additionally divided into 10 parts. The elements of the structure were assumed to be made of steel tubes with yield point fy = 235 MPa and Young's modulus E = 210 GPa; Poisson's ratio v = 0.3. In order to dimension the designed structure, three groups of rods were modeled: meridians, parallels, and diagonals.
The values of cross-sectional forces, load capacity for the most-stressed elements of these groups, and the displacement limits of node 79 are summarized in Table 8. The cross-sections taken for different groups of rods are summarized in Table 9. Figure 16 shows the deformation of the structure.   In the second step, the considered structure was subjected to linear buckling analysis. The lowermost critical load multiplier is equal to µcr = 2.11199. Figure 17 shows the buckling modes for four successive eigenvalues. According to PN EN 1993-1-1/5.2.2 (5)B [3], when the lowermost critical load multiplier µcr < 3.0, it is necessary to carry out a more accurate second order (GNA) analysis for the structure. In the last step, the geometrically nonlinear static analysis was carried out. Nonlinear analysis produced the results of displacements and internal forces while taking into account second-order effects. In GNA, the lowermost critical load multiplier is equal to µcr = 1.562. Exceeding this multiplier results in a peculiar stiffness matrix. GNA is a mathematical tool with which it is possible to generate complete equilibrium paths. Tracking the equilibrium path is inherently related to the analysis of post-critical states. The problem presented in example 2 concerns the design of a real shallow steel lattice structure. These types of structures are susceptible to large displacement gradients, which can only be considered by using geometrically nonlinear relationships. GNA in this case is used to determine the internal forces, displacements, and critical load multiplier. In the case of the analyzed shallow lattice dome, the critical load is related to the phenomenon of the node snap-through. After the critical load is exceeded, only unstable equilibrium states occur. From an engineering point of view, the part of the equilibrium path after reaching the critical load multiplier does not cause significant changes to the structure design process. Figure 18 shows the equilibrium path of the considered bar dome for case 1. Values of the internal forces and the maximum horizontal and vertical displacement of node 54, for the load multiplier µ equal to 1, are presented in Table 10.  Values of the internal forces and the maximum horizontal and vertical displacement of node 54, for the load multiplier µ equal to 1, are presented in Table 10. In the second case, a situation was considered in which the load caused by uneven wind suction was additionally taken into account. Therefore, the considered structure was subjected to: permanent load (G) (self-weight of the steel structure and roofing) and snow load (S) for area 3 of snow load according to PN-EN-1991-1-3 [38] and uneven wind suction load (W) determined on the basis of PN-EN-1991-1-4 [39]. Cross sections of individual groups of members were adopted according to Table 9.
The values of cross-sectional forces, load capacity for the most-stressed elements of these groups, and the displacement limits of node 9 are summarized in Table 11. Figure 19 shows the deformation of the structure. load (G) (self-weight of the steel structure and roofing) and snow load (S) for area 3 of snow load according to PN-EN-1991-1-3 [38] and uneven wind suction load (W) determined on the basis of PN-EN-1991-1-4 [39]. Cross sections of individual groups of members were adopted according to Table  9. The values of cross-sectional forces, load capacity for the most-stressed elements of these groups, and the displacement limits of node 9 are summarized in Table 11. Figure 19 shows the deformation of the structure.   In the second step, the considered structure was subjected to linear buckling analysis. The lowermost critical load multiplier is equal to µ cr = 2.35506. Figure 20 shows the buckling modes for four successive eigenvalues.

Allowable vertical displacement
(mm)-D/300 83.33 Maximum horizontal displacement (mm)-for node 9 4.52 Allowable horizontal displacement (mm)-H/150 6.67 In the second step, the considered structure was subjected to linear buckling analysis. The lowermost critical load multiplier is equal to µcr = 2.35506. Figure 20 shows the buckling modes for four successive eigenvalues. In the last step, the geometrically nonlinear static analysis was carried out. The critical load multiplier is equal to µcr = 1.390. Figure 21 shows the equilibrium path of the considered bar dome for case 2. In the last step, the geometrically nonlinear static analysis was carried out. The critical load multiplier is equal to µ cr = 1.390. Figure 21 shows the equilibrium path of the considered bar dome for case 2. In the last step, the geometrically nonlinear static analysis was carried out. The critical load multiplier is equal to µcr = 1.390. Figure 21 shows the equilibrium path of the considered bar dome for case 2.  Values of the internal forces and the maximum horizontal and vertical displacement of node 9, for the load multiplier equal to 1, are presented in Table 12. The huge impact of uneven wind load on its load bearing capacity in the structure under consideration is clearly visible. During initial verification of cross-sections in linear analysis, an increase in utilization of the elements from the "diagonals" group to 49%, compared to utilization of 8% in case 1, was noticed. Changes in element utilization in the linear analysis for the other groups did not differ significantly and amounted: from 67% to 65% for meridians and from 65% to 67% for parallels, respectively. Larger differences were observed in geometric nonlinear analysis, and for meridians, the utilization increased from 81% to 83%; in the case of parallels, from 85% to 92%; and in the case of diagonals, from 9% to as much as 90%. Differences in cross-sectional forces were also noted. In the linear analysis, the value of the axial force for meridians decreased from 563.694 to 518.05 kN; for parallels, increased from 162.79 to 164.801 kN; and for diagonals, increased from 1.977 to as much as 19.766 kN. Changes in the values of cross-sectional forces follow a similar trend in the case of GNA. For meridians, a decrease from 596.243 to 555.081 kN; for parallels, an increase from 208.020 to 214.545 kN; and for diagonals, an increase from 1.61 to as much as 36.408 kN were observed.
The impact of uneven load is also noticeable on the value of the critical load multiplier determined in LBA. For case 1, it is µ cr = 2.11199, while in case 2 it increases to µ cr = 2.35506. Both values indicate that a second order analysis is necessary. The critical load multiplier values for geometrically nonlinear analysis are, respectively, for case 1, µ cr = 1.562 and for case 2, µ cr = 1.390.

Conclusions
Linear analysis with linear bifurcation analysis works well for the design of steel high-rise single-layer coverings. The situation changes radically when the structure is a low-rise single-layer bar dome. These structures are subjected to large displacement gradients, and the actual configuration has to be taken into account in analysis. During the analysis, a clear decrease in measures describing the stiffness of the structure was observed. This structures should be designed according to geometrically nonlinear analysis.
The nonlinear analysis of bar structures with the finite element method (FEM) involves solving large systems of nonlinear equations. The relation between displacements and loads in a (N+1)-dimensional space is described by the load-displacement path. Stability analysis is closely related to singularities along this path. Controlling the incremental-iterative process to determine the equilibrium path and determine the location of critical points requires the establishment of a number of parameters. Their selection is not always intuitive and is time consuming. It often requires a lot of experience, but it is beneficial. The values of internal forces and displacements determined with two analyses mentioned above show very large differences. This concerns both the ultimate limit state and the serviceability limit state.
The paper focuses on comparing two load cases of low-rise bar domes, classic in accordance with Eurocodes requirements and taking into account uneven wind suction. Based on the analyses carried out, a significant impact of wind suction on the values of cross-sectional forces and displacements of structure nodes was observed. The use of stability analysis methods additionally enhanced the observed effect of the adverse effect of asymmetrical load on the structure. According to the authors, it is advisable to verify the load-bearing capacity of weak structures after taking into account the asymmetrical impact of wind on the structure.