In the previous studies by the authors and others, it was demonstrated that there are two possible defect growth modes and a characteristic material length for any soft material. For a pre-existing defect smaller than the material characteristic length, the energy is dissipated all around the defect as it grows and the critical load for the growth is independent of the defect size. For defects larger than the characteristic length, the growth is by cracking and the energy is dissipated along a plane. Thus, the critical load for the growth is size dependent and can be predicted by fracture mechanics. In this study, we apply the same energy-based argument to the failure of thin membranes, with the focus on the first growth mode that gives the maximum critical load. We assume that strain localization due to damage is the precursor to rupture, and hence, we model the corresponding zone as a through-thickness hole, with its size smaller than the material characteristic length. The defect grows when the elastic energy relaxed by the growth is enough to provide the energy needed for internal microstructure changes. This leads us to the size-independent failure conditions for membranes under the biaxial load. The conditions are expressed in terms of either two principal stretches or two principal stresses for two different types of materials. For verification, we test the theory using the published experimental data on natural and styrene-butadiene rubber. By using the experimental data from equal biaxial loading, we predict the critical principal stretch ratios and critical stresses for different biaxialities. The predictions agree well with the experimental results.
Throughout the literature, it has been postulated that the common source of failure in soft elastic materials is the growth of microscopic defects [1–13]. The defects can be inhomogeneities such as microcavities, microcracks, and weakened or highly strained zones. Therefore, studying the defect growth criteria has been of great interest. However, the main goal of such studies is to predict the failure, which is mainly signaled through a critical state of stress or deformation.
The first approach in studying the mechanical failure gives various classical strength theories [14–16]. These theories usually do not make assumptions on microstructure of the materials and give the size-independent critical load.
However, for soft materials, there is still no widely accepted strength theory. For a neo-Hookean material as an early model for rubbers, an ultimate hydrostatic tensile load of pu = 5E/6 was proposed by Gent and Lindley . Their ultilmate load was obtained by taking the limit of the elastic solution as the cavity elastically expands to infinity. However, the theory did not address the material failure itself, and such elastic limit load does not exist for many other polymers [11,17–20].
Another challenge in strength theory approaches is that although the critical load can be measured for a specific loading type, such as uniaxial loading, it is not easy to be extended to the general multi-axial loading cases. This is especially the case for highly nonlinear elastomers. Some other examples of strength theories are based on ultimate stretch ratios or ultimate stretch tensor invariants such as I1. However, they are not so satisfactory for general multi-axial loadings . The idea that cavitation or failure is not only determined by deformation has been transcribed in the literature and thoroughly reviewed recently by Poulain et al. .
The second type of approach in studying failure leads to the development of fracture mechanics, which gives the size-dependent critical stress. Fracture mechanics introduces an important material property: fracture toughness. Given a pre-cut large crack with a certain size, fracture toughness can be measured experimentally. The critical load for crack propagation is then determined by matching the energy-release rate or J-integral to the fracture toughness [11,22–33]. However, applying the fracture criterion to a real material needs the initial crack size. That is usually not available in practical applications, especially in design processes and in small scale structures.
Continuum damage theory or phase-field models [34–45] might be one way of resolution in defect initiation. For example, energy limiter , local ultimate state variables , and a damage variable  are chosen as failure measures. However, in some cases, the related parameters are hard or impossible to be calibrated by experiments because the evolution in this incubation or incipient stage may not be shown explicitly or externally.
In a recent study by the authors [47,48], we proposed two defect growth modes, omnidirectional growth, and cracking. In the first mode, energy is dissipated all around the defect, and in the second mode, energy is only dissipated in a layer next to the cracking plane. The theory suggests that there is a characteristic length for any soft material. If the defect is smaller than the characteristic length, the defect grows by omnidirectional expansion, and the critical load is size-insensible, as in the strength theories. If the defect is greater than the characteristic length, the defect grows by cracking and size-dependent fracture mechanics theories are applied. The argument is mainly from the geometric consideration and thermodynamics. The idea of the two-stage growth was also discussed in a series of papers recently, including a physics-based theory in forms of material chain scission  and progressive damage models [44,50]. The models used the phase-field–based approach, so that the two growth modes were inherently embedded in: for the small flaw size, the growth is volumetric, and for the large flaw size, the growth is by cracking.
When the defect is very small, we proved in Ref.  that the omnidirectional growth mode is predominant. As the defect grows to a characteristic size of A*, scaled with the ratio Γs/Γp, the growth mode transits from the first mode to the cracking mode. For A > A*, the defect grows by cracking, which is well discussed in fracture mechanics.
The flaw insensibility in small scales was first studied by Gao et al.  as a lesson from nature, where they demonstrated that biological materials select certain structural features in a nanoscale to optimize their strength. The concept of flaw size dependence was revisited by Chen et al.  when studying fatigue. They proposed somewhat same characteristic length when discussing the notch sensitivity. They estimated the typical values of A* to be of ∼1 mm for natural rubber (NR) and ∼10 mm for tough hydrogels . Mao et al.  then gave a clear physical picture using an elastomer chain model and showed that the fracture stretch loses the dependency on the small notch size. Volokh  also estimated a length scale for natural rubber using the work done to failure in their energy limiter theory and the fracture toughness estimated by Rivlin and Thomes . The critical load is insensitive to the defect size if it is smaller than the transitional defect size.
In the above discussion, we assume that the energy barrier for the defect growth is mainly due to the energy dissipation around the defect. Recently, Poulain et al.  and Kumar et al.  postulated and experimentally showed that cavitation is not purely dissipative, i.e., defects in deca-microns can be healed and disappear on unloading, leaving an intact material zone. Thus, besides dissipative energy, the energy reduction in healing process is also a part of the energy barrier. One part of this healing energy is surface energy. As in metal or ceramics, it is proportional to the increment of surface area. In soft materials, however, the kinetic processes involved in healing, such as recombinations, recoils and mass transport, all occur in the zone around the defect with a finite volume. Thus, we define Γp as the total energy needed to drive the unit volume of mass out of its original location next to the cavity surface. It includes the dissipative energy and the volumetric part of the healing energy.
It is novel to see that a size-independent failure criterion similar to classic strength theories can be derived from the Griffith  energy criterion, which is inherently size-dependent. This paper applies the proposed energy-based idea to the failure of thin membrane and derives the explicit failure conditions in terms of remote stresses or stretches. Nait-Abdelaziz et al.  studied the membrane failure under general biaxial loading by assuming an initial defect size and applying fracture mechanics criterion. The fundamental differences between the current study and their approach is that we do not need a specific defect size as they estimated certain sizes. We adopt a different defect growth criterion, and the results are independent of defect size as long as it is smaller than the characteristic size.
The failure of thin soft membrane has been of a great interest in recent years. For example, the natural cavitation and failure in thin membranes of biological tissues and hydrogels [56–60] and dielectric elastomers [61–65] fit into this study. For many soft membranes in practical application, we believe that there is no initial crack or defect that is larger than the characteristic material length. Thus, the critical load for the first growth mode might be more practically relevant than that of the cracking mode, and a size-independent failure criterion might be more useful in the design process.
The objective of this paper is to derive practical and directly usable failure criteria in terms of principal stretches or principal stresses for a membrane under general biaxial loading. Two classes of materials will be examined, one is called I1-materials, whose strain energy density functions mainly depend on the first invariant of Right Cauchy-Green deformation tensor, and the second is called I1-I2-materials, whose energy density functions depend on both invariants, I1 and I2.
In Sec. 2, the first growth mode will be considered for the failure of thin membranes. In the first part of the section, the configurational driving force for the defect growth will be calculated as a function of applied load for equi-biaxial loading. We propose to use the driving force in equi-biaxial case to match the experimental results to get a key parameter in the model, the volume process energy Γp. In the second part of the section, we calculate the driving forces for the general biaxial-loading cases. By using the obtained Γp from equi-biaxial and the driving force of general biaxial loading, we are able to obtain the failure conditions in terms of remote loading, either two principal stretches or two principal stresses. To verify the reliability of the proposed approach, we compare the predictions with the experimental results obtained by Hamdi et al. .
Thin Membranes Under Biaxial Stretching
Consider a thin membrane under general in-plane biaxial loading. We assume that rupture of membrane is first triggered by strain localization, internal microstructural change, and energy dissipation at certain zone, as shown schematically by the shaded area in Fig. 2(a). Physically, it could be due to preexisting defects such as surface defects, inhomogeneities, or the so-called damage zone as described in damage mechanics or phase-field models. As the external load increases, the zone evolves and grows irreversibly. However, it is surrounded by the elastic zone outside. Any bond breaking, chain untangling and sliding in the damaged zone cause strain relaxation in the surrounding elastic zone. Mathematically, for the elastic field of the surrounding zone, which has a relatively large size, the damage zone can be replaced by or modeled as a through-thickness hole of certain size A, as shown in Fig. 2(b). The model is valid as long as the displacements along the dashed lines are same as shown in Figs. 2(a) and 2(b). We assume that the hole is smaller than the material characteristic length A*, in consistence with the first growth mode, which is usually the case of failure initiation in soft materials. Since the critical load for the first growth mode is size independent, the defect size is not needed to calculate the critical load.
Assume that the thin membrane depicted in Fig. 2 is made of an incompressible isotropic hyperelastic material with initial thickness T and side length 2B. It is under remote biaxial Cauchy stresses p1 and p2. The two corresponding principal stretches are λ and λn. n is the ratio of two principal strains and represents the load biaxiality; for example, n = 1 is for equi-biaxial loading and n = 0 for pure shear (second direction constrained). In the thickness direction, the third principal stress is p3 ≈ 0, and the third principal stretch is λ−(n+1) due to incompressibility.
We assume that the hole is circular with the undeformed radius A and in the middle of the sheet. The circular shape is chosen so that only one geometric parameter is needed to describe the configuration of the defect, and also we intend to remove any directional preference other than the external loading itself. Moreover, if there is a real through-thickness and irregular-shaped defect in the membrane, owing to its high stretchability, the stress localization around the defect would soon smooth out as the remote load increases. In fact, the field slightly away from the defect is insensitive to the original shape of the defect. Thus, as a first-order approximation, a strain localized damage zone or physically real through-thickness defect of any irregular shape could be treated as a circular hole.
Equi-biaxial loading is easier to carry out experimentally, so we choose it to determine the volume processing energy Γp from the critical stretch ratio. For mathematical convenience, now we assume the external boundary is a circle of radius B so that we could use polar coordinates. For this axisymmetric problem, we can get the stress and deformation fields by solving an ordinary differential equation (ODE). The equi-biaxial solution also will be used as a benchmark to validate the finite element solution in Sec. 2.2.
Equation (4) is a nonlinear second-order ODE about r(R) for R ∈ [A B]. Instead of matching the stresses or stretches at the remote, we assume the stretch λθ at R = A as a given number λa and integrate and calculate the corresponding remote stresses and stretches. By assigning λa to λθ at R = A, we have the first boundary condition r(A) = λaA. From the traction-free condition at the hole surface and stretch-stress relation, we have , which gives out the second boundary condition . By using the Runge–Kutta method of sixth order , and the values of r(A) and r′(A), we solve Eq. (4) for several different hyperelastic materials.
Figure 3 plots the remote equi-biaxial stress p normalized by C versus λa for a number of material models. Unlike the spherical cavity under hydrostatic load, for this set of membrane materials, the remote stress has no elastic limit when the stretch at the hole approaches to infinity.
As shown in Eq. (7), the normalized driving force for an infinite membrane depends on the remote stretch λ or normalized remote stress p/C. Figure 4 shows the normalized driving force versus the normalized remote tension p/C for several hyperelastic material models. All curves follow the same path for small load and then diverge when p is roughly 2C. The normalized driving forces for I1-materials whose energy functions just depend on I1 (neo-Hookean, Gent  and Seitz et al. ) are very close, unlike the I1-I2-materials, whose energy functions depend on both I1 and I2 (Mooney  and Gent and Thomas ). For the former group, the normalized driving force linearly depends on the remote tension when p/C > 2, with the slope of about π/2.
Consider all the I1-materials in Fig. 4, although these materials give quite different near-field solutions around the hole, the closeness of the normalized driving force indicates that the relaxed strain energy in growth is mainly from the intermediate and far field, in which the normalized fields are similar at moderate or small stretch ratios. The similar result occurs in fracture mechanics, in which the details of the plastic deformation near a crack tip have a weak influence on the energy-release rate of a crack in an infinite body. However, the normalized driving force is quite different for a material whose energy function depends on I2 as well, partly because it needs two material constants to fit its material behavior even at moderate stress level.
The above results from the ODE can be used as a benchmark to cross check the results from the finite element for nonlinear large deformation. In Figs. 3 and 4, the finite element results, for the material model by Seitz et al. , are also plotted (as marks). Results from ODE and FEM are close, indicating that at least for the equal-biaxial loading case, our numerical treatment has no convergence issue.
To show the effect of the shape of the defect, driving force is calculated for a sharp crack of size 2Ac using FEM. To compare with the previous results, we define and normalize the driving force in the same way as in Eq. (1) using the size of an imaginary circular hole, A = 2Ac/π, that gives the same surface area as the crack [47,48]. The material is neo-Hookean. Figure 5 shows the normalized driving force for the crack calculated directly from FEM results (dashed line) together with that of a circular hole of the equivalent size (solid line). For comparison, normalized J-integral of the crack from FEM and the estimation by Yeoh  are also plotted. We expect that the driving force for the defect of other shapes is between the solid line and those dashed lines. For these curves, the normalized driving force linearly depends on the remote tension except when the load is small. Yeoh  estimated the slope to be π/2, which is also shown in the graph and fairly represents the driving force–load relation.
General Biaxial Load.
For the cases of general biaxial loading, we use the finite element method to solve the stress and deformation fields. We follow the same procedure as the equi-biaxial case to obtain the driving force, as described in Appendix.
First, we consider the neo-Hookean material. Figure 6 depicts the normalized driving force versus the normalized remote stress p for the neo-Hookean material under different biaxialities n, where p = p1 is a leading principal stress, the larger one of the two principle stresses. As expected, the equi-biaxial load, n = 1, yields the highest driving force. For loadings with biaxialities between n = −0.5 (uniaxial) and n = 0 (pure shear), the driving force values are close. As shown in the figure, to have the same driving force, loads with smaller biaxiality require a larger leading principal stress p.
For all the other I1-materials in Fig. 4, the driving force–remote stress relation and their dependence on biaxility is similar to the neo-Hookean materials. However, for I1-I2-materials, the relation between driving force and remote load and its dependency on biaxiality are completely different from I1-materials. For example, Fig. 8 shows the driving force versus remote load under different biaxilities for Gum Stock, which fits the Mooney's material model  with C1 = 0.1205 MPa and C2 = 0.0460 MPa. The linear relation between the driving force and the remote stress found in I1-materials does not fit for this material. Also, unlike in the neo-Hookean material where the driving forces for small biaxialities (between −0.5 and 0.5) are close, and for this material, the driving force is more sensitive to biaxiality when it is smaller.
Energy-Based Strength Theory for Membrane
Neo-Hookean and I1-Materials.
We also calculate the failure surfaces for other I1-materials. The results are the same as what depicted in Fig. 9 (two such examples are illustrated in Figs. 13 and 14). The failure surface is close to a straight line in the principal stress plane and a circle in the principal stretch plane. The failure surface in the principal stretch plane is close to . This is very close to the critical I1 criterion proposed by Gent , because is negligible at the failure. However, it is worth noting that the critical I1 criterion is not applicable to general 3D multi-axial loading cases, in which the hydrostatic stress plays a significant role in the rupture but has no effect on the macroscopic deformation.
Unlike I1-material, the relation between the driving force and the remote stress for Mooney's material (Fig. 8) is not close to linear. The failure surface is shown in Fig. 10 for values of Γp/C between 180 and 300, where C = C1 + C2 for Mooney's material. The failure surface in the principal stretch plane (Fig. 10(a)) is completely different from a circle and cannot be fitted by the critical I1 criterion. It also shows that the critical normalized stress or stretch at small values of n are relatively large. This is due to the dependence of energy density function w on I2, which is approximately λ4 for the equi-biaxial case and about 2λ for the uniaxial case when λ >> 1, where λ is the larger one of the two principal stretches. The big drop in energy density from equi-biaxial to uniaxial causes the large raise in the critical load from equi-biaxial to uniaxial, as shown in Fig. 10. Another important observation is that in the principal stress plane in Fig. 10(b), the failure surface again approximately forms a straight line in the principal stress plane.
Natural Rubber and Styrene–Butadiene Rubber.
Hamdi et al.  conducted a series of biaxial experiments on thin membranes made of two elastomers, NR and styrene–butadiene rubber (SBR). They proposed the Yeoh's model  for NR and the Ogden's model  for SBR. The fitting parameters of the materials are listed in Table 2. In the experiments, they measured the critical stretch ratios of the membranes under different biaxialities n: n = 1, 0.88, 0.82, 0.81 for NR, n = 1, 0.88, 0.82, 0.68 for SBR, and n = −0.5 (uniaxial stress state) for both. The measured values of critical stresses by Hamdi et al.  were scattered with a relatively large deviation, indicating the sensitivity to the samples or the conditions of testing.
The rupture in Hamdi et al.'s  experiments might start with the first stage of defect growth, which is triggered by a small-scale defect, i.e., highly strained zone. At this stage, the drastic change might not be observed because the process is localized around the defect and does not induce significant change in the remote stress or stretch. However, when the defect reaches the characteristic size, the transition to cracking is unstable and the membrane ruptures.
By using Yeoh's model for NR and Ogden's for SBR, we calculate the driving force versus remote stretch for the same biaxilities as in the experiments. Figures 11 and 12 show the normalized driving force versus the leading principal stretch λ, the larger one of the two principal stretches, for NR and SBR, respectively. In these two graphs, the measured critical stretches and the corresponding driving forces are marked on the curves for different n values. Different from our ideal model, the critical driving force from the experiments varies with biaxiality for both NR and SBR. This could be due to the statistical variation in the tests or the idealistic assumptions we make in the model. For example, energy dissipation is not exactly uniform around the defect, especially when biaxiality n is small. Also, a slight deviation of the measured stretch ratio in experiment could give large difference in the normalized driving force. However, this should not be a big concern for us, because Γp here is just used as an intermediate value. When we write the failure condition in terms of the principal stretches, as long as Γp is above a certain value, its variation only causes small change in the critical stretch. However, if we express the failure condition in terms of stress, the change of Γp causes the same order of change in pcr/C. This is what was revealed in Hamdi et al.'s experiments that the measured stresses at failure had very large deviation. We suggest to pick the lower measured value as the limit load for practical applications.
It should be pointed out that our numerical results (diamond) and the fitting curve only use one data point from the equi-biaxial test, but obtain the critical conditions close to the experimental results for other biaxialities. In Figs. 13(b) and 14(b), we also note the large deviation of the experimental data in critical stresses. Since the biaxilities in the biaxial tests in the study by Hamdi et al. are between 0.81 and 1 for NR and 0.68 and 1 for SBR, it is hard to tell how good the prediction is for other small biaxialities. One important check point is the extreme case of uniaxial loading, at which the model here is supposed to work worst due to the assumption of even energy dissipation around the hole. Yet, for both materials, we capture the lower experimental values in uniaxial loading tests by using the equi-biaxial data alone. This underestimation is acceptable from the safe design point of view.
We proposed two defect growth modes in Ref. , one is omnidirectional dissipation mode and the other is cracking. There is a characteristic material length that separate the two modes. Defect grows in the first mode when its size is smaller than the characteristic size and grows by cracking when its size is larger. In the both modes, we use the Griffith type energy criterion to determine the conditions for the defect growth. It results in the size-independent critical loading condition for the first growth mode as in the traditional strength theory. For the second growth mode, it follows in the scope of fracture mechanics. The results show that strength theory and fracture theory can be united under one theoretical framework. This paper uses this fundamental principle to study the rupture of membrane under biaxial loading.
This result is also close to the I1 criterion because is negligible at failure. However, we should be aware that the I1 criterion is not applicable for a three-dimensional solid under general multi-axial loading. We also calculate the failure conditions for I1-I2-materials whose energy density depends on both I1 and I12. The results are completely different from the I1 materials. The failure surface cannot be fitted by a circle, and the I1 criterion is not applicable for these materials. The critical stress in uniaxial loading is one order higher than that in equi-biaxial loading.
We used the experimental data for natural rubber and styrene-butadiene rubber to check our theoretical predictions. Only by using the critical stretch or stress in the equi-biaxial loading case, we can predict the critical stretches or stresses for other biaxial loadings with different biaxialities. The predictions are in good agreement with the experimental data on NR and SBR.
This work was supported by the National Natural Science Foundation of China (Nos. 11525210, 11621062, and 91748209), the Fundamental Research Funds for the Central Universities, and PSC CUNY Research Foundation (TRADA-45-324).
Appendix: Normalized Driving Force G¯(λ)
The result for general biaxial loading is similar; the first term of the right-hand side of Eq. (A11) turns into , and the second term is still the normalized work done by the corresponding traction in pulling the hole from the traction-free state to the uniformly stretched state, which is also the energy difference between the two states.