1 INTRODUCTION
Hyperelastic materials are widely used in engineering applications from seismic absorber devices (^{Chiba and Furukawa, 2010}) to aortic valve prosthesis or stents (^{Saab, 1999}). Therefore, much effort has gone into describing the behavior of materials like silicone, latex, rubber (^{Treloar, 1974}; ^{Arruda and Boyce, 1993}; ^{Meunier et al., 2008}), or biological tissues (^{García-Herrera et al., 2011} and ^{2012}); Eilaghi et al., 1984) with the aim of quantifying their ability to undergo very large deformations that can be recovered when the specimen is downloaded without loss of its original abilities (^{Ogden, 1984}).
Several analytical models have been proposed in the literature to describe the deformation energy function W (e.g., see ^{Ogden, 1984}). Ogden, Mooney-Rivlin, Yeoh, Neo-Hookean, or Demiray are the most extensively used models to describe hyperelastic materials. A relevant aspect for dealing with practical applications is to determine the related material constants that realistically describe a specific material's behavior (^{Ogden et al., 2004}). Moreover, to be able to evaluate real situations, numerical modeling is a standard analytical tool. Nowadays, numerical developments to properly deal with very large deformation and the appropriate mechanical characterization of the materials are crucial to get trustworthy predictions. To this end, the simulations also need to be contrasted with experiments. In addition, to analyze the pressurization of hyperelastic membranes, dynamic models have emerged especiallyin applications with rapid loading in thermoforming processes (^{Verron et al., 2001}).
Tensile, compression, shear or inflation tests (^{Treloar, 1974}) are reported in the literature to characterize mechanically the material's behavior, trying to have the best description with a limited number of constants taken from a unique experimental test (^{Arruda and Boyce, 1993}). As an example, silicone behavior is obtained from tensile, compression, shear and inflation test in ^{Meunier et al. (2008}) using image correlation techniques while ^{Skouras et al (2012}) reports on a characterization using the Hart-Smith model via the inflation test. Rubber-like materials are analyzed by ^{Gonzáles et al. (2009}) by fitting the Mooney-Rivlin model with experimental biaxial tests. Uniaxial and biaxial traction tests are presented by ^{Sasso et al. (2008}) to characterize polymers including a finite element (FEM) analysis. The Ogden model is applied in ^{Canseco et al. (2011}) to describe the behavior during compression under the ASTM D 3574 standard of a commercial material named Poron. Mechanical characterization of arteries has been reported by ^{García-Herrera et al. (2011} and ^{2012}), where tensile tests and inflation tests are described. The descending human aorta is studied by ^{García-Herrera et al. (2011}) using the Holzapfel model, while the Demiray model is used to characterize the behavior of a patient-specific human aortic arch in ^{García et al. (2012}).
The present work reports the experimental data obtained from tensile and inflation tests of latex balloons. The mechanical characterization of latex is performed using the tensile test by adjusting the experimental curves with the Yeoh model. Although other models have been used, the best fit was obtained using that model. The inflation test is afterwards taken as a benchmark to assess the performance of the modeling via two different approaches: analytical and numerical. The numerical simulation is performed using a finite element technique able to deal with large deformations that was previously reported and extensively used to analyze other engineering applications (see ^{Celentano et al., 2001}, ^{2009} and ^{2012}). The main objective of this work is threefold: to experimentally characterize the material, to validate the proposed modeling strategy by comparing the results with the experiments and, finally, to extend this analysis to a biomedical problem. The former objective is based on the assumption that the inflation test mimics the pressurization of human eyes. To this end, the study of an eyeball subjected to internal pressure like that reached under glaucoma conditions is presented. The obtained results are compared with those reported in the literature.
The rest of the work is organized as follows: the governing equations are presented in Section 2; the mechanical characterization is described in Section 3; Section 4 reports the inflation of a balloon test including the experimental procedure, measurements and modeling results; Section 5 presents the analysis of an eyeball under glaucoma conditions. Finally, the concluding remarks are summarized.
2 FUNDAMENTAL EQUATIONS
2.1 Governing Equations
In a general context of the mechanics of a continuous medium, the governing equations that model the situation are the continuity and equilibrium equations expressed as:
where ρ_{0} is the density of the material in the initial configuration, ρ is the density in the spatial configuration, and J = det F is the volume relation between the deformed and the initial configuration where F = ∂x/∂X is the deformation gradient tensor (x and X being the spatial and material coordinates, respectively) and u is the displacement vector. The Cauchy stress tensor σ can be expressed by means of a specific deformation energy function W as (^{Holzapfel, 2000}):
where C = F^{T}F is the right Cauchy-Green tensor. In this work, material incompressibility is assumed (i.e., J=1).
2.2 Constitutive Equations
Hyperelastic isotropic materials are characterized by the deformation energy density function expressed as a function of C, e.g., ^{Ogden (1984}) and ^{Gonzáles et al. (2009}), or the invariants of the deformation tensor, e.g., Neo-Hooke in ^{Holzapfel (2000}); Mooney-Rivlin or Yeoh models in ^{Rackl (2015}) and Demiray in ^{Holzapfel (2000}). The Yeoh model is given by:
where I _{1} is the trace of C, and C _{10}, C _{20} and C _{30} are constants of the material to be determined by some fitting method.
Principal Cauchy stress determination for uniaxial tensile samples
The expression of the principal Cauchy stress σ_{i} as a function of the constants of the models and the principal stretchings λ_{i} is:
The incompressibility and isotropic behavior for the case of simple uniaxial traction is given by the following relation:
where sub-indexes 1, 2 and 3 respectively refers to stretchings along the axial, width and thickness directions of the sample.
The principal Cauchy stress computed from Equation (5) including the incompressibility constraint (6) is written as:
Substituting Equation (4) into (7), the analytical expression for the uniaxial Cauchy stress can be determined:
Internal pressure determination for thin-walled hollow spheres
The stress in a thin-walled hollow sphere (i.e., T < R/10, where T is the wall thickness and R is the radius, both in the initial configuration) subjected to a constant inner pressure p_{i} can be determined considering the unique relationship between circumferential λ_{θ} and longitudinal λ _{z} stretching, i.e., λ_{θ} = λ _{z} , and the material's incompressibility. Thus:
The analytical expression for the internal pressure considering the deformation energy given by Equation (4) is:
3 MECHANICAL CHARACTERIZATION
3.1 Tensile Test
The simple tensile test allows setting a uniform stress state in order to be able to determine the mechanical behavior of the material. The effectiveness of this approach has already been verified by other authors in the field of nonbiological, see ^{Canseco et al. (2011}), as well as for biological hyperelastic materials, as in ^{García-Herrera et al. (2011}).
The samples tested in the present work, obtained as shown in Figure 1 from two different principal directions of the sample, have the following dimensions: length 6.3 mm and width 2.29 mm. The thickness in the balloon was found to be variable. In this work, a thickness profile was measured with a micrometer. The average values obtained along the balloon length around 10 mm apart starting from the pole were 0.184 ( 0.005, 0.143 ( 0.011, 0.132 ( 0.008 and 0.115 ( 0.009 mm.
The tensile tests were carried out on a universal testing machine Instron 3342 equipped with a 500 N loading cell at a controlled cross-head displacement speed of 5 mm/min and a constant testing temperature of 25 ºC; see Figure 2.
The curves obtained for the longitudinal and circumferential samples are presented in Figures 3 a) and b), respectively, where the horizontal axis describes the axial stretch λ and the vertical axis indicates the axial Cauchy stress σ (samples were tested for each direction). These plots are computed from the instantaneous axial load P - displacement u (measured between clamps) curves registered during the tensile tests. The Cauchy stress is σ=P/A, A being the instantaneous transversal area of the specimen analytically obtained under the material incompressibility assumptionas A=Ω_{0}⁄L_{f} , where Ω_{0}=L_{0}A_{0} is the initial volume of the specimen, with L_{0} and A_{0} respectively being the initial length and initial transversal area, and the instantaneous specimen length is calculated as L_{f} =L _{0}+u. Considering that λ=L_{f} /L _{0}, the stress-stretch relationship is σ=P/A_{0} λ. The results presented in Figure 3 show that in some samples the curve has discontinuous stretching zones (approximately at λ = 8.0) because during the test some samples undergo partial breakage. In spite of the above, it is seen that for both directions of the tested samples all the tests present the same nonlinear behavioral trend which is characteristic of hyperelastic materials.
3.2 Determination of Material Parameters
The difference between the curves depicted in Figures 3 a) and 3 b) for the longitudinal and circumferential samples is less than 10%, thereby the response of the tested material could be considered as isotropic, i.e., the mechanical behavior does not depend on the orientation of the material. The material constants are obtained by fitting by least squares the average value of the experimental curves reported in Figures 3 with the analytical stress defined by Equation (8). The expression used is:
where σ _{exp} and λ _{exp} are the experimental axial stress and elongation values, respectively, and n is the number of experimental measurements. From Equation (12), we get the material constants by minimizing G. Table 1 summarizes the resulting parameters.
To verify the behavior of the model with the obtained constants, the analytical and measured uniaxial stresses are compared in Figure 4 where it is seen that the Yeoh model presents a satisfactory fit with the experimental behavior for the studied material.
4 INFLATION TEST
4.1 Experimental Design
The pressurization test is performed by injecting air into a latex balloon. The test circuit is illustrated in Figure 5. The objective of this test is to determine the state of the deformations present in the balloon wall during its inflation. The injected mass is at a temperature of 25 ± 0.3 °C controlled with a Proportional Integral and Derivative (PID) device. The pressure is measured by the UNIK 5000 commercial sensor (by Gemeasurement), whose range of operation is from 0 to 0.5 bar with an error of ± 0.2%. The pressure data are captured with a USB acquisition 1608FS card (by MC Measurement Computing).
The procedure for running the test is: first, the tank (1) is pressurized with compressed air at 0.3 bar (inlet pipe, IP), the mass flow of air to be injected into the balloon is adjusted by means of the graduated valve (2), the air is injected into the balloon located in the test chamber (3) by opening a control valve (4) in the outlet pipe (OP), and the pressure inside the balloon is recorded by the pressure sensor (5) installed in the same outlet pipe.
The experiments are registered using two IDS Imaging Development System digital video cameras UI-1545LE-M-GL (using a 75 mm Pentax lens) at their maximum resolution of 1280x1024 pixels and video recorder velocity of 10 fps. Two cameras perpendicular to one another were used to check that the 3D effects were negligible in this test. Later, the video is processed to obtain the balloon profile evolution in terms of the applied pressure. To this end, a code based on the library OpenCV (^{Bradski, 2000}) was developed to fit an ellipse to the contour measured.
4.2 Experimental Results
Figure 6 shows a sequence of pictures as the balloon is inflated which allows us to obtain the evolution of its profile. The dot indicates the center of the ellipse fitted with which the displacement evolutions at the polar and equatorial zones were measured. These displacements were measured with respect to the reference configuration adopted here as that with an inner pressure of 8.3 mbar (see Figure 6) in order to prevent geometric instabilities caused by the self weight of the balloon (the major and minor semi axis in this configuration are a=15.66 mm and b=6.95 mm). As the inflation process proceeds, the balloon tends to a spherical shape.
The pressurization test gives curves of the inner pressure p_{i} versus the displacements at the polar d_{pole} and equatorial d_{equa} zones plotted in Figure 7, where 8 samples were tested (the bars denote the standard deviation). As already commented, these curves are measured with respect to the adopted reference configuration. The experimental curves exhibit a behavior similar to that reported by ^{Verron et al. (2001}) and ^{Verron and Marckmann (2003}), where it is seen that after reaching a maximum pressure, it decreases as the displacements continue increasing. It should be noted that non-uniform strains are developed in the balloon during both the load and unload stages. For instance, for the maximum pressure average value of 70 mbar, the polar and equatorial displacements are 4.27 and 3.66 mm, respectively. The unstable character of the problem is reflected in the very different time intervals associated with the load and unload stages respectively observed in the experiments as 130 s and 2 s.
On the other hand, the experiments have shown that during the inflation of the balloon the material is concentrated in its polar zone, see Figure 8. ^{Voorhies (2003}) reported that such a zone presents the lower risk to develop high stress when the internal pressure increases.
4.3 Modeling the Experiments
As mentioned above, a reference configuration with a small pressure level was chosen to analyze the whole inflation process (see Figure 6). To this end, the estimations of both the initial stress and strain fields together with the thickness profile associated with this configuration must be performed. These estimations are computed in two steps. First, the approach to obtaining a compatible initial stress and strain fields is iteratively tackled by solving the equilibrium equations in the reference configuration (with an internal pressure of 8.3 mbar as indicated in Figure 6) together with the characterized Yeoh model presented above using the finite element method until the condition of a nearly zero displacement field for the whole balloon is fulfilled (^{García-Herrera and Celentano, 2013}). Second, from this simulation, the thickness profile in the reference configuration is computed by completely removing the pressure to check that the thickness distribution at the end of the deflating process corresponds to the measured values reported in Section 3.1. The resulting thickness values applying this procedure obtained along the balloon length around 10 mm apart starting from the pole are: 0.183, 0.136, 0.125 and 0.114 mm. As expected, note that these values are slightly smaller than those listed in Section 3.1.
The study of the inflation process was carried out via two different approaches, analytical and numerical, separately described below.
The analytical results were obtained from the thin-walled hollow sphere case (using Equation (11)) considering a simplified geometry given by an average reference radius of 11.31 mm (assumed as R=(a+b)/2 with a and b being, as mentioned above, the major and minor axis of the reference configuration) and a uniform thickness of 0.140 mm (obtained as the average of the thickness profiles in the reference configuration). The analytical curve shown in Figure 7 qualitatively reproduces the material response during both the load and unload stages. The results in the load stage closely approach the measurements. However, from the maximum pressure onwards, the pressure values are clearly and progressively underestimated.
The numerical analysis is performed in order to more properly take into account the real geometry of the balloon (i.e., ellipsoidal shape with a variable thickness distribution). In the present study, an axisymmetric uniform mesh composed of 12048 four-node elements is considered for the reference configuration (with a height of 45.6 mm), where the thickness is discretized with 8 elements. The element size was determined by a standard mesh convergence analysis (not shown) performed until accurate results with a reduced computation time are obtained. The displacements are fully constrained at the bottom of the balloon while the horizontal displacement is restricted at the top to reproduce symmetry conditions. Figure 9 illustrates the numerical strategy adopted in the present simulation to impose the boundary conditions. First, as shown in Figure 9a, an internal pressure is incrementally imposed on the inner wall of the balloon up to a maximum value of 70 mbar according to the experimental records (note that the air hydrostatic pressure is negligible in this case). Then, in order to circumvent the instability appearing at the instant of maximum pressure that in turn causes numerical convergence problems due to the loss of convexity of the energy function W, the displacement field measured in the experiments is imposed on the curved domain while the straight part of the balloon is loaded with a decreasing pressure according to the experimental values registered. Note that the imposed non-uniform displacement profile (which, as observed in the experiments, has both radial and axial components) is variable along the inflation process. Figures 9b to 9d schematically represent the different displacement boundary conditions during the unload stage. Using this procedure, the numerical internal pressure - displacement curves considering the characterized Yeoh model are depicted in Figure 7. The simulation satisfactorily matches the experiments during the load stage both in pressure and polar and equatorial displacements. In particular, the maximum pressure value is adequately captured. For displacements larger than those corresponding to the maximum pressure, the pressure, obtained as the normal component of the stress tensor at the inner surface of the balloon, exhibits a decreasing trend that agrees with that of the measurements. The predicted slopes of the unload stage are, however, steeper than those of the experimental curves. The excessive deformation that takes place afterwards (particularly in the equatorial zone as shown below) is presumably the cause that hinders the achievement of a converged numerical solution for larger levels of imposed displacement. Although the present numerical analysis is not able to describe the full inflating sequence, it allows, in contrast to the analytical approach, the prediction of non-uniform stress and thickness profiles during the deformation process. These results are shown below.
The computed circumferential stress in terms of equatorial displacement and inflation pressure are respectively plotted in Figures 10a and 10b. The non-uniform stress distribution along the deformation process can clearly be seen. The maximum stresses reached are lower in zones near the top and bottom of the specimen. The equatorial zone (zone 3) presents the greatest circumferential stress for a given equatorial displacement. The stress localization in this zone is apparent when comparing its final stress level at d_{equa}=5 mm (around 1100 kPa) with that corresponding to the thin-walled hollow sphere case, i.e., 784 kPa computed using Equations (10) and (11). Note that although the pressure decreases at the end of the test, the stresses always increase thus denoting the geometrical character of the instability.
Figure 11 shows the thickness ratio - internal pressure curves computed in different balloon zones. As expected, the thickness decreases continuously. The thickness reduction is practically linear before the maximum pressure is reached. The section near the balloon pole (zone 5) is thicker than the others, showing the concentration of mass in that zone. After the maximum pressure is reached, the thickness in the pole remains practically constant, with a nearly uniform value at the end of the test. This result agrees with the qualitative observation presented in Figure 8. In addition, the equatorial zone (zone 3) exhibits the larger thickness reduction. This fact is consistent with the stress localization commented above. Note that the final thickness ratio in this zone (0.38) is clearly lower than that corresponding to the thin-walled hollow sphere case, i.e., 0.48.
5 APPLICATION TO A BIOMECHANICAL ANALYSIS OF HUMAN EYE TISSUES
Once the model and the analytical methodology proposed for the inflation test study had been developed and verified experimentally, the problem of the eyeball subjected to inner pressure was analyzed. Due to the harmful effect of the high intraocular pressure that occurs in the disease called glaucoma (^{Ruiz - Ederra et al., 2005}; ^{Norman et al., 2010}), different experimental (e.g., ^{Schultz et al., 2008}; ^{Elsheikh et al., 2010}; ^{Elsheikh et al., 2008}) as well as numerical (e.g., ^{Belleza et al., 2000}; ^{Lari et al., 2012}; ^{Asejczyk et al., 2011}) studies have been proposed to determine the mechanical characteristics of the sclera (^{Eilaghi et al., 2010}; ^{Elsheikh et al., 2010}; ^{Schultz et al., 2008}) and the cornea (e.g., ^{Elsheikh et al., 2008}) and, in addition, to establish the state of the stress on the walls of the eyeballs subjected to an internal pressure (^{Belleza et al., 2000}).
In the present work, the constants of the Yeoh model are determined from the experimental results of the simple stress tests reported by ^{Elsheikh et al. (2010}) for the sclera and by ^{Elsheikh et al. (2008}) for the cornea, both of them on healthy human tissues. The corresponding curves fitted by least squares to the Yeoh model are presented in Figures 12 and 13. The mechanical characterization of the material is used in the numerical simulation considering the eyeball geometry reported by ^{Stitzel et al. (2002}), which is reproduced in Figure 14. The model is 2D with radial symmetry with respect to the vertical axis, and the uniform finite element grid is composed of 1895 four-node elements. The contour conditions correspond to fixed in the base of the ball (sclera, gray color), free vertical displacement at its cusp (cornea, blue color), and an imposed internal pressure of up to 200 mmHg. The parameters of the materials for the Yeoh model used are given in Table 2.
Model | Parameters | ||
---|---|---|---|
C _{10} | C _{20} | C _{30} | |
Yeoh Sclera | 0.1217 | 0.0058 | 0.0476 |
Yeoh Cornea | 0.0084 | 0.0149 | 0.0112 |
Figure 15 shows the plot of the computed pressure - strain curve, with the strain D_{f} - D _{0}/D _{0}, where D _{0} and D_{f} are the initial and instantaneous diameters of the eyeball at the sclera near the optic nerve zone, respectively. As a qualitative confirmation of the predictions, the experimental data reported by ^{Lari et al. (2012}) during the inflation test of pig eyes is also plotted.
The stress-strain computed curves are shown in Figure 16a together with the results reported by ^{Lari et al. (2012}) for comparison in Figure 16b. Similar trends can be appreciated where in Figure 16a it is seen that the sclera presents a lower degree of deformation (λ_{θ} = 1.34) with respect to the cornea (λ_{θ} = 1.9) under the same stress load, so it is inferred that the sclera shows a more rigid behavior than the cornea.
The circumferential stress-stretch and the circumferential stress-pressure curves for the sclera at zones 1,2,3,4 and 5 are plotted in Figures 17a and 17b, respectively. Similar results are plotted in Figures 18a and 18b for the cornea. From Figure 17, different stress distributions can be seen in different zones of the sclera. This is attributable to the different initial wall thickness, which increases from 0.52 mm (in zone 1) to 1.6 mm (near the optic nerve). The maximum stress value is obtained at the equatorial zone, reaching approximately 0.66 MPa. Figure 18 reports that the stresses in different zones in the cornea (see Figure 14) are very similar to one another because no thickness changes take place.
The thickness evolution for the cornea is presented in Figure 19, showing similar trends between different zones varying from the initial value of 0.52 mm to 0.23 mm. From these results the equatorial zone has a high risk of rupture, as it was also seen in the experimental work reported by ^{Voorhies (2003}), where from the pressurization test results, 80% of eye rupture has occurred in the equatorial plane.
Finally, Figure 20 depicts the thickness evolution with pressure. The thickness decreases continuously, with the rate of reduction greater at the beginning, and this was also reported by ^{Elsheikh et al. (2010}). At high pressure levels the thickness decreases slowly, reaching a constant value of 0.52 mm near the optic nerve (zone 5 in Figure 14). This value is closer to that reported by ^{Bisplinghoff et al. (2009}) for a pressurization test of human sclera up to rupture, registering a final thickness of 0.58 ┴ 0.13 mm.
6 CONCLUSIONS
An experimental and numerical study of a latex balloon pressurization test is presented. To this end, a material characterization was firstly conducted using the hyperelastic Yeoh model. Then, the modeling of this problem was tackled via two different approaches: analytical and numerical. The analytical solution considering a simplified geometry was found to reasonably describe the pressure - displacement curves. On the other hand, the real geometry of the balloon (i.e., ellipsoidal shape with a variable thickness distribution) was taken into account by means of a finite element model proposed to analyze this problem. Although the present numerical analysis has shown a limited performance due to the instability that occurs at the instant of maximum pressure that in turn precludes the analysis of the full inflation sequence, it satisfactorily predicted, in contrast to the analytical approach, the following relevant aspects of the problem: the wall stresses are not uniform, the geometry does not evolve spherically, the equatorial zone presents higher deformations, and the polar zone exhibits material concentration during the pressurization process.
This experiment served as a reference for analyzing the pressurization test in human eyes. The developed methodology was applied to predict the mechanical behavior of a human eyeball subjected to internal pressure, mimicking glaucoma diseases. The study of the eyeball was performed using properties and geometry taken from the literature. The best fit of the material behavior was obtained using the Yeoh model in both materials analyzed. The results obtained in the present work for such a test agree well with the experimental data reported in the references, confirming that agreater stress is developed in the equatorial region where sclera rupture occurs.