Stéphane Miquel Aguirre

Stéphane Avril

An implicit 3D corotational formulation for frictional contact dynamics of beams against rigid surfaces using discrete signed distance fields

Keywords: Contact, 3D Corotational beams, Nonlinear dynamics, Finite Elements, Eulerian Solids, Signed Distance Fields

An implicit 3D corotational formulation for frictional contact dynamics of beams against rigid surfaces using discrete signed distance fields

Introduction

Modeling beam to surface contact has applications in many fields of engineering, such as: civil engineering ( pipeline-soil interaction [START_REF] White | Seabed characterisation and models for pipeline-soil interaction[END_REF], drill-string dynamics in oil extraction [START_REF] Goicoechea | A non-linear Cosserat rod model for drill-string dynamics in arbitrary borehole geometries with contact and friction[END_REF]), offshore and naval engineering (mooring lines and capstains [START_REF] De Cursi | Stress unilateral analysis of mooring cables[END_REF][START_REF] Lee | Coupled analysis method of a mooring system and a floating crane based on flexible multibody dynamics considering contact with the seabed[END_REF], offshore risers sliding on the seabed [START_REF] Neto | Contact between rolling beams and flat surfaces[END_REF][START_REF] Neto | Static analysis of offshore risers with a geometrically-exact 3D beam model subjected to unilateral contact[END_REF]) and mechanical engineering (deployment of flexible antennas [START_REF] Tabarrok | On the dynamics of an axially moving beam[END_REF], belt-drives used in power transmission [START_REF] Wasfy | Effect of bending stiffness on the dynamic and steady-state responses of belt-drives[END_REF], vibration of loosened joints [START_REF] Dumont | Dynamic contact of a beam against rigid obstacles: Convergence of a velocity-based approximation and numerical results[END_REF]). This problem has also attracted the interest in biomedical applications, specifically in the simulation of Minimally Invasive Surgery in which guidewires and catheters are inserted through the opaque patient tissue with limited visual aids [START_REF] Solomon | Real-time cardiac catheter navigation on threedimensional CT images[END_REF][START_REF] Wei | A (Near) Real-Time Simulation Method of Aneurysm Coil Embolization[END_REF][START_REF] Perrin | Patientspecific numerical simulation of stent-graft deployment: Validation on three clinical cases[END_REF][START_REF] Kerrien | Blood vessel modeling for interactive simulation of interventional neuroradiology procedures[END_REF]. In many of the above instances, there is also a need of fast computation, even real-time, such that decisions can be made as the physical process evolves.

The beam to surface contact problem can be accurately solved by combining Total Lagrangian formulations for the beam model, and triangulated or NURBS Lagrangian surfaces for the surface model [START_REF] Temizer | Contact treatment in isogeometric analysis with NURBS[END_REF][START_REF] Konyukhov | On some aspects for contact with rigid surfaces: Surface-torigid surface and curves-to-rigid surface algorithms[END_REF]. The overall process is computationally expensive, due to the nonlinearities involved in the beam formulation but also due to the expensive search algorithms, required to solve a minimal distance problem to find the close point projection [START_REF] Wriggers | Computational contact mechanics[END_REF]. In the pursue of real-time computations alternative approaches have been used to simplify the beam model, on one hand, and simplify the contact interaction on the other hand. Regarding the beam model, simplified approaches have been used such as one dimensional string equations [START_REF] Dumont | Numerical simulation of a model of vibrations with joint clearance[END_REF][START_REF] Liakou | Fast In-Plane Dynamics of a Beam with Unilateral Constraints[END_REF], mass-spring models [START_REF] Basdogan | Virtual environments for medical training: Graphical and haptic simulation of laparoscopic common bile duct exploration[END_REF], rigid multibody links [START_REF] Wang | A virtual-reality simulator and force sensation combined catheter operation training system and its preliminary evaluation[END_REF] or hybrid models including the last two [START_REF] Tang | A stable and real-time nonlinear elastic approach to simulating guidewire and catheter insertions based on cosserat rod[END_REF]. Regarding contact (or collision) detection, this can be simplified by using bounding boxes or spheres surrounding the target object [START_REF] Wang | A robust and fast approach to simulating the behavior of guidewire in vascular interventional radiology[END_REF]. Also, specific boundary conditions can be designed, as in the case of sliding beams in small orifices [START_REF] Vu-Quoc | Dynamics of sliding geometrically-exact beams: large angle maneuver and parametric resonance[END_REF][START_REF] Humer | Dynamic modeling of beams with non-material, deformation-dependent boundary conditions[END_REF][START_REF] Steinbrecher | On the numerical modeling of sliding beams: A comparison of different approaches[END_REF], however this is limited to a limited number of problems. In cardiovascular applications, some authors assume circular cross sections for the vessels, which allows detecting penetration [START_REF] Li | A catheterization-training simulator based on a fast multigrid solver[END_REF]. Each of these approaches provide faster computational speed at the expense of a decrease in accuracy. On top of this, most of the publications use explicit time integration, which limits the maximum time step size allowed for stability reasons. Alternatively, in the search towards an efficient, yet accurate, beam to surface contact algorithm, this paper explores using the following original contributions: 1) development of a consistent implicit 3D corotational formulation for beam against master surface contact and 2) an implicit description of the surface by means of discrete signed distance fields (SDF) defined in a voxel type background mesh. The overall approach can be seen as a Lagrangian beam immersed within an Eulerian (rigid) solid. The use of level sets and/or Eulerian solids in computational contact, has been explored before for 3D solids in [START_REF] Belytschko | A monolithic smoothing-gap algorithm for contactimpact based on the signed distance function[END_REF][START_REF] Bruchon | Using a signed distance function for the simulation of metal forming processes: Formulation of the contact condition and mesh adaptation. From a Lagrangian approach to an Eulerian approach[END_REF][START_REF] Jones | Eulerian solid simulation with contact[END_REF][START_REF] Wolff | Distance fields on unstructured grids: Stable interpolation, assumed gradients, collision detection and gap function[END_REF][START_REF] Chi | A level set enhanced natural kernel contact algorithm for impact and penetration modeling[END_REF][START_REF] Leichner | A contact algorithm for voxel-based meshes using an implicit boundary representation[END_REF] but, to the best of the authors knowledge, this is the first application where it is combined with a Lagrangian beam.

Corotational beams assume a split of the beam deformation into rigid and purely deformational parts and have been largely used in the last decades [START_REF] Rankin | An element independent corotational procedure for the treatment of large rotations[END_REF][START_REF] Rankin | The use of projectors to improve finite element performance[END_REF][START_REF] Crisfield | A consistent co-rotational formulation for non-linear, three-dimensional, beamelements[END_REF][START_REF] Pacoste | Element behavior in post-critical plane frame analysis[END_REF][START_REF] Pacoste | Beam elements in instability problems[END_REF][START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF][START_REF] Felippa | A unified formulation of small-strain corotational finite elements: I. Theory[END_REF][START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF]. This formulation allows modeling all the nonlinearities in the rigid rotation, and using linear constitutive models for the deformational part. In many thin beams applications, this allows for an accurate, yet computationally efficient solution, and makes it one of the preferred choices in real-time simulation [START_REF] Marinkovic | Survey of finite element method-based real-time simulations[END_REF]. Among the challenges faced in corotational beam formulations, comes the parametrization of 3D rotations. Numerous papers use the 'rotational vector' which allows for an additive update of the rotations but it limits its use to rotations under 2π. Instead, Battini and co-workers [START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF] have largely explored the use of a spatial form of the incremental rotation vector, which still allows for additive updates, but only at the level of the iterative corrections. In [START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF], Le and coauthors extended this formulation to dynamics, comparing very well against the total Lagrangian formulation of Simo and Vu-Quoc [START_REF] Simo | A finite strain beam formulation. The three-dimensional dynamic problem. Part I[END_REF][START_REF] Simo | A three-dimensional finite-strain rod model. part II: Computational aspects[END_REF]. Regarding corotational beams and contact, the only contribution in implicit dynamics formulation is found in [START_REF] Cho | Elastoplastic and contact analysis based on consistent dynamic formulation of co-rotational planar elements[END_REF] but it is limited to 2D planar elements. This paper proposes a 3D implicit formulation of frictional dynamic contact by building up from the formulation of Le and coauthors [START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF]. To do so, a friction penalty term is added into the variational formulation, which is integrated along the beam centerline in spatial coordinates and gives rise to a distributed line contact force [START_REF] Meier | A finite element approach for the line-to-line contact interaction of thin beams with arbitrary orientation[END_REF][START_REF] Meier | Geometrically Exact Finite Element Formulations for Slender Beams: Kirchhoff-Love Theory Versus Simo-Reissner Theory[END_REF][START_REF] Meier | A unified approach for beam-to-beam contact[END_REF][START_REF] Meier | Geometrically exact beam elements and smooth contact schemes for the modeling of fiber-based materials and structures[END_REF]. The formulation assumes that the normal gap with respect to the surface is given, as well as its spatial gradient and Hessian (which will be computed via interpolation of a SDF). To facilitate convergence of the Newton Raphson algorithm, the normal contact pressure is regularized via a quadratic contact potential [START_REF] Meier | Geometrically exact beam elements and smooth contact schemes for the modeling of fiber-based materials and structures[END_REF][START_REF] Durville | Contact-friction modeling within elastic beam assemblies: an application to knot tightening[END_REF], whereas the stick-slip behavior is regularized using a square root function [START_REF] Wriggers | Computational contact mechanics[END_REF]. The semi-discrete system of equations is solved in time using HHT-α method [START_REF] Le | A consistent 3D corotational beam element for nonlinear dynamic analysis of flexible structures[END_REF]. In addition to the contact formulation, a new description of the corotational beam kinematics is proposed, traditionally described by means of change of variables between a moving and a global reference frame [START_REF] Crisfield | A consistent co-rotational formulation for non-linear, three-dimensional, beamelements[END_REF][START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF]. To do so, finite-elasticity three dimensional theory is used, as proposed in [START_REF] Auricchio | On the geometrically exact beam model: A consistent, effective and simple derivation from three-dimensional finite-elasticity[END_REF] for geometrically exact beams. By particularizing it to the corotational formulation, interesting links are found with geometrically exact beam kinematics, while making the previous developments in corotational beams still fully valid.

Regarding contact detection, the distance from the beam to the surface is found by interpolating a discrete SDF (level set with gradient of norm one) given in a background voxel-type mesh, to the Lagrangian point of integration. Specifically, the interpolation is carried out using kernel functions as used in Material Point Method [START_REF] Sulsky | A particle method for history-dependent materials[END_REF][START_REF] York | Fluid-membrane interaction based on the material point method[END_REF] and immersed Fluid-Structure Interaction techniques [START_REF] Wang | Extended immersed boundary method using FEM and RKPM[END_REF][START_REF] Liu | Immersed finite element method and its applications to biological systems[END_REF][START_REF] Wang | The immersed/fictitious element method for fluid-structure interaction: Volumetric consistency, compressibility and thin members[END_REF][START_REF] Wang | Interpolation functions in the immersed boundary and finite element methods[END_REF][START_REF] Gil | The immersed structural potential method for haemodynamic applications[END_REF][START_REF] Gil | An enhanced Immersed Structural Potential Method for fluid-structure interaction[END_REF][START_REF] Hesch | On continuum immersed strategies for Fluid-Structure Interaction[END_REF]. Taking advantage of the Cartesian structure of the voxel mesh, a very efficient tensor product of 1D kernels is used, which allows precomputing the interpolation stencils. Specifically, the spline based 1D kernel proposed in [START_REF] Gil | An enhanced Immersed Structural Potential Method for fluid-structure interaction[END_REF] is chosen, which ensures continuity of the first and second spatial derivatives, as required by the beam contact formulation. The paper analyzes numerically the accuracy of the interpolation technique in approximating the SDF, its gradient and Hessian by comparing it against an analytical solution. Next, five examples of beam to rigid surface contact are shown to prove the accuracy and robustness of the proposed formulation.

The remainder of the paper is organized as follows. Section 2 summarizes the applied beam formulation. Section 3 presents the frictional contact penalty formulation and its linearization. Section 4 presents the interpolation used for the discrete SDF and the numerical analysis of its order of accuracy. Section 5 presents the time integration method and computational aspects. Section 6 a series of examples. Finally, section 7 presents the concluding remarks.

Beam formulation

We largely follow the notation in [START_REF] Le | A consistent 3D corotational beam element for nonlinear dynamic analysis of flexible structures[END_REF] in which a fully consistent tangent operator is proposed for the dynamics of co-rotational beams. In what follows, a global reference system is defined in the material configuration, which is given by the triad of Cartesian vectors E i (i = 1, 2, 3). Similarly, a reference system is defined in the spatial (Eulerian) configuration, given by the triad e i (i = 1, 2, 3). As in [START_REF] Ibrahimbegović | Stress resultant geometrically nonlinear shell theory with drilling rotations-Part II[END_REF], both triads are made coincident, i.e. E i = δ ij e j with δ ij being the Kronecker delta. They are denoted with different symbols to differentiate between material and spatial objects.

Parametrization of finite rotations

Rotations are parametrized by means of the rotational vector [START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF][START_REF] Crisfield | Dynamics of 3-D co-rotational beams[END_REF][START_REF] Ibrahimbegović | Finite rotations in dynamics of beams and implicit timestepping schemes[END_REF], defined as

θ = θn, (1)
with n the axis of rotation and θ the rotation angle given by

θ = √ θ • θ. (2)
The relationship between the rotation angle and the rotation matrix is provided by the Rodrigues formula, given by

R(θ) = I + sin θ θ S(θ) + 1 -cos θ θ 2 S(θ)S(θ) = exp(S(θ)), ( 3
)
where I is the 3 × 3 identity matrix and S(θ) is the skew-symmetric matrix of θ, given by

S(θ) =   0 -θ 3 θ 2 θ 3 0 -θ 1 -θ 2 θ 1 0   . ( 4
)
The variation of R(θ) is obtained by first constructing a superimposed infinitesimal rotation [START_REF] Simo | A three-dimensional finite-strain rod model. part II: Computational aspects[END_REF]. This is done by using the exponential map (3) and by noticing that, since R is a two point tensor, the superimposed rotation needs to be a spatial quantity, i.e.

R = exp( S(w))R (5)
and, therefore,

δR = d d R | =0 = S(δw)R, (6)
with δw and S(δw) spatial objects defined as

δw = w i e i , (7)
S(δw) = [S(δw)] ij e i ⊗ e j , (8)
where δw is named the spatial spin vector and are infinitesimal spatial rotations superimposed to the rotation R. Instead of a spatial spin vector, a material one can be used, defined as

δW = δW i E i , (9)
S(δW ) = [S(δW )] ij E i ⊗ E j , (10)
where δW and S(δW ) are the pull back of δw and S(δw), respectively, i.e.

δW = R T δw, (11)
S(δW ) = R T S(δw)R. (12)
Using ( 12) into ( 6) yields [START_REF] Simo | A three-dimensional finite-strain rod model. part II: Computational aspects[END_REF] δR = R T S(δW ). [START_REF] Kerrien | Blood vessel modeling for interactive simulation of interventional neuroradiology procedures[END_REF] The relationship between the spatial spin vector and the variation of the rotational vector are given by

δw = T s (θ)δθ, (14a)
δθ = T -1 s (θ)δw, (14b)
with

T s (θ) = I + 1 -cos θ θ 2 S(θ) + θ -sin θ θ 3 S(θ)S(θ), (15a)
T -1 s (θ) = (θ/2) tan(θ/2) I 3×3 + 1 θ 2 1 - (θ/2) tan(θ/2) θθ T - 1 2 S(θ). (15b)
The spatial form of the angular velocity and acceleration can be obtained from equation [START_REF] Neto | Static analysis of offshore risers with a geometrically-exact 3D beam model subjected to unilateral contact[END_REF] as

S( ẇ) = ṘR T , (16a)
S( ẅ) = RR T + Ṙ ṘT , (16b) with
( •) = d dt (•) and ( •) = d 2 dtdt (•).
Finally, the spatial angular velocity can be calculated from the time derivative of the rotational vector via equation (14b) to give ẇ = T s (θ) θ.

(17)

Corotational beam kinematics

As an addition to the formulation presented in [START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF], the corotational beam kinematics is described using finite-elasticity three dimensional theory as proposed in [START_REF] Auricchio | On the geometrically exact beam model: A consistent, effective and simple derivation from three-dimensional finite-elasticity[END_REF] for geometrically exact beams.

A straight beam with circular cross section is defined with deformed configuration Ω and material configuration Ω 0 . Two additional configurations, Ω0 and Ω are also defined. These configurations contain, respectively, the positions of the beam in the initial and current positions absent of rigid body motions. Both configurations are assumed to be very close, i.e. Ω Ω0 and therefore they are uniquely denoted by Ω. The basis vectors of Ω are Ēi (i = 1, 2, 3) and are made coincident with the material and spatial basis, i.e. Ēi = δ ij E j , Ēi = δ ij e j .

To define the beam position, an orthogonal system T i (i = 1, 2, 3) is placed at the centroid of the beam cross section in the material configuration. This system defines the local coordinates {ξ 1 , ξ 2 , ξ 3 } with the arc-length parameter s = ξ 1 ∈ [0, l 0 ] that coincides with the axis of the beam, defined by T 1 . The cross sections of the beam lie in planes defined by the basis vectors {T 2 , T 3 }, which depend on the arc-length parameter s. Accordingly, the beam is deformed into a current configuration Ω, with an orthonormal frame t i defined also at the centroid of the beam cross section. These vectors are defined through the following transformations

T i (s) = Λ 0 (s) Ēi ; Λ 0 (s) = T i (s) ⊗ Ēi = [Λ 0 ] ij E i ⊗ Ēj , (18)
t i (s, t) = Λ(s, t) Ēi ; Λ(s, t) = t i (s, t) ⊗ Ēi = [Λ] ij e i ⊗ Ēj . ( 19
)
with Λ 0 , Λ being two point tensors that give, respectively, orthonormal triads T i , t i (material and spatial objects, respectively) as a function of the basis vectors Ēi . By comparing relationships ( 18) and ( 19), a map between T i and t i can be defined as

t i (s, t) = R(s, t)T i (s); R(s, t) = Λ(s, t)Λ -T 0 (s) = t i ⊗ T i = [R] ij e i ⊗ E i , (20)
which is the mapping from the material to the current configuration which in fact uniquely describes the rotation of the cross section. Using the above definitions, the material and spatial positions of the beam are given by

X(s, ξ 2 , ξ 3 ) = X 0 (s) + ξ 2 T 2 (s) + ξ 3 T 3 (s), (21) x(s, ξ 2 , ξ 3 ) = x 0 (s) + ξ 2 t 2 (s) + ξ 3 t 3 (s), (22)
where X 0 (s), x 0 (s) are, respectively, the position of the beam centerline in the material and current configurations. In order to remove rigid body rotations, two mappings are defined between Ω and Ω 0 and between Ω and Ω. To do so, two frames that move rigidly with the beam axis are considered [START_REF] Pacoste | Beam elements in instability problems[END_REF]. This is a material frame V i (i = 1, 2, 3) and a spatial frame v i (i = 1, 2, 3 ). In order to compute v i (i = 1, 2, 3), the procedure detailed in [START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF] is used,

v 1 (t) = x 0 (l 0 , t) -x 0 (0, t) x 0 (l 0 , t) -x 0 (0, t) ; v 3 (t) = v 1 (t) × p(t) v 1 (t) × p(t) ; v 2 (t) = v 3 (t) × v 1 (t), (23)
with the auxiliary vector p given by

p(t) = 1 2 (t 2 (0, t) + t 2 (l 0 , t)), (24)
with an identical procedure to obtain, i.e V i = v i (0) (i = 1, 2, 3). The above basis vectors allow defining the following mappings

V i = R 0 e Ēi ; R 0 e = V i ⊗ Ēi = R 0 e ij E i ⊗ Ēj , (25)
v i (t) = R e (t) Ēi ; R e (t) = v i (t) ⊗ Ēi = [R e ] ij e i ⊗ Ēj , (26)
with R 0 e being a two point tensor that goes from Ω to Ω 0 and R e a two point tensor that goes from Ω to Ω. It is worth noticing that R e (0) = R 0 e and that neither R 0 e nor R e depend on the arc-length parameter s, and hence define a rigid body rotation of the beam segment. These mappings allow defining the cross sectional triads in the rigid body free configuration Ω, i.e.

Ti (s) = (R 0 e ) T T i (s), ( 27
) ti (s, t) = R T e (t)t i (s, t). (28)
Equivalently to equation [START_REF] Wang | A virtual-reality simulator and force sensation combined catheter operation training system and its preliminary evaluation[END_REF], a mapping can be defined between Ti and ti as

ti (s, t) = R(s, t) Ti (s), R(s, t) = R T e (t)R(s, t)R 0 e = ti (s, t) ⊗ Ti (s) = [ R] ij Ēi ⊗ Ēj . ( 29
)
where R is a small rotation of the cross sectional vectors Ti in Ω, now absent of rigid body motions. Assuming the material beam as straight, the rotation matrix Λ 0 can be made coincident with R 0 e , this is Λ 0 (s) = R 0 e , while T i = V i and Ti = Ēi . Taking all this into account, and omitting the dependence with respect to t for clarity purposes, the material and spatial positions of the beam centerline are given by

X 0 (s) = X 0 (0) + sV 1 , ( 30
) x 0 (s) = x 0 (0) + sv 1 + u 1 (s)v 1 + u 2 (s)v 2 + u 3 (s)v 3 , (31)
where u i (s) (i = 1, 2, 3) are the positions of the beam centerline respect the rigid frame v i . Given the definition of v 1 in equation [START_REF] Vu-Quoc | Dynamics of sliding geometrically-exact beams: large angle maneuver and parametric resonance[END_REF], the displacements u i (s) follow

u 1 (0) = 0; u 1 (l 0 ) = ū; ū = l n -l 0 , (32a) u 2 (0) = 0; u 2 (l 0 ) = 0, (32b) u 3 (0) = 0; u 3 (l 0 ) = 0. ( 32c
)
By removing rigid body rotations using R 0 e , R e and rigid body translations using the position of the centerline node at s = 0, the position of the beam points in the Ω configuration can be defined as

X(s, ξ 2 , ξ 3 ) = X0 (s) + ξ 2 Ē2 + ξ 3 Ē3 ; X0 (s) = s Ē1 (33) x(s, ξ 2 , ξ 3 ) = x0 (s) + ξ 2 R(s) Ē2 + ξ 3 R(s) Ē3 ; x0 (s) = s Ē1 + u 1 (s) Ē1 + u 2 (s) Ē2 + u 3 (s) Ē3 . (34)
which defines the beam deformational displacement ū = x -X in Ω as

ū(s, ξ 2 , ξ 3 ) = ū0 (s) + ξ 2 ( R(s) -I) Ē2 + ξ 3 ( R(s) -I) Ē3 ; ū0 (s) = u 1 (s) Ē1 + u 2 (s) Ē2 + u 3 (s) Ē3 , (35)
which is assumed to be very small, i.e. x X. Via equation [START_REF] Rankin | An element independent corotational procedure for the treatment of large rotations[END_REF], it is worth noticing that s = X1 , ξ α = Xα .

It is useful now to write the deformation gradient using finite-elasticity three dimensional theory. Following Auricchio and coauthors [START_REF] Auricchio | On the geometrically exact beam model: A consistent, effective and simple derivation from three-dimensional finite-elasticity[END_REF] this is defined as

F (s) = ∂x ∂ξ i ⊗ ∂ξ i ∂X = R(s) + a(s) ⊗ V 1 , (36)
where equation [START_REF] Wang | A virtual-reality simulator and force sensation combined catheter operation training system and its preliminary evaluation[END_REF] has been used alongside with T i = V i . The vector a is defined as

a(s) = γ(s) + ξ α κ α (s), (37)
with γ and κ α defined as (see [START_REF] Auricchio | On the geometrically exact beam model: A consistent, effective and simple derivation from three-dimensional finite-elasticity[END_REF])

γ(s) = x 0 (s) -t 1 (s), ( 38
) κ α (s) = a α (s). ( 39
)
where (•) = d(•) ds . Remark 2.2.1. In comparison with [START_REF] Auricchio | On the geometrically exact beam model: A consistent, effective and simple derivation from three-dimensional finite-elasticity[END_REF], it is worth noticing that R is the rotation from the orthonormal basis V i in the initial configuration (straight, but not parallel to the Ēi axis) to the orthonormal triad t i in the current configuration. In [START_REF] Auricchio | On the geometrically exact beam model: A consistent, effective and simple derivation from three-dimensional finite-elasticity[END_REF], the orthonormal triad in the material configuration coincides with the reference axes, i.e. V i = δ ij Ēj and, therefore, Λ 0 = I and, via equation [START_REF] Wang | A virtual-reality simulator and force sensation combined catheter operation training system and its preliminary evaluation[END_REF], R = Λ.

The deformation gradient [START_REF] Pacoste | Element behavior in post-critical plane frame analysis[END_REF] can be decomposed using a left decomposition to yield

F (s) = R(s)A r (s); A r (s) = I + a r (s) ⊗ V 1 (40)
with a r being, as defined in [START_REF] Auricchio | On the geometrically exact beam model: A consistent, effective and simple derivation from three-dimensional finite-elasticity[END_REF], the rotate-back of a in the material configuration, i.e. a r (s) = R(s) T a(s).

Equation [START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF] shows that the deformation gradient can be decomposed by a pure rotation, on the left, followed by a pure material stretch, on the right. Alternatively, another interesting result can be obtained for the corotational formulation, which links the formalisms described in [START_REF] Auricchio | On the geometrically exact beam model: A consistent, effective and simple derivation from three-dimensional finite-elasticity[END_REF], with the decomposition between rigid body motions and small deformational displacements as described in, for example, [START_REF] Crisfield | A consistent co-rotational formulation for non-linear, three-dimensional, beamelements[END_REF][START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF]. In fact, using the rigid mappings defined in equations ( 25), [START_REF] Li | A catheterization-training simulator based on a fast multigrid solver[END_REF], the following multiplicative decomposition can be carried out

F (s) = R e F (s)(R 0 e ) T (42)
with F being the deformation gradient in the configuration Ω, which is defined by

F (s) = R(s) + ā(s) ⊗ Ē1 (43)
with ā being now the rigid-rotate-back of a to the configuration Ω, i.e.

ā(s) = R T e a(s) = γ(s) + ξ α κα (s), (44) with γ
= x 0 (s) -t1 , (45) κα
= R T e t α (s) = (R T e t α (s)) = ā α = R (s) Ēα . (46)
where the fact that R e is independent of the arc-length coordinate s has been used.

Remark 2.2.2. Differently from a r in (41), ā has been obtained in [START_REF] Cho | Elastoplastic and contact analysis based on consistent dynamic formulation of co-rotational planar elements[END_REF] by removing the rigid body motions via R e and hence the denomination rigid-rotate-back. The relation between the two is given by a r = R 0 e RT ā.

Remark 2.2.3. Equation (40) decomposes the deformation gradient into one material tensor, A r , followed by one rotation, R, both dependent on the arc-length parameter s. On the other hand, equation [START_REF] Simo | A finite strain beam formulation. The three-dimensional dynamic problem. Part I[END_REF] decomposes the deformation gradient into one rotation (R 0 e ) T , followed by one deformation gradient F , followed by another rotation R e . In this case only F is dependent on the arc-length coordinate, which shows why the corotational formulations allow writing the variational equations purely in Ω.

Assuming small rotations in Ω, the rotation tensor R can be parametrized by a first order approximation of (3),i.e.

R(s)

= I + S( Θ(s)) ( 47
)
where Θ is the rotation vector in the configuration Ω. At the same time, assuming also small displacements, it can be proved that the gradient of the displacements ū defined in [START_REF] Crisfield | A consistent co-rotational formulation for non-linear, three-dimensional, beamelements[END_REF] is given by (see Appendix A)

∇ ū = S( Θ) + ā ⊗ Ē1 . (48)
Using equations ( 47) and ( 48) into (43) yields the following approximation of F

F I + ε(∇ ū) + Ω(∇ ū) (49)
which is the additive decomposition of the deformation gradient.

Remark 2.2.4. Equations ( 42) and [START_REF] Durville | Contact-friction modeling within elastic beam assemblies: an application to knot tightening[END_REF] show that, in the corotational approach, the global deformation gradient can be seen as the multiplicative decomposition of a rigid rotation (from Ω 0 to Ω), followed by a small deformation (in Ω), followed by another rigid rotation (from Ω to Ω).

From equation [START_REF] Simo | A finite strain beam formulation. The three-dimensional dynamic problem. Part I[END_REF], the Green-Lagrange strain tensor can be computed as

E = R 0 e Ē(R 0 e ) T , (50) with Ē
= 1 2 ( F T F -I) = Ē11 Ē1 ⊗ Ē1 + Ē1α Ē1 ⊗ Ēα + Ēα1 Ēα ⊗ Ē1 , (51)
and

Ē11 = 1 2 (ā • ā) + (ā • t1 ) , (52a)
Ē1α = Ēα1 = 1 2 (ā • tα ) , (52b)
where equation [START_REF] Simo | A three-dimensional finite-strain rod model. part II: Computational aspects[END_REF] has been used alongside equation [START_REF] Jones | Eulerian solid simulation with contact[END_REF]. The total Green-Lagrange strain tensor depends only on ā, this is the derivative of the local centerline displacement, ū 0 , the local rotation matrix R and its derivative along the centerline R (see equations [START_REF] Cho | Elastoplastic and contact analysis based on consistent dynamic formulation of co-rotational planar elements[END_REF][START_REF] Meier | A finite element approach for the line-to-line contact interaction of thin beams with arbitrary orientation[END_REF][START_REF] Meier | Geometrically Exact Finite Element Formulations for Slender Beams: Kirchhoff-Love Theory Versus Simo-Reissner Theory[END_REF]). Assuming again small strains, via equation ( 49), the Green-Lagrange strain tensor can be written as

Ē = 1 2 ( F T F -I) ε(∇ ū) (53)
which proves that by assuming small strains and small rotations in Ω, Ē approximates to the small strain tensor ε(∇ ū). Using equations (50) [START_REF] York | Fluid-membrane interaction based on the material point method[END_REF], the variation of the internal energy can be computed as

δW int = Ω 0 S : δE dΩ 0 = Ω σ : δ ε d Ω = Ω [σ 11 δ ε11 + σ1α δ ε1α + σα1 δ εα1 ] d Ω ( 54
)
where the notation ε ≡ ε(∇ ū) has been used alongside the definition of

σ = R 0 e T SR 0 e . (55)
Therefore, the internal virtual work can be written in terms of linear elasticity components in Ω, which contributes in the efficiency of the corotational formulation, as all the nonlinearities are taken care of via the elemental rigid body rotation R e . Specifically, the following relationship between stresses and strains in Ω is used

σ = λ tr( ε)I + 2G ε ( 56
)
where λ, G are, respectively, the first and second Lamé parameters. Using the approximation E λ + 2G [START_REF] Auricchio | On the geometrically exact beam model: A consistent, effective and simple derivation from three-dimensional finite-elasticity[END_REF], the stress components contributing into equation ( 54) can be written as

σ11 = (λ + 2G)ε 11 E ε11 , (57a) σ1α = σα1 = 2Gε 1α = 2Gε α1 , ( 57b
)
where E is the Young's modulus. Combining equations ( 54), ( 57), the internal energy of the beam element can be written as

W int = l 0 1 2 E A ε2 11 dA + 1 2 G A (2ε 12 ) 2 + (2ε 13 ) 2 dA ds. ( 58
)
where Ω = l 0 A alongside the fact that the cross sectional area remains unchanged have been used. Finally, and before moving to the Finite Element discretization, it is useful to derive the relationship between the variation of the different rotation tensors involved in the corotational formulation. From equation ( 29), the following relationship is obtained

R(s) = R T e R(s)R 0 e . (59)
The variations of the above rotation tensors are obtained via equation [START_REF] Neto | Static analysis of offshore risers with a geometrically-exact 3D beam model subjected to unilateral contact[END_REF]. Omitting from now on the dependence on s, t for clarity, this is

δR = S(δw)R, ( 60
) δR e = S(δw e )R e , (61)
δ R = S(δ W ) R. ( 62
)
The relationship between these variations is key in the development of the co-rotational formulation. To do so, the variations δw and δw e are firstly rigidly rotated-back to the local configuration, i.e.

δW = R T e δw, (63
)

δW e = R T e δw e , (64)
which after some algebra (see [START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF] for details) and in combination with equation ( 59) yields the following relationship

δ W = δW -δW e . (65)

Finite Element Discretization

The beam centerline is discretized in N initially straight elements of length l

(e) 0 such that l 0 ∪ N l (e)
0 .1 Each of these segments is defined by two nodes a = 1, 2 with material positions X a and spatial positions x a = X a + u a , with u a being the nodal displacement. The cross section rotation at the node is defined by the nodal rotation R a , such that

u 1 = u 0 (0); R 1 = R(0), (66)
u 2 = u 0 (l 0 ); R 2 = R(l 0 ), (67)
(see Figure 1). The rotation matrices R a are parametrized via the rotational vector θ a , while their variation requires the spatial spin vector w a , i.e.

R a = exp(S(θ a )); δR a = S(δw a )R a , (a = 1, 2) (68)
which defines the elemental vector of global displacements and rotations, d, and their variation δd

d = u T 1 θ T 1 u T 2 θ T 2
T ; δd = δu T 1 δw T 1 δu T 2 δw T 2 T . ( 69
)
Following the procedure explained in the previous section, two orthonormal triads V i , v i are build tangent to the axis connecting nodes 1 and 2 in the material and spatial configurations, respectively, which are rigidly rotated using R 0 e , R e (see equations ( 25), ( 26)). These mappings allow to do rigid-rotate-back of the variations δd in the Ω configurations, via

δD = E T δd; δD = δU T 1 δW T 1 δU T 2 δW T 2 ; E = ∂D ∂d = diag(R e ). (70)
Additionally, the nodal centerline displacements and rotations in Ω, ūa and Ra are defined as

ū1 = ū0 (0); R1 = R(0) = R T e R(0)R 0 e = R T e R 1 R 0 e , (71) ū2
= ū0 (l 0 ); R2 = R(l 0 ) = R T e R(l 0 )R 0 e = R T e R 2 R 0 e , (72)
and via equations ( 32), [START_REF] Crisfield | A consistent co-rotational formulation for non-linear, three-dimensional, beamelements[END_REF], the displacement nodal values are given by ū1 =

  0 0 0   ; ū2 =   ū 0 0   (73)
and therefore the only non-zero component is ū, this is, the displacement of node 2 in the X1 direction. At the same time, the rotation matrices Ra are parametrized via the rotational vector Θa , while their variation requires the spin vector Wa , i.e

Ra = exp(S( Θa )); δ Ra = S(δ Wa ) Ra (a = 1, 2). ( 74
)
From the nodal quantities ū, Θ1 , Θ2 the centerline displacements ū0 (ξ), and the cross section rotational vector Θ(ξ) can be interpolated at any point of the centerline. To do so, the Bernoulli hypothesis with linear interpolation for the axial displacement and axial rotation and cubic interpolation for the transverse displacements and rotations is chosen [START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF].

ū0 = u 1 (ξ)E 1 + u t (ξ), (75)
with

u 1 (ξ) = N 2 (ξ)ū; u t (ξ) =   0 u 2 (ξ) u 3 (ξ)   = P 1 (ξ) Θ1 Θ2 . (76)
The rotations are interpolated as

Θ(ξ) = P 2 (ξ) Θ1 Θ2 , (77)
with the expression of P 1 , P 2 given in Appendix B. The above interpolations allow defining a set of local nodal quantities and their variations, which will be used in the discretized variational equation ( 54), i.e.

D = ū ΘT 1 ΘT 2 T ; δ D = δū δ ΘT 1 δ ΘT 2 T . ( 78
)
with the following relationship between δ D and δd

δ D = Bδd, (79)
with B given in Appendix C. On top of this, the dynamic and frictional terms requires the variation and time derivatives of the centerline displacement u 0 as well as those of the centerline spin-vector w 0 . From equation ( 75), the centerline displacement in the current configuration is given by

u h 0 = u 1 N 1 (ξ) + u 2 N 2 (ξ) + R e u t (ξ). ( 80
)
The variation of u h 0 is obtained from the equation above

δu h 0 = R e H 1 E T δd, (81)
and the first and second time derivatives as

uh 0 = R e H 1 E T ḋ, (82)
üh 0 = R e H 1 E T d + R e C 1 E T ḋ, (83)
with H 1 , C 1 given in Appendix D. Regarding the rotational variables, equation ( 65) is used to write

δ W h 0 = δW h 0 -δW e , (84)
which yields the following variations and first and second derivatives of the rotational variables

δw h 0 = R e H 2 E T δd, (85)
ẇh 0 = R e H 2 E T ḋ, ( 86
)
ẅh 0 = R e H 2 E T d + R e C 2 E T ḋ, ( 87
)
where again the expressions of H 2 , C 2 are given in Appendix D.

δW = δW int + δW kin -δW ext = 0, (88)
with W kin being the dynamic potential and W ext the potential of the external forces. The variation of the internal energy is obtained via equation ( 58) and following the procedure described in [START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF]. This variation leads to the elemental internal forces, given by

δW int = (δd) T T int , (89)
with T int the elemental internal forces in Ω, defined as

T int = ∂W int ∂d = B T Tint ≡ F T 1 M T 1 F T 2 M T 2 T , (90)
with B the matrix expressing the relationship between δd and δ D via equation ( 79), F a the nodal force and M a the nodal moment. Finally Tint are the elemental internal forces in Ω defined as

Tint = ∂W int ∂ D ≡ N M T 1 M T 2 T , (91)
with N the elemental axial force and Ma the nodal moments in Ω. As described in [START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF], given the definition of the internal energy for a linear elastic material [START_REF] Gil | The immersed structural potential method for haemodynamic applications[END_REF] and the use small strains in Ω (51) a closed form expression can be obtained for Tint for a given set of displacements D, without resorting to numerical integration. Specifically, given the choice of Finite Element interpolation for the centerline displacements [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF] and rotations [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF], this is given by (see for example [START_REF] Crisfield | A consistent co-rotational formulation for non-linear, three-dimensional, beamelements[END_REF][START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF])

Tint = Kint D, ( 92
) with Kint = 1 l 0          AE 0 0 0 0 0 0 0 GJ 0 0 -GJ 0 0 0 0 4EI 33 0 0 2EI 33 0 0 0 0 4EI 22 0 0 2EI 22 0 -GJ 0 0 GJ 0 0 0 0 2EI 33 0 0 4EI 33 0 0 0 0 2EI 22 0 0 4EI 22          , (93)
where I 22 , I 33 are the moments of inertia respect to the X2 , X3 axes and J is the torsional constant.

In order to obtain the tangent stiffness matrix, the variation of T int is required. This is given by

δT int = K int δd, (94)
with

K int = B T Kint B + Kint (95)
and with K provided in Appendix C. The kinetic energy functional is given by (see for example [START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF][START_REF] Cardona | A beam finite element non-linear theory with finite rotations[END_REF])

W kin = 1 2 l 0 (A ρ ( u0 • u0 ) + ẇ0 • I ρ ẇ0 ) ds. ( 96
)
The variation of the potential W kin in (96) combined with equations ( 81), ( 83), ( 85), ( 86), (87) gives the inertia force vector

δW kin = T T k δd, (97)
with T k given by

T k = E l 0 A ρ H T 1 R T e ü0 + H T 2 Īρ R T e ẅ0 + H T 2 S( Ẇ0 ) Īρ Ẇ0 ds, (98)
with

Ẇ0 = R T e ẇ0 , ( 99
) Īρ = RJ ρ RT . ( 100
)
and with ü0 , ẇ0 , ẅ0 given, respectively, by equations ( 83), ( 86), (87). Finally, to obtain the tangent dynamic matrix, the variation of T k is needed. As in [START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF], this is computed considering only the contributions from the acceleration and velocity, this is, the centrifugal and dynamic matrices

δT k M δ d + C k δ ḋ, (101)
with

M = E l 0 A ρ H T 1 H 1 + H T 2 Īρ H 2 ds E T , (102)
C k = E l 0 A ρ H T 1 (C 1 + C 3 ) + H T 2 Īρ (C 2 + C 4 ) + H T 2 S( Ẇ0 ) Īρ -S( Īρ Ẇ0 ) H 2 ds E T , (103)
with the expressions for C 3 , C 4 provided in Appendix D.

Beam to rigid surface contact formulation

In order to consider the contact against a rigid surface, the weak form (88) is extended using the variation of a penalty potential δW c as

δW = δW int + δW kin -δW ext -δW c = 0. ( 104
)
Following [START_REF] Wriggers | Computational contact mechanics[END_REF], a penalty potential is defined. This potential takes into account the normal and frictional contributions of the beam contact against a rigid master surface Γ s (see Figure 2a)

δW c = l 0 0 (p N δg N + t T • δg T ) ds, (105)
with g N being the normal gap respect the master surface, g T the relative displacement in the tangential direction, p N the contact pressure and t T the frictional contact traction. All values in (105) are evaluated at the beam centerline by assuming that the beam radius is small enough to disregard the couples generated by the frictional force [START_REF] Zavarise | Contact with friction between beams in 3-D space[END_REF]. Given the fact that circular beams are used, g N is defined as

g N (s) = d s (s) -r, (106)
with d s (s) being the distance of a point on the beam centerline to the surface and r the radius of the beam, which is assumed constant. The variation of g N is given by

δg N = ∂g N ∂x • δx 0 = ∂d Γs ∂x • δu 0 = n • δu 0 , (107)
with n being the surface normal. In order to obtain the variation of g T , the centerline displacement is decomposed into its normal and tangential components, respect to the master surface [START_REF] Ortega | Energy-consistent simulation of frictional contact in rigid multibody systems using implicit surfaces and penalty method[END_REF], i.e.

δu 0 = (δu 0 • n) n + δu T = δg N ∂g N ∂x + δg T . ( 108
)
Isolating δg T in the above expression yields In order to facilitate the convergence of the Newton Raphson algorithm, both the contact pressure p N and frictional contact force t T are regularized. The contact pressure is evaluated using the penalty regularization used in [START_REF] Meier | A finite element approach for the line-to-line contact interaction of thin beams with arbitrary orientation[END_REF][START_REF] Meier | A unified approach for beam-to-beam contact[END_REF][START_REF] Meier | Geometrically exact beam elements and smooth contact schemes for the modeling of fiber-based materials and structures[END_REF][START_REF] Durville | Contact-friction modeling within elastic beam assemblies: an application to knot tightening[END_REF]

δg T = I - ∂g N ∂x ⊗ ∂g N ∂x δu 0 , (109) ġT
= I - ∂g N ∂x ⊗ ∂g N ∂x u0 . ( 110
)
p N = p N (g N ) =              pN -ε c g N , g ≤ 0 ε c ḡN -pN ḡ2 N g 2 N -ε c g N + pN , 0 < g N ≤ ḡN 0, g > ḡN (111)
with ε c being the penalty parameter, pN = 1 2 ε c ḡN and ḡN being the normal gap value at which the contact pressure starts increasing (see Figure 3a). As in [START_REF] Meier | A unified approach for beam-to-beam contact[END_REF], ḡN is taken as 10% of the beam radius. The frictional contact traction is evaluated using the regularization proposed in [START_REF] Wriggers | Computational contact mechanics[END_REF]

t T = t T (g N , ġT ) = -µp N (g N ) ġT ġT 2 + ε T . ( 112
)
with ε T being a regularization parameter (see Figure 3b). Replacing equations ( 81), ( 107), ( 109), ( 111) and (112) into the potential (105) yields the following discrete contact force vector

T c = l 0 0 EH T 1 R T e f c ds, (113)
with H 1 given in Appendix D, E given in equation ( 70) and with f c being the contact force, given by

f c = p N (d)G(d; ḋ), (114)
with The contact tangent matrix is obtained by linearizing T c with respect the displacement and velocity contributions (see Appendix E for details),

G(x; d; ḋ) = ∂g N ∂x -µ ġh T ġh T 2 + ε . ( 115
)
δT c = K c δd + C c δ ḋ, ( 116
)
with K c defined as

K c = K 1 c + K 2 c + K 3 c + K 4 c , (117)
and with each of the terms given by

K 1 c = - l 0 0 E Ŝ(H T 1 F c )G T E T ds, (118a)
K 2 c = l 0 0 N 7 l 2 n EA T 1 F c r -EGS(F c )P 1 P E T ds, (118b)
K 3 c = l 0 0 EH T 1 S(F c )G T E T ds, (118c)
K 4 c = l 0 0 EH T 1 R T e K fc ds, (118d)
where the full expression of K fc is given in Appendix E, F c is the contact force rigidly-rotated-back to the Ω configuration, i.e.

F c = R T e f c ( 119
)
and Ŝ is the operator that transforms the 12 × 1 array into a 12 × 3 matrix as

Ŝ(a) =     S(a 1 ) S(a 2 ) S(a 3 ) S(a 4 )     a I =   a 3(I-1)+1 a 3(I-1)+2 a 3(I-1)+3   ; I ∈ {1, 2, 3, 4}. (120)
Finally, the dynamic contribution is given by

C c = l 0 0 EH T 1 R T e C fc ds, (121)
with C fc given in Appendix E.

Interpolation of a discrete Signed Distance Field

The contact formulation presented in the previous section depends on the normal gap g N (x), and its spatial gradient and Hessian, i.e. ∂g N (x) ∂x , ∂ 2 g N (x) ∂x∂x . This section presents the necessary tools to use a discrete SDF to compute the above quantities. To do so, an efficient interpolation based on the tensor product of 1D kernels is presented in section 4.1. Section 4.2 assesses the accuracy of the proposed framework to interpolate the signed distance.

Kernel Functions

At any Gauss point of the beam, represented by the isoparametric coordinate ξ G , which is located in the spatial domain at x G ≡ x 0 (ξ G ), the gap function, first and second derivatives are needed, i.e. from equation (106),

g N (x G ) = d s (x G ) -r, (122) ∂g N (x G ) ∂x = ∂d s (x G ) ∂x , (123)
∂ 2 g N (x G ) ∂x∂x = ∂ 2 d s (x G ) ∂x∂x . ( 124
)
It is assumed that the SDF is provided in a Cartesian grid in a three dimensional prismatic domain

Ω s = [x o 1 , x f 1 ] × [x o 2 , x f 2 ] × [x o 3 , x f 3 ], discretized with N nodes, N = N 1 N 2 N 3 , defined as x A = x A 1 1 , x A 2 2 , x A 3 3 T ; A ≡ (A 1 , A 2 , A 3 ), (125)
with

x A i i = x o i + (A i -1)h i ; A i ∈ {1, N i }; h i = (x f i -x o i )/(N i -1); i ∈ 1, 2, 3 (126a)
and h i being the grid size in the direction i. For simplicity, from now a uniform size will be considered in all directions i.e. h = h 1 = h 2 = h 3 . At each of the voxel grid nodes, the signed distance is given, i.e.

d s (x A ) = d A s ( 127
)
and at any other point within the domain Ω s , the value of d s can be interpolated. Specifically, at the integration point of the beam, i.e.

d s (x G ) = I A (d A s ), ( 128
)
where I is some interpolant function. Given the Cartesian grid chosen for the discrete SDF, one dimensional kernel functions are chosen, as used in Immersed Boundary Fluid Structure Interaction methods [START_REF] Gil | The immersed structural potential method for haemodynamic applications[END_REF][START_REF] Gil | An enhanced Immersed Structural Potential Method for fluid-structure interaction[END_REF][START_REF] Peskin | The immersed boundary method[END_REF]. This is, for any point in x in space, the signed distance d s (x) is interpolated as

d s (x) = A V A d A s ϕ(x -x A ), (129)
where V A = h 3 and φ defined as a tensor product of 1D kernels, i.e.

ϕ(z) = φ z 1 h φ z 2 h φ z 3 h . ( 130
)
This allows computing the first and second derivatives of d s as

∂d s (x) ∂x = A V A d A s ∂ϕ ∂x (x -x A ), (131)
∂ 2 d s (x) ∂x∂x = A V A d A s ∂ 2 ϕ ∂x∂x (x -x A ), (132)
with

∂ϕ ∂z i = 1 h φ z i h φ z j h φ z k h , (133)
∂ 2 ϕ ∂z i ∂z i = 1 h 2 φ z i h φ z j h φ z k h , (134)
∂ 2 ϕ ∂z i ∂z j = 1 h 2 φ z i h φ z j h φ z k h , (135)
where i = j = k. The use of the above interpolation strategy is very advantageous when using Cartesian meshes (as the given voxel mesh), as the stencils involved in the interpolation (129) can be precomputed in advance, and therefore no search is involved in computing the distance d s .

There is a wealth of 1D kernels functions φ proposed in the immersed FSI literature that could be used in the above expressions. Yet, some of the proposed kernels do not satisfy continuity of the second derivative [START_REF] Peskin | The immersed boundary method[END_REF], which would be detrimental for the proposed as second order derivatives of φ are required (see equation 132). Instead, the spline-based kernel proposed in [START_REF] Gil | An enhanced Immersed Structural Potential Method for fluid-structure interaction[END_REF] is used, which has a support of [-3h, 3h] and ensures continuity of the second derivative. This is given by

φ(r) =                        0 r ≤ -
2 < r ≤ 3 0 3 < r (136)
which is designed, by construction, to interpolate exactly linear functions and with second-order accuracy smooth functions [START_REF] Gil | An enhanced Immersed Structural Potential Method for fluid-structure interaction[END_REF][START_REF] Peskin | The immersed boundary method[END_REF].

Numerical analysis of the signed distance interpolation error

In this section the accuracy of the interpolation strategy presented in the previous section is analyzed numerically. To do so, the analytical signed distance function of a sphere is used, which is defined as

d s (x) = (x -P 0 ) • (x -P 0 ) -R 0 , (137)
with P 0 being the location of the center of the sphere and R 0 its radius. For the purpose of the current analysis, these are chosen as P 0 = 0.5 0.5 0.5 T and R 0 = 0.3.

Five sets of voxel meshes are created in the domain Ω s = [-0.5, 1.5] × [-0.5, 1.5] × [-0.5, 1.5] with mesh sizes h = {1/4, 1/8, 1/16, 1/32, 1/64}. At each of the nodes A of a given voxel mesh, the signed distance d A s is stored as d A s = d s (x A ) via equation ( 137), which creates a SDF (see Figure 4a). At the same time, an auxiliary mesh is created with the final purpose of interpolating the SDF values created in the previous step, mimicking what will be done at the beam Gauss point x G . Additionally, the spatial gradient and Hessian are also computed at each of the nodes of this auxiliary mesh. The auxiliary mesh is centered at P 0 , has a side length of l = 1, element size h aux = 1/32 and, in order to avoid overlapping of nodes with the background voxel mesh (which could potentially benefit the global interpolation error), it is rotated 5 degrees respect to the X 1 , X 2 and X 3 axes (see Figure 4b). Additionally, at each of the points of the auxiliary mesh, the exact value of the signed distance, its gradient and Hessian can be obtained from equation (137) which, in turn, allows computing the approximation error. As it can easily be proved, the gradient and Hessian of the signed distance function (137) have a singularity at P 0 . Therefore, in order to avoid numerical issues, the error analysis is limited to the region in which |d s (x)| ≤ 0.8R 0 .

Figure 5 shows the L2-norms for the interpolated signed distance, d h s and the computed spatial gradient and Hessian ∇d h s , ∇ 2 d h s . It can be seen that all values converge to their analytical counterpart with second-order accuracy, i.e. O(h 2 ), which proves the adequacy of the algorithm, initially designed for velocity fields, to interpolate distances and their derivatives. Particularly, it is worth stressing the capacity of this interpolation strategy to compute gradients and Hessians with the same order of accuracy without having the discrete values of these operators initially. Yet, contact against such a smooth surface cannot always be warranted as sharp gradients are present in real-life geometries. In such cases, as proved for the immersed boundary method, first-order accuracy, i.e. O(h) is expected [START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Griffith | On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems[END_REF]. Importantly, the discrete contact force f c in (114) and its spatial derivatives will also approximated at the same order of the distance function d s (x) and its spatial derivatives (i.e. second-order for smooth surfaces, first-order for non-smooth surfaces) which, alongside the beam discretization will determine the contact algorithm order of accuracy.

Time integration and computational aspects

Following the discretization in equations ( 90), ( 98) and (113), and upon assembly, the discretization of the weak form (104) is given by

R(D, Ḋ, D, t) = T int (D, t) + T k (D, Ḋ, D, t) -T c (D, Ḋ, t) -F ext (D, t) = 0 ( 138
)
where R is the residual, T int , T k , T c are, respectively, the assembled vectors of nodal internal, kinetic and contact forces and F ext is the vector of nodal external forces. At the same time, D is the vector of nodal solutions, containing displacements and rotations, and Ḋ, D are respectively, its first and second time derivatives. Following equation ( 69), these are given by

D =       u 1 θ 1 . . . u N θ N       ; ∆D =       ∆u 1 ∆w 1 . . . ∆u N ∆w N       ; Ḋ =       u1 ẇ1 . . . uN ẇN       ; D =       ü1 ẅ1 . . . üN ẅN       , (139)
with N the number of nodes of the beam mesh. It is worth noticing that, in equation (138), T c depends both on D and Ḋ, due to the frictional contact term (see equations ( 114), ( 115)). Equation ( 138) is discretized in time using the HHT-α method [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF], as used in corotational beam dynamics in [START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF][START_REF] Le | A consistent 3D corotational beam element for nonlinear dynamic analysis of flexible structures[END_REF]. This yields with α a user-dependent parameter and where the superscript n (respectively n + 1) indicates that the operator is evaluated at the time step t = n∆t (respectively t = (n + 1)∆t). Finally, the system is solved via a Newton Raphson algorithm using the following linearization

R n+α = (1 + α)(T n+1 int -T n+1 c -F n+1 ext ) -α(T n int -T n c -F n ext ) + T n+1 k = 0 (140)
K∆D = -R n+α (141)
with K given by

K = (1 + α) (K int -K c -K ext ) + 1 β∆t 2 M + γ β∆t (C k -(1 + α)C c ) , (142)
with β = 1 4 (1-α) 2 , γ = 1 2 (1-2α) and K int , K c , M, C k , C c assembled, respectively, from equations (95), (117), (102), ( 103) and (121). Upon obtaining ∆D, the update of the displacement and rotational variables is carried out following the procedure described in [START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF].

At the beginning of each time step, prior to the Newton Raphson algorithm, the predictor proposed in [START_REF] Le | Dynamics of 3D beam elements in a corotational context: A comparative study of established and new formulations[END_REF] is used. Following [START_REF] Meier | A unified approach for beam-to-beam contact[END_REF], the Newton-Raphson algorithm is solved by checking the norm of the increment vector in (139) and the norm of the residual in (140). The iteration finishes once both norms are smaller than their respective prescribed tolerances, R , D , i.e. R ≤ R , ∆D ≤ D . Also, following [START_REF] Meier | A unified approach for beam-to-beam contact[END_REF], if the Newton-Raphson algorithm fails to converge after 10 iterations, the time step is divided by two. Afterwards, and once the algorithm successfully converges in 8 solves, the time step is multiplied by two (this increase stops once the initial time step is recovered).

Regarding the interpolation of the SDF, as described in section 4, this can be efficiently implemented due to the Cartesian structure of the voxel mesh. Specifically, given a Lagrangian point of integration, it is very simple to find in which voxel element it has fallen (if any). Once this is found, the interpolation (129) is carried out using a precomputed stencil. By doing this, the computation of the distance at the Lagrangian point of integration becomes a quick sum over the stencil nodal values, therefore avoiding any search. In case the Lagrangian point of integration falls outside the background voxel mesh, a numerically large value is assigned to d s , which makes p N = 0 via equation (111).

Verification and examples

Each of the SDF involved in these examples has been generated using SDFgen [71]. All examples have been simulated using 3 Gauss points for integrating the dynamic contributions and 5 Gauss points for integrating the contact contributions. The tolerances in the Newton-Raphson iterations are set to to R = D = 1 • 10 -5 . All the examples show quadratic convergence in the Newton-Raphson iterations.

Vibrations of a beam constrained by two end stops

This example is used as a benchmark test and was first published in [START_REF] Dumont | Dynamic contact of a beam against rigid obstacles: Convergence of a velocity-based approximation and numerical results[END_REF][START_REF] Dumont | Numerical simulation of a model of vibrations with joint clearance[END_REF][START_REF] Liakou | Fast In-Plane Dynamics of a Beam with Unilateral Constraints[END_REF]. A beam constrained to move in the X 1 -X 3 plane has a length L = 10 m, moment of inertia I 22 = 1.688•10 -3 m 4 , cross sectional area A = 1.4923 • 10 -5 m 2 , Young's modulus E = 2 • 10 11 Pa, density ρ 0 = 8000 Kg/m 3 and Poisson's ratio ν = 0. The beam has its left end clamped and two rigid stops located at its right end that constrain the displacement above and below X 3 = ±0.1 m. The original results published in [START_REF] Dumont | Dynamic contact of a beam against rigid obstacles: Convergence of a velocity-based approximation and numerical results[END_REF][START_REF] Dumont | Numerical simulation of a model of vibrations with joint clearance[END_REF][START_REF] Liakou | Fast In-Plane Dynamics of a Beam with Unilateral Constraints[END_REF] were obtained by modeling the one dimensional of a string and, therefore, the end-stops where modeled as a singular constrain located exactly at X 3 = 10 m. In this case, a two dimensional beam is modeled and the tip can move freely in the X 1 and X 3 directions. Because of this, the end-stops are created as two rigid prisms of 0.3 × 0.04 × 0.1 m centered at X = 10 0 ±0. [START_REF] Konyukhov | On some aspects for contact with rigid surfaces: Surface-torigid surface and curves-to-rigid surface algorithms[END_REF] T m (see Figure 6a), while non-frictional contact is imposed only at the right end node. In order to mimic the literature results, no radius is considered when computing the normal gap g N , i.e. r = 0 in equation ( 106). A uniform time-varying body force is applied along the beam, with the value b = 0 0 sin(10t) T , which causes the beam right tip to impact against the end-stops and excite high frequency modes. The beam is modeled using three choices for the number of elements, n el = {20, 80, 320} with corresponding time step of ∆t = {5 • 10 -4 , 1.25 • 10 -4 , 3.125 • 10 -5 } s. The penalty parameter is set to ε c = 1 • 10 6 . The end-stops are modeled using a SDF in a voxel-type grid 175 × 45 × 219 nodes of element size h = 0.002 m. Figure 6b shows the isovolume of this SDF.

Figure 7 shows the variation of the u 3 component of the displacement of the beam tip with time for different refinements of the beam mesh. Both the displacement and time are made adimensional using the same procedure as in [START_REF] Liakou | Fast In-Plane Dynamics of a Beam with Unilateral Constraints[END_REF]. It can be seen that high frequency modes are excited with successive impacts and chattering as the beam tip impacts the top and bottom stops. This is highlighted in Figure 8, which zooms in the right end displacement at the first impact. The solution compares very well against that of [START_REF] Liakou | Fast In-Plane Dynamics of a Beam with Unilateral Constraints[END_REF], with some differences when the beam tip is freely moving between the two stops. This difference is probably due to the two-dimensional character of the solution, as the beam tip can freely move in the X 1 and X 3 directions, which is a substantial difference with respect to the string vibration solution obtained in [START_REF] Liakou | Fast In-Plane Dynamics of a Beam with Unilateral Constraints[END_REF]. Figure 9 shows the pathlines of the beam tip in the X 1 -X 3 plane, which shows this two dimensional motion. This plot also shows how the solution converges with successive refinement of the beam mesh.

Beam against sphere

This example compares the accuracy of the discrete SDF versus the use of an analytical signed distance function. A beam of length L = 10 cm and radius R = 0.3 mm has its base point constrained to oscillate back and forth in the X 1 direction (the rest of the beam points are free to move and rotate in any direction). The beam has a Young's modulus E = 5 GPa, Poisson ratio ν = 0.33 and density ρ 0 = 7850. The beam impacts against a sphere centred at P 0 = 5 0.5 8 T cm and radius R 0 = 4 cm. Figure 10a describes the problem setup while Figure 10b shows the prescribed displacement of the beam base point. The sphere surface is defined Figure 12 shows snapshots of the solution at t = 1.25 s, comparing different voxel discretizations with the analytical SDF. Additionally, a video animation is provided in Appendix F in the electronic version of the document. Visually, it is observed that for coarser voxel meshes, the beam undergoes some penetration within the sphere surface. In order to quantify that penetration, the analytical SDF is used to compute the maximum penetration at the centerline. To do so, the centerline position is postprocessed via the interpolation defined in equation ( 80) and using 80 points per element between ξ = 0 and ξ = l 0 . At each of this postprocessed points, the penetration is evaluated via equation ( 106) and the analytical signed distance function (137). Using this procedure, Figure 13 shows the maximum normalized penetration. It can be seen how the maximum penetration is reduced with the refinement of the voxel mesh. Also, the finest voxel mesh gives maximum penetration levels of the order of the analytical solution, which is around 50% of the beam radius.

Finally, the computational efficiency of the proposed distance interpolation algorithm is assessed qualitatively. Figure 14 shows the ratio of computational time (wall clock time) to physical time when using different discretizations of the SDF or the analytical signed distance. It can be seen that when using the discrete SDF the computational time is increased around 50% as compared to the computation when using the analytical signed distance equation. The differences in computational time between the different discretizations are significantly smaller than that, which proves that the interpolation algorithm is barely affected by the size of the background mesh h. This is a significant difference with respect to traditional contact search methods, where the size of the master surface mesh plays a big role in the computational time. Again, this assessment is purely qualitative and a more accurate assessments should be carried out by comparing against standard contact methods. Figure 14 also shows that the ratios of computational to physical time are of the order of 300, approximately. This is due to the fact that the algorithm has been implemented in Matlab using dense matrices and no vectorization and therefore, is open to substantial optimizations. This would allow reducing very significantly the ratio shown in the figure and approach real-time computation (ratio = 1) for this example. Finally, it is worth mentioning that the ratio fluctuations as shown in Figure 14 are due to the adaptive time step strategy described in section 5, which forces a reduction of the time step when contact occurs.

Contact between a cable and a rigid cylinder

This example is taken from [START_REF] Wang | A virtual-reality simulator and force sensation combined catheter operation training system and its preliminary evaluation[END_REF][START_REF] Wang | Dynamic simulation of frictional multi-zone contacts of thin beams[END_REF] and shows the robustness of the contact algorithm in a friction dominated scenario. A circular cable of length 3.0 m is inclined 0.925 rad respect to the axis as shown in Figure 16. The cable is first modeled using 20 beam elements and time step size ∆t = 5 • 10 -4 . The penalty parameter is set to ε c = 1 • 10 6 and three different friction coefficients are used, µ = {0.00, 0.10, 0.30}. Figure 17 shows three snapshots of the solution for µ = 0.1 at different time instants. Figure 18 quantifies the vertical displacement of the beam tip for the different friction coefficients. Additionally, a video animation is provided in Appendix F in the electronic version of the document, where the evolution of the solution is compared between µ = 0.00 and µ = 0.30. As expected, the higher the friction coefficient, the smaller the vertical oscillations, as the cable sticks more onto the cylinder surface. Next, the cable is modeled using three choices for the number of beam elements, n el = {10, 20, 40} with corresponding time steps of ∆t = {1 • 10 -3 , 5 • 10 -4 , 2.5 • 10 -4 } s and µ = 0.30. Figure 19 shows the vertical displacement of the cable tip point for these different discretizations. As it can be seen, the difference in the solution between 20 and 40 elements is very small, in contrast to the Absolute Nodal Coordinate formulation reported in [START_REF] Wang | Dynamic simulation of frictional contacts of thin beams during large overall motions via absolute nodal coordinate formulation[END_REF], which required at least 40 elements to obtain a converged solution.

Finally, Figure 20 shows the normal and frictional force vectors at time t = 3.0 s for a friction coefficient µ = 0.1 s. As it can be seen, the normal force vectors are perpendicular to the isosurface of the discrete SDF, while the frictional force vectors are tangent to that isosurface, proving the accuracy of the SDF interpolation (see section 4) in modeling the gradients of the surface normal, ∂g N ∂x .

Impact of a ring against a rigid surface

This example is designed to test the conservation properties of the algorithm [START_REF] Wriggers | Finite element formulation of large deformation impactcontact problems with friction[END_REF][START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF]. A flexible ring with a centerline radius R = 9.5 m has a circular cross section of radius r = 0.5 m. The ring is made of flexible material with properties E = 100 Pa, ν = 1 • 10 -4 and ρ 0 = 1 • 10 -2 Kg/m3 . The ring has an initial uniform speed of v 0 = √ 2, -√ 2, 0 T m/s with its center initially located at X 0 = [0, 0, 0] T m. A rigid surface is placed at X 3 = -12 m 3 . The penalty parameter is taken as ε c = 5000, the time step size as ∆t = 0.2 s and the ring is discretized using 64 beam elements. The simulation is run using three different friction coefficients, µ = {0.00, 0.15, 0.30}. Figure 22 shows a snapshot of the solution at t = 40.0 s alongside the pathlines for µ = 0.00 and µ = 0.30. It can be seen that the reflection angle increases with increasing friction, as reported in [START_REF] Wriggers | Finite element formulation of large deformation impactcontact problems with friction[END_REF]. Regarding the conservation of energy, Figure 23 shows the evolution of the kinetic, internal and total energy for µ = 0.15. It can be seen that, by using the HHT-α method, the solution shows an overall dissipative trend, and is stable up to very long time steps, not experiencing any energy blow up, which is a crucial test in impact problems [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF].

Figure 24 compares the conservation of the total energy for different friction coefficients. The energy dissipation increases with the friction coefficient, as expected. It can also be seen that the total energy dissipates even in the frictionless case, when no physical dissipation is added into the system. A large amount of this dissipation, of the order of 9 %, occurs during the contact (and sliding) against the rigid surface. This is due to the combination of the proposed contact penalty algorithm with the non-conservative HHT-α time integrator. This result is consistent with other reported results in the literature when a contact spatial discretization is combined with a non-conservative time integrator (see [START_REF] Liakou | Fast In-Plane Dynamics of a Beam with Unilateral Constraints[END_REF][START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF] for different impact examples with numerical energy loss due to a non-conservative time integrator). Briefly explained, this is due to the nonsatisfaction of the persistency condition, i.e. p N ġN = 0, which is one of the fourth Kuhn-Tucker contact conditions, and the only one that involves a rate of the normal gap [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF]. Given that a material time derivative is involved in this condition, the joint spatial discretization and time integration should be designed to ensure its satisfaction. Failing to do this might result in energy blow up, which is not the case for the current algorithm, or dissipation of the total energy during impact (see [START_REF] Dumont | Numerical simulation of a model of vibrations with joint clearance[END_REF] for successive energy losses during consecutive impact), which is due to negative power introduced into the system during contact [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF]. Energy preserving contact time integrators have been proposed in [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF][START_REF] Laursen | Improved implicit integrators for transient impact problems -Geometric admissibility within the conserving framework[END_REF][START_REF] Acary | Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact[END_REF]. It is also worth mentioning that energy preservation in the context of beam dynamics poses also a significant challenge, as time integration of the rotational variables adds additional complexity, and has been explored in, for example, [START_REF] Crisfield | Dynamics of 3-D co-rotational beams[END_REF][START_REF] Ibrahimbegovic | Energy conserving/decaying implicit time-stepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations[END_REF].

Dropping net on the Stanford Dragon

This example is designed to show the robustness of the algorithm in complex contact scenarios. A net of 11 × 13 threads of lengths L 1 and L 2 is placed immediately above the Stanford Dragon [START_REF]The Stanford 3D Scanning Repository[END_REF], covering a larger area of that corresponding to its floor projection. The net geometry and position are detailed in Figure 25 and Table 1. The threads have a circular cross section of radius r = 2 • 10 -4 m, and their material properties are E = 10 GPa, ν = 0.0, ρ 0 = 1500 Kg/m 3 . The net is initially at rest and under the action of gravity, defined as a constant body force b = [0, -9.8, 0] T m/s 2 . Regarding the discretization, 786 beam elements are used to model the net and the initial time step is set to ∆t = 5 • 10 -4 s. The penalty parameter is chosen as ε c = 1 • 10 8 and the friction coefficient as µ = 0.1. The simulation is run from time t = 0.0 s to time t = 19 • 10 -2 s, when the net has covered the whole solid and prior to any self contact of the net. The Stanford Dragon SDF is given in a voxel-type grid of 236 × 176 × 123 nodes with element size h = •10 -3 m (see Figure 26). Figure 27 shows the solution at different time steps. Additionally, a video animation is provided in Appendix F in the electronic version of the document. During the falling process, complex interactions between the net and the different parts of the dragon geometry occur, which are accurately captured by the proposed algorithm. Figure 28 highlights the solution at the final time step, by showing different views. In order to assess the penetration, Figure 29a evaluates the percentage of Gauss points that are in contact or penetrating the body. Gauss points in contact are identified as those that have a normal gap g N ≤ ḡN = 0.1r (see section 3). Further to that, the number of Gauss points at different levels of penetration (represented as a fraction of the beam radius) are also computed. As it can be seen, all Gauss points are either contacting (0.0 ≤ g N ≤ ḡN ) or non-contacting the Dragon surface, with no Gauss point penetrating the surface. Next, the same evaluation is carried out by processing the penetration of the beam points in Figure 29b. To do so, the centerline position is postprocessed in the same way as in example 6.2, but now the distance d s (x) is evaluated using the interpolation technique described in section 4, as no analytical distance function is available. As it can be seen, the majority of the beam length in the contact region (i.e. g N ≤ ḡN ) has a penetration between the range -r ≤ g N ≤ 0.1r (i.e. the majority of the points in the contact region do not penetrate further than the beam radius). In short, the algorithm ensures that no Gauss point penetrates the surface and this results into a small percentage of the beam length penetrating the dragon surface. This is specially important as the voxel mesh element size is significantly larger than the beam radius, i.e. h/r = 5.0, which proves the accuracy of the contact algorithm proposed in section 3. To further strengthen this last point, Figure 30 shows the maximum and mean penetration of the centerline points with g N ≤ 0. It can be seen that the maximum penetration is smaller than 2r (0.4h), while the mean penetration is of the order of 0.4r (0.08h).

Conclusions

This paper pursued the development of a fast, yet accurate and robust, algorithm for beam to surface contact. To do so, a new implicit dynamics formulation for the simulation of 3D beam contact against rigid surfaces has been presented. The formulation builds from an efficient corotational beam dynamics formulation previously reported. The corotational beam kinematics has been rewritten using 3D kinematics description, linking it with geometrically exact beam theory. Frictional contact terms have been added and linearized, providing a consistent implicit dynamics formulation for frictional contact. Contact search is avoided by using a discrete signed distance field to represent the rigid surface, in which the Lagrangian beam is immersed. The contact distance, gradient and Hessian, as required by the beam formulation, are computed at the beam Gauss point by using an interpolation based in a tensor product of 1D kernels, as used in immersed fluid-structure interaction techniques. The use of this interpolation, facilitated by the Cartesian structure of the voxel mesh, is highly efficient as the stencil can be precomputed and no search is required when computing the distance to the surface.

Regarding the numerical results, it has been numerically proved that the interpolated signed distance, spatial gradient and Hessian converge with second-order accuracy to the analytical solution, i.e. O(h 2 ), for smooth surfaces while first-order convergence, i.e. O(h), is expected for sharp surfaces. This fact, alongside the beam discretization, determines the contact algorithm order of accuracy. Next, the robustness and applicability of the overall contact framework has been proven by providing a series of challenging examples. Firstly, the solution matches well against literature results in a beam fast impact problem. Next, an example has been proposed to test the accuracy of the surface discretization by comparing the penetration for different voxel sizes. Using the same parameters for the beam material and discretization, it can be seen that the penetration decreases for finer voxel mesh sizes, as expected. In the same example, the efficiency of the interpolation strategy has been assessed qualitatively, proving that the computational time is barely affected by the refinement of the voxel mesh. A literature example in a friction dominated scenario shows that the proposed algorithm converges faster than previously reported results using Absolute Nodal Coordinate formulation. Using a standard benchmark, it has also been proved that the combination of the proposed frictional contact terms with the HHT-α time integrator provide a globally dissipative solution, avoiding energy blow-ups after long physical time. Finally, the algorithm has been tested in a complex contact scenario by using a SDF of the Stanford Dragon. The penetration levels of the beam centerline are of the order of the beam radius, which is smaller than the voxel mesh size and which proves the robustness of the contact algorithm.

The proposed framework could be improved in different ways. Firstly, the algorithm has been implemented in Matlab in a non-optimized code (dense matrix, no vectorization). Therefore, the current implementation allows for substantial optimization via use of sparse matrices and lower level languages. Once this is done, and keeping in mind that both the corotational beam models and the 1D kernel interpolation have been proven to be effective algorithms in their respective areas of application, the overall framework could be an excellent candidate for accurate real-time contact simulation. Yet, a detailed comparison, in terms of efficiency, of the proposed framework against other standard contact algorithms is required in order to prove that hypothesis. Regarding the time integration, despite the HHT-α method provides a globally dissipative solution, as pointed in [START_REF] Le | A consistent 3D corotational beam element for nonlinear dynamic analysis of flexible structures[END_REF], energy-preserving time integrators could be implemented, which would improve the solution in terms of time accuracy. A significant limitation of the current formulation is that it only considers contact against rigid surfaces. In order to extend it to contact against deformable surfaces, the use of an Eulerian formulation for hyperelastic solids could be explored [START_REF] Jones | Eulerian solid simulation with contact[END_REF]. Finally, the current formulation does not consider contact (or self contact) of the beam elements. To model such phenomena, a unified beam to beam contact approach as the one presented in [START_REF] Meier | A unified approach for beam-to-beam contact[END_REF] could be adapted from a Total Lagrangian to a (Lagrangian) corotational beam formulation as presented in this paper. Overall, this would result into a Lagrangian beam to beam contact formulation combined with the presented SDF beam to surface contact formulation.

Finally, it is worth mentioning that the current algorithm is well-suited to its application in cardiovascular medicine and, specifically, simulation of guidewire navigation and stent-graft deployment [START_REF] Perrin | Patientspecific numerical simulation of stent-graft deployment: Validation on three clinical cases[END_REF][START_REF] Sharei | Navigation of guidewires and catheters in the body during intervention procedures: a review of computer-based models[END_REF]. In these applications efficient or even real time simulations are needed for training purposes and for decision making during the pre-operative and intra-operative phases. Additionally, medical imaging data is traditionally provided in voxel type meshes where SDF can be obtained automatically [START_REF] Jones | 3D distance fields: A survey of techniques and applications[END_REF]. By using the proposed framework, the creation of patient-specific triangular meshes would be avoided, which usually requires experienced user input and takes an important part of the simulation work-flow [START_REF] Neic | Automating image-based mesh generation and manipulation tasks in cardiac modeling workflows using Meshtool[END_REF]. Yet, cardiovascular applications, in general, require taking into account the deformation of the vessel wall [START_REF] Perrin | Patientspecific numerical simulation of stent-graft deployment: Validation on three clinical cases[END_REF]. Therefore, the presented work can be seen as a first step towards this final goal.

At the same time, via approximation (A.2) and by omitting the high order terms, κα can be rewritten as

κ2 = Θ (s) × T2   Θ 1 Θ 2 Θ 3   Ēi ×   -Θ3 1 Θ1   Ēi =   Θ 2 Θ1 -Θ 3 -Θ3 Θ 3 -Θ 1 Θ1 Θ 1 -Θ3 Θ 2   Ēi   -Θ 3 0 Θ 1   Ēi , (A.11) κ3 = Θ (s) × T3   Θ 1 Θ 2 Θ 3   Ēi ×   Θ2 -Θ1 1   Ēi =   Θ 2 + Θ 3 Θ1 Θ2 Θ 3 -Θ 1 -Θ 1 Θ1 -Θ 2 Θ2   Ēi   Θ 2 -Θ 1 0   Ēi . (A.12)
Using the results in equations (A.10), (A.11), (A.12) into (A.6) finally yields

ā =   u 1 -ξ 2 Θ 3 + ξ 3 Θ 2 -ξ 3 Θ 1 ξ 2 Θ 1   Ēi . (A.13) 1 W T 2 T ; δ D * = δū δ W T 1 δ W T 2 T , (C.1)
The relationship between δ D and δ D * given via equation (14a) as

δ D = Bδ D * ; B =   1 0 T 0 T 0 T -1 s ( Θ1 ) 0 3×3 0 0 3×3 T -1 s ( Θ2 )   . (C.2)
Next, the relationship between δ D * and δd is given by

δ D * = B * δd, B * = r P E T , (C.3)
with E given in equation [START_REF] Le | Dynamics of 3D beam elements in a corotational context: A comparative study of established and new formulations[END_REF]. The vector r gives the relationship between the variation δū and the variation of the global displacements, i.e. and the expressions for r and P given by

δū = rδd; r = -v 1 T 0 T v 1 T 0 T . (C.4)
At the same time, using equation ( 65), P gives the relationship between the variation of the spin variables in Ω, δ Wa and δD, i.e.

δ W1 δ W2 = P δD = P E T δd; P = 0 3×3 I 3×3 0 3×3 0 3×3 0 3×3 0 3×3 0 3×3 I 3×3 - G T G T , (C.5)
and Given the above relationship, the tangent stiffness can be obtained via equations ( 90), ( 94)

G T = ∂W e ∂D =   0 0 η/l n η 12 /2 -η 11 /2 0 0 0 -η/l n η 22 /2 -η 21 /2 0 0 0 1/l n 0 0 0 0 0 -1/l n 0 0 0 0 -1/l n 0 0 0 0 0 1/l n 0 0 0 0   (C.6) with η = q 1 q 2 , η 11 = q 11 q 2 , η 12 = q 12 q 2 , η 21 = q 21 q 2 , η 22 = q 22 q 2 . (C.
δT int = δ(B T Tint ) = δ( BT * BT Tint ) = K int δd, (C.9)
with

K int = B T Kint B + K; K = BT * Kh B * + K m (C.10)
with K given in equation ( 93) and K h given by Finally K m is given by

Kh = ∂ BT T int ∂ D * Tint =   0 0 T 0 T 0 Kh1 0 3×3 0 0 3×3 Kh2   ; Kha = ∂T -T s ( Θa ) Ma ∂ Wa Ma , (C.11) Kha = η a Θa M T a -2 Ma ΘT a + ( Θa • Ma )I 3×3 T -1 s ( Θa ) + µ i S( Θa )S( Θa ) Ma ΘT a T -1 s ( Θa ) - 1 2 S( Ma )T -1 s ( Θa ) (C.
K m = D N -EQG T E T + EGar, (C.14)
with N being the first component of Tint (see equation ( 91))), r given in equation (C.4) and D, Q, a given by

D =     D 3 0 3×3 -D 3 0 3×3 0 3×3 0 3×3 0 3×3 0 3×3 -D 3 0 3×3 D 3 0 3×3 0 3×3 0 3×3 0 3×3 0 3×3     ; D 3 = 1 l n (I 3×3 -v 1 v T 1 ), (C.15) Q =     S(Q 1 ) S(Q 2 ) S(Q 3 ) S(Q 4 )     , a =   0 η([ M * 1 ] 1 + [ M * 2 ] 1 )/l n -([ M * 1 ] 2 + [ M * 2 ] 2 )/l n ([ M * 1 ] 3 + [ M * 2 ] 3 )/l n   . (C.16)
In the above definition of Q, Q i are the block components of In order to derive the dynamic force vector as well as the tangent dynamic matrices, the variations of the centerline displacements u h 0 and their time derivatives, as well as the variations of the centerline spin vectors w h 0 and their time derivatives are needed. For completness, the main results of [START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF] are summarized, but for more details the reader is refered to that paper. The variation of the centerline displacement, δu h 0 is obtained via equation ( 80),

P T m = Q T 1 Q T 2 Q T 3 Q T 4 ; m = M T * 1 M T * 2 , (C.17
δu h 0 = δu 1 N 1 (ξ) + δu 2 N 2 (ξ) + δR e u t (ξ) + R e δu t (ξ). (D.1)
The variation δR e is obtained via equations ( 61), [START_REF] Cardona | A beam finite element non-linear theory with finite rotations[END_REF] as

δR e = R e S(δW e ). (D.2)
with (see equations ( 70), (C.6)),

δW e = ∂W e ∂D δD = G T E T δd. (D.3)
The variation of u t is given by

δu t = P 1 δ Θ1 δ Θ2 P 1 δ W1 δ W2 = P 1 P E T δd, (D.4)
where equation (C.5) has been used and where, as in [START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF], it has been assumed that T s ( Θa ) I due to small rotations (see equation 14b). Replacing the above results into (D.1),

δu h 0 = R e H 1 E T δd, (D.5) with H 1 (ξ) = N (ξ) + P 1 (ξ)P -S(u t (ξ))G T . (D.6) with N (ξ) = N 1 (ξ)I 0 N 2 (ξ)I 0 . (D.7)
The variation of H 1 is useful to compute second derivatives and also to obtain the tangent dynamic matrices. This is given by In order to obtain the tangent matrix, the variation of T k in equation ( 98) is needed. As in [START_REF] Le | Corotational formulation for nonlinear dynamics of beams with arbitrary thin-walled open cross-sections[END_REF], this is computed considering only the contributions from the acceleration and velocity, this is, the centrifugal and dynamic matrices This appendix explains in detail the linearization of the discrete contact force vector, which leads to the different terms desribed in equations ( 116), (117), ( 118) and (121). To start with, the contact force vector defined in equation ( 113) is repeated here for convenience,

δH 1 = N 7 l 2 n A 1 rδd -S(δu t )G T , (
δT c = K c δd + C c δ ḋ; K c = K 1 c + K 2 c + K 3 c + K 4 c . (E.3)
and, in order to facilitate their derivation, the following abbreviations are used (see for example [START_REF] Meier | A unified approach for beam-to-beam contact[END_REF]) where equation (D.5) has been used. The above result requires the variation δ ġh T , which is given by

K 1 c δd =
δ ġh T = Aδd + Bδ ḋ; A = A 1 + A 2 + A 3 + A 4 , (E.13)
with Finally, using the results in (E.18) into (E.4), the different tangent terms in (E.3) are given by

A 1 = - ∂g N ∂x • (R e H 1 Ḋ) ∂ 2 g N ∂x∂x R e H 1 E T - ∂g N ∂x ⊗ EH T 1 R T e ∂ 2 g
K 1 c = - l 0 0 E Ŝ(H T 1 F c )G T E T ds, (E.19)

Figure 1 :

 1 Figure 1: Co-rotational element

  (a) Continuum. (b) Discrete voxel mesh with a SDF.

Figure 2 :

 2 Figure 2: Contact of a spatial beam element against a rigid surface.

Figure 3 :

 3 Figure 3: Regularization of the normal and frictional forces.

  (a) SDF in voxel mesh (b) Auxiliary mesh (c) SDF interpolated in auxiliary mesh

Figure 4 :

 4 Figure 4: Interpolation of a SDF in a voxel mesh of h = 1/8

Figure 5 :

 5 Figure 5: Error in the discrete approximation of the signed distance of a sphere surface and its first and second derivatives. Discontinuous line indicates second-order convergence slope.

Figure 6 :Figure 7 :

 67 Figure 6: Vibrations of a beam constrained by two end stops. Initial setup and SDF

Figure 8 :

 8 Figure 8: Vibrations of a beam constrained by two end stops n el = 320. Normalized u 3 displacement of the beam tip.

Figure 9 :

 9 Figure 9: Vibrations of a beam constrained by two end stops, different discretizations. u 1 versus u 3 displacement of the beam tip.

Figure 10 :

 10 Figure 10: Beam against sphere. Initial setup and imposed displacement.

X 3

 3 axis with a lumped mass of 5.0 kg placed at its end. The cable has a radius r = 0.01 m, density ρ 0 = 7919 kg/m 3 and Young's modulus E = 1 • 10 8 Pa. The rigid cylinder has a length of 3.0 m and radius of 0.1 m. The cable has an initial angular speed w = [0, 0, 2] rad/s and is subjected to gravitational force, which is imposed through a constant body force b = ρ 0 [0, 0, -9.8] T m/s 2 . The position of the top end of the cable is set to X o = [-0.8, -0.115, 3.0] m 2 , which accounts for an initial gap of 5 • 10 -3 m (see Figure 15). The surface SDF is given in a voxel-type grid of 25 × 25 × 250 voxels, element size h = 0.0125 m and it is centered at the cylinder longitudinal

(a) 20 ×

 20 20 × 20 voxels (b) 40 × 40 × 40 voxels (c) 80 × 80 × 80 voxels

Figure 11 :

 11 Figure 11: Beam against sphere. Different SDF voxel meshes used for the simulation.

(a) 20 ×

 20 20 × 20 (b) 40 × 40 × 40 (c) 80 × 80 × 80 (d) Analytical

Figure 12 :

 12 Figure 12: Beam against sphere. View on the X 1 -X 3 plane of beam contacting the sphere for different levels of discretization of the SDF as compared to the analytical solution. Time step t = 1.25 s.

Figure 13 :

 13 Figure 13: Beam against sphere. Normalized centerline maximum penetration using the different discretizations of the SDF as compared to using the analytical signed distance function.

Figure 14 :

 14 Figure 14: Beam against sphere. Ratio of the computational (wall clock) time against physical time for different discretizations of the voxel mesh.

Figure 15 :

 15 Figure 15: Contact between a cable and a rigid cylinder. Initial Setup.

Figure 16 :

 16 Figure 16: Contact between a cable and a rigid cylinder. SDF voxel mesh.

  (a) t = 1.0 s (b) t = 1.5 s (c) t = 2.0 s

Figure 17 :

 17 Figure 17: Contact between a cable and a rigid cylinder. Snapshots of the solution at different time steps for µ = 0.1.

Figure 18 :

 18 Figure 18: Contact between a cable and a rigid cylinder. u 3 displacement of the tip point for different friction coefficients.

Figure 19 :

 19 Figure 19: Contact between a cable and a rigid cylinder. u 3 displacement of the beam tip with different discretizations for µ = 0.3.

Figure 20 :

 20 Figure 20: Contact between a cable and a rigid cylinder. Close up of the contact friction vectors at the Gauss points. Time t = 3 s, friction coefficient µ = 0.1.

Figure 21 :

 21 Figure 21: Impact of a ring against a rigid surface. Initial Setup.

  (a) µ = 0.00 (b) µ = 0.30

Figure 22 :

 22 Figure 22: Impact of a ring against a rigid surface. Comparison of the ring pathlines until t = 40.0 s for different friction coefficients.

Figure 23 :

 23 Figure 23: Impact of a ring against a rigid surface. Evolution of the kinetic, internal and contact energies for µ = 0.15.

Figure 24 :

 24 Figure 24: Impact of a ring against a rigid surface. Comparison of the total energy conservation for different friction coefficients.

Figure 25 :

 25 Figure 25: Dropping net on the Stanford Dragon. Initial setup.

Figure 26 :

 26 Figure 26: Dropping net on the Stanford Dragon. SDF voxel mesh.

(a) t = 7 •

 7 10 -2 s. (b) t = 11 • 10 -2 s. (c) t = 15 • 10 -2 s. (d) t = 19 • 10 -2 s.

Figure 27 :

 27 Figure 27: Dropping net on the Stanford Dragon. Position of the grid at different time steps.

Figure 28 :

 28 Figure 28: Dropping net on the Stanford Dragon. Different views of the solution at time t = 19 • 10 -2 s.

  (a) Percentage of Gauss points. (b) Percentage of beam centerline length.

Figure 29 :

 29 Figure 29: Dropping net on the Stanford Dragon. Stacked area plots showing the percentage of Gauss points and length contacting / penetrating the surface.

Figure 30 :

 30 Figure 30: Dropping net on the Stanford Dragon. Maximum and mean normalized penetration.

7 )

 7 Combining equations (C.2) and (C.3) yields the relationship between δd and δ D δd = Bδ D; B = B B * . (C.8)

1 M T * 2 T . (C. 18 )

 1218 ) with M T * a being the nodal moments related to the displacements D * , T * = ∂W int ∂ D * = B * Tint ≡ N M T * Appendix D. Time derivatives and variations of the centerline displacements and rotations

10 )N 7 l 2 nA 1 N 8 = N 5 + N 6 19 )C 2 =

 1021856192 D.8) with N 7 = N 3 + N 4 , (D.9)Via equation (D.5), the first and second derivatives of u h 0 can be computed asuh 0 = R e H 1 E T ḋ, (D.11) üh 0 = R e H 1 E T d + R e C 1 E T ḋ, (D.12)withC 1 = S( Ẇe )H 1 + Ḣ1 -H 1 E (D.13)and Ḣ1 obtained from equations (D.8), (D.4) asḢ1 = r ḋ -S( ut )G T ; ut = P 1 P E T ḋ. (D.14)Regarding the rotation variables, the varition of w h 0 is obtained fom equation (84) and doing a rigid-push-forward with R eδw h 0 = R e H 2 E T δd, (D.15)withH 2 (s) = P 2 (s)P + G T . (D.16)As for H 1 , it is useful to obtain the variation of H 2 . This is given by Via equation (D.15), the angular velocity and angular acceleration are obtained asẇh 0 = R e H 2 E T ḋ, (D.20) ẅh 0 = R e H 2 E T ḋ + R e C 2 ET ḋ, (D.21) with S( Ẇe )H 2 + Ḣ2 -H 2 E. (D.22)

Ũ h 0 H 1 EW h 0 H 2 EẆ h 0 H 2 l 2 nA 1 l 2 nA 2 h 1 , h 2 defined as h 1 = H 1 35 )

 010202212211135 δT k M δ d + C k δ ḋ (D.23)The following variables require linearizationŨ h 0 ≡ R T e üh 0 = H 1 (d)E(d) T d + C 1 (d, ḋ)E(d) T ḋ, (D.24) W h 0 ≡ R T e ẅh 0 = H 2 (d)E(d) T d + C 2 (d, ḋ)E(d) T ḋ, (D.25) Ẇ h 0 ≡ R T e ẇh 0 = H 2 (d)E(d) T ḋ, (D.26)where equations (83), (86) and (87) have been used for the definition of with ü0 , ẇ0 , ẅ0 . The variations of the above quantities are obtained by linearizing only respect to ḋ and d to giveδ T δ d + (C 1 + C 3 )E T δ ḋ, (D.27) δ T δ d + (C 2 + C 4 )E T δ ḋ, (D.28) δ (d)E(d) T δ ḋ. (D.29)In the above equations, C 3 and C 4 are defined asC 3 = -S(h 1 )G T + N 7 ḊrE + S( Ẇe )P 1 P + H 1 F 1 G T , (D.30) C 4 = -S(h 2 )G T + N 8 ḊrE + H 2 F 1 G T , (D.31)with Ḋ taken from equation (70) as Ḋ = E T ḋ = U T Appendix E. Linearization of the contact force vector

  f c given byf c = p N (d)G(d; ḋ); G(x; d; ḋ) = ∂g N ∂x -µ ġh T ġh T 2 + ε . (E.2)

2 c

 2 = E(δH 1 ) T R T e f c , (E.5b)t 3 c = EH T 1 (δR e ) T f c , (E.5c) t 4 c = EH T 1 R T e δf c . (E.5d)Before moving forward with the above linearizations, it is useful to provide some preliminary results. Firstly, f c is rigidly rotated back to giveF c = R T e f c . (E.6)At the same time the operator Ŝ(•) is defined, which transforms a 12 × 1 array a into a 12 × 3 matrix A as A = Ŝ(a) = the operator S as defined in equation (4). In order to obtain t 1 c , the variation δE is needed. This can be obtained from its definition (70) and equation (D.2), to give order to obtain t r c , the variation of δf c is needed. From equation (E.2) this can be obtained asδf c = δp N G + p N δG, (E.9)which in turn requires the variations δp N and δG. The variation of p N depends only on δd and is given asδp N = p N (g N )δg N = p N (g N ) ∂g N ∂x • R e H 1 E T δd , (E.10)with p N obtained from (111) p N = p N (g N ) = δG has contributions in δd and δ ḋ. This can be obtained from (E.2) as

1 F 2 c = EδH T 1 F c = E N 7 l 2 nA 1 1 F 1 F

 1212111 c )δW e = -E Ŝ(H T 1 F c )G T E T δd, (E.18a) t rδd -S(δu t l c (rδd) + EGS(δu t )F c c rδd -EGS(F c )P 1 P E T δd, (E.18b) t 3 c = EH T 1 (δR e ) T f c = EH T 1 (R e S(δW e )) T f c = -EH T 1 S(δW e )F c = EH T 1 S(F c )δW e = EH T 1 S(F c )G T E T δd, (E.18c) t 4 c = EH T 1 R T e δf c = EH T 1 R T e K fc δd + C fc δ ḋ . (E.18d)

1 F 1 S

 11 c r -EGS(F c )P 1 P E T ds, (F c )G T E T ds, e C fc ds. (E.23) Appendix F. Supplementary data

  Using (E.[START_REF] Neto | Static analysis of offshore risers with a geometrically-exact 3D beam model subjected to unilateral contact[END_REF], (E.7), (E.8), (E.15) alongside (D.3), (D.4), (D.8) into (E.5) gives S(δW e ) H T 1 F c 1:3 S(δW e ) H T 1 F c 4:6 S(δW e ) H T 1 F c 7:9 S(δW e ) H T 1 F c 10:12

							
	t 1 c = δEH T 1 F c = E	  
										N ∂x∂x	R e H 1 Ḋ ,	(E.14a)
		A 2 = -I -	∂g N ∂x	⊗	∂g N ∂x	R e S(H 1 Ḋ)G T E T ,	(E.14b)
		A 3 = I -	∂g N ∂x	⊗	∂g N ∂x		N 7 l 2
		A 4 = I -	∂g N ∂x	⊗	∂g N ∂x	R e H 1 Ŝ( Ḋ)G T E T ,	(E.14d)
		B = I -	∂g N ∂x	⊗	∂g N ∂x	R e H 1 E T .	(E.14e)
	Using (E.10), (E.12), (E.13) into (E.9) yields
								δf c = K fc δd + C fc δ ḋ,	(E.15)
	with
	K fc = p N (g N ) G ⊗ (EH T 1 R T e	∂g N ∂x	)
	+ p N	∂ 2 g N ∂x∂x	R e H 1 (s)E T -		µ ġh T	2 + ε	I -	ġh T	1 2 + ε	( ġh
		ġh T	2 + ε		I -	ġh T	1 2 + ε	( ġh

n R e A 1 Ḋr + R e S(G T Ḋ)P 1 P E T , (E.14c) T ⊗ ġh T ) A , (E.

[START_REF] Wriggers | Computational contact mechanics[END_REF]

)

C fc = -p N µ T ⊗ ġh T ) B.

(E.17)

In what follows, the superscript (e) denoting an elemental quantity, is dropped for notational convenience, letter a ∈ {1,

2} is reserved to nodal quantities and letter i ∈ {1, 2, 3} is used for defining dimensions.

The first component of Xo has been decided by the authors, as it was not detailed in[START_REF] Wang | A virtual-reality simulator and force sensation combined catheter operation training system and its preliminary evaluation[END_REF][START_REF] Wang | Dynamic simulation of frictional multi-zone contacts of thin beams[END_REF]

As the main purpose of this example is evaluating the conservation properties of the joint contact algorithm and time integrator, in this case the rigid surface is simply modeled using an analytical signed distance (trivial in the case of a flat surface).

Acknowledgments

The authors are grateful to the European Research Council for ERC grant Biolochanics, grant number 647067.

Appendix A. Kinematics assuming small displacements and rotations

Assuming small rotations in Ω, the rotation tensor R can be approximated as R( Θ) I + S( Θ) (A. [START_REF] White | Seabed characterisation and models for pipeline-soil interaction[END_REF] and its first derivative as R

S( Θ ). (A.2)

At the same time, by assuming small displacement of the centerline in Ω, the derivatives of the transversal displacements are given by u 2 (s) Θ3 (s), (A.3)

The above approximations are used to better analyze the structure of the deformation gradient F , as defined in [START_REF] Simo | A three-dimensional finite-strain rod model. part II: Computational aspects[END_REF], which is repeated for convenience,

with ā given by (see equation ( 44))

and, via equations ( 45), [START_REF] Meier | Geometrically Exact Finite Element Formulations for Slender Beams: Kirchhoff-Love Theory Versus Simo-Reissner Theory[END_REF], γ = x 0 (s) -t1 , (A.7) κα = R (s) Ēα .

(A.8)

Using equations [START_REF] Rankin | The use of projectors to improve finite element performance[END_REF], [START_REF] Crisfield | A consistent co-rotational formulation for non-linear, three-dimensional, beamelements[END_REF] alongside approximations (A.3), (A.4), x 0 can be written as

On the other hand, the displacement ū can be approximated assuming small rotations by replacing approximation (A.1) into equation [START_REF] Crisfield | A consistent co-rotational formulation for non-linear, three-dimensional, beamelements[END_REF], to yield

and the gradient can be obtained by reminding that, in the configuration Ω, X1 = s, Xα = ξ α , which yields

where, again, the small displacement assumption (A.3), (A.4) has been used. Finally, comparing the result above with equation (A.13), yields

The above gradient can be approximated using a sum of symmetric and antisymmetric tensors

)

Appendix B. Isoparametric interpolation

The Finite Element interpolation matrices are

with ξ ∈ [0, l 0 ] and the shape functions given by

Appendix C. Internal force vector, stiffness matrix

The relationship between δd and δ D requires several steps, as explained in [START_REF] Battini | Co-rotational beam elements with warping effects in instability problems[END_REF]. For the sake of completeness, the main results are summarized in what follows. Firstly, the variation of Ra requires the spin vector δ Wa , which does not appear in the variational equation [START_REF] Wang | Extended immersed boundary method using FEM and RKPM[END_REF]. Therefore, an additional set of nodal quantities is considered that inclues the variation of the spin vectors Wa , D * = ū W T Video 1: Beam against sphere. Comparison of the solutions obtained using 3 different voxel mesh discretizations.

Video 2: Contact between a cable and a rigid cylinder. Comparison of the solutions obtained using µ = 0.00 and µ = 0.30.