• Nem Talált Eredményt

Finite Element Analysis of Rubber Bumper Used in Air-springs

N/A
N/A
Protected

Academic year: 2022

Ossza meg "Finite Element Analysis of Rubber Bumper Used in Air-springs"

Copied!
8
0
0

Teljes szövegt

(1)

Procedia Engineering 48 ( 2012 ) 388 – 395

1877-7058 © 2012 Published by Elsevier Ltd.Selection and/or peer-review under responsibility of the Branch Offi ce of Slovak Metallurgical Society at Faculty of Metallurgy and Faculty of Mechanical Engineering, Technical University of Košice

doi: 10.1016/j.proeng.2012.09.530

MMaMS 2012

Finite Element Analysis of Rubber Bumper Used in Air-springs

Tamas Mankovits

a

*, Tamas Szabó

b

aDepartment of Mechanical Engineering, University of Debrecen, Debrecen 4028, Hungary

bDepartment of Robert Bosch Mechatronics, University of Miskolc, Miskolc 3515, Hungary

Abstract

In this paper a solution of a contact problem for large displacements and deformations are analyzed where a rubber bumper applied in air- spring. The nonlinear load-displacement curve is determined. A FEM code written in FORTRAN has been developed for the analysis of nearly incompressible axially symmetric rubber parts. The Hu Washizu type variational principle is formulated for the Mooney-Rivlin material model. Stability and sensitivity analyses are also investigated.

© 2012 The Authors. Published by Elsevier Ltd.

Selection and/or peer-review under responsibility of the Branch Office of Slovak Metallurgical Society at Faculty of Metallurgy and Faculty of Mechanical Engineering, Technical University of Košice.

Keywords: finite element method; rubber; large displacements; contact

1. Introduction

Numbers of papers are devoted to the application of h-version finite element method for the analysis of hyperelastic materials [1-6]. The applications of p-version of the finite element for geometrically and physically nonlinear problems are relatively recent [7-10]. It is well known that the elastomers are regarded as incompressible or nearly incompressible materials. The incompressibility is an auxiliary condition which can be enforced by the so called mixed method. In this paper the traditional penalty method and the mixed formulation of three fields are implemented. The later one is based on the Hu Washizu type variational principle [11]. The displacements and the volumetric change are approximated independently from the hydrostatic pressure. The same order of approximation is selected for the hydrostatic pressure and the volumetric change. The displacement fields are approximated by the tensor product space and its order is higher than the order of space used for the pressure and the volumetric change at each level of polynomial degree p.

In engineering practice high stresses occur in unilateral contact. The contact problem is handled with a simplified penalty approach. The contacting boundary is approximated by a polygon. The edge of the contacting element is forced to have a straight line. In practice, two-node contact elements were implemented.

As a numerical example a rubber bumper applied in air-spring is chosen for this problem. The load-displacement curve for pressure is given by the producer but only when the bumper is loaded between two flat rigid surfaces. The aim of this research is to obtain the load deflection curve up to 20% of compression in realistic working circumstances.

* Corresponding author. Tel.: +36-52-415-155;

E-mail address: tamas.mankovits@eng.unideb.hu

© 2012 Published by Elsevier Ltd.Selection and/or peer-review under responsibility of the Branch Offi ce of Slovak Metallurgical Society at Faculty of Metallurgy and Faculty of Mechanical Engineering, Technical University of Košice

(2)

2. Large displacements and deformations

Using the Lagrangian description in the reference cylindrical coordinate system the motion of continua

r & = r & ( r &

0

, t )

,

in the current configuration at time

t

is given by the reference coordinates

r &

0

of the material points, see in Fig. 1., where

u &

is the displacement.

u0=u

r0 r

r z

t=0

t=t P0(r0,z00)

P(r,z,ϕ)

V0 V

dV0

dV

Fig. 1. The motion of continua

In the subsequent formulation the deformation gradient

r

0

F & r

&

= ∂

(1)

is multiplicatively decomposed into a volumetric changing part

F

V and a volume preserving part

F ˆ F

F

F

V

= ˆ

(2)

where

I J

F

3

1

ˆ =

,

det F ˆ = 1

(3)

and

F J F

V

3 1

=

,

F F J

V

= det =

det

(4)

where

J

is the Jacobian. Using this decomposition we can define the so called right Cauchy-Green tensor and its volumetric preserving part

F F

C =

T ,

C ˆ F ˆ

T

F ˆ J

3

C ˆ

2

=

=

(5)

where Tdenotes the transpose of a tensor. The deformation is expressed by the Green-Lagrange strain tensor

) 2 (

1 C I

E = −

(6)

(3)

where

I

is the unit tensor. The strain energy functions are defined by this specific decomposition

ˆ ) ˆ ( ) ( ˆ ) ,

( J C U J W C

W = +

, (7)

where

U (J )

denotes the strain energy due to volumetric change and

W ˆ ( C ˆ )

is the strain energy determined by volume preserving deformation. The function

U (J )

is often defined in the papers in the simpliest way

)

2

1 2 (

) 1

( J = J

U κ

, (8)

where

κ

is the bulk modulus. The function

W ˆ ( C ˆ )

is given for the Mooney-Rivlin material

) ˆ 3 ( ) ˆ 3 ( ˆ )

ˆ (

01

10

− + −

= I

I

I

II

C

W μ μ

(9)

where

μ

10 and

μ

01 are the material constants,

I ˆ

I and

I ˆ

II are the first and second invariants of the unimodular right Cauchy-Green tensor defined as

33 22

11

ˆ ˆ

ˆ

ˆ C C C

I

I

= + +

,

( ˆ ˆ ˆ )

2

ˆ 1 I C C

I

II

=

I

− ⋅⋅

(10)

In the reference coordinate system, the II. Piola-Kirchhoff stress tensor for the rubber is given by

)

1

( ˆ ˆ

2 +

= ∂ p J C C

C

S W

(11)

where

p

is the hydrostatic pressure. In the current configuration the Cauchy stress tensor can be calculated as

F

T

S F J

T =

1 . (12)

3. Contact kinematics

Let us consider a system which consists of two bodies. One is elastic signed by 1, the other is a rigid body signed by

2see in Fig. 2., where

A

c is the contact region.

V

1

V

2

A

c1

A

c2

n

c

h

u

n1

Q

1

Q

2

Fig. 2. Kinematics in the contact region

(4)

In this work normal contact is assumed, where

n &

c is the normal unit vector of the contact surfaces through two point pairs

Q

1 and

Q

2. The

g

n normal gap can be defined as

h u u g

g

n

=

n

( & ) = −

n1

+

,

A

c

r & ∈

, (13)

where

h

is the initial gap and

u

1n is the normal component of the displacement vector. There can be two cases, one is contact, when

0

0 ≥

=

n

n

p

g

, (14)

there is gap (no contact) between the two bodies if

0

0 =

n

n

p

g

, (15)

where

p

n is the contact pressure. Both cases

= 0

n n

g

p

,

r & ∈ A

c. (16)

4. Variational principles and FEM discretization

The mixed formulation based on the Hu Washizu type variational principle [11], the three fields are independent from each other,

) 2 (

) 1 ( )

( ˆ )

ˆ ( ) , ,

( u J p W C dV U J dV p J J dV g

2

dA

ext

u

A n

V V

hw V

c

&

& = + + − + − Π

Π ³ ³ ³ ³ γ

, (17)

where

J

is the volumetric change,

J

is evaluated by equation (4) making use of the displacements

u &

. The

γ

is the

penalty parameter which can be mentioned as a virtual spring stiffness called as a Winkler-type contact problem.

Since equation (17) is a nonlinear one we should linearize it introducing the increments of the appropriate variables.

Before the FEM discretization we introduce the Green-Lagrange strain tensor (6) and the energy conjugated II. Piola- Kirchhoff stress tensor (11). The increment in the Green-Lagrange strain tensor

Δ E

is decoupled into linear and nonlinear parts

NL

L

E

E E = Δ + Δ

Δ

. (18)

Meridian section is discretized by two-dimensional rectangular elements see in Fig. 3.

ξ

η

1

4 3

2 1

4

3

(1,-1) 2

(1,1) (-1,1)

(-1,-1)

z

r 5

9 8 8

7 7

6 6

5 9

η

Fig. 3. Nine-node two-dimensional isoparametric element

(5)

In case of mixed formulation on each element the displacement fields are approximated by the so called tensor product space using Legendre polynomials [12]

i n i

N

i

u

u & = ¦

=1

&

,

n = 4 + 4 ( p − 1 ) + ( p − 1 )

2, (19)

where

N

i is the shape function,

u &

i

are the displacement parameters,

n

is the number of shape functions which non- linearly depends on the level of approximation

p

. The volumetric change and the hydrostatic pressure are approximated by lower order of polynomials than the displacement

z a r a a p

J ( = 2 ) =

0

+

1

+

2 , (20)

z b r b b p

p ( = 2 ) =

0

+

1

+

2 . (21)

After the discretization we obtain the formula for the Newton-Raphson iteration

f f

f u

K

T

Δ =

ext

int

= Δ

, (22)

where

K

T is the tangential stiffness matrix,

Δ u

is the displacement increment and

Δ f

is the unbalanced load vector. A p-versional finite element code has been developed for the analysis of nearly incompressible materials. The approximation order for the displacement is

p = 2

.

5. Numerical example

A rubber bumper in an air-spring (see in Fig. 4) is analyzed by the FEM code based on the theory discussed above.

Fig. 4. The air-spring

The load-displacement curve for pressure is known from the producer but only when the bumper is loaded between two flat rigid surfaces. A parameter optimization is needed to find the material parameters. The measurement and the code were compared, see in Fig. 5.

(6)

Fig. 5. The load-displacement curve

A numerical analysis is also needed to obtain the optimal input parameters for the finite element analysis such as mesh density and the so called penalty parameter. The effect of the mesh density change was compared at two discrete points (7 and 14 mm compression). The results of the mesh density analysis can be seen in Fig. 6. The 8x8 mesh gives good correlation to the more fined mesh and its time demand is much more cost efficient. The mesh density versus the volume ratio is evaluated, see in Fig. 6.

(a) (b)

Fig. 6. Numerical analyses for (a) mesh density and (b) volume ratio

In examination of a rubber part it is very important to observe how the incompressibility condition is fulfilled which depends on the applied penalty parameter. In Fig. 7. the numerical stability analysis calculated with different penalty parameters at mesh density 8x8 (acceptance limit is 99,9%) and the displacement field calculated by the FEM code can be seen.

(a) (b)

Fig. 7. Results of (a) penalty parameter analysis and (b) the code

(7)

After the numerical stability analysis the code parameters are determined from the results. In Table 1 the most important finite element settings can be seen.

Table 1. Finite element settings

Material parameters

Mooney-Rivlin parameter (

μ

10) 0,63N/mm2

Mooney-Rivlin parameter (

μ

01) 0,1575N/mm2

Finite element code

Penalty parameter (

κ

) 1000

Increment (

Δ u

) 1mm

Supports (top and below) vulcanized

The rubber bumper may contact with two contact regions in working circumstances. The mechanical model of the contact problem and the results of the code for displacement can be seen in Fig. 8a and Fig 8b, respectively.

(a)

z

r RUBBER

STEEL

(b)

Fig. 8. Rubber bumper’s (a) mechanical model and (b) displacement at 13mm compression

The load-displacement curve calculated by the FEM code can be seen in Fig. 9.

Fig. 9. Load-displacement curve of the rubber bumper in working circumstances

(8)

6. Summary

A solution of a contact problem for large displacements and deformations are analyzed where a rubber bumper applied in air-spring is investigated by determining the nonlinear curve of load versus displacements. A FEM code which is able to handle contact written in FORTRAN has been developed for the analysis of nearly incompressible axially symmetric rubber parts. Later a shape optimization problem will be done if the stability analysis and also the comparison with commercial software are adequate.

Acknowledgements

The described work was carried out as part of the TÁMOP-4.2.2/B-10/1-2010-0008 project in the framework of the New Hungarian Development Plan. The realization of this project is supported by the European Union, co-financed by the European Social Fund.

References

[1] Malkus, D.S., 1980. Finite Elements with Penalties in Nonlinear Elasticity, International Journal for Numerical Methods in Engineering 16, p. 121-136.

[2] Swanson, S.R., Christensen, L.W., Ensing, M., 1985. Large Deformation Finite Element Calculations for Slightly Compressible Hyperelastic Materials, Computers & Structures 21, No. 1/2, p. 81-88.

[3] Zienkiewicz, O.C., Qu, S., Taylor, R.L., Nakazawa, S., 1986. The Patch Test for Mixed Formulations, International Journal for Numerical Methods in Engineering 23, p. 1873-1883.

[4] Sussman, T., Bathe, K.J., 1987. A Finite Element Formulation for Nonlinear Incompressible Elastic and Inelastic Analysis, Computers & Structures 26, No. 1/2, p. 357-409.

[5] Simo, J.C., Taylor, R.L., 1991. Quasi-incompressible Finite Elasticity in Principal Streches. Continuum Basis and Numerical Algorithms, Computational Methods in Applied Mechanics and Engineering 85, p. 273-310.

[6] Gadala, M.S., 1992. Alternative Methods for the Solution of Hyperelastic Problems with Incompressibility, Computers & Structures 42, No. 1, p. 1-10.

[7] Düster, A., Hartmann, S., Rank, E., 1991. p-fem applied to Finite Isotropic Hyperelastic Bodies, Computational Methods in Applied Mechanics and Engineering 85, p. 5147-5166.

[8] Nándori, F., Páczelt, I., Szabó, T., 2003. „Analysis of Incompressible Materials with p-version Finite Elements,” microCAD 2003 International Scientific Conference, p. 1-6.

[9] Hartmann, S., Neff, P., 2003. Polyconvexity of Generalized Polynomial-type Hyperelastic Strain Energy Functions for Near-incompressibility, International Journal of Solids and Structures40., p. 2767-2791.

[10] Szabó, Mankovits, T., 2004. „FEM Computations of Hyperelastic Materials,” microCAD 2004 International Scientific Conference, p. 79-84.

[11] Bonet, J., Wood, R.J., 1997. Nonlinear Continuum Mechanics for Finite Element Analysis, Cambridge University Press.

[12] Szabó, B.A., Babuska, I., Chayapathy, B.K., 1989. Stress Computation for Nearly Incompressible Materials by the p-version of the Finite Element Method, International Journal for Numerical Methods in Engineering 28., p. 2175-2190.

Ábra

Fig. 1. The motion of continua
Fig. 2. Kinematics in the contact region
Fig. 3. Nine-node two-dimensional isoparametric element
Fig. 4. The air-spring
+3

Hivatkozások

KAPCSOLÓDÓ DOKUMENTUMOK

In this pape, two specific applications of the Hopfield neural network will be discussed: First, for obtaining the solution of finite element analysis directly by

The design process of the intake manifold system will consist of the usual engineering processes including computer modelling, Finite Element Analysis and finally Computational

In this paper, an e ffi cient method is developed for the formation of null bases for finite element models comprising of rectangular plane stress and plane strain serendipity

The method of finite elements eliminates this difficulty by considering the surface (for simplicity's sake, continua 'will not he treated below) to be diyided into

In order to maintain displacement com- patibility between the beam and the sti ff ened element, a special transformation is used, which includes the coupling of torsional and

Temperature distribution in the wheel 0.257 s after the end of heat input, the maximum surface temperature is 4.45·10 −2 °C (initial temperature: 0

When a finite element stress analysis is used, the three principal stresses and their direc- tions must be calculated and the volume of corresponding elements must

In order to analyze the problems arisen in application of the Preisach model, a two- dimensional field problem is developed for hysteresis motor with determi- nation