Weak coupling of nonlinear isogeometric spatial Bernoulli beams

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. xxx (xxxx) xxx www.elsevier.com/locate/cma

Weak coupling of nonlinear isogeometric spatial Bernoulli beams A.M. Bauer ∗, R. Wüchner, K.-U. Bletzinger Lehrstuhl für Statik, Technische Universität München, Arcisstr. 21, 80333 München, Germany Received 29 July 2019; received in revised form 5 November 2019; accepted 6 November 2019 Available online xxxx

Abstract This paper proposes a coupling scheme for isogeometric elements, which is valid for large rotations and displacements. Two elements can be coupled not only at interpolated control points but also inside the parametric domain of a NURBS patch. Furthermore, each cross section may be arbitrarily oriented in space. The formulation is applicable for all structural elements, where an orthogonal system of base vectors can be derived. An Euler–Bernoulli beam is used as an example in order to illustrate the demands on the coupling conditions in the proposed approach. Moreover, the conditions are chosen such that different types of joints, e.g. scissor joints, can be modeled and such that they are easy and fast to implement. The coupling is incorporated by an additional term in the Principle of Virtual Work which is here computed using a penalty approach. The Euler–Bernoulli beam is summarized in order to introduce the notations in this contribution. This is followed by the proposed coupling conditions in the weak form. Subsequently, alternatives for the coupling conditions are briefly introduced. Eventually, benchmark and demonstrator examples validate and illustrate the proposed coupling methodology. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Isogeometric B-Rep Analysis (IBRA); Isogeometric analysis; Multipatch NURBS; Weak coupling; Nonlinear beam formulation; Large rotations

1. Introduction Hughes et al. [1] introduced Isogeometric Analysis (IGA) in 2005. The properties of Non-Uniform Rational B-Splines (NURBS) as basis functions for the approximation of the solution within Finite-Element analysis (FEA) are beneficial in many cases. However, certain properties lead to challenges, for example the non-interpolating property of the control points with the degrees of freedom (dofs). If the boundaries of the patches match in their discretization, kinematic constraints can be applied to the structural system in order to satisfy the coupling conditions [2–4]. Additionally, some of the well-established isogeometric, dimensionally-reduced elements, such as the Kirchhoff–Love shell by Kiendl et al. [2], do not have rotational dofs to directly couple the rotation. To overcome the definition of complex kinematic constraints, Kiendl et al. proposed the bending strip method in [5], so the beginning of the analysis of isogeometric multipatches was restricted to matching parametrization at a C0 -continuity, which is often at the boundary. However, the matching discretization cannot be guaranteed for complex geometries. A weak imposition independent of the discretization thus becomes necessary. Much research ∗ Corresponding author.

E-mail address: [email protected] (A.M. Bauer). URL: http://www.st.bgu.tum.de (K.-U. Bletzinger). https://doi.org/10.1016/j.cma.2019.112747 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝

Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

2

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 1. Coupling scenarios in 2d space with increasing complexity: (a) matching tangents/cross sectional directors and interpolated control points at patch boundaries, (b) non-matching tangents/cross sectional directors and interpolated control points and (c) non-matching tangents/cross sectional directors and non-interpolated control points.

has been done on coupling of isogeometric surface elements. The weak imposition of boundary conditions was examined for the penalty approach [6,7], Lagrange methods [8–10], Nitsche method [11–17] and Mortar method [9,18,19], among others. Comparisons of several methods can be found in [8]. Many investigations have also been conducted on solids, e.g. in [11,12,20,21], and solid-shell coupling, e.g. in [14]. Nevertheless, not many publications focus on 1D-elements even though many contributions have been recently published on spatial isogeometric beam elements for linear [22–24] and nonlinear [25–29] applications, among others. In general, there are different coupling scenarios in 2d and 3d space (see Figs. 1 and 2). The simplest case is a coupling at interpolated control points with matching tangents and cross section. In 3d space, the cross sections can also be non-matching for parallel tangents (cf. Fig. 2(b)). If the tangents remain parallel during the deformation, this case can be treated like the matching cross sections if the gap can be described by a local twisting angle around the tangent. The next step in complexity is integrating non-matching tangents. The bending strip method [5] can be adopted in 2d space in order to enforce the kink during the deformation. However, the bending strip is not directly applicable in 3d space as it cannot enforce the torsional coupling. The most general case is a coupling inside the parametric domain, where no control points interpolate the curve with non-matching tangents and cross sections. Note that the 2d case is much simpler since the rotational axis is always defined by the normal of the plane. Hence, the rotation can be described by one angle around a fixed axis, whereas 3d space requires a more sophisticated method to uniquely describe the rotation of both sides, which is especially the case for a geometrically nonlinear application. Greco and Cuomo proposed a G 1 -coupling for Euler–Bernoulli kinematics in [30] at interpolated control points. Coupling of Timoshenko beams in the strong form may inherently be given if the degrees of freedom for the cross section are defined globally, e.g. by quaternions in Weeger et al. [29]. In the present paper, we propose an approach to directly deal with the last and the most general scenario for structural elements without these respective dofs. Furthermore, the formulation will be valid for large displacements and rotations in 3d space in the context of a geometrically nonlinear analysis. The present contribution is structured as follows. A brief overview of the specific notation and the structural element formulation of the nonlinear isogeometric spatial Bernoulli beam [26], which will also be used in the numerical examples, is provided in Section 2. This is followed by the newly proposed coupling procedure in 3d space for arbitrary configurations (Section 3). Additionally, some alternatives/variations of the approach are briefly presented. The first part of Section 4 evaluates the quality of the results in comparison with analytical and numerical benchmark solutions. The second part illustrates the possibilities of the coupling formulation for large displacements, scissor joints and complex models. The final section presents a conclusion and an outlook on future research. 2. Structural mechanics The proposed coupling formulation and the beam formulation used are derived from the Principle of Virtual Work. A virtual displacement δu is applied to the system. The performed work of internal and external forces is zero, if the system is in equilibrium [31]: (∫ ) ∫ ∫ δW = δWint + δWext = − S : δE d x + t · δu d x + ρ0 B · δu d x = 0 (1) Ω

ΓN

Ω

Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

3

Fig. 2. Coupling scenarios in 3d space with increasing complexity: (a) matching cross section and interpolated control points at patch boundaries, (b) non-matching cross section, parallel tangents and interpolated control points, (c) non-matching cross section, non-parallel tangents and interpolated control points and (d) non-matching cross section and tangents and non-interpolated control points.

2.1. Beam formulation In this section, the element formulation of Bauer et al. [26] is briefly summarized in order to clarify the requirements to the later proposed coupling formulation for 1d-elements. The applied isogeometric geometrically nonlinear beam accounts for Bernoulli kinematics. The cross sections remain orthogonal to the center line after deformation and there are no changes of the cross sectional dimensions. Torsion according to St. Venant is considered whereas warping effects are neglected. Three translational degrees of freedom and one for a relative rotation around the beam axis are used for the spatial beam. However, the presented coupling formulation is not restricted to this beam formulation but can be used for every element formulation based on a local, orthogonal coordinate system. 2.1.1. Geometric description The continuum of the undeformed and deformed beam, X and x respectively, can be defined by the longitudinal beam axis, represented by a NURBS curve (Xc resp. xc ) and two normalized base vectors Aα , α ∈ {2, 3}, which are perpendicular to the each other and the beam axis A1 (see Fig. 3). Note that upper-case letters will refer to the undeformed and lower-case letters to the deformed configuration. As continuum mechanics is applied, there are three convective contravariant coordinates denoted as θ i with i ∈ {1, 2, 3}. Derivatives with respect to one convective, contravariant coordinate ∂(·) will be written as (·),i . Under the Bernoulli assumption, every kinematic ∂θ i can be described by: ( ) ( ) ( ) ( ) X θ 1 , θ 2 , θ 3 = Xc θ 1 + θ 2 A2 θ 1 + θ 3 A3 θ 1 (2a) ( 1 2 3) ( 1) ( 1) ( 1) 2 3 x θ , θ , θ = xc θ + θ a2 θ + θ a3 θ (2b) The peculiarity of this beam formulation is in the determination of the base vectors Aα . Two mapping matrices based on the Euler–Rodrigues formula are applied to define the undeformed configuration of the beam. Two initial Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

4

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 3. Geometry description of the Euler–Bernoulli beam with position and base vector in undeformed and deformed configuration.

base vectors at the beginning of the beam A0α are mapped to every position of the center line by the mapping matrix Λ(T0 , T) in order to create a unique orientation. Subsequently, they are aligned to the real beam orientation by a rotation around the beam axis by the rotation matrix RT (Ψ ). This results in ( ( )) ( ( )) ( ) (3) Aα θ 1 = RT(θ 1 ) Ψ θ 1 Λ T0 , T θ 1 A0α , where T is the normalized tangent of the beam axis and T0 the tangent at the beginning of the beam. The vectors Ti are the normalized Ai vectors, i.e. represent the orthonormal local coordinate system. We use the same mapping for the base vectors aα of the deformed beam. In contrast to the undeformed 0 configuration, ( )not the initial base vectors at the beginning of the beam Aα are mapped by Λ (T, t) but the base vectors Aα θ 1 of the undeformed configuration. The correction angle ψ is applied by the rotation matrix Rt (ψ) around the axis of the deformed beam t, which is also the fourth degree of freedom. ( ( )) ( ( ) ( )) ( ) ( ) (4) aα θ 1 = Rt(θ 1 ) ψ θ 1 Λ T θ 1 , t θ 1 Aα θ 1 2.1.2. Nonlinear kinematics The Green–Lagrange (GL) strain tensor and the energetically conjugated second Piola–Kirchhoff (PK2) stress tensor are used for the Principle of Virtual Work. Consequently, the covariant base vectors Gi of the continuum are required. They are computed by the differentiation of the position vector X by the respective contravariant coordinate θ i . The same applies for the deformed configuration. G1 = A1 + θ 2 A2,1 + θ 3 A3,1 ,

G 2 = A2 ,

G 3 = A3

(5)

The GL strain tensor is defined as ) ) 1( 1( gi j − G i j Gi ⊗ G j , E i j = gi j − G i j (6) E= 2 2 The curvilinear coordinate system is used for the calculation of the GL strain tensor. However, the strains need to be transferred to the Cartesian coordinate system due to the constitutive law. Entities referring to a Cartesian coordinate system will be denoted with ˜ •. Ei j ˜i j = E (7) ∥Gi ∥2 ∥G j ∥2 In order to simplify the derivation(of the)element for slender beams with h, w ≪ L, the following (( formulation ( 2 )2 ) ) ( 2 3) 3 2 is assumed: quadratic order terms O θ ,O θ , O θ · θ are neglected and ∥G1 ∥2 is approximated by ∥G1 ∥2 ≈ ∥A1 ∥2 . These simplifications imply as well that warping effects vanish. Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

5

Consequently, the strain tensor is left with the following non-zero components: B3

B2

E 11

b3 b2             1 1 2 3 = (g11 − G 11 ) = (a11 − A11 ) + θ (a2,1 · a1 − A2,1 · A1 ) + θ (a3,1 · a1 − A3,1 · A1 ) 2 2 ϵ    κ31 κ21       1 2 3 = (a11 − A11 ) +θ (b2 − B2 ) +θ (b3 − B3 ) 2 Tα

(8)

κβα

tα 1 1       1    E 1α = (g1α − G 1α ) = θ β (aβ,1 · aα − Aβ,1 · Aα ) = θ β (tα − Tα ) 2 2 2 The respective strain in the Cartesian coordinate system is 2 3 ˜11 = E 11 = ϵ + θ κ21 + θ κ31 E 2 2 ∥G1 ∥2 ∥A1 ∥2 E 1α θ β κβα ˜1α = E = , ∥G1 ∥2 ∥Gα ∥2 2∥A1 ∥2

(9)

(10) (11)

where ϵ is the axial strain and καi is the change in curvature in direction of the base vectors aα . καi consists either of bα and Bα denoting the curvature of the respective configuration in the particular direction or tα and Tα denoting the change of the torsion in the structural element formulation of the Euler–Bernoulli beam. The second Piola–Kirchhoff stresses S are then computed with the elasticity tensor C: S=C:E

(12)

2.1.3. Internal work The following term for the internal virtual work results from inserting the non-zero strain components in the equation of the weak form of the equilibrium as defined in Eq. (1), ∫ ∫ ˜ δWint = − S : δE dΩ = − S 11 δ E˜ 11 + ˜ S 12 δ E˜ 12 + ˜ S 13 δ E˜ 13 dΩ (13) Ω⊂R3

Ω⊂R3

This term can be pre-integrated over the cross section assuming that the base vectors Aα are the principal axes of the cross section: ∫ δWint = − E E˜ 11 δ E˜ 11 + G E˜ 12 δ E˜ 12 + G E˜ 13 δ E˜ 13 dΩ Ω⊂R3 ( ) ∫ ( ) E GI 1 1 · Aϵ δϵ + IA3 κ21 δκ21 + IA2 κ31 δκ31 + · − κ32 δκ32 + κ23 δκ23 d L =− 4 2 2 ∥A1 ∥2 2 L ∥A1 ∥2 (14) Young’s modulus E and shear modulus G originate from the constitutive law of a linear elastic material. 3. Coupling formulation A special coupling condition is needed in the context of isogeometric analysis since the isogeometric FE-nodes with the degrees of freedom are in general not interpolating and since there are structural elements without global degrees of freedom for rotations. Consequently, the strong form of the coupling condition is not possible, even for displacements, in the general case if the geometry should remain unchanged. By allowing a change of geometry, it is legitimate to insert as many knots as necessary for a C 0 -continuity in order to generate an interpolating control point at the coupling location (see Fig. 4), which is e.g. used by Greco et al. [30]. The continuity is then maintained by deriving respective kinematic constraints for the dofs on the nodes. This procedure can be bypassed by a weak formulation of the coupling problem as we will see in the following. Generally, there are several methods to apply boundary conditions in a weak sense to a structural system, e.g. penalty approach, (augmented) Lagrange method or Nitsche method. We chose the penalty approach in this contribution as it provides sufficiently good results [32,33] for isogeometric shell models and is a very flexible and straight-forward to implement approach. Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

6

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 4. Two coupled curves with their control points: (a) no interpolation of the curve through a control point at the coupling point and (b) with respective C 0 -continuity at the coupling point for interpolated control points.

Fig. 5. Coupling conditions: (a) same displacement and rotation and (b) same relation between position and base vectors in initial and actual configuration.

3.1. General weak form of the coupling formulation There are basically two conditions which can be used for enforcing coupling. It can be assured that either both sides of the coupling condition follow the same translation and rotation (see Eq. (15)) or that the relation between the both sides (see Eq. (16)) is constant (see also Fig. 5). Both approaches result in an additional term in the virtual work: ( ) ( ) ( ) ( ) δW coup = α disp · um − us · δum − δus + α rot · wm − ws · δwm − δws (15) (( m ) ( m )) (( m ) ( m )) ( ( m s) ( m s )) coup disp s s s s rot δW = α · x − x − X − X · δ x − x − X − X + α · r t , t − δr T , T (16) with α disp being the penalty factor for displacements and α rot for rotations. u and w represent translation and rotation whereas r (Tm , Ts ) quantifies the relation between two vectors on the master and slave sides. A preliminary study on the choice of penalty factor for the proposed coupling conditions was conducted and will be presented in Section 4.1. A more extensive study will be part of future research. coup coup The respective contributions to the right-hand-side Rr and the stiffness matrix K r s are defined similar to the structural element formulations by deriving the additional term with respect to the degrees of freedom uˆ r /uˆ s . This is shown for Eq. (15) since this will eventually be chosen. ( ( ) ∂ (wm − ws ) ) ∂ (um − us ) + α rot · wm − ws · (17) Rrcoup =α disp · um − us · ∂ uˆ r ∂ uˆ r m s m s m s m s ( m ) ∂ 2 (wm − ws ) disp ∂ (u − u ) ∂ (u − u ) rot (w − w ) ∂ (w − w ) rot s K rcoup =α · · + α · · + α · w − w · s ∂ uˆ s ∂ uˆ r ∂ uˆ s ∂ uˆ r ∂ uˆ r ∂ uˆ s (18) The derivation of a coupling condition for large displacements and rotations is shown at the example of a beam since it is one of the most obvious application example. However, it can be transferred to every other formulation, Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

7

Fig. 6. Description of the rotation of base vectors by their tip displacement: (a) one displacement value and respective valid configuration and (b) two displacement values.

such as the Kirchhoff–Love shell, with respective properties. The coupling of different structural elements, e.g. shell and beam, is also possible. 3.2. Coupling conditions A general, nonlinear coupling condition is required to be unique without singularity for every configuration of connection and displacement. As already mentioned, there are in principle two ways of describing a coupling condition. Either the relation between the two triads of the beams is preserved or it has to be ensured that the points including the respective triads on each side deform in the same way. We have chosen a coupling condition which enforces the same displacement and rotation on both sides since the formulation should be general and applicable to rotational supports. By choosing the relative condition, we would need to define a reference triad at a support, which is less intuitive. However, the two conditions are convertible into each other. The respective changes will briefly be shown after the proposed, implemented and tested formulation. Describing large displacements and rotations in a unique way in 3d space can be demanding for arbitrary initial configurations. Therefore, the approach of Lumpe and Gensichen [34] for a nonlinear beam formulation was used and modified. There, the displacements of the tips of the base vectors were used to describe the rotation of the cross section of a nonlinear beam element from the reference configuration to the actual one. These displacements are computed by the difference of the reference and the actual base vector. The displacements of the three base vectors are then split into parts which are oriented along the reference base vectors. Consequently, nine projected displacements evolve for an arbitrary rotation in space. In contrast to [34], the amount of considered projected tip displacements was not reduced to three displacements describing the rotation in space for this coupling condition. We based this decision on the fact that the orthonormal base vectors do not have to be kept orthonormal as they are defined by the structural element and not the rotation components. Furthermore, a singularity can be avoided for very large rotations by taking every tip displacement component into account. Fig. 6 shows a simple rotation around T3 . With the reduction, only the displacement of the tip of T2 in T1 -direction would be taken into account and consequently there are two valid configurations tiα (see Fig. 6(a)) with the same displacement u 2 . If the second displacement component of the tip is also included in the formulation the rotation is uniquely defined since the component of the tip displacement in the direction of the base vector is not equal for the two possible configuration of Fig. 6(a) (see Fig. 6(b)). Note that in our case of 3d rotations of an orthonormal triad, the tip displacements of only two of three base vectors would be sufficient to uniquely describe the rotation. However, all of them are considered in order to have a general implementation. Additionally, the rotations are expressed by a projection into the local coordinate system of the initial beam of the master side in this contribution, which is indicated by w•/Tim . This allows the coupling of selected directions. Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

8

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 7. Rotation of the cross sectional triad and respective tip displacements of the base vectors separated into the components of the initial base vectors.

As the same displacement and rotation should be enforced on both sides, the following coupling conditions can be used, where u is the displacement of the spatial point of the center line and w is the difference of the base vectors in the actual and reference configuration: s um x = ux

wtmm /Tm 1 1

=

w ss,mod m t1 /T1

wtmm /Tm = w ss,mod 2

1

t2

wtmm /Tm = w ss,mod 3

1

u my = u sy

t3

/Tm 1 /Tm 1

wtmm /Tm 1 2

=

wss,mod m t1 /T2

wtmm /Tm = w ss,mod 2

2

t2

wtmm /Tm = w ss,mod 3

2

s um z = uz

t3

/Tm 2 /Tm 2

wtmm /Tm 1 3

=

wtmm /Tm = wss,mod 2

3

t2

wtmm /Tm = wss,mod 3

3

(19)

wss,mod m t1 /T3

t3

/Tm 3 /Tm 3

(20) (21) (22)

with ( ) wti /Tmj = ti − Tim · Tmj

(23)

All considered tip displacements wti /T j are shown in Fig. 7. Since the triads of the initial configuration are not aligned in the general case and the rotation should be described by the tip displacements of the base vectors, one side has to be mapped in a unique way to the other in order to describe the displacement with respect to one local initial coordinate system. This can be done by using the projection rule of the initial coordinate systems also for the actual coordinate systems. A mapping of a base vector tis from the slave side to the master side tis,mod can be computed by the following ( ) ( ) ( ) tis,mod = Tim · Ts1 ts1 + Tim · Ts2 ts2 + Tim · Ts3 ts3 (24) This mapping is unique for arbitrary configurations and preserves the initial relation, which is also illustrated in Fig. 8. Fig. 9 shows an overview of all steps in the coupling formulation, which have to be conducted. The relation between the master and slave point in the reference configuration is used to transform the actual configuration of the slave point. These modified base vectors tis,mod and the actual base vectors of the master point tim are then used for the difference vector with the reference master base vectors Tim . The difference vectors are projected into the reference master coordinate system. These projections wtm/s,mod /Tm are used to evaluate the difference in rotation of i j both sides by comparison. Moreover, coupling conditions can be applied for only one base vector. This corresponds to a scissor joint which can also change its orientation in space during the deformation process (see Fig. 10). The respective conditions for Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

9

Fig. 8. Mapping operation from slave to master side: (a) beam configuration in undeformed and deformed configuration and (b) splitting a base vector from the master side into the components of the slave triad.

Fig. 9. All mapping steps of the proposed coupling condition.

a pin in tm 2 direction would e.g. be: s um x = ux

wtmm /Tm 2 1

=

wss,mod m t2 /T1

u my = u sy wtmm /Tm 2 2

=

w ss,mod m t2 /T2

s um z = uz

wtmm /Tm 2 3

=

wss,mod m t2 /T3

(25) (26)

Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

10

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 10. Scissor joints: (a) Free torsion with wt1 /Tim , (b) free bending moment Mv with wt2 /Tim and (c) free bending moment Mn with wt3 /Tim .

Note that these displacements could also be transformed into angles which would give a direct analogy to a rotational spring. However, this would result in a more complex computation of the derivatives and variations in the implementation. 3.3. Alternative coupling conditions This section briefly shows how the components of the proposed coupling condition can be used in alternative approaches. 3.3.1. Displacements w.r.t. to the local coordinate system The displacements can also be expressed in the local coordinate system of the master. This is simple for the reference local system, which can be used for e.g. supports. m m s m s um Tm = u · Ti = u · Ti = u Tm i

i

(27)

The same holds for the local coordinate system in the actual configuration. This is for example useful for the modeling of shearing force hinges. m m s m s um tm = u · ti = u · ti = u tm i

i

(28)

In contrast to the local reference frame, this is not directly compatible with the implementation of the originally proposed method as additional variations arise with the projection to tim since tim is dependent on the dofs. 3.3.2. Rotations w.r.t. to the global coordinate system The rotation of tim and tis,mod could also be measured in the global coordinate system which results in a small change in the coupling conditions. (m ) m wm (29) ti = ti − Ti ( ) ws,mod = tis,mod − Tis,mod (30) ti s,mod wm ti = wti

(31)

This modification is especially useful for boundary conditions, which are aligned to the global system, if the local coordinate system is not applicable. It implies only a minimal change in the implementation and both options can easily be provided. 3.3.3. Relative condition The relative coupling condition of the displacements is directly convertible to the other one. (( ) ( )) (( ) ( )) δW coup = α disp · xm − xs − Xm − Xs · δ xm − xs − Xm − Xs (( ) ( )) (( ) ( )) = α disp · xm − Xm − xs − Xs · δ xm − Xm − xs − Xs ( ) ( ) = α disp · um − us · δ um − us

(32)

Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

11

Fig. 11. Set-up of a cantilever with tip load consisting of two beams. Source: Adapted from [26].

Fig. 12. Relative error of the u y -displacement of a cantilever under a tip load modeled by two coupled beams with different amount of elements: (a) 16 elements and (b) 64 elements.

The proposed approach can also provide a relative condition for the rotations. Comparing the deformed base vectors tim and tis,mod would be for example a relative condition since the modified slave base vectors are computed such that they reproduce the initial relation between the two coupling partners. tim = tis,mod

(33)

The relative angles between all base vectors could also explicitly be measured and be enforced to be constant. 4. Numerical examples Several benchmarks will be performed in the following in order to prove the validity and generality of the proposed coupling formulation for an Euler–Bernoulli beam. This is followed by some application examples in order to show its applicability for more complex scenarios. 4.1. Cantilever A straight cantilever modeled by two beams is shown as first benchmark (see Fig. 11). The coupling point is in the middle of the total length of the cantilever. It is loaded by a tip load. A linear analysis is performed in order to have analytic results to compare with. Fig. 12 shows the relative error of the displacement for a range of penalty factors for rotation coupling. The penalty factor of displacements is kept constant to α disp = 1010 since it was observed that the result is only affected F by the error which is related to the modeling u disp = αdisp = 100.110 = 10−11 . Furthermore, we varied the number of elements per beam and the polynomial degree. The relative error of the tip displacement is computed with the analytic solution as reference as follows: εu y =

u y − u ref y u ref y

with u ref y =

F L2 3E Iz

(34)

Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

12

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 13. Set-up of a cantilever with tip load consisting of two beams with circular cross section where the second beam is rotated around its axis.

It can be observed that the quality of the results improves in the beginning as expected for an increasing penalty factor. After a certain level, here α rot = 109 , the quality of the result is corrupted by the conditioning of the equation system solved in the analysis. We use a circular cross section as a second benchmark example in order to be able to check different relative orientations of the cross sections against the single beam solution (see Fig. 13). Furthermore, a geometrically nonlinear analysis is conducted in order to check large displacements. Therefore, the cantilever is bent by a single force at the tip. The first beam is kept constant whereas the second is rotated around its axis. Each coupled beam is modeled by a B-Spline with polynomial degree p = 5 and 16 elements, which corresponds to 21 control points. The single beam is modeled by 32 elements with polynomial degree p = 5, which results in 37 control points. The penalty factors are α disp = 1012 and α rot = 109 . We observe that the evolving gap between the two coupled ends correlates to the penalty factor which was chosen to be α disp = 1012 (see Fig. 14(a)). The penalty factor acts like a spring with stiffness α disp between the two beams. The gap opens only in the direction of the load and the length of the gap corresponds to elongation of a spring under the respective load which is transferred in this joint. F 25 = 12 = 2.5 · 10−11 α disp 10 The relative angle between the two beams remains almost constant as can be seen in Fig. 14(b). ∆u = ∥um − us ∥2 = ∆u z =

(35)

εφi = arccos(tim · tis ) − arccos(Tim · Tis )

(36)

The relative error of the absolute displacement at the coupling point umid and at the tip utip is computed using the single patch solution as reference. The results are depicted in Figs. 14(c) and 14(d). εu mid = i

mid mid u i,coupled − u i,single mid u i,single

tip

εu tip = i

tip

u i,coupled − u i,single tip

(37)

u i,single

The coupled beams show very good accordance with the single beam solution, also for varying correlation between the base vectors of the cross section. Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

13

Fig. 14. Results of the cantilever beam with rotating cross section: (a) Evolving gap between the two beams, (b) relative error of the angle between the cross sectional base vector, (c) relative error of the displacements at the coupling point and (d) relative error of the displacements at the tip.

4.2. Circular arch cantilever The next example will test the most complex loading on the coupling, i.e. bending torsion interaction. Therefore, a circular arch cantilever with a point load on the free end is chosen (see Fig. 15). The beams are modeled with the same parametric discretization but the polynomial degree and the number of elements are varied. The linear analytic results of the tip displacement and support forces are computed as follows: ) ∫ π( 2 1 1 tip uz = (−F R · cos α)(−R · cos α) + (F R · (1 − sin α))(R · (1 − sin α)) R dα E I3 GI 0 ( ) F R3π F R 3 3π = + − 2 = 0.138285322802514 m , (38) 4E I3 GI 4 ) ∫ π( 2 1 1 φ tip = (−F R · cos α) · cos α + (F · (1 − sin α)) · sin α R dα E I3 GI 0 F R2π F R2 ( π) =− + 1− = 0.0757632997455722 rad (39) 4E I3 GI 4 Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

14

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 15. Set-up of a quarter circle cantilever. Source: Adapted from [26].

M˜ n = − Fz R = −1.0 kNm M˜ T = Fz R = 1.0 kNm

(40) (41)

The first evaluation of this example compares the relative error for different penalty factors α and polynomial degrees p. The plots in Fig. 16 show the relative error w.r.t. the single-beam-solution in order to separate discretization errors from the investigated coupling errors. The coupled model has 4 or 16 elements per beam, which corresponds to 8 or 32 elements in the single-patch model. Fig. 16 shows a different behavior for polynomial degree p = 2, 3. This convergence to a different result is linked to the still existing differences in the discretization. Fig. 17 shows the control points and weights of two 45◦ -segments, constructed by two independent segments with CAD (see Fig. 17(a)) and with knot insertion in the middle of a 90◦ -segment (see Fig. 17(b)). One can observe in Fig. 17(c) that the weights and the positioning of the Gauss points differ for these parametrizations, where each curve was refined to 4 nonzero knots spans, which in return influences the outcome of the result if the curves are not refined enough. This is confirmed in Fig. 18 where the discretization of Fig. 17(b) is used which is closer to the parametrization of the single beam. Note that there is still a difference due to the continuity in the middle of the knot span. The polynomial degree has no influence on the quality of the result if there are many elements as can be seen in Figs. 18(b) and 18(d). The of-leveling of the error in Figs. 18(a) and 18(c) is therefore not related to the coupling scheme but the discretization. The quality of the result reaches its optimum for a penalty factor α rot = 109 . Fig. 19 shows the convergence behavior for different polynomial degrees p and numbers of elements per beam with penalty factor α rot = 1010 . Note that the discretization of Fig. 17(a) was used for Fig. 19 since this parametrization is the outcome of the “real” modeling in common CAD systems. rot

4.3. Cross junction The next example should check if several beams can be coupled to one point and if the forces are transferred correctly. Non-matching tangents are a further property of this example to be verified. Therefore four beams assembled as a cross with one angular beam are modeled as shown in Fig. 20. Note that the coupling is not restricted to orthogonal configurations. One beam is clamped on the end and the skew one has a tip load. The clamped beam is the master beam for all three coupling conditions. The beams without boundary conditions should ideally remain straight and only be translated and rotated without internal forces. Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

15

Fig. 16. Relative error w.r.t. the single-beam-solution at the tip for: (a) u z for 4 elements per beam, (b) u z for 16 elements per beam, (c) φ for 4 elements per beam, and (d) φ for 16 elements per beam.

Fig. 17. Different discretizations of the coupled circular arch: (a) Control points with respective weights of the two curves, which are constructed by two 45◦ circular segments (dark blue), (b) Control points with respective weights of the two curves, which are constructed by knot insertion into a 90◦ circular segment (lighter blue) and (c) refined control points and Gauss integration points for 4 non-zero knot spans and both discretizations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

16

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 18. Relative error w.r.t. the single-beam-solution and the discretization of Fig. 17(b) at the tip for: (a) u z for 4 elements per beam, (b) u z for 16 elements per beam, (c) φ for 4 elements per beam, and (d) φ for 16 elements per beam.

As can be observed in the deformation plot with respective internal forces, the results match the expected behavior. The unbounded beams remain straight and no internal forces evolve. The internal forces of the other two beams are in equilibrium with the load and the support forces respectively (see Fig. 21). The largest difference ∆φi j measured as angle between two base vectors ti and t j at the coupling point is 6.66e−8rad, which is sufficiently small. ∆φi j = arccos(tim · tsj ) − arccos(Tim · Tis )

(42)

4.4. Mainspring The mainspring benchmark example was chosen in order to test the capability to undergo large rotations. A cantilever with the length L = 20π is again modeled by two beams with flat rectangular cross sections. Every center line is modeled by a B-Spline curve with p = 6 and 10 non-zero knot spans. It is loaded by an incrementally L increased moment M = 4π E I at the tip, which ideally should result in two circles with radius R = EMI = 4π = 5. L Note that the moment cannot be applied directly to the system since the element formulation of the beam is displacement-based and the nodal forces are changing during the simulation. Yet, it can be incorporated by the external work which generates additional terms in the stiffness matrix. The penalty factors for displacement and rotation are both set to α = 107 . The set-up and the corresponding deformation is shown in Fig. 22(a). The comparison of the tip displacements to the analytic reference solution in Fig. 22(b) shows an agreement over the whole load path. Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

17

Fig. 19. Numerical results for the quarter circle cantilever under tip load: (a) u z at the tip, (b) φ at the tip, (c) Mn at the support, and (d) Mt at the support.

Fig. 20. Cantilevering skew cross modeled by four coupled beams.

Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

18

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 21. Inner force of the cantilevering cross: (a) Normal force N , (b) bending moment Mn , (c) bending moment Mv , and (d) torsional moment Mt .

Fig. 22. Mainspring modeled by two coupled beams: (a) Set-up and several deformed states and (b) comparison of the tip displacement with the reference solution.

4.5. Demonstrator examples This section will show some more complex examples where a beam coupling formulation is needed. Each example shows a distinct feature for the simulation of coupled beams and illustrates the applicability. 4.5.1. Tripod The previous examples showed coupling of the ends of mostly straight beams. This is the first example where curved beams are used which are coupled inside their parameter space. Furthermore, this example shows the difference between full coupling with rotations and coupling of the displacements only. Therefore, circular segments are arranged in a tripod configuration (see Fig. 23). The circular segment is defined by the start and end point together with the tangent of the start. The NURBS curve is refined to a polynomial degree of p = 5 and 25 Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

19

Fig. 23. Set-up of a tripod with circular segments (shown with the unrefined control points) as beams subjected to a point load P.

Fig. 24. Deformed tripod: (a) coupling of displacements and (b) coupling of displacements and rotations.

non-zero knot spans. The lower ends are clamped, i.e. all rotations are blocked by the adapted weak coupling formulation for supports, whereas the upper ends are free. A penalty factor of α = 107 is applied for the weak conditions. The structure is loaded by a single load P at the intersection point, which results in a snap-through. Note that none of the base vectors, neither tangents nor cross sectional base vectors, are aligned. This stability problem results in a displacement of u z = −10.55721 of the intersection point for the fully coupled case and u z = −10.95770 for the sole coupling of displacements. Fig. 24 shows the deformed structure after the snap-through. 4.5.2. Gridshell A lattice gridshell is a perfect application example for advanced coupling conditions. Three by three beams are laid over each other and connected at the intersection points. In this case, the formulation for scissor joints is applied, i.e. only t2 is coupled. All ends are supported in z-direction. All beams are loaded with a line dead load in order to bend them into space as shown in Fig. 25(a). Fig. 25(b) shows the close-up of one outer joint. It can be observed that the base vectors t2 at the joint are matching within the tolerance introduced by the penalty approach. The angle between them, which corresponds to the error of the result, is φtm2 ts2 = 0.00198. The initial angle of φTm1 Ts1 = π2 between the tangents of the beam changed to φtm1 ts1 = 1.384. This behavior is perfectly in accordance with the modeling of a scissor joint, which can change its orientation during the analysis and stays local on the beam. Also the internal moments as shown in Fig. 26(a) and Fig. 26(b) are matching the expected result. The laths are bent around their weak axis and torsion occurs only in the outer four laths. Getting the normal forces is more demanding. The continuous basis functions are not able to represent discontinuities as introduced by the point coupling. In order to represent the correct behavior, the length of the tangent would have to change suddenly at the coupling, which is not possible. High oscillations are the result (cf. Fig. 26(c)). This phenomenon can also be observed with other higher-order basis functions. In Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

20

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 25. Gridshell with three by three laths connected by scissor joints: (a) Set-up with supports and loads and respective deformed structure and (b) close-up on one joint with base vectors t2 at every Gauss point of each beam and respective cross sections.

Fig. 26. Inner forces of the gridshell: (a) Bending moment Mn , (b) torsional moment Mt , (c) normal forces N and (d) normal forces N with locally refined curves.

order to restrain the propagation of these oscillations, local refinement, i.e. knot insertion around the discontinuity points, can be applied to the curves, which leads to a proper representation of the normal force as can be seen in Fig. 26(d). However, the displacement of the two different refinements differs only by 0.299%, which implies that the not locally refined curve can be used for the evaluation of the overall structural behavior. 4.5.3. Torus The last example is a torus which is constructed by a hexagonal mesh of individual beams. The beams have a circular cross section with r = 0.1. The Young’s modulus is E = 2.1 · 105 and the Poisson’s ratio is ν = 0. The center line of the torus has a radius of rc = 40. The center line is divided in n c = 28 modules. A module consists Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

21

Fig. 27. Torus with hexagonal beam grid: (a) Parameters for the initial geometry model, (b) load displacement curve and (c) initial and deformed model.

of five circles with ri = 10 around the center line. The circles are subdivided in n i = 13 nodes. These nodes are then used for the construction of the hexagonal grid (see Fig. 27(a)). The model consists of 2184 B-Spline curves, Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

22

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

each with a polynomial degree p = 3 and 4 control points. 4368 coupling conditions are included. The structure is supported in z-direction at four outer opposite points. One point is fully supported whereas the respective opposite point has an additional support in y-direction and a point load P = 7.5 in x-direction (see Fig. 27(c)). The simulation is stable and we can thus state that the coupling formulation is also applicable for large models with many degrees of freedom and a large number of complex coupling scenarios. The displacement of the loading point is given in Fig. 27(b), which reaches its maximum of u = 23.956322 under the load of P = 7.5. 5. Conclusions and outlook A coupling formulation based on two orthonormal triads was derived and implemented in the structural analysis code of the authors’ Chair. It is also accessible from the parametric environment Grasshopper in Rhino. It was applied and systematically tested in combination with a nonlinear isogeometric spatial Bernoulli beam. Several benchmark evaluations were conducted in order to quantify the performance of the proposed formulation. Furthermore, the formulation was applied in more complex examples in regard to the number of coupling points, location, orientation and coupling conditions. The newly proposed formulation can cope with: • • • •

large displacements and rotations coupling at non-interpolated locations arbitrary orientation of the cross sections at both sides, i.e. no aligned base vectors are required modular coupling conditions which only couple certain parts of the cross section, i.e. different types of hinges can be modeled

This weak coupling approach could also be used for other structural elements with constant triads or a combination of different elements, such as beam to shell coupling. Furthermore, other coupling conditions should be tried. A new condition which enforces the relation between two sides would also be beneficial for structural elements where the form of the triad is not preserved as e.g. for a shell. The proposed coupling conditions were implemented with a penalty approach. First investigations, as presented in this paper, showed that a penalty factor in the range of 107 − 1011 returns satisfying results and the concept is validated. However, the choice of the penalty factor with an estimation should be examined more thorough by a future investigation. Moreover, the imposition of the coupling condition by penalty method could be replaced by a Lagrange Multipliers approach. Another interesting point would be an investigation in case of non-linear material behavior. Since the formulation is purely geometric, it can also be used for structural elements with non-linear material law. In case of plasticity directly occurring in the coupling point, the penalty method has to be adapted due the linear elastic behavior of the spring. This could e.g. be realized by a strain-dependent penalty factor. Acknowledgments This work was supported by “Zentrales Innovationsprogramm Mittelstand” of the Bundesministerium f¨ur Wirtschaft und Energie (BMWi) as part of the project “Erforschung der mathematisch-methodischen Grundlagen zur Nutzung der Isogeometrischen Analyse in der Leichtbau-Planung” (ZF4066102BZ6). The support is gratefully acknowledged. References [1] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (39–41) (2005) 4135–4195, http://dx.doi.org/10.1016/j.cma.2004.10.008. [2] J.M. Kiendl, K.-U. Bletzinger, J. Linhard, R. Wüchner, Isogeometric shell analysis with Kirchhoff–Love elements, Comput. Methods Appl. Mech. Engrg. 198 (49–52) (2009) 3902–3914, http://dx.doi.org/10.1016/j.cma.2009.08.013. [3] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, Wiley a John Wiley and Sons Ltd. Publication, Chichester, 2009, http://dx.doi.org/10.1002/9780470749081. [4] D.J. Benson, S. Hartmann, Y. Bazilevs, M.-C. Hsu, T.J.R. Hughes, Blended isogeometric shells, Comput. Methods Appl. Mech. Engrg. 255 (2013) 133–146, http://dx.doi.org/10.1016/j.cma.2012.11.020. [5] J.M. Kiendl, Y. Bazilevs, M.-C. Hsu, R. Wüchner, K.-U. Bletzinger, The bending strip method for isogeometric analysis of Kirchhoff–Love shell structures comprised of multiple patches, Comput. Methods Appl. Mech. Engrg. 199 (37–40) (2010) 2403–2416, http://dx.doi.org/10.1016/j.cma.2010.03.029.

Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

23

[6] M. Breitenberger, A. Apostolatos, B. Philipp, R. Wüchner, K.-U. Bletzinger, Analysis in computer aided design: nonlinear isogeometric B-Rep analysis of shell structures, Comput. Methods Appl. Mech. Engrg. 284 (2015) 401–457, http://dx.doi.org/10.1016/j.cma.2014. 09.033. [7] A.J. Herrema, E.L. Johnson, D. Proserpio, M.C. Wu, J. Kiendl, M.-C. Hsu, Penalty coupling of non-matching isogeometric Kirchhoff– Love shell patches with application to composite wind turbine blades, Comput. Methods Appl. Mech. Engrg. 346 (2019) 810–840, http://dx.doi.org/10.1016/j.cma.2018.08.038. [8] A. Apostolatos, M. Breitenberger, R. Wüchner, K.-U. Bletzinger, Domain decomposition methods and Kirchhoff-Love shell multipatch coupling in isogeometric analysis, in: B. Jüttler, B. Simeon (Eds.), Isogeometric Analysis and Applications 2014, in: Lecture Notes in Computational Science and Engineering, vol. 107, Springer, Cham and Heidelberg and New York, 2015, pp. 73–101, http://dx.doi.org/10.1007/978-3-319-23315-4_4. [9] W. Dornisch, G. Vitucci, S. Klinkel, The weak substitution method - an application of the mortar method for patch coupling in NURBS-based isogeometric analysis, Internat. J. Numer. Methods Engrg. 103 (3) (2015) 205–234, http://dx.doi.org/10.1002/nme.4918. [10] T. Teschemacher, A.M. Bauer, T. Oberbichler, M. Breitenberger, R. Rossi, R. Wüchner, K.-U. Bletzinger, Realization of CAD-integrated shell simulation based on isogeometric B-Rep analysis, Adv. Model. Simul. Eng. Sci. 5 (1) (2018) 19, http://dx.doi.org/10.1186/s40323018-0109-4. [11] M. Ruess, D. Schillinger, A.I. Özcan, E. Rank, Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries, Comput. Methods Appl. Mech. Engrg. 269 (2014) 46–71, http://dx.doi.org/10.1016/j.cma.2013.10.009. [12] V.P. Nguyen, P. Kerfriden, M. Brino, S.P.A. Bordas, E. Bonisoli, Nitsche’s method for two and three dimensional NURBS patch coupling, Comput. Mech. 53 (6) (2014) 1163–1182, http://dx.doi.org/10.1007/s00466-013-0955-3. [13] X. Du, G. Zhao, W. Wang, Nitsche method for isogeometric analysis of Reissner–Mindlin plate with non-conforming multi-patches, Comput. Aided Geom. Design 35–36 (2015) 121–136, http://dx.doi.org/10.1016/j.cagd.2015.03.005. [14] Y. Guo, M. Ruess, Nitsche’s method for a coupling of isogeometric thin shells and blended shell structures, Comput. Methods Appl. Mech. Engrg. 284 (2015) 881–905, http://dx.doi.org/10.1016/j.cma.2014.11.014. [15] W. Jiang, C. Annavarapu, J.E. Dolbow, I. Harari, A robust Nitsche’s formulation for interface problems with spline-based finite elements, Internat. J. Numer. Methods Engrg. 104 (7) (2015) 676–696, http://dx.doi.org/10.1002/nme.4766. [16] Y. Guo, M. Ruess, D. Schillinger, A parameter-free variational coupling approach for trimmed isogeometric thin shells, Comput. Mech. 59 (4) (2017) 693–715, http://dx.doi.org/10.1007/s00466-016-1368-x. [17] Y. Guo, J. Heller, T.J. Hughes, M. Ruess, D. Schillinger, Variationally consistent isogeometric analysis of trimmed thin shells at finite deformations, based on the STEP exchange format, Comput. Methods Appl. Mech. Engrg. 336 (2018) 39–79, http: //dx.doi.org/10.1016/j.cma.2018.02.027. [18] S. Schuß, M. Dittmann, B. Wohlmuth, S. Klinkel, C. Hesch, Multi-patch isogeometric analysis for Kirchhoff–Love shell elements, Comput. Methods Appl. Mech. Engrg. 349 (2019) 91–116, http://dx.doi.org/10.1016/j.cma.2019.02.015. [19] E. Brivadis, A. Buffa, B. Wohlmuth, L. Wunderlich, Isogeometric mortar methods, Comput. Methods Appl. Mech. Engrg. 284 (2015) 292–319, http://dx.doi.org/10.1016/j.cma.2014.09.012. [20] L. de Lorenzis, P. Wriggers, G. Zavarise, A mortar formulation for 3D large deformation contact using NURBS-based isogeometric analysis and the augmented Lagrangian method, Comput. Mech. 49 (1) (2012) 1–20, http://dx.doi.org/10.1007/s00466-011-0623-4. [21] C. Hesch, P. Betsch, Isogeometric analysis and domain decomposition methods, Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 104–112, http://dx.doi.org/10.1016/j.cma.2011.12.003. [22] F. Auricchio, L. Beirão da Veiga, J.M. Kiendl, C. Lovadina, A. Reali, Locking-free isogeometric collocation methods for spatial timoshenko rods, Comput. Methods Appl. Mech. Engrg. 263 (2013) 113–126, http://dx.doi.org/10.1016/j.cma.2013.03.009. [23] L. Greco, M. Cuomo, B-Spline interpolation of Kirchhoff-Love space rods, Comput. Methods Appl. Mech. Engrg. 256 (2013) 251–269, http://dx.doi.org/10.1016/j.cma.2012.11.017. [24] G. Zhang, R. Alberdi, K. Khandelwal, Analysis of three-dimensional curved beams using isogeometric approach, Eng. Struct. 117 (2016) 560–574, http://dx.doi.org/10.1016/j.engstruct.2016.03.035. [25] S.B. Raknes, X. Deng, Y. Bazilevs, D.J. Benson, K.M. Mathisen, T. Kvamsdal, Isogeometric rotation-free bending-stabilized cables: statics, dynamics, bending strips and coupling with shells, Comput. Methods Appl. Mech. Engrg. 263 (2013) 127–143, http://dx.doi.org/10.1016/j.cma.2013.05.005. [26] A.M. Bauer, M. Breitenberger, B. Philipp, R. Wüchner, K.-U. Bletzinger, Nonlinear isogeometric spatial Bernoulli beam, Comput. Methods Appl. Mech. Engrg. 303 (2016) 101–127, http://dx.doi.org/10.1016/j.cma.2015.12.027. [27] E. Marino, Locking-free isogeometric collocation formulation for three-dimensional geometrically exact shear-deformable beams with arbitrary initial curvature, Comput. Methods Appl. Mech. Engrg. 324 (2017) 546–572, http://dx.doi.org/10.1016/j.cma.2017.06.031. [28] O. Weeger, B. Narayanan, L. De Lorenzis, J. Kiendl, M.L. Dunn, An isogeometric collocation method for frictionless contact of cosserat rods, Comput. Methods Appl. Mech. Engrg. 321 (2017) 361–382, http://dx.doi.org/10.1016/j.cma.2017.04.014. [29] O. Weeger, S.-K. Yeung, M.L. Dunn, Fully isogeometric modeling and analysis of nonlinear 3D beams with spatially varying geometric and material parameters, Comput. Methods Appl. Mech. Engrg. 342 (2018) 95–115, http://dx.doi.org/10.1016/j.cma.2018.07.033. [30] L. Greco, M. Cuomo, An implicit G1 multi patch B-spline interpolation for Kirchhoff–Love space rod, Comput. Methods Appl. Mech. Engrg. 269 (2014) 173–197, http://dx.doi.org/10.1016/j.cma.2013.09.018. [31] W. Wunderlich, W.D. Pilkey, Mechanics of Structures: Variational and Computational Methods, second ed., CRC Press, Boca Raton, FL, 2003. [32] M. Breitenberger, CAD-Integrated design and analysis of shell structures (Dissertation), Technische Universität München, München, 2016.

Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.

24

A.M. Bauer, R. W¨uchner and K.-U. Bletzinger / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

[33] A. Apostolatos, R. Schmidt, R. Wüchner, K.-U. Bletzinger, A nitsche-type formulation and comparison of the most common domain decomposition methods in isogeometric analysis, Internat. J. Numer. Methods Engrg. 97 (7) (2014) 473–504, http://dx.doi.org/10.1002/ nme.4568. [34] G. Lumpe, V. Gensichen, Evaluierung der linearen und nichtlinearen Stabstatik in Theorie und Software: Prüfbeispiele, Fehlerursachen, genaue Theorie, Bauingenieur-Praxis, Ernst W. + Sohn Verlag, Berlin, 2014.

Please cite this article as: A.M. Bauer, R. W¨uchner and K.-U. Bletzinger, Weak coupling of nonlinear isogeometric spatial Bernoulli beams, Computer Methods in Applied Mechanics and Engineering (2019) 112747, https://doi.org/10.1016/j.cma.2019.112747.