Efficient computing of the viscoelastic response of helical tendon subunits

Abstract

BACKGROUND AND OBJECTIVE:

Tendons are hierarchical structures, with a viscoelastic, time-dependent mechanical response upon axial loading. Their mechanical behavior strongly depends on the inner composition and on the material attributes of their subunits. With tens of thousands of damaged tendons each year all over the world, an increasing need for robust simulation tools arises, which can assist tendon diagnosis and restoration praxis in the selection of appropriate treatments and substitute materials.

METHODS:

In this paper, we elaborate a numerical model for the computation of the viscoelastic response of tendon fascicles to tensile loads. The model can inherently describe the composite inner tendon structure, accounting for its helical geometric arrangement and for different relaxation behaviors among its constituents. The model’s weak form and finite element implementation are derived making use of the generalized Maxwell-Wiechert model, which inherently accounts for a complex multi-parameter relaxation behavior.

RESULTS:

The model results have been validated with respect to existing homogenization results, with the proposed framework to reproduce homogenization results with an accuracy of 5% over the entire relaxation process. The simulation scheme has been further used to study the viscoelastic behavior of tendon fascicles with different inner structural compositions, namely with different physiologically relevant fiber contents, viscosity values and helical angles, extending the applicability of existent schemes. What is more, its robust formulation has been shown to outperform other simulation methods.

CONCLUSIONS:

Overall, the proposed framework has been shown to constitute a robust simulation tool for the computation of the relaxation behavior of tendon fascicles. It has provided a means to compute the time-evolution of the axial and of the torsional modulus of tendon subunits. The latter has been shown to be significant for typical fascicle structures, with its viscoelastic behavior to well compare with the one recorded for the axial modulus of the tendon’s subunit.

Keywords

Tendon fibers matrix viscoelasticity relaxation efficient computing

1. Introduction

Helical geometries are structural patterns that are extensively encountered in biological members [1]. Tendon sub-units are architected upon distinct inner helical patterns [2, 3], as schematically depicted in Fig. 1. Some of the tendons’ primal mechanical attributes, such as its considerable energy storing capacity or its coupled axial and torsional behavior upon loading have been shown to be directly related to the helical arrangement of its inner sub-units [4].

Figure 1.

Schematic representation of the tendon’s inner structure. Tendon fibrils are encased in fibers which in turn constitute the primal subunits of fascicles forming the tendon unit. Each inner fibrillar structure is immersed in a matrix substance [2].

The geometric configuration of helical structures determines to a large extent their mechanical response [5]. A series of studies have been devoted to the derivation of analytical expressions to characterize their mechanical response upon different loading patterns, such as to axial, torsional, thermal or bending loads [6, 7, 8, 9, 10]. Formulas of the kind can be handily used for small scale calculations, but have been primarily restrained to the elastic material analysis range [11, 12, 13].

Analytical modeling approaches have been complemented by a wide range of numerical schemes, mainly developed for the simulation of helical cable structures [14, 15, 16]. A large part of the numerical schemes presented has been based on finite volume elements, simulating a period evolution of the helical domain. Volume element based models come along with considerable numerical analysis costs, already for a linear elastic analysis [17, 18, 19, 20]. However, incorporating the structural symmetry of the domain so as to simulate the volumetric behaviour through its two-dimensional cross sectional domain has allowed for a substantial computational cost reduction [21, 22]. Note that the simulation of the helix effective volumetric behavior is a prerequisite for its characteristic axial-torsional response to be captured; an attribute that cannot be reproduced by one-dimensional approaches [23] or analytical schemes [24].

Tendons are viscoelastic structures with multiscale stress absorption and relaxation capabilities [25], which are developed over the entire range of their inner hierarchical scales (Fig. 1) [26]. Indicatively, at the tendon’s fascicle scale, relaxation times between 250 and 500 seconds relating to effective viscosity values $\eta$ greater than 125 GPa s have been reported in independent experimental studies [27, 28, 29]. Non-negligible viscoelastic moduli have been accordingly provided for the inner scales of the fiber and fibril (Fig. 1). In particular, already at the fibril scale, relaxation times $\tau$ of several seconds, along with viscoelastic moduli $\eta$ values up to 13 GPa s have been reported [30] and approximated by both single relaxation time [30] and two component Wiechert type relaxation elements [31]. Accordingly, at the fiber scale, single relaxation times lower than a minute and up to several minutes have been experimentally obtained [31, 32]. The elastic fiber modulus has been observed to vary up to an order of magnitude, with reported wet fiber moduli values as low as 200 MPa and up to 1.5 GPa [33]. It has been noted that aging significantly reduces the tendon viscosity at all its inner scales, while it minorly affects its elastic properties [34].

Up to now, different numerical models have been presented for the simulation of fibrous and porous media [35, 36, 37], with a rather limited number of numerical approaches to explicitly address tendon mechanics [20, 38, 39]. Such models are of primal importance, not only for the understanding of the tendon’s inner composition and mechanical attributes, but also for the design of artificial restoration materials [40, 41, 42] and for their use in restoration [43]. Current engineering practice makes use of helically braided scaffold materials, for which dedicated finite element models are required to probe and assist their design specifications [44]. However, most numerical models consider a linear elastic material behaviour, disregarding the tendon’s viscoelastic nature. What is more, models descriptive of the tendons’ time-dependent response consider a simplified composite inner structure with non-undulated fibrillar inner components [45], while no models that inherently account for the viscoelastic, undulated helical formation of the tendon’s inner subunits have been up to now presented. Note that numerical models which account for the composite inner tendon structure come along with considerable numerical costs, necessitating discretized systems of several hundreds of thousands of degrees of freedom for a sufficiently accurate structural representation [46]. Viscoelastic, time-dependent computations can increase the numerical cost by orders of magnitude, rendering inverse engineering approaches practically infeasible [38].

In the current work, we elaborate a numerical framework for the computation of the viscoelastic response of tendon subunits. The latter is computed based on the geometric and mechanical properties of the subunit’s individual microscopic phases, namely of the fibrillar and matrix components. To that scope, in Section 2.1, we introduce the geometric, kinematic and constitutive relations characterizing the inner fibrillar-matrix structure of tendon fascicles. We subsequently provide the constitutive relations that describe the generalized viscoelastic Maxwell-Wiechert model (Section 2.2), upon which we derive the weak, finite element form of the time-dependent viscoelastic behavior of helical fascicles, explicating the model’s algorithmic implementation (Section 2.3). In Section 3.1, we make use of the model to compute the relaxation response of fascicles for the simplified case of a nearly straight inner configuration, which we compare with existing homogenization-based relaxation results. Thereafter, we analyze the role of the helical geometry and material viscosity on the fascicle’s relaxation response deriving useful conclusions on its overall viscoelastic mechanical behavior (Section 3.2). In Section 4, we analyze the model’s computational cost, elaborating its numerical complexity and computing storage requirements. In Section 5, we summarize the principal factors affecting the tendon viscoelastic behavior and conclude.

2. Methods

2.1 Fascicle geometric and constitutive relations

The tendon fascicle is described as a composite helical structure with an undulation angle $\theta$ that is composed of fibers immersed in a matrix substance [38, 46]. The ratio of the area covered by the fibrillar components $A_{f}$ to the total cross sectional $A-A^{\prime}$ fascicle area $A_{\textit{tot}}$ define the fascicle’s fibrillar content, denoted as $C_{f}=A_{f}/A_{\textit{tot}}$ . The fiber content $C_{f}$ is a non-constant model parameter, taking values within the range of 40% and 80%, in accordance with experimental observations [47, 48]. The different fascicle content structures are obtained by means of a dedicated numerical algorithm detailed in Appendix A.3, while the corresponding analysis results are provided in the form of supplementary data. Figure 2 provides a schematic representation of the fascicle geometry with a content value of 60% and a total of fourty fibers (Fig. 2b). It is to note that a minimum of some tens of fibers is required for the simulation of the experimentally observed large Poisson’s ratio behavior of tendon fascicles, a requirement that poses a lower bound on the model size [46, 38].

Figure 2.

Schematic representation of the tendon fascicle geometry (a). A fascicle section on the $\mathbf{n}-\mathbf{b}$ Serret-Frenet basis frame is composed by fibers immersed in a matrix (b). The fascicle cross section is used as the representative domain for the finite element (c) framework elaborated in Section 2.3.

Each point within the fascicle domain is parametrized by the position vector $\mathbf{X}(x_{n},x_{b},\varphi)$ , defined by means of the Serret-Frenet Curvilinear basis, as follows [5]:

$\displaystyle\mathbf{X}(x_{n},x_{b},\varphi)=\mathbf{R}(\varphi)+x_{n}\mathbf{% n}+x_{b}\mathbf{b},\quad\mathbf{R}(\varphi)=[a\cos\varphi,a\sin\varphi,b/\gamma]$ (1)

In Eq. (1), $a$ stands for the centerline of the fascicle, $\varphi$ for the rotation with respect to the loading axis of the fascicle (Fig. 2a), while $x_{n}$ and $x_{b}$ are local coordinates along the normal $\mathbf{n}$ and binormal $\mathbf{b}$ direction vectors attributed to the centerline of the fascicle section (Fig. 2b). The latter are defined as follows [5, 49]:

$\displaystyle\mathbf{n}=\begin{bmatrix}-\cos\varphi\\ -\sin\varphi\\ 0\\ \end{bmatrix}\quad\mathbf{b}=\frac{1}{\gamma}\begin{bmatrix}b\sin\varphi\\ -b\cos\varphi\\ a\\ \end{bmatrix}\quad\mathbf{t}=\frac{1}{\gamma}\begin{bmatrix}-a\sin\varphi\\ a\cos\varphi\\ b\\ \end{bmatrix}$ (2)

In Eqs (1) and (2), $b$ and $\gamma$ are characteristic geometric parameters which directly depend on the helix angle $\theta$ , defined as $b=\alpha\tan\theta$ and $\gamma=\sqrt{\alpha^{2}+b^{2}}$ accordingly [5]. The stress tensor $\sigma^{ij}$ along the fascicle cross section relates to the strain tensor $\epsilon_{kl}$ through the following constitutive expression [49, 50]:

$\displaystyle\sigma^{ij}=C^{ijkl}\epsilon_{kl},\quad\epsilon_{ij}=\frac{1}{2}(% u_{i,j}+u_{j,i})-u_{k}\Gamma_{ij}^{k},$ (3)

where in Eq. (3), superscripts represent contravariant, while subscripts covariant components. The strain tensor of Eq. (3) is defined by means of the Christoffel symbols $\Gamma_{ij}^{k}$ provided in Eq. (A.5) of the Appendix A.1. The elastic material tensor $C^{ijkl}$ is given as a function of the contravariant metric tensor provided in Eq. (A.7) of the Appendix A.1 [49]:

$\displaystyle\mathbf{C}^{el}=C^{ijkl}=\frac{\nu E}{(1+\nu)(1-2\nu)}g^{ij}g^{kl% }+\frac{E}{2(1+\nu)}(g^{ik}g^{jl}+g^{il}g^{jk})$ (4)

In Eq. (4) $\mathit{E}$ and $\nu$ stand respectively for the material’s Young modulus and Poisson’s ratio value. Depending on the position within the finite element mesh of Fig. 2, the previous parameters may refer to the fibrillar components or to the matrix substance accordingly ( $E_{f}$ and $E_{m}$ ).

2.2 Viscoelastic, Maxwell-Wiechert type mechanical response

The time-dependent response of viscoelastic materials is commonly described using the rheological analogy of the generalized Maxwell-Wiechert model [51]. The latter is composed of a linear elastic component with stiffness $E_{\infty}$ in the one dimensional case, which is time independent, positioned in parallel to a series of viscous elements with elastic moduli $E_{i}$ and viscosity $\eta_{i}$ . Each Maxwell element describes a single relaxation time $\tau$ , defined as $\tau_{i}=\frac{\eta_{i}}{E_{i}}$ . Figure 3 provides a schematic representation of the generalized Maxwell-Wiechert model.

Figure 3.

Schematic representation of the Maxwell-Wiechert model.

For an isotropic three-dimensional viscoelastic material the total inner stress $\bm{\sigma}^{n+1}$ developed at a time increment $t_{n+1}=t_{n}+\Delta{t}$ within the Maxwell-Wiechert element is given as a function of the previous timestep $t_{n}$ by the following recursive formula [52, 53, 54]:

$\displaystyle\bm{\sigma}^{n+1}=\bm{\sigma}_{\infty}^{n+1}+\sum_{j=1}^{N}\bm{h_% {j}}^{n+1},\quad\bm{\sigma}_{\infty}^{n+1}=\mathbf{C}_{el}\bm{\epsilon}_{n+1}$ (5)

where in Eq. (5), $\bm{\sigma}_{\infty}^{n+1}$ stands for the elastic stress contribution at the timestep $t_{n+1}$ , while the term $\bm{h_{j}}^{n+1}$ stands for the hereditary stress contribution of an element with a relaxation time $\tau_{j}$ . The hereditary stress contribution at the current time step depends on the one of the previous step $t_{n}$ , as well as on the moduli ratio $M_{r}=E_{i}/E_{\infty}$ , namely on the ratio of the elastic modulus of the viscous to the modulus of the time independent element (Fig. 3), as follows [52, 54]:

$\displaystyle\mathbf{h}_{j}^{n+1}=\exp\left(-\frac{\Delta{t}}{\tau_{j}}\right)% \mathbf{h}_{j}^{n}+M_{j}^{r}{A_{j}}[\bm{\sigma}_{\infty}^{n+1}-\bm{\sigma}_{% \infty}^{n}],\quad A_{j}=\left(\frac{1-\exp(-\Delta{t}/\tau_{j})}{\Delta{t}/% \tau_{j}}\right)$ (6)

In Eq. (6), $A_{j}$ is a constant term which depends on the timestep discretization $\Delta{t}$ and on the relaxation time $\tau_{j}$ of the Maxwell elements. The recursive stress formulas of Eqs (5) and (6) are valid for both for the fiber and matrix components enclosed in tendon fascicle of Section 2.1, though with different viscoelastic parameters for each phase, namely with different $E_{m}$ , $\eta_{m}$ , $E_{f}$ and $\eta_{f}$ (Fig. 2).

2.3 Finite element weak form

In the current section, we make use of the incremental stress definition of Section 2.2 and of the constitutive formulations of Section 2.1 to derive the finite element weak form describing the viscoelastic behavior of tendon fascicles. For the finite element descritization process, we express the strain tensor $\bm{\epsilon}(t)$ entering Eqs (3) and (5) as a function of the displacement vector $\mathbf{u}$ through the linear operator $\mathbf{B}$ , so that $\bm{\epsilon}(t)=\mathbf{B}\mathbf{u}$ . Thereupon, we compute the total internal work produced at increment $n+1$ , using the internal stress definitions of the Maxwell-Wiechert model provided in Eqs (5) and (6), as follows:

$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!W_{\textit{int}}^{n+1}\!=\int_{\Omega}% \!\bm{\epsilon}_{n+1}^{T}\bm{\sigma}_{\infty}^{n+1}dV\!+\!\int_{\Omega}\!\left% (\bm{\epsilon}_{n+1}^{T}\sum_{j=1}^{N}\mathbf{h}_{j}^{n+1}\right)dV\!=\int_{% \Omega}\!\bm{\epsilon}_{n+1}^{T}\mathbf{C}_{el}\mathbf{B}\!\left(\!1+\!\sum_{j% =1}^{N}M_{r,n}A_{j}\!\right)\!\mathbf{u}_{n+1}\!+\!\int_{\Omega}\!\bm{\epsilon% }_{n+1}^{T}\!\sum_{j=1}^{N}e^{-\frac{\Delta{t}}{\tau_{j}}}\mathbf{h}_{j}^{n}dV% \!-\!\!\int_{\Omega}\!\bm{\epsilon}_{n+1}^{T}\mathbf{C}_{el}\mathbf{B}\sum_{j=% 1}^{N}M_{r,n}A_{j}\mathbf{u}_{n}$ (7)

We consider equilibrium conditions by equating the internal work of Eq. (7) $W_{\textit{Int}}^{n+1}=\int_{\Omega}\mathbf{u}^{T}_{n+1}\mathbf{f}^{n+1}_{% \textit{Int}}dV$ with the external work $W_{\textit{Ext}}^{n+1}=\int_{\Omega}\mathbf{u}^{T}_{n+1}\mathbf{f}^{n+1}_{% \textit{Ext}}$ . Thereupon, after separating current ( $n+1$ ) and previous ( $n$ ) step contributions, we arrive at the following discrete form expression:

$\displaystyle\mathbf{K}_{T}\mathbf{u}_{n+1}=\mathbf{f}^{n+1}_{\textit{Ext}}-% \mathbf{f}^{n}_{\textit{Hered}}+\mathbf{f}^{n}_{\textit{Past}}$ (8)

where in Eq. (8), the term $\mathbf{K}_{T}$ stands for the tangent stiffness matrix, while the $\mathbf{f}^{n}_{\textit{Hered}}$ and $\mathbf{f}^{n}_{\textit{Past}}$ terms for the hereditary and past step force contributions are given by the following expressions:

$\displaystyle\mathbf{K}_{T}=\sum_{N_{el}=1}^{N}\mathbf{B}^{T}\mathbf{C}_{el}% \mathbf{B}\left(1+\sum_{j=1}^{N}M_{r,j}A_{j}\right)$ $\displaystyle\mathbf{f}^{n}_{\textit{Hered}}=\sum_{N_{el}=1}^{N}\mathbf{B}^{T}% \sum_{j=1}^{N}\exp\left(-\frac{\Delta{t}}{\tau_{j}}\right)\mathbf{h}_{j}^{n}$ (9) $\displaystyle\mathbf{f}^{n}_{\textit{Past}}=\sum_{N_{el}=1}^{N}\sum_{j=1}^{N}(% M_{r,j}A_{j})\mathbf{B}^{T}\bm{\sigma}_{\infty}^{n}$

In Eq. (2.3), the integral terms have been substituted by discrete summations over the corresponding element contributions. The elastic matrix $\mathbf{C}_{el}$ entering the $\mathbf{K}_{T}$ term of Eq. (2.3) is element specific, while it pertains to the elastic material modulus $E_{\infty}$ of each domain element that is time-independent, leading to a constant algorithmic stiffness matrix. Contrariwise, the hereditary and past force expressions of Eq. (2.3) are not time-independent and need to be computed at each time step. In Fig. 4, we provide a flow chart diagram of the numerical implementation of the Eqs (8) and (2.3):

Figure 4.

Flow chart diagram of the numerical implementation of the viscoelastic finite element relaxation model.

At each time increment the hereditary and past vector force expressions of Fig. 4 need to be computed anew, while the corresponding force vectors needs to be subject to the same boundary conditions employed for the linear model calculations, that is for the helix centerline to be positioned at a distance $a$ (Fig. 2), with its rigid body motions suppressed. For the computations to follow, the ratio of the fascicle’s centerline position to the fiber radius $a/r$ is set to 100, in accordance with experimental observations [55].

The finite element computations summarized in Fig. 4 are carried out on the local fascicle cross sectional plane of Fig. 2b. To that scope, a helix-specific, two-dimensional modeling form of the linear operator $\mathbf{B}$ and external, constant, axial-strain-induced ( $\epsilon_{Z}=\epsilon_{0}$ ) loading term $\mathbf{f}_{\textit{Ext}}^{n+1}$ is employed in Eqs (8) and (2.3), provided in Appendix A.1. The local element stresses are integrated over the area of each element and transformed to the global Cartesian fascicle axis using Eq. (A.2). For the computation of the global force $F_{Z}(t_{n+1})$ at each time increment (Fig. 2), the local normal domain stress are used and transformed after integration to the global tendon loading axis direction (Eq. (A.3)). Accordingly, the evolution of the global moment $M_{Z}(t_{n+1})$ with time is obtained by calculating at each time increment the addition of the local moment $\tilde{M}_{Z}(t_{n+1})$ created along the fascicle section centerline ( $x_{n},x_{b}=0$ in Eq. (1)) to the moment created by the total shear force $F_{b}(t_{n+1})$ , arising from the shear stresses along the binormal $\mathbf{b}$ local helix axis, so that $M_{Z}(t_{n+1})=\alpha F_{b}(t_{n+1})+\tilde{M}_{Z}(t_{n+1})$ .

Different fascicle compositions are simulated by creating fascicle structures of different fiber contents $C_{f}$ , upon the computational process explicated in Appendix A.3. The computational cost of the numerical models are elaborated in Section 4, while in Sections 3.1 and 3.2 the relaxation response of fascicle structures of various physiologically relevant material and geometric parameters is analyzed.

3. Results

3.1 Relaxation response of composite, non-undulated tendon fascicles

In the current Section we compare the stress relaxation behavior of axially loaded tendon fascicles (Fig. 2) obtained from the numerical model of Section 2.3 with the relaxation response obtained from homogenization analysis for the limiting case of fascicles with non-undulated inner components [45]. For the finite element computations, we use an angle of $\theta=$ 89.5 ${}^{\circ}$ . For an angle $\theta\rightarrow{90^{\circ}}$ , the Christoffel symbols of Eq. (A.5) and the curvature and tortuocity paramaters vanish with the fascicle geometry to be approximately parallel to the tendon loading axis of Fig. 2.

For the fibrillar components ( $f$ ), an elastic modulus $E_{f}=$ 200 MPa and a viscoelastic modulus $\eta_{f}=$ 5 GPa s and Poisson’s ratio $\nu_{f}=$ 0.4 is considered, along with a fibrillar content $C_{f}=$ 75% used in [45]. For the embedding matrix substance ( $m$ ) – which incorporates all non-fibrillar components, we make use of an elastic and viscoelastic material modulus that is equal to $E_{m}=$ 4 kPa with $\eta_{m}=$ 0.1 Pa s and $\nu=$ 0.45 [45]. In Fig. 5, we compare the time evolution of the axial modulus $C_{zz}(t)=F_{Z}(t)/A_{t}/\epsilon_{0}$ computed by the finite element scheme $C_{zz}^{\textit{FEM}}$ of Section 2.3 with the axial modulus retrieved by the viscoelastic modulus $C_{zz}^{\textit{HOM}}$ retrieved by homogenization analysis and presented in [45].

Figure 5.

Axial relaxation response $C_{zz}(t)$ of a composite fascicle with a 75% fibrillar content. Black colours pertain to homogenization analysis results, while red colours to the FEM results.

Figure 5 indicates a very good agreement between the relaxation response provided by the homogenization analysis method and by the finite element method. It is to note that the moduli $C_{zz}(t)$ related to the axial direction $Z$ is the primal parameter characterizing the tendon response, as the tendon is primarily loaded in tension, while the material moduli in the $X Y$ direction are affected to a considerably lower extent or remain practically unaffected during the relaxation process [45].

3.2 The role of the fiber helical geometry and material viscosity on the fascicle relaxation response

The finite element scheme of Section 2.3 does not restrain to the limiting simplified case of straight fibrillar components, but it can capture the effect of their undulated inner helical arrangement. The helical geometry of the fascicle’s inner components entails the development of a torsional moment $M_{Z}(t)$ along the tendon loading axis that evolves with time, depending on the viscoelastic paramaters of its inner material constituents. In Fig. 6, we compute the evolution of the effective axial $C_{zz}(t)$ and torsional $C_{\phi{z}}(t)=M_{Z}(t)/(A_{t}{\alpha}\epsilon_{0})$ moduli for the elastic material parameters of Section 3.1. However, we employ different viscous moduli for the fibrillar components, in particular we use viscous moduli values of $\eta=$ 10 GPa s and $\eta=$ 20 GPa s with relaxation times $\tau=$ 100 s and $\tau=$ 200 s accordingly (red and black curves in Fig. 6) [32], retaining the viscoelastic properties of the matrix unaffected and requiring that the remaining elastic stiffness is 50% of the initial elastic stiffness. The relaxation response is computed for two distinct helix angles $\theta$ equal to 70 ${}^{\circ}$ and 75 ${}^{\circ}$ [56, 57] for a fibrillar content $C_{f}$ of 50% [48].

Figure 6.

Axial (left) and torsional (right) moduli relaxation of a composite fascicle with a 50% fibrillar content. Red colors correspond to fascicles with fiber relaxation times $\tau$ of 100 s, while black colours to fascicles with fiber relaxation times of 200 s. Fibrillar angles of 75 ${}^{\circ}$ are depicted with dotted lines.

Comparing Figs 6 and 5 we note that a reduction in the fibrillar content of the fascicle results in a lower initial axial modulus $C_{zz}$ for equal elastic fiber moduli values. The difference is increasing upon decreasing helix angles (Fig. 6 left, $\theta=$ 70 ${}^{\circ}$ ), thus upon a higher fiber undulation. What is more, for a given loading time (e.x. 20 s), higher fiber viscosity and relaxation times $\tau$ relate to fascicles with a higher remaining axial $C_{zz}$ and torsional $C_{\phi{z}}$ moduli, thus to higher force $F_{z}$ and moment $M_{z}$ values. What is more, Fig. 6 suggests that for a given helical fascicle geometry, the torsional modulus $C_{\phi{z}}$ is non-vanishing, contrary to the simplifying case of a non-helical composite fibers (Fig. 5). For typical fascicle geometries, the magnitude of the torsional modulus is non-negligible, being almost one third of the corresponding axial modulus $C_{zz}$ for a fascicle content of $C_{f}=$ 50%, as a direct comparison of the initial time ( $t=$ 0) moduli values of Fig. 6 indicates. In Fig. 7, we compute the relaxation response for the case of fibrillar components with two relaxation times, namely with a low and a high relaxation time of $\tau_{1}=$ 25 s and $\tau_{2}=$ 100 s which relate to viscosity values of 1.25 GPa s and 5 GPa s [31].

Figure 7.

Axial (left) and torsional (right) moduli relaxation of a composite fascicle with a 50% fibrillar content. Fascicles with fibrillar content following two relaxation times of $\tau_{1}=$ 25 s and $\tau_{2}=$ 100 s are depicted. Fibrillar angles of 75 ${}^{\circ}$ are depicted with dotted lines.

Both Figs 6 (right) and 7 suggest that the torsional stiffness modulus value diminishes upon increasing helix angle $\theta$ , thus upon a vanishing helicacy of the tendons’ inner constituents. Furthermore, it may be observed that for a given fiber viscosity value, the time evolution of the torsional modulus $C_{\phi{z}}(t)$ is analogous to the time evolution of the axial modulus. Moreover, fascicles with fibers of high relaxation times retain a higher torsional stiffness value after a certain loading time than their corresponding low relaxation time counterparts. What is more, multiple relaxation times (Fig. 7) lead to relaxation curves with a different curvature profile with respect to the ones computed for a unique relaxation parameter (Fig. 6).

4. Discussion

The viscoelastic behavior of tendon subunits is a complex function of their inner material and geometric parameters. Upon axial loading, their response can well vary depending on a series of inner tendon parameters (Section 3.2). In particular, for a given elastic fiber moduli, the fiber content $C_{f}$ and undulation angles affect both the initial axial forces and moments created, with the difference to approach 20% (Fig. 6) for typical physiological parameter values [32, 56, 57]. What is more, the axial force and moment sustained after a certain time interval $t$ , strongly depends on the viscoelastic moduli and relaxation times of the fascicle’s inner constituents, so that losses as low as 10% (Fig. 6) or up to 50% (Fig. 5) can be obtained, depending on the fascicle’s inner material properties. As a result, it is only for certain sets of fiber-scale mechanical parameters that the upper, fascicle scale, experimentally observed relaxation response can be reproduced [29, 27]; an analysis that exceeds the scope of the current work.

It is to note, that the numerical simulation of fascicle mechanical response requires considerable computational resources, already for a linear elastic analysis [46, 20, 38]. The computational framework elaborated in the current work makes use of a fascicle section of the 3D fascicle structure for the simulation of the effective viscoelastic response (Section 2.3), requiring a minimum of approximately 6000 planar elements for the modeling results to remain invariant with respect to higher domain discretizations to a 5% accuracy. Using three translational degrees of freedom per node, a tangent stiffness matrix $\mathbf{K}_{T}$ with a dimension of approximately 18k $\times$ 18k elements arises for the computation of the iterative process summarized in Fig. 4. For such a system, a memory footprint of around 1.5 GB is required, considering double precision ( bytes per element). Figure 8 depicts the memory requirements (left) and the corresponding computing time requirements (right) for different domain discretizations. As a reference case, we consider the requirements of the simulation case with stiffness matrix of dimension 2k $\times$ 2k. We observe that the memory requirements increase quadratically with the matrix dimension $\left(1.00*\left(\frac{N}{N_{\textit{ref}}}\right)^{1.99}\right)$ . Moreover, the computing time requirements, which are determined by the solution of the linear system, increase approximately cubically $\left(1.02*\left(\frac{N}{N_{\textit{ref}}}\right)^{2.90}\right)$ .

Figure 8.

Memory (left) and run time (right) requirements of the computation model. The reference case corresponds to a linear system of 1983 $\times$ 1983 elements solved on a single-core system using the Eigen high-level C++ library.

The memory requirements presented in Fig. 8(left) are considerably higher if a full 3D domain modeling approach is used, rather than the symmetry-reduced fascicle section formulation presented in the current work, as the minimum element number requirements is substantially higher. Indicatively, full 3D fascicle computational models report the use of 144,635 nodes for a sufficient domain modeling description and a linear elastic analysis [46], while no viscoelastic models of the kind have been up to now elaborated. Considering a minimum of three degrees of freedom (dofs) per node, yields a system size of 433k $\times$ 433k elements, which corresponds to approximately 1.4 TB of memory in double precision for a linear system, thus without the necessary timestep computations for the viscoelastic response of the tendon fascicle to be obtained. Note that for an accurate viscoelastic computation, fine enough time step increments $\delta{t}$ are necessary, which are considerably lower than the relaxation time of the fascicle’s fibrillar components ( $\delta{t}\ll{\tau}$ ). However, the total run time requirements are proportional to the number of time steps. As a result, the elaborated approach manages to reduce by at least two orders of magnitude the overall computational cost, outperforming current state of the art tendon fascicle numerical modeling approaches.

5. Conclusions

In this paper, we have elaborated a modeling framework for the simulation of tendon fascicles, covering a wide range of experimentally reported tendon fascicle structural compositions, thus of fibrillar contents $C_{f}$ and undulations angle values $\theta$ [48, 56]. The framework has been validated with respect to existing viscoelastic homogenization models (Section 3.1), while its capability to compute the time-dependent evolution of the tendon fascicle’s axial and torsional modulus has been showcased for physiologically relevant material and geometric parameters in Section 3.2. It has been shown that the torsional modulus is significant for typical fascicle configurations, while its time-evolution well compares to the one observed for the axial stiffness. What is more, its computational cost has been analyzed and compared with existing modeling approaches, highlighting its merits in the simulation of the viscoelastic, time-dependent relaxation response of tendon fascicles (Section 4). Thereupon, it has been shown that the elaborated numerical scheme outperforms existing modeling approaches, showing promising prospects in the simulation of the complex tendon viscous mechanics (Section 4). As such, we aspire that the proposed framework will provide a research basis in the tendon treatment and restoration praxis, providing a link between the tendon’s inner properties and its macroscopic, experimentally observed mechanical behavior [38, 26].

Footnotes

Acknowledgments

N.K. would like to gratefully acknowledge the support of the Freenovation Grant 2017 and the support of the ETH CSE lab.

Appendix A.1. Helix basis and finite element operators

The helical domain defined by Eqs (1) and (2) pertains to the following contravariant metric tensor $[g^{ij}]$ [22]:

(A.1) $\displaystyle[g^{ij}]=\frac{1}{g}\begin{bmatrix}g+(\tau x_{b})^{2}&-\tau^{2}x_% {n}x_{b}&\tau x_{b}\\ -\tau^{2}x_{n}x_{b}&g+(\tau x_{n})^{2}&-\tau x_{n}\\ \tau x_{b}&-\tau x_{n}&1\\ \end{bmatrix}$

In Eq. (A.1), $g=(1-\kappa x_{n})^{2}$ stands for the determinant of the metric tensor, while $\kappa=\frac{a}{\gamma^{2}}$ and $\tau=\frac{b}{\gamma^{2}}$ denote accordingly the curvature and tortuocity of the helical domain. The parameters $a,b,\gamma$ entering the curvature and tortuocity definitions are characteristic geometric features of the helix explicated in Section 2.1. The Curvilinear basis upon which the strain and stress components of Eqs (3 and (4) are defined, relates to the Cartesian one $\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z}$ through the transformation tensor $\mathbf{F}$ . The latter is given as a function of the previously introduced geometric parameters as follows [22]:

(A.2) $\displaystyle\mathbf{F}=\frac{1}{\gamma}\begin{bmatrix}-\gamma C&bS&x_{n}S+% \gamma\tau x_{b}C-aS\\ -\gamma S&-bC&-x_{n}C+\gamma\tau x_{b}S+aC\\ 0&a&b\\ \end{bmatrix}$

where $C$ and $S$ in Eq. (A.2) stand for the $\cos\varphi$ and $\sin\varphi$ respectively. The centerline position can be taken without loss of generality to be equal to zero ( $\varphi=$ 0). A global (G) Cartesian force or moment vector $\mathbf{f}_{G}$ can be obtained using the corresponding local Curvilinear vector $\mathbf{f}_{C}$ as follows [58]:

(A.3) $\displaystyle\mathbf{f}_{G}=\mathbf{F}\mathbf{f}_{C}$

The linear operator $\mathbf{B}$ entering the weak form of the finite-element model description of Eqs (8) and (2.3) is described as a function of the domain shape functions $\phi$ and incorporates the consideration of a symmetric behaviour with respect to the helix wave guide evolution [21]. The consideration of symmetry leads to a linear operator $\mathbf{B}$ and local Curvilinear strain vector $\bm{\epsilon}$ expression that are a function of the helix geometric attributes as follows [22]:

(A.4) $\displaystyle\mathbf{B}=\begin{bmatrix}\phi_{1,1}&0&0&\ldots\\ 0&\phi_{1,2}&0&\ldots\\ -\Gamma_{33}^{1}\phi_{1}&-\Gamma_{33}^{2}\phi_{1}&-\Gamma_{33}^{3}\phi_{1}&% \ldots\\ -\Gamma_{23}^{1}\phi_{1}&0&\frac{1}{2}\phi_{1,2}&\ldots\\ -\Gamma_{13}^{1}\phi_{1}&-\Gamma_{13}^{2}\phi_{1}&\frac{1}{2}\phi_{1,1}-\Gamma% _{13}^{3}\phi_{1}&\ldots\\ \frac{1}{2}\phi_{1,2}&\frac{1}{2}\phi_{1,1}&0&\ldots\\ \end{bmatrix},\quad\bm{\epsilon}=\begin{bmatrix}0\\ 0\\ b\tau\\ \alpha\tau\\ 0\\ 0\\ \end{bmatrix}\epsilon_{z}=-\begin{bmatrix}0\\ 0\\ b\tau\\ \alpha\tau\\ 0\\ 0\\ \end{bmatrix}\epsilon_{0}$

where in Eq. (A.4) the $\bm{\epsilon}$ stands for the constant strain tensor applied to the local Curvilinear helix domain. The reader is referred to [21] for the derivation details. The Christoffel symbols entering the linear operator $\mathbf{B}$ of Eq. (A.4) are given for the Serret-Frenet basis of Eq. (1) as follows:

(A.5) $\displaystyle\Gamma^{1}_{ij}=\begin{bmatrix}0&0&-\frac{\kappa\tau{x_{b}}}{1-% \kappa{x_{n}}}\\ 0&0&-\tau\\ -\frac{\kappa\tau x_{b}}{1-\kappa x_{n}}&-\tau&\frac{\kappa(\tau x_{b})^{2}}{1% -\kappa x_{n}}+\kappa(1-\kappa x_{n})-\tau^{2}x_{n}\\ \end{bmatrix}$ (A.6) $\displaystyle\Gamma^{2}_{ij}=\begin{bmatrix}0&0&\frac{\kappa\tau{x_{n}}}{1-% \kappa{x_{n}}}+\tau\\ 0&0&0\\ \frac{\kappa\tau x_{n}}{1-\kappa x_{n}}+\tau&0&-\frac{\kappa\tau^{2}x_{n}}{x_{% b}}{1-\kappa x_{n}}-\tau^{2}x_{b}\\ \end{bmatrix},\quad\Gamma^{3}_{ij}=\begin{bmatrix}0&0&-\frac{\kappa}{1-\kappa{% x_{n}}}\\ 0&0&0\\ -\frac{\kappa}{1-\kappa x_{n}}&0&\frac{\kappa\tau x_{b}}{1-\kappa x_{n}}\\ \end{bmatrix}$

For a relaxation experiment, the externally applied force at each element (elem) within the discretized domain $\mathbf{f}^{\textit{Ext}}_{\textit{elem}}=\int_{A}\tilde{\bm{\sigma}}^{\textit% {Ext}}dA$ is computed based on the constant applied strain vector expression of Eq. (A.4). The externally applied stress at each domain element is computed using the constitutive expression of Eq. (3) $\tilde{\bm{\sigma}}^{\textit{Ext}}=\tilde{\mathbf{C}}(\tilde{E})\bm{\epsilon}$ , using the element-specific modulus $\tilde{E}=\Upsilon(t)$ at each timestep. The latter is computed based on the relaxation characteristics ( $\tau_{i}$ ) of the specific element, as elaborated in Section Appendix A.2.

Appendix A.2. Material modulus relaxation kernels

In the current section, we compute the relaxation kernel $\tilde{E}=\Upsilon(t)$ for an element of the finite element domain with area $A$ and relaxations times $\tau_{i},i\in{1\dots{N}}$ , $N$ being the number of its relaxation times that is subject to a constant external strain $\epsilon_{\textit{const}}$ . The externally applied strain is described by means of the Heaviside function $H$ , as follows $\epsilon(t)=\epsilon_{\textit{const}}{H}$ . The total force $f(t)$ developed is computed as the addition of the force contributions of its elastic $f_{E_{\infty}}$ and Maxwell components (elements in parallel), as follows $f(t)=f_{E_{\infty}}+\sum_{j=1}^{N}f_{R_{j}}(t)$ , $R_{j}$ denoting the relaxation components with relaxation times $\tau_{j}=\eta_{j}/E_{j}$ [51]. The derivative of the strain $\epsilon(t)$ within a relaxing element $R_{j}$ is equal to the sum (elements in series) of the strain derivatives of its elastic and viscous contributions, obeying to the following differential expression [51]:

(A.7) $\displaystyle\dot{\epsilon}(t)=\dot{\epsilon}(t)_{E_{j}}+\dot{\epsilon}(t)_{% \eta_{j}}=\frac{\dot{f_{R_{j}}}(t)}{AE_{j}}+\frac{f_{R_{j}}(t)}{A\eta_{j}}% \Rightarrow\bar{f}_{R_{j}}=E_{j}A\frac{{s}}{s+1/\tau_{j}}\bar{\epsilon}$

where at the right hand side of Eq. (A.7), the Laplace transform $L(s)$ of the time dependent strain function has been used ( $L(\dot{\epsilon})=s\bar{\epsilon},L(\dot{f})=s\bar{f}$ ), with the corresponding Laplace transformed force and strain components to be accordingly denoted as $\bar{f}$ and $\bar{\epsilon}$ . In Eq. (A.7), $A$ corresponds to the element area, while $\tau_{j}$ to the element’s relaxation time. For a relaxation test, the Laplace transform of an applied Heaviside type strain is given as $\bar{\epsilon}=\epsilon_{\textit{const}}/s$ . Thus, the Laplace transformed expression of the total force of Eq. (A.7) is given as follows [51]:

(A.8) $\displaystyle\bar{f}=sE_{\infty}A\bar{\epsilon}+\sum_{i=1}^{i=N}\bar{f}_{R_{i}% },\quad i=1,2,\ldots{N}$

The time domain element applied force is given by the inverse Laplace transform $L^{-1}(s)$ of Eq. (A.8), which provides the element moduli relaxation kernel $\tilde{E}=\Upsilon(t)$ , as follows [51].

(A.9) $\displaystyle f(t)=\Upsilon(t)A\epsilon_{\textit{const}},\quad\Upsilon(t)=% \left(E_{\infty}+\sum_{j=1}^{N}E_{j}\exp(-t/\tau_{j})\right)$

Appendix A.3. Fascicle fiber spatial positioning

The composite fiber-matrix fascicle structures (Fig. 2) introduced in Section 2.1 are constructed considering a homogeneous distribution of fibers within a tendon fascicle for different values of contents $C_{f}$ . The positioning of $N$ fibers with radius $r$ within a fascicle cross section $A-A^{\prime}$ is recast to a classical Thompson optimization problem with 2 $N$ unknowns [59], the unknowns representing the $x_{n}$ , $x_{b}$ positions of the fibers, as explicated in [38]. The optimization problem is solved making use of the Covariance Matrix Adaptation Evolutionary Strategy (CMA-ES) method [60]. For completeness, we provide attached the normalized $r/R$ coordinates $x_{n},x_{b}$ of the fiber for contents $C_{f}$ as low as 40% and up to 75%, upon an increment of 5%.

References

Yahia

L.-H.

and Drouin

, Microscopical investigation of canine anterior cruciate ligament and patellar tendon: Collagen fascicle morphology and architecture, Journal of Orthopaedic Research 7(2) (1989), 243–251.

de Campos Vidal

, Image analysis of tendon helical superstructure using interference and polarized light microscopy, Micron (Oxford, England: 1993) 34(8) (January 2003), 423–32.

Maceri

Marino

and Vairo

, An insight on multiscale tendon modeling in muscle-tendon integrated behavior, Biomechanics and Modeling in Mechanobiology 11(3) (Mar 2012), 505–517.

Thorpe

C.T.

Klemt

Riley

G.P.

Birch

H.L.

Clegg

P.D.

and Screen

H.R.C.

, Helical sub-structures in energy-storing tendons provide a possible mechanism for efficient energy storage and return, Acta Biomaterialia 9(8) (August 2013), 7948–56.

Lee

W.K.

, An insight into wire rope geometry, International Journal of Solids and Structures 28(4) (1991), 471–490.

Papailiou

K.O.

, On the bending stiffness of transmission line conductors, Power Delivery, IEEE Transactions on 12(4) (Oct 1997), 1576–1588.

Costello

G.A.

, Theory of Wire Rope, Springer, New York, 1997.

Karathanasopoulos

, Numerical characterization of the structural response of helical constructions to radial and thermal loads, Journal of Computational Methods in Sciences and Engineering 16(4) (2016), 787–800.

Lanteigne

, Theoretical estimation of the response of helically armored cables to tension, torsion, and bending, Journal of Applied Mechanics 52 (June 1985), 423.

10.

Karathanasopoulos

Ganghoffer

J.F.

and Papailiou

K.O.

, Analytical closed-form expressions for the structural response of helical constructions to thermal loads, International Journal of Mechanical Sciences 117 (2016), 258–264.

11.

Sathikh

Moorthy

M.B.K.

and Krishnan

, A symmetric linear elastic model for helical wire strands under axisymmetric loads, The Journal of Strain Analysis for Engineering Design 31(5) (September 1996), 389–399.

12.

Raoof

and Kraincanic

, Prediction of coupled axial/torsional stiffness coefficients of locked-coil ropes, Computers and structures 69(1) (April 1998), 305–319.

13.

Yildirim

, Exact determination of the global tip deflection of both close-coiled and open-coiled cylindrical helical compression springs having arbitrary doubly-symmetric cross-sections, International Journal of Mechanical Sciences, 115 (2016), 280–298.

14.

Jiang

W.G.

Yao

M.S.

and Walton

J.M.

, A concise finite element model for simple straight wire rope strand, International Journal of Mechanical Sciences 41(2) (February 1999), 143–161.

15.

Messager

and Cartraud

, Homogenization of helical beam-like structures: Application to single-walled carbon nanotubes, Computational Mechanics 41(2) (June 2007), 335–346.

16.

Fedorko

Stanova

Molnar

Husakova

and Kmet

, Computer modelling and finite element analysis of spiral triangular strands, Advances in Engineering Software 73 (2014), 11–21.

17.

Nawrocki

and Labrosse

, A finite element model for simple straight wire rope strands, Computers and Structures 77 (2000), 345–359.

18.

Stanova

Fedorko

Fabian

and Kmet

, Computer modelling of wire strands and ropes part II: Finite element-based applications, Advances in Engineering Software 42(6) (June 2011), 322–331.

19.

Stanova

Fedorko

Kmet

Molnar

and Fabian

, Finite element analysis of spiral strands with different shapes subjected to axial loads, Advances in Engineering Software 83 (2015), 45–58.

20.

Reese

S.P.

Ellis

B.J.

and Weiss

J.A.

, Micromechanical model of a surrogate for collagenous soft tissues: Development, validation and analysis of mesoscale size effects, Biomechanics and Modeling in Mechanobiology 12(6) (Nov 2013), 1195–1204.

21.

Karathanasopoulos

and Kress

, Two dimensional modeling of helical structures, an application to simple strands, Computers and Structures 174 (2016), 79–84.

22.

Karathanasopoulos

Reda

and Ganghoffer

J.F.

, Finite element modeling of the elastoplastic axial-torsional response of helical constructions to traction loads, International Journal of Mechanical Sciences 133 (2017), 368–375.

23.

Kmet

and Mojdis

, Postelastic analysis of cable trusses, Journal of Structural Engineering 141(10) (2015).

24.

Karathanasopoulos

and Angelikopoulos

, Optimal structural arrangements of multilayer helical assemblies, International Journal of Solids and Structures (2016), 1–8.

25.

Machiraju

Phan

A.-V.

Pearsall

A.W.

and Madanagopal

, Viscoelastic studies of human subscapularis tendon: Relaxation test and a wiechert model, Computer Methods and Programs in Biomedicine 83(1) (2006), 29–33.

26.

Karathanasopoulos

Arampatzis

and Ganghoffer

J.-F.

, Unravelling the viscoelastic, buffer-like mechanical behavior of tendons: A numerical quantitative study at the fibril-fiber scale, Journal of the Mechanical Behavior of Biomedical Materials 90 (2019), 256–263.

27.

Davis

F.M.

and De Vita

, A nonlinear constitutive model for stress relaxation in ligaments and tendons, Annals of Biomedical Engineering 40(12) (2012), 1–10.

28.

Elliott

D.M.

Robinson

P.S.

Gimbel

J.A.

Sarver

J.J.

Abboud

J.A.

Iozzo

R.V.

and Soslowsky

L.J.

, Effect of altered matrix proteins on quasilinear viscoelastic properties in transgenic mouse tail tendons, Annals of Biomedical Engineering 31(5) (2003), 599–605.

29.

Screen

H.R.C.

, Investigating load relaxation mechanics in tendon, Journal of the Mechanical Behavior of Biomedical Materials 1(1) (2008), 51–58.

30.

Shen

, Tensile Mechanical Properties of Isolated Collagen Fibrils Obtained by Microelectromechanical Systems Technology, PhD thesis, 2010.

31.

Shen

Z.L.

Kahn

Ballarini

and Eppell

S.J.

, Viscoelastic properties of isolated collagen fibrils, Biophys J 100(12) (June 2011), 3008–15.

32.

Sopakayang

De Vita

Kwansa

and Freeman

J.W.

, Elastic and viscoelastic properties of a type I collagen fiber, Journal of Theoretical Biology 293 (2012), 197–205.

33.

Kato

Y.P.

Christiansen

D.L.

Hahn

R.A.

Shieh

S.J.

Goldstein

J.D.

and Silver

F.H.

, Mechanical properties of collagen fibres: A comparison of reconstituted and rat tail tendon fibres, Biomaterials 10(1) (January 1989), 38–42.

34.

Karathanasopoulos

and Ganghoffer

J.F.

, Investigating the effect of aging on the viscosity of tendon fascicles and fibers, Frontiers in Bioengineering and Biotechnology 7 (2019), 107.

35.

Franciosi

Spagnuolo

and Salman

O.U.

, Mean green operators of deformable fiber networks embedded in a compliant matrix and property estimates, Continuum Mechanics and Thermodynamics 31(1) (Jan 2019), 101–132.

36.

Misra

Parthasarathy

Singh

and Spencer

, Micro-poromechanics model of fluid-saturated chemically active fibrous media, ZAMM – Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 95(2) (2013), 215–234.

37.

Giorgio

Andreaus

dell’Isola

and Lekszycki

, Viscous second gradient porous materials for bones reconstructed with bio-resorbable grafts, Extreme Mechanics Letters 13 (2017), 141–147.

38.

Karathanasopoulos

Angelikopoulos

Papadimitriou

, and Koumoutsakos

, Bayesian identification of the tendon fascicle’s structural composition using finite element models for helical geometries, Computer Methods in Applied Mechanics and Engineering 313 (2017), 744–758.

39.

Handsfield

G.G.

Inouye

J.M.

Slane

L.C.

Thelen

D.G.

Wilson Miller

and Blemker

S.S.

, A 3d model of the achilles tendon to determine the mechanisms underlying nonuniform tendon displacements, Journal of Biomechanics 51 (2017), 17–25.

40.

Starly

Yildirim

and Sun

, A tracer metric numerical model for predicting tortuosity factors in three-dimensional porous tissue scaffolds, Computer Methods and Programs in Biomedicine 87(1) (2007), 21–27.

41.

Abdullah

, Usage of synthetic tendons in tendon reconstruction, BMC Proceedings 9(3) (May 2015), A68.

42.

Rawson

S.D.

Margetts

Wong

J.K.F.

and Cartmell

S.H.

, Sutured tendon repair; a multi-scale finite element model, Biomechanics and Modeling in Mechanobiology 14(1) (Jan 2015), 123–133.

43.

Sun

Wang

and Phee

S.J.

, Modeling and motion compensation of a bidirectional tendon-sheath actuated system for robotic endoscopic surgery, Computer Methods and Programs in Biomedicine 119(2) (2015), 77–87.

44.

Laurent

C.P.

Latil

Durville

Rahouadj

Geindreau

Orgéas

and Ganghoffer

J.F.

, Mechanical behaviour of a fibrous scaffold for ligament tissue engineering: Finite elements analysis X-ray tomography imaging, Journal of the Mechanical Behavior of Biomedical Materials 40 (2014), 222–233.

45.

Ganghoffer

J.F.

Laurent

Maurice

Rahouadj

and Wang

, Nonlinear viscous behavior of the tendon’s fascicles from the homogenization of viscoelastic collagen fibers, European Journal of Mechanics – A/Solids 59 (2016), 265–279.

46.

Reese

S.P.

Maas

and Weiss

, Micromechanical models of helical superstructures in ligament and tendon fibers predict large Poisson’s ratios, Journal of biomechanics 43(7) (May 2010), 1394–400.

47.

Svensson

R.B.

Hassenkam

Grant

and Magnusson

S.P.

, Tensile properties of human collagen fibrils and fascicles are insensitive to environmental salts, Biophys J 99(12) (December 2010), 4020–7.

48.

Lavagnino

Arnoczky

S.P.

Frank

and Tian

, Collagen fibril diameter distribution does not reflect changes in the mechanical properties of in vitro stress-deprived tendons, J Biomech 38(1) (January 2005), 69–75.

49.

Chapelle

and Bathe

K.J.

, The Finite Element Analysis of Shells: Fundamentals, Springer, 2011.

50.

Wempner

, Mechanincs of Solids with Applications to Thin Bodies, McGraw-Hill, New York, [1973], 1981.

51.

Christensen

R.M.

, Theory of Viscoelasticity, Academic Press, London, 1982.

52.

Taylor

R.L.

Pister

K.S.

and Goudreau

G.L.

, Thermomechanical analysis of viscoelastic solids, International Journal for Numerical Methods in Engineering 2(1) (1970), 45–59.

53.

Simo

J.C.

, On a fully three-dimensional finite-strain viscoelastic damage model: Formulation and computational aspects, Computer Methods in Applied Mechanics and Engineering 60(2) (1987), 153–173.

54.

Kaliske

and Rothert

, Formulation and implementation of three-dimensional viscoelasticity at small and finite strains, Computational Mechanics 19(3) (Feb 1997), 228–239.

55.

K.-N.

Sun

Y.-L.

and Luo

Z.-P.

, Flexibility of type i collagen and mechanical property of connective tissue, Biorheology 41(3–4) (2004), 239–246.

56.

Järvinen

T.H.

Järvinen

T.L.N.

Kannus

Józsa

and Järvinen

, Collagen fibres of the spontaneously ruptured human tendons display decreased thickness and crimp angle, Journal of Orthopaedic Research: Official Publication of the Orthopaedic Research Society 22(6) (November 2004), 1303–9.

57.

Starborg

Kalson

N.S.

Mironov

Cootes

T.F.

Holmes

D.F.

and Kadler

K.E.

, Using transmission electron microscopy and 3View to determine collagen fibril size and three-dimensional organization, Nature Protocols 8 (Jun 2013), 1433.

58.

Itskov

, Tensor Algebra and Tensor Analysis for Engineers, Springer, 2007.

59.

Thomson

J.J.

, On the structure of the atom: An investigation of the stability and periods of oscillation of a number of corpuscles arranged at equal intervals around the circumference of a circle; with application of the results to the theory of atomic structure, Philosoph Mag Ser 7(39) (1904), 237–265.

60.

Hansen

and Ostermeier

, Completely derandomized self-adaptation in evolution strategies, Evolutionary Computation 9(2) (2001), 159–195.