Sage Journals: Discover world-class research

Abstract

A finite element procedure is presented for the analysis of rubber-like hyperelastic materials. The volumetric incompressibility condition of rubber deformation is included in the formulation using the penalty method, while the principle of virtual work is used to derive a nonlinear finite element equation for the large displacement problem that is presented in a total-Lagrangian description. The behavior of rubber deformation is represented by hyperelastic constitutive relations based on a generalized Mooney-Rivlin model. The proposed finite element procedure using analytic differentiation exhibited results that matched very well with those from the well-known commercial packages NISA II and ABAQUS. Furthermore, the convergence of equilibrium iteration is quite slow or frequently fails in the case of near-incompressible rubber. To prevent such phenomenon even for the case that Poisson's ratio is very close to 0.5, Poisson's ratio of 0.49000 is used, first, to get an approximate solution without any difficulty; then the applied load is maintained and Poisson's ratio is increased to 0.49999 following a proposed pattern and adopting a technique of relaxation by monitoring the convergence rate. For a given Poisson ratio near 0.5, with this approach, we could reduce the number of substeps considerably.

1. Introduction

Hyperelastic rubbers are widely used in industry as sealants, load-bearing parts, and damping parts and in other similar applications. However, the incompressibility and nonlinearity of rubber cause many difficulties, one of which is divergence problem, in analytic and numerical studies like the finite element method [1 –7] as noted in the literature [7]. Furthermore, the problems of self-buckling, follower forces, and contact conditions cause additional difficulties in analyzing rubber. Accordingly, this paper attempts to contribute to solving these problems by proposing a rubber analysis procedure using a numerical differentiation method and Poisson iteration method.

Rubber is a hyperelastic material and its characteristics can be expressed using a strain energy function, which can then be used to derive stresses and constitutive matrices. The strain energy function is a function of the material property and nonlinear strain. There are many models [8] describing the rubber materials such as Mooney-Rivlin model, Neo-Hookean model, Ogden form, and Yeo form. Among them, the general Mooney-Rivlin model is adopted to describe the material characteristics [9 –11]. Yet, the next step for a constitutive matrix involves a quite complex derivation for each Mooney-Rivlin model selected, for example, 2-, 3-, 5-, 6-, or 9-coefficient model, or even more. As such, this requirement for each material model of incompressibility is a formidable process and prone to erroneous properties, especially when the number of Moony-Rivlin coefficients is large. Consequently, a numerical derivation of the constitutive matrices is proposed using a finite difference-type numerical differentiation.

In this paper, a finite element equation is derived for hyperelastic materials based on the theory of virtual work, while considering a large displacement and the deformation characteristics. The incompressibility of rubber is then considered using the penalty function method [12]. The strain energy function involved uses a generalized polynomial of the Mooney-Rivlin model that includes an extended simple Mooney-Rivlin material. The constitutive matrices, numerically obtained from the partial differentiation of the strain energy function, are then included into a nonlinear finite element procedure for an isoparametric 3-dimensional 20-node element. The results obtained with the numerical differentiation method are compared with those obtained with the analytic differential method and the results from commercial packages [8, 13], such as ABAQUS, thereby confirming the validity of the proposed rubber finite element program using the proposed numerical method.

An attempt is also made to analyze the behavior of near-incompressible rubber materials and solve the problem of convergence when Poisson's ratio approaches 0.5. At this point the convergence rate of the equilibrium iteration becomes very slow, and in many cases the iterations diverge. Therefore, a robust method of nondivergent equilibrium iteration is needed, regardless of the number of iterations. Approaching Poisson's ratio to the value of 0.5 is very similar to increasing the load in a conventional nonlinear analysis, as the volumetric strain should approach to zero regardless of the correction load with any magnitude. Thus, when the value approaches 0.5, the rate becomes slower and slower, as in the other cases of severe nonlinearity. Finally, a reduced or selective integration technique can also be used to avoid the volumetric locking phenomenon that can occur when stiffness matrices are integrated in full order.

In this iteration technique, however, it is a critical job to find proper steps of approach without the divergence. Another simple but cumbersome approach is to divide the interval (0.49000∼0.49999) into many small subintervals. This will increase a lot the computational time. Several strategies are possible in Poisson iteration when increasing Poisson's ratio interval from 0.49000 to 0.49999. We propose a method where Poisson's ratio approaches to 0.5 following a pattern that is chosen to avoid most of divergences. In addition, we adopt scale factors to adjust the proposed pattern for proper approaches reflecting current convergence rate. In each equilibrium iteration, we adopt relaxation factors to reduce instability after first substep. Through various trials, we could find proper factors that can lead to a high convergence rate and shorten analysis time.

In short, the main reasons for the development of present procedure are, for free selection of any material model with proper terms, without formidable analytical derivation of material properties, prone to the inclusion of errors. We only need to change the potential function as necessary. Next an analyzer performs a fast computation with the Poisson ratio ν = 0.49, and if necessary, he/she may restart and approach the near-incompressible state (ν≑0.5), by starting Poisson's iteration from 0.49 to 0.4999 or 0.49999. Conventional near-incompressible method can usually handle the case of Poisson's ratio ν < 0.499 (better ν < 0.495). A near-incompressible analysis which can handle the cases of Poisson ratio ν = 0.4999∼0.49999, such as the present method, can represent fully incompressible state (ν = 0.5) almost exactly. Thus present method may be thought to cover the missing range of 0.499 < ν < 0.49999, under which conventional near-incompressible method covers and over which fully incompressible method covers.

2. Finite Element Analysis of Rubber

Based on the theory of virtual work, a finite element equation for hyperelastic materials is described by (1) [6], where $[{}_{0}^{t}{K_{L}}]$ denotes a linear stiffness matrix, $[{}_{0}^{t}{K_{NL}}]$ denotes a nonlinear stiffness matrix at time t referring to the state 0, and ${{}_{0}^{t}F}$ denotes the equivalent nodal force. Material property matrix [D] and second Piola-Kirchhoff stress may be derived from the strain energy potential W and penalty function U [12]; however, if we adopt modified strain energy function equation (3) for near-incompressibility, they are given as in the following:

\begin{matrix} ([{}_{0}^{t + Δ t}{K_{L}}] + [{}_{0}^{t + Δ t}{K_{NL}}]) {Δ U}^{(i + 1)} = {{}^{t + Δ t}p} - {{}_{0}^{t + Δ t}F}^{(i + 1)}, \\ {{}_{0}^{t + Δ t}F}^{(i + 1)} = \int_{0_{V}}^{} {[{}_{0}^{t + Δ t}{B_{L}}]}^{T} [{}_{0}^{t + Δ t}S] d ({}^{°}V), \end{matrix}

(1)

where

\begin{matrix} [{}_{0}^{t + Δ t}{K_{L}}] = \int_{0_{V}}^{} {[{}_{0}^{t + Δ t}{B_{L}}]}^{T} [D] [{}_{0}^{t + Δ t}{B_{L}}] d ({}^{°}V), \\ [D] = [\frac{\partial^{2} W}{\partial ε_{i j} \partial ε_{k l}} + \frac{\partial}{\partial ε_{k l}} (\frac{U}{α} \frac{\partial U}{\partial ε_{i j}})], \\ [{}_{0}^{t + Δ t}{K_{NL}}] = \int_{0_{V}}^{} {[{}_{0}^{t + Δ t}{B_{NL}}]}^{T} [{}_{0}^{t + Δ t}S] [{}_{0}^{t + Δ t}{B_{NL}}] d ({}^{°}V), \\ [{}_{0}^{t + Δ t}S] = [{}_{0}^{t + Δ t}{(\frac{\partial W}{\partial ε_{i j}} + \frac{U}{α} \frac{\partial U}{\partial ε_{i j}})}] . \end{matrix}

(2)

2.1. Strain Energy Function

To express the incompressibility of a hyperelastic material, such as rubber, the penalty term P as a function of penalty function U and penalty coefficient α are added to the strain energy function. The strain energy function W is represented by a generalized Mooney-Rivlin model that includes a simple Mooney-Rivlin model with stress invariants:

W^{'} = W + P = \sum_{r = 0}^{\infty} ‍ \sum_{s = 0}^{\infty} C_{r s} {(I_{1} - 3)}^{r} {(I_{2} - 3)}^{s} + \frac{1}{2 α} U^{2} (I_{3}) .

(3)

For the 2-coefficient Mooney-Rivlin model, a strain energy function with material constants for two coefficients (C₁₀, C₀₁) is used. As such, a modified 2-coefficient strain energy function for incompressibility behavior can then be derived using the penalty function, as in (7). Consider

W^{'} = C_{10} (I_{1} - 3) + C_{01} (I_{2} - 3) - (2 C_{10} + 4 C_{01}) U (I_{3}) + \frac{1}{2 α} U^{2} (I_{3}),

(4)

where the penalty function U may be chosen using third stress invariant I₃,

U (I_{3}) = \frac{1}{2} \ln I_{3},

(5)

the penalty coefficient α is chosen as

\frac{1}{α} = 6 (C_{10} + C_{01}) \times \frac{(1 - ν)}{(1 + ν) (1 - 2 ν)} - 4 (C_{10} + 2 C_{01}),

(6)

and the third terms are included for the stress free condition at the initial state.

When Poisson's ratio (ν) is 0.5, the limit value α approaches zero and the material property matrix diverges. Since U(I₃) → 0, I₃ → 1, α → 0 can effectively express the incompressibility of the material, the second Piola-Kirchhoff stress tensor [5] can then be expressed by the partial differentiation of the strain energy function with respect to the Green-Lagrange strain, as shown in (8):

\begin{array}{l} S_{i j} = C_{10} \frac{\partial I_{1}}{\partial ε_{i j}} + C_{01} \frac{\partial I_{2}}{\partial ε_{i j}} + (\frac{U (I_{3})}{α} - 2 C_{10} - 4 C_{01}) \frac{1}{2 I_{3}} \frac{\partial I_{3}}{\partial ε_{i j}} \\ = 2 C_{10} δ_{i j} + C_{01} (2 I_{1} δ_{i j} - 2 C_{i j}) + (\frac{1}{2 α} \ln I_{3} - 2 C_{10} - 4 C_{01}) \times (\frac{1}{2 I_{3}} e_{i m n} e_{i q r} C_{m q} C_{n r}) . \end{array}

(7)

Plus, the material property tensor D_ijkl can be derived as shown in (10):

\begin{array}{l} D_{i j k l} = D_{i j k l}^{d} + U_{i j k l} = \frac{\partial S_{i j}}{\partial ε_{k l}} \\ = 4 C_{01} δ_{i j k l} - 4 C_{01} δ_{i k} δ_{j l} - (C_{10} + 2 C_{01}) [- A + B] + \frac{1}{2 α} [\frac{A}{2} - \frac{A}{2} \ln I_{3} + \frac{B}{2} \ln I_{3}], \end{array}

(8)

where

\begin{array}{l} A = \frac{1}{I_{3}^{2}} \frac{\partial I_{3}}{\partial ε_{i j}} \frac{\partial I_{3}}{\partial ε_{k l}} \\ = \frac{e_{i m n} e_{j q r} e_{k m^{'} n^{'}} e_{l q^{'} r^{'}} C_{m q} C_{n r} C_{m^{'} q^{'}} C_{n^{'} r^{'}}}{I_{3}^{2}}, \\ B = \frac{1}{I_{3}^{2}} \frac{\partial^{2} I_{3}}{\partial ε_{i j} \partial ε_{k l}} = \frac{4 e_{i k n} e_{j l n} C_{n r}}{I_{3}} . \end{array}

(9)

The derivative of the Cauchy-Green deformation tensor invariant with respect to the Green-Lagrange strain can then be expressed as (11) and (12). Consider

\begin{matrix} \frac{\partial^{2} I_{3}}{\partial ε_{i j} \partial ε_{k l}} = \frac{1}{3} [e_{i m n} e_{j q r} (2 δ_{m k} δ_{q l} C_{n r} + C_{m q} 2 δ_{n k} δ_{r l}) + e_{\sin} e_{p j r} (2 δ_{s k} δ_{p l} C_{n r} + C_{s p} 2 δ_{n k} δ_{r l}) + e_{s m i} e_{p q j} (2 δ_{s k} δ_{p l} C_{m q} + C_{s p} 2 δ_{m k} δ_{q l})] \\ = \frac{2}{3} [e_{i k n} e_{j l r} C_{n r} + e_{i m k} e_{j q l} C_{m q} + e_{k i n} e_{l j r} C_{n r} + e_{s i k} e_{p j l} C_{s p} + e_{k m i} e_{l q j} C_{m q} + e_{s k i} e_{p l j} C_{s p}] = 4 e_{i k n} e_{j l r} C_{n r}, \end{matrix}

(10)

\begin{matrix} \frac{\partial I_{3}}{\partial ε_{i j}} \frac{\partial I_{3}}{\partial ε_{k l}} = \frac{1}{9} (e_{i m n} e_{j q r} C_{m q} C_{n r} + e_{i s n} e_{j p r} C_{s p} C_{n r} + e_{i s m} e_{j p q} C_{s p} C_{m q}) \times (e_{k m^{'} n^{'}} e_{l q^{'} r^{'}} C_{m^{'} q^{'}} C_{n^{'} r^{'}} + e_{k m^{'} n^{'}} e_{l q^{'} r^{'}} C_{m^{'} q^{'}} C_{n^{'} r^{'}} + e_{k s^{'} m^{'}} e_{l p^{'} r^{'}} C_{s^{'} q^{'}} C_{m^{'} q^{'}}) \\ = e_{i m n} e_{j q r} e_{k m^{'} n^{'}} e_{l q^{'} r^{'}} C_{m q} C_{n r} C_{m^{'} q^{'}} C_{n^{'} r^{'}} . \end{matrix}

(11)

2.2. Numerical Differentiation

As shown above, the second Piola-Kirchhoff stress tensor and material property tensor can be derived using an analytic differentiation method. However, there are many models and many orders of material properties to be derived that are formidable job and prone to error. This section suggests a method for deriving the two tensors using a numerical differentiation method [14]. We adopt central finite difference scheme for numerical differentiation. Using the central difference equations, the second Piola-Kirchhoff stress tensor can then be expressed as (13) and the material property tensor can be expressed as (14) and (15). Consider

S_{i j} = \frac{\partial W}{\partial ε_{i j}} = \frac{1}{2} \frac{W (ε_{i j} + Δ ε_{i j}) - W (ε_{i j} - Δ ε_{i j})}{Δ ε_{i j}},

(12)

D_{r r r r} = \frac{W (ε_{r r} + Δ ε_{r r}) + W (ε_{r r} - Δ ε_{r r})}{Δ ε_{r r}^{2}} - \frac{2 W (ε_{r r})}{Δ ε_{r r}^{2}},

(13)

D_{r r s s} = \frac{1}{4} \frac{W (ε_{r r} + Δ ε_{r r}, ε_{s s} + Δ ε_{s s}) + W (ε_{r r} - Δ ε_{r r})}{Δ ε_{r r}^{} Δ ε_{s s}} + \frac{1}{4} \frac{W (ε_{r r} - Δ ε_{r r}, ε_{s s} - Δ ε_{s s})}{Δ ε_{r r}^{} Δ ε_{s s}} - \frac{1}{4} \frac{W (ε_{r r} + Δ ε_{r r}, ε_{s s} + Δ ε_{s s}) + W (ε_{r r} - Δ ε_{r r})}{Δ ε_{r r}^{} Δ ε_{s s}} - \frac{1}{4} \frac{W (ε_{r r} - Δ ε_{r r}, ε_{s s} - Δ ε_{s s})}{Δ ε_{r r}^{} Δ ε_{s s}} .

(14)

2.3. Poisson Iteration Method

When Poisson's ratio approaches 0.5, a finite element analysis is unstable and frequently diverges. Therefore, to prevent this diverging phenomenon, when Poisson's ratio is very close to 0.5, a Poisson iteration method is proposed. First, the problem is analyzed using a smaller Poisson ratio than 0.5 (say 0.49000) for rapid convergence; after the convergence the analysis is restarted, and Poisson's ratio is increased in small increments to the desired value, through iterative approaches from 0.49000 to the desired value very close to 0.5. The Poisson iteration technique is very similar to SUMT (successive unconstrained minimizing technique), which is a sequential unconstrained minimization technique [15]. Poisson's ratios used in this paper are ν₁ = 0.49, ν₂ = 0.4950,…, ν_j = 0.49999 (see Figure 1).

Figure 1:

Approaching routes of load equilibrium iteration (-●- → -●-) and Poisson iteration (-●- → -■-).

The convergence for Poisson's ratio of 0.49000 can be easily achieved for each load increment as in other usual nonlinear problems. After the equilibrium of loads at each node, finite element data are saved and used later for restart. From this state the Poisson iteration is started and performed by increasing the ratio from 0.49000 to 0.49999 in small steps. There will be many equilibrium points for each increased Poisson ratio. However, data for these intermediate equilibrium states are usually not saved without any specific purpose. When the equilibrium is reached finally for the ratio of 0.49999, all the data are saved and we go back to the starting point of Poisson iteration. Data at the former equilibrium point with Poisson's ratio of 0.49000 had been stored, and we perform the restart for next load increment based on the saved status.

The equilibrium point with ν = 0.49000 works as a restarting point of the overall procedure which also is a starting point of Poisson iteration and also behaves as another restarting point for next load increment. We would get two routes of deformed configuration: one for ν = 0.4900 and the other for ν = 0.49999. If the near-incompressibility is required at a specific level of load only, one may proceed to the level for ν = 0.49000 without the Poisson iteration to reach the equilibrium with ease. From this point, then, the Poisson iteration is started to end up with the equilibrium with ν = 0.49999.

The number of iterations for each load increment is usually 3∼7, depending upon the specific problem and on the required convergence criteria. The number of iterations for the Poisson method is, however, greater than that of the load increments, usually of order 10 to 100, depending upon the severity of problem. The severe distortion of elements will cause slow convergence, and follower forces where the loads depend upon the deformed configuration of solid will also cause much slower convergence or even a divergence without so much smaller increment of Poisson's ratio. Therefore, an effective way of increasing the magnitude of Poisson's ratio is needed to be proposed for fast and reliable convergence.

The flow chart of finite element procedure with two iterations, load equilibrium iteration and Poisson iteration, is shown in Figure 2.

Figure 2:

Flow chart for the finite element analysis of near-incompressible rubber with Poisson iteration technique.

2.4. Verification of Proposed Method

To verify the proposed method, a single element of an isoparametric 3-dimensional 20-node rubber element was used, as shown in Figure 3, and the numerical results from the analytic differentiation method and numerical differentiation method compared with each other and the results of commercial packages [8, 13].

Figure 3:

Three-dimensional geometry of single element model under tension (P = 150 N).

The comparison between the results from the analytic differentiation and those from the numerical differentiation is shown in Figures 4, 5, and 6 for 2-coefficient, 5-coefficient, and 6-coefficient Mooney-Rivlin model, respectively. All three models showed a good agreement between the results from the numerical differentiation method and analytic differential method. Plus, the results from the numerical differentiation method showed a good agreement with those from NISA II, also shown in Figure 7.

Figure 4:

Comparison of deflections obtained from numerical differentiation and analytic differentiation (2-coefficient Mooney-Rivlin model: C₁₀ = 0.18299; C₀₁ = 0.00342 MPa; ν = 0.49000).

Figure 5:

Comparison of deflections obtained from numerical differentiation and analytic differentiation (5-coefficient Mooney-Rivlin model: C₁₀ = 0.179; C₀₁ = 0.009; C₂₀ = 0.0018; C₁₁ = − 0.0000846; C₀₂ = − 0.0000371 MPa; ν = 0.49000).

Figure 6:

Comparison of deflections obtained from numerical differentiation and analytic differentiation (6-coefficient Mooney-Rivlin model: C₁₀ = 0.179; C₀₁ = 0.009; C₂₀ = 0.0018; C₁₁ = − 0.0000846; C₀₂ = − 0.0000371; C₃₀ = 0.0000445 MPa; ν = 0.49000).

Figure 7:

Comparison of deflections obtained from numerical differentiation and NISA II (2-coefficient Mooney-Rivlin model: C₁₀ = 0.18299; C₀₁ = 0.00342 MPa; ν = 0.49000).

The comparison between the results from the numerical differentiation method and those from ABAQUS (4-digit data) is shown in Figures 8 and 9 for 2-coefficient and 5-coefficient Mooney-Rivlin model, respectively. The results from the numerical differentiation method and ABAQUS showed a good agreement. All the comparisons are summarized in Table 1.

Table 1:

Comparison of displacements for single element model.

Number of Mooney-Rivlin coefficients	Poisson's ratio	Analytic differentiation	NISA II	ABAQUS	Numerical differentiation	Error (%)

2	0.49000	56.2268	N.A.	N.A.	56.2267	0.00018
2	0.49000	N.A.	56.2270	N.A.	56.2267	0.00053
2	0.49900	N.A.	N.A.	56.21	56.21	0.0
5	0.49900	N.A.	N.A.	46.60	46.60	0.0

Figure 8:

Comparison of deflections obtained from numerical differentiation and ABAQUS (2-coefficient Mooney-Rivlin model: C₁₀ = 0.18299; C₀₁ = 0.00342 MPa; ν = 0.49900).

Figure 9:

Comparison of deflections obtained from numerical differentiation and ABAQUS (5-coefficient Mooney-Rivlin model: C₁₀ = 0.179; C₀₁ = 0.009; C₂₀ = 0.0018; C₁₁ = − 0.0000846; C₀₂ = − 0.0000371 MPa; ν = 0.49900).

The proposed numerical differentiation method was also compared with existing literature [12], and the 3-dimensional geometry and finite element models of rubber under tension used by Peeken are shown in Figure 10. The maximum deflection and strains of the element are compared in Figure 11 and Table 2.

Table 2:

Comparison of results with Peeken and numerical differentiation methods.

Load (MPa)	Peeken		Numerical differentiation		Error (%)
Load (MPa)	Deflection (mm)	Strain	Deflection (mm)	Strain	Error (%)

1.5	400.84	1.20042	400.84	1.20042	0.0

Figure 10:

Three-dimensional geometry and finite element models of rubber under tension (Peeken's model).

Figure 11:

Comparison of strain obtained from numerical differentiation with Peeken's results for four-element model (C₁₀ = 0.370; C₀₁ = 0.015 MPa; ν = 0.48000).

3. Examination of Convergence for Different Poisson Ratios

In this section we will test the convergence of the finite element method for different Poisson ratios using 3-dimensional 20-node elements. First example is simple extension of single element for easy evaluations of volume changes. Second example is the case which is very severe to converge for near-incompressible rubber because of severe distortion of the model under following loads, affected by the applied pressure and inflated shapes.

3.1. Examination for Isoparametric Three-Dimensional 20-Node Rubber Element Model

To examine the convergence for different Poisson ratios, an isoparametric 3-dimensional 20-node rubber element model was examined under loads P₁ and P₂ simultaneously, as shown in Figure 12, and results of the nonlinear analyses under loads P₁ = 300 N and P₂ = 150 N are summarized in Table 3. In this case, the Mooney-Rivlin coefficients are C₁₀ = 0.15 and C₀₁ = 0.05. When the Poisson iteration method was not used, the convergence failed and the incompressibility condition (ΔV/V = 0) could not be approached. However, when using the Poisson iteration method, convergence was achieved, even with Poisson's ratio of 0.49999.

Table 3:

Test of convergence with various Poisson ratios.

Number of Mooney-Rivlin coefficients	Poisson's ratio	Conventional method	Poisson iteration method

2	0.49000	Convergence	Convergence
2	0.49500	Convergence	Convergence
2	0.49900	Convergence	Convergence
2	0.49990	Divergence	Convergence
2	0.49995	Divergence	Convergence
2	0.49999	Divergence	Convergence

Figure 12:

Isoparametric 3-dimensional 20-node rubber element model.

In addition, the volume change (elongation 96.8%) was also examined for various Poisson ratios using the same model with a uniaxial tension: P₁ = 600 N and P₂ = 0 N. The results are shown in Figures 13 and 14, while Table 4 shows comparisons of the volume changes for various Poisson ratios.

Table 4:

Comparison of volume change with various Poisson ratios.

Number of Mooney-Rivlin coefficients	Poisson's ratio (ν)	Deflection (mm)		Convergency of solution	Volume change (%) ((ΔV/V) × 100)
Number of Mooney-Rivlin coefficients	Poisson's ratio (ν)	Analytic differentiation	Numerical differentiation	Convergency of solution	Volume change (%) ((ΔV/V) × 100)

2	0.49000	9.68506	9.68506	Convergence	0.63
2	0.49500	9.68479	9.68479	Convergence	0.31
2	0.49900	9.68456	9.68456	Convergence	0.063
2	0.49950	9.68452	9.68452	Convergence	0.031
2	0.49990	9.68449	9.68449	Convergence	0.0061
2	0.49995	9.68448	9.68448	Convergence	0.0030
2	0.49999	9.68448	9.68448	Convergence	0.00047

Figure 13:

Volume change (elongation 96.8%) with various Poisson ratios under uniaxial tension (P₁ = 600 N and P₂ = 0 N).

Figure 14:

Volume change (elongation 96.8%) with various Poisson ratios under uniaxial tension (detailed view: ν = 0.49950∼0.49999).

3.2. Examination by Numerical Example: 64-Element Rubber Pad Analysis

A rubber pad inflation analysis, as reported in literatures [16, 17], was used as a numerical example for verification. The rubber pad was 200 × 200 mm with a thickness of 12.5 mm, and it was modeled with 64 elements of 3-dimensional 20-node hexahedron. The four sides of the rubber pad were fixed and supported, while pressure was applied to the bottom upwardly. The Mooney-Rivlin coefficients were C₁₀ = 0.55152 MPa and C₀₁ = 0.13788 MPa. Table 5 shows the convergence for different Poisson ratios (ν = 0.4900, 0.4950, 0.4990, and 0.4999) when using the conventional method and proposed Poisson iteration method, while Table 6 shows the displacement results for each pressure when ν = 0.4900 and ν = 0.4999. Since no theoretical results were found for this example, the results were compared with those from commercial packages [8, 13] and existing literature [16], where 450 of 8-node brick elements have been used and compared with the experiment to show a good agreement. Figures 15 and 16 show a comparison of the displacements obtained from the various rubber pad inflation analysis methods; Figure 17 shows the deformed and undeformed models of the rubber pad, analyzed by the proposed techniques using follower force, using I-DEAS as pre- and postprocessor.

Table 5:

Convergency in rubber pad inflation analysis.

Number of Mooney-Rivlin coefficients	Poisson's ratio	Conventional method	Poisson iteration method

2	0.4900	Convergence	Convergence
2	0.4950	Convergence	Convergence
2	0.4990	Divergence	Convergence
2	0.4999	Divergence	Convergence

Table 6:

Comparison of displacement results in rubber pad inflation analysis.

Pressure (MPa)	NISA II (ν = 0.49)	Present (ν = 0.4999)	Difference (%)

0.1	43.5679	38.1005	14.35
0.2	63.7776	60.2794	5.90
0.3	85.0884	83.2739	2.51
0.4	114.348	108.587	6.48

Figure 15:

Comparison of center displacements obtained using various rubber pad inflation analysis methods.

Figure 16:

Comparison of center displacements in rubber pad inflation analysis with different Poisson ratios.

Figure 17:

Deformed and original shapes of rubber pad used in inflation analysis.

4. Approaching Techniques in Poisson Iteration Method

We adopt a pattern technique, expressing Poisson's ratio increase in a certain function, and it was named as a pattern approach.

4.1. Linear Approaching Technique

Poisson's ratio in linear approaching technique increases linearly with approaching 200 steps. The value of it in current step is represented as in (16) and shown in Figure 18:

ν = 0.49 + k \times [\frac{(0.4999 - 0.4900)}{K}],

(15)

where ν is Poisson's ratio, k is number of current steps, and K is number of total steps.

Figure 18:

Poisson's ratio as a function of steps in linear approaching technique.

Using the linear approaching technique, we got a divergence at Poisson's ratio of 0.49801 (number of step: 197) for a given load in the example of Section 3.2, and the number of total substeps was 528.

4.2. Pattern Approaching Technique

Poisson's ratio in pattern approach increases following our proposed pattern that is expressed in (17) and shown in Figure 19. We found it after many trials of pattern functions.

ν = \frac{{1 - [1 / (50 + (4950 \times k / K))]}}{2},

(16)

where ν is Poisson's ratio, k is number of current steps, and K is number of total steps.

Figure 19:

Poisson's ratio as a function of steps in pattern approaching technique.

4.3. Relaxation Factors x

The instability of numerical analysis can be improved effectively using relaxation factor [18] when it is apt to divergence. The instability of rubber analysis can also be reduced with relaxation factors x that are given by analyzer. The ratio is a value that divides present unbalanced load by previous one. When the ratio is lower, then the relaxation factor is x = 1, which means no relaxation. If the ratio is bigger than 1, we adopt some value of x less than 1. The relaxation factors are multiplied to the nonbalanced iterative loads which have been originated from the change of Poisson's ratio approaching to 0.5. The unbalanced load of substep 2 is most critical to the convergency of numerical procedure in other problems like postbuckling [18]; however, whole substeps in Poisson iteration are critical to the convergency. Consider

\begin{matrix} Ratio < 1 : x = 1, \\ Ratio > 1 : x < 1 . \end{matrix}

(17)

Next, we also introduced scale factor S which scales the number of total steps from the initial value of K = 200 based on the number of subiterations as follows:

\begin{array}{l} number of subiteration \leq 4 : sf = S, \\ number of subiteration = 5 : sf = 1, \\ number of subiteration \geq 6 : sf = 1 / S, \\ Number of total steps = current number of total steps \times sf . \end{array}

(18)

When the scale became S = 0.75, the iteration number is shown in Table 7 for a given load in the example of Section 3.2. When the relaxation factors are x = 0.1 and 0.2, the total subiteration number of 235 has resulted. Comparing the optimal relaxation factor of 0∼0.3 in postbuckling analysis [18], the recommended factor of relaxation of 0.1∼0.2 in near-incompressible analysis falls in similar range. Here, we can see that the severe near-incompressible rubber analysis problem better needs the relaxation method with proper factor as in other severe problems in convergence like the postbuckling analysis.

Table 7:

Comparison of the total subiterations for relaxation parameters with S = 0.75.

X	0.1	0.2	0.3	0.4	0.5

Total subiterations	235	235	352	352	Diverged

5. Conclusion

A finite element procedure was derived for hyperelastic materials, such as rubber, using isoparametric 3-dimensional 20-node rubber elements, one of analytic differentiation and numerical differentiation methods, and the Poisson iteration method. In addition, a pattern approaching method in the Poisson iteration technique was proposed to suppress the diverging phenomenon in rubber problem even for the case that Poisson's ratio is very close to 0.5. To resolve this problem, Poisson's ratio of 0.4900 is used first to get an approximate solution without divergence. After convergence, the analysis is restarted and Poisson's ratio is increased to 0.4999 following a proposed pattern and adopting a technique of relaxation and scaling.

The contents of this study can be summarized as follows.

(i) A finite element procedure for analyzing rubber problems was developed using analytic differentiation method and numerical differentiation method to obtain the constitutive matrices. In addition, a Poisson iteration method was proposed to increase the convergence rate when dealing with near-incompressibility.

(ii) When analyzing compressible rubber (ν = 0.4900), the results from the numerical differentiation method matched quite well with those from the analytic differentiation method and with other commercial packages, such as NISA II and ABAQUS. The use of the numerical differentiation method minimized the change in the finite element program when modifying the rubber model (hyperelastic potential function, penalty function, etc.).

(iii) A pattern approaching method in the Poisson iteration technique was proposed. A scale factor for the total number of steps and relaxation factors for effective subiteration were applied to the method. With these factors, the proposed Poisson iteration method was very effective in resolving the divergence problem and obtaining a fast convergence rate in the near-incompressible rubber problem analysis (ν is close to 0.5, i.e., 0.4999) and produced a rather different solution compared to that for compressible rubber (ν = 0.4900). A maximum discrepancy of 14% was displayed between these displacements, indicating that special care should be provided when dealing with most cases of incompressible rubber with mostly confined boundary.

(iv) In respect of the subiteration number, linear approaching method resulted in a divergence after 528 subiterations for a certain case. The present method showed a fast convergence rate and gave the result for Poisson's ratio of 0.4999 in 235 total subiteration steps. The optimum relaxation factor x with the scale factor S = 0.75 which is 0.1∼0.2 is recommended, which is much similar to the other severely nonlinear cases like postbuckling problem, for the fast convergence.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Footnotes

Acknowledgments

This work was supported by the National Research Foundation of Korea (NRF) Grant funded by the Korean Government (MOE) (NRF-2012R1A1A2008903). The authors are partially supported by BK21 plus Project of the National Research Foundation of Korea.

References

Oden

J. T.

and Key

J. E.

, “Numerical analysis of finite axisymmetric deformations of incompressible elastic solids of revolution,” International Journal of Solids and Structures, vol. 6, no. 5, pp. 497–518, 1970.

Scharnhorst

and Pian

T. H. H.

, “Finite element analysis of rubber-like materials by a mixed model,” International Journal for Numerical Methods in Engineering, vol. 12, no. 4, pp. 665–676, 1978.

Jankovich

Leblanc

Durand

, and Bercovier

, “A finite element method for the analysis of rubber parts, experimental and analytical assessment,” Computers and Structures, vol. 14, no. 5–6, pp. 385–391, 1981.

Sussman

and Bathe

K.-J.

, “A finite element formulation for nonlinear incompressible elastic and inelastic analysis,” Computers & Structures, vol. 26, no. 1–2, pp. 357–409, 1987.

Weiss

J. A.

Maker

B. N.

, and Govindjee

, “Finite element implementation of incompressible, transversely isotropic hyperelasticity,” Computer Methods in Applied Mechanics and Engineering, vol. 135, no. 1–2, pp. 107–128, 1996.

Bathe

K.-J.

, Finite Element Procedures, Prentice Hall, Upper Saddle River, NJ, USA, 1996.

Borst

Crisfield

M. A.

Remmers

J. J. C.

, and Verhoosel

C. V.

, Non-Linear Finite Element Analysis of Solids and Structures, John Wiley & Sons, London, UK, 2nd edition, 2012.

ABAQUS User's Manual, version 6.13, 2014.

Treloar

L. R. G.

, “Stress-strain data for vulcanised rubber under various types of deformation,” Transactions of the Faraday Society, vol. 40, pp. 59–70, 1944.

10.

Haines

D. W.

and Wilson

W. D.

, “Strain-energy density function for rubberlike materials,” Journal of the Mechanics and Physics of Solids, vol. 27, no. 4, pp. 345–360, 1979.

11.

Attard

M. M.

and Hunt

G. W.

, “Hyperelastic constitutive modeling under finite strain,” International Journal of Solids and Structures, vol. 41, no. 18–19, pp. 5327–5350, 2004.

12.

Peeken

Döpper

, and Orschall

, “A 3-D rubber material model verified in a user-supplied subroutine,” Computers and Structures, vol. 26, no. 1–2, pp. 181–189, 1987.

13.

“NISA II User's Manual,” EMRC, version 18.0, 2005.

14.

Kwon

Y.-D.

Rho

K.-T.

, and Won

D.-W.

, “A study on the convergency of the finite element analysis of rubber using numerical differentiation method,” in Proceedings of the AMDP Conference vol. 1, pp. 353–358, 1999.

15.

Arora

J. S.

, Introduction to Optimum Design, McGraw-Hill, New York, NY, USA, 1989.

16.

Shin

S.-J.

Lee

T.-S.

, and Oh

S.-Y.

, “Three-dimensional finite element analysis of rubber pad deformation—part II: applications,” Transactions of Korean Society of Mechanical Engineers A, vol. 22, no. 1, pp. 121–131, 1998 (Korean).

17.

Kwon

Y.-D.

Kwon

H.-W.

Kim

W.-J.

, and Yeo

S.-D.

, “Estimation of rubber material property by successive zooming genetic algorithm,” Journal of Solid Mechanics and Materials Engineering, vol. 1, no. 6, pp. 815–826, 2007.

18.

Kwon

, and Lim

, “Nonlinear analysis using load-displacement control,” Structural Engineering and Mechanics, vol. 19, no. 2, pp. 153–172, 2005.

A Finite Element Procedure with Poisson Iteration Method Adopting Pattern Approach Technique for Near-Incompressible Rubber Problems

Abstract

1. Introduction

2. Finite Element Analysis of Rubber

2.1. Strain Energy Function

2.2. Numerical Differentiation

2.3. Poisson Iteration Method

2.4. Verification of Proposed Method

3. Examination of Convergence for Different Poisson Ratios

3.1. Examination for Isoparametric Three-Dimensional 20-Node Rubber Element Model

3.2. Examination by Numerical Example: 64-Element Rubber Pad Analysis

4. Approaching Techniques in Poisson Iteration Method

4.1. Linear Approaching Technique

4.2. Pattern Approaching Technique

4.3. Relaxation Factors x

5. Conclusion

Conflict of Interests

Footnotes

Acknowledgments

References