Consistent coupling of positions and rotations for embedding 1D Cosserat beams into 3D solid volumes

This article proposes a mortar type finite element formulation for consistently embedding curved, slender beams, i.e. 1D Cosserat continua, into 3D solid volumes. A consistent 1D-3D coupling scheme for this problem type is proposed, which enforces both positional and rotational constraints. Since Boltzmann continua exhibit no inherent rotational degrees of freedom, suitable definitions of orthonormal triads are investigated that are representative for the orientation of material directions in the 3D solid. The rotation tensor defined by the polar decomposition of the deformation gradient is demonstrated to represent these material directions in a L2-optimal manner. Subsequently, objective rotational coupling constraints between beam and solid are formulated and enforced in a variationally consistent framework. Eventually, finite element discretization of all primary fields results in an embedded mortar formulation for rotational and translational constraint enforcement. Based on carefully chosen numerical test cases, the proposed scheme is demonstrated to exhibit a consistent spatial convergence behavior and to offer the up-scaling potential for studying real-life engineering applications such as fiber-reinforced composite materials.


Introduction
Embedding fibers or beams, i.e. solid bodies that can mechanically be modeled as dimensionally reduced 1D structures since one spatial dimension is much larger than the other two, into a 3D matrix material is a common approach to enhance the mechanical properties of a structure. Fiberreinforced structures can be found in many different fields, e.g. in form of steel reinforcements within concrete structures, lightweight fiber-reinforced composite materials based on carbon, glass or polymer fibers in a plastic matrix, or additively manufactured components allowing for a very flexible and locally controlled reinforcement of plastic, metal, ceramic or concrete matrix materials [32,33,42]. At a different length scale, fiber embeddings play a key role for essential processes in countless biological systems, e.g. in the form of embedded networks (e.g. cytoskeleton, extracellular matrix, mucus) or bundles (e.g. muscle, tendon, ligament) [2,21,31,41]. Most of these applications are characterized by geometrically complex embeddings of arbitrarily oriented, slender and potentially curved fibers. Computational models predicting the response of such reinforced structures are essential for a time-and cost-efficient design and development of technical products, but also to gain fundamental understanding of biological systems at length scales that are not accessible via experiments. In the context of computational modeling, as considered in the following, the embedded 1D structures will be referred to as fibers or beams, respectively, and the 3D matrix as solid.
One common modeling approach for this physical beamto-solid volume coupling problem is based on homogenized, anisotropic material models for the combined fiber-matrix structure [1,70]. This widely used approach is appealing since, e.g. no additional degrees of freedom are required to model individual fibers, and existing simulation tools can be used as long as they support anisotropic material laws. However, such models cannot give detailed information about the interactions between fibers and surrounding matrix as, e.g. required to study mechanisms of failure. Moreover, the fiber distribution in the solid has to be sufficiently homogeneous and a separation of scales is required, i.e. the fiber size has to be sufficiently small as compared to the smallest dimension of the overall structure. Eventually, when modeling new fiber arrangements, the homogenization step inherent to these continuum models requires sub-scale information, e.g. provided by a model with resolved fiber geometries.
Another modeling approach consists of fully describing the fibers and surrounding solid material as 3D continua. This leads to a surface-to-surface coupling problem at the 2D interface between fiber surface and surrounding solid. In the context of the finite element method, these surfaces can be tied together by either applying fiber and solid discretizations that are conforming at the shared interface or via interface coupling schemes accounting for non-matching meshes, such as the mortar method [45,[47][48][49]. Alternatively, extended finite element methods (XFEM) [40] or immersed finite element methods [29,55] can be used to represent 2D fiber surfaces embedded in an entirely independent background solid mesh. While such fully resolved modeling approaches allow to study local effects with high spatial resolution, the significant computational effort associated with these models prohibits their usage for large-scale systems with a large number of slender fibers.
The class of applications considered here typically involves very slender fibers. In this regime it is well justified, and highly efficient from a computational point of view, to model individual fibers as beams, e.g. based on the geometrically exact beam theory [9,12,14,24,30,36,38,51,52,[59][60][61], which is known to combine high model accuracy and computational efficiency [5,53]. Based on the fundamental kinematic assumption of undeformable cross-sections, such beam models can be identified as 1D Cosserat continua with six degrees of freedom defined at each centerline point to describe the cross-section position (three translational degrees of freedom) and orientation (three rotational degrees of freedom). Thus, the problem of beams embedded in a 3D solid volume can be classified as mixed-dimensional 1D-3D coupling problem between 1D Cosserat continua and a 3D Boltzmann continuum. A variety of 1D-3D coupling approaches exist in the literature, however most of them involve truss / string models, i.e. 1D structural models account only for internal elastic energy contributions from axial tension, e.g. [13,17,20,25,26,43,50]. Work on the 1D-3D coupling between beams, i.e. full Cosserat continua, and solids is much rarer. In [16], collocation along the beam centerline is applied to couple beams with a surrounding solid material. A mortar-type coupling approach is proposed in the authors' previous work [63], where a Lagrange multiplier field is defined along the beam centerline to weakly enforce the coupling constraints. The 1D-3D coupling between beams and a surrounding fluid field, as relevant for fluid-structure-interaction (FSI) problems, has been considered in some recent contributions [22,68].
All the aforementioned 1D-3D beam-to-solid coupling schemes have in common that only the beam centerline positions, but not the cross-section orientations, are coupled to the solid, which will be denoted as translational 1D-3D coupling. In such models, an embedded fiber can still perform local twist/torsional rotations, i.e. cross-section rotations with respect to its centerline tangent vector, relative to the solid. While this simplified coupling procedure can reasonably describe the mechanics of certain problem classes where such relative rotations will rarely influence the global system response, e.g. embedding of straight fibers with circular cross-section shape, for most practical applications a more realistic description of the physical problem requires to also couple the rotations of beam and solid.
In a very recent approach by [27] the full 1D-3D coupling problem involving positions and rotations has been addressed for the first time. The coupling of the two directors spanning the (undeformable) beam cross-section with the underlying solid continuum together with the coupling of the cross-section centroids results in a total of nine coupling constraints. One specific focus of this interesting contribution lies on a static condensation strategy, which allows to eliminate the associated Lagrange multipliers and the beam balance equations from the final, discrete system of equations. The requirement of a C 1 -continuous spatial discretization of the solid domain, as resulting from the proposed condensation strategy, is satisfied by employing NURBS-based test and trial functions.
The present work proposes a full 1D-3D coupling approach based on only six, i.e. three translational and three rotational, coupling constraints between the cross-sections of 1D beams, modeled according to the geometrically exact beam theory, and a 3D solid. The finite element method is employed for spatial discretization of all relevant fields. Consistently deriving the full 1D-3D coupling on the beam centerline from a 2D-3D coupling formulation on the beam surface via a first-order Taylor series expansion of the solid displacement field would require to fully couple the two orthonormal directors spanning the (undeformable) beam cross-section with the (in-plane projection of the) solid deformation gradient evaluated at the cross-section centroid position. It is demonstrated that such an approach, which suppresses all in-plane deformation modes of the solid at the coupling point, might result in severe locking effects in the practically relevant regime of coarse solid mesh sizes. Therefore, as main scientific contribution of this work, different definitions of orthonormal triads are proposed that are representative for the orientation of material directions of the 3D continuum in an average sense, without additionally constraining in-plane deformation modes when coupled to the beam cross-section. It is shown that the rotation tensor defined by the polar decomposition of the (in-plane projection of the) deformation gradient appears as a natural choice for this purpose, which even represents the average orientation of material directions of the 3D continuum in a L 2 -optimal manner. Moreover, several alternative solid triad definitions are investigated that potentially allow for a more efficient numerical evaluation.
Once these solid triads have been defined, objective (i.e. frame-invariant) rotational coupling constraints in the form of relative rotations are formulated for each pair of triads representing the beam and solid orientation. Their variationally consistent enforcement either based on a penalty potential or a Lagrange multiplier potential, with an associated Lagrange multiplier field representing a distributed coupling moment along the beam centerline, is shown. Eventually, finite element discretization of the Lagrange multiplier and relative rotation vector field along the beam centerline results in an embedded mortar-type formulation for rotational constraint enforcement. In combination with a previously developed mortar-type formulation (BTS-TRANS) for translational 1D-3D coupling [63], this results in a full 1D-3D coupling approach denoted as full beam-to-solid volume coupling scheme (BTS-FULL). Finite element discretization of the solid and the embedded (potentially curved) beams inevitably results in non-matching meshes, which underlines the importance of a consistently embedded mortar-type formulation as proposed in this work. Based on elementary numerical test cases, it is demonstrated that a consistent spatial convergence behavior can be achieved and potential locking effects can be avoided if the proposed BTS-FULL scheme is combined with a suitable solid triad definition. Eventually, real-life engineering applications are considered to illustrate the importance of consistently coupling both translational and rotational degrees of freedom as well as the upscaling potential of the proposed formulation to study complex mechanical systems such as fiber-reinforced composite materials, containing a large number of curved, slender fibers with arbitrary orientation embedded in a matrix material.
The remainder of this work is organized as follows: In Section 2, we state the fundamental modeling assumptions of the proposed BTS-FULL scheme. Specifically, the importance of enforcing both rotational and translational coupling conditions is demonstrated, and the general implications of a 1D-3D coupling approach are discussed. In Section 3, we give a short summary of the theory of large rotations as required to formulate rotational coupling conditions. In Section 4, the governing equations for the solid and beam domains are presented, and objective rotational coupling constraints are defined and enforced in a variationally consistent manner, either based on a penalty or a Lagrange multiplier potential. In Section 5, we propose different definitions of orthonormal triads that are suitable to represent the orientation of solid material directions in an average sense. In Section 6, discretization of the coupling conditions based on the finite element method is considered, once in a Gauss point-to-segment manner and once as mortar-type approach along with a weighted penalty regularization. Finally, numerical examples, carefully selected to verify different aspects of the proposed formulation, are presented in Section 7.

Motivation and modeling assumptions
In Section 2.1, the main modeling assumptions generally underlying 1D-3D coupling schemes will be discussed. Subsequently, in Section 2.2, the importance of a full position and rotation coupling (BTS-FULL) will be motivated for general application scenarios, and special cases will be discussed, where also a purely translational coupling (BTS-TRANS) can be considered as reasonable approximation.

Modeling assumptions underlying the 1D-3D coupling
The considered class of 1D-3D coupling schemes is based on the assumption that the fiber material is stiff compared to the solid material, and local fiber cross-section dimensions are small compared to the global solid dimensions. Thus, the solid may be discretized without subtracting the fiber volume, formally resulting in overlapping solid and fiber domains. While consistent 2D-3D coupling on the fiber surface would allow for high-resolution stress field predictions in the direct vicinity of the 2D fiber-solid interface, such approaches require an evaluation of coupling constraints on a Fig. 1 Plane coupling problem of a single fiber cross-section with a solid finite element mesh -full 1D-3D coupling (left) vs. 2D-3D coupling (right).
2D interface and a sufficient discretization resolution of the solid with mesh sizes smaller than the fiber cross-section dimensions, thus in large parts deteriorating the advantages provided by a reduced dimensional description of the fibers.
In truly 1D-3D coupling approaches, the coupling conditions are exclusively defined along the beam centerline, thus preserving the computational advantages of the dimensionally reduced beam models. Of course, such approaches inevitably introduce a modeling error as compared to the 2D-3D coupling, i.e. the surface tractions on the 2D beam-solid interface are approximated by localized resultant line forces and moments acting on the beam centerline. This has a significant impact on the analytical solution of the problem, as line loads acting on a 3D continuum result in singular stress and displacement fields, cf. [19,44,66]. Thus, convergence of the 1D-3D solution towards the 2D-3D solution is not expected. However, in the realm of the envisioned applications, we are rather interested in global system responses than in local stress distributions in the direct vicinity of the fibers. Thus, practically relevant solid element sizes are considered that are larger than the fiber cross-section dimension. In this regime of mesh resolutions, this inherent modeling error of 1D-3D approaches can typically be neglected.
To verify this statement, consider a plane problem of a beam cross-section, loaded with a moment, that is coupled to a solid finite element as depicted in Figure 1. As long as the cross-section diameter is smaller than the solid finite element mesh size, the resulting discrete nodal forces F S acting on the solid are independent of the employed coupling approach, i.e. either 1D-3D coupling with associated coupling moment M (Figure 1, left) or 2D-3D coupling with associated coupling surface load τ (Figure 1, right). Obviously, this is an idealized setting, but it illustrates that 1D-3D coupling approaches can be considered as valid models for solid mesh sizes larger than the cross-section diameter, which will also be verified in Section 7. For a more detailed discussion of this topic the interested reader is referred to our previous publication [63], specifically to Figure 15 in [63], which depicts an analogous scenario for the coupling of translational degrees of freedom.

Motivation for full translational and rotational coupling
To differentiate the scope of validity of the proposed BTS-FULL scheme (coupling of positions and rotations) and of existing BTS-TRANS schemes (coupling of positions only), two application scenarios are discussed.
As first scenario, systems are considered (i) that contain only transversely isotropic fibers (e.g. circular cross-section shape and initially straight) and (ii) whose global system response is dominated by the axial and bending stiffness of the fibers, i.e. the torsional contribution is negligible. As demonstrated in [63], BTS-TRANS schemes can be considered as a reasonable mechanical model in this case, since local (twist/torsional) rotations of the fibers with respect to their straight axes will rarely influence the global system response. Torsion-free beam models [37] represent an elegant mechanical description of the fibers for such applications.
As second scenario, systems are considered that contain transversely anisotropic fibers (e.g. non-circular crosssection shape or initially curved). First, it is clear that twist rotations of the fiber cross-sections with respect to the centerline tangent (even if not possible in their simplest form as rigid body rotations) will change the global system response, since such fibers exhibit distinct directions of maximal/minimal bending stiffness or initial curvature. Second, due to the inherent two-way coupling of bending and torsion in initially curved beams [37], bending deformation will inevitably induce torsion in such application scenarios, i.e. the global system stiffness is approximated as too soft if these torsional rotations are not transferred to the matrix by a proper coupling scheme. Thus, a unique and consistent mechanical solution for this scenario can only be guaranteed by BTS-FULL schemes.
Remark 2.1 In fact, both aforementioned application scenarios might lead to non-unique static solutions if neglecting the rotational coupling. However, for transversely isotropic fibers the non-uniqueness only occurs at the local fiber level, i.e. the twist orientation of the fibers is not uniquely defined, which does not influence the global system response. The locally non-unique fiber orientation is typically only an issue from a numerical point of view (e.g. linear solvers), and can be effectively circumvented by employing, e.g. torsion-free beam models not exhibiting the relevant rotational degrees of freedom. For transversely anisotropic fibers, such local twist rotations will change the global system response. This gives rise to non-unique static solutions on the global level and, thus, has significant implications from a physical point of view.

Large rotations
This section gives a brief overview on the mathematical treatment of finite rotations as required for the formulation of rotational coupling constraints. For a more comprehesive treatment of this topic, the interested reader is referred to [8,12,24,38,52,60]. Let us consider a rotation tensor where SO 3 is the special orthogonal group and the base vectors g i form an orthonormal triad, that maps the Cartesian basis vectors e i onto g i . In the following, a rotation pseudovector ψ is used for its parametrization, i.e. Λ = Λ(ψ). The rotation vector describes a rotation by an angle ψ = ψ around the rotation axis e ψ = ψ/ ψ . The parametrization can be given by the well-known Rodrigues formula [3] Λ(ψ) = exp (S (ψ)) where exp(·) is the exponential map. Furthermore, S ∈ so 3 is a skew-symmetric tensor, where so 3 represents the set of skew-symmetric tensors with The inverse of the Rodrigues formula (2), i.e. the rotation vector as a function of the rotation tensor, will be denoted as ψ(Λ) = rv(Λ) in the remainder of this work. In practice, Spurrier's algorithm [62] can be used for the extraction of the rotation vector. Two triads Λ 1 (ψ 1 ) and Λ 2 (ψ 2 ), with their respective rotation vectors ψ 1 and ψ 2 , can be related by the relative rotation Λ 21 (ψ 21 ). The relative rotation is given by with the identity Λ T = Λ −1 for all elements of SO 3 . Thus, the (non-additive) rotation vector ψ 21 = rv (Λ 21 ) = ψ 2 − ψ 1 describes the relative rotation between Λ 1 and Λ 2 .
In a next step, the variation of the rotation tensor shall be considered, which can be expressed either by an (infinitesimal) additive variation δψ of the rotation vector or by a (infinitesimal) multiplicative rotation variation δθ, which is also denoted as spin vector: With the relation above and the definition of S, the variations of the triad basis vectors δg i read The (infinitesimal) additive and multiplicative rotation vector variations can be related according to where the transformation matrix T (ψ) [61] is defined as In [34], the objective variation δ o of a spatial quantity defined in a moving frame Λ 1 is defined as the difference between the total variation and the variation of the base vectors of the moving frame. In the context of rotational coupling constraints this will be required when expressing the objective variation of a relative rotation vector ψ 21 : For a detailed derivation of this expression for the objective variation the interested reader is referred to [34].
Remark 3.1 Via right-multiplication of (8) with the rotation vector ψ it can easily be shown that ψ is an eigenvector (with eigenvalue 1) of T and also of T T , i.e. T ψ = ψ and T T ψ = ψ. This property will be beneficial for derivations presented in subsequent sections. Every vector parallel to ψ is also an eigenvector of T . This can be interpreted in a geometrical way: If the additive increment δψ to a rotation vector ψ is parallel to the rotation vector, i.e. δψ = δψe ψ and ψ = ψe ψ , the resulting compound rotation ψ + δψ = (ψ + δψ) e ψ is still defined around the rotation axis e ψ . In this case, the rotation increment is a plane rotation relative to Λ(ψ), and the multiplicative and additive rotational increments are equal to each other, δψ = δθ.

Remark 3.2
In addition to Λ, also the symbol R will be used in the following to represent rotation tensors.

Problem formulation
We consider a 3D finite deformation full beam-to-solid volume coupling problem (BTS-FULL) as shown in Figure 2. All quantities are refereed to a Cartesian frame e 1 , e 2 , e 3 . For simplicity, we focus on quasi-static problems in this work, while the presented BTS-FULL method is directly applicable to dynamic problems as well. The principle of virtual work serves as basis for the proposed finite element formulation.
Contributions to the total virtual work of the system can be split into solid, beam and coupling terms, where the solid and beam terms are independent of the coupling constraints, i.e. well-established modeling and discretization techniques can be used for these single fields, cf. [63].

Solid formulation
The solid body is modeled as a 3D Boltzmann continuum, defined by its domain Ω S,0 ⊂ R 3 in the reference configuration, with boundary ∂Ω S,0 . Throughout this work, the subscript (·) 0 indicates a quantity in the reference configuration. A solid material point can be identified by its reference position X S ∈ R 3 . The current position x S ∈ R 3 relates to X S via the displacement field u S ∈ R 3 , i.e.
The domain and surface of the solid in the deformed configuration are Ω S and ∂Ω S , respectively. Virtual work contributions δW S of the solid are given by where δ denotes the (total) variation of a quantity, S ∈ R 3×3 is the second Piola-Kirchhoff stress tensor, E ∈ R 3×3 is the work-conjugated Green-Lagrange strain tensor,b ∈ R 3 is the body load vector andt ∈ R 3 are surface tractions on the Neumann boundary Γ σ ⊂ ∂Ω S,0 . The Green-Lagrange strain tensor is defined as For the compressible or nearly incompressible solid material, we assume existence of a hyperelastic strain energy function Ψ (E), which allows to determine the second Piola-Kirchhoff stress tensor according to S = ∂Ψ (E) ∂E .

Geometrically exact beam theory
The beams are modeled as 1D Cosserat continua embedded in 3D space based on the geometrically exact Simo-Reissner beam theory. Thus, each beam cross-section is described by six degrees of freedom, namely three positional and three rotational degrees of freedom. This results in six deformation modes of the beam: axial tension, bending (2×), shear (2×) and torsion. The cross-section centroids are connected by the centerline curve r(s) ∈ R 3 , where s ∈ [0, L] =: Ω B,0 ⊂ R is the arc-length coordinate along the beam centerline Ω B,0 in the reference configuration, and L the corresponding reference length. The displacement of the beam centerline u B (s) ∈ R 3 relates the reference position r 0 to the current position r via r(s) = r 0 (s) + u B (s).
The orientation of the beam cross-section field is described by the following field of right-handed orthonormal triads Λ B (s) := [g B1 (s), g B2 (s), g B3 (s)] = Λ B (ψ B (s)) ∈ SO 3 , which maps the global Cartesian basis vectors e i onto the local cross-section basis vectors g Bi (s) = Λ B e i for i = 1, 2, 3. Therein, ψ B ∈ R 3 is the rotation pseudo-vector chosen as parametrization for the triad. Moreover, the triad field in the reference configuration is denoted as Λ B,0 (s) := [g B1,0 (s), g B2,0 (s), g B3,0 (s)] = Λ B,0 (ψ B,0 (s)), and the relative rotation between the triads in reference and current configuration is denoted as R B := Λ B Λ T B,0 . According to the fundamental kinematic assumption of undeformable cross-sections, the position of an arbitrary material point within the beam cross-section either in the reference or in the current configuration can be expressed as follows: x B (s, α, β) = r(s) + αg B2 (s) + βg B3 (s), where α and β represent in-plane coordinates. Based on a hyperelastic stored-energy function according to can be derived. Here, Γ ∈ R 3 is a material deformation measure representing axial tension and shear, Ω ∈ R 3 is a material deformation measure representing torsion and bending, and C F ∈ R 3×3 and C M ∈ R 3×3 are cross-section constitutive matrices. Eventually, the beam contributions to the weak form are given by with the virtual work δW B ext of external forces and moments.

Full beam-to-solid volume coupling (BTS-FULL)
In the proposed BTS-FULL method, the pointwise six degrees of freedom associated with the beam centerline positions and cross-section triads are coupled to the surrounding solid, i.e.
Herein, Γ c = Ω S,0 ∩ Ω B,0 is the one-dimensional coupling domain between the beam centerline and the solid volume, i.e. the part of the beam centerline that lies within the solid. The rotational coupling between beam cross-section and solid as presented in this section is in close analogy to the generalized cross-section interaction laws proposed in [34]. The rotation vector ψ SB describes the relative rotation between a beam cross-section triad Λ B and a corresponding triad Λ S associated with the current solid configuration, Opposite to Λ B , which is well defined along the beam centerline, there is no obvious or unique definition for Λ S in the solid domain. In Section 5, different definitions of the solid triad Λ S are presented and investigated. However, for the derivation of the coupling equations, it is sufficient to assume the general form Λ S = Λ S (F ), i.e. formulating the solid triad as a general function of the solid deformation gradient in the current configuration. The formulation of the constraint equations along the beam centerline brings about an advantageous property of the BTS-FULL method: the translational (18) and rotational (19) coupling constraints are completely decoupled. Therefore, the rotational coupling equations (19) can be interpreted as a direct extension to the BTS-TRANS method, which only couples the beam centerline positions to the solid as derived and thoroughly discussed in [63]. In what follows, two different constraint enforcement strategies for the rotational coupling conditions will be presented.
Remark 4.1 In Section 7, we compare the BTS-FULL method to a full 2D-3D coupling approach that enforces constraints at the 2D beam-solid interface. The governing equations, as well as the discretized coupling terms for this 2D-3D coupling scheme are stated in Appendices A.2 and A.3.

Penalty potential
We consider a quadratic space-continuous penalty potential between beam cross-section triads and solid triads defined along the beam centerline: with the cross-section coupling potential π θ = π θ (s) and the symmetric penalty tensor c ∈ R 3×3 . Variation of the penalty potential leads to the following contribution to the weak form: Therein, δ o ψ SB is the objective variation of the rotation vector ψ SB . Making use of (9), the variation of the total potential becomes, cf. [34], where δθ S and δθ B are multiplicative variations associated with the solid and beam triad, respectively. Here, we consider penalty tensors of the form c = θ I with a scalar penalty parameter θ ∈ R + with physical unit Nm/m. With this definition and the identity T T (ψ)ψ = ψ (cf. Remark 3.1) the variation of the penalty potential simplifies to It is well-known from the geometrically exact beam theory that the (multiplicative) virtual rotations δθ B are workconjugated to the moment stress resultants. Therefore, θ ψ SB can be directly interpreted as the (negative) coupling moment acting on the beam cross-section.

Lagrange multiplier potential
Alternatively, the Lagrange multiplier method can be employed to impose the rotational coupling constraints. A Lagrange multiplier field λ θ = λ θ (s) ∈ R 3 is therefore defined on the coupling curve Γ c . For now, this field is a purely mathematical construct in the sense of generalized coupling forces associated with the coupling conditions (19). The Lagrange multiplier potential for the rotational coupling is Variation of the Lagrange multiplier potential again leads to a constraint contribution to the weak form, i.e.
Therein, δW λ θ and δW C θ are the variational form of the coupling constraints and the virtual work of the generalized coupling forces λ θ , respectively. With (9) the virtual work of the generalized coupling forces becomes Since the multiplicative rotation variations δθ B are workconjugated to the moment stress resultants of the beam, the term −T T (ψ SB )λ θ can be interpreted as a distributed coupling moment acting along the beam centerline.

Remark 4.2
For a vanishing relative rotation ψ 21 = 0, as enforced in the space-continuous problem setting according to (19), the identity −T T (ψ SB ) = I holds true and the rotational Lagrange multipliers exactly represent the coupling moments along the beam centerline. However, for the discretized problem this is only an approximation.

Objectivity of full beam-to-solid volume coupling
As indicated above, the solid triad field depends on the solid deformation gradient F . It can easily be shown, that the presented solid triad definitions STR-POL, STR-AVG and STR-ORT, in Section 5 are objective with respect to an arbitrary rigid body rotation R * ∈ SO 3 , i.e.
The geometrically exact beam model employed in this contribution is also objective [36,38], i.e.
Equations (28) and (29) inserted into the definition of the relative beam-to-solid rotation vector according to (20) gives the rotated relative rotation vector, where the identity rv(R * ΛR * T ) = R * rv(Λ) has been used. Thus, the rotational coupling conditions (19) in combination with the proposed solid triad definitions and the employed geometrically exact beam models are objective. As shown in [34], in this case also an associated penalty potential of type (21) or an associated Lagrange multiplier potential of type (25) is objective. The previous considerations show objectivity of the proposed (space-continuous) 1D-3D coupling approaches. However, in the realm of the finite element method, cf. Section 6, it is important to demonstrate that objectivity is preserved also in the discrete problem setting. It is well known that the discretized deformation gradient, as required for the definition of solid triads, is objective as long as standard discretization schemes (e.g. via Lagrange polynomials) are applied to the displacement field of the solid. Also the employed beam finite element formulation based on the geometrically exact beam theory is objective, even though this topic is not trivial and the interested reader is referred to [36,38]. Therefore, it can be concluded that the proposed 1D-3D coupling schemes are objective for the space-continuous as well as for the spatially discretized problem setting.

Remark 4.3
Objectivity is the main reason for formulating the rotational coupling constraints (19) based on the relative rotation vector, i.e. ψ SB = 0, cf. [34]. As alternative choice for the rotational coupling constraints the difference between the beam and solid triad rotation vectors, i.e. ψ B − ψ S = 0, could be considered. However, such coupling constraints would result in a non-objective coupling formulation [34].

Definition of solid triad field
One of the main aspects of the present work is the definition of a suitable right-handed orthonormal triad field Λ S in the solid, which is required for the coupling constraint (19). This is by no means a straightforward choice, and different triad definitions will lead to different properties of the resulting numerical coupling scheme. In the following, a brief motivation will be given for the concept of solid triads before different solid triad definitions will be proposed.

Motivation of the solid triad concept
If the embedded beam is considered as a 3D body, a consistent 2D-3D coupling constraint between the 2D beam surface and the surrounding 3D solid can be formulated as Therein, Γ 2D-3D is the 2D-3D coupling surface, i.e. the part of the beam surface that lies withing the solid volume. In the following, let X S,r denote the line of material solid points that coincide with the beam centerline in the reference configuration, i.e. X S,r = r 0 . Furthermore, the orthonormal triad Λ S,0 = [g S1,0 , g S2,0 , g S3,0 ] shall represent material directions of the solid that coincide with the beam triad in the reference configuration according to The corresponding quantities in the deformed configuration are denoted as x S,r and Λ S . Let us now expand the position field in the solid as Taylor series around x S,r , i.e.
where F is the deformation gradient of the solid according to (12). The 1D-3D coupling strategy underlying the proposed BTS-FULL scheme relies on the basic assumption of slender beams, i.e. R L, where R is a characteristic cross-section dimension (e.g. the radius of circular cross-sections). This assumption allows to truncate the Taylor series after the linear term as long as small increments ∆X = αg S2,0 + βg S3,0 , with α, β ≤ R, are considered: which results in an error of order O(R 2 ). Here, the directors g S2 and g S3 , which are not orthonormal in general, represent the push-forward of the solid directions g S2,0 and g S3,0 : It follows from (34) and (15) that the 2D-3D coupling conditions (31) between the beam surface and the expanded solid position field are exactly fulfilled if the following 1D-3D coupling constraints are satisfied: Coupling constraints of the form (37) enforce that the material fibers g S2 and g S3 of the solid remain orthonormal during deformation, thus enforcing vanishing in-plane strains of the solid at the coupling point x S,r = r. In Section 7, it will be demonstrated that constraints of this type lead to severe locking effects when applied to finite element discretizations that are relevant for the proposed BTS-FULL scheme, i.e. solid mesh sizes that are larger than the beam cross-section dimensions. It will be demonstrated that such locking effects can be avoided if the solid triad field is defined in a manner that only captures the purely rotational contributions to the local solid deformation at x S,r = r without additionally constraining the solid directors in the deformed configuration. As will be demonstrated in the next sections, the rotation tensor defined by the polar decomposition of the deformation gradient is an obvious choice for this purpose, but also alternative solid triad definitions are possible. Table 1 gives an overview of the solid triad variants proposed in the following. All of these solid triad definitions Λ S = [g S1 ,g S2 ,g S3 ] will be a function of the solid deformation gradient F , i.e. Λ S = Λ S (F ). Moreover, all solid triad definitions will be constructed in a manner such that the associated orthonormal base vectorsg S2 andg S3 represent the effective rotation of the non-orthonormal directors g S2 and g S3 in an average sense. Thus, it will be required thatg S2 andg S3 lie within a plane defined by the normal vector in the following denoted as the n-plane. Eventually, in the examples in Section 7, two desirable properties of the solid triad field for the proposed BTS-FULL method are identified: (i) The solid triad should be invariant, i.e. symmetric/unbiased with respect to the reference in-plane beam cross-section basis vectors g B2,0 and g B3,0 . (ii) The resulting BTS-FULL method should not lead to locking effects in the spatially discretized coupled problem.
These properties will be investigated for the following solid triad definitions.

Polar decomposition of the deformation gradient (STR-POL)
Based on polar decomposition, the deformation gradient of the solid problem can be split into a product F = vR S = R S U consisting of a rotation tensor R S ∈ SO 3 and a (spatial or material) positive definite symmetric tensor v or U , respectively, which describes the stretch. An explicit calculation rule for the rotation tensor, e.g. based on v, can be stated as: As mentioned above, it is desirable that the orthonormal base vectorsg S2 andg S3 of the solid triad Λ S lie in a plane with normal vector n according to (38). It can easily be verified that the rotation tensor R S associated with the total deformation gradient F according to (40) will in general not satisfy this requirement. Thus, a modification will be presented in the following to preserve this property.

Construction of STR-POL triad
Since the sought-after solid triad shall be uniquely defined already by the two in-plane directors g S2 and g S3 , a modified version of the deformation gradient will be considered, which consists of the projection of the total deformation gradient F into the n-plane extended by the additional term n ⊗ g S1,0 . This modified deformation gradient ensures that the two relevant in-plane basis vectors are correctly mapped, i.e. g S2 = F n g S2,0 and g S3 = F n g S3,0 , while the third basis vector, which is not relevant for the proposed coupling procedure, is mapped onto the normal vector of the n-plane, i.e. n = F n g S1,0 . This specific definition of a deformation gradient allows for the following multiplicative split: where R n describes the (pure) rotation from the initial solid triad Λ S,0 onto a (still to be defined) orthonormal intermediate triadΛ = [ḡ 1 ,ḡ 2 ,ḡ 3 ], whose base vectorsḡ 2 andḡ 3 lie within the n-plane, and F 2D represents a (quasi-2D) in-plane deformation betweenḡ 2 andḡ 3 and the non-orthonormal base vectors g 2 and g 3 . Now, by applying the polar decomposition only to the in-plane deformation, i.e.
a solid triad can be defined from the initial triad Λ S,0 as: Once an intermediate triadΛ is defined, the required rotation tensors R 2D S and R n can be calculated as follows: The last remaining question is the definition of the triad Λ. It can be shown that the choice of this triad is arbitrary and does not influence the result, since a corresponding in-plane rotation offset would be automatically considered/compensated (in the sense of a superposed rigid body rotation) via the rotational part R 2D S of the in-plane polar decomposition (43). For example, a simple choice is given byḡ 1 = n,ḡ 2 = g S2 / g S2 andḡ 3 = n ×ḡ 2 , which coincides with the solid triad definition later discussed in Section 5.3.1.
S R n is fulfilled for quasi-2D deformation states, e.g. for pure torsion load cases where the beam axis remains straight during the entire deformation (see example in Section 7.5). In this case, the (simpler) polar decomposition of the total deformation gradient F according to (40) can exploited.

Properties of STR-POL triad
In contrast to alternative solid triad definitions that will be investigated in the following sections, the definition according to (44), referred to as STR-POL or by the subscript (·) POL , is not biased by an ad-hoc choice of material directors in the solid that are coupled to the beam. Instead, the rotation tensor R S describes the rotation of material directions coinciding with the principle axes of the deformation (i.e. it maps the principle axes from the reference to the spatial configuration), which has two important implications: First, the choice of material directions that are coupled depend on the current deformation state and will in general vary in time. Second, the principle axes represent an orthonormal triad per definition, and, thus the coupling to the beam triad will not impose any constraints on the local in-plane deformation of the solid. Consequently, this solid triad variant fulfills both requirements (i) and (ii) as stated above.
Eventually, a further appealing property of the STR-POL triad shall be highlighted. Let θ 0 ∈ [−π, π] represent the orientation of arbitrary in-plane directors in the reference configuration defined to coincide for solid and beam according to Their push-forward is given by g S (θ S (θ 0 )) = F n g S,0 (θ 0 ) for the solid and g B (θ B (θ 0 )) = R B g B,0 (θ 0 ) for the beam, where the angles θ S ∈ [−π, π] and θ B ∈ [−π, π] represent the corresponding in-plane orientations in the deformed configuration (see Appendix A.1 for a detailed definition). Since in-plane shear deformation is permissible for the solid but not for the beam, the orientations θ S (θ 0 ) and θ B (θ 0 ) cannot be identical for all θ 0 ∈ [−π, π] and arbitrary deformation states. However, as demonstrated in Appendix A.1, when coupling the beam triad to the STR-POL triad according to (44), the beam directors g B (θ B (θ 0 ) represent the orientation of the solid directors g S (θ S (θ 0 ) in an average sense such that the following L 2 -norm is minimized: In conclusion, STR-POL is an obvious choice for the solid triad with many favorable properties, e.g. it represents the average orientation of material solid directions in a L 2 -optimal manner. However, it requires the calculation of the square root of a tensor, and more importantly, for latter variation and linearization procedures also the first and second derivatives of the tensor square root with respect to the solid degrees of freedom. This results in considerable computational costs, since this operation has to be performed at local Gauss point level. Therefore, alternative solid triad definitions will be proposed in the following that can be computed more efficiently, while still being able to represent global system responses with sufficient accuracy.

Alternative solid triad definitions
All solid triad variants considered in the following rely on the non-orthonormal solid directors g S2 and g S3 according to (35), their normalized counterparts and the corresponding normal vector n according to (38). Based on these definitions, three different variants will be exemplified in the following.

Fixed single solid director (STR-DIR 2/3 )
In the first variant, denoted as STR-DIR 2/3 , the orientation of one single solid director, either g S2 or g S3 , is fixed to the solid triad, cf. Figure 3(b). The choice which solid material direction to couple is arbitrary. Therefore, two variants will be distinguished: Since the variant STR-DIR 2/3 does not fulfill the requirement (i) as stated above, it will only be considered for comparison reasons in the 2D verification examples in Section 7.

Fixed average solid director (STR-AVG)
In order to solve this problem, i.e. to define a solid triad that is symmetric with respect to the base vectors g 2 and g 3 , an alternative variant denoted as STR-AVG is proposed, which relies on the average of the directors g 2 and g 3 , cf.
With this average vector the solid triad can be constructed as: The rotation tensor R (−(π/4)n) in (50) represents a "backrotation" of the constructed reference triad Λ S,AVG,ref by an angle of −π/4 to ensure that the resulting solid triad aligns with the beam triad in the reference configuration according to (32). In Section 7, it will be shown numerically that the variant STR-AVG, similar to the variant STR-POL, fulfills both requirements (i) and (ii) stated above.

Remark 5.2
Theoretically, an additive director averaging procedure such as (49) can result in a singularity if the underlying vectors are anti-parallel, i.e. g S2 = −g S3 . However, since the associated material directors are orthogonal in the reference configuration, i.e. g T S2,0 g S3,0 = 0, and shear angles smaller than π/2 can be assumed, this singularity will not be relevant for practical applications.

Fixed orthogonal solid material directions (STR-ORT)
In the last considered solid triad definition, both material directors g S2 and g S3 are coupled to the solid triad simultaneously. This variant enforces that the directors g S2 and g S3 remain orthogonal to each other, and thus it is denoted as STR-ORT, indicated by a subscript (·) ORT . The STR-ORT variant is realized by applying the rotational coupling constraints (19) twice, once with Λ S,DIR2 according to (47) and once with Λ S,DIR3 according to (48).
Opposed to the other triad definitions in this section, this version additionally imposes a constraint on the solid displacement field by enforcing all shear strain components to vanish at the coupling point. In Section 7, it will be demonstrated that this over-constrained solid triad definition can lead to severe shear locking effects, i.e. requirement (ii) from Section 5.1 is not satisfied. Thus, also this variant will only be considered for comparison reasons in the 2D verification examples in Section 7.

Variation of the solid rotation vector
In the coupling contributions to the weak form (24) and (26) the multiplicative rotation vector variation δθ S (spin vector) of a solid rotation vector ψ S arises. The spin vector is workconjugated with the coupling moments, i.e. it is required to calculate the virtual work of a moment acting on the solid in a variationally consistent manner. In contrast to the beam spin vector δθ B , which represents the multiplicative variation of primal degrees of freedom in the finite element discretization of the geometrically exact Simo-Reissner beam theory and is discretized directly, no such counterpart exists for the solid field. Therefore, it is assumed that the solid spin vector can be stated as a function of a set of generalized solid degrees of freedom q (which will later be identified as nodal position vectors in the context of a finite element discretization) and their variations δq. The additive variation of the solid rotation vector ψ S (q) then reads The multiplicative and additive variations are related via (7), which gives the spin vector associated with the solid triad as a function of the generalized solid degrees of freedom: Remark 5.3 Alternatively, the solid spin vector can be expressed by the variations of the corresponding solid triad basis vectors g Si and their variations δg Si , cf. [36,38]: ∂q δq.
This formulation for the solid spin vector is equivalent to the one in (52), but only contains the solid triad basis vectors and their variations. Therefore, this definition of the solid spin vector is better suited for solid triads constructed via their basis vector. Especially in the implementation of the finite element formulation, it is advantageous to avoid the computation and inversion of the transformation matrix in (52). Nonetheless, in the remainder of this contribution, the solid spin vector as defined in (52) is used to improve readability of the equations.

Spatial discretization
In this work, spatial discretization of the beam, solid and coupling problem will exclusively be based on the finite element method. In the following, a subscript (·) h refers to an interpolated field quantity, superscripts (e) and (f ) indicate that the quantity is defined for a solid element e and a beam element f , respectively. Accordingly, (e, f ) refers to coupling terms between the solid element e and beam element f . The global element count is made up of n el,S solid finite elements and n el,B beam finite elements.

Solid and beam problem
For the solid domain an isoparametric finite element approach is used to interpolate position, displacement and virtual displacement field within each solid element Ω Therein, N (e) ∈ R 3×n (e) dof is the element shape function matrix, which depends on the solid element parameter coordinates ξ S , η S , ζ S ∈ R. Furthermore, X S(e) ∈ R n (e) dof , d S(e) ∈ R n (e) dof and δd S(e) ∈ R n (e) dof are the element reference position vector, element displacement vector and element virtual displacement vector, respectively. Each solid element has n (e) dof degrees of freedom. The beam finite elements used in this work are based on the Simo-Reissner formulation presented in [35,38].    Fig. 4 Degrees of freedom for a single beam element used in this work. All quantities related to the beam centerline position are depicted in blue, all cross-section orientation related quantities are depicted in red.
freedom describing the beam centerline position. The interpolated position of the beam centerline is with the beam position shape function matrix H (f ) ∈ R 3×12 , the beam centerline reference position vector X B(f ) ∈ R 12 and the beam centerline displacement vector d B(f ) ∈ R 12 . Furthermore, ξ B ∈ R is the parameter coordinate along the beam centerline.
A triad interpolation scheme based on three element nodal rotation vectorsψ 3 is utilized [14]. The third node is placed in the middle of the element and carries no translational degrees of freedom, only rotational ones. The three local nodal rotation vectors serve as primal degrees of freedom for the interpolated rotation field along the beam centerline. Each local rotation vector has 3 degrees of freedom, thus resulting in a total of 9 rotational degrees of freedom per beam finite element. The interpolation of the beam cross-section triad along the beam centerline is a nontrivial task and requires an orthonormal interpolation scheme for the interpolated triad Λ (f ) B,h (ξ B ) to guarantee that the interpolated triad field is still a member of the rotational group SO 3 . Furthermore, objectivity of the discrete beam deformation measures has to be preserved by the interpolation, which is a challenging task if rotational degrees of freedom are involved. In this contribution we will refer to the interpolated triad field Λ as an abstract nonlinear function of the beam parameter coordinate and the nodal rotation vectors. The corresponding interpolated field of multiplicative rotation vector increments ∆θ (f ) B,h (ξ B ) has been consistently derived in [14] and reads: Therein,Ĩ (f ) i ∈ R 3×3 are generalized shape function matrices for the multiplicative nodal rotation increments ∆θ are the corresponding element-wise assembled quantities. It should be pointed out thatĨ (f ) i are nonlinear functions of the beam parameter coordinate and the nodal rotation vectors of the beam element, i.e. these rotational shape functions are deformationdependent. To avoid this nonlinearity in the discretized spin vector, i.e. the virtual rotation field δθ (f ) B,h , which would require additional linearization contributions to calculate a consistent tangent, the beam finite elements employed in this work follow a Petrov-Galerkin discretization approach. Therein, standard Lagrange shape functions are used to interpolate the discretized nodal spin vectors: Here, L In what follows, all coupling terms are evaluated on the beam centerline. This requires the projection of points along the beam centerline parameter space into the solid element parameter space, which in turn is achieved by solving the set of nonlinear equations X S(e) h ξ S , η S , ζ S = r (f ) 0,h ξ B , for a given ξ B . To improve readability, the superscripts indicating the beam and solid elements will be omitted from now on. They will however be stated in the virtual work contributions and the integration domains in order to clearly indicate pair-wise values. Additionally, any dependency on element parameter coordinates will not be stated explicitly.
Remark 6.1 While the C 1 -continuous centerline representation of the employed beam elements [35] is not mandatory for the considered beam-to-solid volume coupling problem, it offers significant advantages in problems additionally involving beam-to-beam [39] or beam-to-solid contact interaction, which will be addressed in our future research.

Gauss point-to-segment coupling of cross-section rotations
Evaluating the variation of the total coupling potential (24) based on the discretized solid position field and beam crosssection rotation field as presented in the last section yields the discrete variation of the coupling potential: S,h is the discretized coupling domain between beam element (f ) and solid element (e). The integral in (59) is evaluated numerically via a Gauss-Legendre quadrature, resulting in a Gauss point-to-segment (GPTS) coupling scheme. From a mechanical point of view this can be interpreted a weighted enforcement of the rotational constraints at each integration point along the beam, i.e. a Gauss point-to-segment type coupling: where n GP is the number of Gauss-Legendre points,ξ B i is the beam element parameter coordinate for Gauss-Legendre point i with the corresponding weight w i . Again, in order to improve the readability of the remaining equations in this subsection, the explicit indication of the evaluation at the Gauss-Legendre points will be omitted in the following. The previous equation can now be stated in matrix form as Therein, the abbreviations f B c,GP ∈ R 9 and f S c,GP ∈ R n (e) dof for the generalized Gauss point coupling forces on the rotational beam degrees of freedom and the generalized solid element degrees of freedom, respectively, have been introduced: Furthermore, r B c,GP ∈ R 9 and r S c,GP ∈ R n (e) dof are the beam and solid coupling residual vectors. Employing a Newton-Raphson algorithm to solve the global system of nonlinear equations, a linearization of the residual vectors with respect to the element degrees of freedom is required, which reads: Therein, the transformation matrix T (ψ B,h ) appears, since the linearization is performed with respect to the multiplicative rotation increments ∆θ B,h . Furthermore, the generalized shape function matrixĨ follows from the interpolation of the multiplicative rotation increments, cf. (57). The previously derived matrices and vectors are all defined on beamto-solid element pair level. Since no additional degrees of freedom are introduced, the pair-wise contributions can simply be assembled and added to the global linear system of equations. The Gauss point-to-segment coupling approach is presented here to illustrate how the rotational coupling conditions can be enforced in a point-wise manner. However, in [63] it has been shown that a Gauss point-to-segment coupling approach leads to spurious contact locking for embedded one-dimensional beams in three-dimensional solid volumes. Therefore, this approach will not be investigated further in the remainder of this contribution, but a mortartype coupling is proposed instead.

Mortar-type coupling of cross-section rotations
Employing a mortar-type coupling approach, the rotational Lagrange multiplier field λ θ introduced in Section 4.3.2 is also approximated with a finite element interpolation, cf. [7,46,71]. The rotational Lagrange multiplier field is defined along the beam centerline and accordingly its finite element approximation is defined along the beam finite element and reads as follows: where n (f ) λ is the number of Lagrange multiplier nodes on beam element (f ), Φ (f ) i is the shape function for the local node i and λ is the elementwise assembled Lagrange multiplier shape function matrix for a beam element and λ is the vector with all corresponding discrete rotational Lagrange multiplier values per beam element. As indicated by the dependency on beam parameter coordinate ξ B , the Lagrange multiplier field is defined along the beam centerline. However, there is no requirement that the Lagrange multiplier shape functions are identical to the beam centerline shape functions, or even that the number of beam nodes matches the number of Lagrange multiplier nodes. A more thorough discussion on the choice of Lagrange multiplier shape functions is given at the end of this section.

When inserting the finite element interpolations, the discretized variation of the coupling constraints (26) reads
δW (e,f ) Therein, the abbreviations g c,λ θ ∈ R n (f ) dof,λ and r c,λ θ ∈ R n (f ) Therein, the abbreviations f B c,λ θ ∈ R 9 and f S c,λ θ ∈ R n (e) dof represent the integrand of the beam and solid element coupling forces, i.e.
Furthermore, r B c,λ θ ∈ R 9 and r S c,λ θ ∈ R n (e) dof are the beam and solid coupling residual vectors, respectively. Again, a linearization of the residual contributions with respect to the discrete beam-to-solid pair degrees of freedom is required for the Newton-Raphson algorithm. The linearization is: Therein, the abbreviations q (·) for the stiffness matrices of the pair-wise coupling terms have been introduced, i.e.
As in the GPTS case, the previously derived vectors and matrices are all defined on beam-to-solid element pair level. However, in this case additional unknowns have been introduced, i.e. the rotational Lagrange multipliers λ θ on pair level. In practice, all derivatives explicitly stated in (66) and (68) are evaluated using forward automatic differentiation (FAD), cf. [28], using the Sacado software package [56], which is part of the Trilinos project [67].
At this point it should be pointed out that all coupling integrals are evaluated numerically using so-called segmentbased integration, cf. [18,63]. Therein, the beam finite element parameter space is divided into subsegments at points where the beam crosses a solid finite element face. Each subsegment is subsequently integrated using a Gauss-Legendre quadrature with a fixed number of integration points. This leads to a highly accurate numerical integration procedure and allows for the resulting finite element coupling method to pass classical patch tests in surface-to-surface problems as well as constant stress transfer tests in beam-to-solid problems, cf. [18,63].
The choice of proper Lagrange multiplier basis functions is important for the mathematical properties of the resulting finite element discretization. The Lagrange multiplier shape functions must fulfill an inf-sup condition to guarantee stability of the mixed finite element method [11]. This is a wellstudied topic in the context of classical surface-to-surface mesh tying or contact. However, as pointed out in [63], beam-to-solid coupling problems diverge from the standard surface-to-surface case in some aspects. First, the discretization along the beam centerline with Hermite polynomials is unusual compared to standard (i.e. Lagrange polynomialbased) finite element discretizations. Also, the 1D-3D coupling can be classified as a mixed-dimensional embedded mesh problem, since there is no explicit curve in the solid domain to match the beam centerline, which can lead to stability issues [57]. Additionally, in this contribution we deal with rotational coupling, which also differs from the standard surface-to-surface case. A deep mathematical analysis of these properties is beyond the scope of the present contribution. However, we will build upon the extensive studies and findings from [63], where it has been shown that a linear interpolation of the Lagrange multipliers combined with a penalty regularization leads to a stable finite element formulation of the coupled problem. Instabilities might only occur if the beam finite elements become shorter than the solid finite elements. However, this is not a mesh size relation that is within the envisioned applications for the BTS-FULL method. Alternative approaches to avoid the mentioned instabilities are available in the literature, such as Nitsche's method [15,23,58], or discontinuous Galerkin formulations [23,57]. Another appealing approach is the so-called vital vertex method introduced in [6], where discrete Lagrange multipliers are inserted at intersection points between the coupled meshes, i.e. the intersections between a beam finite element centerline and the solid elements in the presented beam-to-solid case.

Beam-to-solid volume coupling (BTS-TRANS)
In Sections 6.2 and 6.3, GPTS and mortar-type coupling discretizations for the rotational coupling between the beam cross-section and the solid volume (19) have been presented. The translational coupling of the beam centerline (18), however, is entirely based on the mortar-type beam-to-solid volume coupling (BTS-TRANS) method previously introduced in [63]. The resulting linearized system of equations for the centerline translation coupling with saddle point structure reads: Therein, K S is the solid tangent stiffness matrix, D S and ∆D S are the global solid displacement vector and its increment, respectively, and R S is the residual of the solid degrees of freedom. In (70), the beam terms are split into centerline and rotation contributions indicated by (·) r and (·) θ , respectively. At this point it should be noted that due to the employed Petrov-Galerkin method the beam stiffness matrices are non-symmetric, as will be the case for the rotational coupling contributions to the global stiffness matrix. Furthermore, D B r and ∆D B r are the global beam centerline displacement vector and its increment, and ∆D B θ represents the global vector collecting the multiplicative rotation increments associated with the nodal triads of the beam finite elements. The update of the rotation state has to be preformed according to [38]. The BTS-TRANS coupling is represented by the discrete mortar matrices D and M and the centerline Lagrange multipliers Λ r . The structure of the global stiffness matrix in (70) illustrates that only the beam centerline degrees of freedom (and not the rotational degrees of freedom) are coupled to the solid degrees of freedom in the previously proposed BTS-TRANS coupling scheme.

Combined mortar-type coupling of translations and rotations (BTS-FULL)
The global system of equations for the mortar-type BTS-FULL method is the combination of the mortar-type coupling for the centerline positions BTS-TRANS (70) and the mortartype coupling of the beam cross-section rotations (68). The resulting global system of equations becomes: where Q (·) are the globally assembled pair-wise stiffness matrices q (·) , and R (·) c,λ θ are the globally assembled residual vectors of the local rotational coupling contributions r (·) c,λ θ . Additionally, Λ θ are the globally assembled rotational Lagrange multipliers λ θ . Therefore, the size of the global system of equations of the uncoupled system is extended by the total number of translational and rotational Lagrange multipliers.

Remark 6.2
The structure of the global system of equations for BTS-FULL (71) illustrates the direct coupling of the rotational degrees of freedom of the beam with the solid degrees of freedom, i.e. Q Bλ θ and Q λ θ B . Disregarding all other advantages of our BTS-FULL method, this motivates its from a pure numerical point of view, as possible rigid body rotations of straight embedded fibers around their centerline are constrained, which is not the case for the BTS-TRANS method (70).

Remark 6.3
In the BTS-TRANS method, the mortar-type coupling matrices D and M only depend on the reference configuration, i.e. they only have to be calculated once and can be stored for the entire simulation. In the BTS-FULL method, the (rotational) coupling terms Q (·) depend on the current configuration, i.e. the coupling terms have to be re-evaluated in each Newton-Raphson step. However, this should not be viewed as a drawback of BTS-FULL scheme, rather as a simplification of the BTS-TRANS variant, which results from neglecting the rotational coupling terms.

Penalty regularization
In the present mortar-type coupling case (BTS-FULL) the constraint equations are enforced with the Lagrange multiplier method, thus resulting in a mixed formulation. However, a direct solution of the global system (71) might introduce certain drawbacks, such as an increased system size compared to the uncoupled system and a generalized saddle point structure. In [63] the constraint equations have therefore been enforced using a well-known penalty regularization, which means that a relaxation of the translational coupling constraints −MD S + DD B r = R c,r = 0 in the form of Λ r = r V −1 r R c,r is introduced. Therein, r ∈ R + is a scalar penalty parameter and V r is a scaling matrix to account for non-uniform weighting of the constraint equations [63,72]. The numerical examples in [63] show that for reasonably chosen penalty parameters the resulting violation of the constraint equations due to their relaxation does not have any impact on the accuracy of the BTS-TRANS method. Therefore, the constraint enforcement of the new rotational coupling equations (65) is also carried out with the penalty method. The constraint relaxation is achieved through again with a scalar penalty parameter θ ∈ R + and a global scaling matrix for the rotational Lagrange multipliers V λ θ . The global scaling matrix is assembled from the nodal scaling matrices κ With the introduction of the constraint relaxation (72), the Lagrange multipliers Λ θ are no longer independent degrees of freedom of the system, but a function of the beam rotations and solid displacements. Therefore, they can be eliminated from the global system of equations (71), which results in the condensed linear system of equations Therein, the following abbreviations have been introduced for improved readability:

Examples
The following numerical examples are set up using the beam finite element pre-processor MeshPy [64] and are simulated with our in-house parallel multi-physics research code BACI [4].

Single element moment test
The first problem setup is depicted in Figure 5. A straight beam is embedded inside a solid cube (E = 1 N/m 2 , ν = 0) and the beam is loaded with a distributed torsion moment in e 1 direction, which is constant along the beam centerline. This example is used to investigate how a moment on a beam is transferred to solid nodal forces. The cube is modeled with a single eight-noded hexahedral element and all solid degrees of freedom are fixed. A single Simo-Reissner beam finite element is used to discretize the beam. No Dirichlet boundary conditions are applied on the beam and the coupling between the beam and the solid is realized with our novel mortar-type BTS-FULL method. Thus, the only interaction between the beam and the solid is the transfer of the external moment. The resulting nodal reaction forces for the different solid triad definitions introduced in Section 5 are depicted in Figure 6. Therein, the results for the STR-POL, STR-AVG and STR-ORT variants, cf. Figures 6(a), 6(d) and 6(e), match up to machine precision. In general, however, the solid coupling reaction forces may differ for different definitions of the solid triad, as visible for the variants STR-DIR 2/3 in Figures 6(b) and 6(c). This observation can be explained by the fact that the representation of a moment via nodal forces is non-unique, i.e. there is an infinite number of possible force pair combinations to achieve this. However, from a mechanical point of view, the force pairs resulting from the STR-POL, STR-AVG and STR-ORT variants seem more natural than the ones for the STR-DIR 2/3 variants. Moreover, the former three variants result in the (unique) force pair solution if the moment is applied as a constant shear stress on the beam surface, cf. Section 2. Additionally, it can be observed for the STR-DIR 2/3 variants that the choice which local solid direction is coupled to the solid triad drastically affects the result for the nodal forces.

Shear test
The next elementary test case is illustrated in Figure 7. The problem geometry is the same as in the previous example. cube, as depicted in Figure 7. No boundary conditions are applied to the beam. This problem illustrates how the specific solid triads affect the shear stiffness of the solid element and will be studied in two steps.
In a first step, we want to investigate the impact of the local stiffening effect the beam cross-section has on the surrounding solid material. To do so, a reference solution is created by applying a full 2D-3D beam-to-solid coupling  Figure 8 illustrates the shear stress in the solid, with and without the embedded beam, for an exemplary beam radius r = 0.1 m. As expected, the solution is uniform in the entire solid volume for the pure solid variant. In the 2D-3D beam-to-solid coupling variant, the embedded beam affects the solid stress and displacement fields. The overall displacement of the solid is smaller than for the variant without a beam, thus demonstrating the stiffening effects of the beam cross-section. In agreement with our fundamental modeling assumption of overlapping beam and solid domains (see Section 2), the solid shear stress inside the beam domain is zero. Outside of the beam domain, the solid shear stress field shows slight fluctuations due to the local constraints enforcing the 2D-3D coupling at the beam surface. However, close to the boundaries of the cube, these fluctuations become negligible and the shear stress field is quite homogeneous and therefore very similar to the pure solid shear problem.
In a second step, this problem is simulated with one single solid finite element to investigate potential shear locking effects. The coupling between beam and solid is now realized with our novel BTS-FULL method and a rotational penalty parameter of θ = 100 Nm/m. In Figure 9, the deformed solid element and the resulting coupling reaction forces on the solid nodes are depicted for the different solid triad definitions and again for the problem without embedded beam. Due to the orthogonality constraints in the STR-ORT variant, no shear mode remains in the solid finite element, i.e. it is rigid with respect to shear deformations (in fact, small deformations can be observed due to the penalty regularization). In this example, all other solid triad definitions result in a solid displacement field matching the variant without embedded beam up to machine precision. Table 2  The results show that the presented solid triads lead to either no shear stiffening effects in the solid (STR-POL, STR-DIR 2,3 and STR-AVG) or to severe locking resulting in a complete constraining of all shear modes (STR-ORT). To assess which variant resembles best the resolved 2D-3D coupling scheme, the relative L 2 -displacement error is compared. In the results presented in the following, the reference solution is the solution obtained with a fine solid mesh and a 2D-3D coupling. Figure 10 illustrates e L2,rel for different beam diameter to solid cube length ratios d/h. The relative error for the STR-ORT variant is almost constant 1 for all beam diameters ratios, i.e. even for beam cross-section sizes similar to the cube dimensions a full constraining of all shear modes does not accurately describe the physical coupling. For all other variants the behavior of the relative error is the same, since none of them constrain the shear deformation mode in the solid, i.e. the beam cross-sections rotate with the solid without constraining it. For small ratios of beam radius to solid cube length the error is close to zero. For larger ratios of beam radius to solid cube length, the error increases as there is a real physical stiffening effect due to the embedded beam cross-section in the 2D-3D problem that is not captured by the 1D-3D coupling schemes. However, in the entire range of practically relevant solid mesh sizes (relative to the beam cross-section size) as illustrated in Figure 10, the solid triad variants that do not constrain the in-plane deformation of the solid result in a better approximation of the physical system behavior as compared to the STR-ORT triad. In this example, the results obtained with the BTS-FULL method and different solid triads will be compared with a spatially converged reference solution, where the coupling between the beam surfaces and solid volume is discretized in a surface-tovolume (2D-3D) manner, i.e. the beam surface instead of the centerline is fixed to the solid, cf. Appendices A.2 and A.3. The resulting shear stresses are visualized in Figure 12. In the full 2D-3D model, there are stress concentrations at the interface between the beam surfaces and the solid. It is important to point out that the BTS-FULL method (1D-3D), is not able to capture these stress concentrations, regardless of the employed solid triad. However, this has not been the intention of the BTS-FULL method in the first place, but instead we want to make sure that the far field stress in the solid is represented accurately. Figure 12 illustrates the shear stress results obtained with the STR-POL, STR-DIR 2/3 and STR-AVG solid triads. In the reference solution the in-plane shear stress is positive at the top and bottom of the cube and negative in the middle. The results obtained with the STR-POL, STR-DIR 2/3 and STR-AVG solid triads are similar to the ones obtained with 2D-3D coupling. However, the results with the STR-ORT solid triads clearly exhibit drastic shear locking effects due to the (over-) constraining of orthogonal solid directions. Table 3 provides the displacement at the top right corner of the cube for the 2D-3D reference solution and different types of solid triad fields, as well as the relative error. The error for the STR-ORT solid triad is six times larger than for all other solid triads. This again illustrates the unwanted locking effects introduced by the STR-ORT solid triads variant.
At this point a short recap of the first three examples for each of the investigated solid triad constructions is given to summarize their applicability in the context of our BTS- Therefore, these variants will not be employed in the following. However, for comparison purposes they will be included in the spatial convergence example, cf. Section 7.5. STR-AVG All basic consistency tests are fulfilled by the averaged solid triad and the results are very close to the ones obtained via the STR-POL variant, while being less expensive from a computational point of view. This variant is used in the remaining examples of this contribution. STR-ORT This variant leads to considerable shear locking in the range of coarse solid mesh resolutions, which is exactly the range of interest for the proposed 1D-3D coupling schemes. Therefore, this variant will not be used in the remainder of this contribution.

Transfer of constant torque
This example serves as a consistency test for the BTS-FULL method and its ability to transfer a constant torque. It is an extension of the constant stress transfer problem for the BTS-TRANS method previously presented in [63]. The example is inspired by classical patch tests, which are well-established tools to investigate the consistency of finite element formulations [65]. The constant torque test is depicted in Figure 13. It consists of a solid block Ω S with two embedded beams Ω B1 and Ω B2 . The two beams occupy the same spatial position. The solid is fixed at the lower surface and no external loads are applied. One beam is loaded with a torsion load m, and the other beam with a torsion load −m, both acting along their axial direction. The magnitude of the torsion load is 10 Nm/m. Based on the space-continuous problem description, the opposing loads on the two beams cancel out each other, and in sum the two beams transfer no loads to the solid. This gives the trivial solution u S = 0 for the displacement field in the solid, cf. [63], and a constant solution for the beam rotations along their axis. In this test it shall be verified that this solution can also be represented in the spatially discretized setting using an arbitrarily coarse discretization. Both beams are coupled to the solid via the BTS-FULL method. There is no direct interaction between the two beams, but all interactions are transferred through the solid domain. The geometry and material parameters are taken from [63]. The dimensions of the solid block are 1 m × 1 m × 2 m and a Saint Venant-Kirchhoff material model (E = 10 N/m 2 , ν = 0.3) is employed. The block is discretized with 4 × 4 × 7 eight-noded, first-order hexahedral elements. The circular cross-sections of the two beams have a radius of 0.05 m, and the beam material parameters are E = 100 N/m 2 and ν = 0. The beams B1 and B2 are discretized with 5 and 7 Simo-Reissner beam finite elements, respectively. This results in a non-matching discretization between the two beams as well as between the beams and the solid. Coupling between the beams and the solid is realized with a linear interpolation of both the translational and rotational Lagrange multipliers. The STR-AVG solid triads are employed in this example, cf. Section 5.3.2. The penalty parameters are r = 100 N/m 2 and θ = 100 Nm/m. Figure 14 illustrates the results of this test. The stress in the solid and the curvature in the beam are indeed zero up to machine precision, thus matching the expected analytical solution.
This example illustrates the ability of the BTS-FULL method to exactly represent a constant torsion state along the beam and the consistency of the coupling terms despite the fact that arbitrary non-matching meshes are involved.

Spatial convergence
This numerical example investigates the spatial convergence properties of the BTS-FULL method under uniform mesh refinement. The problem is depicted in Figure 15  thus resulting in a total torque of 1.65885 · 10 −2 Nm. This example can be interpreted as an adapted version of the spatial convergence problem in [63] to verify the scenario of rotational coupling. A similar problem is also investigated in [27]. The spatial convergence behavior of the BTS-FULL method will be analyzed with respect to a spatially converged reference solution obtained with a 2D-3D coupling discretization, as described in Appendix A.2. To compare the results, the L 2 displacement error in the solid is calculated via Here, V 0 = 1 m 3 is the solid volume in the reference configuration. It should be pointed out that the 2D-3D coupling problem does not have the same analytical solution as the BTS-FULL problem, because the 1D-3D coupling results in a singularity in the analytical solution, cf. Section 2.1. Therefore, spatial convergence of the BTS-FULL method towards the reference solution is not expected all the way towards the asymptotic limit of arbitrarily small solid element sizes, but only in the practically relevant regime of solid mesh sizes that are larger than the beam cross-section radius. In this regime, the singularity, i.e. the difference between the 1D-3D and 2D-3D models can not be fully resolved by the finite element solution space. This fact can be exploited to obtain reasonably accurate results with our BTS-FULL (i.e. 1D-3D) method for the envisioned applications and practically relevant mesh resolutions. The solution to the presented problem has a point symmetry around the e 1 axis. Therefore, the STR-POL and STR-AVG solid triad variants coincide and give the same numerical results up to machine precision. Similarly, the results obtained with the STR-DIR 1/2 variants match up to machine precision. Figure 16 shows the convergence plot for different types of solid triads as well as for the 2D-3D coupling approach. For coarse discretizations an excellent convergence behavior can be observed for the BTS-FULL variants, slightly below the convergence rate of the reference 2D-3D method, but with a significantly reduced computational cost. All BTS-FULL convergence plots exhibit a kink at a certain solid mesh resolution: the STR-DIR 1/2 variants at around h solid = 0.07 m, and the STR-POL and STR-AVG variants at around h solid = 0.06 m. Figuratively speaking, the difference between the 1D-3D and 2D-3D coupling model becomes dominant at the kink position, since the solid element size to beam cross-section diameter ratio becomes smaller. Nevertheless, the kink only occurs when the solid element size is already smaller than the cross-section radius, which is far away from the envisioned geometric relations for the BTS-FULL method anyways. The results confirm that for solid element sizes larger than the beam cross-section diameter, i.e. our desired and practically relevant discretization case, the results obtained with the BTS-FULL method (1D-3D) exhibit excellent spatial convergence properties and thus give a very good approximation of the 2D-3D coupling problem.

Plane cantilever bending
In this example we consider a cantilever structure modeled as a solid continuum subject to a moment load. The problem is illustrated in Figure 17 This moment is chosen such that, if the cantilever were modeled using 1D beam theory, it should bend exactly to a quarter circle, due to a pure bending deformation in the region between the Dirichlet boundary and the applied moment. Directly imposing a conservative moment load on a solid, i.e. a Boltzmann continuum, which exhibits no rotational degrees of freedom is a non-trivial task. Standard approaches would require to model the moment as a (deformation-dependent) load/traction field distributed across an arbitrarily chosen sub-volume of the solid. In this example we impose the external moment on the solid structure by defining a solid triad (STR-AVG) at the application point of the moment. The nodal external forces effectively acting on the solid are obtained by projecting the moment to the solid finite element space via the discrete version of (52). The cantilever is discretized with 20 × 4 plane fournoded, first-order quadrilateral elements. In Figure 17(b) the deformed cantilever is illustrated. The global displacement behavior is as expected, i.e. the cantilever bends to a quarter circle. Of course, the local strain state close to the point where the external moment is applied is not meaningful in a continuum mechanics sense, since we impose a singular moment at that point. However, according to Saint Venant's principle, a linear stress distribution across the beam height, as expected for the pure bending of a slender beam-like structure, can be observed at a sufficient distance from the point where the moment is induced. This example illustrates that the presented rotational coupling approach is not limited to the coupling of beam cross-section orientations, but can also be used as a stand-alone feature to impose moments onto a solid domain in a variationally consistently manner. It should be pointed out that this example has only been carried out in 2D for reasons of simplicity, while the illustrated capability is available in 3D problems, too.

Plate with embedded beam
In this example a beam is only partially embedded inside a solid plate and loaded with a tip force. Two different geometry variants of the embedded beam are considered, cf. Figure 18. The coupling of beam and solid is realized with our novel BTS-FULL method and compared to the BTS-TRANS method from [63], i.e. the one without rotational coupling.
First-order Lagrange polynomials are employed to discretize the translational and rotational Lagrange multipliers. The penalty parameters are r = 100 N/m 2 and θ = 100 Nm/m. The solid plate is modeled with 1 × 10 × 10 eight-noded solid-shell elements [10,69], while the entire beam is discretized with six Simo-Reissner beam finite elements. The resulting global finite element model has 807 degrees of freedom. A full 3D model, also resolving the beam with threedimensional solid finite elements and consisting of 90,190 second-order tetrahedra (tet10) elements, serves as a comparison. The discretization of the full 3D model has been chosen such that mesh convergence is guaranteed. Consequently, the full 3D model consists of 270,570 degrees of freedom.
The results for variant A are shown in Figure 19. It can be seen that the the full 3D model and the new BTS-FULL method exhibit the same overall behavior, while the beam experiences much larger deformations and the solid smaller ones in the BTS-TRANS model without rotational coupling. This is due to the fact that in the full 3D problem a considerable portion of the external load is transferred from the beam to the solid via shear stresses on the beam surface, which are represented by moments in the reduced-dimensional model. Only the new BTS-FULL method is able to capture these coupling moments. Figure 20 shows the results for variant B. In this case a solution for the purely translational BTS-TRANS method (i.e. only centerline position coupling) does not even exist within a quasi-static framework, since the beam has an unconstrained rigid body rotation mode around its axis of the embedded part. Again, the displacement results of the full 3D problem and the BTS-FULL model are very close to each other. A more detailed comparison of the different variants in given in Figure 21. Therein, the displacements along the curve indicated in Figure 18 are visualized. Now it also becomes clear quantitatively that the displacement results obtained with the BTS-FULL method are very close to the ones obtained with the full 3D problem. Considering that the former reduces the number of degrees of freedom by a factor of about 330 as compared with the latter, this is a remarkable result and showcases the efficiency of the new BTS-FULL method for challenging applications.

Twisted plate
In this final example we consider a plate, with complex, spatially distributed fiber reinforcements in 3D, cf. to the e 2 axis to make the example more challenging and represent general 3D fiber-solid element intersection scenarios. The embedded fibers have a cross-section radius r = 0.01 m and the material parameters are E = 400 N/m 2 and ν = 0. The coupling of fibers and solid is realized with the BTS-FULL method (STR-AVG solid triad, r = 100 N/m 2 and θ = 100 Nm/m). First-order Lagrange polynomials are employed to discretize the translational and rotational Lagrange multipliers. The solid plate is modeled with 10 × 35 × 2 eight-noded solid-shell elements, while each fiber is discretized with four Simo-Reissner beam finite elements, thus resulting in a total of 92 beam finite elements. The displacement controlled twisting deformation of the plate is applied within 100 quasi-static load steps. At this point it should be mentioned that this example could not be solved with the BTS-TRANS method, since the rigid body rotation modes of the straight fibers lead to a non-converging Newton-Raphson algorithm in the very first load step. This underlines the advantages of the mechanically consistent coupling provided by the BTS-FULL method. Figure 23 illustrates the deformed structure at different load steps. Until load step 75, the reinforced plate exhibits a more or less homogeneous twist along the e 2 axis. From load step 75 to load step 100, the reinforced plate folds around the e 2 axis. To assess the non-linear behavior of this structure and evaluate the global impact of the fiber-reinforcements, the fiber-reinforced plate is compared to a simple plate (same material) without any fibers. Figure 24 depicts the reaction moment M 2 around the e 2 axis at the fully clamped surface of the plate with and without fiber-reinforcements. Until load step 70, the structures behave similarly. However, as expected the fiber-reinforcements lead to an increased reaction moment for the same twist angle φ, i.e. to a stiffer structural response. Both structures exhibit a limit point with an unstable post-critical solution, i.e. the structures would collapse if the twist is applied in a load-controlled manner. The fiber-reinforcements affect the critical point of the structure such that the instability occurs at a smaller twist angle and the critical moment is increased. This illustrates the complex influences that fiber-reinforcements may have on the global non-linear behavior of a structure. Figure 25 illustrates the final configuration and shows a close-up view of the deformed embedded fibers. The maximum normal stresses in the fibers resulting from axial and bending deformations can be estimated for this example as ≈ 15 N/m 2 and ≈ 26 N/m 2 (not visualized in the figure), respectively. In the solid the maximum principal Cauchy stress is 0.578 N/m 2 (not visualized in the figure). As expected, the stresses in the stiff fibers are much larger than in the relatively soft solid matrix. To further investigate the influences of the different deformation modes of the fibers, Figure 26 depicts the tension, shear, torsion and bending contributions to the total internal elastic energy of the fibers over the course of the   simulation. In the first few load steps, the main contributors to the internal elastic energy of the system are bending and torsion deformations, cf. right part of Figure 26. This can be attributed to the fact that in the beginning of the simulation the deformations of the plate mainly take place in in e 3 direction, which predominantly causes bending and torsion deformations of the fibers. As the plate is twisted further, geometrically non-linear effects materialize especially on the outer edges of the plate. The edges form helix like curves. Due to the constrained displacements in e 3 direction at the clamped surfaces, the outside edges of the plate are stretched in e 3 direction. This causes axial tension in the fiber semicircles at the outside. Starting at approximately load step 25, the main contribution to the internal elastic fiber energy comes from axial deformations. In the post-buckling state, the bending deformation of the plate, and therefore also of the fibers increases. This causes an increase in the internal elastic bending and torsion energy of the fibers. Moreover, shear deformations only have a minor contribution to the total internal energy of the fibers, which is expected due to the slenderness of the embedded fibers, thus motivating a future use of the BTS-FULL method in combination with shear stiff Kirchhoff-Love beam theories [36,37]. The considerable contributions of bending and torsional energy to the internal elastic energy of the fibers demonstrates the importance of consistently representing these modes and coupling them to the background solid material as done by our proposed BTS-FULL scheme. For this example, this would not be the case if simplified models for the fibers (e.g. modeled as strings without bending stiffness) or for the fiber-solid coupling (e.g. BTS-TRANS) were applied.
This example also showcases the maturity of the implemented BTS-FULL method from an algorithmic point of  view. The chosen solid mesh, in combination with the tilted fiber semicircles results in complex 3D intersections between the beam finite elements and the solid finite elements, thus illustrating the robustness of the employed numerical integration algorithm. As a final example, a more complex model of a fiber-reinforced plate is considered. Therein, the dimensions of the plate are repeated 5 times in e 1 and e 2 direction and 3 times in e 3 direction. The pattern and size of the fiber-reinforcements is similar to the one illustrated in Figure 22, however, in this case there are 3 layers of fiberreinforcements over the thickness of the plate. This results in a total of approximately 53,000 solid finite elements and 1,800 fibers with 4 beam finite elements each, i.e. the problem size is scaled by a factor of approximately 75 compared to the previously considered plate. The deformed configuration of the plate is visualized in Figure 27 . This further illustrates the robustness and scalability of the presented BTS-FULL method for large-scale problems.

Conclusion
In this work we have proposed a 1D-3D coupling method to consistently embed 1D Cosserat beams into 3D Boltzmann continua (solids). Six constraint equations act on each point along the Cosserat beam centerline, namely three transla- tional constraints and three rotational constraints, thus resulting in a full, mechanically consistent coupling between the 1D beams and the 3D continuum. Deriving the full 1D-3D coupling on the beam centerline from a 2D-3D coupling on the beam surface via a Taylor series expansion of the solid displacement field would require to fully couple the deformed solid directors with the undeformable beam crosssection triad. It is demonstrated that such an approach, which suppresses all in-plane deformation modes of the solid at the coupling point, might result in severe locking effects in the practically relevant regime of relatively coarse solid mesh sizes. Therefore, a suitable triad field has to be defined in the 3D Boltzmann continuum that only represents solid material directions in an average sense without constraining them. It has been shown that the rotational part of the polar decomposition of the (in-plane projection of the) solid deformation gradient is a natural choice, since it represents the average orientation of material directions of the 3D continuum in a L 2 -optimal manner. Additionally, several other solid triad definitions have been presented, which allow for a more efficient numerical evaluation. Future work will focus on the extension of the proposed beam-to-solid volume coupling approach to beam-to-solid surface coupling as well as to beam-to-solid surface contact. Especially in the latter case, the proposed rotational coupling constraints are essential to capture effects such as frictional contact between beam and solid. Another topic of interest for future research is the combination of 1D-3D and 2D-3D coupling within a unified beam-to-solid coupling approach. This would allow to use 2D-3D coupling along with a refined solid mesh only in domains where high resolution of solid stress fields is of interest, and using the proposed, highly efficient 1D-3D coupling approach in the remaining problem domain. Moreover, also a combination of the developed schemes with concepts allowing for a consistent coupling of the beam ends with the solid domain, cf. [54], are considered as promising future research direction.

A Appendix
A.1 Proof of L 2 -optimality of STR-POL triad In the following, a proof shall be given for (45). First, an angle θ 0 ∈ [−π, π] is defined that represents the orientation of arbitrary in-plane directors in the reference configuration defined to coincide for solid and beam according to g S,0 (θ 0 ) = g B,0 (θ 0 ) = cos (θ 0 ) g B2,0 + sin (θ 0 ) g B3,0 . The push-forward to the spatial configuration is given by g S = F n g S,0 (θ 0 ) for the solid and g B = R B g B,0 (θ 0 ) for the beam. As stated in Section 5.1, it is a desirable property of the (to be defined) solid triad and therefore also of the beam triad that the base vectors g B2 and g B3 lie in the n-plane, i.e. the plane spanned by the solid base vectors g S2 and g S3 . Thus, in analogy to (42) as stated for the solid, it is assumed that the total rotation of the beam cross-section R B is split in a multiplicative manner into two successive rotations where R n describes the 3D rotation from Λ B,0 toΛ and R 2D B the quasi-2D rotation fromΛ to Λ B . Thus, after push-forward to the intermediate configuration defined by R n the corresponding directors of solid and beam still coincide: g S (θ 0 ) =ḡ B (θ 0 ) = R n g B,0 (θ 0 ) = cos (θ 0 )ḡ 2 + sin (θ 0 )ḡ 3 .
In the following, the material and spatial principle axes associated with the polar decomposition (43) of the in-plane deformation gradient F 2D are denoted as G P 2 and G P 3 as well as g P 2 = R 2D S G P 2 and g P 3 = R 2D S G P 3 . Since the principle axes G P i and the orthonormal base vectorsḡ i are related by quasi-2D rotations with respect to the normal vector n, the directors in (80) can alternatively be stated as g S (θ 0 ) =ḡ B (θ 0 ) = cos (θ 0 ) G P 2 + sin (θ 0 ) G P 3 , whereθ 0 = θ 0 − θ diff is defined via the constant offset value θ diff describing the rotation fromḡ i to G P i . The final beam cross-section triad follows from the second (quasi-2D) rotation R 2D B = R(θ 2D B n) fromḡ i to g Bi described by the scalar rotation angle θ 2D B . In a similar fashion the (quasi-2D) rotation R 2D S = R(θ 2D S n) from G P i to g P i is described by the scalar rotation angle θ 2D S . Due to the 2D-nature of these rotations, the beam directors g B in the spatial configuration can be derived from (81) according to: Let the principle stretch ratios associated with the in-plane deformation gradient F 2D be denoted as λ 2 and λ 3 . Then, the solid directors g S = F n g S,0 (θ 0 ) in the spatial configuration can be derived according to: S [cos (θ 0 ) G P 2 + sin (θ 0 ) G P 3 ] = λ 2 cos (θ 0 ) g P 2 + λ 3 sin (θ 0 ) g P 3 . (83) Here, the relation g P i = R 2D S G P i and the diagonal structure v 2D = λ 2 g P 2 ⊗ g P 2 + λ 3 g P 3 ⊗ g P 3 of the spatial stretch tensor has been exploited. From (82) and (83), the orientation angles of the spatial beam and solid directors g B (θ 0 ) and g S (θ 0 ) relative to the spatial principle axis g P 2 can be identified according to: θ S (θ 0 ) = arctan λ 3 sin (θ 0 ) λ 2 cos (θ 0 ) .
By exploiting the property θ S (−θ 0 ) = −θ S (θ 0 ) of (85), it can easily be shown that (87) results in the requirement: This means that the beam directors g Bi have to coincide with the principle axes g P i and, thus, the total beam triad has to satisfy Λ B = R 2D S R n Λ B,0 , which is identical to the solid triad definition STR-POL according to (44) with the initial condition Λ B,0 = Λ S,0 . By checking the second derivative, it can easily be confirmed that this solid triad choice indeed results in a minimum of the L 2 -norm in (86).

A.2 Full 2D-3D coupling
In the example section we compare the 1D-3D (i.e. BTS-FULL) method to reference solutions obtained with a 2D-3D coupling approach. For the sake of completeness, we state the kinematic coupling constraints for the employed 2D-3D coupling approach. The coupling constraints read Therein, Γ 2D-3D is the 2D-3D coupling surface, i.e. the part of the beam surface that lies within the solid volume. Furthermore, r CS ∈ R 3 is the cross-section position vector, i.e. the vector that points from the cross-section centroid to the cross-section perimeter. The cross-section position vector can be expressed by the current beam triad basis vectors g B2 and g B3 , or via the cross-section rotation tensor Λ B and the Cartesian basis vectors e 2 and e 3 , i.e.
Therein, α ∈ R and β ∈ R are the beam cross-section coordinates, i.e. they parametrize the beam cross-section. In the following, two methods to enforce the 2D-3D coupling conditions (89) are presented, once with a Lagrange multiplier method and once with a quadratic penalty potential.

A.2.1 Penalty potential
The quadratic penalty potential reads Here, 2D-3D ∈ R is a scalar penalty parameter. Variation of the penalty potential gives the following contributions to the weak form: Therein, the coupling force f 2D-3D , acting on the beam centerline and solid, can be identified. Furthermore, m 2D-3D is the coupling moment acting on the beam cross-section. This demonstrates the projection of purely positional coupling constraints (on the surface of the beam) onto the beam centerline, and illustrates the arising rotational coupling terms in a 2D-3D coupling approach.

A.2.2 Lagrange multiplier potential
The 2D-3D coupling conditions (89) can also be enforced with a Lagrange multiplier method. A Lagrange multiplier vector field λ 2D-3D ∈ R 3 is therefore defined on the coupling surface Γ 2D-3D . The total Lagrange multiplier potential for the 2D-3D coupling reads The variation of the total Lagrange multiplier potential gives the following contributions to the weak form: Again, this showcases the projection onto the beam centerline, in this case of the Lagrange multiplier field λ 2D-3D , i.e. the coupling surface tractions on the beam surface.
A.3 Gauss point-to-segment approach for full 2D-3D coupling Evaluating the variation of the total coupling potential (92) on the basis of the discretized solid position field and beam cross-section rotation field yields the discrete variation of the 2D-3D coupling potential: Therein, Γ 2D-3D,h is the discrete beam surface. It is important to point out that the beam surface is not directly discretized. It is an analytical surface defined by the discretized beam centerline, the beam crosssection orientations and the beam cross-section geometry. Equation (95) e 3 e 2 e 1 r(ξ B j ) g B2 g B3 Fig. 28 Illustration of the discrete coupling points for 2D-3D coupling along a single cross-section [63].
can be stated in matrix form as Furthermore, r S c,2D-3D , r B c,r,2D-3D and r B c,θ,2D-3D are the local residual vectors. Again, a linearization of the residual contributions with respect to the discrete beam-to-solid pair degrees of freedom is required for the Newton-Raphson algorithm. The linearization reads: The local contributions (96) and (98) to the global residual and the stiffness matrix, respectively, can be assembled in a straightforward manner and will not be stated here for the sake of brevity. As in the BTS-FULL mortar-type coupling, all derivatives explicitly stated in (98) are evaluated using forward automatic differentiation (FAD).
In practice, all integrals presented in this section are evaluated using a GPTS approach as illustrated in Figure 28, cf. [63]. At each Gauss-Legendre pointξ B i along the beam centerline, multiple equally spaced coupling points (illustrated with the symbol '×' in Figure 28) are defined along the circumference of the corresponding cross-section. Mechanically speaking, each coupling point is tied to the underlying solid via a linear penalty constraint.