Slight

Finite Element simulations of rubbers and biological soft tissue usually assume that the material being deformed is slightly compress-ible. It is shown here that in shearing deformations the corresponding normal stress distribution can exhibit extreme sensitivity to changes in Poisson’s ratio. These changes can even lead to a reversal of the usual Poynting eﬀect. Therefore the usual practice of arbitrarily choosing a value of Poisson’s ratio when numerically modelling rubbers and soft tissue will, almost certainly, lead to a signiﬁcant diﬀerence between the simulated and actual normal stresses in a sheared block because of the diﬀerence between the assumed and actual value of Poisson’s ratio. The worrying conclusion is that simulations based on arbitrarily specifying Poisson’s ratio close to 1 / 2 cannot accurately predict the normal stress distribution even for the simplest of shearing deformations. It is shown analytically that this sensitivity is due to the small volume changes which inevitably accompany all deformations of rubber-like materials. To minimise these eﬀects, great care should be exercised to accurately determine Poisson’s ratio before simulations begin.


Introduction
When solid rubber-like materials are assumed to be incompressible so that only isochoric (i.e. volume preserving) deformations are allowed, then the theory of non-linear elasticity is known to be extremely effective in providing analytical predictions for their mechanical response when both the original geometry and boundary conditions are simple. These analytical solutions include shear, torsion, bending, straightening and inflation (see, e.g., Rivlin [2], Green and Zerna [3], or Ogden [4]). To analyze more complex problems, use of Finite Element (FE) simulations is usually necessary. However, in most commercial finite element codes, incompressibility is not assumed ab initio in order to prevent element locking (see, e.g., [5,6,7]). Thus in FE analyses, slight compressibility is usually assumed and the solution for fully incompressible materials is obtained in a limit process if necessary (see, e.g., [7,8,9,10]). Indeed, because all solid rubbers are to some extent compressible, models for slight compressibility have been investigated independently of any FE considerations (see, e.g., Ogden [4], Horgan and Murphy [13,14,15,16] and references cited therein).
Rubber components are often subjected to shearing deformations in applications and biological soft tissue is often sheared in vivo (see, for example, LeGrice et al. [17] for a discussion on shear strain in the left ventricular myocardium or Horgan and Murphy [18] on the simple shearing of tissues). It is therefore essential that accurate, reliable and efficient models be available to predict shear behaviour. In a recent paper, Gent et al. [1] conducted some numerical simulations using the commercial code ABAQUS of the simplest possible shearing deformation, where one face of a block, modelled as a compressible neo-Hookean solid, is displaced relative to the parallel face. This shearing deformation will be called here experimental simple shear. One set of their results is puzzling: in their Figures 12 and 13, an infinitesimal change in Poisson's ratio (ν) leads to a large change in the predicted values of the normal stress components, with, in some cases, a small change in Poisson's ratio leading to a change in the sign in the components, leading effectively to an 'inverse' Poynting effect. These numerical results have been confirmed by our own numerical simulations ( §2) for which the adopted procedure mirrors, as closely as possible, the earlier work.
In §3 this extreme sensitivity of the normal stresses to changes in ν is explained by combining a realistic mathematical model of experimental simple shear that incorporates the possibility of small volume changes with the constitutive models of slight compressibility used in FE simulations. It will be shown that the predicted normal stresses for the slightly compressible neo-Hookean material obtained from this combination display exactly the sensitivity to changes in the Poisson ratio as was first observed by Gent et al. [1]. The shearing strains considered here are only moderate in range, with a maximum shear strain of 100% being imposed. It is likely that the extreme sensitivity noted here is more pronounced for larger strains.
We suggest that the results discussed above have significant implications for the numerical simulation of materials traditionally considered as being incompressible. We note that in practice, when simulating these materials, experimental values of the corresponding Poisson ratio are not known. Instead, most engineers adopt the pragmatic solution of specifying instead a value of Poisson's ratio 'close' to 1/2, assuming implicitly that the specific value of Poisson's ratio adopted will not have a significant impact on the numerical results and, in particular, that small percentage changes in the value of Poisson's ratio result in corresponding small percentage changes in the stress distribution. Most simulations of biological soft tissue, for example, are conducted on this basis. However, the simulations of Gent et al. [1] and our own numerical experiments ( §2) all suggest that, with the models currently used, a small change in the value of Poisson's ratio chosen (say from 0.499 to 0.495) will result in significant changes in the stress distribution. There is almost certainly a difference between the assumed and actual values of Poisson's ratio and this should generally result in a significant difference between predicted and actual stress distributions. Our conclusion is that an accurate determination of Poisson's ratio is essential if FE simulations are to yield accurate normal stresses in shearing deformations (See also Horgan and Murphy [19] for a discussion of the results of Gent et al. [1] relevant to the present paper.) It is finally noted that the seemingly obvious solution to this extreme sensitivity to Poisson's ratio of simply simulating perfect incompressibility is not a valid approach. No material is perfectly incompressible and assuming this idealisation (denoted in our results by ν = 1/2) will result in a difference between the idealised value assumed and the actual value of Poisson's ratio. This difference will lead to a significant percentage error in the predicted normal stresses.

Modelling slight compressibility
Many different constitutive models have been proposed to reflect slight deviations from incompressibility on assuming that the material is homogeneous, isotropic and hyperelastic. We refer to the recent papers [13,14,15,16] for background, references to the pertinent literature and, very importantly, a summary of the fit of these models, and, in particular, the models used here, with experimental data. Here we briefly describe the usual form of such slightly compressible (or almost incompressible) strain-energy functions used in the FE simulations.
Assume that an incompressible strain-energy function W = ψ (λ 1 , λ 2 , λ 3 ), where ψ is a symmetric function of the principal stretches λ i , has been obtained that gives an excellent fit to the experimental data collected from some set of material characterization experiments. Then the corresponding slightly compressible form implemented in FE models is usually of the form where F is some specified polynomial function of its indicated argument, and thē λ's are the so-called 'deviatoric stretches' (see the ABAQUS [11] and ANSYS [12] manuals). On assuming that ψ (1, 1, 1) = 0 to ensure zero strain-energy in the reference configuration for the incompressible material, it follows from (1) that we need to assume that for the same property to hold for the slightly compressible material. To ensure zero stress in the reference configuration, it is sufficient to require that Finally, to ensure compatibility with the classical linear form of the elastic strain-energy function on restriction to infinitesimal deformations, it is required that where µ and κ are the infinitesimal shear and bulk moduli, respectively. Recall that they are connected to Poisson's ratio ν through Motivated by these relations, it is now assumed that the ψ term in (1) is proportional to µ and that F is proportional to κ. On setting i 3 = 1 in (1), one recovers the original strain-energy function for the corresponding incompressible material. It is important to emphasise that the main motivation for strain-energy functions of the form (1) seems to be mathematical convenience: a simple additive term in i 3 is included to account for the compressibility so that, on restriction to infinitesimal deformations, ψ and F are proportional to the shear modulus µ and the bulk modulus κ, respectively (see e.g., Horgan and Murphy [16] for further details).
Many FE codes are based on models of the form (1). The initial moduli are usually assumed to be such that κ/µ 1 in order to model slight incompressibility. In the variational formulation, the bulk modulus κ is called the penalty number (see, e.g., [5,20,7]) and is chosen to be large. From (6), we find that for ν = 0.49, 0.499, 0.4999, we have κ/µ 50, 500, 5000, respectively. The default setting in Abaqus/Explicit is quite low, as it assumes that κ/µ = 20, corresponding to a Poisson ratio of ν 0.475. Equivalently, the model can be written in terms of the Cauchy-Green strain invariants, defined in terms of the principal stretches λ i , i = 1, 2, 3, as follows, The corresponding strain-energy function that is used in FE simulations to model almost incompressible behavior then has the equivalent form 3 . (8) Further details can be found, for example, in Crisfield [10].
The general representation of the stress-strain relations for compressible, non-linear elasticity can be found, for example, in Atkin and Fox [21] and is given by where Here T is the Cauchy stress, I is the identity tensor, and B is the left Cauchy-Green strain tensor, defined in terms of the deformation gradient tensor F through the relation B = F F T . The stress-strain relation for the almost incompressible material, (8), therefore has the form (9) with the following coefficients: The stress-strain relation for almost incompressible materials can thus be put in the form with α 1 , α 2 given in (11) Thus the stress is always the sum of two parts: a hydrostatic term proportional to the volumetric function F (i 3 ), and therefore of the order of the infinitesimal bulk modulus κ, and another term of the order of the infinitesimal shear modulus µ. The numerical simulations of Gent et al. [1] are based on the slightly compressible neo-Hookean specification of this strain-energy function, readily implemented in ABAQUS [11]. It has the form where the constitutive volumetric function F is the quadratic form which is the simplest form to satisfy (3), (4), and (5) 2 . For this material the constitutive coefficients simplify tô Predictions based on the neo-Hookean model have been found to yield a good fit with data from a variety of experiments on rubber-like materials for small to moderate strains. Since moderate strains are only considered here, the results obtained in our analyses are therefore likely to be obtained from simulations based on other forms of the deviatoric part of the strain-energy function.
The form (14) was assumed in our repeating, and furthering, of the numerical simulations of Gent et al. [1]. We entered their physical and geometrical parameters into ABAQUS 6.7 to try to replicate their results for a sheared neo-Hoohean block. We obtained the stresses and strains throughout the thickness. We varied the mesh resolution, the time step duration and the number of increments used and found virtually no difference with their simulations. Figure 2 shows an example of the stress field distribution in the sheared block, displaying clearly an almost homogeneous field except near the edges. Further pictures obtained for the stress and strain distributions are identical to those of Gent et al. [1], see Figure 1.
We also carried out two re-meshes in addition to the original mesh [1] of 20 (thickness) by 100 (length), i.e. 2,000 rectangular elements. The first re-mesh was a square mesh of 20 × 200 (i.e. 4,000 elements) and the second re-mesh was a square mesh of 40 × 400 (i.e. 16,000 elements). Each new mesh was run for the two following values of Poisson's ratio: ν = 0.49 and ν = 0.5. These re-meshes had no effect on the original results of Gent et al. [1]. Similarly, we varied the time step and incrementation, from their original [1] values of time step = 4 seconds, 400 increments (giving a shearing rate of 2 mm/s). Thus we halved (time=2, increments=200, shearing rate=4mm/s) and then doubled (time=8, increments=800, shearing rate=1mm/s) the time steps and increments. These time changes also has no effect.
In conclusion, we found, just like Gent et al. [1], that for the compressible neo-Hookean model (14)- (15) infinitesimal changes in the value of Poisson's ratio can lead to finite changes in the simulated stress distribution. Our simulations suggest that this problem is especially acute for ν ∈ [0.495, 0.5] and that that the problem of extreme sensitivity of the normal stresses to Poisson's ratio increases the closer the value of ν is to 1/2. These features are evident in Figure 1 and, perhaps more emphatically, in a comparison of Figures 2 and 3. In Figure 2, where perfect incompressibility has been assumed, the stress is almost homogeneous, except for more complicated behaviour near the corners, and tensile throughout a block with aspect ratio 1:10 sheared by a unity amount. Figure 3 displays the normal stress component T 22 distribution with Poisson's ratio ν = 0.495. In contrast with Figure  2, the stress is no longer homogenous: it can be seen there that T 22 is tensile and non-uniform near the slanted faces, but is compressive towards the centre of the specimen.
It is natural to see if the problems in simulating slightly compressible materials are reflected in analytical approaches to simple shear because their specific cause might then be isolated. To investigate these problems analytically, the classical mathematical approximation to experimental simple shear is generalised in the next section and the corresponding stress distribution obtained for both the general slightly compressible model (12) and for the slightly compressible neo-Hookean model (17) used in FE simulations.

Analytical approaches to simple shear
The classical mathematical model of simple shear assumes that a block of rubber is deformed according to where (X 1 , X 2 , X 3 ) and (x 1 , x 2 , x 3 ) denote the Cartesian coordinates of a typical particle before and after deformation, respectively, and γ is an arbitrary dimensionless constant called the amount of shear. Within the context of nonlinear elasticity, this problem was first considered by Rivlin [2] and has been much discussed and analyzed since. The reference works of Green and Zerna [3], Truesdell and Noll [22], and Ogden [4], for example, all consider various aspects of this problem.
Intuitively it seems reasonable that this homogeneous deformation should accurately model experimental simple shear, especially within the bulk of the material and for modest amounts of shear. The corresponding in-plane stress distribution for the general strain-energy function adopted in FE simulations to model almost incompressible rubbers is easily obtained from (11), (12) and is given by Modelling experimental simple shear using the intuitive deformation (18) therefore results in the in-plane stress field for almost incompressible materials of the form (1) being independent of the volumetric function F . This is not surprising given that classical simple shear is an isochoric deformation (i 3 = 1). Therefore, for a fixed value of µ, the in-plane stress field is independent of the bulk modulus κ, or equivalently, of the Poisson ratio ν. This is in conflict with the numerical predictions discussed earlier.
This conflict can be explained by the fact that using the isochoric approximation (18) to model experimental simple shear is too restrictive and does not capture an essential feature of every deformation of rubbers, i.e., that an infinitesimal volume change always occurs and therefore occurs in experimental simple shear in particular. To get an idea of the order of the change, we rely on the classical experimental data of Penn [23]. In a series of elegant experiments on peroxide vulcanizates of natural gum rubber using a dilatometer technique, Penn measured the volume change in simple tension for axial stretches of the order of 2. The volume change, calculated as i 3 − 1, was found to be of the order of 0.0001. No comparable data exist for simple shear. We now investigate the effect, if any, of a change in volume.
We model the infinitesimal volume change with the following simple generalization of the classical approach (18), where the 's are infinitesimal quantities. This deformation seems intuitively to capture the main features of the deformation of a block when one face is displaced relative to a parallel face. It corresponds to a simple shear of amount γ superposed on a triaxial extension, where the 's are the infinitesimal elongations.
Retaining only the first-order terms in the 's gives i 3 1 + , where ≡ 1 + 2 + 3 is the volume change. It follows that, to the first-order in , the constitutive function F (i 3 ) occurring in the stress-strain relation (12) for almost incompressible materials can be approximated as where we used (4) and (5) 2 . The corresponding in-plane stresses are given by Rajagopal and Wineman [24] as where λ i ≡ 1 + i are the stretch ratios (here it is understood that their powers are linearized with respect to the 's). Compare these general expressions with those for simple shear given in (19): the shear stress T 12 is again independent of the bulk modulus and volume change, but now both normal stresses T 11 and T 22 contain a new κ term. This is in agreement with the effects observed in the numerical simulations, showing that the normal stresses, not the shear stress, are sensitive to volume changes (see Figures 1  and 4(a)). Now we try to mirror the numerical simulations: we specialize the stressstrain relation to that of the slightly compressible neo-Hookean solid (17); we impose a parallel displacement for the moving face with respect to the fixed face (so that λ 2 = 1); and we impose a two-dimensional formulation (so that λ 3 = 1). From [24] we find that the first principal invariant is I 1 = 3+γ 2 +2 , where we linearized with respect to = 1 . For this specialisation the normal stress component T 22 , for example, is easily computed and we find that where we used the connection (6) between the bulk modulus, the shear modulus and Poisson's ratio. To compare this stress with the simulated stress given in Figure 2, the volume change = i 3 − 1 is needed. In Figure 4(b) i 3 , as measured by ABAQUS, is plotted as a function of the amount of shear γ for different values of Poisson's ratio: ν = 0.499, 0.495, 0.49. In Table 1, we collected in particular results obtained when the amount of shear is γ = 1.0, as in the simulations. The volume change i 3 − 1 is displayed in the second column and the normalized Cauchy stress component T 22 /µ, found from (23), in the third column. It follows that (23) predicts normal stresses that are very close to the average values of T 22 predicted by the numerical simulations given in Figure  1(b). Therefore infinitesimal changes in Poisson's ratio close to the limiting   Table 1: Variation of volume change with Poisson's ratio value of ν = 1/2 may indeed result in dramatic changes in the value for the normal stress, even leading to an inverse Poynting effect, due to the volume changes which accompany every deformation.

Conclusion
It has been shown that the usual practice of arbitrarily choosing a value of Poisson's ratio when numerically modeling rubbers results in widely different normal stress distributions when the third decimal of ν ∈ [0.495, 0.5] is modified. The sign of the Poynting effect may even be reversed. It is thus essential for the accurate modeling of rubbers to experimentally determine Poisson's ratio to the utmost precision. This is not as hard as it seems, especially when the measure of ν relies on ultrasonic wave speeds evaluation. There, the contrast between the speed of a compression wave v L and the speed of a shear wave v T is so large for nearly incompressible solids, that the connection gives excellent precision. For instance, Wood and Martin [25] find ν = 0.49986 using acoustic waves. Similarly, Gennisson et al. [26] measure the following body wave speeds for an Agar-Agar gel: v L = 1500 m/s and v T = 1 m/s, giving ν = 0.49999978, according to (24). The issues raised here ought to be addressed by anyone wishing to use a commercial FE code for large deformation behavior. For instance the Abaqus/Explicit manual [11] indicates that the default setting for slight compressible models is ν = 0.475, a value which allows volume changes sufficient to inverse the Poynting effect according to Figure 1. It goes on to state that if we "are defining the compressibility rather than accepting the default value, an upper limit of 100 is suggested for the ratio of κ/µ", which corre-sponds to ν 0.495, also within the range of extreme sensitivity in shear. Furthermore, the manual goes on to say that slight compressibility "does not warrant special attention for plane stress", a statement clearly at odds with our findings.