Modeling of the Achilles Subtendons and Their Interactions in a Framework of the Absolute Nodal Coordinate Formulation

Experimental results have revealed the sophisticated Achilles tendon (AT) structure, including its material properties and complex geometry. The latter incorporates a twisted design and composite construction consisting of three subtendons. Each of them has a nonstandard cross-section. All these factors make the AT deformation analysis computationally demanding. Generally, 3D finite solid elements are used to develop models for AT because they can discretize almost any shape, providing reliable results. However, they also require dense discretization in all three dimensions, leading to a high computational cost. One way to reduce degrees of freedom is the utilization of finite beam elements, requiring only line discretization over the length of subtendons. However, using the material models known from continuum mechanics is challenging because these elements do not usually have 3D elasticity in their descriptions. Furthermore, the contact is defined at the beam axis instead of using a more general surface-to-surface formulation. This work studies the continuum beam elements based on the absolute nodal coordinate formulation (ANCF) for AT modeling. ANCF beam elements require discretization only in one direction, making the model less computationally expensive. Recent work demonstrates that these elements can describe various cross-sections and materials models, thus allowing the approximation of AT complexity. In this study, the tendon model is reproduced by the ANCF continuum beam elements using the isotropic incompressible model to present material features.


Introduction
The Achilles tendon (AT) is the strongest tendon in the body and serves an important function during locomotion. It can reach loads up to four times the body weight while walking and approximately 10 times while running, with the upper border as much as 9 kN [1]. At the same time, it is vulnerable to traumatic injuries due to chronic or acute overloading [2], with the determinants of good recovery not well understood. The AT possesses a complex structure whereby three subtendons, having subject-specific crosssections, arise from the soleus and the lateral and medial heads of the gastrocnemius muscles. The three heads twist around each other, counterclockwise in the right AT and clockwise in the left. With the complex structure comes functional consequences. Studies have revealed nonuniform displacements within healthy AT [3], whereas the displacement can be more uniform in an injured tendon [4]. Furthermore, loading from the three different muscles causes nonhomogeneous longitudinal strains, compression, and transverse strains in the AT [5]. Longitudinal and transverse strains have also been reported in human studies [6,7]. These nonhomogeneous strains are likely linked to the architectural structure of the tendon [8]. The AT can endure the large forces transmitted axially, and they are most often studied, whereas less attention is placed on shear forces. Because of its significant role in the human musculoskeletal system, a better understanding of AT can provide valuable information for diagnosis and treatment.
What is the significance of the complex geometry and twist of AT? Previous research has shown that region-specific susceptibility to strain injury changes with the amount of tendon twist [9] and that the changes in AT stress are more sensitive to volumetric tendon shape rather than material properties [10]. The appropriate model can help study stress and strain distributions within the AT and improve understanding of tendon function in health, disease, and rehabilitation.
One of the most popular methods in modern biomechanical research for studying tendons is the finite element method (FEM). Using finite solid elements in a framework of the nonlinear FEM helps to comprehend the tendon's complex geometry and materiality. However, AT modeling with the solid finite elements leads to a significant number of degrees of freedom (DOFs) [11][12][13], which results in a long computation time to obtain a solution. Solid elements can approximate almost any shape and provide reliable results, but require dense discretization in all three dimensions. Therefore, other types of finite elements are necessary to decrease the computational cost.
This research introduces a new approach to the deformation analysis of the Achilles subtendons and their interactions. The main idea is to use one-dimensional finite element discretization over the subtendon's length (in the longitudinal direction) to decrease computational costs. It can be achieved by considering the tendon as a beam-like structure using ANCF-based continuum beam elements with specific descriptions for geometrically complicated deformable cross-sections. This approach leads to the finite element discretization over the subtendon's length and makes it possible to consider material laws based on the continuum mechanics. Recent studies [14][15][16][17] show that this transition is possible without significant losses in the quality of the results within the ANCF framework. For example, the work [16] considers the deformations of the beam-like structures described with the ANCF continuum beam elements and from the various soft material models. In [17], the approximation of the rat Achilles tendon experiment with the ANCF elements is given and verified against experimental results. Ref. [15] provides the approximation way for arbitrary cross-sections, which is suitable for the continuum-based ANCF beams. For example, the provided technique approximates one of the Achilles subtendons. In the case of multibeam construction, the question of mutual interaction between beams arises, i.e., the so-called contact problem. In the work [14], the methods for solving contact problems between beams with arbitrary cross-sections are presented.
In this study, we explored the method given in [15] and modeled the whole AT with the ANCF continuum beam elements. The subtendons' cross-sections are obtained with the integration scheme proposed in [15]. Then, the obtained beams are pretwisted, one around the other. In this study, the neo-Hookean material model describes the soft tissue response [10,18]. There are also possibilities to use others' material models in a way given in [16]. However, Annaidh et al. [19] demonstrated that using anisotropic material models within FEM can have inconsistencies. The possible contact between subtendons can be described via the surface-to-surface procedure, thus, taking into account the complicated border interactions between two bodies.

ANCF Beam Element
This section provides the geometrical setup for the continuum-based ANCF beam element. The idea behind this element type is to use the slope vectors for defining the cross-section orientation and deformation. The advantages of these finite elements are discussed in [20][21][22].
There are various types of ANCF elements, and one can divide them into several groups and subgroups; for more details, the reader is referred to Nachbagauer et al. [23], Obrezkov et al. [24], Patel and Shabana [25]. Here, the higher-order three-nodded element with the second-order interpolation in longitudinal and thickness directions denoted 3363 is used. It does not require any modifications to demonstrate good performance even for complicated loading cases [26] and allows the use of all material laws based on the 3D elasticity.

Kinematics of the ANCF Continuum Beam Elements
Let r = r(x, y, z) ∈ R 3 be the position vector field of any particle in the current configuration. The position in the initial configuration is denoted as r [15], see Figure 1. The connection between the two vectors is where u h is a displacement vector. Hence, the body motion is where N m is a shape function matrix and q is a vector of nodal coordinates. q contains the position of the nodes as well as their derivatives. Therefore, we accept the following notation for i th node: Illustration of a three-nodded beam element with an arbitrary particle p in current and p in reference configurations. The three nodes are denoted by r (i) and r (i) , respectively, i = 1, 2, 3 [17].
As mentioned above, in this work, we consider 3363 beam elements [26]. The vectors of nodal coordinates related to this element are presented as follows: Accordingly, the vector of displacements u h has the form where u is a vector of nodal displacements. The element is isoparametric. Here, we introduce a new local coordinate system ξ = {ξ, η, ζ} with the range for the local coordinates Here, l x , l y and l z are the physical dimensions of the element. The substitutions are made to deal with the Gaussian integration procedure [15]. Now, we have Then, the form of the shape function matrix is where I is a 3 × 3 identity matrix and components of N m are For further investigation, it is necessary to define the deformation gradient F. From (1) and (2), it can be written as The determinant of F defines the volume ratio of the element, we assume

Cross-Section Geometry Description
The standard Gaussian quadrature formula for the integration of any function f (x, y) in the general form can be written as follows, where 2n − 1 is the polynomial exactness degree of function f over one of the axis lines, and w is the weight of the point. For simple cross-sections (circular, etc.), we send our readers to [27], where weights and points in a binormalized coordinate system can be found. Below we present the method for more complicated domains, which can also be found in [15]. Let us consider a closed domain Ω, which has a piecewise border ∂Ω with points V i on it: . Subsequently, the "cumulative chordal" formula parametrization is recalled: Then the cubature formula with the 2n − 1 polynomial exactness degree over the Ω domain has the form where and w λ , η λ and ζ λ are: Thus, only τ n i k , w n i k and Ξ need to be defined. Ξ is an arbitrary straight line The choice of Ξ does not have any influence. However, it is necessary to obtain the nodes and weights. τ

Equilibrium Equation
Our task involves many subroutines, each of them contributing to the energy balance and equilibrium of the whole system. The common approach for calculating is to use the variational formulation. The variations can be grouped as inertia, external, contact, and internal: δΠ inert can be written as where ρ is the mass density, and V is the volume of the element in the reference configuration. The mass matrix is M = V ρN T NdV. In the case of the static problem, which is the concern of this work, δΠ inert = 0. The variation of Π int with respect to the nodal coordinates is [16] S is the second Piola-Kirchhoff stress tensor, and its form depends on the material model, which will be presented in Section 4. E is the Green-Lagrange strain tensor The last is the variation of the contact force work δΠ con , which will be explained.

Approximation of the Tendon Tissue
The elastic properties of the Achilles tendon tissue can be presented in different ways. One can find examples in [10,12,13,28,29], where the Helmholtz free energy function Ψ describes elastic features for such material models. In the isotropic case, Ψ depends only on the right Cauchy-Green tensor C = F T · F, therefore, Ψ = Ψ(C). In the case of anisotropy, the additional structural tensor A can be added to define the preferable deformation direction Ψ = Ψ(C, A). Models describing AT are usually incompressible. The common approach to deal with it is to split the deformation gradient F into dilational (volumetric) and distortion (isochoric) parts. Here again, we want to send our readers to the work [19], where the authors point out the possible problems associated with the decomposition of the anisotropic material models. We have This leads to the follow representation of the right Cauchy-Green tensor: Thus, after the decomposition, we have where Ψ vol (J) = k(J − 1) 2 , k is a penalty coefficient to guarantee the incompressibility. The part Ψ iso might be reformulated in the terms of the Cauchy-Green deformation tensor invariants, where I 1 and I 2 have the forms The second Piola-Kirchhoff stress from (14) is formulated as follows: Using (18), it can be expressed as ∂J ∂C The corresponding volumetric part has the form In this study, we consider one type of material model: the neo-Hookean model. The isochoric part of the neo-Hookean model is with the expression for the second Piola-Kirchhoff stress tensor:

Contact Formulation
Working with an assembled structure consisting of two or more bodies, the question of interaction between the substructures appears. That problem requires the solution of a contact task. In this study, we are concerned with the description of the bodies of nonstandard forms, such as only surface-to-surface contact formulation, which can describe this contact [14]. Let us describe the task of two contacting beams (denoted as A and B) in the terms of the distances between the two closest position vector fields r A and r B . Then, assuming that along the contact surface there is no penetration, the minimum distance problem in the most general case can be formulated as follows: The nonpenetration condition is defined via the so-called gap function, which in this work is given as follows, The subscript c denotes the orthogonal projection of the point on beam A on the beam B obtained from (27), r B (ξ B c ) η,ζ=0 is the point projection r B c ξ B c , η B c , ζ B c on the beam centerline. In the model, we assume that two bodies are closely placed to each other, and only sliding is allowed. Therefore, the nonpenetration condition is g = 0. Then the variation reads as follows, where Ω c is the contacting surface between A and B beams, and p n is the penalty parameter. The weak form of contact energy (28) presented in Section 3 can be expressed in the discrete form as follows, where In (29), n i is the amount of Gauss points in the A beam, along the ξ direction, w i are their corresponding weight, η k and ζ k are the Gauss points coordinates along the η and ζ directions parameters. ξ B c , η B c and ζ B c are the parameters of the closest projected point r(ξ A j , η A k , ζ A k ) on B, n is a normal vector from the B to A beam elements' surfaces.

Numerical Examples
Previous studies found that the Achilles tendon consists of three subtendons with each having a complicated cross-section shape [30,31]. Additionally, there are three common types of AT with varying subtendon regions and torsion [30]. In this work, we consider the AT of Type III due to its relatively simple cross-section form. We extracted the geometrical description of subtendons from [30]. Although the exact geometrical data are not presented, we use CAD software to obtain the positions of the points, as in [15]. Here, we also considered the pretwist of the tendon about the centroidal axis (line, where all three subtendons are connected) from 0 • at x = 0 to ψ degrees at x = L. The centroidal axis of the beam remains straight. See Figure 3. The representations of the Gauss points for all three subtendons are given in Figure 4a-c.
The length of the tendon was set at L = 0.07 m [29]. The geometrical results based on the approximation are 16.31 mm 2 for the soleus subtendon, 15.98 mm 2 for the medial and 19.57 mm 2 for lateral subtendons, with total area equaling to 51.86 mm 2 . That slightly exceeds the average female tendon cross-section 51.2 mm 2 and is smaller than the average male cross-section 62.1 mm 2 [29].   We used the neo-Hookean material model with three shear modulus equal to c 10 = 103.1 MPa for soleus, c 10 = 143.2 MPa, and c 10 = 226.7 MPa for medial and lateral subtendons, respectively [29]. We considered three different pretwisted designs: ψ = 0, ψ = 15, and ψ = 45. The choice is based on the work [28], where the optimal value of twisted is found between 15 and 45 degrees. Then, the soleus subtendon was subjected to forces along the longest direction and applied at the last node, the maximum applied tensile load is 400 N. The applied force exceeds four times the loading conditions given in [15,29], allowing the demonstration of the nonlinear deformations, about 10% of the initial length. On the other edge, r = 0 from (3) is fixed at the first node, and this condition forbids the displacement, but allows the cross-sectional contraction.
The results presented in Table 1 are consistent with the ones given in [15], where the pretwisted subtendons show higher elongations under the same load in comparison to straight subtendons. The elongations for other subtendons are near zero, which indicates that there is sliding between the subtendons as in Section 5 holds. Table 2 presents the converge tests, wherein the elongation results for a number of mesh refinements for the straight and pretwisted soleus subtendon of Type III from the neo-Hookean material model subjected to N = 400 N tensile force are given. The deformed shapes for straight pretwisted ψ = 15 the tendons are given in Figure 5a, where the shapes are discretized by four ANCF-based continuum beam elements at each subtendon.

Conclusions
This work uses the continuum-based ANCF beam element to describe the human Achilles tendon's deformation due to elongation. In the study, the AT is presented as a combination of three substructures, pretwisted and sliding one around the others. The contact between them is described with the segment-to-segment algorithm. The Gauss-Green cubature integration formula captures the sophisticated cross-section form of each subtendon. The neo-Hookean isotropic material model describes the pure elastic response. The results show that the model is feasible, but more careful verification is necessary. That can include models built with conventional 3D solid elements and the comparison with experimental data.
Additionally, the work possesses certain limitations. For example, the cross-sectional area is taken to be the same for all subtendons along their longitudinal axes. That is a substantial simplification, but there is no available geometrical data to approximate such variation.