Comparative study on the sealing performance of packer rubber based on elastic and hyperelastic analyses using various constitutive models

Evaluation of the sealing performance of the packer rubber varies according to specific simulation models. This paper aims at revealing the difference between elastic and hyperelastic analyses using the finite element mothed (FEM). The study extracts the hyperelastic parameters of the neo-Hookean, the Mooney-Rivlin, and the Yeoh models from a uniaxial tensile test. Then, the setting process of a mechanical packer is simulated by elastic and hyperelastic calculations. We compare the deformed configuration and the contact stress given by these models. Our results show that the Yeoh model produces the minimum residual sum of squares (RSS) among the hyperelastic models for hydrogenated nitrile butadiene rubber (HNBR). The mooney-Rivlin model has a negative parameter, making the calculation unstable. The linear elastic model fails to simulate the setting process, while the neo-Hookean model overestimates the contact stress. Despite the similar stress distribution, the nonlinear elastic model provides a 17.8% higher contact stress on average than the Yeoh model. A parametric study based on the Yeoh model points out that the sub-thickness of the packer rubber needs an elaborate design. Reducing the sub-thickness could increase the contact stress but decrease the seal length in the force control mode. From an engineering perspective, this study demonstrates that it needs to pay more attention when selecting an appropriate material model and a sound analysis method to evaluate the sealing performance of an oil packer.


Introduction
An oil packer is set in the annular between the tubing and casing or between the casing and formation in a wellbore, isolating produced fluids and standing the pressure difference [1,2]. The primary function of the packer is to protect the tubing or casing, including preventing corrosion from the formation fluid, supporting the tubing weight, and constraining the tubing movement [1,3]. The rubber sleeve is the core component of the packer, which serves as the principal sealing element [4,5]. With the increasing growth of unconventional oil and gas exploration, the downhole condition becomes more and more adverse, imposing an enormous challenge on the packer rubber [6][7][8]. For instance, the wellbore pressure is higher than 70 MPa, and the formation temperature is over 150°C in high pressure and high temperature (HPHT) wells, often encountered in deep-water fields [9]. This extremely harsh environment needs a robust seal in the annular that depends on the mechanical behavior of the rubber. Therefore, evaluating the sealing performance is essential in the packer design.
Analytical and numerical models are valuable to assess the seal quality of the packer before manufacture, device selection, and putting it into use [5]. Most of these models calculate the contact stress at the interface between the rubber and casing (or formation) as a sealing indicator. Analytical models usually regard the setting Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. process as a rubber sleeve compressed confined by two metal tubes or by a metal tube and the surrounding formation [10][11][12][13]. They pursue a closed-form solution describing the deformation and stress of the rubber. Al-Kharusi et al [10] proposed a mechanical model to study the performance evaluation of swelling elastomer seals. They assumed that the elastomer was linear elastic and the axial displacement had a parabolic distribution. The axial stress component in elastomer was expressed in terms of pressures at the two ends and the annular geometry. Using similar assumptions, Al-Hiddabi et al [11] gave a more sophisticated formula for pressure in the seal. Zhang et al [13] took advantage of the Love displacement function to determine the axial displacement for compression packers. The derived formula calculated the maximum shear stress in the rubber as a mechanical failure criterion. These analytical models agree with numerical simulations, but they all rely on linear elastic constitutive equations.
One prominent advantage of analytical models is their explicit formulae that could quickly give an estimation. However, their limitations due to restricted assumptions (linear elastic or small deformation) constrain their use, especially for complex geometries and finite deformation. Numerical models could conquer this difficulty by a nonlinear analysis. The finite element method (FEM) is one of the most popular numerical techniques applied to engineering problems [14]. Numerous FEM simulations on the packer setting have been reported [15][16][17][18][19][20]. Alzebdeh et al [15] conducted a FEM simulation to investigate the effects of the seal length and thickness, the compression ratio, and boundary conditions on the contact pressure for an expandable elastomer. They demonstrated that the contact pressure would reach a stable value when the seal length was longer than 200 mm. Hu et al [17] analyzed the influence of rubber materials on the seal quality for compression packers by FEM calculations. In their conclusion, the Ogden model was better to characterize the mechanical behavior of the investigated rubber. By setting up a three-dimension FEM model, Ahmed et al [18] indicated that the most critical factor affecting the performance of an elastomer seal is the annular fit, followed by the Poisson's ratio and elastic modulus of the rubber, and the compression ratio. Their model is linearly elastic and validated by the analytical formula from Al-Hiddabi and his co-workers' results [11]. Zheng et al [19] presented a nonlinear FEM simulation on the sealing behavior of the packer rubber with the Mooney-Rivlin constitutive model. They optimized the total thickness, the sub-thickness, the height, and the material hardness of the rubber by parametric studies.
There still exist some issues in assessing the packer seal from analytical and numerical models. Packer rubber is usually made of hyperelastic materials such as nitrile butadiene rubber (NBR), hydrogenated nitrile butadiene rubber (HNBR), ethylene propylene diene monomer (EPDM), etc [21,22]. These materials have polymer structures whose mechanical behavior needs hyperelastic constitutive models to describe [23]. Analytical models neglect the nonlinear response of hyperelastic polymers to simplify the solution, which distorts the predicted values in some cases [24]. On the other hand, numerical models may offer diverse or even contrary results by different constitutive models and analysis methods. For example, the previous study showed that the contact stress maximized at the upper end of the rubber and decreased near the lower end [19]. At the same time, another FEM analysis pointed out that the contact stress was evenly distributed along with the packer [25]. Therefore, it is necessary to identify the difference between simulations by various constitutive models and analysis methods.
This study compares FEM simulations using different modeling methods, including linear elastic, nonlinear elastic, and hyperelastic analysis. The hyperelastic parameters are extracted from a uniaxial test for HNBR. We calculate the contact stress at the sealing interface and the reaction force during the setting to illustrate the difference between analysis strategies. Moreover, a parametric study reveals the effects of the height, the total thickness, and the sub-thickness of the rubber. Results are conductive to evaluate the sealing performance of the packer rubber accurately.

Constitutive models for hyperelastic materials
It is essential to use an appropriate constitutive model to predict the mechanical behavior of rubber-like materials. Continuum mechanics supplies a fundamental theory to describe the constitutive relationship for hyperelastic materials [26]. The proposed model is very complicated for generic materials. Fortunately, for isotropic and incompressible materials, e.g., polymers, the constitutive equation has a simplified form as [26]: where σ represents the Cauchy stress tensor, b is the left Cauchy-Green strain tensor, Ψ is the strain-energy function, I denotes the unit tensor, and p is a Lagrange multiplier that can be identified as hydrostatic pressure. The remained terms I 1 and I 2 are the invariants of b, which are related to the stretch ratio as [26]: are the stretch ratios in three mutual orthogonal directions.
Equation (1) shows that the strain energy function determines the constitutive relationship between the stress and strain. Many models have been proposed to characterize Ψ. The simplest one is the neo-Hookean model [27]: where μ is the shear modulus. It often carries out a uniaxial test to plot stress-strain curves and then fits the curve using equation (3) to estimate μ. For instance, under uniaxial tension, the strain invariants become: where λ is the uniaxial stretch ratio. Substituting equation (3) and (4) into equation (1) gives the uniaxial stress: where σ uni is the stress in the uniaxial test. Once the σ uni -λ curve is available, μ can be easily extracted from equation (5). It should be noted that σ uni is the true stress in the uniaxial test, which can be expressed in terms of the engineering stress and strain as [28]: where σ 0 and ε 0 represent the engineering stress and strain, respectively, and the stretch ratio λ equals 1+ε 0 .
Although the neo-Hookean model is simple, it becomes awkward for complicated polymers. More sophisticated models have emerged to characterize the mechanical behavior of hyperelastic materials. Among them, the Mooney-Rivlin and Yeoh models are two relatively concise ones. The mooney-Rivlin model introduces a linear dependence on I 2 in the strain energy function [27]: where C 10 and C 01 are two parameters having the same unit with μ. Equation (7) makes equation (5) into the following form: Unlike the Mooney-Rivlin model, the Yeoh Model assumes that the strain energy function is independent of I 2 [27]: Hassani et al [29] conducted uniaxial tensile tests for the HNBR compound at a 50 mm min −1 constant speed for different temperatures. The results show that HNBR has a stable tensile property without the aging effect. Using their data at 75°C, we compare the performance of the neo-Hookean, mooney-Rivlin, and Yeoh models, as illustrated in figure 1. The mooney-Rivlin and Yeoh models outperform the neo-Hookean model since they are more consistent with the experimental data. It is difficult to identify which is better between the Mooney-Rivlin and Yeoh models only through observation in figure 1. Then, we list the model parameters and the residual sum of squares (RSS) for different fitting curves [30]. The RSS of the Yeoh model is smaller than the one of the Mooney-Rivlin model. Besides, the Mooney-Rivlin model has another inconvenience in that C 01 is negative and has the same order as C 10 , which would make the numerical simulation unstable. Hence, this study selects the neo-Hookean and Yeoh models to describe the mechanical behavior of rubber and compare the results with the elastic model. Table 1 also presents an elastic modulus of E, which is estimated by the Yeoh model for incompressible materials [31]: In the following sections, we will use E to simulate the setting process of the packer by elastic analysis.

Mechanical model for packer setting
A typical unit in the conventional mechanical oil packer is composed of the rubber, push ring, thrust ring, and the annular between the central pipe and casing, as sketched in figure 2. When the packer reaches the desired position in the wellbore, the push ring moves down and compresses the rubber. As the axial length shortens, the  rubber expands in the radial direction until it contacts the casing inwall. A prescribed compression induces enough contact stress between the packer and casing to provide a robust seal. Boundary conditions include loading and displacement pre-setting. This study considers two loading modes: displacement control and force control. The displacement control applies a downward vertical displacement (see v c in figure 2) on the push ring to compress the rubber. In the force control, we set a compression load (with the same value of F r but has an opposite direction in figure 2(b)) on the push ring. The thrust ring, central pipe, and casing are fixed during the compression. The rubber and other metal parts connect by surface-to-surface contacts.
To simplify the simulation for the packer setting, we make some primary assumptions as follows: (1) The model is axisymmetric due to the symmetry of the geometry and boundary conditions. A cylindrical coordinate system (r, z) locates in the center of the wellbore with the z-axis parallel to the well depth, as shown in figure 2(a).
(2) The packer rubber is isotropic, incompressible, and elastic or hyperelastic with the material parameters listed in table 1.
(3) Since the metal parts (the rings, central pipe, and casing) are much harder than the rubber, they are deemed as rigid bodies.
Geometry dimensions involved in figure 2 contain the height (h), thickness (t), and outer radius (R 1 ) of the rubber, the outer radius (R 0 ) of the central pipe, and the inner radius (R 2 ) of the casing. This study also considers the sub-thickness (t s ) and fixes the chamfer angle as 45°. We undertake the commercial FEM software ABAQUS to simulate the setting process of the packer. This software has built-in elastic and hyperelastic material constitutive models, and the required input parameters are listed in table 1. The rubber is meshed with 4-node bilinear axisymmetric quadrilateral and hybrid elements (CAX4H) to accommodate the incompressibility of the material, as illustrated in figure 2(c). This element has four nodes, each having two degrees of freedom. As the deformation at the shoulders of the rubber is more significant than the other parts, the two ends of the packer have meshed with refined elements. The model constrains all degrees of freedom in the thrust ring, the central pipe, and the casing. A vertical displacement v c applies to the push ring to compress the rubber (see figure 2(b)). Surface-to-surface contacts establish between the rubber and metal parts with a frictional coefficient of 0.2. The reaction force between the rubber top and the push ring F r and the contact stress σ c at the rubber-casing interface are recorded during the simulation. The contact stress characterizes the sealing performance of the packer, and the reaction force reveals the relationship between the set-down load and the contact stress.
As a supplement to the above model, we present the equilibrium equation of stress components for axisymmetric problems [32]: Figure 3. The contact stress calculated by different element sizes.
where σ rr , σ θθ , and σ zz are the stress components in the radial, hoop, and axial directions, respectively, and τ rz denotes the tangential stress. By the current coordinates in figure 2(a), the magnitude of σ rr at the right side of the rubber is the contact stress. It needs to emphasize again that the stress in Eq. (12) is the Cauchy stress, i.e., the true stress. Equation (1) shows that the true stress should be in terms of the left Cauchy-Green strain tensor, which is measured in the deformed configuration. This measurement is the fundamental difference between the hyperelastic analysis and the linear elastic theory. The latter often supposes that the deformed configuration coincides with the reference frame due to the small deformation assumption. In the next section, we will illustrate this difference in detail.

Base case study
We select a typical geometry of the packer for the base case study, as listed in table 2. Before comparing the results of elastic and hyperelastic models, we first investigate the dependence of the simulation on the element size. Figure 3 shows the contact stress σ c estimated by different element sizes with the Yeoh model. As the element size decreases, the curves in the middle of the contact interface (about-15∼15 mm) are close. However, the stress at the two ends of the packer still disperses. Extreme fine elements would produce a reliable FEM simulation result with a long computational time price. Table 3 lists the maximum contact stress obtained by increasing the element number for the current case. More refined elements result in minor contact stresses, but the maximum  does not decrease monotonically. Thus, we make a tradeoff to set the element number as 614 with 1456 degrees of freedom, which can make a convergent estimation and save computation time.
After fixing the element size, the simulation repeats with four analysis strategies: the linear elastic model, the nonlinear elastic model, the neo-Hookean model, and the Yeoh model. Figure 4 illustrates the final deformed configuration when the push ring moves down 9.5 mm. The linear elastic analysis indicates that the right side of the packer does not reach the casing ( figure 4(b)), which means a failed setting. It infers that the linear elastic model makes the rubber too rigid to be compressed. In contrast, the nonlinear elastic, neo-Hookean, and Yeoh models mimic the setting process successfully in that the rubber contacts the casing after compression. The deformed configurations of the nonlinear elastic model ( figure 4(c)) and the Yeoh model ( figure 4(e)) are like each other. Figure 4(d) shows that the neo-Hookean model induces milder deformation, making the rubber more rigid than the other two models.
As illustrated in figure 4, the linear elastic model cannot simulate the setting process. Then, we compare the contact stress estimated by the nonlinear elastic, neo-Hookean, and Yeol models, as shown in figure 5. The similarity between the three models is the stress distribution that σ c increases from the bottom to the top. There is a stress drop at both ends of the packer, and the maximum contact stress occurs near the upper end. Previous studies have reported a similar stress distribution [19]. As expected, the neo-Hookean model produces higher σ c since it makes the rubber stiffer. In particular, the maximum contact stress estimated by the neo-Hookean model is 4.03 MPa for v c =8.5 mm ( figure 5(a)). For the same compression displacement, the nonlinear elastic and Yeoh models give 2.81 MPa and 2.57 MPa, respectively. The maximum σ c increases to 7.49 MPa, 6.58 MPa, and 5.75 MPa by the neo-Hookean, nonlinear elastic, and Yeoh models, respectively, when v c is 9.5 mm  ( figure 5(b)). It is worth noting that the stress curves from the nonlinear elastic and Yeoh modes are nearly parallel, which indicates that the difference between these two models is mainly in magnitude. It may attribute to the elastic modulus E obtained from the Yeoh model (equation (11)). The nonlinear elastic model overestimates the contact stress, and its prediction is 11.1% and 17.8% higher on average than the Yeoh model for 8.5 mm and 9.5 mm compressions, respectively. As figure 5 shows, more compression will widen the gap between the nonlinear elastic analysis and the Yeoh model. This discrepancy suggests that the Yeoh model is more suitable to describe the mechanical behavior of rubber for evaluating the performance of a packer.
We focus on the reaction force (F r ) at the push ring to present the difference between the elastic and hyperelastic models more clearly. Figure 6 plots F r as a function of the compression displacement from the nonlinear elastic and Yeoh models. It also identifies a boundary of v c (8.18 mm), over which the relative deviation δ between these two models will surpass 5%. This boundary corresponds to a compression ratio of 15.6%. The subgraphs in figure 6 are the rubber configurations during the setting process. They show that before the packer finishes setting, δ has been greater than 5%, which reconfirms the competence of the Yeoh model.

Sensitivity analysis
Section 4.1 has demonstrated that the nonlinear elastic model would overvalue the contact stress, which may give an illusion to engineers that the seal is strong enough. Conversely, the Yeoh model is more conservative and provides a safe estimation for engineering applications. This section will reveal the impact of the packer geometry on the contact stress based on the Yeoh model. The parametric study sets two analysis modes: one fixes the compression displacement v c as 9.5 mm, and the other controls the reaction force F r as 42.0 kN. The magnitude of F r is extracted from the base case study (see figure 6). Figure 7 tells that shorter rubber could generate greater contact stresses in both control modes. For v c =9.5 mm ( figure 7(a)), as the rubber height  reduces from 62.5 mm to 42.5 mm, the maximum contact stress will increase from 2.36 MPa to 12.29 MPa. The maximum σ c stays almost stable for different rubber heights in the force control mode. The contact stress in the middle of the rubber increases with the reduction of h. Shorter rubber would undergo more deformation with the same compression displacement or force, leading to more significant contract stress. Moreover, the rubber height becomes a more sensitive parameter in the displacement mode.
Thicker rubber would induce higher contact stress both in the displacement and force control modes, as shown in figure 8. Under displacement control, as t increases from 16.2 mm to 18.2 mm, the maximum σ c goes from 0.68 MPa to 14.37 MPa ( figure 8(a)). Note that a 0.68 MPa contact stress will not provide a reliable seal between the packer and the casing. 1 mm thicker rubber than the base case would enhance σ c by 20.2% in the force control mode ( figure 8(b)). It must be pointed out that increasing the thickness may make the packer problematic to run in the hole since the gap between the rubber and the casing decreases, as sketched in figure 2. Besides, thicker rubber also costs more, then increasing t is hardly a wise decision.
Generally, the annular gap between the tubing and casing is fixed once the well has been built, constraining the total thickness of the rubber. The other parts of a packer also limit the rubber height. Hence, adjusting the sub-thickness of the rubber is more practical than changing its height or total thickness. Like the effects of the rubber height and total thickness, figure 9 indicates that the contact stress is more sensitive to the sub-thickness in the displacement mode. Unlike h and t, σ c increases with t s under the displacement control but decreases with t s under the force control. For instance, the maximum σ c grows from 1.83 MPa to 11.10 MPa when t s varies from 8 mm to 12 mm in the displacement mode ( figure 9(a)). Meanwhile, the maximum σ c reduces from 6.03 MPa to 5.67 MPa for the same change of t s in the force mode ( figure 9(b)). A smaller t s needs a more insignificant force to be compressed to a certain amount, leading to less contact stress. Oppositely, under the force control, reducing t s makes the rubber compressed more to provide greater contact stress.  In engineering practice, a mechanical oil packer is often subjected to a prescribed compression load. Then, it is reliable and economical to decrease the sub-thickness of the rubber to ensure a good seal. Nevertheless, modifying t s may cause other problems. Figure 10 plots the radial stress component of the rubber in the deformed configuration for different sub-thicknesses under the force control. Reducing t s would decrease the contact area in that the seal length is shorter in figure 10(a) than the one in figure 10(c). On the other hand, the longer t s would induce extrusion of the upper shoulder, leading to stress concentration and rubber damage [33]. Therefore, it needs an elaborate consideration on t s .
To illustrate the cons and pros of t s thoroughly, we define a seal length, l seal , as marked in figure 2(b). l seal is estimated from the FEM simulation as the following procedure. Extract the nodal displacements and the original coordinates of the elements on the potential contact area, as labeled in figure 2(a). Let u re and u ze denote the radial and axial displacements, respectively, and r e and z e are the undeformed radial and axial coordinates, respectively. It identifies that the corresponding element contacts the casing if:

+ >
The factor on equation (13) right hand is slightly less than 1 to accommodate numerical errors. Then the seal length is given by: Using the maximum contact stress and the seal length, figure 11 shows a comprehensive assessment of the effects of t s . The maximum σ c decreases nonlinearly with t s , while l seal presents a linear increase. It recommends an optimal value of t s as 9.6 mm, which could guarantee high sealing stress with a large contact area.

Conclusions
This paper has presented a comparative study on the sealing performance of the packer rubber based on different constitutive models and analysis strategies. The considered hyperelastic constitutive equation includes the neo-Hookean, mooney-Rivlin, and Yeoh models. The analysis methods are linear elastic, nonlinear elastic, and hyperelastic. Remarkable conclusions from the FEM simulation results are as follows: (1) Parameter fitting for a uniaxial tensile test of HNBR shows that the Yeoh model is more reliable and accurate. It produces the minimum RSS among the three investigated hyperelastic constitutive equations.
The mooney-Rivlin model induces a negative material parameter that would result in an unstable FEM simulation.
(2) The linear elastic model makes the rubber too stiff to simulate the setting process of the packer. The neo-Hookean model gives higher contact stress than the nonlinear elastic and Yeoh models. Although the nonlinear elastic analysis and the Yeoh model provide a similar stress distribution, the former would overestimate the contact stress 17.8% higher on average than the latter. (3) Parametric studies on the effects of the rubber geometry from the Yeoh model indicate that the contact stress is more insensitive to the height, thickness, and sub-thickness under the force control. Tuning the sub-thickness is more practical than adjusting the other two sizes. The study suggests an optimal subthickness, which could supply great contact stress and a large seal area.