A hyperelastic or green-type elastic material [13–15] is a type of material that follows constitutive models where the stress–strain behavior is defined by a strain energy density function [16]. They are considered to be truly elastic as they store energy during loading and dissipate equal amounts of energy during the unloading process. These materials experience large strains that are mostly recoverable [17–19]. To define hyperelasticity, many mathematical models have been developed and different aspects of non-linear elastic material behavior have been explained [20, 21]. These types of models have been very successfully applied to soft materials and tissues.

Polyvinyl alcohol cryogel (PVA-C) samples of 5%, 10%, and 15% concentrations are modeled as soft tissue mimics to examine their non-linear material properties. Experimental force–displacement (F–D) data were used as an input parameter for the Ogden hyperelastic model.

A hyperelastic material is a type of constitutive model where the presence of a strain energy density function is assumed, and the stress–strain relationship is derived from a strain energy density function. The Ogden hyperelastic model for isotropic material can be obtained from strain as follows:

where μ, α, and K are presented as constitutive parameters, J is the determinant of the strain tensor, and λ is known as the principle stretch. For incompressible material deformation, J = 1 leads the second term in the aforementioned equation to vanish. For the first-order Ogden model, n = 1. Thus, the model will have only two unknown parameters of shear modulus (μ) and strain hardening exponent (α).

For a uniaxial compression test, the nominal stress (σ) is represented as a function of the stretch ratio λ. The first-order Ogden material model can be presented in the form as given in equation (2):

The aforementioned two equations can take only positive values. The experimental force–displacement (F–d) data are used as an input parameter for the hyperelastic model. By applying curve fitting to the experimental (σ – λ) data, the first-order Ogden parameter is extracted.

Commercial FE software package Abaqus was used to create an axisymmetric model for examining indentation. PVA-C was used as a soft tissue mimic, and a flat-ended cylindrical indenter was modeled in this simulation. The contact surface between the indenter and the soft tissue was set as frictionless. Soft biological tissue is generally considered to be incompressible [22]. For optimal numerical accuracy, high mesh density was adapted underneath the indenter, and convergence criteria were verified. Indentation tests were simulated under linear and non-linear hyperelastic model assumptions. Numerical analysis was performed by FEM using Abaqus version 6:13–4 (2013). Cylindrical indenters and rigid flat-ended cylindrical indenters were used for numerical analysis. An indenter of radius 4 mm was indented on a soft tissue-mimicking sample of PVA-C with 5%, 10%, and 15% (w/w) concentrations. The sample size was L = 19 mm and B = 12 mm. The soft tissue sample was modeled as homogeneous, isotropic, and nearly incompressible.

An axisymmetric model was developed and meshed with (CAX4RH) and (CAX3H) linear quadrilateral elements. A very fine mesh was made at the contact zone and comparatively, the coarse mesh was utilized outside the contact zone. The boundary conditions for the three sides of the sample were fixed (U1 = U2 = U3 = UR1 = UR2 = UR3 = 0), and one side is variable (U1 = U3 = UR2 = 0). During the indentation test, a 4N vertical load was applied, and the contact between the indenter and the sample was considered frictionless. The Ogden hyperelastic material model was considered for PVA-C samples undergoing the indentation test, the data of which are provided in Ref. [23]. The simulation was completed in 13 steps with a step size of 0.01. To observe the indentation responses, the hyperelastic FE model was analyzed including a sample with finite thick and infinite thickness made of PVA-C with 5% and 15% concentrations.

A two-dimensional axisymmetric cylindrical indentation model was developed by using Abaqus version 6:13–4 (2013). The Ogden first-order strain energy density function was used in the numerical simulation. In total, 7818 (CAX4RH) nodes were generated in the meshing process, as shown in Table 1.

An axisymmetric FE model was developed by using Abaqus, and PVA-C 5% and 15% 2FTC experimental data obtained from Ref. [23] were used as an input parameter for the first-order Ogden hyperelastic model. The load–displacement graph for both PVA concentrations is shown in Figures 1, 2.

Based on previous research, it has been found that the hyperelastic parameters determined using cylindrical or spherical indentation testing are not always the same as those obtained from experimental uniaxial compression tests [4, 6, 10, 20]. To investigate the sources of such disagreement, a parametric study was conducted to examine the influence of (μ, α) on the shape of the simulated data. An optimization algorithm combined with MATLAB was used to minimize the sum of the squared difference between the experimental measurements and the FE-simulated Ogden hyperelastic model. This allowed for the unknown parameters (μ, α) to be determined. A sensitivity analysis was then conducted to verify the accuracy and robustness of the parameters. For cylindrical uniaxial indentation obtained from inverse analysis shown in Table 2.

The inverse FE method was used to determine the hyperelastic material parameters of soft tissue specimens from indentation responses [24, 25]. Based on Hadamard's definition, an inverse problem is posed if one of the three conditions is violated: (i) existence, (ii) uniqueness, and (iii) stability. Many optimization algorithms are commonly used as such least square fitting power low, Levenberg–Marquardt (LM) Trust Region Algorithm, and Kalman Filter to solve the inverse problems.

The inverse analysis is introduced to minimize an objective function with respect to unknown constitutive material parameters (μ, α) that match the experimental data [10, 26]. The Levenberg–Marquardt (LM) method was used in this dissertation to extract the unknown parameters based on the inverse analysis. The LM method is defined as the minimization of the error function Φ with respect to a vector P^. The error function is represented as follows:

Here, P88 is a vector that contains unknown constitutive parameters  P^T= {μ, α}, and n is the number of measurements. The vector r^  is defined as follows:

where t^* and t^ are the model-predicted and experimental data.

A detailed flowchart of the inverse optimization process to characterize the tissue hyperelastic parameters used in the current study is shown in Figure 3.

First, experimental (F–D) data were used as an input parameter for the numerical model (Abaqus). Simulated (F–D) data were then used as an input parameter for inverse analysis (MATLAB). An object junction was introduced to minimize the quadratic difference between simulated- and model-predicted data. Through the optimization process, a numerical convergence was achieved, and an optimized unknown parameter was obtained. These parameters were used as an input parameter of Abaqus for numerical validation.

Although the implicit method with Newton–Raphson iterative solver is enough to obtain the converged solution, mesh optimization through adaptive meshing was adopted for numerical accuracy. Optimized μ and α through mesh convergence are shown in Figures 4, 5.

The Ogden model is used to identify the hyperelastic material properties. The Levenburg–Marquardt (LM) algorithm was used to minimize the difference between experimental- and model-predicted data. Here, a novel approach is introduced to determine and validate material hyperelasticity. This method includes the following steps:

Through inverse analysis, extracted Ogden parameters are readily available to use in simulation and also to verify the stability of the solution. The detailed process is shown in the flowchart of Figure 3.

A new technique has been developed to characterize the biomechanical properties of nonlinear material using Ogden and hyperelastic models. Experimental results are compared with the novel inverse technique that can be further investigated to develop patient-specific artificial organs.

FE-simulated cylindrical load–displacement data were used in the robust optimization algorithm through inverse analysis. The entire procedure is shown as a flowchart in Figure 3. A validation exercise was conducted using Ogden parameters (μ), which were varied while α was kept constant and vice versa as shown in Figures 6, 7 [10].

The results obtained from the load–displacement curve indicated that there is a significant effect of (μ), as this parameter explains the strength of the material and α is the strain hardening coefficient. In the process of validation, μ and α values were compared with that of Ref. [10] and found good agreement.

A detailed flowchart of the inverse optimization process to characterize the tissue hyperelastic parameters used in the current study is shown in Figure 11.

First, experimental (F–D) data are used as an input parameter for the numerical model (Abaqus). Simulated (F–D) data were then used as an input parameter for inverse analysis (MATLAB). An object junction was introduced to minimize the quadratic difference between simulated and model-predicted data. Through the optimization process, a numerical convergence was achieved and an optimized unknown parameter was obtained. These parameters were used as an input parameter of Abaqus for numerical validation.

The extracted Ogden parameters from the inverse optimization algorithm are plotted against various concentrations of PVA-C 5%, 2FTC and PVA-C 15%, 2FTC as shown in Figures 12, 13.

It can be hypothesized that (μ, α) values have a greater dependency on thickness. As the thickness increases, both μ and α increase.

The freeze and thaw technique is a part of the stability testing that determines whether any formulation will remain stable under various conditions. The extracted Ogden parameters from the inverse optimization algorithm are plotted against concentration and freeze–thaw cycles are shown in Figures 14, 15.

This led us to conclude that there is a proportional relationship between (μ, α) values and PVA-C concentration and freeze-thaw cycles. With the increase in PVA-C concentration, (μ, α) values increase. With the increase of the freeze–thaw cycle time (FTC), the values of (μ, α) also increase.

To overcome the stability problem, sensitivity analysis was conducted by adding noise (±1% and ±2 %) to the solutions. This noise modulation with the input data will enable us to investigate whether there is any extraneous influence on material and equipment. Fellay et al. [6, 13] considered this problem as another minimization approach.

After determining optimal material parameters, the sensitivity analysis of the solution was analyzed by adding noise (±1% and ±2%) to the iterative solutions. Based on Hadamard's [9] article, a solution to any inverse problem continuously depends on the stability of the data [28]. The sensitivity of the solution indicates that there was no external influence in the solution. This again proved the uniqueness of the solution.

The objective of the study was to characterize the nonlinear behavior of soft tissue phantom. The combination of the inverse method in conjunction with the FE method enables the identification of unknown material parameters. The LM optimization algorithm was used to optimize material properties by minimizing the sum of squared differences between the model-predicted and experimentally measured load–displacement data, which provided benchmarks for accurate Ogden parameters (μ, α). The accuracy, effectiveness, and robustness of such procedures were validated through FE-simulated data, which is cross-checked by the LM optimization technique and finally, compared with published results.

First, the simulated data were used as an input parameter for the MATLAB optimization algorithm. The robust optimization technique was performed as expected, which confirmed the uniqueness of the solution. Second, the Ogden parameters were plotted at various thicknesses. This proved that μ and α have a significant effect on specimen thickness. Moreover, μ and α were plotted against concentration and freeze–thaw cycles. The simulated result confirmed the dependency of μ and α on PVA concentration and freeze–thaw cycles. This means that the value of μ and α increases with the increase in concentration and freeze–thaw cycles with the little exception of PVA at higher concentrations.

R-squared values for 5% and 15% PVA-C, 2FTC were obtained through the parameter estimation techniques, and their fitness coefficients were 0.99999 and 0.9998, respectively. This led us to conclude that the optimization algorithm provided accurate results in the determination of the Ogden parameter. Finally, PVA-C 5% and PVA-C 15% have undergone a sensitivity test by adding y ± 2% noise to the solution. No extraneous effect was observed. Thus, supporting the uniqueness of the solution. This robust technique will be proposed as a “gold standard” for future biomedical research.

Finally, some perspectives and future works of this study have been summarized in the Conclusion section.