ISOGEOMETRIC BEAM ELEMENT EXTENDED FOR GRIDSHELL ANALYSIS

Abstract. The design of gridshell structures is a complicated process, as the resulting shape of structure depends on initial undeformed grid geometry as well as on the history of boundary conditions used to form the final structure. In practice, the physical model is often used to determine the shape of the structure and the initial grid at the same time. Introducing Isogeometric analysis into a design of gridshells simplifies the design process as the problem can be easily recalculated when initial grid or boundary conditions change and the resulting shape can be immediately illustrated. The presented paper discusses possibilities in isogeometric gridshell modeling and proposes possible solutions of identified problems.


Introduction
Gridshell is a light-weight structure gaining its strength from a doubly-curved geometry.Such a structure is typically made of wood and it is constructed from initially straight grid of laths.The laths are mutually connected by joints allowing only for the in-plane rotation.The final shape of the structure is reached by subsequent introduction of a load or a prescribed displacement at selected nodes.Once the desired shape is reached, the joints are tightened in order to fix the in-plane rotation and stiffen the structure.See Fig. 1 for an example of a gridshell structure.
It is important to design the structural shape properly to reach high low-bearing capacity.The design process of a gridshell can be divided into three main parts: form-finding, analysis, and construction.Formfinding is typically performed by an architect who is seeking for the optimal shape of the gridshell taking into the account the aesthetics and structural behavior of the structure.Either physical modeling or computer simulation is used for form-finding, nevertheless these methods (and available tools) do not provide knowledge about structural behavior.When the desired shape is defined the structure is assessed in a Finite Element Method (FEM) solver, however the analysis considers only the final shape of the structure.There are no available tools which would take into account the history of the construction steps and would allow to assess all the intermediate construction stages.Currently the construction itself relies on trial-error methods supported by a human intuition and experience.
Our goal is to develop a complex tool for analysis assisted form-finding which would allow to find the optimal shape of gridshell and which would help to propose sequence of construction steps.The idea is to use the isogeometric approach which eases incorporation of the analysis within available design tools and which may bring higher accuracy into the analysis of doubly-curved geometry thanks to the more suitable geometry description.The presented paper discusses the advantages of isogeometric gridshell modeling and also points out some issues connected to its application.

Isogeometric analysis of beams
Isogeometric analysis has been introduced by Thomas Hughes in 2005 [2] with the aim to bridge the gap between the Computer-Aided Design (CAD) and the Finite Element Analysis (FEA).This goal has been achieved by using the same CAD basis functions for geometry description and unknowns approximation in the analysis.Thus one model can be shared between the design and the analysis.

NURBS-based analysis
The most wide-spread technology in CAD industry are NURBS (Non-Uniform Rational B-Splines).The NURBS curve is obtained by linear combination of Cartesian coordinates of the control points P using NURBS functions Each NURBS function is generated by weighting Bspline functions N p i with a given weight w i associated with control point B-spline functions form a special subset of NURBS functions (corresponding to constant w i ) and are derived recursively starting with a piecewise constant functions where p is the degree of the B-spline function, ξ i is the coordinate of the i th -knot and ξ ∈ 0, 1 is a parametric coordinate.The computational domain in Isogeometric analysis (IGA) is firstly divided into patches, which are further divided into knotspans (often referred to as elements in IGA).Understanding of knotspans is similar to elements in standard FEM, nevertheless higher continuity up to C p−1 between knotspans can be achieved naturally using NURBS, unlike C 0 in standard FEM.Moreover, the NURBS basis functions are generally non-iterpolatory.In the knots (points which are dividing patch into knotspans), the continuity can be locally decreased up to C 0 by knot multiplication, which is a standard technique in NURBS technology.This procedure imposes the interpolatory behavior of the functions at the particular point.
For better understanding of NURBS see Fig. 2, where mapping between real and parametric geometry is illustrated.For more information the reader could refer to [3].

Beam element formulation
The formulation of the three-dimensional beam element [4] is derived in local coordinate system (t, n, b), where t, n, b are tangent, normal and binormal vectors and s is the curvilinear coordinate measured along the centerline of the beam.Presented formulation is based on Timoshenko beam theory and therefore, in three-dimensional space, there are six independent unknowns: tangential displacement u t , normal displacement u n , binormal displacement u b and rotations θ t , θ n , θ b (see Fig. 3).
A strain-displacement matrix B is defined as where Curvature κ is given as and torsion where r(s) is the position vector.The material stiffness matrix D is given by constitutive relations where E is the Young's modulus, G is the shear modulus, I k is the torsional moment, I n , I b are the moments of inertia, A is the cross-sectional area, A b , A n are the reduced shear areas.Finally, the element stiffness matrix K is defined as where the integral is evaluated using Gaussian quadrature, which is not exact for NURBS functions, but it has a sufficient accuracy.

Numerical locking
Well-known problem of numerical locking of Timoshenko beam elements in standard FEM occurs also for the isogeometric formulation.The same approximation order for displacements and for rotations leads to the shear strains of higher order than bending moments, while it should be vice-versa.Moreover the constant strains cannot be approximated along the entire patch due to the field inconsistency in the shear terms.
The reduced integration, popular locking-treatment in standard FEM for its computational efficiency, is not convenient for IGA, as the optimal reduced integration rule depends on the choice of approximation order and continuity and such a rule is hard to be determined.On the contrary, Discrete shear gap (DSG) method [5,6] or B-method [5] can be used independently of the choice of approximation, the latter being used in the presented implementation of the element.
The main idea of the B-method is to project membrane strain ε m and shear strains γ n , γ b onto a basis of lower order.For example, the contribution of membrane strain ε m to the B matrix is where for example and where N j is the j th basis function of the original basis and Ñi is the i th basis function of the lower order basis.The membrane contribution to stiffness matrix K m is where M is the element "mass matrix" calculated in The analogical procedure can be followed also for the shear strains γ n and γ b .

Application to gridshells
Isogeometric formulation of beam element allows to perform the analysis directly within the CAD system, however the modeling of gridshells with presented isogeometric beam element formulation is not straightforward.There are two possible approaches to modeling gridshells: a) multiple patches per one lath, b) single continuous patch per lath.

Lath modeled by multiple patches
This approach is similar to the use of standard FEM, where the joint is usually modeled by inserting a node.
In IGA, it means to assign at least one patch to each lath segment between two joints leading to the interpolatory degrees of freedom at particular points (joints).This approach is simple and does not require any additional implementation, nevertheless it leads to loss of some advantages of isogeometric formulation, namely loss of higher inter-element continuity along the entire lath and loss of possibility to handle one lath as a single entity, which might be highly convenient for the development of a complex design tool.A simple test has been performed to verify the "multiple-patches" approach.Two laths mutually connected in the middle, each modeled by two patches with one knotspan and cubic B-spline approximation, have been subjected to the transverse force placed at the joint.The model and resulting shear force γ n and bending moment χ b are shown in Fig. 4. The results overlap the exact solution and prove the possibility to use this modeling technique for gridshell analysis.

Lath modeled by a single patch
Isogeometric approach offers a possibility to model each grid lath by a single patch.Such a model requires to add additional constraints to enforce the compatibility of displacements and rotations in joints as the compatibility cannot be achieved by sharing particular degrees of freedom due to the non-interpolatory nature of isogeometric basis functions.A special boundary condition based on the Lagrange multipliers method has been implemented to support this approach.The potential energy is enhanced by an additional term yielding the required continuity where λ k is vector of Lagrange multipliers for k th joint, rk i and rk j are the vectors of constrained unknowns at the k th joints on the i th and j th beams, respectively.Note, that due to the non-interpolatory nature of NURBS basis functions, the vectors rk i and rk j have to be expressed as a linear combination of the control points values using NURBS basis functions evaluated at the joint position.The local coordinate ξ i (ξ j ) of the joint on the beam i (j) has to be calculated according to the parametrization.The physical meaning of Lagrange multipliers λ corresponds to the forces and moments required to enforce the compatibility and thus such an approach results in concentrated force and moment loadings at the joints positions.The problem of such a formulation is the inability of NURBS basis functions to represent exact solution with discontinuities in strains and internal forces corresponding to concentrated loadings.This inability is a consequence of the higher continuity of NURBS basis along the entire patch.
This problem is documented on the similar test as the one used in Section 3.1 (see Fig. 4).Unlike in previous section, each lath is now modeled by a single patch using just a single span with cubic Bspline approximation, which results in C 3 -continuity at the joint position causing the oscillations in internal forces, see Fig. 5a for the results.The continuity can be lowered to C 2 by placing a knot at the particular position, the continuity can be further reduced up to C 0 by increasing the knot multiplicity.The resulting internal forces for all choices of continuity are shown in Fig. 5.The oscillations in the numerical solutions are sufficiently removed only in the case of C 0 -continuity.Note, that the procedure of decreasing the continuity to C 0 results in the same model as in Sec.3.1 and thus the advantages of isogeometric analysis are eliminated.
The main idea, how to preserve the higher continuity while avoiding oscillations in numerical solutions, is to enrich the NURBS basis by the special basis functions tailored to match the exact solution.In this paper, the derivation of the additional basis functions is illustrated on the case of two-dimensional straight beam subjected to concentrated force.
From the knowledge of the exact solutions of bending moment M ex and shear force Q ex under the unit force where β is a rotation and v is a transverse displacement, the additional basis functions N * β , N * v can be derived See Fig. 6 for the resulting basis functions.The final approximation for the enhanced solution is The integration constants were determined with respect to the boundary conditions of a simply supported beam.The performance of the enhanced approximation has been tested by means of two simple test.Firstly, simply supported beam subjected to the concentrated force load was considered, in the second test the concentrated force has been combined with continuous loading.The resulting internal forces for both load cases overlap exact solutions, see Fig. 7.

Conclusions
Two possible approaches to gridshell modeling have been presented.The simple FEM-like approach, where the joints are modeled by connection of two patches, provides good results, however the advantages of isogeometric approach are lost.In order to keep the higher continuity along the entire lath, it has to be modeled by a single patch.In such a case, the enforcement of continuity of displacements and rotations in joints by Lagrange multipliers inevitably leads to application of concentrated forces and moments.It has been shown, that enhancement of NURBS basis functions by a set of the special functions is sufficient to support concentrated loadings.
The paper documents the first necessary steps of the development of a complex design tool for gridshells.The presented idea of enhancement of the NURBS basis has to be combined with locking treatment to provide robust element formulation, moreover due to the large deformations, the formulation has to be extended for geometrically non-linear analysis.

Figure 2 .
Figure 2. Description of B-spline (NURBS) finite element geometry.The geometry is modeled by one patch consisting of two knotspans with quadratic NURBS approximation.

Figure 3 .
Figure 3. Local coordinate system and degrees of freedom of a three-dimensional beam element.

Figure 4 .
Figure 4. Simple gridshell structure modeled by four patches with cubic approximation.Joint is modeled by connection of patches.Resulting internal forces overlap exact solution.

Figure 5 .
Figure 5. Resulting internal forces (blue continuous line) and exact solution (red dash-dot line) for different continuity at the position of concentrated force load.

Figure 6 .
Figure 6.Additional basis functions to support case of concentrated force load on 2D straight beam with NURBS approximation.

Figure 7 .
Figure 7. Resulting internal forces for simply supported beam subjected to a) a concentrated force, b) a concentrated force and continuous force loading.NURBS approximation enriched by the functions tailored to support case of concentrated force has been used.The results overlap the exact solution.