Benchmark Problems of Hyper-Elasticity Analysis in Evaluation of FEM

The paper proposes benchmark problems on exact solutions of hyper-elastic analysis, which can be used to evaluate analysis capabilities of rubber-like materials provided by a finite element program or other approximate solution methods. Special attention was concentrated on analysis and derivation of the exact solutions for the thick-walled rubber cylinders under internal pressure and axial extension, the thick-walled rubber balloons under internal pressure and the rubber cylinders under torsion or tension-torsion. Deformation and stress analysis on the above three cases were conducted to provide equations and methods for data processing. Exact standard solutions of the problems combined with the strain energy function of generalized high-order polynomials are given. Numerical examples and evaluation results of two commercial packages that are in common use (ABAQUS and ANSYS) are presented. Good agreements are found in the comparisons between the present exact standard solutions and the simulation results.


Introduction
Hyper-elasticity is a highly non-linear material property. Structural analysis of hyper-elastic materials is always combined with large deformation. Hence analytical solutions in this field were rarely found and numerical approximate methods, such as FEM (finite element method), become musts in this field. However validation of the approximate formulation and program implementation has to be done using some benchmark problems before practical application of an approximate method [1,2]. Many commercial FEM codes, such as reference [3,4], have provided extensive analysis options and incorporated many complicated demonstration examples regarding hyper-elasticity. However, users may like to understand the reliabilities of the analysis capabilities and correctness of the answers by typical benchmark problems with exact solutions.
In addition, the key issue of an appropriate strain energy function for a specific material is still uncertain. Many research works have been devoted in evaluation of the existing ones [5][6][7] or creation of new ones [8][9][10][11]. Thirdly, new methods of solution are still under development [12][13][14] and need to be implemented in commercial codes and tested by a benchmark problem before application.
Therefore, this paper tries to propose some typical problems of which exact numerical solutions can be obtained and can serve as tools to evaluate availabilities of hyper-elasticity analysis in FEM packages or other approximate solutions. Some of the simplest deformation modes, such as uniaxial extension, pure shear and biaxial extension, etc., which were traditionally used to do material property calibration, can no longer be used as standard solutions to check numerical methods. This paper tries Materials 2020, 13 to present some other problems of which solutions can be found by elementary analysis or given in literature. The strain energy functions of a generalized high-order polynomial will then be combined to give exact solutions and numerical examples. Finally the evaluation results of FEM packages ABAQUS and ANSYS will be reported.

Constitutive Relation of Hyper-Elastic Materials
The theory of constitutive relation of hyper-elastic materials has been developed extensively and a great amount of literatures can be found. Reference [15,16] indicated that ABAQUS and ANSYS, etc. have implemented many forms of the well-known strain energy functions, or potentials in short, in their material constitutive relation libraries for users to choose from. This paper has no attempt to introduce the achievements in this field but pays attention to conduct numerical calculations of some typical problems. Therefore the stress-strain relation from Rivlin [17] is cited and re-written into a form convenient for later derivations and computer programming. and where t ∆ -True stress; ∆ = 1→2→3→1-cyclical index indicating the principle directions of strain and stress; Λ = Principle stretch; W = Strain energy density function; W ,1 , W ,2 = Derivatives of strain energy density function to invariants of right Cauchy-Green strain tensor; t = Hydrostatic pressure.
There are many different forms of the strain energy density function. As examples, only the high-order polynomial functions will be considered hereafter: where C is a vector containing the material constants: while F is a vector containing the base functions of the polynomial: In which the invariants of the right Cauchy-Green strain tensor are: Taking derivation of Equation (3) with respect to the strain invariants result in with  (1) and (2) gives: and

Thick-Walled Rubber Cylinder
A thick-walled rubber cylinder under internal pressure and axial extension is shown in Figure 1. The dimensions and loads are denoted also. The symbol P and p denotes axial pulling force and internal pressure. L, R 0 , R i and T represent the length, outer radius, inner radius and thickness of the cylinder; R is an arbitrary radius; D is the distance from an arbitrary point to the inner wall. A fundamental assumption is going to be made: The well-known plane section assumption. Therefore, the deformation mode is quite simple: Uniform axial extension and uniform radial expansion. Besides, under the consideration of fully incompressibility, the deformation is completely defined by only two parameters: Axial extension of the cylinder length and radial expansion of the inner or outer radius.

Thick-Walled Rubber Cylinder
A thick-walled rubber cylinder under internal pressure and axial extension is shown in Figure  1. The dimensions and loads are denoted also. The symbol P and p denotes axial pulling force and internal pressure. L, R0, Ri and T represent the length, outer radius, inner radius and thickness of the cylinder; R is an arbitrary radius; D is the distance from an arbitrary point to the inner wall. A fundamental assumption is going to be made: The well-known plane section assumption. Therefore, the deformation mode is quite simple: Uniform axial extension and uniform radial expansion. Besides, under the consideration of fully incompressibility, the deformation is completely defined by only two parameters: Axial extension of the cylinder length and radial expansion of the inner or outer radius.

Deformation Analysis
Apparently, the deformation in the radial direction is non-uniform. Examine the deformation at an arbitrary radius R. Let ξ be the dimensionless coordinate: Let ri, r0 be the inner and outer radius after deformation. Due to full incompressibility, there exists: where l denotes length of the cylinder after deformation. Hence, or Figure 1. A thick-walled cylinder.

Deformation Analysis
Apparently, the deformation in the radial direction is non-uniform. Examine the deformation at an arbitrary radius R. Let ξ be the dimensionless coordinate: or Let r i , r 0 be the inner and outer radius after deformation. Due to full incompressibility, there exists: where l denotes length of the cylinder after deformation. Hence, An arbitrary radius deforms into: where t is the thickness of the cylinder after deformation.

Strain Analysis
The principle stretch, numbered 1, in axial direction is: The principle stretch in hoop direction, numbered 2, is: It is worth noting that the hoop strain varies with ξ. According to full incompressibility, the principle stretch in the radial direction, numbered 3, is: (20) which is also a function of ξ. The invariants are:

Stress Under Pre-Specified Deformation
The stress can be calculated if l and r i , are known. According to Equations (18)- (21), all of the strain quantities can be obtained. In addition, the stress calculation can start from the radial equilibrium equation. It is well known that in the small deformation regime the equation is [18]: where σ denotes nominal stress in hoop and radial directions. The symbol R denotes radius before deformation. In hyper-elasticity analysis, large deformation must be taken into account. The equilibrium relation is then expressed by true stresses and radius after deformation, as shown in Figure 2, where r denotes an arbitrary radius after deformation, t 2 and t 3 represent the hoop stress and the radial stress respectively. The equation becomes: By substitution of Equation (10) with ∆ = 2, ∆+1 = 3 into Equation (23), we have: Materials 2020, 13, 885

of 16
Taking integration in the whole thickness to calculate the internal pressure: Taking integration in the thickness direction to calculate the radial stress: where By substitution of Equation (27) into Equation (10) with ∆ = 2, ∆ + 1 = 3, we have: By substitution of Equation (27) into Equation (10) with ∆ = 3, ∆ + 1 = 3, we have: Finally the axial reaction force is computed by axial equilibrium: Taking integration in the thickness direction to calculate the radial stress: By substitution of Equation (27) into Equation (10) with Δ = 2, Δ + 1 = 3, we have: By substitution of Equation (27) into Equation (10) with Δ = 3, Δ + 1 = 3, we have: Finally the axial reaction force is computed by axial equilibrium:

Numerical Example
In order to verify the correctness of the analytical solution, a finite element model of thick-walled rubber cylinder was established. The cylinder, parameters of which are shown in Table 1, deformed uniformly under internal pressure and axial extension.

Numerical Example
In order to verify the correctness of the analytical solution, a finite element model of thick-walled rubber cylinder was established. The cylinder, parameters of which are shown in Table 1, deformed uniformly under internal pressure and axial extension. The problem is computed by Equations (27), (29)-(31). Very fine intervals are employed when conducting numerical integrations in the radial direction. Figure 3 shows the stress distribution in the radial direction.
Materials 2020, 13, x FOR PEER REVIEW 6 of 17 The problem is computed by Equations (27), (29)-(31). Very fine intervals are employed when conducting numerical integrations in the radial direction. Figure 3 shows the stress distribution in the radial direction. The radial expansion and axial extension are measured in percentage of the thickness and length, respectively. These curves provide a graphics approach to solve the rubber cylinder problem under action of forces rather than pre-specified displacements. To this end the curves with finer intervals of the axial extension can be drown following the same approach. As an example, a loading case with axial force P = 100 kN and internal pressure p = 0.03 MPa is considered. The first step is to draw lines in Figures 6 and 7 to get all possible deformation data under the load P or p, as shown in Figures 6 and 7. The obtained data are listed in Table 2.  The radial expansion and axial extension are measured in percentage of the thickness and length, respectively. These curves provide a graphics approach to solve the rubber cylinder problem under action of forces rather than pre-specified displacements. To this end the curves with finer intervals of the axial extension can be drown following the same approach.
Materials 2020, 13, x FOR PEER REVIEW 6 of 17 The problem is computed by Equations (27), (29)-(31). Very fine intervals are employed when conducting numerical integrations in the radial direction. Figure 3 shows the stress distribution in the radial direction. The radial expansion and axial extension are measured in percentage of the thickness and length, respectively. These curves provide a graphics approach to solve the rubber cylinder problem under action of forces rather than pre-specified displacements. To this end the curves with finer intervals of the axial extension can be drown following the same approach. As an example, a loading case with axial force P = 100 kN and internal pressure p = 0.03 MPa is considered. The first step is to draw lines in Figures 6 and 7 to get all possible deformation data under the load P or p, as shown in Figures 6 and 7. The obtained data are listed in Table 2.  As an example, a loading case with axial force P = 100 kN and internal pressure p = 0.03 MPa is considered. The first step is to draw lines in Figures 6 and 7 to get all possible deformation data under the load P or p, as shown in Figures 6 and 7. The obtained data are listed in Table 2.       In the second step, two lines are drawn based on the data of Table 2, as shown in Figure 8. The answer is located at the intersection point of the two lines, which satisfies the loading P = 100 kN and p = 0.03 MPa simultaneously.
In the second step, two lines are drawn based on the data of Table 2, as shown in Figure 8. The answer is located at the intersection point of the two lines, which satisfies the loading P = 100 kN and p = 0.03 MPa simultaneously.

Evaluation of ABAQUS and ANSYS
The above results are then used to evaluate FEM codes: ABAQUS and ANSYS. Taking account of uniformity of axial deformation and axisymmetry, the FEM mesh is as shown in Figure 9. In the axial direction, only one element is enough; in the hoop direction, the 2-D axisymmetric elements are adopted. Pre-specified axial displacements are applied on all nodes and pre-specified radial displacements are applied on the two nodes on the inner surface.  Table 3 lists the comparison of the stresses from exact, ABAQUS and ANSYS. The data proved that the two commercial packages provide answers of excellent accuracy.

Thick-Walled Rubber Balloon
A thick-walled rubber balloon under internal pressure is shown in Figure 10. The dimensions and loads are denoted also. The symbols p denotes internal pressure. R0, Ri and T represent the outer radius, inner radius and thickness of the rubber balloon; R is an arbitrary radius; D is the distance

Evaluation of ABAQUS and ANSYS
The above results are then used to evaluate FEM codes: ABAQUS and ANSYS. Taking account of uniformity of axial deformation and axisymmetry, the FEM mesh is as shown in Figure 9. In the axial direction, only one element is enough; in the hoop direction, the 2-D axisymmetric elements are adopted. Pre-specified axial displacements are applied on all nodes and pre-specified radial displacements are applied on the two nodes on the inner surface. In the second step, two lines are drawn based on the data of Table 2, as shown in Figure 8. The answer is located at the intersection point of the two lines, which satisfies the loading P = 100 kN and p = 0.03 MPa simultaneously.

Evaluation of ABAQUS and ANSYS
The above results are then used to evaluate FEM codes: ABAQUS and ANSYS. Taking account of uniformity of axial deformation and axisymmetry, the FEM mesh is as shown in Figure 9. In the axial direction, only one element is enough; in the hoop direction, the 2-D axisymmetric elements are adopted. Pre-specified axial displacements are applied on all nodes and pre-specified radial displacements are applied on the two nodes on the inner surface.  Table 3 lists the comparison of the stresses from exact, ABAQUS and ANSYS. The data proved that the two commercial packages provide answers of excellent accuracy.

Thick-Walled Rubber Balloon
A thick-walled rubber balloon under internal pressure is shown in Figure 10. The dimensions and loads are denoted also. The symbols p denotes internal pressure. R0, Ri and T represent the outer radius, inner radius and thickness of the rubber balloon; R is an arbitrary radius; D is the distance  Table 3 lists the comparison of the stresses from exact, ABAQUS and ANSYS. The data proved that the two commercial packages provide answers of excellent accuracy.

Thick-Walled Rubber Balloon
A thick-walled rubber balloon under internal pressure is shown in Figure 10. The dimensions and loads are denoted also. The symbols p denotes internal pressure. R 0 , R i and T represent the outer radius, inner radius and thickness of the rubber balloon; R is an arbitrary radius; D is the distance from an arbitrary point to the inner wall. The deformation mode is quite simple: Uniform radial expansion.
Besides, under the consideration of full incompressibility, the deformation is completely defined by only one parameter: Radial expansion of the inner or outer radius.
Materials 2020, 13, x FOR PEER REVIEW 9 of 17 from an arbitrary point to the inner wall. The deformation mode is quite simple: Uniform radial expansion. Besides, under the consideration of full incompressibility, the deformation is completely defined by only one parameter: Radial expansion of the inner or outer radius. Figure 10. A thick-walled balloon subjected to internal pressure.

Deformation Analysis
Let ri, r0 be the inner and outer radius after deformation, respectively. Due to full incompressibility there exists: Thus, if r0 is already know: Or, if ri is already know: At an arbitrary radial position, the radius before and after deformation are: Therefore, the hoop stretch in both meridional and circumferential directions is: In the radial direction the principle stretch is: The strain invariants are:

Deformation Analysis
Let r i , r 0 be the inner and outer radius after deformation, respectively. Due to full incompressibility there exists: Thus, if r 0 is already know: Or, if r i is already know: At an arbitrary radial position, the radius before and after deformation are: and Therefore, the hoop stretch in both meridional and circumferential directions is: In the radial direction the principle stretch is: The strain invariants are:

Stress Under Pre-Specified Deformation
Suppose the balloon is under action of a pre-specified expansion of the inner radius: R i →r i . All the strain quantities can then be calculated based on Equations (34)-(39).
Similar to the cylinder problem the stress calculation can start from examining the radial differential Equation. Figure 11 shows an infinitesimal slice cut out from the balloon. Examining the radial equilibrium results in the following Equation: Materials 2020, 13, x FOR PEER REVIEW 10 of 17

Stress Under Pre-Specified Deformation
Suppose the balloon is under action of a pre-specified expansion of the inner radius: Ri→ri. All the strain quantities can then be calculated based on Equations (34)-(39).
Similar to the cylinder problem the stress calculation can start from examining the radial differential Equation. Figure 11 shows an infinitesimal slice cut out from the balloon. Examining the radial equilibrium results in the following Equation: Figure 11. Radial equilibrium of a thick-walled balloon.
By integration along the whole thickness, the internal pressure inducing the pre-specified expansion is: In which, Figure 11. Radial equilibrium of a thick-walled balloon.
Substituting the constitutive Equation (41) into (43), one can obtain: By integration along the whole thickness, the internal pressure inducing the pre-specified expansion is: In which, Materials 2020, 13, 885

of 16
By integration along a part of the thickness, the true radial stress is obtained.
In which, Combining Equation (47) with (41), the hoop or tangential stress is

Numerical Example
A thick-walled rubber balloon under internal pressure was established to verify the correctness of the analytical solution. As an example, the rubber balloon with the following parameters (Table 4) were computed: The problem is computed by Equations (45), (47) and (49). Figure 12 shows the stress distribution in the thickness direction. Figure 13 shows the curve of the internal pressure versus radial expansion. It can also serve as a graphic tool calculating radial expansion of the balloon under action of known internal pressure. By integration along a part of the thickness, the true radial stress is obtained.
In which, Combining Equation (47) with (41), the hoop or tangential stress is

Numerical Example
A thick-walled rubber balloon under internal pressure was established to verify the correctness of the analytical solution. As an example, the rubber balloon with the following parameters (Table 4) were computed: The problem is computed by Equations (45), (47) and (49). Figure 12 shows the stress distribution in the thickness direction. Figure 13 shows the curve of the internal pressure versus radial expansion. It can also serve as a graphic tool calculating radial expansion of the balloon under action of known internal pressure.

Evaluation of ABAQUS and ANSYS
The example has been calculated by ABAQUS and ANSYS. The adopted mesh is quite similar to the one used in the cylinder problem studied in the last section. But the cross section is now a trapezoid rather than a rectangle, as shown in Figure 14. All the nodes are fixed tangentially. Pre-

Evaluation of ABAQUS and ANSYS
The example has been calculated by ABAQUS and ANSYS. The adopted mesh is quite similar to the one used in the cylinder problem studied in the last section. But the cross section is now a trapezoid rather than a rectangle, as shown in Figure 14. All the nodes are fixed tangentially. Pre-specified radial displacements are applied on the two nodes on the inner surface.

Evaluation of ABAQUS and ANSYS
The example has been calculated by ABAQUS and ANSYS. The adopted mesh is quite similar to the one used in the cylinder problem studied in the last section. But the cross section is now a trapezoid rather than a rectangle, as shown in Figure 14. All the nodes are fixed tangentially. Prespecified radial displacements are applied on the two nodes on the inner surface. Comparison of the calculated stresses is made in Table 5. It can be seen that both the software packages possess strong capabilities in hyper-elasticity analysis.

Shear Stress and Torque Under Simple Torsion
When a cylinder, either hollow or solid, is subjected to twist, the main deformation mode is shear. Let us examine a segment of a cylinder with unit length as shown in Figure 15. An arbitrary material point located on a circle of radius r will move a distance u in the hoop direction, see Comparison of the calculated stresses is made in Table 5. It can be seen that both the software packages possess strong capabilities in hyper-elasticity analysis.

Shear Stress and Torque Under Simple Torsion
When a cylinder, either hollow or solid, is subjected to twist, the main deformation mode is shear. Let us examine a segment of a cylinder with unit length as shown in Figure 15. An arbitrary material point located on a circle of radius r will move a distance u in the hoop direction, see Figure 15a. Hence a rectangular element on the cylindrical surface will develop the deformation mode of well known "simple shear", as shown in Figure 15b. The expression of shear stress can then be copied from the relevant references.
Materials 2020, 13, x FOR PEER REVIEW 13 of 17 15a. Hence a rectangular element on the cylindrical surface will develop the deformation mode of well known "simple shear", as shown in Figure 15b. The expression of shear stress can then be copied from the relevant references. Let the relative twist angle be denoted by θ then from Figure 15a: By integration, the torque is: By substitution of Equations (7) and (8) into (53) with M = m 4 m 4 2m 6 θ 2 2m 6 θ 2 2m 6 θ 2 · · · m k = 1

Poynting Effect
It must be pointed out that there is also normal stress acting in the plane of the shear stress due to the so called "Poynting effect". An elaborate analysis of all the problems can be found in, for instance, Reference [20,21]. The final forms of the solutions are as follows.
Normal force under simple torsion: Torque and normal force under combined tension-torsion: with Substitution of Equations (7) and (8) into (56) and (58) yields the following final expressions. Normal force under simple torsion: Torque under tension-torsion: Normal force under tension-torsion:

Numerical Example and Evaluation of ABAQUS
Examine a solid or hollow shaft with the following given data (units are in mm, MPa.): Outer radius 10, inner radius 0 or 2, length l = 1, material constants are the same with that used in examples of the forgoing sections, relative twist angle = 5 • . The problem was calculated by the analytical solution expressed by Equations (54), (60) and ABAQUS. For ABAQUS the asymmetric element CGAX4H was used to model the shaft. The computed shear stress, torque and normal force are listed in Table 6, which exhibit excellent agreement between the two solutions. The solid cylinder was also computed under a combined tension-torsion with axial stretch equal to 1.5. Results of the torque and normal force are listed in Table 7, which also show good correlation between the analytical solution and answers from ABAQUS.

Other Forms of the Strain Energy Function
Even though only the polynomial function is discussed in the foregoing sections, the benchmark problems combined with other forms of the strain energy function can be solved in the same manner. First, we should notice that some other functions are nothing but subsets of the polynomial form. In other words, they are incomplete polynomial expressions with some terms missing. This is summarized in Table 8.
By combining the table with Equations (3)-(6), as well as all the related subsequent expressions, the formulation of other strain energy functions can be derived.

Conclusions
The benchmark problems, thick-walled rubber cylinders under axial extension and/or radial expansion, thick-walled rubber balloons under radial expansion and solid or hollow rubber cylinders under torsion or tension-torsion have been studied. For the strain energy function of generalized polynomial form, the paper gives exact solutions of the problems. Numerical examples and evaluation results of two commercial packages that are in common use are presented. The results were then used to evaluate FEM analysis capabilities of ABAQUS and ANSYS in the field of hyper-elasticity. Good agreements are found in the comparisons between the present exact standard solutions and the simulation results, which can serve as tools to evaluate availabilities of hyper-elasticity analysis in FEM packages or other approximate solutions.