Numerical Simulation and Experimental Validation of the Inflation Test of Latex Balloons

Experiments and modeling aimed at assessing the mechanical response of latex balloons in the inflation test are presented. To this end, the hyperelastic Yeoh material model is firstly characterized via tensile test and, then, used to numerically simulate via finite elements the stress-strain evolution during the inflation test. The numerical pressure-displacement curves are validated with those obtained experimentally. Moreover, this analysis is extended to a biomedical problem of an eyeball under glaucoma conditions.

Latin American Journal of Solids and Structures 13 (2016) 2657-2678 Arruda and Boyce, 1993;Meunier et al., 2008), or biological tissues (García-Herrera et al., 2011 and2012); 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 (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 and2012), 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., 2001Celentano et al., , 2009Celentano et al., 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 Latin American Journal of Solids and Structures 13 (2016) 2657-2678 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.

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 ρ is the density of the material in the initial configuration, ρis the density in the spatial configuration, and det is the volume relation between the deformed and the initial configuration where is the deformation gradient tensor (x and X being the spatial and material coordinates, respectively) and is the displacement vector.The Cauchy stress tensor can be expressed by means of a specific deformation energy function as (Holzapfel, 2000): where is the right Cauchy-Green tensor.In this work, material incompressibility is assumed (i.e., J=1).

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 I1 is the trace of C, and C10, C20 and C30 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 σ as a function of the constants of the models and the principal stretchings λ is: Latin American Journal of Solids and Structures 13 (2016) 2657-2678 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: (8)

Internal pressure determination for thin-walled hollow spheres
The stress in a thin-walled hollow sphere (i.e., , where is the wall thickness and is the radius, both in the initial configuration) subjected to a constant inner pressure can be determined considering the unique relationship between circumferential and longitudinal stretching, i.e., , and the material's incompressibility.Thus: (9)

2
(10) The analytical expression for the internal pressure considering the deformation energy given by Equation ( 4) is: 3 MECHANICAL CHARACTERIZATION

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).
Latin American Journal of Solids and Structures 13 (2016) 2657-2678 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.   and Structures 13 (2016) 2657-2678 Lf =L0+u.Considering that λ=Lf/L0, the stress-stretch relationship is σ=P/A0λ.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.

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 σ and 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 .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.

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.

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.66mm 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 versus the displacements at the polar dpole and equatorial dequa 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.

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 dequa=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 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 addi-Latin American Journal of Solids and Structures 13 (2016) 2657-2678 tion, 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.

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 byElsheikh 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.    Figure 15 shows the plot of the computed pressure -strain curve, with the strain , where D and D 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 in- creases 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.

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.

Figure 3 :
Figure 3: Uniaxial tensile test results for a) longitudinal and b) circumferential specimens.

Figure 6 :
Figure 6: Balloon inflation sequence during the pressurization test.

Figure 7 :
Figure 7: Inflation pressure -displacement curves in the pressurization test.

Figure 8 :
Figure 8: Material concentration zone in the pressurization test at 50 mbar.

Figure 9 :
Figure 9: Boundary conditions for the FEM simulation.

Figure 10 :
Figure 10: a) Circumferential stress -equatorial displacement and b) circumferential stress -inflation pressure computed curves in different zones of the balloon.

Figure 11 :
Figure 11: Thickness ratio -inflation pressure computed curve in different zones of the balloon.

Figure 12 :
Figure 12: Experimental results of data fitting with the Yeoh model for the human sclera.

Figure 13 :
Figure 13: Experimental results of data fitting with the Yeoh model for the human cornea.

Figure 14 :
Figure 14: 2D geometry of the human eye for FEM simulation (the boxed numbers indicates the zones of interest in this study).
Figure 16: a) Stress -stretch curve for human sclera and cornea (average values from the present computation) and b) stress -strain curve for pig sclera and cornea taken from Lari et al. (2012).

Figure 18 :Figure 19 :
Figure 18: a) Stress -circumferential stretch and b) stress -pressure computed curves in different zones of the cornea.

Figure 20 :
Figure 20: Thickness -pressure curves in different zones of the sclera.

Table 1 :
Material parameters obtained by curve fitting from experimental data.

Table 2 :
Material parameters obtained by curve fitting of experimental data.