Insight into numerical characteristics of embedded finite elements for pile‐type structures employing an enhanced formulation

A number of problems encountered in the field of computational geotechnics requires the modelling of numerous pile‐type structures embedded in a three‐dimensional soil continuum. Traditional modelling approaches to this problem result in computational costs that must be regarded as unbearable for most practical purposes. As an attractive alternative, we propose an enhanced embedded FE model with implicit interaction surface (EB‐I) where coupling between the contacting domains is realized by an implicit surface‐to‐volume (2D‐to‐3D) coupling scheme. As a novelty, the latter implements a non‐linear interface constitutive model that allows for explicit consideration of endpoint interaction, but does not represent a constraint for the solid mesh generation. As the slender structure is discretized employing the Timoshenko beam theory, shear deformability is explicitly considered, as opposed to earlier EB‐I‐type models reassessed in this paper. The credibility of the proposed EB‐I is numerically validated on the basis of comparative studies. It is found that the 2D‐to‐3D coupling scheme generally improves the well‐posedness of the resultant global stiffness matrix, making the proposed EB‐I computationally competitive to geometrically simplified line‐to‐volume coupling schemes. Future lines of research are carefully addressed throughout this work and include the normal stress recovery technique as well as applications to large‐scale simulations.


INTRODUCTION
Embedding slender structures into a matrix material represents a versatile strategy tailored for the development of sustainable and cost-efficient solutions to numerous technical problems. 1In the context of computational modelling, as considered in the following, the former will be referred to as sub-scale elements (SSE), the matrix material as macroscale element (MSE) and the composite material as embedded reinforcement model (ERM).][8][9][10] Likewise, this concept can be found in other fields than solid mechanics 11,12 ; for example, countless biological systems comprise biopolymer networks of highly slender filaments, whereas the latter govern crucial processes, including cell migration and cell division. 13In a very recent publication, 14 focus is placed on the electro-chemical analysis of fibrous electrodes where discrete fibres are embedded in an electrolyte matrix, a multi-functional material that can be used for energy storage and harvesting.6][17] In this view, pile-soil systems may be regarded as ERM.At variance to analytical and empirical methods traditionally utilized to design ERM involving pile-type structures, such as the p-y-method originally proposed by Reese et al. 18 or the strain wedge method first applied in the work of Norris 19 that are based upon more or less significant model assumptions, the finite element method (FEM) allows for the general and simultaneous consideration of complex geometries, material non-linearities, soil-structure interaction (SSI) phenomena, non-trivial loading situations and inhomogeneous geological environments.Therefore, the FEM possesses a relatively higher potential to provide a time-and cost-efficient design of related problems. 20][23] Considering the examples presented in the previous paragraphs, it appears almost logical that the evolution of suitable finite element modelling techniques has attracted widespread interest from the scientific community.The vast majority of associated developments has been introduced on a case-by-case basis resulting in cross-disciplinary lines of research, some of which overlap.In a more general sense, Goudarzi and Simone 24 follow a simplistic classification scheme into two broad modelling categories, namely explicit and implicit approaches.In the former case, frequently referred to as meanfield 25 or homogenization approach, 26,27 the collective effect of randomly distributed SSE on the global ERM response is taken into account in an indirect manner by means of anisotropic material models (Figure 1A).This entails the assumption that the analysis domain is macroscopically homogeneous and uniform in properties, albeit being composed of disjointed sub-regions. 28From a computational perspective, implicit modelling strategies, such as deployed in Kattis et al., 29 are particularly appealing, since they are applicable to a variety of simulation programs offering the opportunity to apply anisotropic material laws.Essentially, this allows to restrain the computational effort to stay beyond practical limits, as no additional degrees of freedom (DOF) are in demand to resolve the SSE. 30However, implicit modelling approaches require sub-scale information, presumably provided by ERM enclosing numerous resolved SSE geometries, to determine material properties.Furthermore, they are not capable of providing detailed insight into local effects, 31 and inconvenient for optimization purposes and the assessment of failure mechanisms, a crucial output of finite element analyses (FEA). 32onsequently, implicit modelling approaches should be restricted to global response analyses.
The second modelling approach explicitly resolves SSE either employing full three-dimensional (3D) domains comprising solid finite elements (FE) or dimensionally reduced discretizations, also known as discrete FE models (D-SSE) 33 and embedded FE models (EB-SSE), 34 respectively (Figure 1B,C).6][37] With regard to the modelling of piles through D-SSE, solid FE corresponding to the pile domain may be additionally circumscribed by contact formulations, such as isoparametric contact elements [38][39][40][41][42] in order to account for SSI.In the context of computational geotechnics, this expanded D-SSE modelling approach is often termed standard FE approach (SFEA). 43onventional D-SSE are subjected to mesh conformity constraints, that is, nodes coupled at the contact interface assume the same coordinates.Recently, FEM-based strategies that allow this restriction between contacting SSE-MSE domains to be lifted have emerged; prominent examples encompass the mortar FEM, 44,45 immersed FEM 46 and extended FEM, 47 some of which have been applied to geotechnical problems. 48A deep theoretical discussion is beyond the scope of the present paper, but it is important to point out that these approaches invoke numerically expensive cutting procedures to implicitly model interfaces, and constitute no general remedy against meshing complexities inherent to D-SSE. 49lthough D-SSE allow for a detailed description of the complete ERM geometry, and may be equipped with advanced constitutive models to describe inherent sources of non-linearity, they necessitate high levels of mesh refinement close to SSE-MSE contacts in order to capture the sought field variables with sufficient accuracy.This inevitably induces a substantial increase of computational burden. 50Considering practical computational limits in terms of memory capacity and total runtime duration, D-SSE are, therefore, not generally favoured and rarely applied to study the behaviour of largescale ERM involving a high number of SSE. 51In the realm of computational geotechnics, this holds particularly true for dynamic SSI problems. 52dopting simplifying geometrical and physical assumptions, 17,49 EB-SSE have evolved as attractive alternatives to overcome aforementioned obstacles as they substantially alleviate computational costs compared to D-SSE. 53In the range of our modelling assumptions, numerous numerical examples [54][55][56][57] demonstrate that related formulations yield very accurate results, are well-posed and exhibit optimal spatial convergence; therefore, EB-SSE are increasingly used to examine SSI phenomena in a general manner for different kinds of geotechnical problems. 58On the one hand, relative merits compared to D-SSE are attributed to the elimination of the conformal mesh generation constraint, 24 on the other hand, EB-SSE reduce the number of unknowns to be solved substantially. 59rom a numerical perspective, the development of EB-SSE is non-trivial to handle as the description of SSE-to-MSE interaction boils down to a mixed-dimensional coupling scheme linking dimensionally reduced 1D structural FE to 3D solid FE.Since the different types of involved FE are endowed with different discretizations, mixed-dimensional coupling schemes of EB-SSE need to keep track of a set of variables which reside on the boundary between domains of different dimensions; to this end, they require the simultaneous consideration of equations in different dimensions which adds complexity to their development. 60Although EB-SSE have gained increasing popularity within the field of computational geotechnics, commercially available implementations continue to suffer from numerical pitfalls, some of which are addressed in a very recent contribution. 61nspired by recent developments in the field of computational geotechnics, 62,63 the present paper introduces the generic framework of an enhanced EB-SSE with implicit interaction surface.The proposed EB-SSE is general in a sense that it is able to reproduce key features of the SFEA, such as plastic frictional tangential slip along the pile shaft 42 as well as a mechanically consistent redistribution of coupling forces at the base to realistically capture the response of pile-type structures.Moreover, its variational framework explicitly accounts for shear deformation occurring inside the SSE, making it more suitable for idealizing the behaviour of short piles in stiff soils. 16,64he remainder of this contribution is organized as follows: Section 2 reconsiders important landmarks in the development of EB-SSE based on a thorough literature review.To the authors' knowledge, this is the first attempt to categorize the latter in a rigorous manner with explicit crosslinks to the field of computational geotechnics.In Section 3, we state the fundamental modelling assumptions and implementation details of an enhanced EB-SSE, which builds on the existing EB-SSE 51 currently implemented in the FE code PLAXIS 3D. 65Section 4 is devoted to numerical experiments where we carefully study its numerical performance with particular focus on the coupling scheme and computational performance; we further shed light on the importance of a stress dependent formulation of the embedded interface constitutive model.Section 5 closes with the conclusions of this work and highlights further research aspects.

EMBEDDED FE MODELS FOR MULTI-POINT COUPLING OF 1D STRUCTURES WITH 3D CONTINUUM DOMAIN: AN OVERVIEW
From a numerical point of view, EB-SSE considered in this section are composed of 1D structural FE constituting the SSE, such as truss, bar or beam FE, and mixed-dimensional coupling schemes with high-dimensional gap, 60 whereas the dimensionality refers to the number of dimensions used to parametrize the coupled FE types. 66In view of the envisaged ERM application class, this concerns pile-type structures embedded in a background solid mesh representing the MSE.
Principally, the development of numerical coupling schemes with high-dimensional gap is more challenging compared to those designed for low-dimensional gap models, as in the latter case, the mutual interaction between contacting bodies defined across their contact surfaces is better understood in both mathematical and physical terms. 42,60Typically, there is no explicit geometry object in the 3D solid mesh to match the centerline of dimensionally reduced 1D structural FE, which renders high-dimensional gap problems naturally more prone to numerical stability issues. 67Furthermore, EB-SSE entail several modelling errors, including the artificial continuity assumption 61 imposed on the solid domain (Ω  ) and the spurious superposition of material points that share the same spatial position. 68,69The latter aspect constitutes a general limitation inherent to the family of embedded overlay elements 70,71 in which SSE have to be understood as being superimposed on the MSE domain rather than being surrounded by it.
Over a period of more than five decades, numerous researchers have successfully contributed to the gradual evolution of EB-SSE, with some enhancements having also found application in the realm of computational geotechnics.In spite of remaining numerical pitfalls reported in the previous paragraph, these have paved the way to the ubiquitous deployment of EB-SSE in geotechnical FEA.As a consequence, the analysis aim has been increasingly shifted from global response analysis to the investigation of complex boundary value problems (BVP) that provide local insight into the mechanical behaviour at critical locations; for example, see Refs. 56,57,72The purpose of this section is to organize relevant advances in the development of EB-SSE in a comprehensible manner, many of which have been integrated in the enhanced variational EB-SSE framework discussed in Section 3. In addition, the authors firmly believe that the understanding of the theoretical concepts is crucial to assess the scope of validity for geotechnical FEA involving EB-SSE, and exploit the capabilities of state-of-the-art EB-SSE in future research works.

Embedded FE models with implicit interaction line
The pre-dominant share of EB-SSE documented in the literature is equipped with multi-point line-to-volume (1D-to-3D) coupling schemes that enforce constraint equations among the relevant DOF.This means that the mutual interaction is mathematically described at the centerline position of the 1D structural FE.Early work on so-called embedded FE models with implicit interaction line (EB-L) has been carried out in the context of reinforced concrete, where they have been found eligible for the modelling of reinforcing agents, 73,74 with the restriction that the reinforcements have to align with the parametric coordinate system of the solid discretization.In particularly relevant publications, 2,75,76 this approach has been extended to allow for an arbitrary orientation of straight 1D structural FE relative to the solid mesh, a key feature of modern EB-SSE.Further work documented in Elwi and Hrudey 77 has pursued the generalization of EB-L to curved 1D geometries.
In view of the suitability for geotechnical FEA, however, it must be noted that the basic idea of those mentioned previous works is to literally add stiffness contributions stemming from an EB-L to the global stiffness matrix of the ERM, that is, no DOF are explicitly solved for the EB-SSE.As a consequence, this EB-L class is subject to the same obstacles as it is the case with implicit modelling approaches.With reference to Kang et al., 78 we will denote this class as semi-discrete EB-L (Figure 2A).More specifically, semi-discrete EB-L require perfect adherence conditions to be imposed on all coupled nodal pairs, and imply that the kinematics of 1D structural FE can be appropriately described by the solid displacement field, commonly expressed by means of Boltzmann continua 31 incorporating standard C0-continuous element functions 79 ; given the C1-continuous centerline representation of conventional 1D structural FE, 80 however, this assumption is not generally warranted.As semi-discrete EB-L demand the implementation of a relatively expensive post-processing procedure to obtain stress resultants, 62 they must also be regarded as prohibitive for numerous structural design tasks.To some extent, fully discrete 78 EB-L resolve aforementioned issues by introducing independent translational DOF at discrete coupling points (Figure 2B), resulting in the addition of unknowns to the discrete system of equations.Adopting the classification scheme for the treatment of contact constraints after Wriggers, 42 variationally consistent methods to enforce the coupling constraints of fully discrete EB-L include, but are not restricted to, the penalty method, 20 Lagrange multiplier method, 81 mortar method, 49 Nitsche method 82 and the formulation of constitutive equations at the interaction line. 51Similar to semi-discrete EB-L, fully discrete EB-L can be superimposed on Ω  irrespective of the solid mesh. 76rom a geotechnical perspective, the development of fully discrete EB-L, in lieu of semi-discrete EB-L, has been motivated by an increasing demand for modelling techniques that allow for a reliable prediction of local field variables (for example, stresses and displacements) in the vicinity of pile-type structures, account for non-linear SSI effects, and offer the advantage of providing direct insight into stress resultants as crucial ingredients to the structural design. 58In the realm of computational geotechnics, discrete EB-L can, therefore, be regarded as state-of-the-art, and are generally preferred by modern geotechnical FE codes 58 ; in turn, semi-discrete EB-L are favoured in FEA concerning the global system response of composite materials comprising a high number of randomly dispersed glued or moulded SSE at microscopic length scales, such that all modes of relative motion are suppressed and full adherence conditions prevail. 28part from the distinction between semi-or fully discrete EB-L, a reasonable classification scheme designed for EB-SSE has to distinguish between the various types of 1D structural FE that can be used to discretize the SSE.In the case of reinforced composite materials, 2,3,10 bar elements may be well suited as the internal elastic energy contribution mainly infers from axial compression or tension 83 ; we refer to members of this EB-SSE class as bar EB-L.While convenient for the modelling of geotechnical problems dominated by axial deformation, such as micropiles subjected to axial tension, 84 bar EB-L are prohibitive for application scenarios with significant bending deformation.Beam EB-L, 51,85,86 in contrast, are conceptually different to embedded bar EB-L proposed in the 1990s for the modelling of steel cords in large strain FEA of rubber composites 71 or reinforcing agents in reinforced concrete structures 87 ; the salient difference concerns the bending part of the elastic strain energy which is explicitly regarded in the internal virtual work equations of beam EB-L. 83his extension ensures a mechanically consistent consideration of the structural resistance stemming from the bending deformation mode; pile bending imposed by lateral ground pressure is a case in point where the bending component is expected to dominate the internal potential energy contribution of SSE. 15 In comparison with bar EB-L, beam EB-L can, therefore, be considered as more general version of EB-SSE.A detailed discussion on the role of beam theories adopted in the context of EB-L can be found in Steinbrecher et al. 49 From a continuum mechanics perspective, the definition of realistic contact conditions is crucial to capture the contact kinematics between two contacting bodies with high accuracy. 42,69This holds particularly true for geotechnical problems experiencing SSI, as typically encountered in the envisaged field of EB-SSE applications.In the course of EB-SSE evolution, a paradigm shift from simplified linear stick EB-L 73 to non-linear stick-slip EB-L described by a regularization of Coulomb's law 51,85 could be observed, whereas the latter has been found eligible for the analysis of pile-type structures subjected to axial loads. 51,58Research on the performance of EB-L under inclined and horizontal loading, potentially resulting in the formation of a gap between the contacting bodies, is much rarer.However, recent studies 68,88 indicate that the performance of EB-L under mixed-mode loading, as it may occur in the case for horizontally loaded piles, is yet sub-optimal.
Recent work concerning the generalization of EB-L has followed two lines of research, namely, (i) the introduction of weak discontinuity enrichment functions at solid element nodes crossed by 1D structural FE using the partition of unity enrichment method. 89The so-called strain discontinuity enriched embedded FE model, 24 herein referred to as discontinuous EB-L, mitigates some shortcomings of traditional continuous EB-L.This includes the artificial continuity constraint 61 inside Ω  and numerical oscillations. 58The authors, however, state that discontinuous EB-L are not yet competitive as they require the generation of a background integration mesh, at present a computationally expensive process that impairs their deployment in traditional engineering problems. 24ii) the supplementary enforcement of rotational coupling constraints between 1D structural FE and Ω  , that is, full coupling EB-L (Figure 2C).31,90 Unlike translational coupling EB-L, full coupling EB-L allow not only for the coupling of translational DOF between the contacting bodies to describe the cross-section positions of 1D structural FE, but also to enforce rotational constraints at mutual contacts in order to couple cross-section orientations to the solid.In this way, EB-L gain the ability to transfer torsional loads in a mechanically consistent manner; for detailed information regarding application scenarios where full coupling EB-L are generally preferred one is referred to Steinbrecher et al. 31 Although full coupling EB-L appear promising, it should be noted that their validation has so far been restricted to linear material models and conceptually simplified test scenarios; hence, future work is in demand to prove their suitability to geotechnical BVP considering more realistic contact and material behaviour, that is, inherent sources of non-linearity.
In principal, the abovementioned EB-L configurations involve truly dimensionally reduced 1D-to-3D coupling strategies formulated in a weak variational sense on the basis of the virtual work principle, 91 in which the coupling terms are described by balance, compatibility and constitutive equations that are combined into one partial differential equation for displacement 1,92 ; from a numerical point of view, related problems consist in finding the equilibrium state within a discrete Ω  that is subject to distinct interaction load components, namely, inner line loads acting at the centerline position and inner point loads imposed at the proximal end. 61Adopting the nomenclature of Podio-Guidugli and Favata, 92 these load components give rise to the juxtaposition of problems similar to the Kelvin problem.As a consequence, the exact solution to be approximated contains singularities in the sought field quantities close to the positions of interaction, possibly yielding a badly conditioned problem that disqualifies the use of EB-L.Despite the existence of well-established solution strategies to overcome deduced numerical obstacles, 43,86 the 1D-to-3D coupling scheme impairs the general applicability of EB-L in geotechnical FEA. 61

Embedded FE models with implicit interaction surface
A more natural choice for mathematically describing the coupling formalism is to expand the interaction geometry from an implicit interaction line to an implicit interaction surface.This concept is adopted by embedded FE models with implicit interaction surface (EB-I), as conceptually illustrated in Figure 2D.In the context of computational geotechnics, fundamental research on EB-I can be found in Refs. 62,63,79,93,94Since there is no explicit interaction surface in Ω  to impose the coupling constraints on, these models incorporate frame-invariant mapping schemes that map the relevant DOF of the contacting bodies to a well-defined interaction surface, over which coupling constraints are imposed in a weak sense.
In contrast to 1D-to-3D coupling schemes associated with EB-L, this strategy induces a surface-to-volume (2D-to-3D) Generic overview of sub-scale element modelling approaches with precise categorization of embedded FE models.
coupling problem where an implicit interaction surface, defined by the geometrical properties of 1D structural FE, is embedded into the background solid domain; in this way, 1D structural FE interact with all solid FE along their surface, instead of one central solid FE located at the centerline position only.It should be pointed out that 1D-to-3D versus 2D-to-3D coupling terms have a different physical dimensionality, and therefore different units; while the former are associated with a line load, the latter refer to a surface load. 49Increasing the physical dimensionality of the coupling scheme allows the high-slenderness restriction of EB-SSE 49 to be lifted.Conceptually, EB-I can, therefore, be considered suitable for analysing geotechnical structures other than micropiles, as initially intended with respect to EB-L. 85ecently published comparative studies demonstrate principal merits of EB-I compared to EB-L, some of which are considered relevant for geotechnical FEA. 58,61,62,68,88,95A considerably important advantage is that EB-I predict the near-field displacements close to the physical contact interaction with significantly higher accuracy (Figure 2E); EB-I can, therefore, be regarded as generally superior in terms of calculation fidelity.This observation, together with EB-L related considerations stressed in Section 2.1, are put into context in the form of an EB-SSE classification scheme that assesses the individual numerical performance in a qualitative manner based on findings reported in the literature (Figure 3).To the authors' best knowledge, the latter represents the first attempt to structure a dynamically growing number of developments related to EB-SSE in the context of computational geotechnics, thereby adopting mechanically justified categories combined with a semantically consistent terminology.It is hoped that this endeavour sheds light on the fundamental principles of EB-SSE, and paves the way to generally applicable formulations that reduce the computational demands of D-SSE while retaining their solution fidelity.
Above considerations have motivated the development of an enhanced EB-I that explicitly accounts for non-linear interface behaviour along implicit interaction surfaces covering both, the shaft and the base (combined coupling).The reader should notice that this feature represents an enhancement compared to the EB-I versions highlighted in Figure 3 where coupling is restricted to perimetral positions along the shaft (shaft coupling); since the proposed EB-I is additionally equipped with a Coulomb-type failure criterion, it can be classified as stick-slip EB-I.In the next section, insight into relevant implementational details is provided.

EMBEDDED BEAM WITH INTERACTION SURFACE: MODEL DEVELOPMENT
Extending the mechanical framework of the EB-I models detailed in Refs., 62,63 we present an enhanced EB-I formulation for embedding three-noded shear flexible Timoshenko beams 96 into a 3D solid domain, that is, shear deformability is explicitly taken into account. 83Contrary to previous publications concerning the capabilities of earlier versions of the proposed EB-I model in a practical context, 58,61,68,84,88 the present paper provides detailed information about the underlying implementational framework developed in the FE code PLAXIS 3D as well.Without loss of generality, coupling between the contacting domains is described at implicit 61 interaction surfaces Γ  where displacement compatibility and stress equilibrium conditions are satisfied.Essentially, the EB-I interface behaviour is numerically described not only by means of contact conditions established along the shaft, but also across the base surface.In addition, a non-linear interface constitutive model is introduced which is able to account for the soil stress state in the vicinity of Γ  .The significance of these features in geotechnical FEA will be demonstrated in Section 4.2.
Adopting continuum mechanics fundamentals, non-contacting boundaries of the beam domain Ω  and Ω  are decomposed into Dirichlet and Neumann boundaries, Γ  and Γ  , which satisfy where Γ is the complete set of boundaries over which boundary conditions are defined (Figure 4).Subsequently, the mathematical description of the proposed EB-I model is provided in symbolic notation where scalars are denoted by means of Greek and Italic letters, for example,  and R. Lowercase bold letters are used to denote vectors (e.g., ũ ).Likewise, uppercase bold letters represent tensors of higher order (e.g.,   ).Components of vectors and higherorder tensors are defined considering a right-handed Cartesian basis, with a fixed set of orthonormal base unit vectors.Unless stated otherwise, the normal pressure is assumed positive in compression.

Preliminaries and kinematical description
Figure 4 illustrates the basic layout of the proposed EB-I model, which comprises the following ingredients and assumptions:  error is proportional to the beam volume fraction and the relative stiffness ratio between beam and solid. 49In view of the envisaged application classes, the beam stiffness is considerably higher compared to the solid stiffness; hence, the influence of the modelling error can be considered insignificant.• At present, the proposed EB-I model is unable to keep track of gap formation mechanisms as a consequence of lateral loading; it should therefore be applied to geotechnical problems where the beam is assumed to remain in contact with the surrounding soil throughout the analyses.Future EB-I model extensions may address its generalization to more complex loading situations where gapping phenomena are likely to occur.• The variational fundamentals concerning coupling constraints between Ω  and Ω  read: which are fulfilled in the reference configuration, that is, the coupling operator is constant and requires no evaluation of the normal vector  on Γ  in the current configuration.In Equation (3), t denotes a motion parameter.It should be noted that this formulation of the coupling conditions is only valid for contact problems which do not include frictional sliding, as the friction process is dissipative, that is, the solution becomes path-dependent.As will be shown in Section 3.4, however, this hypothesis represents no limitation of the EB-I model, which is completely general from a constitutive point of view.This includes, but is not restricted to the constitutive description of the non-linear interface behaviour.

Variational fundamentals
Under the assumption of small deformations, contributions to the total potential energy of the ERM can be split into beam Π  (•), solid Π  (•) and coupling terms Π  (•), whereas the beam and solid terms are independent of the coupling constraints, that is, traditional discretization and modelling techniques can be used for Ω  and Ω  .Applying the principle of virtual work, the Lagrangian functional (•) reads where   ,   ∈ ℝ 3 denote beam and solid displacement vector fields, respectively.Equilibrium of the mechanical ERM system is achieved by minimizing (•), 97 hence we search for the solution to the problem: where   ×   denotes the admissible solution space.With regard to the first and the second term in Equation ( 4) representing the virtual work equations of the beam and Ω  , well-established discretization and modelling techniques can be adopted without modifications; since this is not the main aspect of the present work, the interested reader is referred to Refs. 83,91The last term in Equation ( 4) constitutes the potential energy related to the kinematic coupling constraints in Equation ( 3), which are enforced at x : Herein,   ∈ ℝ + represents a scalar-valued penalty parameter used to penalize any violation of the kinematic coupling constraints.Adopting the penalty method, it can be shown that Equation ( 3) is exactly enforced for   → ∞; compare Wriggers. 42The stationary point of the potential energy function formulated in Equation ( 4) is found by setting its variation equal to zero: In Equation ( 7),   (ũ  − ũ ) can be interpreted as discrete coupling forces mobilized at x .In order to mathematically explore the load transfer mechanism of the EB-I model in more detail, we can re-formulate the coupling conditions defined in Equation (3): x (, , ) =   () +  ⋅ cos () ⋅   1 () +  ⋅ sin () ⋅   2 () ⏟⎴⎴⎴⎴⎴⎴⎴⎴⎴⎴⎴⎴⏟⎴⎴⎴⎴⎴⎴⎴⎴⎴⎴⎴⎴⏟   (,,) (9)   where x , x ∈ ℝ 3 denote the current positions of coupled solid and beam material points initially sharing the same location at Γ  in the reference configuration (Figure 5B),   () ∈ ℝ 3 defines the beam centerline curve connection cross-section centroids (Figure 6A), and   (, , ) ∈ ℝ 3 represents cross-section position vectors that point from the cross-section centroid to Γ  , defined by means of the polar coordinates  ∈ ℝ + and  ∈ [0, 2); compare Figure 4.The rate form of Equation ( 8) in the reference configuration is defined by the material time derivative: which may be interpreted as enforcing zero relative velocities between the contacting bodies at x .By inserting Equations ( 8) and ( 9) into Equation ( 6), an equivalent relationship is obtained for the penalty potential: Eventually, variation of Equation ( 11) allows for a direct inspection of the load transfer mechanism invoked by the proposed EB-I model: Therein,    signifies an infinitesimal rotation variation, Γ , represents the discrete interaction surface segment (Figure 4) belonging to x , and  ∈  3 denotes a skew-symmetric tensor, with  3 representing a set of skew-symmetric This reformulation in Equation ( 12) leads to a direct expression of coupling forces  ,x  acting on both, the beam centerline and the solid.In addition, Equation ( 12) pinpoints the coupling moment  , acting on the beam crosssection.Obviously, this demonstrates the projection of purely translational coupling constraints enforced on Γ  onto the beam centerline, and highlights arising beam rotational coupling terms.Interested readers may refer to the work of Granitzer and Tschuchnigg 61 which highlights engineering problems in which this aspect may govern the structural response considerably.

Finite element approximation
The numerical treatment of coupling constraints is crucial in the development of EB-SSE since they may add severe nonlinearity to the problem at hand.This includes the constitutive description of the interface behaviour, the numerical enforcement of coupling constraints as well as the discretization of Γ  .With regard to the proposed EB-I model, coupling is incorporated by means of kinematical constraints established over the beam and solid displacement field as previously stated in Equation (3).More specifically, at each beam node   ( ξ) along the beam centerline, multiple equally spaced coupling points x are defined along Γ , as well as at Γ , covering the virtual beam geometry.Therefore, all displacements and tractions describing the coupling behaviour are evaluated at the barycenter of individual Γ , (Figure 6A).From a mechanical point of view, each x is tied to the underlying Ω  through linear penalty constraints.Coupling contributions are numerically integrated, stipulating that where   ,   and   denote the number of x specified over the perimeter, the length and the base of respective beam elements (Figure 6B).In this context, it should be noted that an Γ , (x  ) does not follow from a direct discretization of Γ  , instead it is analytically described via the beam cross-section orientation, the beam cross-section geometry and the discretized beam centerline.As a consequence, the proposed EB-I does not require an additional surface mesh to be integrated in Ω  , hence it maintains the mesh-independency feature of EB-L models, as discussed in Section 2.1.
Adopting the Galerkin method, 91 we can mathematically express the beam constituent in Equation ( 3) on the basis of Lagrangian element functions    occupied by the ith node of consecutive three-noded beam FE, and related translational and rotational DOF û , θ ∈ ℝ 9 : Therein,   , denotes a scalar projection of the Euclidian base unit vector    in the direction of the Euclidian base unit vector   ; see Section 3.1 for a specification of reference systems.Assembling û and θ in a single vector ̂û  ∈ ℝ 18 , we can re-write the discrete beam displacement field at x in compact form introducing the mapping function matrix  ∈ ℝ 318 : In view of the solid constituent in Equation ( 3), the most natural choice for estimating the solid displacement field at Γ , is to interpolate within the corresponding solid FE: ũ,ℎ (x  ) = (x  ) ⋅ û (17)   where  ∈ ℝ 3× is the element function matrix of the -noded solid FE, and û ∈ ℝ  denotes the solid FE nodal displacement vector.Depending on the mesh topology, a solid FE may include zero, one or multiple x within its boundaries.For this purpose, a search algorithm is implemented which allows for the determination of corresponding pairs of target solid FE (TSFE) and x (Figure 6A) based on the inverse mapping of isoparametric solid element functions, similar to Ninić et al. 20 To limit the computational effort, this procedure is executed only once during the mesh discretization.
The coupling constraints stated in Equation ( 8) are locally enforced through a point collocation method, meaning that the differential equation in Equation ( 10) is explicitly solved for all x using the penalty method.In this way, coupling contributions to the global stiffness matrix as well as the global force vector are computed in a point-wise manner; see also Remark 3.1.Describing the variation of the total coupling potential in Equation ( 12) on the basis of discretized beam and solid displacement fields as defined in Equations ( 16) and ( 17), we can derive the local stiffness matrix K,x  and the resisting force vector f,x  with respect to x in single tensorial equations: whereas point-wise information on nodal quantities related to individual x is stored in vector form  x = [ ̂û  , û ]  .For practical reasons, ̂û  is explicitly solved, hence they are kept in the discrete weak form of the BVP; relative merits of this approach have been addressed in Section 2.1 in the context of semi-discrete ERM.The local contributions can readily be assembled to the global system of equations using standard procedures, hence will not be stated here for the sake of brevity.
Remark 3.1.The proposed coupling approach naturally captures the three-dimensionality of the problem, circumventing the need to incorporate additional kinematical hypotheses, such as incorporated on the basis of a projection operator in Turello et al. 63 Without loss of generality, this aspect reduces the implementational complexity of the proposed EB-I model.

Incorporation of non-linear interface constitutive model
Conceptually, the proposed EB-I is not necessarily treated as direct problem of compatibility on the basis of Equation ( 3), but can be extended by incorporating constitutive relationships that relate coupling forces to slip and compression.Inspired by the EB-L formulation 51 currently employed in the commercial software code PLAXIS 3D, the proposed EB-I model with weakly enforced coupling constraints is augmented with a non-linear interface constitutive in order to account for axial slipping and separation between Ω  and Ω  .This yields the advantage that the gradual mobilization of skin resistance can be tracked with higher accuracy, while allowing for reasonable estimates of the ultimate bearing capacity, 43 as will be shown in Section 4.2.Contrary to the abovementioned case, in which perfect adherence coupling constraints stated in Equation ( 3) are enforced along Γ  by adding the penalty term in Equation ( 12) to the weak form in Equation ( 7), the coupling term to be used here is governed by a constitutive relationship that is locally established at each x : The latter is equivalent to the virtual work produced by the traction vector t ∈ ℝ 3 (Figure 6B) and the virtual relative displacements ũrel , omitting virtual work contributions stemming from any external loads and body forces.Adopting the methodology first employed by Tschuchnigg and Schweiger, 51 an elasto-plastic relationship between ũrel and t in both contact tangential directions ẽ2 and ẽ3 is defined.This means that the hypothesis regarding the coupling constraints documented in Equation ( 3) is no longer enforced, that is, depending on the loading history relative displacements between Ω  and Ω  may contribute to the dissipation of internal energy induced by friction processes. 98he energy-conjugated components of ũrel = [ũ rel,1 , ũrel,2 , ũrel,3 ]  and t = [ t1 , t2 , t3 ]  defined by the local axes x 1 x2 x3 with the origin placed at x are related using 1D constitutive laws: One can easily see that the introduction of the above constitutive equations can be interpreted as an uncoupled nonlinear penalty formulation 42 that describes the constitutive behaviour during the pre-failure loading stage; compare Remark 3.2.Therein,   ,  cir and  ax represent the embedded interface stiffness in normal, circumferential and axial directions with respect to Γ , (x  ) specified over the perimeter of the virtual beam geometry,   denotes the embedded interface stiffness in both lateral directions as well as in axial direction with respect to Γ , (x  ) specified at the base (Figure 7A), and  ∈ ℝ 3 represents the identify matrix.  and   denote the diagonal matrices comprising the embedded interface stiffness parameters computed at the perimeter and the base, respectively, which are defined as (see also Remark 3.3) In Equations ( 22)-( 24),    is the average shear modulus at the stress points belonging to the TSFE corresponding to x , 51   is the interface Poisson's ratio, which has a default value of 0.45 43 and R is the virtual beam radius.Γ 1 , Γ 2 and Γ 3 as well as Δ 1 and Δ 2 are scalars assigned with default values of 1 and 0, respectively; these optional parameters are incorporated for scientific purposes, for example, to study the influence of the embedded interface stiffness magnitude on the well-posedness of the global stiffness matrix in order to ensure a high level of numerical robustness of the proposed EB-I model. 42mark 3.2.The proposed EB-I model adopts well-validated relations from Tschuchnigg 43 where    ,  and   have been found suitable for describing the non-linear interface behaviour.In case, the behaviour of Ω  is described using a constitutive model that accounts for the stress-dependent stiffness of the soil,  𝑎𝑣  is updated at the beginning of each load step based on the solid stress field state at the end of the previous load step.In this way, the non-linearity of the coupling behaviour is naturally included in the proposed EB-I model, omitting the tedious definition of fictitious constitutive parameters on a case-by-case basis followed by existing approaches for describing non-linear interface behaviour in combination with EB-I models; for example, see Turello et al. 63 Remark 3.3.In each time step [  ,  +1 ], a quasi-static implicit time integration is applied, in which the traction vector t,+1 ∈ ℝ 3 at x is calculated based on the elastic predictor-plastic corrector concept.Correspondingly, the traction vector update can be rewritten in incremental form: As can be inferred from Equations ( 25) and ( 26), the proposed EB-I performs an incremental traction update at each x since the coupling behaviour is described on the interaction surface.This differs from the underlying concept of the existing EB-L model 51 where the coupling behaviour is described on a line, resulting in an incremental update of the nodal force vectors Δ  residing on the jth beam node r, .However, with some abuse in notation, the latter can be related knowing the traction vectors Δ t at each time step residing on the same cross-section: Regarding the constitutive description, it should be noted that this requires dimensional changes of the proposed constituents in   and   , defined in Equations ( 22)-( 24), compared to the existing EB-L model, which are, therefore, defined as the original relations divided by the perimeter (2 ⋅  ⋅ ) and the base area ( 2 ⋅ ), respectively.Remark 3.4.In the scope of envisaged application classes described in Section 3.1, publications 61,84,88 showcase the suitability of the proposed interface constitutive model for describing the mechanical behaviour of slender structures embedded in different soil types, while ensuring a high level of numerical robustness.However, the authors believe that it is somewhat premature to claim the general applicability of the proposed relations in Equations ( 22)- (24) given the wide field of potential EB-SSE application cases. 58Hence, further validation studies are in demand to prove the generality of the above relations, or detect limitations for practical use.In the latter case, Γ 1 , Γ 2 and Γ 3 can readily be modified to achieve the desired system behaviour.
As a novelty, the interface constitutive model of the proposed EB-I is equipped with a limit criterion which couples the ultimate skin resistance to both the stress state and the soil parameters of the TSFE.The response along the perimeter is, therefore, locally controlled at each x by means of a penalty regularization of Coulomb's friction law, in combination with the Kuhn-Tucker-Karush condition 42 in order to exclude the occurrence of tensile stresses perpendicular to the interaction surface: In analogy to Tschuchnigg et al., 51 this approach is referred to as layer-dependent (LD) option, indicating that the ultimate skin resistance EB-I is a result of the analysis; see Tschuchnigg and Schweiger 51 for alternative options.Therein,  ′ and  ′ are the effective shear strength parameters of the TSFE associated with x .It should be noted that the effective normal stress σ′  considered in Equation ( 29) is coaxial with the normal traction constituent t1 of t with its origin placed at x .However, it is pointed out that they are not necessarily equal in magnitude.While t1 is directly related to the relative displacement field, as described in Equation ( 21), σ′  is obtained from the known effective stress tensor at the Gauss points  ′ ,GP belonging to the TSFE; see Remark 3.5.As schematically illustrated in Figure 7B, the mapping procedure includes four operations.At first, the stress components computed at the Gauss points are projected from the solid Gauss points to the solid nodes of the TSFE using a nodal extrapolation scheme 99 : where σ′ ,  represents the effective stress tensor extrapolated to the th node belonging to the TSFE,  GP denotes the number of Gauss points occupied by the TSFE,  ′ ,  ′ and  ′ are the stretched intrinsic coordinates of the TSFE, while   GP and  ′ , GP constitute the element function and the stress tensor associated with Gauss point  GP .Subsequently, the effective stress tensor at x can be determined using the isoparametric interpolation concept, similar to Equation ( 17): (32)   with σ′  being the nodal stress vector obtained from Equation (31).Since the stress quantities are evaluated in the global coordinate system  1  2  3 of Ω  , the mapping procedure needs to consider two rotations in three-dimensional space, which are defined by where ,   ∈ ℝ 33 are the (transposed) orthogonal rotation matrices that describe the sequential transformation from the global to the local coordinate system of the beam axis, and from the local coordinate system at the beam axis to the local coordinate system at Γ , .Finally, σ′  can be expressed as the Euclidean norm of the stress component vector: where  is the unit outward normal vector with respect to Γ , .
The definition of the ultimate base resistance  max that can be mobilized across Γ , conceptually differs from the LD option described above, as it is provided as direct input to the analysis.The selection of this approach is attributed to the wide range of possible values, which may be well in excess of the ultimate structural strength of Ω 100 ; hence, in many cases,  max is not of great importance to design and performance considerations, as it is very often the base resistance mobilization rate that governs the numerical performance under both serviceability and ultimate limit state conditions, which in turn is highly controlled by the definition of the embedded interface stiffness at the base expressed by Equation (24).Viewed in the context of ultimate limit state design, however,  max can be used as a safeguard in FEA, and allows to define limiting values based on empirical data ranges, such as documented in German Geotechnical Society, numerical benchmark analyses 58 or experimental results. 101We note that in the limit case,  max is internally averaged over Γ , , satisfying Remark 3.5.After initialization of an EB-I, t1 is locally computed at each x based on the incremental change in ũrel,1 representing the relative displacement between Ω  and Ω  in the perpendicular direction to Γ , ; compare Figure 7A.Initially, the material points belonging to Ω  and Ω  are assumed to share the same position across the interaction surface.With respect to the activation phase of a vertically oriented EB-I, it follows for ∀ x ∈ Γ , : From Equation ( 36), it can be inferred that zero-valued normal traction constituents t1 tend to underestimate the mobilized confining stresses along the perimeter, giving rise to overly conservative estimates of the shear resistance along the interaction surface throughout the simulation.To overcome this shortcoming, σ′  considered in the Coulomb-type limit criterion formulated in Equation ( 29) is directly obtained from the stress state in TSFE, which in turn naturally accounts for both the initial stress field and the gradual evolution upon loading.It is, therefore, concluded that σ′  addresses the normal stress history at Γ , in a more realistic way.

Numerical implementation
The proposed EB-I model is general in the sense that the mapping procedure described in Section 3.3 allows for the modelling of arbitrarily oriented structures, without occurrence of "geometrical" inaccuracies, such as explored in Turello et al. 53 This is attributed to the point collocation method used to assemble the global system of equations, with which local coupling contributions, defined in Equations ( 18) and ( 19), stemming from x that lie outside Ω  can be omitted by the subroutine used to loop over Γ  in order to form the linearized system of equations.In this paper, the formulation is implemented for 10-noded tetrahedral solid FE with C0-continuous element functions of second order.It is noted that this is not a limiting factor inherent to the presented EB-I framework, that is, it can be applied with different types of solid FE.Due to the introduction of discrete coupling forces localized inside the solid FE, the exact solution to BVP involving embedded FE models may contain singularities; hence, optimal spatial convergence towards the exact solution is not generally expected. 61A simple way to overcome this issue is to define an elastic zone which comprises the virtual beam geometry as well as the half sphere adjoining the pile base. 51,86In these regions, solid FE stress points are assigned with a stiffness similar to the soil stiffness and forced to remain elastic, yielding local stress field approximates that may be regarded as unrealistic.However, this is not considered critical within the scope of envisaged application classes since the mechanical response is well covered by the Timoshenko beam.
In this manuscript, the global system of non-linear equations that can be formally written as is iteratively solved using the multi-core solver PICOS 102 which achieves parallelism through the domain decomposition technique.In Equation (37),  is the global stiffness matrix, Δ û is the increment vector describing the iterative change in sought nodal field quantities,  +1  is the external force vector applied at the current load step and    is the internal reaction force vector obtained from the previous load step.A detailed description of the incremental-iterative solution process including stopping criteria is beyond the scope of the present manuscript and can be found elsewhere. 65 I G U R E 8 (A) Analysis domain and boundary conditions, (B) mesh discretizations, (C) cross-sectional view of shaft and base coupling point patterns.Due to limited space, only one half and one quarter of the FE models are described in (A) and (B), respectively, although full 3D models are used.

VERIFICATION AND NUMERICAL VALIDATION
The next section uses a single-pile problem to analyse the numerical characteristics of the EB-I model proposed in this work.Although single piles are rarely constructed in isolation, it is useful to consider their analysis for the performance assessment of novel EB-SSE; for example, see Refs. 62,63,86It should be noted, however, that the proposed formulation can be applied to a variety of technical problems comprising ERM; see Section 1.The parametric studies are designed to gain insight into two key model features of the proposed EB-I, namely the combined coupling scheme and the nonlinear interface constitutive model, and to assess its computational efficiency.All FEA are carried out using a research calculation kernel of the FE code Plaxis 3D Connect Edition V22.01.00.452, installed on a conventional personal computer with 64-bit Windows 11 Pro OS, Intel Core i7-1165G7 2.8 GHz processor, 16 GB of available RAM and a solid-state drive.

Investigated scenario and model description
In the course of this section, the well-documented Alzey Bridge static pile load test 103 executed in slightly overconsolidated Frankfurt clay serves as reference scenario to study the numerical performance of the proposed EB-I model.Figure 8A depicts the model geometry and boundary conditions, together with the groundwater table that is located 3.5 m below the ground surface.In an additional series of FEA, the pile domain Ω  is discretized employing both the EB-L model and the SFEA model described in Tschuchnigg and Schweiger 51 for purposes of comparison.Standard zero-thickness interface elements 104 are introduced at the pile-soil interface with the same properties of the soil.In order to reduce the effect of singular plasticity points developing close to the pile edge, these elements are extended beyond the physical pile boundary; compare van Langen and Vermeer. 105As discussed in Section 2, this SFEA model represents a D-SSE modelling approach and is, therefore, regarded as numerical benchmark; compare Figure 3.The pile material parameters are listed in Table 1.
As shown in Figure 8B, Ω  is discretized with a varying number of 10-noded tetrahedral solid FE; unless otherwise stated, the results are obtained with 8329 (EB-L, EB-I) and 84,148 (SFEA) solid FE.For the sake of numerical consistency, all simulations are carried out using full 3D models, instead of axisymmetric or reduced 3D models for the SFEA.This allows for a direct comparison between both EB-SSE and the SFEA in terms of computational performance.Likewise, all simulations are carried out employing identical calculation phase sequences; these include the initial stress field generation employing the K0 procedure, followed by the activation of Ω  as well as displacement-driven and force-driven head-down loading, respectively.Figure 8C provides a visual representation of the standard EB-I coupling point patterns considered along Γ , and across Γ , .While these coupling point patterns comply with respect to the number of perimetral coupling points x (•), , the base coupling point pattern additionally incorporates a central coupling point x  .In order to avoid co-local x (•), at the perimeter of the base interaction surface (i.e., pairs of x , and x , sharing the same initial spatial coordinates), x , are specified at an offset ratio r∕R = 0.75.
Ω  is modelled with the Hardening Soil Small model (HSS) with non-associated flow rule and small-strain stiffness overlay, 106 an extension of the Hardening Soil model (HS). 107Table 1 lists the HSS material parameters; while the calibrated HS parameters and drainage conditions are defined in accordance with Engin et al., 86 HSS-specific model parameters, namely the initial shear modulus,  0 , and the threshold shear strain for stiffness degradation,  0.7 , are specified on the basis of empirical correlations documented in Refs. 108,109

Pile resistance mobilization
Despite the remarkable developments in the field of EB-SSE explained in Section 2, the coupling formalism at the pile base has received little attention to date, with the notable exception of Tschuchnigg 43 where this aspect is investigated in the context of EB-L.As a crucial enhancement of the proposed EB-I, its combined coupling scheme additionally accounts for SSI evolving between Ω  and Ω  at the base in the form of x , placed across Γ , (Figure 8C), that is, Γ  = Γ , ∪ Γ , .It should be pointed out that this feature extends existing EB-I models with purely shaft coupling documented in Refs. 62,63here Γ  = Γ , ; see also Figure 3.While this simplified coupling procedure may reasonably describe the mechanics of related problem classes where the base response will rarely influence the global pile response (e.g., laterally loaded piles or friction piles), it undermines the EB-I calculation fidelity in geotechnical problems with a potentially significant base resistance mobilization.This is demonstrated in Figure 9A, which shows the load-settlement curves under vertical loading for the cases of shaft and combined 2D-to-3D coupling, respectively.While the results obtained with the shaft coupling scheme in the first stage of loading are in reasonable agreement with the numerical benchmark, it becomes apparent that this method fails to accurately represent the pile response at load levels   > 1800 kN with dominant pile resistance mobilization rates.The proposed combined coupling scheme, in contrast, captures the numerical benchmark throughout the complete calculation phase sequence; apparently, the incorporation of Γ , extends the EB-I applicability, such that it is no longer limited to loading situations with insignificant base resistance mobilization.For instance, this allows for more realistic estimates of load-sharing mechanisms traditionally of interest in the design of pile-supported foundations; for example, see Nguyen et al. 110 Figure 9B studies the development of skin tractions along the pile length computed with the combined coupling scheme.In order to obtain the equivalent results in the SFEA model, the stress field at Ω  is post-processed by means of integration at the pile cross-section and filtered to remove unwanted noise triggered by stress gradients at both pile ends. 104At   = 1500 kN, the EB-I achieves a reasonable agreement with the numerical benchmark, with approximately linear skin traction increase along the pile shaft despite an almost constant distribution of relative displacement in axial direction (Figure 9C).This behaviour is a favourable consequence of the stress-dependent formulation of the interface constitu- tive model formally described in Equations ( 21)- (24), in which t is directly coupled to    computed at the stress points belonging to the TSFE; compare Figure 6A.From the power law exponent  = 1.0 specified in Table 1, it follows that these coupled quantities have to develop linearly with depth. 43,107Skin traction oscillations observed with the EB-I represent an intrinsic feature of EB-SSE, which are attributed to the non-smoothness of the relative displacement field; the interested reader may refer to Refs., 24,95 where this aspect has been discussed in detail.The EB-I skin traction distribution computed at a relative pile head displacement u∕D = 3% is automatically limited by the Coulomb-type formulation of the limit criterion, in which σ′  is recovered from TSFE; compare Figure 7B.Apparently, the EB-I behaves similar to the numerical benchmark where standard zero-thickness interface elements are used to limit the skin resistance.The corresponding relative displacement profile in axial direction shown in Figure 9C implies the transition from stick to slip behaviour at all shaft coupling points x , , that is, the embedded interface stiffness in axial direction adopts zero values.

Coupling point pattern analysis
The choice of the coupling point pattern is important for the mathematical properties of the discretized system, since it influences the connectivity between the DOF occupied by the beam FE and solid FE, respectively.To the best of the authors' knowledge, this aspect along with the sparsity structure of K, a crucial factor that may significantly influence the performance of iterative solvers, 59 will be investigated for the first time in the context of EB-I models.
Let us consider the shaft coupling point pattern illustrated in Figure 8C, in which eight x , are placed at Γ , .As can be inferred from Equations (25) to (27), this coupling point configuration results in Δ  being linked to a maximum of eight solid FE, depending on the mesh discretization; compare Figure 2E.In extreme cases of high occupation density (Figure 10A), in which each potential TSFE is occupied by at least one x , , the element stiffness matrix describing the coupling between three-noded beam FE (with three translational and rotational DOF per   ) and eight 10-noded TSFE (with three translational DOF per solid node   ) becomes    ∈ ℝ 738×738 , because 3 (  ) ⋅ [6 (DOF∕  ) + 3 (DOF∕  ) ⋅ 10 (  ∕TSFE) ⋅ 8 (TSFE∕  )] = 738 DOF.In the presence of both loose and dense arrangements with varying occupation density, however, the size of    may vary between consecutive three-noded beam FE, thereby adding complexity to the global stiffness matrix assembly procedure.To overcome this obstacle and simplify the assembly process, we pursue a point  collocation method where the point-wise contributions obtained from Equations ( 18) and ( 19) are separately assembled to  as well as the internal reaction force vector.
Figure 10B shows the number of non-zero entries (nnz) in  for different mesh discretizations and a varying number of x , indicating the occupation density.nnz around the main diagonal imply interactions between different nodes, either belonging to Ω  or Ω  .Consequently, counting the number of nnz provides a direct picture of the nodal connectivity. 111dopting this analysis strategy, the results highlight two interesting characteristics with respect to the sparsity structure of .First, the number of nnz increases with increasing mesh refinement.This observation can be intuitively explained by the fact that the number of DOF in the ERM increases with increasing mesh density, resulting in an increase of nnz as each DOF is associated with one or more nnz in .The second observation is that the number of nnz (or equivalently the sparsity pattern) converges to a certain state beyond which the number of nnz remains constant, regardless of the number of x , .This refers to the situation where all potential TSFE are occupied by at least one x , .With reference to the coarse and medium mesh, this state is reached at #x  , = 64 and #x  , = 128, respectively.With reference to the fine mesh, Table 2 summarizes selected numerical values describing the sparsity structure of  including the bandwidth. 38he latter produces constant values irrespective of #x  , , indicating a similar level of memory consumption for storing .Although not shown in this paper, the same holds true for the coarse and medium mesh.A more detailed analysis of the sparsity structure of , together with an improvement of the iterative solver performance on the basis of improved deflation and preconditioning methods is at the core of ongoing research and will be presented in future work.
The choice of the coupling point pattern could become critical with regard to the potential occurrence of unwanted contact locking effects triggered by an over-constraining of the ERM, that is, coupling point configurations with relatively high occupation density may produce inaccurate results that are regarded too stiff.In the context of 1D-to-3D coupling,  this is a thoroughly studied topic that can be mathematically attributed to a violation of the discrete inf-sup condition. 49,112igure 10C provides insight into the influence of the coupling point pattern on the mechanical pile response employing the combined coupling scheme in combination with the LD option.The results show an insignificant increase of pile stiffness with increasing number of x , , demonstrating the numerical robustness of the proposed EB-I model under vertical loading.As of now, however, the potential occurrence of contact locking phenomena under more general loading situations cannot be excluded and may be revisited in future research.

Assessment of computational performance and condition number
The aim of this study is to assess the computational efficiency of the proposed EB-I. Figure 11A compares the total runtime duration recovered from EB-I simulations with different coupling point configurations.It should be noticed that the presented total runtime comprises the setup time, required for setting up the pre-conditioner, and the iteration time for solving the system of equations, and is obtained by running the simulations until the out-of-balance force obtained in the loading phase with u∕D = 3% is less than the default tolerated global error of 1%. 65The presented speedup ratio, defined as the total runtime of the EB-I model divided by the total runtime consumed by the SFEA model, serves as comparative measure to assess the runtime efficiency relative to the numerical SFEA benchmark.The reader should note that this measure is also deployed in Table 3 to assess the runtime efficiency of the EB-L model.It can be inferred from the figure that an increase in #x  , leads to an increase of total runtime and, inversely, a decrease of the computed speedup ratio.With reference to the coarse mesh, a maximum speedup ratio of 8.1 is reported for #x  , = 4, while a minimum value of 4.5 is achieved with #x  , = 256.A similar trend with reduced values is found for the fine mesh.Succinctly, this observed tendency is not surprising as it can be explained by an increasing computational effort required to keep track of the coupling variables between Ω  or Ω  , which is, in turn, influenced by the coupling point configuration.For example, this concerns the setup time period consumed by the search algorithm to detect the TSFE belonging to the individual x , followed by the solution of a local system of equations required for the projection of the point on the beam centerline into the solid FE parameter space; compare Section 3.3.With reference to Section 4.3, an increase in #x  , may lead to assembly procedures that must be regarded cumbersome for practical problems.The computational cost may, therefore, be considered as roughly proportional to #x  , .
Table 3 additionally compares the computational efficiency of the proposed EB-I with the EB-L; it is important to note that deviations in the tabular data are caused by the different coupling schemes as it is the only variable in the sensitivity analysis.Interestingly, the EB-I produces generally smaller iteration numbers as well as similar speedup ratios close to the lower bound of considered #x  , , despite the relatively more expensive coupling formalism, compare Section 2.2.To investigate the numerical origin of this aspect in more detail, we evaluate the condition number  2 describing the wellposedness of K.This is because the performance of iterative solvers is closely linked to  2 , whereas smaller values indicate more efficient simulations and faster convergence to the exact solution. 24,97Since K is symmetric positive definite, 113  2 () =   ()   () where   and   are the maximum and minimum eigenvalues of .To facilitate the numerical calculation of  2 , K is recovered from the coarse mesh.Preliminary analyses of K indicate that  2 does not significantly change as the load level changes.Obviously, the EB-I reduces  2 by more than one order of magnitude compared to the EB-L, that is, the conditioning of K is significantly improved.This is a beneficial feature of the EB-I as it speeds up the total runtime compared to the EB-L and, to some extent, compensates for the more expensive coupling scheme.From a solid mechanics perspective, relative merits of the EB-I with respect to  2 can be reasonably interpreted by re-assessing the mechanical problems illustrated in Figure 11, for which analytical solutions, valid under the premise of elasticity, exist. 92With regard to the EB-L, the definition of the coupling forces in the 1D-to-3D coupling scheme is a problem similar to the Kelvin problem (Figure 11B), while the 2D-to-3D coupling scheme employed with the EB-I can be interpreted as a generalization of the plane Kelvin problem (Figure 11C).It should be noted that the exact solution to both mixed-dimensional coupling problems exhibits local singularities in the displacement and stress fields close to the position of concentrated contact interaction, that is, beam centerline for EB-L and interaction surface for EB-I.One of the main consequences of the distributed mobilization of coupling forces encountered with the EB-I is that it reduces the influence of the singularities, which is reflected in the reduced condition number with respect to K. In this context, it should be pointed out that relative merits of the EB-I are further amplified due to the proposed changes in the coupling formalism at the pile base, where the redistribution of coupling forces suppresses unwanted coupling force localization effects.
Although our numerical evidence seems to imply that the computational efficiency is optimal for small #x  , , however, one could argue that coupling point patterns with low occupation density may lead to the non-fulfilment of basic consistency tests, such as described in Steinbrecher et al. 49 ; needless to say that spurious predictions may have adverse consequences in the interpretation of the results.Of course, related studies concerning the accuracy of the solutions and convergence rates are generally favoured by a regular subdivision of Ω  , 114 a condition that is not met in the simulation framework used in this manuscript where we generally deploy an unstructured mesh with tetrahedral solid FE.As a consequence, we can only suggest to produce similar studies to that presented in Figures 10C and 11A.Nevertheless, as a first empirical recommendation deduced from the results of numerous simulations involving the EB-I, some of which are presented in this paper, we suggest #x  , = 8.

CONCLUSIONS
This work has been motivated by the growing demand for slender structure modelling concepts that combine a high level of computational efficiency with a reasonable level of numerical fidelity, with particular emphasis on embedded FE models (EB-SSE) where pile-type structures are discretized by 1D structural FE.The key scientific contributions of this work can be split into three parts.First, significant contributions related to EB-SSE are theoretically discussed and put into context in the form of an EB-SSE classification scheme.To the authors' knowledge, the latter represents the first approach to structure-related developments with clear reference to the field of computational geotechnics, and pinpoints a number of relevant limitations that should be considered in future research.
The second part introduces the theoretical framework of an enhanced EB-SSE formulation, in which the interaction between the Timoshenko beam FE and surrounding solid FE is established by means of an implicit surface-to-volume (2Dto-3D) coupling scheme comprising non-linear coupling constraints enforced at multiple coupling points.As a novelty of the proposed embedded FE model with implicit interaction surface (EB-I), the combined 2D-to-3D coupling formalism additionally accounts for endpoint interaction, a crucial feature considering the potentially significant mobilization of base resistance encountered in many geotechnical problems.Relevant components of the proposed EB-I framework are mathematically described and visually explained, whereas particular focus is placed on the variational principles, mapping scheme and formulation of the non-linear interface constitutive model.
In part three, we investigate the performance of selected EB-I model features on the basis of numerical studies concerning a single pile problem under vertical loading.The first set of analyses showcases the beneficial consequence of the enriched coupling scheme, which allows for more realistic estimates of the mechanical pile behaviour at load levels fairly above the ultimate pile skin resistance.Moreover, it is shown that the LD normal stress recovery technique considered by the Coulomb-type limit criterion enforced along the shaft interaction surface captures the response of the numerical benchmark solution with reasonable accuracy.The second set of analyses investigates the role of the coupling point configuration with regard to the global stiffness matrix K properties and mechanical pile response.The results indicate that the proposed 2D-to-3D coupling scheme is beneficial in terms of global stiffness matrix conditioning.This is a particularly relevant feature of the EB-I as it improves the computational performance of the iterative solver compared to the classical EB-L with geometrically simplified 1D-to-3D coupling schemes, making the EB-I competitive for practical applications despite the computationally more expensive coupling formalism.Although it is observed that an increasing number of coupling points has detrimental effects on the sparsity structure of K as well as the total runtime duration, the determination of a generally optimal coupling point configuration has been out of reach, as it depends on the problem at hand.Nevertheless, within our studies, we have observed a relatively small influence of this variable on the numerical predictions, which makes the proposed EB-I model suitable for a wide range of geotechnical problems.Further understanding of the numerical performance under more complex loading situations will be instrumental to moving the proposed EB-I forward to eventually harness its capabilities for large-scale 3D FE simulations involving a high number of slender structures.Furthermore, future work will investigate the capabilities of alternative normal stress recovery procedures considered in the ultimate skin resistance limit criterion, as well as the extension of the 2D-to-3D coupling scheme to explicitly account for non-circular cross-sections.

A C K N O W L E D G E M E N T S
The first author would like to thank Laurin Hauser for his interest and constructive comments.Open access funding enabled and organized by the Graz University of Technology.

C O N F L I C T O F I N T E R E S T S TAT E M E N T
The authors declare no conflicts of interest.

D ATA AVA I L A B I L I T Y S TAT E M E N T
Data supporting findings of this work are available on request from the corresponding author.

F I G U R E 1
Schematic representation of sub-scale element modelling techniques -(A) implicit, (B) explicit-discrete and (C) explicit-embedded (Note: For clarity, 3D solid FE nodes are not shown).

F
I G U R E 2 (A-D) Conceptual layout of embedded FE models with implicit interaction line/surface and (E) concept of independent displacement field of 1D structural FE (Ω 1 ) and 3D solid FE (Ω 3 ).

F I G U R E 4
Basic notation used in the development of the proposed EB-I model with implicit interaction surface Γ  including reference systems of interest.
3 ∈ [0,  beam ].The position vector of a generic local point at the beam axis is defined as   =    ⋅    .• The proposed EB-I coupling formalism enforces coupling constraints on multiple coupling points x distributed over the perimeter Γ , as well as the base Γ , of the implicit coupling surface.With reference to Equations (1) and (2), this can conveniently be expressed by Γ  = Γ , ∪ Γ , .Depending on the beam loading mode, this results in the distributed mobilization of discrete coupling forces; see Figure 5A.The corresponding orthonormal base unit vectors are denoted as {ẽ 1 , ẽ2 , ẽ3 }.In Section 3.3, we will introduce a consistent mapping procedure to describe the interaction behaviour in a discrete manner by means of beam and solid displacement fields ũ , ũ ∈ ℝ 3 computed at x , assuming that the intersection boundary of both domains initially occupies the same position; see Figure 5B.• Due to inherent limitations of EB-SSE addressed in Section 2, Ω  is discretized without subtracting the beam volume from the solid volume, resulting in the spurious overlapping of material points at identical positions.This modelling F I G U R E 5 (A) Mobilization of coupling forces over the perimeter of the implicit interaction surface; (B) schematic representation of displacement vectors at a selected coupling point.

F
I G U R E 6 (A) Spatial description of interaction surface segment Γ , as well as (B) illustration of complete set of Γ , and specification of coupling point patterns.tensors that satisfies (  )(  − x ) =   × (  − x ), ∀   ,   , x ∈ ℝ 3

F
I G U R E 7 (A) Definition of embedded interface stiffness parameters and (B) mapping procedure considered in layer-dependent option to recover σ′  from the stress field of the solid Gauss points.

F
I G U R E 9 (A) Influence of base coupling surface Γ , on load-settlement response; (B) skin resistance mobilization and (C) relative displacements  ,3 at different loading conditions.In (C),  ,3 is the mean value of ũ,3 averaged over x , belonging to the same cross section.

F
I G U R E 1 0 (A) Schematic representation of coupling point-to-solid arrangements; (B) influence of coupling point pattern on (B) nnz in  obtained at u∕D = 3% and (C) load-settlement response.TA B L E 2 Global stiffness matrix properties of the Alzey Bridge model (fine mesh).

F
I G U R E 1 1 (A) Mesh-dependent computational performance; embedded (B) point load and (C) line load acting on infinite solid, representing the Kelvin and plane Kelvin problem, respectively.TA B L E 3 Comparison of computational costs as function of sub-scale element modelling approach.
2 ,  3 } are described by means of orthonormal base unit vectors { 1 ,  2 ,  3 }; by adopting the Einstein notation, the position vector of a generic global point is defined as  =   ⋅   .Likewise, it is useful to define Overview of constitutive model parameters considered in Alzey Bridge model.
TA B L E 1