Isogeometric analysis applied to 2D Bernoulli-Euler beam model: imposition of constraints by Lagrange and Penalty methods

The isogeometric analysis (IGA) consists of using the same shape functions, usually employed on Computer-Aided Design (CAD) technologies, on both geometric modelling and approximation of the fields of physical models. One issue that concerns IGA is how to make the connection or apply general constraints in the connection of structures described by different curves and surfaces (multi-patch structures), particularly when the shape functions are not interpolatory at the selected point for the imposition of the constraint or the desired constraint is not related directly to degrees of freedom, which may be an issue on Kirchhoff-Love shells and Euler-Bernoulli beams, since usually no rotational degrees of freedom are employed. In this context, the present contribution presents an isogeometric 2D curved beam formulation based on Bernoulli-Euler assumptions. An approach about the implementation of multi-patch structures enforcing constraints, such as same displacement or same rotation among neighbor paths, is developed based on Penalty and Lagrange methods. The applicability of the methods is verified by examples of application.

The paper starts with a brief overview about B-Splines and NURBS, showing the main aspects of the CAD technologies and their relation with finite element related analysis. Then, the 2D isogeometric beam formulation is presented. Subsequently, the Penalty and Lagrange Methods approaches are developed. Finally, examples of application are presented.

B-SPLINES AND NURBS
This section presents briefly the main concepts for self-containing of herein developed contributions. Further information about NURBS and isogeometric analysis can be found on [20] and [1], respectively.
Knot vector is a set of non-decreasing ordinates of the parametric space from which B-Spline basis functions are generated. It is a vector of the form Ξ = � 1 , 2 , … , + +1 �, where ∈ ℝ is the i-th knot, n is the number of basis functions, and p is the polynomial order of the B-spline basis.
When the knots are equally spaced, the knot vector is called uniform, otherwise, the knot vector is called nonuniform. The knots also may be repeated, and when the first and last knot appears + 1 times, the knot vector is said to be open.
An open knot vector generates B-spline basis functions that are interpolatory at ends of the parametric space interval � 1 , + +1 �, that is, at the first end of the interval, 1 , 1, is unitary while the other basis functions are null; and at the last end of the interval, + +1 , , is unitary while the other basis functions are null. This feature makes the basis functions analogous to the Lagrange polynomials used in the traditional Finite Element Method, which are interpolatory at the ends of the element.
B-spline basis functions are generated recursively according to a knot vector, starting from piecewise constant polynomials: An interesting aspect about B-splines lies on the continuity of the basis functions. If there are no repeated knots, B-spline basis functions of order have ( − 1) continuous derivatives. If an internal knot appears times, the number of continuous derivatives at the knot is equal to ( − ). When the multiplicity of an internal knot coincides with , the basis functions are interpolatory at the knot. Figure 1 illustrates a second-order B-spline basis created by a uniform knot vector, Ξ 1 = {0,1,2,3,4,5,6}. Since there are no repeated knots, the basis functions have the first derivative continuum, which is indicated by the smoothness of their shapes. While Figure 2 exemplifies a second-order B-spline basis generated by an open knot vector, Ξ 2 = {0,0,0,1,2,3,3,4,5,6,6,6}, which also is composed by an internal knot (3) with multiplicity equal to the polynomial order. So, Basis functions are interpolatory at ends of the parametric space, 1 and 6, and at the repeated knot 3. 3) Non-negativity: , ≥ 0, ∀ .
B-spline curves, in turn, are generated by a linear combination of control points in ℝ 2 or ℝ 3 ( ) employing B-spline basis functions such that: If there are no repeated knots or control points, the B-spline curve has ( − 1) continuous derivatives. If a knot or control point appears times, the number of continuous derivatives is equal to ( − ). Figure 3 illustrates the B-spline curve obtained through the combination between the second order basis functions from Figure  The dashed lines are called control polygon and correspond to the linear interpolation between the control points. Note that the curve first derivative is not continuous at control point 5 . This fact is due to the repetition of knot 3, whose multiplicity is responsible for the generation of the interpolatory basis functions at = 3 (see Figure 2).
The rational basis functions can be defined as: So, the NURBS curve is also given by: The rational basis functions have the same properties described above for the B-spline basis, such as, partition of unity, piecewise polynomials, non-negativity and continuity conditions related to the multiplicity of knots. Thus, if all of the weights are equal to unity( = 1, 1 ≤ ≤ ), the rational basis functions are equal to the B-spline basis functions. So, B-Splines are a particular case of NURBS.
One interesting aspect about NURBS, in detriment of B-splines, consists on the additional local shape control due to the weights. Both NURBS and B-splines can be shaped by changing the position of the control points, but NURBS also are composed by the weights, which control the proximity of the curve to the control points. It is important to state that any change in either a weight ( ) or a control point ( ) affects only the portion of the curve comprehended in the interval � , + +1 �, which is the support of the rational ( ) or B-spline basis function ( ) they are related to (see Eq. (3) and (4)). Figure 4 illustrates local shape control associated with the weight 2 . The second order NURBS curves are generated by the same open knot vector Ξ 3 = {0,0,0,1,1,1} and control points C 3 = �{1,0}, {1,1}, {0,1}�, and the initial curve, which is shown in black, is created by a set of unitary weights = {1,1,1}. So, the initial curve is also a B-spline curve. Decreasing the value of 2 to √2 2 ⁄ , the NURBS curve gets far from control point 2 , matching a quarter circumference, represented by the red dashed curve. Increasing the value of 2 to 1.5, the curve, drawn in dashed blue, approaches the control point 2 . So, one can observe that increasing the value of weight, the curve approaches the respective control point, and vice-versa. The aim of the isogeometric analysis is the approximation of the element displacements using the same basis functions employed on the construction of the geometry, that is, the rational or B-splines basis functions (in the present work). Sometimes, even when the basis functions are able to perfectly construct the geometry, they might be unsatisfactory to approximate displacements, yielding bad results. However, there are some kinds of refinements in which the basis functions can be enriched, without changing the curve geometrically or parametrically, and, so, the results for the displacements can be improved. In isogeometric analysis, there are two types of refinements that are analogous to the traditional finite element method: h-and p-refinements. In h-refinement, knots are inserted into the knot vector. For each knot inserted into the knot vector, one basis function is created. So, inserting knots into the knot vector leads to a more numerous basis, which might conduct to a better approximation. In p-refinement, the order of the basis functions are elevated, the functions get more complex and the approximation is improved. There is one more refinement procedure that can be applied in isogeometric analysis that is called k-refinement, which consists of first order elevating the basis functions (p-refinement) and then inserting knots to the knot vector (h-refinement). The k-refinement takes advantage of the non-commutativity of order elevating and knot insertion to construct a higher continuity basis. Further information about refinements can be found at [1] and [20].

BERNOULLI-EULER CURVED BEAM MODEL
In this section, a curved beam model based on geometrically-linear Bernoulli-Euler kinematic assumptions is presented. The beam axis, which is described by a 2-D B-spline or NURBS curve, passes through the centroid of the successive constant cross sections. As established on Bernoulli-Euler kinematic assumptions, the cross sections remain plane and orthogonal to the deformed beam axis. So, kinematics assumptions leads to neglecting the shear strain energy. Figure 5 illustrates the curved beam model. Coordinate describes the beam axis, and corresponds to the tangent and transverse local axes, respectively, to the beam axis, 1 and 2 are the global axes, and corresponds to the radial and tangential loads, respectively.
= ( + ) + ; Then, the longitudinal strain of a point at an arbitrary distance from the beam axis is related to the change of the beam curvature ( ), as follows: Adopting a linear elastic constitutive model, in which is the Young modulus, the longitudinal stress of the beam is given by: As the beam axis is curved, the tangential displacement affects the rotation of the cross sections, which is represented in the Figure 7. So, the change of the beam curvature is given by: where is the rotation of the cross sections.
Where and are the internal and external virtual works, respectively. The symbol is used, in the present work, in front of variables to denote that they are virtual quantities.
The internal virtual work is given by: where V is beam volume. One may develop: where L is the beam length and A is the cross-sectional area. Since the beam axis passes through the centroid of cross sections: So, Eq. (16) reduces to: Eq. (19) also can be written as the sum of two contributions, as follows: where: The external virtual work, in turn, is given by the following equation: where and are the radial and tangential loads distributed along the beam axis, respectively, , and are the radial and tangential forces and the moment, respectively, applied at an arbitrary coordinate value 0 of the beam axis.
As the geometry of beam axis is described, in this model, by a B-spline or a NURBS curve, transverse and tangential displacements are approximated by a linear combination between B-spline or rational basis functions, { 1 , … , }, and the unknown parameters, { 1 , … , } and { 1 , … , }. Thus, Galerkin Theorem is applied, that is, the virtual displacements, and , are approximated by the same functions used on the approximation of the real displacements, { 1 , … , }. The approximation for the transverse displacement is showed below, the other ones can be found analogously: Rearranging Eq. (24) using vectors, we have: Where: Being the Jacobian of transformation between coordinates and , we can find the following equations of the first and second derivatives of in relation to : Latin American Journal of Solids and Structures, 2020, 17(1), e248 9/27 = , (28) Rewriting Eq. (30) and (31) using Eq. (26) and (27), we have: where: ( )′ and ( )′′ are the first and second derivatives with respect to , respectively. Then, substituting the approximations for the real and virtual displacements in Eq. (21), (22) and (23), we have: where: and 1 , 2 and 3 are the ordinates from the parametric space in which the radial and tangential concentrated forces and concentrated moment are applied, respectively. The suppressed calculations are detailed in section "Development of the equations of the Bernoulli-Euler beam model" of the Appendix A.
So, we find an analogous form to the traditional finite element method, since corresponds to the stiffness matrix, is the vector of the displacements and is the vector of external forces.
In addition, according to [14], the Jacobian of the transformation between coordinates and , and the radius of curvature are given by: where 1 and 2 correspond to the projections of the curve on the axis 1 and 2 , respectively.

CONSTRAINTS
When the independent variables of the isogeometric Bernoulli-Euler beam solution patch, from now on named "element", and , are approximated by a basis generated by an open knot vector, the basis functions are interpolatory at the ends of the element, and it is not a hard task to impose displacement constraints on and at the ends of the element. Since the rotation of the cross sections is not an independent variable in the isogeometric Bernoulli-Euler beam element, that is, is function of the derivative of the transverse displacement with respect to and the tangential displacement , it is not so straightforward to impose the rotation-related constraints at the ends of the element. However, this problem can be solved by using the Penalty or the Lagrange method. In both methods, contributions are added into the equation of the PVW in order to impose the kinematic constraints, as follows: where is the contribution of the constraints to the model equations, which is calculated differently for each method: = − for Lagrange method, and = for penalty method. The approaches presented in this section may be used to impose displacement constraints on any point of the element (patch) and also establish the connection between different elements (patches).
For the development of the contributions to the model PVW of a given desired constraint, we assume a scenario involving two isogeometric Bernoulli-Euler beam elements, numbered 1 e 2, with the external constraints applied only to the element 1 and the internal constraints for the linking of both elements. Any other situation can be derived analogously.
According to Lagrange Method, the potential contribution �Π � to enforce a desired constraint given by � 1 − � = 0 may be written as follows: where is a scalar Lagrange parameter (or Lagrange multiplier) and 1 = is the desired constraint we want to enforce, where 1 is a given component of displacement for an element 1 and is the prescribed displacement.
The contribution of the constraint term to the model weak form (PVW) is obtained by the variation of Eq. (43): According to Lagrange Method, the potential contribution �Π � to enforce a desired constraint given by � 1 − 2 � = 0 may be written as follows: where is a scalar Lagrange parameter (or Lagrange multiplier) and 1 = 2 is the desired constraint we want to enforce, where 1 is a given component of displacement for an element 1 and 2 is a given component of displacement for an element 2.
The contribution to the model weak form (PVW) is obtained by variation of Eq. (45): The external constrained displacements can be written in terms of their approximations by B-splines or rational basis functions (see Eq. (25) and (30)), as follows: , are the coordinates from parametric space that generate the basis functions of element 1 in which the external constraints are applied on 1 , 1 and 1 , respectively, is the vector that contains the basis functions of element 1, and are the vectors of the unknowns of element 1, 1 and 1 are the radius of curvature and the Jacobian of element 1, respectively.
Analogously, the internal constrained displacements can also be written in terms of their approximations by B-splines or rational basis functions, as follows: Where: The suppressed calculations can be found in "Development of the equations of the Lagrange Method" of the Appendix B.
So, according to Eq. (56), the contributions of the internal and external constraints calculated through Lagrange method are added to the stiffness matrix of the structural model constructed with the isogeometric Bernoulli-Euler beam elements, and the structural analysis can be accomplished by solving the forthcoming linear system, already encompassing the original degrees of freedom, but now also containing the Lagrange multipliers as unknowns.
The approach for the Penalty method starts in the same way as done for Lagrange method, that is, it begins with the potentials of the constraints: where Π is the potential of the external constraint, Π is the potential of the internal constraint, and are the penalty constants, which corresponds to real arbitrary values, usually taken as big as possible, in order to better impose the desired constraint.
The contribution to the model weak form (PVW) is obtained by variation of Eq. (69) and (70): where δW and δW are the virtual works of the external and internal constraints, respectively. Substituting the approximations for the external constrained displacements, Eq. (47) to (49), into Eq. (71) leads to: Where: The suppressed calculations can be found in "Development of the equations of the Penalty Method" of the Appendix C.
Analogously, the approximations for the internal constrained displacements, Eq. (50) to (55), are substituted into Eq. (72), as follows: Where: The suppressed calculations can be found in "Development of the equations of the Penalty Method" of the Appendix C.
Remembering Eq. (37), (38) and (42), Eq. (73) to (87) can be organized to compose the following global system: Where: So, through the approaches based on the Penalty and the Lagrange method presented in this section, constraints can be imposed to the isogeometric Bernoulli-Euler beam elements in order to compose a great sort of scenarios involving 2-D curved beams.

NUMERICAL EXAMPLES
In this section, we show examples to demonstrate the proposed isogeometric beam formulation and the imposition of the constraints using the Penalty and the Lagrange methods. All the numerical results here shown were obtained by an implementation using Mathematica™.
In the numerical examples, the integrations were made using the standard Mathematica´s numerical integration function, which usually uses adaptive algorithms that recursively subdivide the integration region as needed. The method and rule options were set to "automatic". The block of zeros in the stiffness matrix produced by the Lagrange Method (for cases with constraint imposition), by its turn, was handled by a linear solver standard function from Mathematica´s library, which can work even with sparse matrices.
1. Quarter circumference arch and Penalty methods The first example consists of a cantilever quarter circumference arch with a vertical force applied at the free end ( Figure 8). The problem data are listed on Table 1. The geometry is described by a fourth order NURBS curve, whose knot vector is open and corresponds to {0,0,0,0,0,0.2,0.4,0.6,0.8,1,1,1,1,1}, and the weights and control points are organized on Table 2. As the problem is simple, the external constraints can be imposed exploring rational basis functions properties on the following equations: Since the knot vector is open, basis functions are interpolatory at the ends of the parametric space, and Eq. (93) and (94) can be simplified as: Another feature of the basis functions generated by the open knot vector is that, at = 0, only 1 ′ and 2 ′ are non-null. Thus, Eq. (95) is simplified as: So, the external constraints of this problem can be applied by only making 1 , 2 and 1 equal to zero. Then, the obtained linear system is solved, and the result for the transverse displacement at the free end is given by: According to [21], the analytical solution for this problem is given by: A 0.49% error is verified. So, we conclude that the problem is well represented by the isogeometric Bernoulli-Euler beam element and the approximations for the displacements by means of nine rational functions, that is, nine degrees of freedom, are satisfactory for engineering purposes.
2. Tudor arch The second example corresponds to a structure named Tudor arch, which is composed, geometrically, by two circular arches and two straight lines. Figure 9 illustrates the geometry and the loads applied on the structure, and the problem data are listed on Table 3.   The structure's computational model is composed by four isogeometric elements: two circular arches and two straight bars. At the ends of the arches, there are two external constraints, one on w and the other one on u, for each arch. The arches are linked to the straight bars by three internal constraints, on , and . The straight bars, in turn, are connected only by translational displacements and , the continuity of the rotations are not imposed between them.
For the approximation of the displacements of the straight beams, it was used a fourth order B-spline basis generated by the open knot vector {0,0,0,0,0,0.5,1,1,1,1,1}. As they are straight beams, their curvature radii are infinite, which simplifies some of terms in the model equations. For instance, for straight beams, Eq. (37) and (38) are simplified to: Firstly, the problem was solved using Lagrange method, and the resultant Lagrange multipliers are listed on Table 5. Then, the results obtained for the displacements were used to calibrate Penalty parameters. Table 6 shows a comparison between the resultant displacement � = � 1 2 + 2 2 � at the top of the structure obtained with the Lagrange � = 4.1372 × 10 −4 � and Penalty Method, varying penalty parameters in a range comprehended between 10 12 to 10 18 . One can observe that for 10 15 ≤ ≤ 10 16 , the solutions obtained through Penalty method are practically the same as the solution obtained using Lagrange method, considering four digits of precision.  The Lagrange multiplier can be also used as a hint about the magnitude of the penalty parameter to be employed. For instance, λ e u,1 is the Lagrange multiplier associated to the external constraint on the axial displacement u of the left arch, and it is also the axial force on the end of the arch. Based on the Lagrange multiplier and an expected violation on the constraint ∆ , an estimative of the penalty parameter is given by: In fact, the violation on the external constraint above mentioned using 1.0E+12 as the penalty parameter is about 1.6E-07, and the Lagrange multiplier is about 1.6E+05, so: The Lagrange multiplier may be understood as the force or moment required to impose a constraint, and, comparing the Lagrange multipliers from Table 6 with the diagrams from Figure 10, this fact is exemplified. For instance, ,1 and ,4 correspond to the axial forces on the ends of the arches (160.7 kN and 138.6 kN, respectively), and ,12 and ,34 are the bending moment on the linkage of the arches with the straight bars (29.5 kNm and 21.2 kNm, respectively), and so on.
3. Bridge application The last example corresponds to a 2D beam structure inspired on an arch bridge, which is composed, geometrically, by a circular arch and a straight line. Figure 11 illustrates the geometry and the loads applied on the structure, and the problem data are listed on Table 7.  The structure's computational model is composed by four isogeometric elements: two circular arches and two straight bars. At the ends of the arches, there are three external constraints, on , and , for each arch. The arches are joined themselves by three internal constraints, on , and . At the ends of the straight bars, there are two external constraints, on and , for each bar. The straight bars are joined themselves by three internal constraints, on , and . The arches are linked to the straight bars only by translational displacements, and , the continuity of the rotations are not imposed among them.
The geometry of the arches are given by fourth order NURBS curves, whose knot vector is open and corresponds to {0,0,0,0,0,0.2,0.4,0.6,0.8,1,1,1,1,1}, and the weights and control points are listed on Table 8. For the approximation of the displacements of the straight beams, it was used a fourth order B-spline basis generated by the open knot vector {0,0,0,0,0,0.5,1,1,1,1,1}. Firstly, the problem was solved using Lagrange method, and the resultant Lagrange multipliers are listed on Table  9. Then, the results obtained for the displacements were used to calibrate Penalty parameters. Table 10 shows a comparison between the resultant displacement � = � 1 2 + 2 2 � at the top of the arch obtained with the Lagrange � = 3.5402 × 10 −4 � and Penalty Method, varying penalty parameters in a range comprehended from 10 11 to 10 16 . One can observe that for 10 15 ≤ ≤ 10 16 , the solutions obtained through Penalty method are practically the same as the solution obtained using Lagrange method, considering four digits of precision.    The Lagrange multiplier may be understood as the force or moment required to impose a constraint, and, comparing the Lagrange multipliers from Table 9 with the diagrams from Figure 12, this fact is exemplified. For instance, ,1 and ,2 correspond to the axial forces on the ends of the arches (89.3 kN), and ,12 and ,34 are the bending moment on the linkage of the arches and the straight bars (28.9 kNm and 14.9 kNm, respectively), and so on.

CONCLUSION
In the present paper, a 2D isogeometric curved beam element based on Bernoulli-Euler assumption is presented. In order to establish connections and aiming the continuity of rotation between elements, two different approaches, according to Lagrange and Penalty methods, of imposing constraints on the beam elements are developed.
Three examples of application are presented. The first one corresponds to a single patch example and is dedicated to verify the concordance of the isogeometric curved beam element with analytical solution. The other ones are multi-patch and they are solved by using the Penalty and Lagrange methods to impose the displacements constraints. The results are compared with each other and the calibrating of the penalty parameters is discussed.
Attention should be given to some features about Lagrange and Penalty methods in order to not compromise the problems solutions. The Lagrange method brings new unknowns to the problem, which are the Lagrange multipliers, and also blocks of zeros to the structure global stiffness matrix. These blocks of zeros may be carefully handled, since the resulting linear system of equations may not be solved by a simple algorithm. The Penalty method, in turn, does not introduce unknows to the problem, but the precision of the solution is highly dependent of the arbitrary parameter β.

APPENDIX A: DEVELOPMENT OF THE EQUATIONS OF THE BERNOULLI-EULER CURVED BEAM MODEL
The development of the equations of the Bernoulli-Euler curved beam model is detailed in this appendix. Substituting the approximations for the real and virtual displacements (see Eq. (24) to (35)) in Eq. where 1 , 2 and 3 are the ordinates from the parametric space in which the radial and tangential concentrated forces and concentrated moment are applied, respectively.
Creating the vectors: we find an analogous form to the traditional finite element method, since corresponds to the stiffness matrix, is the vector of the displacements and is the vector of external forces.