Next Article in Journal
Assessment of Mineralogical Characteristics of Clays and the Effect of Waste Materials on Their Index Properties for the Production of Bricks
Next Article in Special Issue
Organic Anode Materials for Lithium-Ion Batteries: Recent Progress and Challenges
Previous Article in Journal
Effects of Heat Treatment and Severe Surface Plastic Deformation on Mechanical Characteristics, Fatigue, and Wear of Cu-10Al-5Fe Bronze
Previous Article in Special Issue
On the Stability of Complex Concentrated (CC)/High Entropy (HE) Solid Solutions and the Contamination with Oxygen of Solid Solutions in Refractory Metal Intermetallic Composites (RM(Nb)ICs) and Refractory Complex Concentrated Alloys (RCCAs)
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Leonid P. Obrezkov
1,2,*,
Taija Finni
1 and
Marko K. Matikainen
2
1
Faculty of Sport and Health Sciences, University of Jyväskylä, 40014 Jyväskylä, Finland
2
Mechanical Engineering, LUT University, 53850 Lappeenranta, Finland
*
Author to whom correspondence should be addressed.
Materials 2022, 15(24), 8906; https://doi.org/10.3390/ma15248906
Submission received: 10 November 2022 / Revised: 7 December 2022 / Accepted: 9 December 2022 / Published: 13 December 2022
(This article belongs to the Special Issue Mechanics and Analysis of Advanced Materials and Structures)

Abstract

:
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.

1. 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 cross-sections, 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.

2. 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.

2.1. 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
r = r ¯ + u h ,
where u h is a displacement vector. Hence, the body motion is
r ( x , y , z , t ) = N m x , y , z q t ,
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:
r , x i = r i x , r , y i = r i y , r , z i = r i z ,
r , y y i = 2 r i y 2 , r , y z i = 2 r i y z , r , z z i = 2 r i z 2 .
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:
q i = [ r i , r y i , r z i , r y y i , r z z i , r y z i ] .
Accordingly, the vector of displacements u h has the form
u h ( x , y , z , t ) = N m x , y , z u t ,
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 1 , 1 , where ξ = 2 x l x , η = 2 y l y , ζ = 2 z l z . 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
r ( ξ , η , ζ , t ) = N m ξ , η , ζ q ( t ) , u h ( ξ , η , ζ , t ) = N m ξ , η , ζ u ( t ) .
Then, the form of the shape function matrix is
N m ( ξ , η , ζ ) = N 1 I   N 2 I   N 3 I . . . N 18 I ,
where I is a 3 × 3 identity matrix and components of N m are
N 1 = 1 2 ξ ( ξ 1 ) N 2 = 1 4 l y ξ η ( ξ 1 ) N 3 = 1 4 l z ξ ζ ( ξ 1 ) N 4 = 1 8 l z l y ξ η ζ ( ξ 1 ) N 5 = 1 16 l y 2 ξ η 2 ( ξ 1 ) N 6 = 1 16 l z 2 ξ ζ 2 ( ξ 1 ) N 7 = 1 ξ 2 N 8 = 1 2 l y η ( 1 ξ 2 ) N 9 = 1 2 l z ζ ( 1 ξ 2 ) N 10 = 1 4 l z l y η ζ ( 1 ξ 2 ) N 11 = 1 8 l y 2 η 2 ( 1 ξ 2 ) N 12 = 1 8 l z 2 ζ 2 ( 1 ξ 2 ) N 13 = 1 2 ξ ( ξ + 1 ) N 14 = 1 4 l y ξ η ( ξ + 1 ) N 15 = 1 4 l z ξ ζ ( ξ + 1 ) N 16 = 1 8 l z l y ξ η ζ ( ξ + 1 ) N 17 = 1 16 l y 2 ξ η 2 ( ξ + 1 ) N 18 = 1 16 l z 2 ξ ζ 2 ( ξ + 1 ) .
For further investigation, it is necessary to define the deformation gradient F . From (1) and (2), it can be written as
F = r r ¯ = r ξ r ¯ ξ 1 = I + u h ξ r ¯ ξ 1 .
The determinant of F defines the volume ratio of the element, we assume
J = det F > 0 .

2.2. 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,
Ω f x , y d Ω = i = 1 n j = 1 n f x i , y j w i w j ,
where 2 n 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:
V i = ( α i , β i ) , i = 1 , . . , φ ,
Ω = [ V 1 , V 2 ] [ V 2 , V 3 ] . . . [ V φ , V 1 ] .
The lines [ V i , V i + 1 ] also have several additional “control” points, such as P i 1 = V i , P i 2 , . . . , P i m i = V i + 1 , or in the binormalized coordinates as P i 1 ξ = V i ξ , . . . , P i m i ξ = V i + 1 ξ . Subsequently, the “cumulative chordal” formula parametrization is recalled:
[ α i j ξ , β i j ξ ] = 0 , j = 1 m i 1 Δ t i j , Δ t i j = P i j + 1 ξ P i j ξ , j = 1 , . . . , m i 1 .
Then, each line [ V i ξ , V i + 1 ξ ] is tracked by a spline curve S i ( t ) = ( S i 1 ( t ) , S i 2 ( t ) ) degree of p i , where p i m i 1 , see Figure 2.
Then the cubature formula with the 2 n 1 polynomial exactness degree over the Ω domain has the form
I 2 n 1 = λ Λ 2 n 1 f ( η λ , ζ λ ) w λ ,
where
Λ 2 n 1 = { λ = ( i , j , k , h ) : 1 i φ , 1 j m i 1 ,
1 k n i , 1 h n } ,
and w λ , η λ and ζ λ are:
η λ = S i 1 ( q i j k ) Ξ 2 τ h n + S i 1 ( q i j k ) + Ξ 2 ,
ζ λ = S i 2 ( q i j k ) ,
w λ = Δ t i j 4 w k n i w h n ( S i 1 ( q i j k ) Ξ ) d S i 2 ( t ) d t t = q i j k ,
q i j k = Δ t i j 2 τ k n i + t i j + 1 + t i j 2 , Δ t i j = t i j + 1 t i j ,
n i = n p i + p i / 2 , p i is even , n p i + ( p i + 1 ) / 2 , p i is odd .
Thus, only τ k n i , w k n i and Ξ need to be defined. Ξ is an arbitrary straight line
Ω R 2 = [ a , b ] × [ c , d ] , Ξ ( η ) [ a , b ] , η [ c , d ] .
The choice of Ξ does not have any influence. However, it is necessary to obtain the nodes and weights. τ k n i , w k n i are the nodes and weights, respectively, of the Gauss–Legendre quadrature formula of the exactness degree 2 n i 1 on [ 1 , 1 ] .

3. 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:
δ Π e x t δ Π i n t + δ Π i n e r t δ Π c o n = 0 .
δ Π i n e r t can be written as
δ Π i n e r t = q ¨ T V ρ N T N d V · δ q ,
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 N d V . In the case of the static problem, which is the concern of this work, δ Π i n e r t = 0 . The variation of Π i n t with respect to the nodal coordinates is [16]
δ Π i n t = V S : δ E d V = V S : E q d V · δ q .
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
E = 1 2 F T · F I .
The last is the variation of the contact force work δ Π c o n , which will be explained.

4. 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
F = J 1 3 F ¯ , J = det F > 0 .
This leads to the follow representation of the right Cauchy–Green tensor:
C ¯ = F ¯ T · F ¯ .
Thus, after the decomposition, we have
Ψ = Ψ v o l ( J ) + Ψ i s o ( C ¯ , ) ,
where Ψ v o l ( J ) = k ( J 1 ) 2 , k is a penalty coefficient to guarantee the incompressibility. The part Ψ i s o might be reformulated in the terms of the Cauchy–Green deformation tensor invariants,
Ψ i s o = Ψ i s o ( I 1 ¯ , I 2 ¯ ) ,
where I 1 ¯ and I 2 ¯ have the forms
I ¯ 1 = tr C ¯ , I ¯ 2 = 1 2 tr C ¯ 2 + tr 2 C ¯ .
The second Piola–Kirchhoff stress from (14) is formulated as follows:
S = 2 Ψ C = 2 Ψ C ¯ : C ¯ C .
Using (18), it can be expressed as
S = 2 Ψ ¯ C ¯ C ¯ C + 2 Ψ v o l J J C = 2 k Ψ ¯ I ¯ k I ¯ k C ¯ C ¯ C + Ψ v o l J J C 1 ,
J C = 1 2 J C 1 .
The corresponding volumetric part has the form
S v o l = d J 1 J C 1 .
In this study, we consider one type of material model: the neo-Hookean model. The isochoric part of the neo-Hookean model is
Ψ ¯ = c 10 I 1 ¯ 3 ,
with the expression for the second Piola–Kirchhoff stress tensor:
S ¯ = 2 c 10 J 2 3 I 1 3 I ¯ 1 C ¯ 1 .

5. 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:
d = r A r B .
The nonpenetration condition is defined via the so-called gap function, which in this work is given as follows,
g ( ξ A , ξ B , η A , η B , ζ A , ζ B ) = r A ( ξ A ) η , ζ = 0 r B ( ξ c B ) η , ζ = 0 ( r A ( ξ A , η A , ζ A ) r A ( ξ A ) η , ζ = 0 + r c B ξ c B , η c B , ζ c B r B ( ξ c B ) η , ζ = 0 ) ,
where g ( ξ A , ξ B , η A , η B , ζ A , ζ B ) 0 . The subscript c denotes the orthogonal projection of the point on beam A on the beam B obtained from (27), r B ( ξ c B ) η , ζ = 0 is the point projection r c B ξ c B , η c B , ζ c B 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,
δ Π c o n = p n Ω c g δ g d Ω ,
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,
δ Π con = δ u A T p n i = 1 n i j = 1 n j k = 1 n k g ξ i A , η j A , ζ k A N A T n i j k w i w j w k + δ u B T p n i = 1 n i j = 1 n j k = 1 n k g ξ c B , η c B , ζ c B N B T n i j k w i w j w k ,
where
n i j k = n ξ c B ( ξ i A ) , η c B ( ξ i A , η j A , ζ k A ) , η c B ( ξ i A , η j A , ζ k A ) ,
N A T = N T ( ξ i A , η j A , ζ k A )
N B T = N T ξ c B ( ξ i A ) , η c B ( ξ i A , η j A , ζ k A ) , ζ c B ( ξ i A , η j A , ζ k A ) .
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. ξ c B , η c B and ζ c B are the parameters of the closest projected point r ( ξ j A , η k A , ζ k A ) on B, n is a normal vector from the B to A beam elements’ surfaces.

6. 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.

7. 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.

Author Contributions

Conceptualization, M.K.M. and T.F.; methodology, M.K.M. and L.P.O.; software, Matlab; validation, M.K.M. and T.F., and L.P.O.; formal analysis, L.P.O.; investigation, L.P.O.; writing—original draft preparation, L.P.O.; writing—review and editing, M.K.M. and T.F. and L.P.O.; visualization, L.P.O.; supervision, M.K.M. and T.F.; project administration, M.K.M. and T.F.; funding acquisition, T.F. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by The Academy of Finland (Decisions No. 299033 and 323168).

Data Availability Statement

Not applicable.

Acknowledgments

We would like to thank the Academy of Finland (Decisions No. 299033 and 323168) for funding.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
ATAchilles tendon
FEMFinite Element Model
ANCFAbsolute Nodal Coordinate Formulation
DOFDegrees of Freedom

References

  1. Komi, P.V.; Fukashiro, S.; Järvinen, M. Biomechanical loading of Achilles tendon during normal locomotion. Clin. J. Sport Med. 1992, 11, 521–531. [Google Scholar] [CrossRef]
  2. Järvinen, T.A.; Kannus, P.; Maffulli, N.; Khan, K.M. Achilles Tendon Disorders: Etiology and Epidemiology. Foot Ankle Clin. 2005, 10, 255–266. [Google Scholar] [CrossRef] [PubMed]
  3. Slane, L.C.; Thelen, D.G. Non-uniform displacements within the Achilles tendon observed during passive and eccentric loading. J. Biomech. 2014, 47, 2831–2835. [Google Scholar] [CrossRef] [Green Version]
  4. Khair, R.M.; Stenroth, L.; Péter, A.; Cronin, N.J.; Reito, A.; Paloneva, J.; Finni, T. Non-uniform displacement within ruptured Achilles tendon during isometric contraction. Scand. J. Med. Sci. Sport. 2021, 31, 1069–1077. [Google Scholar] [CrossRef] [PubMed]
  5. Obuchowicz, R.; Ekiert, M.; Kohut, P.; Holak, K.; Ambrozinski, L.; Tomaszewski, K.; Uhl, T.; Mlyniec, A. Interfascicular matrix-mediated transverse deformation and sliding of discontinuous tendon subcomponents control the viscoelasticity and failure of tendons. J. Mech. Behav. Biomed. Mater. 2019, 97, 238–246. [Google Scholar] [CrossRef] [PubMed]
  6. Farris, D.J.; Trewartha, G.; McGuigan, M.P.; Lichtwark, G.A. Differential strain patterns of the human Achilles tendon determined in vivo with freehand three-dimensional ultrasound imaging. J. Exp. Biol. 2013, 216, 594–600. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Obst, S.J.; Renault, J.B.; Newsham-West, R.; Barrett, R.S. Three-dimensional deformation and transverse rotation of the human free Achilles tendon in vivo during isometric plantarflexion contraction. J. Appl. Physiol. 2014, 116, 376–384. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  8. Bojsen-Møller, J.; Magnusson, S.S. Heterogeneous Loading of the Human Achilles Tendon In Vivo. Exerc. Sport Sci. Rev. 2015, 43, 190–197. [Google Scholar] [CrossRef] [Green Version]
  9. Shim, V.; Handsfield, G.; Fernandez, J.; Lloyd, D.; Besier, T. Combining in silico and in vitro experiments to characterize the role of fascicle twist in the Achilles tendon. Sci. Rep. 2018, 8, 13856. [Google Scholar] [CrossRef] [Green Version]
  10. Hansen, W.; Shim, V.B.; Obst, S.; Lloyd, D.G.; Newsham-West, R.; Barrett, R.S. Achilles tendon stress is more sensitive to subject-specific geometry than subject-specific material properties: A finite element analysis. J. Biomech. 2017, 56, 26–31. [Google Scholar] [CrossRef]
  11. Kinugasa, R.; Yamamura, N.; Sinha, S.; Takagi, S. Influence of intramuscular fiber orientation on the Achilles tendon curvature using three-dimensional finite element modeling of contracting skeletal muscle. J. Biomech. 2016, 49, 3592–3595. [Google Scholar] [CrossRef] [PubMed]
  12. Morales-Orcajo, E.; Souza, T.R.; Bayod, J.; Barbosa de Las Casas, E. Non-linear finite element model to assess the effect of tendon forces on the foot-ankle complex. Med. Eng. Phys. 2017, 49, 71–78. [Google Scholar] [CrossRef] [PubMed]
  13. taş, R.A.; Lucaciu, D.O. Finite Element Analysis of the Achilles Tendon While Running. Acta Med. Marisiensis 2013, 59, 8–11. [Google Scholar] [CrossRef]
  14. Bozorgmehri, B.; Obrezkov, L.P.; Harish, A.B.; Matikainen, M.K.; Mikkola, A. A contact description for continuum beams with deformable arbitrary cross-section. Finite Elem. Anal. Des. 2023, 214, 103863. [Google Scholar] [CrossRef]
  15. Obrezkov, L.; Bozorgmehri, B.; Finni, T.; Matikainen, M.K. Approximation of pre-twisted Achilles sub-tendons with continuum-based beam elements. Appl. Math. Model. 2022, 112, 669–689. [Google Scholar] [CrossRef]
  16. Obrezkov, L.P.; Matikainen, M.K.; Harish, A.B. A finite element for soft tissue deformation based on the absolute nodal coordinate formulation. Acta Mech. 2020, 231, 1519–1538. [Google Scholar] [CrossRef]
  17. Obrezkov, L.P.; Eliasson, P.; Harish, A.B.; Matikainen, M.K. Usability of finite elements based on the absolute nodal coordinate formulation for the Achilles tendon modelling. Int. J. Non-Linear Mech. 2021, 129, 103662. [Google Scholar] [CrossRef]
  18. Weiss, J.A.; Maker, B.N.; Govindjee, S. Finite element implementation of incompressible, transversely isotropic hyperelasticity. Comput. Methods Appl. Mech. Eng. 1996, 135, 107–128. [Google Scholar] [CrossRef]
  19. Annaidh, A.N.; Destrade, M.; Gilchrist, M.D.; Murphy, J.G. Deficiencies in numerical models of anisotropic nonlinearly elastic materials. Biomech. Model. Mechanobiol. 2013, 12, 781–791. [Google Scholar] [CrossRef]
  20. Escalona, J.L.; Hussien, H.A.; Shabana, A.A. Application of absolute nodal co-ordinate formulation to multibody system dynamics. J. Sound Vib. 1998, 214, 833–851. [Google Scholar] [CrossRef]
  21. Maqueda, L.G.; Bauchau, O.A.; Shabana, A.A. Effect of the centrifugal forces on the finite element eigenvalue solution of a rotating blade: A comparative study. Multibody Syst. Dyn. 2008, 19, 281–302. [Google Scholar] [CrossRef]
  22. Shen, Z.; Li, P.; Liu, C.; Hu, G. A finite element beam model including cross-section distortion in the absolute nodal coordinate formulation. Nonlinear Dyn. 2014, 77, 1019–1033. [Google Scholar] [CrossRef]
  23. Nachbagauer, K.; Gruber, P.; Gerstmayr, J. A 3D Shear Deformable finite element based on the absolute nodal coordinate formulation. Multibody Dyn. 2013, 28, 77–96. [Google Scholar] [CrossRef]
  24. Obrezkov, L.P.; Mikkola, A.; Matikainen, M.K. Performance review of locking alleviation methods for continuum ANCF beam elements. Nonlinear Dyn. 2022, 109, 531–546. [Google Scholar] [CrossRef]
  25. Patel, M.; Shabana, A.A. Locking alleviation in the large displacement analysis of beam elements: The strain split method. Acta Mech. 2018, 229, 2923–2946. [Google Scholar] [CrossRef]
  26. Ebel, H.; Matikainen, M.K.; Hurskainen, V.V.; Mikkola, A. Higher-order beam elements based on the absolut nodal coordinate formulation for three-dimensional elasticity. Nonlinear Dyn. 2017, 88, 1075–1091. [Google Scholar] [CrossRef]
  27. Abramowitz, M.; Stegun, I.A.; Romer, R.H. Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables. Am. J. Phys. 1988, 58, 958. [Google Scholar] [CrossRef] [Green Version]
  28. Handsfield, G.G.; Greiner, J.; Madl, J.; Rog-Zielinska, E.A.; Hollville, E.; Vanwanseele, B.; Shim, V. Achilles Subtendon Structure and Behavior as Evidenced From Tendon Imaging and Computation Modeling. Front. Sport. Act. Living 2020, 2, 70. [Google Scholar] [CrossRef]
  29. Yin, N.Y.; Fromme, P.; McCarthy, I.; Birch, H. Individual variation in Achilles tendon morphology and geometry changes susceptibility to injury. eLife 2021, 10, e63204. [Google Scholar] [CrossRef]
  30. Edama, M.; Kubo, M.; Onishi, H.; Takabayashi, T.; Inai, T.; Yokoyama, E.; Hiroshi, W.; Satoshi, N.; Kageyama, I. The twisted structure of the human Achilles tendon. Scand. J. Med. Sci. Sport. 2015, 25, e497–e503. [Google Scholar] [CrossRef]
  31. Finni, T.; Bernabei, M.; Baan, G.C.; Noort, W.; Tijs, C.; Maas, H. Non-uniform displacement and strain between the soleus and gastrocnemius subtendons of rat Achilles tendon. Scand. J. Med. Sci. Sport. 2018, 28, 1009–1017. [Google Scholar] [CrossRef] [PubMed]
Figure 1. 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].
Figure 1. 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].
Materials 15 08906 g001
Figure 2. An arbitrary domain in initial and local coordinate systems.
Figure 2. An arbitrary domain in initial and local coordinate systems.
Materials 15 08906 g002
Figure 3. The Type III tendon representation. (a) The Achilles sub-tendons’ cross sections. (b) The pretwisted underformed Achilles tendon discretized by four ANCF-based continuum beam elements at each subtendon with ψ = 45°.
Figure 3. The Type III tendon representation. (a) The Achilles sub-tendons’ cross sections. (b) The pretwisted underformed Achilles tendon discretized by four ANCF-based continuum beam elements at each subtendon with ψ = 45°.
Materials 15 08906 g003
Figure 4. Integration approximations of the three subtendons by the Gauss–Green cubature formula. (a) Soleus. (b) Medial gastrocnemius. (c) Lateral gastrocnemius.
Figure 4. Integration approximations of the three subtendons by the Gauss–Green cubature formula. (a) Soleus. (b) Medial gastrocnemius. (c) Lateral gastrocnemius.
Materials 15 08906 g004
Figure 5. The deformed shapes of the pre-twisted Achilles tendons when the soleus is loaded and discretized by four ANCF-based continuum beam elements at each subtendon. (a) ψ = 15 (b) ψ = 45 .
Figure 5. The deformed shapes of the pre-twisted Achilles tendons when the soleus is loaded and discretized by four ANCF-based continuum beam elements at each subtendon. (a) ψ = 15 (b) ψ = 45 .
Materials 15 08906 g005
Table 1. The elongation test results in [mm] for the straight and pretwisted soleus subtendon of Type III tendons from the neo-Hookean material model.
Table 1. The elongation test results in [mm] for the straight and pretwisted soleus subtendon of Type III tendons from the neo-Hookean material model.
Elongation [mm] of the Soleus Sub-Tendon
Applied LoadVariation of ψ
[N] ψ = 0 ψ = 15 ψ = 45
100.1430.1460.168
200.2860.2890.313
300.4300.4330.457
400.5740.5780.602
450.6470.6500.675
600.8640.8680.893
801.1571.1601.186
901.3041.3071.333
1001.4511.4541.481
1502.1972.2012.228
2002.9582.9612.989
3004.5244.5274.557
4006.1516.1556.187
Table 2. Elongation results in [mm] for several mesh refinements for the straight and pretwisted soleus sub-tendon of Type III from the neo-Hookean material model under N = 400 N tensile force.
Table 2. Elongation results in [mm] for several mesh refinements for the straight and pretwisted soleus sub-tendon of Type III from the neo-Hookean material model under N = 400 N tensile force.
Elongation [mm] of the Soleus Sub-Tendon
Element NumberVariation of ψ
per Sub-Tendons ψ = 0 ψ = 15 ψ = 45
n Sol × n MG × n LG
1 × 1 × 1 6.0516.0556.123
2 × 2 × 2 6.0896.0936.123
4 × 4 × 4 6.1516.1556.187
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Obrezkov, L.P.; Finni, T.; Matikainen, M.K. Modeling of the Achilles Subtendons and Their Interactions in a Framework of the Absolute Nodal Coordinate Formulation. Materials 2022, 15, 8906. https://doi.org/10.3390/ma15248906

AMA Style

Obrezkov LP, Finni T, Matikainen MK. Modeling of the Achilles Subtendons and Their Interactions in a Framework of the Absolute Nodal Coordinate Formulation. Materials. 2022; 15(24):8906. https://doi.org/10.3390/ma15248906

Chicago/Turabian Style

Obrezkov, Leonid P., Taija Finni, and Marko K. Matikainen. 2022. "Modeling of the Achilles Subtendons and Their Interactions in a Framework of the Absolute Nodal Coordinate Formulation" Materials 15, no. 24: 8906. https://doi.org/10.3390/ma15248906

APA Style

Obrezkov, L. P., Finni, T., & Matikainen, M. K. (2022). Modeling of the Achilles Subtendons and Their Interactions in a Framework of the Absolute Nodal Coordinate Formulation. Materials, 15(24), 8906. https://doi.org/10.3390/ma15248906

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop