The time-dependent properties of rubber-like synthesized and biological materials are crucial for their applications. Currently, this behavior is mainly measured using axial tensile test, compression test, or indentation. Limited studies performed on using multi-axial loading measurements of time-dependent material behavior exist in the literature. Therefore, the aim of this study is to investigate the viscoelastic response of rubber-like materials under multi-axial loading using cavity expansion and relaxation tests. The tests were performed on PVA hydrogel specimens. Three hyperelasitc models and one term Prony series were used to characterize the viscoelastic response of the hydrogels. Finite element (FE) simulations were performed to verify the validity of the calibrated material coefficients by reproducing the experimental results. The excellent agreement between the experimental, analytical and numerical data proves the capability of the cavity expansion technique to measure the time-dependent behavior of viscoelastic materials.
Rubber-like materials are known for their viscoelastic behavior. Elaborate characterization of this behavior is crucial to numerous applications in a variety of engineering fields and industries, in addition to the biomedical field. The time-dependent properties of these materials control their mechanical response under thermal conditions [1,2], and large deformations [3,4]. Thus, these properties have been the focus of numerous studies that utilized different mechanical characterization techniques.
In the biomedical field, studying the viscoelastic properties of biological tissues indicated that their mechanical properties change based on their pathological conditions. For instance, fibrosis was reported to change the viscoelastic behavior of liver tissues [5,6]. This change occurs at early stages of fibrosis mainly due to degeneration of fatty tissues . The mechanical quantification of such diseases at early stages is very critical for treatment processes. Currently, the measurement of the mechanical properties of diseased and healthy tissues in-vivo is commonly achieved using transient elastography (TE) [8,9]. Although this technique is capable of measuring the mechanics of tissues, it is mainly used to quantify their Young’s modulus when small deformations are applied as demonstrated in [10,11]. On the other hand, TE proved to be efficient in measuring the viscous (relaxation) response of these tissues in-vivo .
In general, a number of theories have been developed to describe the time-dependent response of different materials. Those theories describe materials as linear viscoelastic (LVE) , quasi-linear viscoelastic (QVE) , and non-linear viscoelastic (NVE) [15,16]. The LVE description is the most common due to its simplicity in addition to the validity of the superposition principal it utilizes. Moreover, it has been thoroughly examined mathematically , and experimentally . It is also available in commercial finite element (FE) packages. Therefore, it is the focus in this study.
The implementation of the LVE theory is based on combining the instantaneous elastic response of a given material with its viscous behavior. This behavior is commonly measured by the creep test , or the relaxation test , including in-vivo characterization in cases of biological tissues . The elastic behavior is measured commonly by using axial loading and indentation [21–23]. These conventional techniques are used to measure the material response under different deformation rates. Although uniaxial loading is very common when studying the mechanics of viscoelastic materials, to our knowledge, little attention has been paid to viscoelasticity characterization using multi-axial loading. Thus, in this work, we aim to investigate the viscoelastic properties of a rubber-like material under different rates by using cavity-based multi-axial loading (pressure + equi-biaxial tension). This technique is based on inducing an expanding cavity within the structure of soft materials by using a balloon and an incompressible fluid. This technique was found to be efficient in measuring the hyperelastic properties of rubber-like materials [24,25]. Therefore, it represents a good option to investigate the rate-dependent behavior of such materials under multi-axial loading condition.
A hyperelastic material that exhibits the LVE behavior is PVA hydrogel as will be demonstrated in this study. These hydrogels are considered excellent alternatives to biological tissues due to their high hydrophobicity [26,27], biocompatibility , mechanical strength [29,30], physical integrity under large deformation [24,31], viscoelastic properties [32,33], thermal stability and non-toxicity . Thus, they are the most suitable material for this study. Cavity expansion tests and simple shear relaxation tests are conducted to measure the instantaneous and long-term behaviors of PVA hydrogels, respectively. The data obtained from these tests will be used to calibrate the material constants of hyperelastic models (Yeoh, Arruda-Boyce and Ogden), and the Prony series viscoelastic model coefficients. Finally, FE simulations are performed to reproduce the time dependent behavior of PVA hydrogels.
This study is based on using the cavity expansion and the relaxation tests to measure the instantaneous and time-dependent behaviors of the PVA hydrogels, respectively; thus, twelve cylindrical samples (40 mm in height and 36 mm in diameter) were prepared. The synthesis process was based on physically crosslinking PVA solution. The solution was made by mixing 99+% hydrolyzed PVA (molecular weight of 146000–186000 g.mol-1 with deionized water (10% w/w ratio) at 90 oC by using standard flask/column combination. Thereafter, the solution was poured into cylindrical molds. The hydrogels were formed by freezing the solution at -20 oC for 3 hours and then thawing at room temperature (~21 o C) for 3 hours. This freeze-thaw cycle (FTC) was repeated two times. Thereafter, three gel cylinders were cut into smaller cylinders (5 mm in height) then punched in the middle to create rectangular cuboid samples with dimensions of 20 mm in length, 10 mm in width and 5 mm in thickness. Fig 1 shows the samples’ geometry for both cavity expansion and shear tests.
The procedure followed in  was adopted to perform the cavity test. However, few changes were made to the system that was used: pressure sensor (Model PRESS-S-000, PENDOTECH, USA) and pressure reader (Model PMAT-S, PENDOTECH, USA) were used to observe the pressure due to the pressurized incompressible fluid (water). In addition, an Instron machine (Model 4465; Canton, MA, USA) was used to introduce the water into the needle-balloon tool. A new design needle-balloon tool was used in this study, the balloon is silicon based with smaller dimensions and has softer mechanical behavior than the one used in the previous study. To ensure no air was entrapped inside the system, its elements (syringe, needle-balloon tool and Y-shaped tubes) were filled with water and submerged under water until they sank and assembled together. Afterwards, the tubes were connected to the pressure sensor with a water inlet and outlet, which allows depleting any remaining air in the system. Nine hydrogel samples were investigated by this test at three different deformation rates: 5 μl/s, 20 μl/s and 50 μl/s. A schematic diagram for the system is shown in Fig 2. In addition, the balloon response was evaluated separately using the Instron machine.
This test was performed on the gel cuboids; they were subjected to simple shear relaxation test. Acrylic platens (50 mm × 20 mm × 2.5 mm) were machined with spacers glued to the platens to ensure a good alignment as well as to avoid overstressing of the gel specimens. The top and the bottom surfaces of the cuboid specimens were placed and glued between the platens with a thin layer of a fast-acting adhesive. The assembly was then mounted on the Instron machine and loaded at a rate of 0.1 s-1 . The relaxation test was performed at various shear strains: 10%, 20%, 30%, 40%, 50%, 60%, 70% and 80% to investigate the behavior of the hydrogels in a step-like strain at various magnitudes. Fig 3 shows the loading configuration for the simple shear relaxation test.
In this work, the time-dependency of the PVA hydrogel behavior is modelled as LVE. This description is considered accurate for the tested gels as will be shown in section 4. A general description of the LVE material stress in large deformations at time (t) can be expressed as reported in 
The first term represents the hyperelastic stress function (σh.el). The second term represents Prony series expression of stress relaxation function (σv.el).
At time t+Δt, and by using “N” number of terms in the Prony series, Eq 1 can be expressed as
σh.el (t+Δt) can be evaluated by using material parameters of a calibrated strain energy function (SEF), and, while (t+Δt) can be expressed as
Where gi and τi are the relaxation parameters. See  for full derivation.
A number of strain energy functions (SEFs) are available in literature. In this work, we focus only on three SEFs, namely: Yeoh, Ogden and Arruda-Boyce. These models are available in most of the commercial FE software packages; in addition, they proved to be efficient in capturing the non-linear elastic response of rubber-like materials subjected to uniaxial loading, biaxial loading , and cavitation . Yeoh, Ogden, and Arruda-Boyce models are expressed, respectively, as
Where μi, αi, Ci0 and Ci are material constants; λm is the extensibility limiter.
In the cavity expansion test, the pressure applied on the cavity wall generates radial stress normal to the wall and hoop circumferential stresses along the wall. These stresses can be expressed in Cauchy hyperelastic form as
Where W1 and W2 are the SEF derivatives with respect to I1 and I2 , respectively. Eq 8 can be used to calculate σh.el (t+Δt).
The cavity tests were performed on cylindrical specimens at each injection rate () of 5 μl/s, 20 μl/s and 50 μl/s up to a volume of 110 μl in order to investigate the behavior of the gels at particular loading rates, as shown in Fig 4. The pressure (MPa) and time (s) data were measured directly by the system. The volume data shown in Fig 4 were calculated through multiplying the time data by their corresponding injection rate (t ˑ). It was observed that, at the initial stage of the balloon inflation, the gels did not show any resistance for time periods of 1.96 s, 0.45 s and 0.16 s at injection rates of 5 μl/s, 20 μl/s and 50 μl/s, respectively. Multiplying these short periods of times by their corresponding volume rates will result in initial volumes (Vi ) of 8.635 μl, 8.23 μl and 8.24 μl, respectively. The lack of resistance to these volumes can be due to the displacement of the channel wall created by the needle during the insertion process, which relieves any constraints against balloon inflation at early stages. Afterwards, the gels start to resist deformation as the balloon inflates spherically. During the test, the hydrogels were internally deformed to a volume of 110 μl. At this volume, the balloon material started to participate in resisting the inflation. Thus, to avoid the effect of the balloon’s silicon material, all tests stopped at 110 μl. The balloon response was investigated separately at the lowest rate (5 μl/s) by the Instron machine, using the built-in linear variable differential transformer (LVDT) and load cell, see S1 Fig. The effect of the friction between the rubber-stopper and the syringe plunger was also investigated, see S2 Fig
The cuboid specimens were loaded at various shear strain levels (10% - 80%) and held for 300 seconds to record the relaxation force as stresses showed no tangible relaxation after 300 seconds of loading. The force data was then converted into stress form using the cuboid dimensions, as shown is Fig 5A. The long-term shear stresses and their correspondent strains were plotted, and it was found that the gels exhibited a linear long-term behavior, see Fig 5B, which indicates the validity of the LVE assumption. It is worth mentioning that stress-strain data in the simple shear test represented plane stress-plane strain behavior, i.e., the initial area (20 mm × 10 mm) that was used to calculate the shear test from the force data remains unchanged after loading the specimens.
The experimental data obtained from the cavity expansion tests were in pressure (P)–volume (V). The deformation term (V) was converted into a hoop stretch term (λ), which was calculated as . “R” values were found to be 1.399 mm, 1.382 mm and 1.383 mm at injection rates of 5 μl/s, 20 μl/s and 50 μl/s, respectively. These magnitudes were calculated using the concept of Equivalent Volume Diameter (EVD)  as follows:
Vi can be evaluated as
Where Vn-b is cylindrical volume (5 mm in length and 0.85 mm in diameter) of the balloon region before inflation; Vp are the volumes of water introduced into the balloon, 8.635 μl, 8.23 μl and 8.24 μl (calculated in the Results section), before the hydrogel started to resist balloon inflation at each of injection rates, 5 μl/s, 20 μl/s, 50 μl/s, respectively.
The magnitude of “r” was calculated as
Where Vapp is the injected volumes of water during the cavity expansion test.
Fig 6A shows the experimental data in the form of P–λ of the three injection rates. The instantaneous rate of change of the three curves was calculated. Before the commencement of the calibration process, it is worth mentioning that due the local nature of the cavity expansion test, Eq 9 can be expressed as
This indicates that the radial stress at the cavity wall, where the material parameters are calculated, is equivalent to the applied pressure.
The resulting theoretical curves exhibited strong agreement with experimental data. Arruda-Boyce model showed the least fit, yet it is still valid, (Fig 6B, 6C and 6D). The material parameters are summarized in Table 1.
The initial shear moduli of the PVA hydrogels can also be estimated from the calibrated material parameters at each deformation rate. These moduli are summarized in Table 2.
In this work, Abaqus® software package was used to evaluate the Prony series coefficient using its built-in calibration tool. The long-term shear stress data, shown in Fig 5, was converted into shear moduli data, and then normalized and input in Abaqus® to run the calibration process. One term Prony series provided an excellent fit to the relaxation response of the hydrogel, see Fig 7. The Prony series coefficients were g = 7.837e-02 and τ = 47.845 s.
Most of the reported work on the viscoelasticity of PVA hydrogels was describing the material as QLV [39,40]. Thus, it was necessary to validate the predictable behavior of the LVE description adopted in this study by using numerical simulations, an analytical solution, and measured experimental data.
In order to reproduce the experimental results, FE simulations were performed using Abaqus® by implementing the model described in , and the material parameters obtained in the previous section. A mass density of 1030 kg/m3 (±65 kg/m3) was used for the model’s elements (CAX4H). This magnitude was estimated by using Archimedes principal. The observed pressure during the experiments was applied in the FE simulations.
The FE simulations predicted numerical pressure (radial stresses at the cavity wall) that was in a good agreement with experimental pressure. In addition, the simulations provided numerical hoop stresses that were matching their analytical counterparts calculated using Eq 2, see Fig 8. This equation can be solved using Matlab® , see S1 Appendix for an example solution using Yeoh model for 50 μl/s data.
The numerical results showed an excellent agreement with the experimental data. In addition, they predicted numerical hoop stresses that are in robust agreement with the analytical hoop stresses. These two agreements at each of the deformation rates are good indications to the strong validity of the calibrated material parameters of the implemented SEFs.
The data obtained from the cavity tests were experimental P-λ and analytical σθ—εθ. The SEFs’ parameters, used to fit experimental data as well as to predict the σθ—εθ data analytically, were verified by using FE. The root mean square (RMS) error was calculated to evaluate the difference between the realistic data (experimental and analytical) and the reproduced data (numerical). The error values are summarized in Table 3
The characterization of PVA hydrogels using cavity expansion at different deformation rates allows for measuring the rate-dependent response of rubber-like materials in multi-axial loading. In this research, the properties of PVA hydrogels subjected to cavitation loads have been characterized up to 80% hoop logarithmic strain at deformation rates of 5 μl/s, 20 μl/s and 50 μl/s. Force relaxation experiments in simple shear were also performed at various strain magnitudes (10% - 80%).
The parameters of the three SEFs were calibrated through Eq 9, while the time-dependent Prony coefficients were calibrated using Abaqus’ calibration tool. Overall, the good agreement between the numerical, experimental and analytical results indicates that the three SEFs are suitable to model the viscoelastic behavior of the PVA hydrogels up to 80% strain.
The initial shear moduli of the PVA hydrogels were estimated by the SEFs to be ≈ 2.5 kPa, 5.5 kPa and 7 kPa at deformation rates of 5 μl/s, 20 μl/s and 50 μl/s, respectively. While the three models predicted similar initial shear moduli at each deformation rate, the overall performance of Arruda-Boyce model in data fitting was the least sufficient. This is due to its limited number of parameters unlike the other two phenomenological models. However, Arruda-Boyce model has the capacity to predict the material response in other forms of loading using the same material parameters due to its micromechanical nature (use of eight Langevin chains network) [41,42], an advantage that does not exist in other empirical models such as Yeoh and Ogden.
In this study, the PVA hydrogels were considered incompressible. This assumption was based on the very low compressibility nature of this type of gels. Nafo and Al-Mayah  reported that there is no significant difference between the compressible and incompressible assumptions in simulating the behavior of PVA hydrogels subjected to cavity loads. In addition, Chen et al  reported a relatively high Poisson’s ratio for this type of hydrogels (≈0.5); thus, we considered the gels as incompressible materials in this study. It is also worth mentioning that the volumetric behavior of most rubber-like materials has very little or no time-dependency . Therefore, only the relaxation parameters (gi and τ) were calibrated and considered for the viscous response of the hydrogels.
To investigate the syringe-balloon response as well as the injection system at various injection rates, we also collected data from the machine’s load cell. The friction force between the rubber-stopper and the inner wall of the plunger increased when the injection velocity was increased as shown in S2 Fig. Overall, the friction response is divided into two stages, static and dynamic. The dynamic stage seems very steady when the injection velocity was low. However, when the injection velocity was increased, the force data exhibited relatively volatile behaviors. The two friction phases observed by the load cell makes it difficult to estimate Vp, particularly, at high velocities, which may lead to fallible calculation of “λ” and inaccurate characterization of material properties. Therefore, observing the pressure from an independent source (pressurized water that flows into the system) will result in measurements that are more accurate. There are limited data, to authors’ knowledge, available on testing the viscoelastic response of PVA hydrogels under multi-axial loading condition; however, the test proved to be capable of measuring the rate-dependent response of the hydrogels. Attributed to its multi-axial nature and simplicity, the technique has the potential to provide enough information about in-vivo viscoelastic behavior of biological tissues in different loading axes, which will be addressed in a different study.
There are a number of models to describe linear viscoelasticity [45,46]. However, the most commonly used model is the generalized Maxwell model. This model is usually represented by a Prony series , which is used to fit the experimental data to evaluate the discrete times (τi) and the long term shear moduli (gi ). Although this representation is feasible and the fitting process is relatively simple, Prony series is prone to challenges related to the coefficients evaluation process, particularly, when the fitting process renders negative coefficients, which result in oscillations in the modelled material behavior . However, this challenge can be addressed using several fitting approaches that have been developed over the past few decades [49–51]. On the other hand, the use of Prony series is efficient analytically and numerically as shown in this study, if the coefficients were determined correctly. This is mainly because Prony series permits analytical integration of the relaxation equations.
Measuring the rate-dependent response of rubber-like materials commonly performed using uniaxial loading. This investigation introduced an alternative that allows for multi-axial measurement to the viscoelastic properties of these material. The following conclusions can be made from this study: