1 Introduction

Soft biological tissues are ensembles of cells of various types and their extracellular matrixes that together have a specific functionality. Biological tissues are among the most complex materials in terms of structural heterogeneity. This complicates the characterisation of their mechanical behaviour that depends strongly on the morphological features and mechanical behaviours of the underlying tissues and cell types. A human soft tissue consists of skin, adipose (fat), muscle, veins, each with its own substructure. The result of this complexity is an anisotropic mechanical response [1] varying throughout the tissue with the local morphology. Soft tissues are known to exhibit highly nonlinear [2] and time-dependent behaviours [3]. Due to high water concentration, soft tissues can be considered incompressible. Typically, soft tissues are characterised by their J-shaped stress-strain curves [3] with stress remaining (mostly) proportional to strain at initial small deformation, followed by a steep increase after reaching a certain deformation limit.

Previous studies into mechanics of soft tissues tried to obtain these curves experimentally and characterise the soft tissue both in vivo [4] and in vitro [5]. Since there is no standard sample-preparation protocol for such tests, both intrinsic and extrinsic factors can greatly influence the testing results of in vitro studies. In such studies, the tissues are usually analysed in extreme conditions exploring their maximum stretch or compression that are impossible to study in vivo. In in vivo experiments, indentation is often employed as a non-invasive method to measure the soft tissue’s response [6]. It is a common practice for an indentation depth to be limited by a pain tolerance threshold for in vivo experiments [7], with cases of pain tolerance explored as well [8]. Although the reported pain levels vary significantly among participants due to the individuals’ characteristics, similar regions exhibit similar behaviours.

The main challenge of studies of human tissues in vivo is the necessity of involvement of human participants. Mostly, conventional experimental equipment is not suitable to carry out the measurements needed for the studies, so custom testing rigs and sensors are usually developed [4, 9, 10]. This equipment needs to produce accurate results while complying with strict safety regulations, so that the participants are not in risk during the examination. Ethical regulations for the methodology also should be followed. Hence, in vivo studies can be a very challenging process. Computational modelling can be a useful substitute as it provides tools to test hypotheses and predict the mechanical behaviour without the need for human involvement, thus eliminating the risk for accidents. It also constitutes an ideal method for parametric studies based on multiple repetitions needed to study the effects of various parameters and their combinations.

To predict the soft-tissue behaviour, a finite-element method (FEM) has been extensively used. Both 2D and 3D models were developed to capture spatial distributions of deformations and pressure in the tissues [11]. Taking advantage of a smaller number of elements and the lower computational cost, 2D models were used to introduce various properties of the analysed tissues. The assumptions employed in this simplified modelling approach can significantly affect the modelling strategy and the obtained results thus requiring further detailed investigations [12].

Taking advantage of the computational power that become available in the recent years, bigger and more complex models were developed, including full-limb models with patient-specific anatomical and morphological features. These developments were supported by the use of magnetic resonance imaging techniques for better accuracy of the measurements and identification of individual tissue’s behaviour [13,14,15]. Stress distributions on the surface and inside the tissues were studied with the FEM together with interactions between tissues and external parts such as prosthetic sockets [16, 17].

Although the model geometry is becoming more and more complicated and patient-specific, still little is known about the material behaviour of soft tissues. Although describing their mechanical behaviours is an active field of research, with many variants of constitutive models reported in the literature. Initial approaches that used linearly elastic material models were quickly phased out due to their limited suitability for cases with large deformation and material’s nonlinearity [18]. Various hyperelastic models were proposed for more adequate description of tissue’s deformations. In these models, the strain energy of a solid is expressed as a function of the invariants or principal stretches, with the stresses calculated with the proper derivation. Such models are introduced in Table 1.

Table 1 Strain-energy functions of hyperelastic models [19]
Full size table

Identifying the in vivo material parameters for soft tissues is an important but challenging aspect of experiments. When testing rigs are available, an inverse engineering approach can be employed to fit the computed results to the respective experimental data [20, 21]. Despite a significant number of studies carried out in this field, there is still no established universal model developed for soft tissues, while the respective material parameters characterising their mechanical behaviour have significant variations. As a result, the existing models are mostly limited to subject-specific studies and local behaviour. A summary of the material parameters reported in the literature for limb’s soft tissues is presented in Table 2.

Table 2 Hyperelastic models used to model soft tissue in different schemes and values of their parameters as well as experimental data for elastic properties
Full size table

Various methods were employed to model the soft tissue. Initially, single-domain models (i.e. with one bulk tissue throughout its volume) were common, usually with linear elastic or Neo-Hookean material formulation of their mechanical behaviour [17, 66]. The schemes with multiple domains, presenting skin, fat and muscle tissues separately, were also developed; they mostly considered the skin and fat tissue as a joint domain due to their similar morphological character in clinical images. More recent studies investigated the necessity of including further features in the model, such as fascia, veins and individual muscle compartments, but there was no (or very limited) information about the response of these individual materials [23, 67].

This study, focused on a numerical approach, aims to investigate and identify a necessary size of a 3D FE model for simulation of indentation of a human lower limb. In addition, the effect of the chosen levels for material parameters was investigated after the optimal model size was identified. The effect of the loading direction was also assessed.

2 Development of FEA model

One of the challenges that arise when modelling human soft tissues is unique anatomical characteristics of every individual, as mentioned above. To keep analysis consistent, in this study, a contour of a human thigh was manually segmented from a transverse-cut medical image (Fig. 1) and used for all subsequent assessments. Contours of a bone and a dermal tissue of the modelled thigh were also identified and separately segmented using NX 11.0. A skin and an adipose tissue were segmented as single part, as reported in previous studies [67]. A final sketch of the limb including three types of tissue (bone, muscle, dermal), each as a separate domain, was then imported into ABAQUS 2018 CAE software for analysis.

The size of the thigh was scaled according to the average size of male’s European human thigh circumference as in Table 3, and the transversal cross section was extruded along the main (bone) axis, as in the previous approach [47]. The length of the extrusion was determined based on requirements of the subsequent analysis. Each tissue’s region (domain) was partitioned to introduce the respective material properties. The mesh was created using 4-node linear tetrahedral elements with hybrid formulation (C3D4H), and mesh convergence was achieved with an approximate node spacing of 3 mm at the outer site and 5 mm at the bone region. A cylindrical, flat-ended indenter with diameter of 10 mm, with an edge radius of 1 mm, was modelled as analytical rigid surface. To implement the indentation simulation, the bone section was fully constrained, while the indenter’s flat end was positioned in contact with the skin and allowed to translate axially towards the centre of the bone. The selected friction coefficient was 0.5 and used for all the models developed (see also Sect. 2.2).

2.1 Constitutive model

The soft tissues’ nonlinear behaviour in our simulations was described using the Ogden hyperelastic model for incompressible materials [72]. In an isotropic, first-order Ogden model, the strain energy W is given as a function of the principal stretches (see Eq. (2) in Table 1). In this equation, N represents the number of terms used in the strain energy function, \(\mu _{p}\) are components of shear moduli and \(\alpha _p\) are dimensionless constants. The classical shear modulus \(\mu \) can be obtained in the reference configuration from

$$\begin{aligned} 2\mu =\sum \limits _{p=1}^{\mathrm{N}} \mu _{p} \alpha _{p},\quad \mu _{p}\alpha _{p}>0, p=1,\ldots ,N \end{aligned}$$
(8)

while the three principal Cauchy stresses \(\sigma _{\alpha }\) are given by

$$\begin{aligned} \sigma _{\alpha }=-p+\sum \limits _{p=1}^N {\mu _{\alpha }\lambda _{\alpha }^{\alpha _{p}}} , \quad a=1, 2, 3 \end{aligned}$$
(9)

Material parameters selected for the two hyperelastic tissues used in this study are highlighted in bold font in Table 2 describing the average reported material stiffness. Linear elasticity was used for the bone, with the Young’s modulus of 20 GPa and the Poisson’s ratio of 0.3. All the modelled tissues were assumed to be homogeneous and isotropic materials.

Fig. 1
figure 1

Anatomy of human femur [68]

Full size image

A preliminary numerical study was undertaken to assess the effect of friction at the indenter/limb interface in a broad range of the coefficient of friction. Eventually, the magnitude of 0.5 was chosen for this parameter in the main study to emulate conditions of the dry-skin examination.

2.2 Study of model size

As well known, a type of boundary conditions has a significant impact on the accuracy of simulation results. While simulating indentation of a human limb can be captured relatively accurately with a larger size of the limb, it may not be necessarily viable for multi-model and/or parametric studies. Therefore, in this study, indentation characteristics at two different sites—Site A (with soft-tissue thickness of 47 mm) and Site B (87.5 mm) (Fig. 2b)—were evaluated with models with varying length: from 20 to 220 mm (Fig. 2a). The simulations were performed with both an unconstrained and a fully constrained boundary conditions applied at both transversal ends (full cross sections) of the limb. Convergence was assumed when the difference in the reaction force between the models with the unconstrained and fully constrained boundary conditions was smaller than 5%.

For Site A, with tissue thickness of 47 mm, convergence in the resulting reaction force was achieved for a model length of 80 mm, while for Site B (thickness 87 mm) this happened at 150 mm (Fig. 3)., A ratio of 1.7 between the minimum model length and the tissue thickness was found for both cases. Considering this ratio, the main study was implemented with two different models, each with the minimum identified size.

Table 3 Average size of thigh’s circumference of healthy adults for various geographical regions
Full size table
Fig. 2
figure 2

a Models with length of 20, 60, 150 and 220 mm. b Indentation sites used in model corresponding to minimum and maximum tissue thicknesses

Full size image
Fig. 3
figure 3

Reaction forces for models with various lengths and constrained and free ends for Site A (a) and Site B (b) (The arrow shows the convergence point)

Full size image

3 Methodology

3.1 Effect of loading direction

When performing and simulating indentation, ideally, the indenter’s axis should pass through the bone’s centre. Obviously, small misalignment can occur in real-life experiments as a result of the unique feature of individuals or a human error. A set of indentations was simulated for various loading directions on Sites A and B, to examine the sensitivity of the indentation protocol to the misalignment of the indenter and the significance of the internal geometry of the limb. Effects of both linear and angular misalignments were investigated. An angular misalignment of \(10^{\circ }\) was considered for both clockwise and counter-clockwise directions. On the thinner Site A, the indentation axis was still meeting the bone tissue for the entire range of this misalignment, while on the thicker Site B, the axis was off the bone. To study the effect of linear misalignment, the indenter was repositioned at \(\pm \,5\), \(\pm \,10\) and \(\pm \,20\) mm parallel to the original position. The original and new indentation directions are shown in Fig. 3. For this study, the models optimal identified for each site (as discussed above) were used.

Fig. 4
figure 4

Angular misalignments on Site A (a) and Site B (b). Linear misalignments on Site A (c) and Site B (d)

Full size image
Fig. 5
figure 5

Sets of material parameters used for muscle domain in simulations (Set R represents the reference values; \(\mu \) is in kPa.)

Full size image

3.2 Effect of material parameters

In this study, the soft-tissue response was investigated for various sets of material parameters for the first-order Ogden hyperelastic model (see Sect. 2.2). A broad range in properties’ magnitudes available in the literature necessitated the assessment of the respective effect. In our numerical simulations, changes in material properties were limited to the domain of muscles that constituted the largest part of the modelled soft tissue, whereas those of the skin and adipose were kept at the same level. In reference to the initial set of material parameters—Set R in Fig. 5, additional sets were created to explore the changes in the limb’s response to indentation. Two additional levels were introduced for \(\mu \), the first equal to one half of the reference value and the second was three times higher than the reference. For \(\alpha \), a nearly fivefold decrease of the reference magnitude was used, with the total range being from 4.5 to 21.4. Permutation of these magnitudes resulted in five additional sets of material parameters presented in Fig. 5. So, Sets 2 and 4 have the same level of \(\alpha \) as Set R, while differing in the magnitude of \(\mu \): Set 2 has the lowest one, while Set 4—the highest. Sets 1, 3 and 5 share the same—lowest—magnitude of \(\alpha \), with Set 3 retaining the reference level of \(\mu \), while Sets 1 and 2 having the same lowest level of \(\mu \) and Sets 4 and 5—highest \(\mu \).

Levels of the reaction force at the indenter tip, obtained in numerical simulations, were compared for all the sets, with their ratios used to evaluate the contribution of each parameter. The aim of this part of the study was to identify the range of the expected reaction forces for the specific choice of material parameters and to assess the ability of this formulation to capture the tissue response during indentation.

3.3 Axisymmetric model

An axisymmetric finite-element model was also developed, with the same geometrical features of the limb tissue. A rectangular domain with width of 40 mm and height of 47 mm, with the latter coinciding with the symmetry axis, was created for simulations of indentation at Site A. Its bottom edge was fixed in all degrees of freedom, representing the tissue attached to the bone, while the edge of the tissue parallel to the symmetry axis was left free. Indentation to the depth of 20 mm was modelled, using an indenter identical with that in the 3D limb model. To elucidate the mechanical behaviour of the soft tissue, it is necessary to assess the contribution of the bone to the stresses so that the measurements present the data for the material of interest rather than a combined data for the limb. The minor principal strain was chosen for comparison as it describes compressive strains of the material under the indenter tip.

Fig. 6
figure 6

Evolution of reaction force against compressive stretch ratio (ratio of indentation depth to tissue thickness): effect of angular misalignment for Site A (47 mm) (a) and Site B (87.5 mm) (b); effect of linear misalignment for Site A (c) and Site B (d)

Full size image

4 Results

4.1 Effect of angular misalignment

The results calculated for the cases of angular misalignment presented significant differences. In both sites examined, the difference in reaction forces for the misaligned directions and the reference exceeded the 5% threshold used throughout the simulations. Specifically, a \(\pm 10^{{\circ }}\) misalignment of the indenter tip at Site A resulted in 25.3% (\(+\) 10%) and 14.4% (− 10%) reductions in the peak reaction force, respectively, with a slight difference as a result of non-symmetric limb profile. For Site B (thicker), the resulting force was 8.6% (\(+\) 10%) lower and 1.4% (− 10%) higher than the original value, (Fig. 6a, b). The results indicate that the reaction force measured in indentation is sensitive to the alignment and structure of the underlying tissue, while the thicker layer of the tissue helps to mitigate some of the misalignment effect.

Table 4 Difference in reaction force (%) for various indentation directions at maximum indentation depth 20%
Full size table

4.2 Effect of eccentricity

The study of the effect of eccentricity to the indentation direction demonstrated the growing distance from the ideal direction (i.e. to the bone centre) caused a monotonous decrease in the reaction force. The maximum decline in this force measured for a 20 mm offset at Site A exceeded 20% for both directions (Fig. 6c). For smaller offset positions, the positive offsets, at \(+\) 5 mm and \(+\) 10 mm, seem to have less effectiveness compared with the negative offsets.

The extent of reduction in the maximum reaction force at site B was generally lower. Besides, the diminishing trend was observed only for the positive offset distance, while for the negative one it was inversed (Fig. 6d and Table 4). The difference between these two cases was attributed to a rapid change in the thickness and curvature of the soft tissue as the offset changed from − 20 to \(+\) 20 mm.

4.3 Effect of material parameters

Since both constitutive parameters of the first-order Ogden model affect directly the shear modulus, an increase in either of them resulted in an increased reaction force, as was demonstrated by numerical simulations. Still, as a result of their different roles in the constitutive equation (see Eq. 9), contribution of \(\mu \) was noticeable throughout the entire modelled stretching. On the other hand, parameter a contributed more to the response after the compressive stretch ratio exceeded 0.1 for both Site A and Site B (Fig. 7). Thus, \(\mu \) defined the slope of force-stretching curves at the early stage of the process.

Fig. 7
figure 7

Evolution of reaction force for studied material sets for Site A (47 mm) (a) and site B (87.5 mm) (b)

Full size image
Fig. 8
figure 8

Comparison of force–displacement curves for 3D model (Site A) and corresponding axisymmetric model for reference material parameter set (Set R: \(\mu = 15.6\) kPa; \(\alpha _1=21.4\))

Full size image

Analysis of the results for Site A (Fig. 7a) demonstrate that changes in \(\mu \) resulted in almost proportional variations in the reaction force for compressive stretch ratio above 0.15. Specifically, at compressive stretch ratio of 0.2, the ratio of reaction forces for Set 2 and Set R (i.e. the same \(\alpha \)) was 0.58, while that for Set 1 and Set 3 (i.e. the same, but lower \(\alpha \)) was 0.61. For the same loading conditions, the ratio of reaction forces for Set 4 and Set R (at constant \(\alpha \)) was 2.56, while that for Set 5 and Set 3 was 2.68. A change in the magnitude of \(\alpha \) from 4.5 and 21.5 for the same level of \(\mu \) returned continuously a reaction force around 1.5 times larger at stretch ratio of 0.2.

Similarly, the relationship between the reaction forces for various sets of material parameters for Site B (Fig. 7b) followed the same trends as in Site A. The ratio of reaction forces for Set 2 and Set R was 0.53 at compressive stretch ratio of 0.2, while that for Set 4 and Set R was 2.7. For the lower level of \(\alpha \), the ratios for forces for Set 1–Set 3 was 0.55 and for Set 5 to Set 3 was 2.65, respectively. At compressive stretch ratio of 0.2, the variation of \(\alpha \) resulted again in the approximate ratio of 1.5 for forces in materials with the higher level of \(\alpha \).

4.4 Comparison with axisymmetric model

A detailed comparison between the 3D limb model and an axisymmetric one was made using the reference set (Set R) of material parameter from Fig. 5. For direct comparability, both models had the same ratio of the length to the tissue thickness (equal to 1.7) identified earlier in the study of the model size (Sect. 2.2). The boundary conditions at the edges of the models were free for all six degrees of freedom since the stress concentration disappeared fully at these dimensions.

A good agreement was revealed, notwithstanding differences in geometry, boundary conditions and domain shapes. Good accordance for evolution of the reaction force with indentation depth was observed for the entire loading process as shown in Fig. 8. In addition, distribution of a principal compressive strain at 10 mm and 20 mm indentation displacement were also compared (Fig. 9), with good compliance both for the character of distribution and the magnitude. So, despite the existence of a continuous dermal layer in the three-dimensional limb model, the material response could be captured highly accurately with the simplified axisymmetric model during the entire process of indentation.

Fig. 9
figure 9

Comparison of principal logarithmic strains for 3D full-limb and axisymmetric models for Site A (47 mm) for various indentation depths: a 10 mm; b 20 mm

Full size image

5 Discussion

It is evident that a small misalignment of the indenter can result in significant differences in the resulting force–displacement curves, and, therefore, influence the quality of experimentally obtained data. So, when performing an indentation experiment using a flat indenter tip to extract force–displacement data with an inverse analysis, it is important to simulate the exact position and direction of the indenter, instead of using idealised boundary and loading conditions. A recent study on rat tibia [27] concluded that tests with a spherically shaped indenter along various indentation directions altered the strain distribution within a defined volume of tissue.

However, reaction forces of the primary measurement were not reported, so further experimental investigation is necessary to assess differences between their testing protocol and the idealised modelling.

As demonstrated, the changes in material parameters of the first-order Ogden hyperelastic model significantly influenced the evolution of reaction force with the indentation depth. The extent of the effect of material properties varied for specific regions. In order to elucidate the contribution of each parameter, the ratios of reaction forces for the same sets of \(\mu \) (or \(\alpha \)) with different levels of \(\alpha \) (or \(\mu \)) were compared (Fig. 10) for both Site A and Site B using the sets of material parameters introduced in Sect. 3.2 (Fig. 5). The effect of changing the shear modulus was evident throughout the entire deformation process, starting at around 1.4 at the initial loading stage, but quickly saturated to around 2.6 for the ratio \(\mu = 3\) (1.6 for \(\mu = 2\)) at larger indentation depths (Fig. 10a, solid lines). In contrast, the contribution of \(\alpha \) was mostly negligible up to the indentation depth of some 10% of the initial thickness of the soft tissue (Fig. 10a, dashed lines). Interestingly, it was also revealed that these ratios were mostly unchanged for Site B, suggesting thickness-independent tissue responses.

Fig. 10
figure 10

Ratio of reaction force for various sets of material parameter (see notation in Fig. 5) at different indentation depths: a Site A; b Site B. Solid lines are used for comparisons of sets with the same \(\alpha \); dashed lines—with the same \(\mu \)

Full size image

It is important to recall at this point that since only the muscle material parameters were changed, the case with Set 1 can effectively be considered as single-material tissue since the muscle properties were almost identical to those of the layer of skin and adipose. Although the reference material, Set R, is stiffer than that corresponding to Set 2 (Fig. 7), no significant differences in the behaviour was observed as a result of introduction of the skin/adipose domain. The response followed a similar trend although with a different magnitude.

The model was further compared with specially performed experimental indentation tests of a human calf with the same range of tissue thicknesses (its detailed experimental protocol and results will be discussed in a separate paper). Although this comparison is suitable only for preliminary conclusions, it was evident that the model can reproduce the experimentally observed behaviour of soft tissues at large indentation depths (Fig. 11). A significant extent of variability in mechanical behaviours of lower-limb soft tissue is obvious, especially for the cases of tissue with the same thickness (58.8 mm) in different persons, proving the difficulty of studying these materials. The comparison of the obtained numerical results with the experimental response suggested that many modelling assumptions resulted in overestimation of the soft-tissue stiffness. Considering more recent studies focusing on measuring mechanical properties [57,58,59,60,61,62,63,64,65], a lower range of elastic moduli was reported (see Table 2). It is thus important to consider lower values for material parameters that correspond to lower elastic moduli for thick muscles. These values should be also considered as a reference point in future modelling attempts.

Fig. 11
figure 11

Comparison of experimental data for indentation of calves with different thickness (solid lines; the legend provides the soft-tissue thickness in mm) and simulation results (lines with markers)

Full size image
Fig. 12
figure 12

Distribution of normalised stress in indentation direction along tissue thickness

Full size image

Although the indentation experiments were used extensively to characterise the in vivo soft tissue behaviours of the lower limb, there were only a few studies that investigated the distributions of stress and strain in the tissue directly beneath the indenter. Also, there was no information regarding the reaction contribution of the rigid support to the tissue’s response in the indentation process. Furthermore, computing hyperelastic material properties based on the assumption of a quasi-statically, uniformly deformed tissue would be an oversimplification. To elucidate the stress development directly beneath the indenter tip, the normalised (with its peak magnitude) stress in the direction of indentation for different indentation depths (10 mm, 20 mm and 30 mm) was plotted along the normalised tissue thickness for Site A in Fig. 12. The results demonstrated a significant variation of stress along the indentation direction. For indentation depth of 10 mm and 20 mm, less than half of the thickness was affected by stress in excess of 50% or the peak stress concentrated at the top of the tissue. On the contrary, for a higher indentation depth (30 mm), the stress concentration band at the upper tissue was much larger, with some 90% of the thickness exposed to higher stresses (i.e. exceeding 50% of the peak stress). This is an indication that the properties measured at this depth were affected by the rigid support (bone), suggesting that the maximum indentation depth for studying this tissue should be limited to 0.45 of its original thickness. A further study on the link between the material stiffness and the indentation depth is currently under progress.

6 Conclusions

In this paper, various modelling strategies for numerical simulations of indentation on human lower-limb soft tissue were suggested and evaluated. A consistent minimum size of the model was established based on the ratio between the model length and the tissue thickness, so that the simplified boundary conditions in 3D and axisymmetric models can be used. For analysis of hyperelastic materials and the indentation depth up to 20% of the initial thickness of soft tissue, a ratio of 1.7 between the model length and the tissue’s thickness can be assumed to satisfy sufficiently the Saint-Venant’s principle.

It was also evident that the indentation reaction force was sensitive to the indentation direction and the indenter’s initial position with respect to morphological features of the tissue (e.g. the bone shape). A further experimental study is needed to accurately quantify this effect; obviously, the precise alignment of the indenter in experimental procedures is advisable to avoid potential errors.

The findings of this study further revealed that the mechanical response to indentation of the lower-limb soft tissue was affected differently by the two parameters of the first-order Ogden hyperelastic material model. The shear modulus is proved to contribute significantly at all the stages of the material’s deformation, while \(\alpha \) was mostly contributing to the nonlinear behaviour for deeper indentations. The ratios material parameters for different sets indicated that the difference in reaction force due to changes in shear modulus remained constant once the indentation depth reached around 10% of the soft-tissue’s initial thickness. For the inverse analysis, it is suggested that hyperelastic materials described by this formulation could have the material parameters identified individually, reducing the computational cost of multi-parametric optimisation. The shear modulus can be identified by shallow indentation tests (\(<10{\%}\)) for a given \(\alpha \). Once the slope of the curve is identified, the level of \(\alpha \) can be calibrated for the obtained shear modulus. Assessment and comparison of the calculated distributions for compressive stress and strain under the indenter for 3D and axisymmetric models suggested that hyperelastic materials can be characterised reasonably accurately for stretch ratios up to 0.45, without being affected by reaction of the rigid support (bone).

In conclusion, this study indicated the parameters necessary for modelling the indentation processes of human soft tissue. The results—also compared with experimental data—agree with the literature but also indicate that further studies should be carried out to advance the existing modelling strategies. There are also opportunities to further enhance numerical models by employment of material formulations for soft tissues, accounting directly for the effect of fibres (see e.g. [74,75,76]).