A mixed 3D-1D finite element formulation for analysis of geomaterial structures with embedded curvilinear inclusions: application to load transfer in mooring anchor systems

A finite element model for structural analysis of media with embedded inclusions is presented. The “embedded element concept” is adopted to model the contact interaction of two medium components along the contact interface considering a mixed 3D-1D formulation. The Mohr-Coulomb interface model is employed to define the bond-stress and bond-slip relation and strains associated with bond-slip are assumed to remain infinitesimal along the interface. Nonlinear analysis is performed with a corotational kinematics description introduced in the context of embedded approach. The problem of load transfer in mooring anchor systems was investigated and reasonable results were obtained using the present model.


INTRODUCTION
Many problems of structural engineering as well as geotechnical and petroleum engineering involve linear or curvilinear inclusions embedded in a solid material matrix. Reinforced concrete is the most common example with reinforcing steel bars modeled as cable elements and incorporated into the finite elements referring to the concrete material (see for instance Manzoli et al. (2008), Oliver et al. (2008) or Figueiredo et al. (2013), to cite some recent references among the numerous contributions in the field). Similar approaches, such as those implemented by Zhou et al. (2009) or Maghous et al. (2012), may also be used in geotechnical applications, especially for load evaluations in anchoring systems employed in offshore oil platforms, where soil-mooring line interaction can be characterized using a mixed 3D-1D finite element formulation. In this case, the kinematic and constitutive descriptions of the interface phenomena play a key role.
It is underlined that when the medium consists in a homogeneous matrix reinforced by several linear inclusions that are arranged periodically, the homogenization method or multiphase modelling appear as alternatives to handling the matrix/inclusion interaction problem (see for instance Bernaud et al. (1995), de Buhan and Sudret (1999), Sudret and de Buhan (2001), Bennis and de Buhan (2003), de Buhan (2006), de Buhan et al. (2008), Bernaud et al. (2009), Hassen et al. (2013), to cite a few).
Mooring systems are utilized in the offshore petroleum industry to maintain floating platforms attached to the exploitation site. A mooring system is basically composed by mooring lines and anchors, which is submitted to hydrodynamic and aerodynamic loads applied on the floating structure and transferred to the mooring lines through fairlead points located on the platform. Owing to friction forces developed along the mooring line, the load applied on the anchoring device is not the same as that observed on the platform. Taking into account that the mooring system depends totally on the strength of the anchor, dissipation produced by friction forces acting on the soil-cable interface along the buried segment of the mooring lines must be determined in order to properly evaluate the load applied on the anchoring device, considering that the friction forces acting along the line segment immersed in the sea water are determined. Previous works on this subject may be found in Degenkamp and Dutta (1989), Neubecker and Randolph (1995), Yu and Tan (2006) or Wang et al. (2010aWang et al. ( , 2010b. The structural response of one-dimensional structures immersed in a three-dimensional solid medium requires a special system. Although a numerical model based on contact mechanics for multibody systems can be adopted for the analysis (Wriggers, 2006;Laursen, 2010), the use of the so-called "embedded element concept" is much more efficient, considering that the diameter of the embedded one-dimensional structure is significantly smaller than the typical size of the surrounding solid. The embedded formulation was introduced by Phillips and Zienckiewicz (1976) to analyze concrete structures. Since then, improvements have been made regarding the description of mechanical behavior and search algorithms for localization of the reinforcement within the solid matrix. Chang et al. (1987) modified the original formulation in order to allow straight reinforcement elements embedded in arbitrary direction with respect to the local axes of the solid matrix element. Balakrishnan and Murray (1986) presented an embedded formulation with bond-slip interface, while Elwi and Hrudey (1989) extended the two-dimensional embedded formulation by using reinforcements with curved shape. Barzegar and Manddipudi (1994) proposed a general model for spatial modeling of straight segments of embedded reinforcement using inverse mapping and a search algorithm for intersection points between reinforcement and solid elements. Owing to the use of straight elements, a refined mesh is necessary in the matrix solid field in order to define reinforcements with curved shape. Ranjbaran (1996) proposed a numerical formulation for embedded reinforcements in 3D brick elements, where full bond between the solid matrix and reinforcement are assumed. Gomes and Awruch (2001) extended the search algorithm proposed by Barzegar and Manddipudi (1994) to account for curved elements embedded in hexahedral finite elements with quadratic interpolation functions. More recent works have improved the description of the solid-structure system including flexural rigidity in order to represent the one-dimensional structure using beam models, which can be utilized to simulate soil-structure interaction problems involving piles (e.g., Sadek and Shahrour;Ninic et al., 2014).
In the context of mixed 3D-1D formulation, special emphasis must be given to properly describe the interface phenomena. A Mohr-Coulomb interface model can be adopted to define the constitutive relationship between bond-stress and bond-slip, where an elastic-plastic formulation is utilized in order to define the bond-stress evolution based on the slip motion along the tangential direction of the interface, which may be reversible or irreversible. Unlike matrix and embedded inclusion, which may undergo large strains, the bond-slip is assumed to remain infinitesimal along matrix/inclusion interface. This assumption may be viewed as first approach to assessing the kinematical description of the relative motion on the interface. Clearly enough, a more comprehensive modeling is still to be developed in order to properly take into account large relative displacements along the interface.
In this context, it is noted that the objective of this contribution is to formulate a mixed 3D-1D finite element method for analyzing solid-structure interactions in geomaterials with embedded curvilinear inclusions. In terms of engineering application, the proposed finite element modeling is applied to evaluate the load attenuation at the anchoring point of a typical mooring line utilized in anchoring systems of offshore platforms, considering the solid-structure interaction occurring along the buried segment of the mooring line.

FUNDAMENTALS OF THE EMBEDDED FORMULATION IN THE INFINITESIMAL FRAMEWORK
This section is intended to provide the basic elements of the 3D-1D mixed formulation used for implementation of the so-called "embedded model".

Geometric and kinematic setting
A basic assumption to be defined in embedded formulations is related to the strain field of the material involving a one-dimensional structure. Let us assume that the motion of any material point located on the contact surface between a curvilinear body, such as a bar or a cable, and the surrounding solid medium can be described using a local coordinate system attached to that curvilinear structure. In this system, the tangent direction is aligned with the longitudinal axis of the cable and the normal direction is related to any radial direction defined on the plane of its cross section (see Fig. 1). Depending on the physical problem, one can also assume that the contact along the normal direction is permanent and the normal component of relative displacement between the solid matrix and the cable is null, while the relative displacement in the tangential direction is not. This condition is called bond-slip and is valid when the surrounding material is sufficiently compact such that only tangential movements are permitted on the contact interface. A limit condition can be also adopted when no slip is allowed and the immersed curvilinear body adheres perfectly to the surrounding material. In the case of perfect adherence between reinforcement and the solid medium, the axial strain in the immersed cable εc can be obtained according to the following compatibility condition: . .
where li, mi and ni represent the direction cosines referring to the cable local axes (l, m, n) with respect to the global coordinate system (X1, X2, X3). When relative movement between the immersed cable and the surrounding material is considered, the displacement field in the solid matrix domain must be described using a composition of solid matrix displacements and relative displacements at material points located on the solid/cable interface. In the present work, this composition is only valid along the tangential direction of the immersed cable, whereas full bond between matrix and cable is assumed along the normal direction. Denoting by uc the cable displacement and by um that of the solid matrix, the displacement jump along the matrix-cable interface may be expressed as: T  T   1  2  3  c,1  c,2  c,3  m  1  2  3  m,1  m,2  m,3   . , . c u l l l u u u u l l l u u u   In infinitesimal framework, the weak form of the equilibrium condition applied to the mechanical system constituted by the solid matrix and embedded inclusions, takes the form: stands for the traction in axial direction of the embedded inclusion. In order to avoid volume superposition between matrix and embedded cable when rigidity and geometrical characteristics are locally similar, integration of terms referring to the volume of the matrix material may be performed considering the following spatial domain:

Constitutive equations for matrix and embedded inclusion
Relations between stress and strain are first presented to describe the constitutive equations of the matrix material considered in the present formulation. The state equation for the soil matrix phase is formulated within the framework of finite plasticity. Explicit rate-form expression involving the Jaumman derivative of the stress tensor as well as the associated corotational description is detailed in section 3.1. For the sake of clarity, the main features of the soil matrix constitutive behavior are provided in this section, restricting the description to the context of infinitesimal plasticity.
In the elastic range, the Cauchy stress tensor σm and the small strain tensor εm referring to the matrix material are related according to the following equation: where D e is the fourth-order elastic constitutive tensor for an isotropic material, which may be described as: where δij denotes the components of the Kroenecker delta and K and G are the bulk and shear moduli, respectively. The incremental form of Eq. (8) is given as: In the elastoplastic range, the constitutive equation referring to the matrix material is described in the incremental form as follows: with: and: where κ is the hardening parameter, m f is the yield function and m g is the plastic potential. The total strain rate is decomposed into elastic and plastic parts using the additive decomposition, that is: The flow rule is written as: where   is a positive scalar representing the plastic multiplier. The complete description of the elastoplastic formulation requires the yield function m f , the plastic potential m g and a hardening law to be prescribed. In the present work, the generalized formulation presented in references (Nayak and Zienkiewicz 1972;Owen and Hinton, 1980;Souza Neto et al., 2008) is adopted to describe the material models utilized here for the matrix material.
As regards the embedded inclusions, an elastic behavior is considered in the analysis, with account for geometric nonlinearities. The specific state equation for inclusions will be described in section 3.2. Meanwhile, it takes the following form in the context of infinitesimal elasticity: where σc and εc are the axial components of the Cauchy stress and small strain tensors, respectively, and Ec is the Young modulus of the constitutive material. The incremental form of Eq. (15) may be expressed as

Constitutive behavior for interface material
We address hereafter the elastoplastic constitutive formulation referring to the contact interface. The differential form of the total relative displacement on the interface is decomposed into reversible and irreversible parts, which are described according to the local coordinate system defined in terms of tangential and normal components, that is: where superscripts el and ir indicate the reversible and irreversible parts of the corresponding jump displacement components. Notice that the total relative displacement in the normal direction is zero owing to the assumption of bond-slip contact at the interface. In order to define the interval of stress states that correspond to reversible relative displacements, a yield criterion must be defined. In the present work, the Mohr-Coulomb yield criterion is adopted (see Fig. 2): where c is the cohesion and μ is the friction coefficient of the material interface. Coefficient μ is usually defined from the internal friction angle ϕ: tg    . The stress vector acting upon a current point of the interface is defined in terms of the tangential (τint) and normal (σint) stress components. The tangential stress component is updated according to the elastic-plastic formulation, where the yield criterion described above indicates if the stress increment associated with the relative displacement on the interface is obtained using an elastic constitutive equation or an elastic-plastic constitutive equation.
We deal now with the crucial issue related to how the normal vector v to cylindrical interface can be defined in the plane of inclusion cross-section (see Fig. 3). It is emphasized that this key point is inherent to the mixed 3D-1D formulation when the interface behavior is accounted for. The underlying question is: how interaction forces between a 1D structural element and the surrounding 3D continuous body are possibly modeled from a rigorous mechanical viewpoint? An approximate approach to deal with this challenging issue has been proposed in Figueiredo et al. (2013). In the subsequent analysis, we develop a heuristic approach in which interaction efforts are defined from minimization of the shear component of the stress vector generated on all the facets of the soil/inclusion interface by the stress tensor in the soil. The stress vector acting along the interface is first computed from the stress tensor σm in the solid matrix: The normal stress is simply the projection of t on the normal direction: From a computational viewpoint, tensor σm is evaluated at macroscopic points of the matrix element containing the inclusion element, coinciding with the integration points of the inclusion element. Unit normal vector v is defined in the global coordinate system (X1, X2, X3) by its orientation θ in the inclusion cross-section plane and the direction cosines li, mi and ni of the local frame (l, m, n): The solution to Equation (22) that leads to higher value of int  is hence retained.
The direction of relative plastic displacement is derived from the flow rule given as follows: where d is the plastic multiplier and g is the plastic potential. The latter is assumed under form: where    is the dilatancy angle.
The consistency condition employed in the present formulation may be expressed as: The last term of the equation above is related to material hardening when this feature is considered in the mechanical description of the material. A linear isotropic strain hardening model is assumed here to describe the cohesion evolution after yielding, which may be expressed as follows: where h is the hardening modulus. Notice that the cohesion evolution does not explicitly depend on the normal component of the stress state along the interface.
Recalling that the normal component of the total relative displacement n du is zero, one obtains: The elastic constitutive equation describing the behavior of interface written in terms of differential stress components and differential relative displacements reads: On the other hand, making use of the consistency equation together with constitutive equations (32) leads to: Finally, the plastic multiplier may be rewritten as follows: It stems from the constitutive equation of the interface following the tangential direction that: Taking into account the particular form (28) of potential g and expression (35) of plastic multiplier, the elastoplastic constitutive equation of the matrix-reinforcement interface is given as follows: Therefore, the incremental form of the constitutive equations referring to the contact interface reads: in the elastic range, and: in the plastic range. When 0   (no dilation), condition 0 n du  can be ensured by choosing /~1/ n s k k tg .

FINITE STRAIN APPROACH
Nonlinear problems are usually analyzed employing the incremental approach, where strain increments are evaluated based on the last incremental displacement field obtained from the solution of the equilibrium equations. At each iterative step, the respective stress updates are obtained according to the material behavior. When large strains are involved, the latter can be described in the corotational coordinate system adopting objective stress rate.
Although the matrix and the inclusions may undergo large strains, a fundamental assumption of the present modeling is that the strains along matrix/inclusions interface remain infinitesimal.

Corotational description of matrix particles motion
A basic issue to the corotational description of the motion of a continuum body is the decomposition of the reference frame into base and corotated configurations Haugen, 2005a, 2005b). The base configura-tion usually refers to the undeformed configuration, while the corotated one is obtained through a rigid body motion from the base configuration. The coordinate system of the corotated configuration follows the material motion and can be easily related to the global coordinate system through rotation operations. In the finite element approach, each element possesses its own set of coordinate systems located at the respective quadrature points. The starting point of the analysis is the decomposition of geometric transformation of matrix into a rigid body motion and pure deformation performed at the element level in the corotational coordinate system. Such a corotational description maintains orthogonality of the attached reference frame as illustrated in Figure 4, thus providing an effective setting for rate form formulation of the constitutive equations in the context of large strains. Note that rotation R transforming the current configuration into corotational configuration can be defined from the rotation component in polar decomposition of deformation gradient (Espath et al., 2014).
Assuming that all kinematical variables are known at the previous configuration t = tn of the matrix, the displacement field at the end of the current load step can be obtained from integration of the strain rate tensor along time interval [tn, tn+1]. The whole steps of the reasoning are performed in the corotational coordinate system, where notation  shall be used to express any field  in the corotational frame. The strain rate tensor in the corotational system is defined as: where def m v represents the velocity field associated with the deformation part of the motion expressed in the corotational system. The strain increment is computed using the mid-point integration proposed in Hughes and Winget (1980). In this procedure the velocity is assumed to be constant within the time interval and the reference configuration is attached to the intermediate configuration at t = tn+1/2 in the corotational system. Accordingly: where def m u is referred to the deformation part of the total displacement increment m u in the corotational system and n 1 2  x is the intermediate geometric configuration of the matrix element in the corotational system, which can be computed as: where n+1 2 R is the orthogonal transformation tensor performing rotation from the global coordinate system to the corotational coordinate system defined locally at the intermediate configuration 1/2 n t  of the matrix element.
The displacement increment m u referring to time interval [tn, tn+1] is decomposed at element level according to: In the above local relationship, def m u and rot m u denote respectively the deformation and rotation parts of the displacement increment defined in the global coordinate system. The increment of deformation displacements expressed in the corotational system is defined by: since the strain rate tensor is evaluated at the intermediate configuration Coordinates n x and n 1  x corresponding to geometric configurations in the corotational system at n t t  and 1 n t t   are obtained from the following transformations: where n R and n+1 R are orthogonal transformation tensors performing rotations from the global coordinate system to the corotational coordinate system defined locally at n t t  and 1 n t t   , respectively. Vectors n x and n+1 x refer to geometric configurations defined in the global coordinate system. Omitting subscripts   , 1/ 2, 1 n n n   referring to considered time, the components of the transformation R are given by: with: x r x r r r r r r r r   The Cauchy stress tensor in the global system is obtained from objective tensor transformation as follows: where m  and m  are the Cauchy stress tensors evaluated in the global and corotational system, respectively. To formulate the constitutive behavior of matrix constitutive material, a set of assumptions are stated. First, the elastic part of the deformation gradient of matrix particles is assumed to remain infinitesimal, meaning that large strains are of irreversible (plastic) nature. In addition, the constitutive material is considered as elastically isotropic and that elastic properties are not affected by the plastic strains. Under these conditions, it can be shown that the state equation formulated in rate form relates a rotational time derivative of stress tensor (Jaumman derivative) and the strain rate tensor (see for instance Dormieux and Maghous, 1999;Bernaud et al., 2002;Bernaud et al., 2006) through: where p m d denotes the plastic strain rate and m  is the spin tensor defined in the corotational system. The above equation is used for stress updates in the corotational system in the context of large strains. It is observed that the corotational spin tensor has also to be integrated over the time interval n n t t  following the same mid-point rule adopted in (41).

Description of deformation in embedded inclusion
As adopted for surrounding matrix, the kinematics of the embedded inclusion is described referring to configuration n t . With respect to this configuration, the Green-Lagrange axial strain writes: where c u is the displacement vector at any point along the embedded inclusion. Derivations are taken with respect to coordinate t l x x  along the tangential direction, as illustrated in Figure 5.
It is recalled that the displacement jump at the surrounding matrix/inclusion interface is purely tangential (see section 2.1): where l is the unit vector along the cable. Accordingly, the components of c u in the local coordinate system may be expressed from s u and the displacement components in the global coordinate system um T = (um,1, um,2, um,3) of the geometrically coinciding matrix particle: The axial strain c E is computed using the matrix notation (52): The embedded curvilinear inclusion is discretized by means of succession of linear elements. Each element is embedded within a brick element (parent element) associated with a constant tangent vector l and a constant matrix transfer R*. When perfect bonding is considered at the interface solid matrix/inclusion, the first term in the right hand-side of (53) should be removed (i.e., Keeping in mind that the approach is based on the updated Lagrangian scheme, the axial strain c E is defined referring to the last equilibrated configuration of the embedded inclusion. The latter configuration is defined in the corotational coordinate system referring to the embedding matrix element (parent element) as sketched in Assuming an elastic behavior for the inclusion constitutive material, the stress increment during time interval 1 [ , ] n n t t  is related to axial strain by means of the uniaxial linear relationship:  is the second Piola-Kirchhoff stress tensor and c E is the elastic stiffness of the inclusion constitutive material.

Finite element discretization
In the subsequent analysis the eight-node hexahedral finite element formulation with one-point quadrature is adopted to discretize the displacements of matrix particles. At the element level, the displacement vector of matrix particles between n t and 1 n t  is approximated, using Voigt notation, as follows: where  m u is the vector of nodal displacements referring to eight-node hexahedral element and m N is the 3 24  matrix aimed to generate the symmetric part of the displacement gradient at element level (see Eq. (41)). It is should be observed that matrix m B must be also evaluated considering the differential operator in the corotational coordinate system, that is operating with m B . Finally, similar procedures are used for the finite element approximation of strain rate tensor and spin tensor.
Regarding the embedded inclusions, their geometry is discretized into piecewise linear elements. Instead of nodal displacements of the inclusion elements, the approach operates with nodal displacements jump. Along each two-node element, the tangential relative displacement between matrix and inclusion particles defined in (4) is approximated from the nodal values: where  s u is the tangential displacement jump vector related to the linear finite element nodes of the embedded inclusion. Matrix Φ, whose dimension is 1 2  , contains the shape functions of the two-node linear element.
Recalling that (51) and (52)), the discretized expression of axial Green-Lagrange strain between n t and 1 n t  reads in each inclusion element: is the 1 2  matrix containing the derivatives of shape functions with respect to tangential coordinate along the inclusion element,  T is the 1 6  matrix introduced in (2). Matrix G , whose dimensions are 3 26  , is defined by: Vector  d contains 26 components of nodal displacements associated with the matrix element (brick element) and embedded linear inclusion element: Rearranging (58) yields: If perfect bonding is considered at the matrix/inclusion interface, the latter expression reduces to where matrix G  is obtained from matrix G by removing its last two columns. The subsequent developments deal with the finite element formulation of the equations governing the response of the mechanical system to a prescribed external loading. It is recalled that the mechanical system refers to the solid matrix and embedded inclusions in mutual interaction. Discretization of the weak form of balance momentum expressed at current configuration of the mechanical system results in a set of nonlinear equations. These equilibrium equations must be iteratively satisfied using the incremental approach (see Bathe, 1996), since both the stiffness matrices and the internal force vectors are functions of the current element configuration. A linearization procedure based on Newton-Raphson method and Taylor series expansion of the general internal force vector within the time interval 1 [ , ] n n t t  leads to the following global system: where subscript 1 n  denotes the current position in the time marching, while subscripts i and 1 i  refer respectively to current and previous iterative steps in the Newton-Raphson procedure applied over the time n n t t  . tan K is the tangent stiffness matrix, ext F and int F are the external and internal force vectors, respectively. U represents the global vector whose components are node values of the matrix displacement and the interface relative displacement between n t and 1 n t  . Vector U is obtained by assembling procedure to incorporate the contribution of all element nodal displacements  m u or  d .
At each iterative step, the stiffness matrix and the internal force vector are evaluated in the corotational coordinate system. In particular, the terms related to the elementary contribution of matrix (matrix element without embedded inclusion) take the following expressions: where ˆe m  stands for the current configuration of the matrix element volume in the corotational system. ep D is the stress-strain constitutive matrix defined in (49) where R is the rotation transformation matrix defined in Eqs. (45) and (46).
The finite element formulation is based on reduced integration, where hourglass control techniques are employed in order to avoid numerical instabilities, such as volumetric locking and/or shear locking. Detailed description of the stabilization procedures adopted in this work may be found in Awruch (2008, 2013).
As regards the contribution of embedded inclusions to tangent stiffness matrix and force vector, the work performed by the internal axial force  in any virtual elastic evolution c  E should be analyzed. Denoting by ˆe c  the last equilibrated configuration of the inclusion element with respect to the corotational system of the surrounding matrix element (parent element), the internal force vector related to inclusion element is defined as: The elementary tangent stiffness matrix due to inclusion contribution is obtained by linearizing the above expression of internal force vector with respect to generalized displacement vector where c E is the longitudinal elastic stiffness of inclusion material and c  is the axial Cauchy stress tensor along the inclusion element. For a single element with embedded inclusion segment, the discretized weak form of equilibrium equation writes   n 1,i n 1 n 1,i 1 n 1,i 1 After assembling all the elements, the global system takes the form given by Eq. (63), which is solved iteratively with implementation of Generalized Displacement Control Method (e.g. Yang and Shieh, 1990).

VERIFICATION TEST OF NUMERICAL PROCEDURE
In order to check the finite element implementation of the proposed model as well as related accuracy, the numerical simulation of a kind of inclusion pull-out test is analyzed in the sequel. The structure sketched in Fig. 7 consists of a reinforced cantilever beam with length L and rectangular b x h cross-section. The loading is defined by the following conditions: • body forces are neglected; • both matrix and inclusion are clamped along the plane x = 0; • a force x Pe is applied at the free end of the reinforcing inclusion, while the matrix is free of stresses along the plane x = L.
The matrix constituent material is elastic perfectly plastic with a von Mises condition for plastic yielding. A Mohr-Coulomb elastoplastic constitutive law with a non associated plastic flow is considered for the matrix/inclusion interface. The corresponding model data are given in Table 1. The finite element grid used for the simulation is shown in Fig. 7 and consisted of 110x20x20 hexahedral eight-node elements, resulting from a preliminary mesh sensitivity (convergence) analysis. The influence of interface properties on the structure response is investigated by varying the value of two typical parameters, namely the interface elastic stiffness s k and the interface cohesion c .   Table 2 summarizes the numerical predictions derived from the formulation implemented in this paper, together with the finite element solutions obtained from ANSYS software. The latter has been used considering hexahedral 20-node elements (SOLID 95) for both the matrix and inclusion domains, while quadrilateral 8-node elements (CONTA 174 and TARGE 170) have been used for the contact interface. The structural response is evaluat- at the same point, and the reaction force R applied to the inclusion at its fixed (clamped) end ( 0 x  ). As it can be observed from this comparison, a good agreement is obtained from the two distinct numerical approaches, thus providing a first validation of the proposed numerical formulation. The accuracy of the latter can be illustrated by evaluating the maximum relative difference observed between the two approaches considering the 36 simulations that are reported in

NUMERICAL APPLICATIONS
The numerical model is applied to analyze loads transfer in undrained conditions along a mooring line embedded in soil and used in petroleum exploitation units.

Problem description
Mooring systems are basically composed by mooring lines and anchors, which are responsible for maintaining floating platforms attached to the exploitation site. Special anchorage devices, such as the torpedo anchor (see Wodehouse et al. (2007) and Souza et al. (2011) for a more detailed information), have been developed for oil extraction in deep waters, where fixed anchoring points are utilized to connect the mooring lines to the anchors. In catenary mooring configuration, the suspended length of mooring line in the seawater assumes the typical shape of a free hanging line, reaching the touchdown point on the seabed surface horizontally. The final segment of the mooring lines lies embedded in the soil, comprehending the touchdown point and the anchoring point located at some position on the anchorage device. This embedded segment is usually represented considering an inverse catenary curve.
From a mechanical viewpoint, load transfer in anchorage system may be summarized as follows. As the platform drifts horizontally with the wind action and ocean current, loads are induced to the mooring system and finally transferred to the anchor (Fig. 12). Although many anchoring devices are available, torpedo anchors have proved to be a very attractive alternative for deep water anchorage. Torpedo anchors are introduced into the seabed by free fall from a previously determined altitude immersed in the seawater. The anchorage depth is determined by soil conditions and drop height, which usually varies from 30 m to 150 m. The holding capacity of a typical torpedo anchor is basically defined by the penetration depth and soil properties. The present investigation shall be restricted to load action along a mooring line segment defined between the touchdown point located on the seabed surface and the anchoring point located at the top of a torpedo anchor. Aerodynamic and hydrodynamic loads, denoted symbolically by p F , acting on the floating structure are transferred to the mooring system through fairlead points located on the platform. However, the anchor is subjected to lower loads owing to friction forces acting throughout the length of the mooring line. Therefore, the main objective of the present study is to determine the load t f applied to the anchor after dissipation produced by friction forces acting on the soil-cable interface along the buried segment of the mooring line, considering that dissipation produced in the seawater results in a load b f at the touchdown point (see Fig. 12), which is lower than that applied at the platform level, i.e., b p  f F . It should be emphasized that evaluating the load t f acting on the torpedo to load transfer in mooring anchor systems anchor is an important issue in design of such anchors, since its value is directly related to the undrained load capacity.

Geometry discretization aspects
The geometry model is shown in Figure 13. It typically consists of a parallelepiped volume of soil whose upper face lies on the seabed plane 0 y  . The finite element discretization is defined by means of hexahedral finite elements. The shape of the mooring line segment utilized in the present analyses is defined according to an inverse catenary curve, which corresponds to the geometrical form usually observed from field measurements and numerical predictions. The final form of the curve depends on the anchoring depth. An equivalent bar with circular cross-section is adopted here to represent mooring lines made up of chains (Fig. 12b). The inverse catenary curve, located on the vertical mid-plane as shown in Fig. 13a, is approximated using Lagrange polynomial interpolation, taking into account that a minimum number of points to define the embedded leg of the mooring line is available. As the mooring line is modeled using linear cable elements, all intersections between the discretized mooring line and the finite element mesh referring to the soil material must be obtained employing a numerical algorithm for intersection detection. The numerical scheme adopted in this work to define the cable elements may be summarized as follows: • A specific number of points is chosen along the cable line, which are obtained using the shape of the mooring line determined from field measurements or numerical predictions. Then, Lagrange polynomials are adopted to approximate the shape of the mooring line, where the polynomial degree is specified according to the number of interpolation points initially chosen. Now, a large number of points must be created in order to discretize the interpolation curve using small straight line segments. These segments are employed to define all the intersections between the interpolation curve and element faces corresponding to the matrix material, whereas the embedded cable and matrix material are discretized using linear finite elements. • Once the points on the interpolated curve representing the embedded cable are defined, the elements of the matrix material containing the end points of the mooring line segment are to be determined. The intersection search between the interpolation curve and the faces of matrix elements begins with that matrix element containing the first point of the interpolation curve. • Taking into account the element containing the first point of the discretized curve, the subsequent points are covered in order to find the intersection between the interpolated curve and some face of the current matrix element. The subsequent points are tested to define whether the point is contained or not in the present element of the matrix material. The first point outside the spatial domain of the matrix element and that point immediately before are chosen to construct a straight line from which the intersection point can be obtained. • After the first intersection point is obtained, the procedure described above is repeated considering the element of the matrix material in contact with the element face of the preceding element where the intersection occurs. The first point of the straight segment representing a cable segment is now considered as the last intersection point. The points initially defined on the interpolation curve are again covered from the point subsequent to the point outside the previous matrix element. The first point outside the present element is chosen as the last point of the straight line representing a cable segment. After all the points on the interpolation curve are covered, all the straight segments representing the mooring line will be defined. dimensions were chosen in order to avoid the influence of boundary effects over the present results. The finite element mesh consists of 44800 eight-node hexahedral elements, which are uniformly distributed over the computational domain. Similar characteristics have been considered for the finite element meshes referring to the anchoring points at h=20 m and h=25 m using 72000 and 109200 eight-node hexahedral elements, respectively.
Smooth-wall boundary conditions are imposed on the bottom and lateral surfaces of the domain geometry. On the upper surface, a hydrostatic load w q that stands for the seawater pressure acting at the seabed level (depth H=2135 m from sea level) is applied. As mentioned previously, the curvilinear bar representing the embedded segment of the mooring line is located on the intermediate vertical plane (z = 3.5 m) of the computational domain. The anchoring point at the end of the embedded mooring line segment is modeled considering that the anchor is a rigid body with all degrees-offreedom restricted. Hence, the anchor and the soil-anchor interface are not considered in the present approach. The lower end of the cable is positioned according to the anchoring depths analyzed here, where all the displacement degrees-of-freedom are restrained to simulate the connection between the mooring line and the torpedo anchor. On the other hand, the upper end is located on the seabed surface, where an axial load is applied representing the aerodynamic and hydrodynamic forces acting on the floating platform, and transmitted from the fairlead to the touchdown point through the immersed segment of the mooring line in the seawater. The inverse catenary shapes referring to the different anchoring depths analyzed in this work can be approximated with the following polynomial functions, according to the anchoring depth h: These functions are obtained from least square fitting applied on typical data obtained from mooring systems installed in Brazilian offshore fields.

Model data
A main purpose of the analysis is to investigate the mechanical response of soil/mooring cable system under undrained loading conditions. In the context of total stress constitutive behavior, the soil is assumed to be elastic perfectly plastic with isotropic physical properties evolving with depth. A Tresca-like plasticity model is adopted for plasticity with undrained shear strength profile: which means that u S linearly increases with depth. Parameter  [Pa/m] can be evaluated from field tests like Cone Penetration Test. As regards the soil elastic properties, the Young modulus is considered to evolve with depth according to:  being a dimensionless parameter (stiffness ratio) whose typical value ranges from 100 to 500 for normally consolidated soils. The behaviors of mooring cable and related soil contact interface have been described in sections 2.2 and 2.3. The equivalent bar with circular cross-section, which is adopted in the present analysis to model the mooring line (see Fig. 12b), is considered as elastic with axial rigidity defined by the product c c E A  . The elastic-plastic behavior of the soil/inclusion interface is defined by normal and tangential stiffness moduli n k and s k , together with cohesion c , friction angle  and dilatancy angle  .  An important aspect of the analysis is related to the initial stress distribution in the soil, i.e., the stress state prior to the application of structural load applied by the floating structure and transferred to anchor through the buried mooring cable. A comprehensive approach would theoretically rely on the analysis of whole torpedo anchor penetration process in the soil, which falls beyond the scope of the present investigation. Instead, it is assumed that the structure is loaded by the soil stress state undisturbed by the anchor installation process. More precisely, the following simplified expression for initial stress distribution is considered: are the horizontal and vertical stress components of the initial stress state.
Parameter 0 K is similar to the coefficient of lateral earth pressure at rest, and value 0 1 K  is considered for most of subsequent numerical simulations. The last component of loading mode considered in the present analysis is the structural load originating from the floating platform and transferred to the anchor system defined by soil -embedded mooring cable -interface. The load is modeled by means of an axial force b f applied to the upper extremity of buried mooring cable, which corresponds to the touchdown point located on the seabed surface ( 0 y  ) (see Fig. 13a). A typical load record referring to mooring lines installed in Brazilian offshore fields shown in Fig. 14a will be considered in the numerical simulations. The corresponding spectral frequency density function is displayed in Fig. 14b, showing that the main frequencies associated with the load record are ranging around 0.002 Hz . Note also that the mean load value is approximately b 4000 kN mean  f . Figure 15: Schematic representation of applied load and associated structural response in terms of axial displacement.
The spectral density function obtained from the numerical simulations carried out herein is shown in Figure   16. It is observed from this analysis that the main frequencies of the mechanical system are around 20 Hz , thus indicating that the main frequencies associated with the dynamic load applied to the system are much lower than those of the system cable-surrounding soil-contact interface.
This preliminary analysis suggests that although the mooring system is subjected to dynamic forces, the mechanical response of the system can be evaluated by means of numerical simulations disregarding inertial effects, i.e., considering that loads are applied under quasi-static conditions.
Clearly enough, this conclusion has been drawn in light of calculations performed with a fixed model data. Further simulations based on a parametric study by varying several problem parameters would be necessary to corroborate this result. role of this parameter in the structure response. It should be kept in mind, however, that this parameter is not commonly available from field or experimental measurements.

Further simulations
Referring to model data of Table 3, simulations with different anchoring point depths are performed to evaluate the load t f applied to the torpedo anchor, as well as associated load attenuation induced by friction along the soil/cable interface. Three configurations corresponding to anchoring points located at depths h=15 m, 20 m and 25 m , are investigated. Axial load b f obtained from spectral filtering applied to the load record presented in Fig. 14a is applied at the free end of the mooring cable. This procedure leads to an equivalent load record with a smaller time interval.  Table 4 presents the corresponding load attenuation values obtained from the numerical simulations. As expected, the magnitude of load attenuation increases with anchoring depth.  The effect of stress yielding induced by the presence of embedded cable on surrounding soil can be qualitatively assessed by visualization of the von Mises equivalent stress contours (Fig. 19). It is observed from these figures that the higher is the anchoring depth, the lower is the surrounding soil yielding (i.e., decreasing von Mises equivalent stress). The structural response is also evaluated through axial displacement and soil/cable relative displacement at the free end of the cable. For illustrative purposes, Figures 20a and 20b show the evolution in time of these quantities in the particular case of anchoring point at h=15 m . Some statistical values referring to the structural response obtained in the present analysis are summarized in Table 5.

CONCLUSIONS
A mixed 3D-1D finite element formulation for structural analysis of solids with embedded curvilinear inclusions with account for interface behavior has been presented in the context of elastoplasticity. The analysis focused on the specific situation of a geomaterial surrounding bar-like inclusions that are assumed to take only tensile-compressive forces. The embedded approach together with a corotational kinematics description is adopted for nonlinear analysis of such geo-structures.
The numerical implementation of the model to investigate the problem of load transfer in mooring anchor systems has demonstrated its aptitude to deal with and macroscopically capture essential features of deformation in a complex structure. Although the simulations have mainly concerned quasi-static analyses, the current version of finite element tool can be adopted for dynamic analyses. However, such analyses handled in 3D setting require developing more efficient computational strategies including remeshing procedure and parallel computing.
From both theoretical and numerical viewpoints, several fundamental issues should be foreseen in the future. In particular: • Development of a mechanically consistent method to properly model the interaction forces at the soil/cable interface, in the spirit of the multiphase approach proposed in Figueiredo et al. (2013). • Formulation of the interface constitutive behavior in the context of large strains, that is to address large relative displacements between embedded inclusion and surrounding matrix. It is well know that such a task is complex regarding the finite element implementation. • Extension of the model formulation to the situation of porous media. At the macroscopic scale, the hydromechanical coupling in soils, and more generally in geomaterials, can be account for in the framework of infinitesimal or finite poroplasticity. This is an essential aspect for addressing problems connected with geotechnical or petroleum engineering. This issue is actually the object of current developments.
As mentioned above, dynamic analyses in 3D nonlinear setting involve large computational cost. This includes parallel implementation of the finite element model with specific data structure storage and iterative solver.