Optimization of the design of a barbed suture for flexor tendon repair using extended finite element analysis

Cross-sectional area of the flexor tendons

The main motivation behind varying the aspect ratio of the suture arises from the elliptical nature of the flexor tendons cross-section. It is anticipated that better conformity between the suture and the tendon is achieved when both possess similar aspect ratios. Widths and heights of Flexor Digitorum Superficialis and Profundus tendons (FDS and FDP respectively) for the different fingers were obtained from the literature[23,24], in order to evaluate their aspect ratios and gain insight on possible design values for “ρ” [Table 1].

The dimensions of the smallest tendon should be taken as the reference design such that the developed suture is suitable for use on all flexor tendons. As Table 1 shows, the FDS of the little finger has the smallest major and minor axes; a = 1.5 mm and b = 0.858 mm. On average, the aspect ratio “ρ” is 2 and varies between a maximum of 2.88 for the FDS of the middle finger and a minimum of 1.46 for the FDP of the ring finger.

Based on these findings, the values of “ρ” were varied between 1/3 and 3. The effect of the aspect ratio on the strength of the suture was investigated for constant cut angle and cut depth. Also, three different cut depths and cut angles: 0.07-mm, 0.12-mm, and 0.19-mm as depths and 150 154 and 160 as angles were studied for a constant “ρ”.

Furthermore, the radius of the suture (or the equivalent ellipse axes) was increased from 0.3 mm to a larger diameter that is suitable with both the dimensions of the smallest tendon (FDS of little finger) and the optimum aspect ratio. It is desired that the optimized suture accommodates about 30% of the tendon’s cross-sectional area. The diagram in Figure 4 clarifies the scheme that was followed.

Figure 4. Flowchart clarifying the plan followed throughout the analysis

The values for diameter, cut angle, and cut depth were set for a circular cross-sectional area (CSA). When considering elliptical profiles, the major and minor axes were determined such that the area of the ellipse is equal to that of the circle. Furthermore, while the cut angles remain the same, the cut depths should be adapted for the ellipses in order to retain the same amount of intact area after inducing a cut in the sutures. The performed calculations are detailed below.

Finding major and minor axes:

(CSA)_ellipse = (CSA)_{circle and πab = πr}2

and for an ellipse,

ρ = a/b

therefore,

a = r × ρ^1/2 and b = r/(ρ^1/2)

To find equivalent cut depth, the remaining area of the ellipse [Figure 5] can be found using the following relation:

Figure 5. Diagram of an ellipse with a cut on its major axis

(RA)ellipse = a × b × [π/2 + sin^-1 (d/a)] + [b × d × (a² - d²)^1/2]/a (Equation 1)

Where “a” and “b” are the major and minor axes respectively, and “d” is the distance from the centre of the cut (d = a - cut depth).

Since the remaining area of the ellipse should be equal to that of the circle, Equation (1) was used with a = b = r (radius of the circle) along with the desired cut distance “d” in order to determine the remaining area of the circle. Once the later was known, it was set as the remaining area of the ellipse and the equation was solved for the equivalent “d”. The corresponding cut depth is simply found by subtracting the major axis from the cut distance; that is, cut depth = a - d [Table 2].

For simplification purposes, a single barb was included in the model with five different geometries used for this experimental analysis [Figure 6].

Figure 6. Different geometries used in modeling

Material modeling

The material properties of polypropylene which correspond to the material usually used for flexor tendon repair were obtained after performing tensile tests to failure on 12 different suture specimens. This allowed us to determine polypropylene’s short-term linear elastic behavior, which requires the definition of an instantaneous Young’s Modulus and a Poisson’s ratio [Table 3].

The time-dependent behavior, on the other hand, was specified as previously described, via inputting into ABAQUS normalized shear relaxation data with respect to the initial stress along with the respective time instants as obtained from the relaxation tests. ABAQUS then performs a curve fitting analysis to determine the prony series coefficients. Computed values for the normalized shear g_i and bulk k_i moduli, in addition to the corresponding relaxation times t_i [Table 4] and a plot of the normalized stress relaxation behavior and its curve were calculated [Figure 7].

Figure 7. Shear relaxation curve from experimental data (red) and ABAQUS curve fitting (blue)

It is noteworthy to point out that the relaxation times gotten from experimentation were scaled by 100 when inputted into ABAQUS such that the total time becomes 36 s instead of 3600 s (60 min). This was done in order to avoid lengthy simulations and to be able to capture polypropylene’s relaxation behavior within the model.

As for damage initiation, the maximum principal stress was chosen[24-28]. Since this criterion defines the nucleation of a defect within the body, and since the geometric model at hand has a crack already defined, the value for Maxps is not of high significance for this crack propagation analysis. The calculated true ultimate strength was used as the limiting maximum principal stress. The conversion of the engineering tensile strength obtained from tensile experiments to its true equivalent was done using the following relation: σ_true = σ_eng (1 + σ_eng).

Applied to polypropylene [Table 3]: σ_true = 353.147 (1 + 0.39145) = 491.39 MPa, the value of Maxps was set at 491.39 MPa specified at the crack tip. For a mesh that is not sufficiently refined in the vicinity of the crack tip, the default centroidal approximation of the Maxps may not be accurate. It is recommended that the stresses be extrapolated to the crack tip[29-36]. As for tolerance, the default value of 0.05 was kept.

For a pre-existing crack, damage evolution is much more important for the analysis. Based on the principles of LEFM, a power law mixed mode energy-based criterion was selected with the power set as 1. The critical Mode I, Mode II, and Mode III energy release rates (G_Ic, G_IIc, and G_IIIc respectively) should be specified. These values were obtained as critical stress intensity factors K_Ic and K_IIc form the literature[37] and converted to values for GIc and GIIc assuming plane strain conditions using the following equations:

G_Ic = [K_IC² (1 - v²)]/E and G_IIc = [K_IIC² (1 - v²)]/E

Where E is the Young’s Modulus and v is the Poisson’s ratio.

For K_Ic = 5.5 MPa (m)^1/2 , and K_IIc = o.44 K_Ic = 2.42 MPa (m)^1/2:

G_Ic = [K_IC² (1- v²)]/E = [(5.5)² (1 - 0.36878²)]/1668.67 = 15.663 MPa.mm

and

G_IIc = [K_IIC² (1- v²)]/E= [(2.42)² (1 - 0.36878²)]/1668.67 = 3.0323 MPa.mm

Or alternatively,

G_IIC = (0.44)² G_IC = 3.0323 MPa.mm

For Mode III, it was assumed that the critical energy release rate G_IIIc was equal to that of mode II due to the lack of fracture toughness data for this out-of-plane shear. As such, G_Ic, G_IIc, and G_IIIc were set to be equal to 15.663 MPa.mm, 3.0323 MPa.mm, and 3.0323 MPa.mm respectively.

For damage stabilization, viscous regularization was used to specify viscosity coefficients[25-30]. When the model exhibits material softening and stiffness degradation, severe convergence problems arise. A viscous regularization scheme helps improve convergence by stabilizing the model during damage. The value for the viscosity coefficient should be chosen such that stabilization does not significantly influence the final results. To make sure that this is the case, the viscous dissipation output (ALLVD) should be small compared to the strain energy (ALLSE). Throughout this analysis, viscous regularization was chosen to be 1E-005. ALLVD and ALLSE for a circular suture having a cut angle 154 and cut depth 0.19-mm are shown in Figure 8. ALLVD is nearly zero and is very small compared to ALLSE. Similarly, the selected value was checked for all other models to ensure that it did not affect the solution.

Figure 8. ALLSE and ALLVD for the whole model with respect to time

Crack modeling

For the definition of the crack in ABAQUS[27-33], the suture geometry was selected as the domain while the barb created as a cell partition was specified as the crack [Figure 9]. An XFEM crack growth interaction was created and crack propagation was allowed in the initial step.

Figure 9. Defining the crack using XFEM

As aforementioned, XFEM based on VCCT was implemented as the fracture criterion to model crack propagation. A contact interaction property was created with energy-based mixed mode power law. Maximum tangential shear stress was chosen to be the direction of crack growth. A tolerance of 0.05 and viscosity of 1E-005 were specified. Tolerance for unstable crack propagation was toggled on and kept at its default value of 0.2.

Mode I, Mode II, and Mode III critical energy release rates were given the same values as those previously specified for the damage evolution criterion. As for the powers a_m, a_n and a_o [Table 5], they were found in the literature to be 1, 2, and 1 respectively[27].

In order to visualize a crack when using XFEM analysis, PHILSM and PSILSM should be requested as output. PSILSM is required to view the initial crack front, while PHILSM is necessary to observe the location of the crack. Also, the output STATUSXFEM should be chosen in order to perceive the status of the enriched elements. A value of 0 for this output means that the element is undamaged, a value of 1 means that the element has been completely cut by the crack, and a value between 1 and 0 signifies that the element is damaged but some traction forces remain.

Boundary conditions and loads

The performance of the barbed suture was investigated under two loading conditions. The first case studies the strength of the suture and simulates a tensile test whereby one end is held fixed while the other is subjected to a finite displacement. The second case investigates the strength of a single barb. The boundary conditions are similar to those of the first load case, but in this analysis, the surface of the barb is constrained in the direction of application of the displacement[34-37].

Three loading situations were performed. In loading case 1, the right extremity was fixed by an ENCASTRE boundary condition applied in the initial step and propagated to the following Static and Visco steps. As for the left edge, a displacement/rotation boundary condition was specified at the surface. The displacement in z-direction, U3 was set to 0 initially. However, it was modified in the first Static step to a ramped displacement of magnitude 0.5-mm and was propagated to the subsequent Visco step [Figure 10].

Figure 10. Boundary conditions for load case 1

Average stresses and strains in the body of the suture below the crack tip were of interest at the onset of failure. For this reason, an element set, Elem1, was created and associated with a Field Output Request to report strains and stresses in the z-direction, NE33 and S33 respectively [Figure 11].

Figure 11. The position of Elem1 set highlighted in red

In loading case 2, the left extremity was fixed by an ENCASTRE boundary condition applied in the initial step and propagated to the following Static and Visco steps. As for the right edge, a displacement/rotation boundary condition was specified at the surface. The displacements in x, y, and z-directions, U1, U2, and U3 respectively as well as the rotations in these directions, UR1, UR2, and UR3 respectively were set to 0 initially. However, U3 was modified in the first Static step to a ramped displacement of magnitude 0.3-mm and was propagated to the subsequent Visco step. Furthermore, a third boundary condition was applied at the crack which was constrained in the direction of application of the displacement; roller in z-direction. Since the created partition to represent the barb was defined as a crack, loads and boundary conditions could not be defined exactly at the surface. For this reason, a new partition was created by offsetting the plane of the crack by 0.005-mm. The roller boundary condition was applied at that surface [Figure 12].

Figure 12. Boundary conditions for load case 2

Since the main concern from applying this test case is to investigate the strength of a single barb, average stresses at the surface of the crack were monitored at the onset of failure. For this reason, an element set, BabrElem, was created and associated with a Field Output Request to report Von-mises stresses [Figure 13].

Figure 13. The position of BarbElem set highlighted in red

Finally, loading case 3 was set in order to investigate the behavior of the suture when loads are applied to both, the body of the suture and the surface of the barb, so a third test case which combines the former test conditions was applied. The significance of this study lies in its resemblance to the real-life suture-tendon interaction. These loading conditions were considered for the final optimized geometry only.

After assessing the forces exerted in the different repaired FDS & FDP tendons during the five different protocols of rehabilitation [Table 6], the highest median forces were noted during FDP maneuver and can reach a mean peak of 25 N[38,39].

A first factor of safety of 18% was added in order to withstand the tendon softening and the decrease in the repair strength during the early days after tendon repair and a second factor of safety of 30% was considered to compensate the risk of gap formation in the zone of repair. So the force of 25 N becomes 37 N and this was comforted by the literature which advices an ultimate repair strength of at least 30 to 55 N with a mean force of 40 N when an active rehabilitation maneuvers are involved[40].

So in loading case 3, the left extremity was subjected to a displacement/rotation boundary condition. The displacements in x, y, and z-directions, U1, U2, and U3 respectively as well as the rotations in these directions, UR1, UR2, and UR3 respectively were set to 0 initially. However, U3 was modified in the first Static step to a ramped displacement of magnitude 0.096-mm equivalent to 40N, as obtained from the load-displacement curve generated for the model on ABAQUS [Figure 14]. The right extremity was fixed by an ENCASTRE boundary condition applied in the initial step and propagated to the following Static and Visco steps. Furthermore, a third boundary condition was applied at the crack. A ramped displacement of magnitude 0.05-mm was defined perpendicular to the offset surface that was partitioned in the second test case [Figure 15]. This boundary condition was created in the Visco step such that the body of the suture has been preloaded with a value of 40 N.

Figure 14. Force-Displacement curve as obtained from ABAQUS

Figure 15. Boundary conditions for load case 3

Meshing

The impact of the quality of the mesh on the correctness of the generated solution can never be overemphasized. A high-quality mesh is one that serves the objective of the model; obtain a solution with the needed accuracy by using only as many degrees of freedom as necessary.

The mesh is refined until a critical result, such as the maximum stress in a specific location, converges; that is, it does not change significantly with further refinement (h-method).

A mesh convergence analysis was performed for this study for the two loading cases. The investigation was carried out for the circular suture (ρ = 1) having a cut angle of 154, and a cut depth 0.19-mm.

In loading case 1, meshing was done using linear tetrahedral elements (C3D4). A mesh sensitivity analysis was carried to determine the optimum mesh size among the tested values: 0.015, 0.02, 0.025, 0.027, 0.03, 0.031, 0.035 and 0.037. The corresponding average strains at onset of failure for the set Elem1 for each mesh size are summarized in Table 7.

No significant trend was observed and convergence was obviously not attained. An error of 31% between the maximum and minimum values of the strain was noted. To resolve this problem, a local structured mesh of brick elements with reduced integration (C3D8R) was generated in the critical region near the crack tip while the global mesh type was kept as tetrahedral elements C3D4 [Figure 16].

Figure 16. Mesh for load case 1 consisting of structured brick elements near the crack

Again, a mesh sensitivity analysis was performed for the mesh sizes: 0.015, 0.022, 0.024, 0.025, 0.027, 0.03, 0.035, 0.04 and 0.05. For each mesh size, the corresponding average strain at onset of failure was recorded with the results summarized in Table 8.

Finally, a converged mesh has been obtained with a maximum error of 12% between the maximum and minimum values of the strain. A mesh size of 0.025 (124,261 elements) was chosen, taking into consideration both accuracy of the results and the required simulation time.

In loading case 2, in order to avoid the difficulties faced in loading case 1 particularly the possibility of nodes connections forming new elements instead of being deleted at the crack site, linear tetrahedral elements (C3D4) were used to mesh the whole model. For the mesh sizing, a mesh convergence study was carried out considering the following sizes; 0.02, 0.025, 0.03 and 0.035.

For each mesh size, the corresponding average maximum principal stress at the surface of the barb when a crack started to propagate (BarbElem) was recorded [Table 9].

A converged mesh has been attained with a maximum error of 8% between the maximum and minimum values of the maximum principal stress and a mesh size of 0.025 (108,708 elements) was chosen for the same reasons as in loading case 1 [Figure 17].

Figure 17. Mesh for load case 3 consisting of linear tetrahedral elements

Results with variant aspect ratio, constant cut angle & constant cut depth

At first, the aspect ratios were varied for a constant cross-sectional area equivalent to that of a circle having a radius of 0.3-mm, a constant cut angle (CA) of 154 and cut depth (CD) of 0.19-mm for both loading conditions. The second step was to investigate the effect of varying CA and CD for a constant.

With constant cut angle and cut depth, in the loading case 1, an ENCASTRE and a displacement/rotation boundary conditions were applied to the suture geometry as thoroughly described in the boundary conditions and loads section. Stresses and strains in the z-direction, S33 and NE33 respectively, in the body of the suture (Elem1) were reported at the onset of crack propagation and compared among the different cross-sections [Table 10].

So, changing the aspect ratio can affect the strength of the suture. The highest strains and stresses are observed for aspect ratio 1/3. As such, ρ = 1/3 is the strongest, followed by 1/2, 3, 2, and finally 1. A circular suture is found to be the weakest among all configurations. However, since the differences in the magnitudes of NE33 and S33 are not significant, the results cannot be considered conclusive. No clear trend is detected.

In this circular cross-sectional model, the direction of propagation of the crack is almost perpendicular to the applied displacement and the rupture failure is spotted. Similar behaviors are observed for all other aspect ratios [Figures 18 and 19].

Figure 18. Initial situation of the suture under load case 1

Figure 19. A propagating crack in a circular suture subjected to load case 1

In loading case 2, the results reveal an evident trend [Table 11]. As the aspect ratio decreases below 1, the barb strength increases. The opposite is noticed for increasing the aspect ratio beyond 1. Therefore, the strongest geometry belongs to ρ = 1/3, followed by 1/2, 1, 2, and finally 3. Notice that the area of the barb remains constant despite the variations in both major and minor axes of the ellipses. This is anticipated since the criterion for determining the dimensions of the ellipses was to have the same remaining area, or alternatively, equal areas removed.

The direction of crack propagation for a circular cross-section model is being almost parallel to the applied displacement. Under such boundary conditions, peeling rather than rupture failure is detected. All other aspect ratios behave in a similar manner [Figure 20].

Figure 20. A propagating crack in a circular suture subjected to load case 2

The results obtained for constant cross-sectional area, constant cut angle and cut depth for both loading cases 1 & 2, cannot offer a decisive conclusion regarding the best cross-sectional configuration. A mere 7% increase in %NE33 is observed when comparing the weakest suture (ρ = 1) with the strongest one (ρ = 1/3). Also, it is not clear which condition leads to an increased suture strength, decreasing the aspect ratio below one or increasing it. On the other hand, a 30% rise in the barb force is detected for LC2 when changing the aspect ratio from 3 to 1/3.

So since the aspect ratio ρ of the tendon is 2, a suture with an aspect ratio of 1/2, given the findings up till this point, appears to be most suitable. Furthermore, since the smallest major and minor axes among the tendons are a = 1.5-mm and b = 0.858-mm respectively, the radius (or the equivalent ellipse axes) of the suture can be increased from 0.3 mm to a larger diameter that is suitable with both the dimensions of the smallest tendon and the optimum aspect ratio. It is desired that the optimized suture accommodates about 30% of the tendon’s cross-sectional area. Consequently, a suggested geometry has the following characteristics: ρ = 1/2, a = 0.25-mm and b = 0.5-mm.

Results with constant aspect ratio, variant cut angle & variant cut depth

The results of a constant suture aspect ratio , with different cut angles and cut depths are reported below. A circular cross-section was used with a radius of 0.3536-mm equivalent to an ellipse having a = 0.25-mm and b = 0.5-mm. While the cut angles remain the same, the cut depths were adapted for the ellipse [Table 12].

In loading case 1, the boundary conditions were applied to a cylinder with a circular cross-section of r = 0.3536-mm for the cut angles 150, 154, and 160, and cut depths 0.07-mm, 0.12-mm and 0.19-mm. The outputs of interest are the same as before for the same load case [Table 13].

For a constant cut angle, decreasing the cut depth leads to stronger sutures that fracture at higher average strains. This result is expected since a smaller cut depth is synonymous with a smaller crack. On the other hand, for a constant cut depth, increasing the cut angle results in stronger geometries. From this analysis, the best configuration is CA = 160 and CD = 0.07-mm.

In loading case 2, the analysis was performed for a circular suture of r = 0.3536-mm for the cut angles and cut depths 150, 154, and 160, and 0.07-mm, 0.12-mm and 0.19-mm respectively. The outputs of interest are the same as before for the same load case [Table 14].

For a constant cut angle, decreasing the cut depth results in a weaker barb that can withstand lower forces. On the other hand, for a constant cut depth, increasing the cut angle results in stronger barbs. This is expected, for when both CA and CD increase, the barb area increases as well, allowing it to endure higher forces. From this analysis, the best configuration is CA = 160 and CD = 0.19-mm

Optimized geometry

With respect to aspect ratios, it was concluded that ρ = 1/3 presents the strongest design as compared to other aspect ratios. Taking into consideration the tendon’s cross-section, it was decided that ρ = 1/2 is most suitable.

Concerning cut angles and cut depths, it was found from the load case 1, that the best configuration corresponds to CA = 160 and CD = 0.07-mm. For the second test case, CA = 160 and CD = 0.19-mm were found to offer the best results. As such, a cut angle of 160 is evidently superior. However, a decision on the cut angle is not direct. A middle value between 0.19-mm and 0.07-mm was chosen; that is, 0.12-mm.

Therefore, the optimum geometry has an aspect ratio ρ = 1/2, CA = 160 and CD = 0.12-mm. Load cases 1 and 2 were performed again for this optimized geometry [Table 15].

Comparing the results with those obtained for a circular cross-section [Table 16] for the same cut angle and cut depth confirms the conclusion that aspect ratio 1/2 is superior. The average strains as well as the forces at the surface of the barb at onset of crack propagation increase by 12%, and 11% respectively.

Since the values for both cut angle and cut depth were varied for a circular cross-section, the observed trend was to be validated for ρ = 1/2. The test was performed once for a constant CA of 160 and two different values for CD: 0.19-mm and 0.12-mm, and another time for a constant CD of 0.19-mm while changing the CA from 160 to 154.

For a constant cut depth, decreasing the cut angle resulted in a weaker suture and a barb that can tolerate lower forces (decrease of 14% in %NE33 and 33% in the force), as observed for ρ = 1 [Table 17]. For a constant cut angle, decreasing the cut depth resulted in a stronger suture (24% higher %NE33) but a weaker barb (120% lower barb forces). This was also noted for a circular cross-section.

A loading case 3 was performed for the optimized cut angle of 160and cut depth of 0.12-mm, but for ρ = 1, since generating a working mesh for ρ = 1/2 was problematic under the prescribed boundary conditions.

Reporting the findings for this test case, a small difference was noticed in comparison with the second loading condition. A force of 6.26 N was obtained as opposed to 6.93 N (10% difference). The same trend is expected for ρ = 1/2 but most likely, the percent variation differs. However, this decrease in barb strength would be accounted for in safety factors when determining the number of barbs required to be anchored within the tendon [Table 18].

As aforementioned, an estimation of the number of barbs required to hold a tendon force equivalent to 40 N can be simply found by dividing this force by the one withstood by a single barb.

For ρ = 1/2, CA = 160, CD = 0.12-mm, the force at the barb surface was recorded to be 7.71 N. This value represents the force at which crack propagation initiates. Certainly, the barb should be loaded at forces lower than the threshold. And since for the third loading case, a value lower than 7.71 is anticipated, the force at the barb surface is assumed to be 5 N, 35% less. Taking into consideration a factor of safety of 3, the required number of barbs to be anchored in the tendon to withstand 40 N becomes:

Hence, 24 barbs are required. Note that this number is a mere estimate.