A mortar-type finite element approach for embedding 1D beams into 3D solid volumes

In this work we present a novel computational method for embedding arbitrary curved one-dimensional (1D) fibers into three-dimensional (3D) solid volumes, as e.g. in fiber-reinforced materials. The fibers are explicitly modeled with highly efficient 1D geometrically exact beam finite elements, based on various types of geometrically nonlinear beam theories. The surrounding solid volume is modeled with 3D continuum (solid) elements. An embedded mortar-type approach is employed to enforce the kinematic coupling constraints between the beam elements and solid elements on non-matching meshes. This allows for very flexible mesh generation and simple material modeling procedures in the solid, since it can be discretized without having to capture for the reinforcements, while still being able to account for complex nonlinear effects due to the embedded fibers. Several numerical examples demonstrate the consistency, robustness and accuracy of the proposed method, as well as its applicability to rather complex fiber-reinforced structures of practical relevance.

behavior of engineering structures. In many cases, the reinforcements can be considered as being one-dimensional (1D), i.e. one dimension is much larger than the other two. Applications can be found in different fields, such as civil engineering, where steel reinforcements are embedded into concrete to improve its low tensile strength. In mechanical engineering, fiber-reinforced composites take advantage of fibers with high stiffness by embedding them inside a softer matrix material. This results in lightweight structures that are used in various applications, such as spacecrafts, boats, or sports equipment. Last but not least, also nature exploits the benefits of fiber-reinforced materials, as can be seen for example in arterial wall tissue with collagen fibers. Numerical simulation of such engineering and biomechanical structures is of high importance during the development and design phase, but it is also quite challenging.
Different modeling techniques exist to create a numerical model of the reinforced materials, almost all of them being based on the finite element method. From a mechanical point of view the matrix surrounding the beams is a three-dimensional (3D) continuum, which we will refer to as solid. In this work, we will denote the combined problem of arbitrary curved beams being embedded inside the solid volume as a beam-to-solid volume coupling problem. Figure 1 illustrates different beam-in-solid modeling techniques on the basis of the same physical problem of three fibers being embedded inside a material matrix, with the modeling complexity increasing from left to right. In the model shown in Figure 1(a), the stiffness contributions from the fibers and matrix are homogenized, thus resulting in an anisotropic material law for the combined volume [1,43]. In this case, the fibers are not explicitly modeled, and therefore, this is the most simple case of the models shown in Figure 1 regarding modeling effort and computational complexity. The main complexity in this approach lies in the accurate homogenization of the fibers and the matrix material. arXiv:1912.04153v1 [cs.CE] 9 Dec 2019 Figure 1(c) shows a model with the fiber volume explicitly cut out of the solid volume. This yields a fully 3D surfaceto-surface mesh tying problem between the fiber surfaces and the corresponding solid surfaces. In the shown model the coupling is realized by discretizing the solid and fibers with matching meshes. Alternatively, the interfaces could also be tied together with coupling methods for non-matching meshes, e.g. mortar finite element methods [29,[31][32][33]. Even in the non-matching case, the creation of the finite element mesh with explicit boundaries at the interface between beam surface and solid can be a non-trivial task. The extended finite element method (XFEM) [26] or immersed finite element methods [19,36] have been used to overcome this issue by implicitly defining the interface between beam surface and solid. Therefore, a very simple, in many cases even structured Cartesian finite element mesh can be employed. The drawback of those approaches are the numerically expensive cutting procedures required to implicitly model the interface. An approach as in Figure 1(c) is the closest to the real physical beam-to-solid volume coupling problem and is expected to provide very accurate solutions also close to the interface between fiber and matrix. Yet, it results in a complex model and an expensive numerical simulation, since resolving the fibers as 3D continua increases the system size by several orders of magnitude. This limits the applicability for large-scale engineering structures. Figure 1(b) shows a model with explicitly modeled fibers embedded into the matrix. In this case, the 1D reinforcements are modeled with a beam theory, which provides accurate and efficient numerical models for the fibers [21,24,25,35,38,39]. All kinematic fields of the beams are defined along the 1D centerline of the beam. The solid is modeled, and in particular it is meshed, without subtracting the beam volume from the solid volume, thus resulting in overlapping volumes. This introduces a modeling error, since in the physical problem no two material points can share the same spatial position. This modeling error is proportional to the fiber volume fraction as well as the stiffness ratio of fiber and matrix. The high fiber stiffness compared to the matrix stiffness in the considered cases reduces the influence of this modeling error. The new beam-to-solid volume coupling approach we present in this work follows the modeling ideas from Figure 1(b) and will exclusively use 1D beam formulations to model the fibers. The resulting beam-in-solid model boils down to a mixeddimensional 1D-3D coupling problem. Early work on 1D-3D coupling of structures has been carried out in the context of reinforced concrete in [27], with the restriction that the reinforcements have to align with a parameter coordinate of the solid element. In [7], this approach was extended to straight reinforcements with arbitrary directions relative to the solid elements, and in [10,13,34] also curved reinforcements are considered. All of those mentioned previous works do not introduce additional degrees of freedom for the reinforcements, but instead incorporate the beam stiffness contributions into the stiffness matrices of the solid elements. Alternatively, the beam degrees of freedom can be kept in the discrete system, which introduces the need for kinematic coupling constraints acting on the beams and solid [2,9,14,15,46]. A collocation method is used in [9] to couple 1D beams into a 3D matrix. In [15], a CutFEM approach is employed to embed 1D structural elements without bending stiffness into a 3D matrix material. The application of 1D-3D coupling can also be found in other fields than solid mechanics [8,16,17]. For example, vascular tumor growth is simulated in [16] by coupling the 1D vasculature to the surrounding 3D tissue. The approach recently presented in [18] combines the techniques from Figure 1(b) and 1(c) by using a 3D representation of the beams in zones of interest and 1D structural models otherwise.
In this work, we use C 1 -continuous geometrically exact beam finite elements [24]. Moreover, we propose an embedded 1D-3D mortar-type approach to model the coupling interaction between beam and solid finite elements. Specifically, a Lagrange multiplier field, representing a line load, is defined along the beam centerline to enforce the coupling constraints, similar to [16]. The coupling constraints are therefore formulated in a weak variational sense. The definition of interaction forces between beam and solid as a line load is a problem similar to the plane Kelvin problem of a line load acting on an infinite solid [12,28,40], which is illustrated in Figure 2. The exact solution to the Kelvin problem contains singularities in the stress and displacement fields close to the point of action of the line load. This has a major impact on the well-posedness and applicability of the proposed 1D-3D coupling method, and to the best of the author's knowledge, this aspect along with the resulting spatial convergence behavior will be discussed in detail for the first time. In the range of our modeling assumptions, i.e. relatively high beam stiffness compared to the solid stiffness as well as relatively small beam cross-section dimensions compared to the solid finite element sizes, the presented beam-to-solid volume coupling method is well-posed, yields very accurate results and exhibits optimal spatial convergence. In comparison to the available modeling techniques for thin fibers being embedded into a background material, this allows for an extremely efficient and simple model of the solid phase, while still being able to account for complex nonlinear effects due to the embedded fibers represented by 1D beam formulations.
The remainder of this paper is organized as follows: In Section 2, we derive the weak form of the quasi-static equilibrium equations for the beam-to-solid volume coupling problem via the principle of virtual work. This is done by combining the individual contributions from 3D solid structures, 1D beams and, in particular, the coupling/interaction terms between them. In Section 3, the finite element method is used to spatially discretize the weak form of the equilibrium equations. Further, the choice of suitable Lagrange multiplier basis functions as well as numerical integration techniques are discussed. The final discrete linearized system of equations is then derived by enforcing the coupling constraints in a weighted node-wise manner and by introducing a penalty regularization of the mortar method. Finally, numerical examples are given in Section 4. The examples are designed to assess the impact of modeling choices on the quality of the results, as well as to show the applicability of the presented methods to real-life engineering applications.

Problem formulation
We consider a 3D finite deformation beam-to-solid volume coupling problem as shown in Figure 3. For both the beam and the solid, a Cartesian frame {e 1 , e 2 , e 2 } is employed as a fixed frame of reference. The principle of virtual work (PVW) serves as basis for the employed finite element method. 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 constraint. Therefore, well-established formulations for the solid as well as the beam can be used without modifications. Without loss of generality, only quasi-static problems are considered in this work. This only impacts the virtual work contributions from the solid and the beam, but the coupling terms for the beam-to-solid volume coupling problem hold also for time-dependent problems. The fundamentals of both formulations as well as their individual contribution to the virtual work will be outlined in the next two sections. Finally the coupling between beam and solid will be described in detail in Section 2.3.

Solid formulation
The solid is modeled as a 3D continuum, represented by the open set Ω S 0 ⊂ R 3 in the reference configuration and by Ω S ⊂ R 3 in the deformed configuration. The reference surface ∂Ω S 0 can be divided into the Dirichlet and Neumann boundary surfaces, Γ u and Γ σ , respectively. In the current configuration they are denoted as γ u and γ σ . In the reference configuration, a material point on the solid can be identified by its reference position X S . The current position x S is related to the reference position through the displacement field u S via The variational formulation of the quasi-static balance equations serves as basis for the finite element method, resulting in the solid contribution δW S to the total virtual work. A Lagrangian formulation is used, i.e. all field variables refer to the reference configuration. Hence, the integration of the field variables is performed over the reference volume Ω S 0 and its boundary ∂Ω S 0 . Since the variation along the Dirichlet boundary Γ u vanishes, the only remaining surface integral in the variational formulation is over the Neumann boundary Γ σ . The virtual work δW S of the solid is given by Fig. 3 Notation of the finite deformation beam-to-solid volume coupling problem.
where δ denotes the variation of a quantity, S the second Piola-Kirchhoff stress tensor and E the energy-conjugate Green-Lagrange strain tensor. Contributions to the external virtual work δW S ext result from the prescribed body loadb and surface tractiont, both in defined the reference configuration. The Green-Lagrange strain tensor E is given as with F = ∂x S ∂X S being the material deformation gradient and I the 3D second-order identity tensor. For simplicity, we assume a hyperelastic material with the strain energy function Ψ(E), which relates to the second Piola-Kirchhoff stress tensor as follows: All the subsequent examples in Section 4 employ a hyperelastic material model for the solid, but this is by no means a requirement of beam-to-solid volume coupling, which for example can also be used for elasto-plastic solids.

Beam formulation
The beams used in this work are based on the geometrically exact beam theory, which in turn builds upon the kinematic assumption of plane, rigid cross-sections. Figure 4 shows the reference and current configuration of the beam without any additional kinematic assumptions. For illustration purposes, the reference configuration shows a straight beam, but unless stated otherwise, the presented beam theories can also be applied to beams with initial curvature. The complete beam kinematics can be defined by a centerline curve r(s) ∈ R 3 , connecting the cross-section centroids, and a field of right-handed orthonormal triads Λ(s) := The kinematic quantities X B , x B , u B ∈ R 3 , i.e. reference position, current position and displacement of an arbitrary point within the cross-section, are functions of the centerline coordinate s as well as the cross-section coordinates α, β ∈ R: x B (s, α, β) = r(s) + αg where u B r = r − r 0 is the displacement of the beam centerline.
In this work, three different geometrically exact beam theories are employed. Their basic kinematic assumptions as well as the corresponding internal elastic energy Π int,(·) will be stated in the following subsections. The beam contribution to the global virtual work reads where the term ∫ also differ for the three employed beam theories. A consistent and objective handling of the rotational variations contained in δW B ext,(·) is a non-trivial task. Since it is not the main aspect of the current work, the interested reader is referred to [24]. To improve readability of the following equations, a derivative with respect to the beam centerline coordinate s will be represented by (·) := ∂(·)/∂s throughout this section.

Simo-Reissner beam theory
Of the three beam theories considered in this work, the Simo-Reissner (SR) beam theory is the most general one, as it does not introduce additional kinematic constraints on the beam. This results in shear-deformable beams capturing six modes of deformation: axial strain, two bending modes, torsion and two shear modes. The cross-section kinematics can be described with six degrees of freedom: the spatial position of the cross-section r(s) and its rotation vector ψ(s) ∈ R 3 , which defines the cross-section triad Λ(s) = Λ(ψ(s)) based on the well-known Rodrigues formula [24]. The internal elastic energy of the beam is given as where axial tension and shear strains are represented by the material deformation measure Γ := Γ(ψ, r ) = Λ T r − e 1 ∈ R 3 , while torsion and bending are represented by the material curvature vector Ω ∈ R 3 , which in turn follows from Ω × a = Λ T Λ a ∀ a ∈ R 3 . Using the rotation vector parameterization of the triad field Λ(ψ(s)) as discussed above, the resulting curvature vector can be formulated as a function of ψ and ψ , i.e. Ω = Ω(ψ, ψ ). The constitutive matrices C F and C M are defined as where E is the Young's modulus, G the shear modulus, A the cross-section area, A 2 and A 3 the effective shear areas, and I T , I 2 , I 3 are the polar and planar second moments of area, respectively.

Kirchhoff-Love beam theory
The Kirchhoff-Love (KL) theory introduces an additional kinematic constraint, restraining the shear deformation of the beam. This is equivalent to the requirement that the first cross-section basis vector g 1 is parallel to the centerline tangent r , or While the position of the cross-section is described in the same manner as for the Simo-Reissner beam, the additional constraints reduce the number of independent rotations to one, thus a total of four degrees of freedom remain to fully describe the cross-section. The sole remaining rotational degree of freedom ϕ(s) ∈ R describes the twist rotating of the cross-section around the tangent vector r measured with respect to a properly defined reference triad Λ ref (r ), such that the cross-section triad can be described as a function of the centerline tangent and the twist, Λ(s) = Λ(r (s), ϕ(s)).
A detailed overview, how to parametrize the twist degree of freedom, can be found in [24]. The curvature of the beam centerline is described with the Frenet-Serret vector which only depends on the beam centerline. The definition of the curvature contains second derivatives of the beam centerline, therefore resulting in the smoothness requirement of C 1 continuous centerlines. Defining the material curvature vector Ω identical to the Simo-Reissner case above, it can be formulated as a function of ϕ, ϕ , r and κ for the Kirchhoff-Love case [24], i.e. Ω = Ω(ϕ, ϕ , r , κ). The internal energy for the Kirchhoff-Love beam is Therein, = r − 1 is the axial tension of the beam.

Torsion-free beam theory
The torsion-free (TF) beam formulation considered in this work was first proposed in [22] and extended in [23]. It represents a special case of the Kirchhoff-Love beam theory. For certain properties of the problem, i.e. straight undeformed beams with axisymmetric cross-sections and no external torsional moments, it can be shown that the static equilibrium configurations resulting from the Kirchhoff-Love beam theory are characterized by (exactly) vanishing torsion [22]. The fact that these requirements are fulfilled in many practically relevant systems, and also in most of the examples considered in this work, justifies and motivates the application of this type of beam element formulation. Compared to the Kirchhoff-Love beam, the twist degree of freedom is not present anymore and the beam can be completely described by its centerline position, i.e. three degrees of freedom per cross-section. Since the discrete representation and algorithmic treatment of large rotations is the main complexity of geometrically nonlinear beam theories, the employed torsion-free beam theory, which can completely abstain from any rotational degrees of freedom, is particularly appealing and easy to handle. The internal energy of the torsion-free beam reads with the scalar curvature κ = κ .

Beam-to-solid volume coupling
In the beam-to-solid volume coupling problem shown in Figure 3, the beam is embedded inside the solid volume. The most natural choice for the coupling conditions is to couple the beam surface ∂Ω B 0 to the solid volume Ω S 0 . However, there is no explicit surface in the solid domain, to define the coupling conditions on. Therefore, this is a surface-tovolume (2D-3D) coupling problem, i.e. the beam surface is embedded into the background solid volume. The coupling constraints are formulated in the reference configuration and read with Γ 2D-3D c = ∂Ω B 0 being the coupling surface. The Lagrange multiplier method is employed to impose the coupling constraint. A Lagrange multiplier vector field λ 2D-3D (s, α, β) ∈ R 3 is defined on Γ 2D-3D c , which can be interpreted as the negative interface tractions acting on the beam surface. Contributions to the total virtual work are the coupling interface contribution and the variational form of the coupling constraints This leads to a saddle point-type weak formulation of the 2D-3D beam-to-solid volume coupling problem: The integrals in equations (16) and (17) are evaluated on the coupling surface Γ 2D-3D c , which requires a computationally expensive numerical integration of δW 2D-3D c and δW 2D-3D λ . For the inherent assumption in this work, that the cross-section dimensions of the beam are small compared to the other dimensions of the beam-to-solid volume coupling problem, we can approximate the surface integrals as line integrals along the beam axis Ω B L,0 . These line integrals can be evaluated very efficiently. The approximation changes the physical coupling dimensionality applied to the beam-tosolid volume coupling model from surface-to-volume (2D-3D) to a line-to-volume (1D-3D) coupling. The new coupling Since this is a significant change in the mathematical description of the mechanical model, the implications of this choice will be discussed in several remarks at the end of this section. The approximated variational coupling terms read Here λ 1D-3D (s) ∈ R 3 is a new Lagrange multiplier field defined along the beam centerline. We point out that λ 2D-3D and λ 1D-3D have different physical dimensionality and, accordingly, also different units: the first one is a surface load, while the latter one represents a line load. The final PVW for the 1D-3D beam-to-solid volume coupling problem reads For improved readability, the superscript 1D-3D for the lineto-volume coupling terms will be omitted from now on.

Remark 2.1
In the previous considerations, it was assumed, without loss of generality, that the beam consists of a single fiber which lies completely inside the solid. The derived equations also hold if the beam sticks out of the solid volume.
In this case the coupling integrals are not evaluated on the complete beam domain, but instead only on the portion of the beam centerline inside the solid. The impact on the numerical integration will discussed in Section 3.3.

Remark 2.2
With the definition of the line-to-volume coupling terms in (19) and (20), the coupling is now exclusively formulated through the beam centerline position, which decouples the cross-section rotations Λ from the solid deformations. In particular, relative rotations between the crosssection and the solid around the tangent vector r are not restrained. At first glance, this might be considered as a rather coarse approximation for certain physical systems such as fiber-reinforced composite materials, where fibers are e.g. molded / glued into a matrix such that all modes of relative motion are blocked. However, in our target applications, the main contributions to the internal energy of the beams and the mechanical resistance of the overall structure stem from bending and axial tension of the fibers, therefore justifying the choice to neglect the coupling of cross-section rotations. Additionally, this will lead to coupling terms, which only contain the centerline degrees of freedom and are independent of the actual beam theory. This allows for an easy adaptation of our beam-to-solid volume coupling method to different beam theories. As a further consequence, the rotation of beam fibers around their centerline might be unconstrained, possibly yielding a singular linear system to solve. In practice, this repairable deficiency is limited to static analyses and the undeformed configuration. As a remedy, one either imposes Dirichlet boundary conditions on at least one of the twist degrees of freedom or uses standard linear solvers with deflation capabilities to properly exclude such nullspace modes. The problem is cured as soon as the beam centerlines have deformed, i.e. usually after the first Newton step. This discussion also underlines an advantage coming with the torsionfree beam theory: The corresponding beam finite elements do not have any rotational degrees of freedom, and consequently, such rigid body modes cannot occur.

Remark 2.3
Another aspect to be addressed, when switching from 2D-3D to 1D-3D coupling, is the introduction of singular solutions. From a mechanical point of view, the line-to-volume coupling is equivalent to a line load inside the solid. This is a generalized version of the Kelvin problem [12,28,40], which consists of an infinite solid loaded with an embedded line load. The analytical solution for the Kelvin problem has a singularity at the line load point of action, not only in the stress field, but also in the displacement field. This has a significant impact on the spatial convergence of the finite element discretization and will be discussed in detail in Section 4.2.

Spatial discretization and numerical integration
An isoparametric finite element discretization is employed to approximate the continuous fields for geometry, displacement as well as virtual displacement. The interpolation of the positions and displacements in the solid domain is given by and Here, N k ∈ R is the finite element shape function for the solid node k, x S k ∈ R 3 and d S k ∈ R 3 are the nodal reference position and displacement, respectively. The total number of solid nodes is n S . The variables ξ S , η S and ζ S are the 3D coordinates of the solid finite element parameter space.
Kirchhoff-Love and torsion-free beam elements require a C 1 -continuous centerline interpolation, which is realized with third-order Hermite polynomials [24,41]. An objective and path-independent interpolation of the rotational field Λ (e) h (ξ B ) along the beam centerline is a non-trivial task and will not be discussed here, since the rotations do not appear in the coupling terms anyway. A comprehensive overview on this topic can be found in [24]. The resulting beam element has two centerline nodes with six degrees of freedom per node, i.e. three positional and three tangential degrees of freedom. Due to its superior numerical properties, this discretization scheme is also used for the Simo-Reissner beam element, as derived in [20]. The beam centerline reference position and displacement are interpolated by and where H r l ∈ R and H t l ∈ R denote the Hermite shape functions for the positional and tangential degrees of freedom for the beam node l. Both shape functions are a function of the scalar beam centerline parameter coordinate ξ B . The discrete vectors x B,r l , x B,t l ∈ R 3 are the reference position and tangent, respectively. The discrete degrees of freedom d B,r l , d B,t l ∈ R 3 denote the nodal displacements and tangent increments. The total number of beam centerline nodes is n B . To improve readability of the following derivations, the beam centerline displacement is redefined in the following way The discrete nodal displacement vector d B l now contains all centerline degrees of freedom for the node l.
Employing a mortar-type coupling approach, the Lagrange multipliers are also approximated with a finite element interpolation [3,30,44]. The continuous Lagrange multiplier field λ is defined along the beam centerline. Therefore, the Lagrange multiplier interpolation is defined along the 1D beam elements. All subsequent integration is performed on the domain Γ B c,h , which is the projection of the beam centerline domain Ω B L,0 onto the beam finite element function space. In the nomenclature of classical contact mechanics, the beam would be considered the slave side, and the solid the master side. The approximated Lagrange multiplier field reads where Φ j ∈ R is the shape function for the discrete Lagrange multiplier vector λ j ∈ R 3 at node j. The total number of discrete Lagrange multiplier nodes is n λ , which is not necessarily equal to n B . The shape function Φ j is a function of the scalar beam centerline parameter coordinate ξ B . Note that even though the Lagrange multipliers are defined along the beam centerline domain, the displacement shape functions H l will not be used to interpolate the Lagrange multiplier field. An adequate choice of Lagrange multiplier shape functions will be discussed in Section 3.1. The nodal discrete unknowns d S k , d B l and λ j are assembled into the global degrees of freedom vectors d S , d B and λ.
Insertion of the finite element approximations (23), (25) and (29) into the variational form of the coupling constraints (20) gives where h defines a suitable projection from a point on the beam centerline to the corresponding point in the solid volume. In the previous equation, two local matrices with mass matrix-like structure can be identified: There, D (j,l) describes the coupling between the Lagrange multiplier node j and the beam node l and M (j,k) describes the coupling between the Lagrange multiplier node j and the solid node k. They can be assembled into global, so called mortar matrices D ∈ R 3n λ ×6n B and M ∈ R 3n λ ×3n S , which both are rectangular in general. A similar expression containing D and M can also be derived for the virtual work δW c,h of the coupling forces. All in all, the coupling contributions to the weak form can now be stated in global matrix form .
Here, f S c and f B c are the vectors with the discretized coupling forces acting on the solid and beam degrees of freedom, respectively. The vector g c contains the discretized constraint equations and its entries can be interpreted as the relative displacement between beam centerline and solid weighted with the Lagrange multiplier shape functions. Inserting all discretized variables into (21) gives the discrete nonlinear system of equations for the quasi-static beam-to-solid volume coupling problem: Here, f S int and f B int are the internal force vectors of the solid and beam, respectively. The Newton-Raphson algorithm is used to obtain solutions to the system of nonlinear equations. Therefore, a linearization of equations (35) to (37) with respect to the global unknowns d S and d B has to be derived. The linearized system of equations with saddle point structure reads: where are the stiffness matrices of the solid and beam, respectively.

Lagrange multiplier shape functions
The choice of Lagrange multiplier shape functions is important for the mathematical properties of the discretized system, since the discrete Lagrange multiplier bases, i.e. shape functions, must fulfill an inf-sup condition with the displacement field [5]. In the context of surface-to-surface contact or mesh tying in solid mechanics, this is a well studied-topic. However, in the considered beam-to-solid volume coupling problem, we employ Hermite polynomials as primary shape functions for the slave side, i.e. the beam, which is unusual compared to the standard surface-to-surface case. Additionally, we deal with an embedded 1D-3D coupling, i.e. there is no explicit curve representation in the solid mesh to match the beam centerline, which can lead to stability issues [37]. The numerical experiments in Section 4.2 and 4.3 carefully evaluate the influence of different Lagrange multiplier bases on the numerical properties of the beam-to-solid volume coupling problem.
Since the Lagrange multipliers are defined on the beam centerline, a natural choice in the spirit of the mortar method would be to use the same shape functions as for the beam elements, i.e. third-order C 1 -continuous Hermite polynomials. However, for neighboring beam elements with equal length the integral over the Hermite shape functions associated with the tangential degrees of freedom becomes zero. This can lead to numerical difficulties in the constraint enforcement. Therefore, in this work, standard Lagrangian shape functions are used to interpolate the Lagrange multiplier field. Three different types of shape functions will be compared: linear, quadratic and cubic. In surface-to-surface mortar methods, the use of stable lower order interpolations for the Lagrange multipliers compared to the displacement interpolation order was already successfully explored in [31,33].

Remark 3.1
The previous derivations are given for the case, where the constraint equations are fulfilled in a truly weak (variational) sense. In Section 4, this mortar-type coupling will be compared to a classical Gauss-point-to-segment (GPTS) coupling approach. In the GPTS coupling, the strong form of the constraint equations (15) is fulfilled at each Gausspoint along the beam, i.e. a discrete 3D Lagrange multiplier vector λ GPTS j is defined at each Gauss-point in the sense of a collocation method. However, GPTS coupling can also be interpreted as a special case of the mortar coupling, namely if the Lagrange multiplier field is interpolated as Here, δ is the Dirac delta distribution with the property ∫ . The position and weight of the j-th Gauss-point are denoted withξ B j andw j , respectively.

Enforcement of constraint equations
The constraint equations (15) are discretized with a mortar coupling approach using Lagrange multipliers, this resulting in a mixed formulation. However, due to certain drawbacks, e.g. an increased system size compared to the uncoupled problem and a saddle point structure, (38) will not be solved directly here to obtain solutions to the beam-to-solid volume coupling problem. Instead, the penalty method is used to obtain approximate solutions of (38). This results in a formulation that is purely displacement-based and does not contain any additional variables. The main idea behind this well-known penalty regularization of the mortar method is to allow a relaxation of the discretized coupling constraints g c = 0 in the form Herein, the Lagrange multipliers are no longer independent variables, but well-defined functions of the beam and solid displacements. They can subsequently be removed from the global system of equations. In (40), ∈ R + is the penalty parameter and it is clear that for → ∞, (40) becomes equivalent to (37). The entries in the weighted relative displacement vector g c are proportional to the support of the corresponding Lagrange multiplier shape function, i.e. they depend on the beam element length. If unaccounted for, this dependency would result in a violation of the basic patch tests presented in Section 4.1. To resolve this problem, the relaxation of the constraints in (40) is additionally multiplied with the inverse of the diagonal nodal scaling matrix κ, similar to the approach in [45]. The local scaling matrix for the Lagrange multiplier node j is defined by and is assembled into the global scaling matrix κ. With the penalty approach, the coupling forces f S c and f B c can be stated as, With this the final global system of equations (38) becomes: Here, the number of global unknowns is the same as in the uncoupled case. An additional effect of the penalty-regularized version of the mortar method is the elimination of the saddle point structure in the stiffness matrix. However, there are some drawbacks of the penalty approach. The constraint equations are violated by definition, which only can be reduced with higher penalty parameters, but this in turn leads to an ill-conditioned tangential system matrix. Therefore, it is desirable to choose a penalty parameter that results in a sufficiently accurate solution of the constraint equations, but also limits unwanted numerical effects. The influence of the penalty parameter in practice will be discussed in detail in Section 4.3.

Numerical Integration
The beam-to-solid volume coupling contributions to the global system of equations are all calculated via integration over the beam domain in the reference configuration, cf. (31) and (32). Numerical integration, namely a Gauss-Legendre quadrature, is used to evaluate the coupling matrices D and M and the scaling matrix κ during the finite element simulation. An accurate numerical evaluation of the coupling integrals is absolutely essential to pass basic consistency tests, such as the patch tests in Section 4.1. The integrands in D and κ solely contain fields defined along the beam centerline, namely the beam displacements and the Lagrange multipliers. If the Jacobian ∂ r 0,h /∂s along the beam element is constant, the integrand is of polynomial form and the numerical integration is exact, if enough quadrature points are used. In the cases considered in this work, the maximal polynomial degree of the integrand in D and κ is 6, i.e. third-order beam shape functions and third-order Lagrange multiplier shape functions. Therefore, 4 Gauss-Legendre points are needed for the numerical integration to be exact. The integrand of M contains fields defined along the beam centerline as well as the solid volume. In Figure 5, it can be seen that the evaluation of the solid shape functions along the beam centerline results in a general nonlinear function which contains socalled weak discontinuities, i.e. kinks at the points where the beam crosses between solid elements, and strong discontinuities, i.e. jumps at points where the beam sticks out of the solid volume. Moreover, the continuous parts of the integrand in M are not of polynomial degree. To still guarantee high accuracy of numerical integration for the integrand in M, two different algorithms will be investigated and compared, cf. Figure 6. They are illustrated in Figure 6. Element-based integration uses a fixed number of Gauss-points per beam element. The only exception occurs at strong discontinuities, where the integration is only performed for the part of the beam element inside the solid volume. In segment-based integration, the integration domain along the beam element is split into multiple segments, such that the integrand in the individual segments does not contain any kinks. Each segment is then integrated with a fixed number of Gauss-points.
The global coupling matrices only depend on the initial configuration of the beam-to-solid volume coupling problem, i.e. they remain constant over the course of the simulation. From a computational point of view, it makes sense to evalu-ate the coupling matrices D, M and κ once and store them for subsequent Newton iterations and time steps. Nevertheless, it is important to address the impact of the different numerical integration schemes with regard to computational performance and accuracy. Independent of the integration scheme used, each Gauss-point evaluation requires the solution of a local nonlinear system of equations, i.e. the projection of the point on the beam centerline into the solid finite element parameter space. For element-based integration, the evaluation time for the coupling terms is more or less proportional to the number of Gauss-points used. Since the integrand contains kinks, a relatively high number of Gauss-points is necessary to obtain a sufficiently accurate numerical integration. On the other hand, segment-based integration requires calculation of the intersections of the beam elements with the solid surfaces. This intersection operation also requires the solution of local nonlinear systems. The total number of intersections, which have to be calculated, depends on the mesh configuration and cannot be quantified in a general manner. The advantage of the segment-based integration is that the integrands over a segment are smooth, see the left part of Figure 5, and an acceptable integration error can be obtained with a reasonable number of Gauss-points. Unless stated otherwise, all the examples in this work use 6 Gauss-points per integration segment. A direct comparison of the two integration schemes regarding evaluation time is difficult, as the times depend on the mesh configuration of the individual problem. In [11], an elaborate comparison of different numerical integration algorithms for mortar methods is given. It should be stated that, in general, due to the non-polynomial integrand in M both integration schemes cannot integrate M exactly. Nevertheless, the segment-based integration has clear advantages: the accuracy of its numerical integration is independent of the beam-to-solid element length ratio and a higher accuracy can be achieved with the same global number of Gauss points.

Examples
The following examples are chosen to evaluate the different beam-to-solid volume coupling methods proposed in this work, cf. Table 1, and to demonstrate their accuracy and robustness for the simulation of challenging engineering applications. All simulations are performed with our in-house parallel multi-physics research code BACI.

Patch Tests
Patch tests are a well-established tool to investigate the consistency of finite element formulations. In the realm of solid mechanics, they are typically used to check that the finite element method is able to exactly represent a constant stress   [32], where a constant traction across non-conforming interfaces can exactly be represented via a mortar mesh finite element approach. However, the choice of a suitable patch test for beam-to-solid volume coupling is not straightforward, as this is technically not a discretization of two domains with non-matching meshes at the interfaces, but rather an embedded problem of a 1D curve (the beam) lying inside a 3D volume (the solid). Figure 7 shows the first patch test presented here. It consists of a solid cuboid Ω S with two embedded straight beams B1 and B2, where Ω B1 and Ω B2 , the domains of both beams, occupy the same spatial domain. No surface loads or body forces are applied on the solid, while constant line loads with a magnitudet act in opposite directions ±e 3 on the beams.

Beams inside a solid volume
Therefore, the opposing loads on the two beams cancel each other out and in sum the two beams transfer no loads to the solid. This gives the trivial solution for the solid displacement field u S = 0 and the constant solution −u B1 = u B2 = e 3t /ε for the beam displacements, where ε is the contact stiffness between the beams. This patch test uses the proposed beamto-solid volume coupling method to couple both beams to the solid. By doing so, all interactions between the beams are transfered via the solid domain, resulting in a patch test like problem for beam-to-solid volume coupling. This test case will be used to assess the influence of discretization and integration error on the performance of the proposed beam-to-solid volume coupling method. The solid is discretized with 4 × 4 × 7 eight-noded, first order hexahedral elements (hex8). Simo-Reissner beam elements are used to represent both beams B1 and B2, which are discretized with 5 and 7 equidistant elements, respectively. Mortar coupling is applied between the beam centerline and the solid, with a linear interpolation of the Lagrange multiplier field λ along the beam elements. To circumvent numerical problems in the solution of the resulting linear system of equations, the solid is constrained such that all six rigid body modes of the system are eliminated. Additionally, any rotation of the first nodes of the two beams is constrained to prevent a rigid body rotation of the beams around their axes. The magnitude of the line loads on the beams ist = 5 N/m. For the given geometry, there is no discretization error, since the chosen shape functions for the beams and the solid are able to exactly represent both geometry and numerical solution. To assess the numerical integration error, the problem is solved once with element-based integration of the mortar coupling terms and once with segment-based integration. Figure 8(a) shows the result obtained with element-based integration and 6 Gauss points per beam element. Clearly, the solution is not exact, as the solid is not stress-free and the deformation of the beams is not constant, thus resulting in non-vanishing curvatures along the beams. Figure 8(b) shows the results obtained with segment-based integration of the mortar coupling terms, where each segment is integrated with 6 Gauss points. In this case, the numerical results exactly match the analytical solution up to machine precision, which confirms the vanishing integration error for segment-based integration.
A second patch test is set up similar to the first one, with the straight beams being replaced by two helix-shaped beams. The helix has the following geometrical parameters: A radius of 0.45 m, three turns with a pitch of 9/5 m and a right handed screw type. In this case, the beams B1 and B2 are discretized with 23 and 31 elements, respectively. The employed C 1 -continuous Hermite polynomials used for the beam centerline interpolation can not represent the helix geometry exactly, which results in two slightly different geometries of the beams and different arc lengths of the two helices, thus introducing a discretization error. In order for the two beams to be in equilibrium, the loadt on beam B2 is scaled with a factor of 0.999318, to correct for the different beam lengths. In this case, the beams can not perform a rigid body motion when coupled to the solid. Therefore, only the six rigid body modes of the solid are constrained. All other parameters are equal to the previously described patch test. Figure 9(a) shows the results with element-based integration of the mortar coupling terms. Similar to the previous scenario, one can see non-vanishing stresses in the solid and curvature oscillations in the beams. In this case, also the result with segment-based integration, shown in Figure 9(b), does not match the analytical results up to machine precision, because of the previously described discretization error. However, when comparing the quantitative results, one can see that the influence of the numerical integration error for element-based integration is about one order of magnitude larger than the discretization error, which confirms that element-based integration introduces a significant additional integration error.
The presented results were all calculated with our mortar beam-to-solid volume coupling approach and first-order interpolation of the Lagrange multipliers. Quantitatively, the results change only slightly if a GPTS (1D-3D) approach is used or if a different interpolation scheme for the Lagrange multipliers is applied. Therefore, the conclusions obtained from the shown examples, i.e. the importance of an accurate numerical integration of the beam-to-solid volume coupling terms and the superiority of segment-over element-based integration, can be applied to all aforementioned cases.

Strong discontinuities in beam-to-solid volume coupling
To check the ability of the proposed methods to handle strong discontinuities, i.e. a beam sticking out of a solid domain, two more patch tests are introduced. Both problems consist of a solid cube and a straight beam which starts inside of the cube and ends outside of it. In the first case, the beam intersects a face of the solid, in the second one it intersects an edge. All solid degrees of freedom are constrained and a constant line load −t e 3 = −1 N/m e 3 is applied only to the part of the beam inside of the cube. Similar to the previous patch test the analytical solution for the beam displacement is u B = −t/ε e 3 . Furthermore, the beam can only be in equilibrium if the coupling interface traction is λ =t e 3 = 1 N/m e 3 . For reasons of simplicity, the solid cube is discretized with a single hex8 element, the beam with a single Simo-Reissner element. In this case, segment-and element-based integration are identical to each other, as both schemes will result in the same integration points and weights. Segmentation has to be performed at the point where the beam exits the solid volume. The patch tests are analyzed once with a Gauss-point-to-segment approach and once with mortar coupling using a linear interpolation of the Lagrange multipliers. Figure 10 shows the results. For the Gauss-point-to-segment method, the coupling forces at the integration points are illustrated and it can be observed that they are exact up to machine precision. The same holds true for the Lagrange multiplier interface tractions in the mortar case. The discrete Lagrange multipliers should not be confused with discrete nodal loads on the beam element, as the Lagrange multiplier field is only integrated on the beam segment that resides inside the solid. This underlines the importance of segmentation at solid surfaces, i.e. proper treatment of strong discontinuities.

Spatial Convergence
The following numerical example investigates the spatial convergence properties of the beam-to-solid volume coupling method as well the validity of the evaluation of the coupling terms along the beam centerline instead of the beam surface, i.e. our fundamental 1D-3D modeling assumption. The considered problem is shown in Figure 11.  The spatial convergence behavior of the different coupling methods will be analyzed with respect to the 2 displacement error error is computed relative to a reference finite element solution obtained with 2D-3D (surface-to-volume) coupling. The 2D-3D coupling is realized with a GPTS approach as illustrated in Figure 12. At each Gauss-Legendre pointξ B j along the beam centerline, multiple equally spaced coupling points (illustrated with the symbol '×' in Figure 12) are inserted along the circumference of the corresponding cross-section. The coupling points are constrained to the circumference of the cross-section. All coupling points along a single crosssection ran rotate around r (ξ B j ), i.e. the resulting moment around r (ξ B j ) vanishes, which corresponds to the kinematic assumptions of the torsion-free beam theory. The coupling points are tied to the underlying solid mesh via a linear penalty constraint. In all of the following results obtained with 2D-3D coupling, 6 Gauss-Legendre point in axial direction and 128 integration points in circumferential direction are used. This ensures a sufficiently accurate numerical evaluation of the 2D-3D coupling terms in order to hold as reference solution, and the chosen penalty parameter does not lead to unwanted stiffening effects.
The solid block is meshed with first-order hex8 solid elements, with an element size h solid . The rod is discretized with torsion-free beam finite elements with a length of h beam = 2.5h solid . The penalty parameter for all 1D-3D coupling methods is 100 N/m 2 , for 2D-3D GPTS coupling it is 100 N/m 3 . Additionally to the previously described integration rule for 2D-3D coupling, all 1D-3D coupling schemes in this example are evaluated with segment-based integration and 6 Gauss points per segment. These values are chosen according to Section 4.3 in order to avoid unwanted contact locking effects. For models purely consisting of either first-order solid or third-order beam elements, the expected convergence rate of the 2-error is O(h 2 ) and O(h 4 ), respectively. The expected convergence rate for the coupled problem is thus the lower of the two, i.e. O(h 2 ). Figure 13 shows the convergence plot of the coupled structure with different coupling methods. The 2D-3D GPTS coupling scheme exhibits the expected convergence rate of O(h 2 ) for the entire dataset. All 1D-3D coupling schemes behave very similar to each other. For coarse meshes, the expected optimal convergence order O(h 2 ) can be observed. At around h solid = 0.12 m the convergence behavior of all 1D-3D coupling methods has a kink, and for smaller element sizes the error does not decrease any further, it even slightly increases. The bottom part of Figure 13 illustrates the solid mesh size compared to the beam cross-section at different points in the convergence plot. In the case of a coupling along the beam surface (2D-3D), the beam interacts with all solid elements along its surface. For beam-to-centerline coupling (1D-3D), the beam only interacts with the solid elements along its centerline. For finer discretizations, the influence of the different interactions becomes more evident, which materializes in the kink in the convergence plot. This result gives rise to a very important finding, namely that our 1D-3D coupling is valid down to a certain element size, i.e. up to the kink in the convergence plot. Exemplarily, the beam tip displacement of the 2D-3D reference solution is 0.19009 m. With our 1D-3D beam-to-solid volume coupling method the tip displacement at the kink in the convergence plot (mesh B from Figure 13) is 0.18895 m, which amounts to a relative error of approximately 0.5%. The coupling interactions of the 2D-3D and 1D-3D schemes are shown in Figure 14. The critical solid element size, i.e. up to which the 1D-3D coupling is accurate, depends on a number of different parameters and can not be given in closed form. However, for the problems considered in this work, i.e. rather stiff beams and soft solids, a rule of thumb can be given: the solid element size should not be smaller than the beam cross-section diameter. Keeping in mind the envisaged applications, one can conclude that this does not pose any restrictions on our beam-to-solid volume coupling methods, but is perfectly in line with their modeling goal.
Remark 4.1 Consider a plane problem of a beam crosssection coupled with a solid finite element mesh as depicted in Figure 15. As long as the cross-section fully lies within a single solid element, i.e. the cross-section diameter is smaller than the solid element size, the resulting nodal forces on the solid nodes should be independent of the used coupling scheme -as long as the resultants of the 1D-3D and 2D-3D coupling are equivalent. Obviously this is an idealized setting, but this still underlines and nicely illustrates the validity of our 1D-3D coupling approach down to a solid element size of about the cross-section diameter.

Influence of the penalty parameter
In this example, the influence of the penalty parameter as well as the beam-to-solid element length ratio is investigated. The analyzed problem is the same as in Section 4.2, now with a fixed solid element length h solid = 0.25 m. The model is simulated with different penalty parameters and beam-to-solid element length ratios. To quantify the differences between results obtained with different parameters, the 2-errors relative to the same reference solution as used in Section 4.2 are compared. The results are shown in Figure 16. Each of the four plots represents a fixed beam-to-solid element length ratio. The penalty parameter is plotted on the abscissa. The line style identifies the employed coupling scheme. For both element and segment-based integration, 6 integration points are used per element and segment, respectively. The desired behavior for an increasing penalty parameter is a convergence towards the exact fulfillment of the constraint equations, i.e. the solution of (38). In the presented plots, this corresponds to a horizontal line for high penalty parameters. For  all element length ratios, the GPTS scheme with segmentbased integration exhibits an increasing error for increasing penalty parameters. The GPTS scheme with element-based integration behaves better for high element length ratios, but as the beam length gets closer to the solid element size, the same behavior can be observed. This effect is sometimes referred to as contact locking and occurs due to an overconstraining of the system, i.e. too many discrete coupling constraints are enforced, and as a result, the coupling discretization becomes too stiff. Mathematically, this is related to a violation of the discrete inf-sup condition [5]. This effect is especially distinct for GPTS schemes, where each Gauss point represents three coupling constraints, i.e. the number of discrete coupling constraints depends on the integration scheme used. A smaller number of Gauss points can usually improve the contact locking properties for GPTS schemes, but this in turn can lead to the non-fulfillment of the patch tests given in Section 4.1. The beam-to-solid volume coupling mortar schemes behave better: for element length ratios of 10 and 5 no locking can be observed at all. For smaller element length ratios the schemes with quadratic and cubic interpolation also show signs of contact locking. Linear interpolations of the Lagrange multipliers do not show such behavior for the considered element length ratios. By using a lower order interpolation of the Lagrange multipliers, the number of constraints is reduced, which explains the better behavior of the lower-order Lagrange multiplier interpolations regarding contact locking. The employed numerical integration scheme does not affect the contact locking behavior of mortar beam-to-solid volume coupling methods, as the number of coupling constraints is independent of the number of Gauss points used.
The results show that a GPTS-based coupling discretization tends to be prone to spurious contact locking effects. A linear interpolation of the Lagrange multipliers within a mortar-based coupling discretization, as suggested in this contribution, is the most robust coupling scheme regarding the choice of the penalty parameter.

Fiber-Reinforced Composite Plate
In this example, a fiber-reinforced composite plate is modeled with the proposed beam-to-solid volume coupling method and the results are compared to a homogenized approach. Figure 17 shows the problem setup of a two-layer composite plate. The plate has a length and width of 2 m and 1 m, respectively. The layer buildup is asymmetric: it consists of two layers with fiber directions of 45 • and −45 • , each with a thickness of 0.02 m. A hyperelastic Saint Venant-Kirchhoff material model (E = 10 N/m 2 , ν = 0.3) is used to model the matrix material. The fibers are modeled as torsion-free beams (E = 1000 N/m 2 , ν = 0) with circular cross-sections (radius r = 0.045 m). Figure 17 shows the fiber placement in the layers, which results in a fiber volume ratio of 0.25. At one of its short ends the plate is clamped in e 1 and e 3 direction, and a surface Neumann load p in e 1 direction of 2.5 N/m 2 is applied at the other short end. The matrix is modeled with 288 eight-noded solid-shell elements [4,42] and the fibers with 1498 torsion-free beam elements, respectively. On average, the beam-to-solid element length ratio is about 2.5. Mortar coupling with linear interpolation of the Lagrange multiplier shape functions and a penalty parameter of 1000 N/m is used to couple the beams to the solid. Segment-based integration is used to evaluate the coupling terms. All boundary conditions are exclusively applied to the solid-shell elements.  Figure 18(a) shows the deformed plate, where for illustration purposes only three quarters of the solid elements are visualized. Due to its asymmetric layer buildup, the plate deforms out of the e 1 − e 2 plane, even tough all applied loads and boundary conditions are exclusively in-plane. Figure 18(b) shows only the beam elements and a vector plot of the negative discrete nodal values of the coupling tractions calculated with (40). The largest coupling tractions occur at the boundary of the plate, especially at the corners. These coupling tractions will be used to gain insight on fiber pull-out and related composite damage phenomena in future research. Such information cannot be obtained at all from a homogenized theory.
The same plate is also modeled using a homogenized approach. Each layer is modeled with a transversely isotropic material, this representing a homogenization of the fibers and matrix in that layer. As is common practice, the material properties for the transversely isotropic material are calculated according to a homogenization approach for linear strains, cf. [43]. For the nonlinear simulation of the plate, a combination of a purely isotropic hyperelastic material and a transversely orthotropic hyperelastic material is employed, cf. [6]. Each layer is modeled with 288 eight-noded solidshell elements, thus resulting in a total of 576 finite elements for the homogenized model. In Figure 19, the deformations of the mid-plane at the right end (e 1 = 2 m) are compared to the results obtained with our new beam-to-solid volume coupling method. Only for larger loads, there is a tiny discrepancy between the different methods, which can be attributed to a number of factors, e.g. the different strain measurements used in the beam and the homogenized solid, or small scale effects in the composite that can not be resolved by the continuum model. Nevertheless, the results are in excellent agreement with each other, which underlines the general applicability of our beam-to-solid volume coupling method to fiber-reinforced composites.

Remark 4.2
The presented beam-to-solid volume coupling model of the composite plate consists of 1,950 solid degrees of freedom and 10,992 (torsion-free) beam finite element degrees of freedom. This example can also be modeled with Kirchhoff-Love beam elements, which yields the same numerical results up to machine precision, due to exactly vanishing torsion [22]. However, the number of beam degrees of freedom for the Kirchhoff-Love model increases by about 30% to 14,322, thus justifying and encouraging the application of torsion-free beam element formulations if the underlying assumptions are met.

Fiber-Reinforced Pipe
The final numerical example is a fiber-reinforced pipe under pressure. The problem setup, illustrated in Figure 20(a), consists of a pipe modeled with a Neo-Hookean material law (E = 10 N/m 2 , ν = 0.3). The pipe is 2 m long and has an inner and outer radius of 0.9 m and 1 m, respectively. It is reinforced with Simo-Reissner beams (E = 1000 N/m 2 , ν = 0) as also shown in Figure 20(a). The cross-section radius of the beams is 0.04 m. The inner surface of the pipe is loaded with a Neumann surface pressure p of up to 2.5 N/m 2 . At the top and bottom, symmetry boundary conditions are applied to the pipe as well as to the beams. For further symmetry reasons, only a quarter of the depicted pipe is actually simulated with the following element numbers referring to the quarter model. Coupling between the beams and the solid is realized with our mortar approach and linear Lagrange multiplier shape functions with a penalty parameter of 1000 N/m and segment-based integration. The pipe is discretized with 225 C 1 -continuous isogeometric solid elements (based on second-order NURBS) and 45 Simo-Reissner beam elements. Figure 20(b) shows the deformed configuration of the pipe. The expected stiffening effects of the beams onto the structure can clearly be seen. In-between the beam reinforcements, the relatively soft pipe exhibits larger displacements. Although only qualitative in nature, this example could obviously not be modeled with a homogenized approach and illustrates a very interesting problem class for the new beamto-solid volume coupling method. Even tough all the previous derivations and examples used first-order interpolation of the solid finite elements, this example also showcases the straightforward applicability of our beam-to-solid volume coupling method to higher-order and even C 1 -continuous solid interpolations. This allows for a coupling of beam and solid fields with the same order of interpolation continuity.

Conclusion
In this work, we have proposed new modeling techniques for the coupling of 1D continua embedded inside full 3D continua. Two different finite element-based coupling schemes have been introduced, a Gauss-point-to-segment (GPTS) and an embedded mortar-type approach. The resulting constraint equations of both schemes are enforced via a penalty method, and in the case of the mortar-type approach, the penalty regularization was performed in a weighted node-wise manner. For the mortar-type method, different discrete Lagrange multiplier bases were investigated. Moreover, different numerical integration methods of the coupling terms were compared. Several numerical experiments have been conducted to assess the behavior of the different schemes regarding the choice of the penalty parameter and the numerical integration of the coupling terms. For relevant physical application scenarios of the beam-to-solid volume coupling method, i.e. relatively slender and stiff fibers compared to the surrounding matrix material, the validity of the fundamental modeling assumption of 1D-3D coupling has been verified, and its optimal spatial convergence behavior has been shown numerically. Furthermore, the results underline the importance of an accurate numerical integration of the coupling terms as provided only by carefully chosen segmentation schemes. Overall, the embedded mortar-type discretization with linear interpolation of the discrete Lagrange multiplier basis emerges as the better modeling choice due to its superior robustness regarding the choice of the penalty parameter, the  beam element to solid element length ratio and its optimal spatial convergence properties. Future work will focus on the extension of the presented beam-to-solid volume coupling approach to beam-to-solid surface coupling as well as beam-to-solid surface contact. Another topic of interest for further research is to make use of the improved local resolution of the numerical solution close to the fiber-matrix interface, compared to homogenized approaches, for analyzing progressive damage and failure phenomena in fiber-reinforced materials, such as fiber pullout and the onset of delamination.