ISOGEOMETRIC COLLOCATION MIXED METHODS FOR RODS

Isogeometric collocation mixed methods for spatial rods are presented and studied. A theoretical analysis of stability and convergence is available. The proposed schemes are locking-free, irrespective of the selected approximation spaces.

1. Introduction. Isogeometric Analysis (IGA), introduced by Hughes et al. [8,12], is a powerful analysis tool aiming at bridging the gap between Computational Mechanics and Computer Aided Design (CAD). In its original form IGA has been proposed as a Galerkin method where the geometry is represented by the spline functions typically used by CAD systems and, invoking the isoparametric concept, field variables are defined in terms of the same basis functions used for the geometrical description.
The high regularity properties of the employed functions not only lead in many cases to superior accuracy per degree of freedom with respect to standard FEM (cf., e.g., [5,9,13,15,16]), but opens the possibility to employ the IGA concept within the collocation scheme framework (see [2,3]), Within this context, a comprehensive study on the advantages of isogeometric collocation over Galerkin approaches is reported in [16].
Aim of this paper is to report our results regarding the IGA collocation scheme when applied to elastic rod problems, as detailed in [4].
The outline of the paper is as follows. In Section 2, we briefly review the basic idea of IGA, after having introduced the Non-Uniform Rational B-Splines (NURBS) description of a spatial curve. Section 3 presents the rod model, while in Section 4 the discrete scheme is introduced together with our error estimates. Section 5 shows some numerical tests supporting the theoretical convergence results. Finally, a few conclusions are provided in Section 6.
2. NURBS-based isogeometric analysis. We begin by recalling that B-Splines are smooth functions constructed by piecewise polynomials. A B-spline curve in R d is composed of linear combinations of B-spline basis functions and coefficients (B i ). These coefficients are points in R d , referred to as control points.
To define such functions we introduce a knot vector as a set of non-decreasing real numbers representing coordinates in the parametric space of the curve where p is the order of the B-spline and n is the number of basis functions (and control points) necessary to describe it. For p = 0 (piecewise constants): for p ≥ 1 : where, in (3), we adopt the convention 0/0 = 0.
If internal knots are not repeated, B-spline basis functions are C p−1 -continuous. If a knot has multiplicity m the basis is C k -continuous at that knot, where k = p−m. In particular, when a knot has multiplicity p, the basis is C 0 and interpolates the control point at that location. We define To obtain a NURBS curve in R 3 , we start from a set B ω i ∈ R 4 (i = 1, ..., n) of control points ("projective points") for a B-spline curve in R 4 with knot vector Ξ. Then the control points for the NURBS curve are where [B i ] k is the k th component of the vector B i and ω i = [B ω i ] 4 is referred to as the i th weight. The NURBS basis functions of order p are then defined as The NURBS curve α α α is defined by As usual, we denote the support of the curve α α α by Γ(α α α), hence Γ(α α α) ⊂ R 3 . In addition, we suppose that the map α α α : [0, 1] → Γ(α α α) is smooth and invertible, with smooth inverse denoted by α α α −1 : Γ(α α α) → [0, 1]. Following the isoparametric approach, the space of NURBS vector fields on Γ(α α α) is defined, component by component, as the span of the push-forward of the basis functions (6): We finally note that the images of the knots through the function α α α naturally define a partition of the curve support Γ(α α α) ⊂ R 3 , called the associated mesh M h , h being the mesh-size, i.e., the largest size of the elements in the mesh. 3. The rod mixed model. We here briefly present a variant of the (scaled) model of [1], where, in order to include the NURBS parameterization directly into the analysis, the parameterization used to describe the rod is not assumed to be the curvilinear abscissae. Therefore, we assume that the rod axis is defined by a NURBS curve α α α(ξ), with ξ ∈ [0, 1], see (7). Accordingly, the tangent vector is α α α (ξ), and we notice that α α α (ξ) is not necessarily of unit length. In order to develop the analysis, we require the minimal regularity assumption α α α ∈ C 1 [0, 1]. In the considered mixed formulation, the unknowns are the displacements v(ξ), the rotations ϕ ϕ ϕ(ξ), and a variable, τ τ τ (ξ), associated to internal forces. Assuming for simplicity, but without loss of generality, clamped boundary conditions, the problem to be solved is the following: In (9), the vector field q ∈ C 0 [0, 1] represents the load acting on the rod, the tensor fields E ∈ C 1 [0, 1] and A ∈ C 0 [0, 1] are uniformly positive definite and symmetric, and they are associated to the given material law and section geometry, while d is a slenderness parameter (cf. [1,7]). Moreover, we notice that data q, E, A are suitably scaled with respect to the local length factor of the rod parametrization. For the collocation method purposes, we rewrite system (9) as: Using the variational approach of [7] and standard regularity results, one gets the following proposition. (9) (and, therefore, also to Problem (10)). Moreover, it holds: 4. The collocation scheme. We now present the collocation method for the rod model introduced previously.
, used for the rotation approximation, and associated with the knot vector The knot vector (12) will be used for each of the three components of the approximated rotation field. Accordingly, we set (cf. (4) and (8)): Analogously, we define the NURBS space for the displacement approximation, and associated with the knot vector (used component-wise) Finally, we define the NURBS space for the internal force approximation, and associated with the knot vector (used component-wise) We notice that it holds Remark 1. We remark that the three knot vectors above induce, in principle, three different meshes: with corresponding mesh-sizes h ϕ , h v , and h τ . In the sequel, we will set h = max {h ϕ , h v , h τ }. However, we notice that, in practical applications, the three meshes most often coincide.
Remark 2. We remark that, in principle, one might also think of using different knot vectors for the different components of the approximated fields. However, this latter choice does not seem to be of practical interest.
In the sequel, we will also use the spaces of first and second derivatives: whose dimensions are given by dim (18). Furthermore, we introduce suitable sets of collocation points in [0, 1]: We notice that it holds 3( We are now able to present the proposed scheme. Given the finite dimensional spaces defined in (13), (14), and (16), together with the collocation points introduced in (20), the discretized problem reads as follows.
We finally present the following fundamental assumption on the collocation points.
A comment on possible practical collocation points satisfying Assumption 4.1 can be found in Section 4.1.

Error estimates.
We now state the main theoretical result proved in [4], concerning the convergence of the scheme (22).
and where the constant C is independent of the knot vectors and the thickness parameter d.
We remark that the theoretical results establish error estimates in the W 2,∞norm while the convergence plots in Section 5 are reported in terms of L 2 -norm errors, which are more relevant in engineering applications. However, we point out that, since the L 2 -norm is bounded from above by the W 2,∞ -norm, our theoretical error estimates hold for the L 2 -norm as well.
One could extend the above results to the case of less regular loads q, obtaining a lower convergence rate β. Moreover, approximation results in higher order norms can also be derived by using inverse estimates. We also remark that equation (24) should not be intended as a recipe to find the optimal balancing among p ϕ , p v and p τ since the provided estimates are not sharp.
The optimal selection of points for interpolation of one-dimensional splines is addressed in various papers. The only choice proven to be stable (i.e., satisfying Assumption 4.1) for any mesh and degree are the so-called Demko abscissae, see for instance [10,11]. A different approach proposed in the engineering literature [14] is to collocate at the Greville abscissae. We refer to [2] for a deeper investigation and comparison between the Demko and Greville choices.
Remark 3. Theorem 4.1 yields a converge estimate, uniform in the thickness parameter, without requiring any particular compatibility condition among the three discrete spaces Φ Φ Φ h , V h , and Γ Γ Γ h . Therefore, the proposed method is locking-free regardless of the chosen polynomial degrees and space regularities. Even different meshes can be adopted among the three spaces. Such result is surprising, at least in comparison with Galerkin schemes. Indeed, in typical Galerkin approaches the discrete spaces Φ Φ Φ h , V h , Γ Γ Γ h must be carefully chosen, in order to avoid the locking phenomenon and the occurrence of spurious modes.
5. Numerical tests. In this section, some numerical experiments are shown to test the accuracy of the method.
We alway collocate the relevant equations at the physical images of the Greville abscissae of the unknown field knot vectors as described below (cf. [6]).
The Greville abscissae related to a spline space of degree p and knot vector {ξ 1 , . . . , ξ n+p+1 } are points of the parametric space defined by: In addition, we also need the Greville abscissae related to the k-th derivative space, which are defined as: We consider a 90 • circular arch, of radius r and thickness t, clamped on one end and subjected to an out-of-plane concentrated load on the other end, as sketched in lation. Since the system is statically determinate, stress resultants can be obtained directly from equilibrium considerations. Figure 2 shows the stress resultants as functions of the angle θ, expressed in terms of an intrinsic basis [ζ 1 , ζ 2 , ζ 3 ]. The analytical solution for the displacements in z-direction is therefore: This reference solution is used to compute the L 2 -norm of the error of displacements. The convergence studies are performed for two different slenderness ratios, namely t/r = 10 −1 and t/r = 10 −4 , and for polynomial degrees ranging from 3 to 8. The L 2 -norm of the error of displacements is plotted versus the total number of collocation points in Figures 3-4. Figure 3 shows the convergence plots using equal orders for displacements, rotations, and internal forces. A very good convergence behavior is observed independently of the slenderness. These results confirm the locking-free properties of mixed collocation methods. This characteristic is independent of the choice of spaces for the three fields, i.e., the spaces for displacements, rotations, and internal forces can be chosen freely without any inf-sup-like condition to be fulfilled by the spaces. To further test this, the example is repeated with the following "exotic" choice of polynomial degrees: p v = p ϕ − 1 = p n − 2. The results are plotted in Figure 4 and confirm the behavior described above.
It should also be noted that in Figures 3 and 4, zig-zagging of the curves due to converged results is really at machine precision, for both the thick and the thin cases, which proves that the mixed formulation does not cause conditioning problems.  Figure 3. Circular arch with out-of-plane load, solved by a mixed collocation method with p v = p ϕ = p n = p. L 2 -norm of displacements error versus number of total collocation points, for t/r = 10 −1 (left) and t/r = 10 −4 (right). Dashed lines indicate reference orders of convergence. 6. Conclusions. In this work we have presented the application of isogeometric collocation techniques to the solution of spatial rods, as detailed in [4]. The obtained collocation scheme has been numerically tested in order to assess their accuracy and efficiency. In particular, it is interesting to highlight that the considered mixed formulations appeared to be locking-free for any choice of the discrete spaces for  displacements, rotations, and internal forces; such a remarkable behavior has also been analytically proven in the second part of the paper. These results propose isogeometric collocation methods as a viable and efficient alternative to standard approximation methods for curved beams.