Coupling 2D continuum and beam elements: a mixed formulation for avoiding spurious stresses

This paper presents a novel approach to coupling beam and solid elements. Connecting standard beam elements with solid elements results in a situation where the solid part of the model covers the cross-sectional deformation, while at the beam part assumptions only consider the rigid body movement of the cross-section. Therefore, kinematic assumptions do not include cross-section deformations such as warping and contraction. Due to this restriction, spurious stresses occur at the transition zone. To circumvent this problem, this contribution introduces a mixed hybrid transition element based on an extended Hu–Washizu functional. The extension allows the cross-section to contract and warp. Compared with other approaches that include these phenomena, precomputation of the warping function is unnecessary. The present approach considers the cross-sectional deformation using local variables. On the solid interface, there are only displacement degrees of freedom (depending on the solid discretization), whereas on the beam interface, there are two displacement degrees of freedom and one rotation. This allows beam-like parts of the solid model to be replaced with standard beam elements without affecting the overall accuracy of the model. Besides solid–beam coupling, the proposed formulation can also be used to apply beam-like boundary conditions as well as stress resultants to the solid interface. For the sake of simplicity and without restricting generality, the underlying element formulation is introduced for the 2D linear elastic case. Numerical examples demonstrate that the transition element causes no spurious stress distributions on the solid part. The results of full solid models are compared with mixed beam–solid models.


Introduction
In many cases, transitioning from a continuum theory to a beam theory is useful for reducing the number of degrees of freedom of a model. Furthermore, in applications such as homogenization considering macroscopic beam or shell kinematics [5,7] and multiscale modelling considering solid and structural elements [19], questions arise on how to apply the boundary conditions without introducing spurious stresses at the boundary of the solid model. This transition is difficult because of the different kinematic assumptions of a beam and a solid element.
B Simon Klarmann klarmann@lbb.rwth-aachen.de 1 Over the past decades, different approaches have been proposed for achieving the transition between solids and structural elements have been proposed. These approaches can be categorized into multipoint constraints (MPC) and transition elements; see for example [6,12]. The present contribution proposes a transition element; therefore, a brief overview of the existing relevant works is provided. The first transition element for 2D solids and beams was proposed in [13], which modified the element geometry and associated shape functions. This approach was extended to shell-solid coupling in [14] and to non-linear deformation in [15]. In [4], a transition element was introduced to couple 3D solid elements with beam elements, again using a special element geometry to describe the displacement field. To avoid locking inside the element, the authors in [3] introduced a formulation based on a Hellinger-Reissner functional. The kinematic assumption did not include warping or contraction of the cross-section. Later, the authors in [18] proposed a transition element for coupling shell structures with beam elements. The element considered contraction of the crosssection, by allowing the movement of each cross-sectional point on a straight line toward the center of gravity of the cross-section. This formulation was extended in [1] to thermoelastic analysis. Furthermore, the authors introduced a warping function that needed to be evaluated a priori. In [2], a mesh-geometry-based solution was introduced by adapting the mesh to cross-sections' geometry and using beam elements to achieve the coupling. A transition element based on asymptotic beam theory was proposed in [12]. The aforementioned approaches either do not include warping, require an a priori calculation of a warping function or require the assumption of linear elastic material behavior.
Another approach is multipoint constraint coupling. Here, constraint equations are formulated to connect the solid nodes to either beam or shell kinematics. The method requires assumptions on how the load is distributed over the solid face. The coupling 3D solids and 1D beams was presented in [10], which was extended to more general 1D, 2D, and 3D coupling in [9,11]. To find the correct stress distribution on the solid side, [8,19] have analyzed the force distribution in separate numerical analyses. With the first contributions in 1979 and the newest ones in 2022, the topic is already quite old but not entirely solved. Summing up, there are two limit cases to achieve the coupling between solids and beams. One assumes continuity in the displacement field, allowing the cross-section only to perform a rigid body translation; the second assumes continuity of the stresses but is limited to constant and linear stress distribution over the cross-section.
All papers discussed before either do not contain crosssectional deformation at the interface; their resulting stress distributions do not fulfill the stress boundary conditions; require a separate calculation of the stress distribution, or require an additional warping function. In contrast, our element allows for a general cross-section deformation without needing other precalculated cross-section parameters or functions. This is one of the key features of the method presented in this paper.
This contribution aims to couple a continuum displacement field with the displacement field of a 2D beam by reducing the continuum displacement field to two displacements and one rotation. The proposed element does not evoke spurious stresses, and the resulting stress state is equivalent to the one of a long and slender solid model.
As two displacements and one rotation are the general degrees of freedom of a beam element, either a Timoshenko or a Bernoulli beam can be connected to the 2D solid model. The assumptions made, lead to a formulation that is closer to a solid element than a beam element. The novel aspects and key features of the proposed formulation are as follows: • No a priori calculation of the stress distribution or cross-section deformation is required, e. g. cross-section warping; • General cross-section deformations are determined by the numerical analysis method and are essential to avoid spurious stresses at the interface; • Although the model considers arbitrary cross-section deformations represented by local variables, no additional global degrees of freedom are necessary for the resulting finite element.
The resulting finite element formulation has only displacement and rotational degrees of freedom, has a symmetric stiffness matrix, and requires only the input parameters of a classical 2D solid element.
To focus on the main idea of the formulation, which is introduced in Sect. 3, this paper addresses the 2D, homogeneous, linear elastic case and thoroughly investigates it.
The remainder of this paper is structured as follows. Section 2 presents the kinematic assumptions and motivates the modification of the weak form. Section 3 introduces the weak from of equilibrium based on an extended Hu-Washizu functional to allow cross-sectional deformations. The weak form contains additional correction strains. Section 4 deals with the derivation of their specific shape functions. Section 5 presents details on the implementation of the proposed element formulation. Finally, Sect. 6 thoroughly investigates the performance of the proposed element formulation using various numerical examples.

Kinematic assumptions
In this section, kinematic assumptions for the transition element are derived. The goal is to connect the displacement field u S ∈ R 2 of a 2D solid at the face S with those of a 2D beam consisting of displacements u B0 ∈ R 2 and rotation β ∈ R. It is assumed that the structure has a constant width b in the out-off-plane direction. This constant part b is omitted hereinafter. This approach deals with a 2D geometry described by the position vectors X S , X B , X B0 , and X ∈ R 2 ; the local basis A i ∈ R 2 ; the global basis e i ∈ R 2 ; and the displacement vectors u S , u B , u B0 , and u G ∈ R 2 . Index S represents the 2D solid quantities, index B the beam quantities, and index B0 quantities on the beam reference axis. Vector X is the position vector inside the element and u G is the displacement vector inside the element in the global coordinate system e i (see Fig. 1). The kinematic assumptions allow a deformation of the 2D solid face S while the beam assumptions state that the beam face B stays plane and does not contract. Furthermore, small displacements and rotations are assumed as well as that the faces S and B being initially parallel. The impact of the restriction on the solution, namely To connect the beam with the 2D solid kinematics, Fig.  1 presents the initial structure. Establishing the connection between the 2D solid part (gray: S ) and the beam (line: B ) is achieved using a transition zone (yellow part: I ) of length which consists of the 2D solid kinematic assumptions on the left face S and the beam assumptions on the right face B . The beam node does not have to coincide with the center of gravity of the cross-section. The center of gravity is not required to be known. Furthermore, each point of the 2D solid face S is described by position vector X S , on the beam face B with position vector X B and inside the transition zone with position vector X. The connection point to the beam has the position vector X B0 . Like in classical beam theories, a local orthonormal coordinate system A i (i = 1, 2) with the coordinates ξ i ∈ R is attached to the beam node. The vector A 1 is the normal of the faces S and B , A 2 is parallel to these faces. The faces are assumed to be plane and represent the cross-section of the beam with a height of h. Furthermore, in the undeformed configuration the faces S and B are the same except that B is shifted by the distance in the A 1 -direction.
The position vectors X S (ξ 2 ) on S are explicitly given, while those on B are described by The beam parametrization in Eq. (1) means that each point on the face B is described by the distance ξ 2 in direction A 2 from the reference point X B0 . Compared with Fig. 1, it is clear that ξ 2 ranges over the height h of B . Using X S and X B , a position X inside the transition zone is given by a linear interpolation over the length of the element: On the 2D solid face S , the global displacement vector of a point P S is u S (ξ 2 ) and on the beam face B of a point P B it is u B (ξ 2 ). With respect to beam kinematics and considering small rotations, u B (ξ 2 ) reads as follows: In Eq. (3), u B0 and β are the displacement vector and the rotation at X B0 , respectively. Assuming a linear interpolation using the displacements u S and u B leads to the global displacements at X: Using the orthogonal transformation R = e i ⊗ A i to project global displacements given in Eq. (4) into the local coordinate system leads to Equation (5) finally represents the displacement field in the local coordinate system A i . With the assumption of small deformation, the strain field is given using Equation (6) describes continuum strains in the local coordinate system A i of the element. Herein (. . .) ,i , i = 1, 2 represents the derivatives w.r.t. local coordinates ξ i . It is important to note that the derived strains contain the assumption of a plane and an inextensible cross-section at the face B .
As the kinematics are composed of 2D solid and beam kinematics, it is essential to make a connection to pure beam kinematics. This motivates the extension of the weak form in Eq. (17) and in later Sect. 4. To achieve this goal, the 2D solid displacements u S are compared with averaged displacement u C and a rotation β C of the face S . They are defined as follows: In Eqs. (7) and (8) the displacements u S (ξ 2 ) are functions of ξ 2 , while u C and β C are constant through the face S . They describe the movement of S in terms of beam kinematics. With these assumptions, the displacements u S can be written as follows: Equation (9) means that the movement of a point on S can be described by using the beam kinematics enriched by a general deformation of the cross-section ω C (ξ 2 ). Inserting Eqs. (9) into (5) leads to the following equation Here, the quantitiesū andβ are given as follows: Using Eqs. (10) with (11) and inserting it into Eq. (6) leads to the following strains: With Eq. (12), it is clear that within the kinematic assumptions, the classical Timoshenko beam strains (stretch ε, curvature κ, and shearing γ ) are included: The matrix E 1 in Eq. (13) represents the relation between the beam strains and continuum strains. The additional strain parts due to the deformable 2D solid surface S are collected inε. It is important to note thatε contains the assumption of the inextensible and plane beam cross-section of face B . Precisely this part is relaxed in Eq. (17) and motivates the construction of the correction strains in Eq. (19). Based on these kinematical considerations, the weak form of equilibrium can be derived.

Weak form of equilibrium-Mixed formulation
In this section, the weak form of equilibrium is introduced.
To derive the weak form of equilibrium for the transition element, the weak form of the whole structure is split into three parts first: In Eq. (14), the three parts are given by g S , g B ,and g I , namely the 2D solid part, beam part, and transition zone, respectively. Regarding g S and g B , it is assumed that only continuity in displacements is required; note that for g B , the displacements are composed of displacements of the reference line and rotations. Therefore, classical beam and 2D solid formulation can be used for these two parts. As mentioned in Sect. 2, the kinematics assume that the cross-section is inextensible and plane at the beam side B . To remove this constraint, a mixed formulation is used. The starting point is the standard Hu-Washizu functional for the transition zone of the model: The initial weak form of equilibriumḡ I of the transition element in Eq. (15) depends on the displacements u or geometric strains ε according to Eq. (6); the physical strains ε p and physical stresses σ p ; as well as their virtual counterparts δu, δε p , and δσ p . The constitutive relation between the physical strains and stresses is expressed by using the matrix C, which in case of plane-stress assumptions is as follows: In Eq. (16), E is the Young's modulus, ν is Poisson's ratio, and b is the width of the structure (perpendicular to the e 1 , e 2 -plane).
Furthermore, the term δε T σ p in Eq. (15) is responsible for transferring the forces from the beam node to the 2D solid surface. To obtain the correct distribution of stresses on the 2D solid surface, σ p must have the correct shape. As σ p is related by the constitutive model to ε p through the term C ε p −σ p , this can only be achieved if the physical strains ε p have the correct shape over the height of the transition zone. This leads to the important part * of Eq. (15), which enforces the balance between the physical strains ε p and geometric strains ε. For the geometric strains ε, the relation according to Eq. (6) is used. Therefore, the assumption of an inextensible and plane cross-section on B is still included. The idea to solve this problem is the enrichment of the geometric strains ε with additional correction strains ε c in a manner that the physical strains ε p have the correct shape. This leads to the following modified version of Eq. (15): Performing an integration by parts of Eq. (17) results in the strong form: The strong form [Eq. (18)] consists of the equilibrium neglecting the body forces, the relaxed kinematics, and the constitutive relations inside the domain I . The last statement says that all higher-order stresses related to the correction strains ε c should be zero in the domain ξ ∈ [− , 0]. In the modified weak form in Eq. (17), special care needs to be taken for the correction strains ε c which is discussed in the next section.

Design of correction strains
Due to the extension of the weak form of equilibrium from Eqs. (15) to (17), special care must be taken regarding ε c . If ε c is assumed to be an arbitrary function, then it would accumulate all of the strains, which would result in ε p = 0 and σ p = 0. As a result, the beam node can move independently of the solid surface. To determine how to reduce the solution space of ε c , it is crucial to examine the relation of the kinematic assumption; see Eqs. (8) to (13). The assumption of an inextensible and plane cross-section is included with an additional deviation of ω C from S . The face B can only perform a rigid translation and rotation according to Eq.
(3). To achieve the correct deformation state inside the element, the correction strains ε C must remove the constraint of an inextensible and plane cross-section on B . To construct these strains, the same approach is used as that in [16,17]. Following this approach with some slight modifications, the continuum strains are represented using the following equation: In Eq. (19), the term ε c = E 2 ε cω is used to remove the constraint of an inextensible and plane cross-section on the beam side B . Like the expression for the classical beam strains, E 2 ε cω is composed of functions over the crosssection arranged in E 2 , which are scaled by scalar values in ε cω that represent higher-order beam strains. The functions in E 2 are constructed using the general deformation functions ω 1 and ω 2 , where ω 1 represents cross-sectional warping in the local A 1 -direction and ω 2 represents the contraction of the cross-section in the local A 2 -direction: In the definition of E 2 in Eq. (20), it is clear that the strains ε c are introduced as independent strains and are not defined through additional displacements on the beam surface, despite the fact that ω 1 and ω 2 can be interpreted as superposed displacements on B . For ε c to not accumulate all of the strains, the general displacement functions ω 1 and ω 2 must fulfill the following orthogonality conditions: The constraints in Eq. (21) result from the interpretation of ω 1 and ω 2 as superposed displacements on B and compared with Eq. (9). Within the reformulated 2D solid displacements [Eq. (9)], the rigid body translation and rotation are described using the beam kinematics, and therefore they must be removed from ω 1 and ω 2 . As demonstrated in [16,17], the representation of ε c = E 2 ε cω with the constraints in Eq. (21) allow the enriching of classical beam strains to represent a complete continuum strain state.

Finite element formulation
This section derives the matrices of the transition element. Due to the orthogonality conditions [Eq. (21)] of the general displacement functions ω 1 and ω 2 , their shape functions consist of a local and a non-local part. This will result in a non-standard finite element; the composition of the finite element is depicted in Fig. 2. In the addressed 2D case, the element is subdivided into regular sub-elements e Si based on

Approximation of the geometry
Due to the different kinematic assumptions of 2D solid and beam kinematics as well as the regular geometry, separating the interpolation of geometry and solution is reasonable. To achieve this, a distinction is made between vertices and nodes. Vertices represent points with specific coordinates, whereas nodes represent the degrees of freedom associated with either a vertex or an edge. It is sufficient to interpolate the geometry with linear edge elements on the solid surface by employing the shape functions: Using the locally attached parametric coordinate system η 1 , η 2 and the shape functions in Eq. (22), a coordinate on the surface S is described as follows: According to Eq. (23), X h S is the surface coordinate approximated by interpolating the two coordinates X e i j , j = 1, 2 of the associated edge element e i . If the n S element vertex coordinates are arranged inX ∈ R 2(n s +1) , vertex coordinates of an edge X e i j can be extracted using The matrices a e i X in Eq. (24) are comparable to the assembly matrix, which here is used to extract the vertex coordinates. With the approximated surface coordinate X h S on the 2D solid surface, the corresponding coordinate on the beam side can be evaluated using In Eq. (25) is the length of the transition element and A 1 is equivalent to the normal vector on the 2D solid surface S . As each sub-element is assumed to have the shape of a rectangle, the normal vector A 1 in Eq. (25) can be computed using the vertex coordinates X e i of an arbitrary edge e i in Eq. (24): With Eqs. (23) and (25), a general point of sub-element e Si is given as follows: Furthermore, with Eqs. (25) and (26), the local coordinate ξ 2 reads as follows: In Eq. (28) X B0 is the coordinate of the beam vertex n B . Finally, based on Eq. (27), the Jacobean with respect to ξ i evaluates to Here, l e i is the length of the edge e i in Eq. (29). Using Eq.
(29), the derivatives of the shape functions with respect to the local coordinates ξ 1 , ξ 2 can be computed as follows: As the local and parametric coordinate axes coincide, only the derivatives specified in Eq. (30) exist, and the remaining derivatives are zero. Furthermore, the parameter i specifies the shape function number.

Approximation of the displacement field
Next, the approximation of the displacement field is addressed. Nodes with degrees of freedom v i = u x , u y T i are assigned to the geometric objects (vertices and edges) and interpolated inside the sub-elements. Vertex nodes on the solid surface are associated with the linear shape functions, as in Eq. (22). The higher-order shape functions are achieved by hierarchically extending the approximation order by assigning nodes to the edges and associating them with Lobatto shape functions: Using this hierarchical extension, a sub-parametric concept for the element formulation is employed. On the beam vertex n B , two nodes are defined -one with degrees of free- representing the displacements and one with v B2 = [β, 0] T B2 representing the rotation. Distributing n N nodes on the surface S results in a vector and v B2 are the beam nodes with the displacement and rotational degrees of freedom, respectively. Defining the polynomial order of the shape functions of the displacement field on h S as n p leads to n N = n p · nel S + 1 surface nodes. The n p + 1 nodes of an edge and the two beam nodes are extracted as follows: Like a e i X , the matrix a e i S extracts the relevant degrees of freedom from V. With Eq. (32) and the shape functions Eqs. (22) and (31), the approximation of the displacement field at the surface h S is as follows: On the beam side, the displacement approximation, using Eq. (28) to evaluate ξ 2 , is as follows: Finally, with Eqs. (33) and (34), the displacement approximation inside one sub-element e Si is derived using With Eq. (35), the approximation of the displacement field is complete.

Approximation of the geometric strains
With the representation of the displacement field in Eq. (35), the geometric strains in Eq. (6) can be written as follows: The positions of the entries in B in Eq. (36) depend on the order of the nodes in V as well as the position in the cross-section, for example the current sub-element during the numerical integration. Only the entries B B1 and B B2 are always present in the same position. Considering Eq. (6), B I evaluates to and B B1 and B B2 evaluate to With Eqs. (37) and (38), the representation of ε h in Eq. (36) is complete.

Approximation of the correction strain field
To derive the shape functions for the correction strain field ε c , the functions ω 1 and ω 2 are interpreted as displacement fields. Thus, the same distribution of degrees of freedom as for the 2D solid displacement field on S is assumed. Furthermore, rigid body translations and rotation are already included in the kinematic assumptions and must be restricted for ω 1 and ω 2 . This leads to the following non-local warping shape functions: To fulfill the orthogonality condition in Eq.
It is necessary to highlight that the index i in Eq. (40) associates the functions H i with a specific node and does not indicate the shape function order. The order of the local shape functions N i are the same as those of the displacement field. Additionally, the shape functions only vary over the thickness of the cross-section and are constant in the length-direction ξ 1 of the transition element. To evaluate the unknowns c 1i and c 2i , the system of equations must be solved for each shape function. It is important to note that the integration of the left-hand side in Eq. (41) is performed over the whole cross-section and not only locally inside one sub-element e Sj . The vector of contraction functions ω 2 is approximated with As there is only one constraint equation for shape functions ω 2 in Eq. (39), they are extended with a single constant parameter. Therefore, these functions are given byĤ i = N i +ĉ i . Looking at Eq. (20), only the derivatives w.r.t. ξ 2 are required, causing the constant extension to drop out. As a result, local shape functions in each sub-element can be used to represent the functions ω h 2 andĤ i = N i is chosen. By inserting Eqs. (39) and (42) into (20), the matrix E 2 is defined as follows: Noteworthily, the first and third rows of E h 2 in Eq. (43) always have n N −2 entries, while the second row has n N −1 entries, because the linear dependent equations must be removed due to constraints in Eq. (41). With the matrix E h 2 , the correction strain field reads as follows: With Eq. (44), the approximation of the correction strains is completed and only the approximations of the physical strains and stresses are missing.

Approximation of the physical strains and stresses
As the displacement field is approximated linearly in the length direction of the transition element, the physical strains ε h p and stresses σ h p are approximated constant in this direction. In the thickness direction, the components ε h x and σ h x are interpolated using the same shape function order n p as for the displacement field on the 2D solid surface S , while the remaining parts are interpolated with order n p − 1. The approximation reads as follows: In Eq. (45),ε p andσ p are the degrees of freedom of the physical stresses and strains, while A e Si ε and A e Si σ include the associated shape functions with respect to the current sub-element e Si . If the degrees of freedom of the current sub-element are extracted withε e Si p = a e Si εε p andσ e Si p = a e Si εσ p , the matrices read as follows: It is advantageous to use Legendre shape functions L i , i = 0, . . . , n p for Eq. (46). These shape functions are as follows: The reason for this can be found in Eq. (49) as the integration of the term A T ε A σ leads to a diagonal matrix if Legendre functions are employed.

Discretized weak form of equilibrium
To evaluate the weak form of equilibrium of the element Eq. (17), the same shape functions are chosen for the associated virtual quantities. Therefore, they read as follows: With these at hand, the element stiffness matrix resulting from Eq. (17) reads as follows: The integration of Eq. (49) is performed using the sum over each sub-element and a local Gauss quadrature. After the integration over the transition zone I , the full element stiffness matrix reads as follows: The matrices G, H, L, M result from the integrals in Eq. (49). As only continuity in displacements and rotations is required, the local variables can be condensed in Eq. (49). Considering the system of equations As the matrices M and H are both invertible, and H is also a diagonal matrix due to the properties of the Legendre polynomials, Eq. (55) can be reformulated tô Inserting Eq. (56) into (54) and solving it forε c leads tô Finally, inserting Eqs. (56) and (57) into (51) leads to the condensed element stiffness matrix: Notably, the resulting element stiffness matrix k e according to Eq. (58) has exactly three non-zero eigenvalues, independent of the number of degrees of freedom of the transition element. One can think of these three non-zero eigenvalues as the transfer of the normal force, shear force, and moment from the beam point to the 2D solid surface. A numerical investigation is performed in Sect. 5.7.

Stiffness entries at the beam point and eigenvalues of the element
To obtain insight into the possible performance of the proposed element, this section investigates entries of the stiffness matrix with a focus on the part that represents the beam nodes. At first, the entries of the stiffness matrix k e of the proposed transition element [Eq. (58)] are compared with the stiffness matrix k e Timo of a Timoshenko beam element. Using a reduced integration with linear shape functions, said stiffness matrix has the following entries for the comparable degrees of freedom: Entries in Eq. (59) consist of the element length L e , which is equivalent to the transition element length , tension stiffness E A, shear stiffness G A S considering a suitable shear correction factor, bending stiffness E I , and the distance of the center of gravity of the cross-section to the beam reference line s y . To compare the results, a height of h = 1, width of b = 1, Young's modulus E = 100, Poisson's ratio ν = 0.3, a shear correction factor of κ q = 5/6, and s y = −2 are assumed. Assuming a transition element length of = 0.001, which is equivalent to the element length L e = 0.001, the relative error, considering 40 sub-elements over the height, evaluates to Except for the terms containing the shear stiffness, the error is only affected by numerical inaccuracy.
As already mentioned, the element only has three nonzero eigenvalues. It is clear that rigid body modes cannot be related to the zero eigenvalues as in a classical finite element; instead, the three non-zero eigenvalues represent the number of forces transferred from the beam to the 2D solid surface. As a result, it is possible to relate the eigenvalues to the stiffness properties E A/ , E I / , and G A S / as depicted in Fig. 3. The eigenvalues are first multiplied with the element length and then normalized by the associated cross-sectional parameter. All eigenvalues converge against their related value of Timoshenko's theory. For all examples, the geometry of the transition elements is interpolated with linear shape functions, while the displacements are approximated with bi-quadratic shape functions. The 2D solid part is modeled with a displacement-based formulation. For the beam part, a Timoshenko beam element with linear shape functions is chosen. Like the 2D solid element, the beam element is based on a pure displacement formulation, but reduced integration is employed to avoid shear locking.
For a brief overview, the examples and their purposes are listed in Table 1. The first and third examples consist of two pure 2D solid models (Tension, shear, and bending test-Benchmark test; and Arc). Their aim is to demonstrate the quality of the resulting stress distribution over the crosssection where the proposed transition element is applied.
In these examples, the transition element is used to impose boundary conditions on the 2D solid models as well as to connect two 2D solid models with only three degrees of freedom. To demonstrate the suitability of the proposed element in mixed models, consisting of 2D beam and 2D solid elements, the second example is a simple frame structure. The  Simple frame Application of the proposed element to sub-structure modelling.

Arc
Testing the limits due to the assumptions of an initially straight transition zone.

Remark 1 (Plot)
It is important to note that in plots of deformed meshes, the proposed transition element is represented assuming the kinematic assumptions introduced in Sect. 2. This means, warping, and contraction on the beam side h B of the transition element are not considered. In addition, the stresses are not calculated for representation inside the transition element. This allows an easier distinction between 2D solids and the transition element in the figures as the latter is depicted in gray.

Tension, shear, and bending test -Benchmark test
In this example a simple cantilever beam is chosen. The aim is to test the different beam deformation modes, such as tension, bending, and shearing of the beam. The focus is on the area where the transition element connects the beam with the 2D solid model. The systems are depicted in Fig. 4. To distinguish between the solid part and the transition zone, they are identified by S and I , respectively. The length L S and transition element length together lead to the total length of the structure L, resulting in a total length of L S + in the classically clamped structures sys. a and sys. b. Applying the supports as well as the loading using two transition elements results in sys. c with length L S + 2 . The structure has a height of h and the beam node has a distance of s y to the center of gravity. If not otherwise specified, meshes of the solid part consist of elements with a width and height of e. Examples are chosen to demonstrate the performance of the transition element in achieving the connection between solid kinematics and beam kinematics without introducing spurious stress distributions over the cross-section. The material parameters for all of the examples are Young's modulus E = 100 and Poisson's ratio ν = 0.3. The same applies to the geometry parameters, for which the total length L is chosen as 10, which means either L = L S + = 10 if the left side is fully clamped or L = L S + 2 = 10 if the transition element is applied to both sides. Furthermore, the height and width of the structure are set as h = 1 and b = 1. If not specified otherwise, the distance to the center of gravity is s y = 0. In case of sys. a and sys. b, the displacement boundary conditions are directly applied on the 2D solid elements. In case of sys. a, the cross-section can contract on the left side, but vertical displacements are only supported through a single node, which introduces a stress concentration at that node. Using sys. b, the left side can no longer contract, but vertical forces of the support are distributed over the height of the cross-section. F x = 1, s y = 0 First, a tension rod is investigated. This means that the geometry is chosen according to Fig. 4, the beam node is at the center of gravity (s y = 0), and the structure is loaded in the x-direction. The load F x = 1 is applied to the beam node on the right side of the structure.

Tension test -Load
In Fig. 5, the numerical result is compared with the analytical result. The analytical solution can be evaluated with u ana. = F x L/E/b/h according to the beam theory. To observe the impact of the proposed element formulation, a numerical investigation is performed with the three different boundary conditions. The first diagram (Fig. 5) depicts the relative displacement error of the beam node with varying lengths of the proposed transition element. As expected for case ν = 0, the fully clamped version sys. b has a much higher error, whereas in case of sys. a and sys. c only numerical errors are visible. The reason for this is that in sys. b, the cross-section cannot contract on the left side due to the boundary conditions. Furthermore, the length of the transition element has no impact on the solution.
Next, the mesh refinement is investigated, and the resulting relative error is depicted in Fig. 6. It is important to note that the error is depicted over the element size e. Consequently, the horizontal axis ranges from fine (left) to coarse meshes (right). For sys b. a slight convergent behavior can be observed, when reducing the element size e. However, compared with the analytical solution the error remains high due to the boundary effects with ν = 0.3. When using the transition element on both sides (sys. c) or allowing the contraction of the cross-section on the left side (sys. a), the error is in the range of the numerical accuracy. This is due to the fact that with a reduced element size e, the number of degrees of freedom increases, resulting in a large equation system and thus a larger numerical error.
To highlight the properties of the proposed transition element, the structure is reduced to a length of = h = 1 and a load of F x = 1, and the deformed geometries with all stress components inside the 2D solid elements are depicted in Fig.  7. To create the plots, the transition element length is set as = 0.1 and the element size is set as e = 0.05. To enable directly distinguishing between the transition element and the 2D solid elements, the stresses are only plotted inside the solid part and the transition zone is presented in gray. As is clearly visible, the stress state inside the solid part is not affected by the boundary conditions nor by the proposed transition element. As expected, the normal stress σ x = 1 and the other two stress components are σ y = σ xy = 0.
By contrast, Fig. 8 depicts the stress and deformation state if sys. b is used. On the clamped side, the cross-section cannot contract, leading to distortions in the normal stress σ x and the non-zero stresses σ y and σ xy .
To investigate the impact of the length L of the system, it is varied from L = 1 to L = 15 in Fig. 9. Therefore, an L/h-ratio ranging from 1 to 15 is investigated. In case of sys. a and sys. c the resulting displacement error compared with the analytical beam solution is numerically zero. In case the structure is fully clamped on the left side (sys. b) the error remains high but is reduced with increasing length L.

Tension test, off-centered -Load F x = 1, s y = h
As mentioned in the theory part, the distance of the beam vertex to the center of gravity must not be known. To determine whether the effects of an off-centered load are covered correctly, the system is modified such that s y = h (see Fig.  4). Here, only sys. c is considered as the impact of boundary conditions on the numerical solution has already been investigated for a comparable load case. The material parameters are the same as before and the length of the transition element is = 0.1. As the deformation involves only ten-sion and bending, the analytical solution of the Bernoulli and Timoshenko beam theories is the same in this example. The reference solution u ana. for the horizontal displacement u x , vertical displacement u y , and rotation β according to the beam theory is as follows: As expected with sys. c, the results in Fig. 10 indicate that the error is only affected by numerical accuracy. All three results, namely the tip deflections u x and u y and tip rotation β, are evaluated accurately compared with the analytical solution. Errors grow due to reduced element sizes and therefore an increasing number of degrees of freedom. This indicates that only numerical accuracy has an impact on the solution.
In Fig. 11, the solution is compared with the analytical one when the distance to the center of gravity is modified. For this investigation an element size of e = 0.25 and a transition element length of = 0.01 are chosen. Again, the error is in the range of 10 −12 and 10 −10 due to numerical accuracy.
The deformed structure with normal stresses σ x is displayed in Fig. 12. For the plot, the structure has a total length of L = 10, height of h = 1, and distance of beam vertices to the center of gravity of s y = 2.5 h. In the figure, the beam vertices are depicted as black points. The lower-left Fig. 10 Tension test off-centered -Relative error with fixed distance to the center of gravity (s y = h) and varying element sizes for sys. c Fig. 11 Tension test off-centered -Relative error between the FE solution and analytical solution with varying s y for sys. c Fig. 12 Tension test off-centered -Deformed structure with normal load and s y = 2.5 h. The contour plot depicts normal stresses σ x Fig. 13 Tension test off-centered -Stresses over the cross-section one is completely clamped, whereas the lower-right one is loaded by a single load F x = 1. In the deformed configuration, the lower right vertex is moved to the right due to the stretch of the structure and rotation of the cross-section. As the off-center normal load introduces an additional moment, the vertex moves upwards. The resulting normal stresses are not affected by the boundary conditions. Finally, Fig. 13 presents the stress distribution over the height ranging from the bottom −0.5 to the top 0.5. The cut is performed at the connection line between the 2D solids and the transition element on the loaded side. As expected, the stresses σ y and σ xy are zero, while the distribution of σ x is linear over the cross-section. Note: Kinks in the σ x stresses result from the stress projection on a linear geometry element at the boundary.

Shear and bending test -Load F y = 1
After investigating the tension load case, the load direction is changed to a load F y = 1 perpendicular to the beam axis. The material and geometric properties are maintained as before. Unless otherwise specified, the transition element length is set as = 0.1. The analytical solution for the vertical displacement u y using Timoshenko's beam theory evaluates to: To compare the results, the analytical solution is evaluated using Timoshenko's beam assumptions with a shear correction factor of κ = 5/6. As in the tension case, the solutions of sys. a and sys. c are the same with respect to numerical accuracy; at first, the In later examples, it becomes clear that the length of the transition element can be interpreted as a mesh parameter. Therefore, for the comparison of the different boundary conditions, it is set as = 0.01 to reduce its influence on the solution quality.
In Fig. 14 the relative error of the vertical displacement is plotted over the element size e. Compared with the pure tension test, in which the results of sys. a and sys. c were the same, here the solution of sys. a no longer converges against the analytical solution. The reason is that the support of sys. a introduces a singularity. As a result, the solution starts to converge against the analytical solution at first, but ends up overestimating it as soon as the mesh is too fine. Only sys. c converges against the analytical solution. This demonstrates that neither sys. a nor sys. b can represent the boundary conditions according to the beam theory. Therefore, in the following, only sys. c is considered.
Next, the convergence behavior due to mesh refinement in the solid domain-and therefore a refinement over the height of the cross-section-is investigated. The results are presented in Fig. 15. It can be observed that the convergence behavior is of an order of 4 over the elements' edge length. This behavior was expected as the displacements are approximated with quadratic shape functions. Furthermore, using a relatively long transition element with = 0.1, the error stays above 10 −7 . Lowering the transition element length leads to lower overall errors in the solution when using fine meshes in the solid domain.
To investigate the influence of the transition element length on the solution, Fig. 16 depicts the resulting error due to varying the length with different solid element sizes. As is clearly visible, the solution converges when the transition The optimal convergence behavior is achieved when = e is used. In this case, the rate of convergence is approximately 3, or more precisely 3.35. Due to the linear interpolation over the transition element length , a convergence order of 4 cannot be achieved. As a result, the length of the transition element can be considered a mesh parameter and should be set as = e.
Thus far, we have fixed the geometric and material parameters to demonstrate that the results converge against the analytical solution. In Fig. 17, the impact on the solution is investigated when different Young's moduli E are chosen. The analysis is conducted for different element sizes e and a transition element length of = e. When the element size is reduced from e = 0.1 to e = 0.02, the displacement error of course is reduced. From the diagram, it is more crucial to note that material properties do not affect the quality of the solution. Only for the finest mesh with e = 0.02, slightly different results are observed as numerical accuracy is reached.

Shear and bending test, clamped -Load F y andˇ= 0 on the right side
In this example, the focus is on the mixed loading case, which consists of a load perpendicular to the beam and a moment due to the constraint β = 0 on the right side. The left side of the beam remains fully clamped. Only sys. c is investigated.
A mesh convergence study is performed and the results are presented in Fig. 18. The study is performed for transition element lengths = 0.1, 0.01, ,0.001. The results are comparable to those in Fig. 15. This indicates that mixed  As presented in Fig. 19, a convergence study for the transition element length is performed with different 2D solid element sizes e. The results converge as expected with a decreasing transition element length . As the interpolation order in the length direction of the transition element is linear, a convergence order of 4 cannot be achieved. The highest convergence order is reached when using = e. Again, this demonstrates that the transition element length can be interpreted a mesh parameter and should be set as = e. Figure 20 depicts the resulting shear stresses on a deformed mesh for a short beam with L/h = 2. As the transition element is applied to both sides of the structure, the resulting shear stresses exhibit no distortions and are in perfect agreement with Timoshenko's beam assumptions.
The resulting shear stresses over height are depicted in Fig.  21. The cut is performed in the 2D solid domain at the right solid face S . Normal stresses σ x are linearly distributed over the height, shears stresses σ xy are distributed with a quadratic parabola, and normal thickness stresses σ y are zero. Figure 22 depicts the relative error of the vertical tipdisplacement over the beam length L. The ratio L/h is varied from 0.5 to 10. Even for very thick beams, the error compared with Timoshenko's beam theory is very small and becomes even smaller for thinner beams. This is expected as for thick beams shear deformation becomes increasingly important (see Sect. 5.7).

Simple frame
In this example, a simple frame structure is considered. To demonstrate the performance of the proposed transition element, a full model consisting of 2D solid elements is compared with a mixed model consisting of 2D solid and beam elements as well as a pure beam model. In case of the mixed model, only areas in which beam theory is no longer valid are represented by 2D solid elements.
The geometry of the system with dimensions and loads is depicted in Fig. 23 for the full and mixed models. Beams are represented by lines ending at the distance (transition element length) of the 2D solid parts. The thickness perpendicular to the drawing is set as b = 1. A homogeneous, isotropic material with Young's modulus E = 100 and Poisson ratio ν = 0.3 is chosen. For the beam elements based on Timoshenko's theory, a shear correction factor of κ = 5/6 is employed. The loads on the structure are q x = 1 and q y = 1 as specified in Fig. 23. The 2D solid part is modeled with perfect square elements. In the case of the pure beam model, the 2D solid parts are removed and only the lengths of already presented beams are adjusted to connect the beams and reach the supports.
A mesh convergence study is depicted in Fig. 24. The full model (full) and mixed model (mixed) converge against the same solution. Using the mixed model, the same solution  quality can be achieved with approximately half the number of degrees of freedom compared with the full model. By contrast, the pure beam model exhibits behavior that is overly soft. Its corner displacement is computed by adding the vertical distance to the corner multiplied by the rotations to the corner displacement of the beam model. Even if the rotational part is neglected, the displacement value evaluated by the beam model is approximately 3% too high. The reason for this deviation in the results of the beam model is that the beam theory is neither valid in the corners nor at the supports compared with the following stress contour plots in Fig. 25.
In Fig. 25, the resulting stress contour plots in the mixed model are presented on the deformed structure. To compare the deformation state, the deformed mesh of the full model is also depicted. The deformation behavior of the full and mixed models is the same. Furthermore, no spurious stresses due to the transition element can be observed. The boundary conditions are applied using the proposed transition element.
In Fig. 26, the stresses over the thickness at the cut A-A are depicted (see Fig. 23). Even though the cut is performed directly at the transition element, the results of the mixed model agree with those of the 2D solid model.

Arc
In this example, a circular arc is investigated. The aim is to investigate the impact of the proposed transition element when applied to curved structures. It is important to note that in this case, the assumption made in Sect. 2 of the faces S and B being initially parallel is violated.
The structure's geometry and load are illustrated in Fig.  27. It is loaded by a single load F in the horizontal direction at the outer position of the arc under 45 • to the negative xaxis. To compare the solution, the displacement u of a node in the x-direction at the outer position under 45 • to the xaxis is observed. The radius is r = 100 and the thickness is t = 5. The material properties are a Young's modulus E = 100 and a Poisson's ratio ν = 0.3. Boundary conditions are applied using the proposed transition element, Fig. 25 Simple frame -Stress contour plots (σ x , σ y , and σ xy ) of the coupled model with the proposed transition element on a deformed structure underneath the deformed mesh of the full model. The transition element length for these plots is set as = 0.01 such that both ends of the arc are clamped. To determine the impact on solution quality, the arc is subdivided into several smaller sections. Each segment-segment connection is modeled with two transition elements. The parameter s specifies the number of sections that the model is built of. In Fig. 27, the model is depicted for s = 5 subdivisions. The transition element length is set as = 0.01.
In Fig. 28, the relative error between the full and subdivided models for different numbers of subdivisions is presented. To compute the relative error, the displacements u S of the subdivided model and u f of the full model are calculated using regular meshes with the same element size e. Thus, it can be expected that the error due to the discretization with 2D solid elements is the same in both solutions and only the impact of the transition element is extracted. Even though the curvature of the structure is not considered in the transition element formulation, the resulting deviation from the full model is very small but increases when more subdivisions are used. Only for very coarse meshes can a convergence of the coupled model toward the full model be seen. For finer meshes, the deviation of the solution of the subdivided model from that of the full model converges against a certain value depending on the number of subdivisions.
In Fig. 29, stress contour plots on a deformed mesh are presented. The thickness of the depicted structure is t = 5 with a total of 18 sub-structures connected by the proposed transition element. The transition element nodes are depicted with black dots. In the plots, no deterioration of the stresses can be observed due to the transition element.
Next, the same investigation is performed for a thicker arc. The results in Fig. 30 indicate that the relative error in displacements compared with a pure solid model became higher, even though it is still small with a magnitude of 10 −5 -10 −3 .
A detailed insight into the stress state at the transition zone is illustrated in Fig. 31. The top row presents the subdivided structure while the bottom row presents the results of the full model. The transition element can be identified by the vertical black line with the dot. While the normal stresses σ x and shear stresses σ xy are in very good agreement, the stresses in the thickness direction σ y (depicted in the center) exhibit some distortion in the transition zone.
These spurious stresses are two to three orders of magnitude smaller compared with the stresses illustrated in Fig. 32. Nevertheless, it is important to note that two solid surfaces with 64 quadratic elements over the thickness are coupled by only three degrees of freedom in this example or, in other words, 258 degrees of freedom on the left side are coupled with 258 degrees of freedom on the right side by only two translations and one rotation. Furthermore, the transition element itself is assumed to be straight and does not consider the curvature of the structure. In Fig. 33, the impact of the ratio r /t is depicted. For large ratios r /t, the relative deviation of the resulting displacement of the subdivided model from the full model is very small with a value of 10 −7 . With a decreasing ratio (increasing thickness t), the relative deviation increases. Furthermore, this example indicates that with a ratio r /t ≥ 10, the proposed element can be used in curved areas despite the geometrical assumption of initially parallel faces S and B .

Conclusions
The element formulation proposed in this paper allows for the coupling of 2D solid elements with beam elements. The numerical investigations demonstrated that there were no spurious stresses on the solid part even when the length of the transition zone was close to zero. The resulting stress distributions in a section of the solid, located directly in the transition zone, were in very good agreement with the analytical solution. An interpretation of the proposed element as a solid element is possible. Therefore, the length can be interpreted as a mesh parameter with converging results for → 0. The element contains general deformations of the cross-section without the necessity of a precomputed warping function. Even though the element assumes a straight transition zone, with initially plane and parallel faces at the solid and beam parts, applying it to moderate curved, slender structures is possible. The element possesses only one rotational degree of freedom and displacement degrees of freedom. The resulting stiffness matrix is symmetric and, most crucially, it leads to stress distributions on the solid face of the transition element that are in accordance with the results of a full 2D solid model. This makes it perfectly applicable in situations where it is critical that the stress distribution on the solid surface is not affected by spurious stresses. A further application field of the presented approach is homogenization, in which the beam kinematics are coupled to a representative volume element.  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.