Temperature-dependent modelling of a HNBR O-ring seal above and below the glass transition temperature

Materials and Design 156 (2018) 1–15

Contents lists available at ScienceDirect

Materials and Design journal homepage: www.elsevier.com/locate/matdes

Temperature-dependent modelling of a HNBR O-ring seal above and below the glass transition temperature J. Troufflard a , H. Laurent a, * , G. Rio a , B. Omnès b , S. Javanaud b a b

Univ. Bretagne Sud, UMR CNRS 6027, IRDL, Lorient F-56100, France Cetim, 74 route de la Jonelière, Nantes CS850814, France

H I G H L I G H T S

G R A P H I C A L

A B S T R A C T

• Identification of the parameters of a hyper-visco-hysteresis model simulating an O-ring test for low temperature cycles. • Taken the loss of contact force occurring at low temperatures into account. • Decomposition into volumetric/ distortional and viscous/non-viscous parts for developing the model. • The hyperelastic part can be identified directly by performing anisothermal stress relaxation tests.

A R T I C L E

I N F O

Article history: Received 18 April 2018 Received in revised form 7 June 2018 Accepted 8 June 2018 Available online 15 June 2018 Keywords: Elastomer Sealing Environmental performance Low temperature Hyperelasto-visco-hysteresis model

A B S T R A C T Sealing products can undergo many different atmospheric conditions and hence large variations in the in-service temperatures. These products, typically made of elastomer, are known to have heat sensitive mechanical characteristics. Since the reliability of sealing systems can be affected by cyclic temperature conditions, it is essential to predict seal performances depending on the temperature in order to prevent the occurrence of interfacial leaks. The aim of the present study is to simulate the thermo-mechanical behavior of a HNBR O-ring seal during a temperature cycle ranging from room temperature to low temperatures down to −34 ◦ C. The main objective is to present the construction and the thermal dependency of a hyperelasto-visco-hysteresis (HVH) model so as to be able to predict the mechanical behavior of the elastomer in function of temperature. This phenomenological model is designed in the form of a combination of hyperelastic, viscoelastic and elasto-hysteretic stress contributions. The relevant parameters are identified using a material database based on the results of classical homogeneous isothermal relaxation tests and anisothermal relaxation tests. The predictions of the HVH model obtained shows good agreement with the experimental contact force during a 20% squeezed O-ring test under cyclic thermal loading below the glass transition temperature. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction

* Corresponding author. E-mail address: [email protected] (H. Laurent).

https://doi.org/10.1016/j.matdes.2018.06.016 0264-1275/© 2018 Elsevier Ltd. All rights reserved.

Most of the elastomeric seals designed for static or pseudo-static high-pressure applications in the oil and gas industry are of the “squeeze” type. The most common of these seals is the O-ring, which

2

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

is set in a rectangular housing and “squeezed” by initial compression to give its initial sealing capability. This initial squeeze is then increased by applying system pressure to deform the seal. The overall sealing force thus created results from the combination between this initial sealing force and the system pressure applied [1]. In many technical applications, it is essential for the tightness of this kind of systems to be able to resist a large range of temperatures. Ensuring reliability of systems exposed to low temperature conditions is therefore an important industrial issue. In order to determine the low temperature limits of sealing applications, several methods have been frequently used to define the efficiency of the tightness [2] and to establish suitable standards [1,3]. The temperature retraction procedure, the Gehmann test and compression set measurements have been classically used for this purpose, for example [4-6]. The tightness of elastomer O-ring seals at low temperatures has also been frequently investigated by measuring the gas leakage rate and the sealing force during specific thermal cycles [7]. In this case, leakages can occur at low temperatures because the contact pressure decreases due to the thermal contraction and a loss of elasticity of the elastomer occurs during the cooling phase. Homogeneous Dynamic Mechanical Analysis (DMA) tests provide another means of measuring the sealing characteristics of material at low temperatures [8,9]. Parallel to these experimental campaigns, the numerical modelling of sealing systems at low temperatures is still a challenge because of the complexity of the heat-dependent structural changes involved [10-12]. High temperatures can result in leakage because of the decomposition of the material, whereas at low temperatures, the rubber-glass transition (abbreviation: glass transition, noted Tg ) results in even more significant parametric changes which affect the efficiency of the seal [13]. The transition to the glassy state is accompanied by the complete loss of rubber elasticity [14-16]. This change is due to a conspicuous switch in the behavior of the structure from entropic elastic to strongly viscous behavior below the glass transition temperature. This long-term viscosity induces an irreversible change in what can be called the quasi elasto-plastic behavior of the material. The springback and recovery processes become very slight near the glass transition temperature Tg . The stiffness of the elastomer increases as the temperature decreases, and the glass transition amplifies this trend. However, the reversibility of the glass transition process is another important characteristic which has to be taken into account. Although a leakage can occur at low temperatures, all the performances of the seal can be recovered if a sufficiently high temperature is restored. This transition from a highly reversible material towards a non-reversible one has to be modelled to ensure accurate finite element analysis of sealing products’ reliability. The main problem here is how to model the mechanical behavior of a material in which the properties change considerably depending on the temperature. However, the classical hyperelastic and viscoelastic models can be still used for the numerical modelling of O-ring seals as a function of the temperature [17,18]. Li et al. [19] have presented a novel phase-transition-based viscoelastic model including the time factor above and below the glass transformation temperatures separately, based on multiplicative thermoviscoelasticity principles. Martinez et al. [11] have used a statistical approach to model a wide range of frequency and strain amplitudes as a function of the temperature. Akulichev et al. [20], after a short review of the findings in the literature on elastomeric seals at low and high pressures and depending on the temperature, have presented a model in which the longterm hyperelastic and viscoelastic (time-dependent) responses and the finite compressibility of the materials are taken into account. The latter authors have established that seal failures are mainly due to detachment of elastomer seals from their mated sealing parts because of thermal contraction processes depending on the adhesion bond.

The aim of this study is to present a new phenomenological approach to simulate the finite-strain relaxation behavior of a hydrogenated nitrile butadiene rubber (HNBR) seal under ambient and sub-ambient conditions below the Tg . In a previous study [21], the experimentally measured loss of contact of a HNBR O-ring at low temperatures was described. In this paper, the authors reported that a one order Ogden material model combined with Prony series and time-temperature equivalence method did not efficiently simulate the loss of contact which occurs during the cooling phase. Using a hyperelasto-visco-hysteresis (HVH) model, previously adopted in the context of a FKM elastomer at room temperature and at temperatures above Tg [22], it is therefore proposed to study the response of this model under the experimental conditions described in [14] and [21], in which the rubber is under constant strain throughout the glass transition phase. The general framework of the HVH model consists of several contributions, each of which is defined as a power related to the total macroscopic strain. This is a parallel rheological model which involves an additive decomposition of the stress. From the practical point of view, the experimental stress can be analyzed in terms of these contributions in the case of a uniaxial test. The total strain is used as a macroscopic experimental datum without making any assumptions about its decomposition. The temperature dependent parameters of the model is used in order to transform the nature of the elasticity from hyperelasticity to elasto-plasticity near Tg . In addition, the fact that the viscosity of the material increases near Tg has to be taken into account when modelling the thermal evolution of the corresponding parameters. It is also important to mention that the ageing of the material in high temperature, which leads to residual deformation in compression is not taken into account in this study. This residual deformation corresponds to the “compression set” (CS) and this parameter highlights the sealing capability of the seal. The paper is composed as follows. In Section 2, the HVH model is described in terms of the constitutive equations for the mechanical partitioning contributions and the type of decomposition choices used on the volumetric and shape contributions. In Section 3, the experimental database presented in [14] is described and used to meet the parameter identification requirements. The identification procedure and the material parameters of the HVH model are described in Section 4. In Section 5, it is reported how the validity of the model was established in experimental tests such as O-ring relaxation tests under cyclic temperature conditions. The study is summarized and the results are discussed in Section 6 along with some suggestions for future studies.

2. Material model The HVH model, which is based on a phenomenological approach, consists of several stress contributions, each of which is based on a mechanical power and defines a different qualitative mechanical aspect of the overall behavior of the material. For further details about the origin of this model and its compatibility with the thermodynamic dissipation principle, see [23-25]. This 3D model is implemented in the in-house finite element software program Herezh++ [26]. This program is interfaced with Abaqus [27] but in the present study, all the simulations were performed using only Herezh++. Prior to this phenomenological decomposition, a more general mechanical assumption was adopted in order to split the deformation locally into a volumetric part and a shape (isochoric) part, as originally proposed by Flory [28] and subsequently done in the context of isothermal finite strain elasticity (e.g. [29] and [30]). The volumetric changes depend on the behavior of the material under purely spherical loading, as in oedometric tests, whereas the changes in shape depend on the distortional strain. This partitioning into two

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

parts is classical in the case of the simple shear tests, which involves distortion of the material and no change in the volume. It requires an appropriated choice of strain measure and invariants in order to prevent the occurrence of undesirable spherical-deviatoric coupling. In the following sections, the basic stress contributions are first listed, giving some details about how they will be used. A more general point of view is then presented, in which the advantages of the volume-shape partitioning procedure are described and how it improves the modelling of the material behavior. The last subsection describes the assumptions and arbitrary choices made in order to simplify the parametric identification procedure. The aim of this study is to design a very flexible and powerful model and to provide suitable coupling between the parameters to fit the O-ring thermal cycle test, using a non-exhaustive database of classical mechanical tests. 2.1. Constitutive equations 2.1.1. General framework The general framework of the HVH model results from a balance between the mechanical powers, where the internal power is decomposed as follows: dPext = dPint = s : e˙ dv = (se + sv + sh ) : Ddv

(1)

where dPext is the external power in the case of an elementary volume dv, dPint is the internal power, s is the Cauchy stress tensor, e is the Almansi total strain tensor and D is the strain rate tensor. e˙ denotes the temporal derivative of the strain tensor in an objective frame. It has been classically established using the Lie two times covariant derivative of the Almansi strain tensor that e˙ ij = Dij (see for example [29]) where eij are the components in the dual frame g i of  /∂ hi where hi is a set of material coordithe material frame gi = ∂ OM nates of the position point M. The total stress tensor s is divided into three parts: a fully reversible hyperelasticity part se , a viscous part sv and the hysteresis part sh , which introduces a non-reversibility in the model. These contributions are included in the HVH model by means of an additive law. The parallel rheological model is given by Eq. (1) and the total stress is equal to ¯ e + sv + sh s = Se + s

(2)

¯ e denote the spherical and deviatoric parts of the where Se and s hyperelasticity, respectively. The following sections give further details about the partitioning of the mechanical behavior and the basic constitutive equations for the various contributions. A more comprehensive description of the HVH model can be found, for example in the case of elastomeric material, in [22]. Before presenting the model in detail, a general hyperelastic framework has to be defined. All the hyperelastic potentials can be defined relatively to the initial or deformed configuration. The following general equation gives the hyperelastic stress tensor obtained from any hyperelastic potential: √ √ 1 ∂W 1 ∂ (w g) g0 ∂ W se = √ = √ = g g ∂e V ∂e ∂e

(3)

where w and W are the energy density in the current and initial √ √ configurations, respectively. g and g0 are the roots of the determinants of the Jacobian matrix in the current and initial states, respectively. V is the volumetric change, which is equal to the ratio √ √ between the current volume v = g and the initial volume v0 = g0 v  according to V = v . In the current state, for instance g = det(gi .gj ) 0

3

√ and dv = gdhi . In the context of the finite element method, the material parameter hi selected were simply the coordinates of any point M in the reference finite element. It is worth noting that the strain tensor e is the Almansi measure in Eqs. (1) and (3) because the temporal derivative of the Almansi strain tensor components applies to the current configuration. This derivation therefore yields the Cauchy stress tensor s. 2.1.2. Spherical hyperelasticity The spherical part of the total strain tensor is related to the volumetric changes. The compressibility is governed by a hyperelastic potential defined as follows: we =

K 2 ln (V) 2

(4)

where K is the bulk modulus. Based on Eq. (3), the hydrostatic pressure can be expressed in the more clearly understandable form: −p = V

∂ we + we ∂V

(5)

where −p = 13 trace(s) is the hydrostatic pressure. This potential gives the spherical part Se of the hyperelastic contribution in Eq. (2). More specifically, it defines the whole spherical part of the total stress tensor s because the rest of the HVH model is purely deviatoric, as explained below. v−v It is worth noting that the relative volumetric change v 0 is close 0 to 0 in the case of a quasi-incompressible material. Given the potenv−v0 tial (4), Eq. (5) can be rewritten with the variable V = v + 1, and 0 developed as a Taylor series in the region of 0. Neglecting the second order terms, the spherical part selected can be related to the physical meaning of the bulk modulus K, defined as follows: −p ≈ K

v − v0 v0

(6)

This comment should help to understand the mathematical expression for the spherical potential. 2.1.3. Deviatoric hyperelasticity The Hart-Smith potential [31,32] is used to define the entropic ¯ e is concerned. elasticity of the elastomer. Only the deviatoric part s This potential is written using the original decomposition of the strain energy function We proposed by Hart-Smith in [31] as follows: ⎧ ∂ We ⎪ ⎪ ⎨ ∂ J1 ∂ We ⎪ ⎪ ⎩ ∂ J2

  = C1 exp C3 (J1 − 3)2 =3

C2 J2

(7)

where J1 , J2 are the first and second modified invariants of the left Cauchy-Green tensor B and C1 , C2 , C3 are material coefficients. By definition, the invariants of B are sensitive to the volumetric changes. To remove this volume variation, the volumetric variation independent invariants J1 , J2 have been classically used with potentials in the framework of elastomer models such as the Mooney-Rivlin and Hart-Smith models. Their definitions are −1/3

J1 = I1 • I3

−2/3

J2 = I2 • I3

(8)

(9)

4

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

where I1 , I2 , I3 are the invariants of B expressed as ⎧ IB = I:B = trace(B) ⎨ I1 = I2 = I(B • B) = I : (B • B) = trace(B • B) ⎩ I3 = I(B • B • B) = I : (B • B • B) = trace(B • B • B)

(10)

and I is the identity tensor. Along with the spherical potential described above, the Hart¯ e requires in order to fully Smith potential gives the deviatoric part s define the hyperelastic contribution se in Eq. (1). 2.1.4. Deviatoric viscosity The generalized Maxwell model is used to introduce time dependency into the HVH model. This rheological model consists of a parallel assembly of several branches composed of an elastic spring and a viscous dashpot in series. A purely spherical test on the elastomer is assumed to involve negligible viscosity, and no viscous hysteresis will therefore be observed in oedometric tests. The constitutive equation corresponding to the generalized Maxwell model for any branch i is written as follows: ˙¯ v + 1 s ¯ = 1 s ¯v D 2Gi gi

(11)

¯ and s ¯ v are the deviator tensors of D and sv , respectively. where D ˙¯ v is the Jaumann time derivative of the deviatoric part of the stress s ¯ v . Gi is the shear modulus of the linear elastic spring, which tensor s is classically defined in terms of the Young’s modulus Ei and the Poisson ratio mi . gi is the material parameter associated with the dashpot which are linked by the characteristic time ti defined by gi = 2Gi ti

(12)

In the subsequent application to the HNBR O-ring, five Maxwell branches were selected. In Eq. (11), only the deviatoric parts of the strain rate and stress tensors are used. The viscous contribution sv in Eq. (1) therefore does not have a spherical part and does not affect the compressibility of the overall model. In addition, the value of the Poisson ratio has no effect because deviatoric behavior combined with an exact decomposition into volume and distortional deformations results in a perfectly incompressible viscous part. The underlying value of mi could be 0.5. Only the parameters ti and Gi still have to be identified. 2.1.5. Elasto-hysteresis model The elasto-hysteresis model [33] is not a classical and it is less well known. This model is the integral form of Saint-Venant’s rheological model, which consists of an infinite parallel assembly of spring and slider branches [34]. It can be used for modelling metal plasticity and shape memory alloys [35], but is also suitable for use with polymers and elastomers [36]. In particular, it involves discretely stored loading inversions which can be used to capture the Mullins effect under cyclic loading conditions. The incremental constitutive equation in the elasto-hysteresis model is defined as follows: ˙¯ h = 2l • D ¯ + b • 0 • Dt s ¯ s R h

(13)

The two following relations give the expressions for b and 0: b=

−2l (w Q0 )np (QDs¯ h )2−np

¯ − ¯h :D 0 = DtR s

(QDs¯ h )2 w ˙ • 2l w

(14)

(15)

¯ h and 0 is a non-reversible where b is related to the intensity of DtR s ¯ h describes the evolution of the deviintrinsically dissipated rate. DtR s atoric part of the hysteretic stresstensor between a reference state r ¯ h : DtR s ¯ h . The parameter w is and current state t and QDs¯ h = DtR s the Masing similarity function, which was taken here to be equal to 2. l and Q0 are the initial slope and the yield limit of the hysteresis stress-strain curve, respectively, under shear loading conditions. np is the Prager’s parameter which makes the non-linear transition for reach the stress saturation Q0 . The elasto-hysteresis model is a useful tool for modelling a metalplasticity like behavior with only three parameters and lends itself to the simulation of loading cycles. During these cycles, the inversion and crossing points are accounted for using the intrinsic dissipation rate function 0 defined above. This value is referred to a volume element and must always be positive. The state at time t is therefore an inversion point which occurs when the function 0 becomes negative. But the present study focus on the ability of the elasto-hysteresis model to introduce non-reversibility and stiffness at low temperatures. The irreversibility, which is given by sh , can be said to be an approximation of the very long term viscosity. 2.2. Mechanical partitioning The HVH model is based on a two levels partitioning procedure. On the one hand, the volumetric part is decoupled from the isochoric part of the transformation. This transformation is described in terms of the strain tensor. This mathematical object must be well chosen and is strongly linked to the time-integration of the strain rate tensor. On the other hand, there is a phenomenological partitioning, which is clearly expressed by Eq. (1) and depends on how the mechanical origin of the stress is interpreted. These two parts have a practical purpose since they determine the type of mechanical tests to be performed and the experimental data used in the identification process. The transition from an initial volume unit to a deformed state can be regarded as a combination between a purely volumetric variation and an isochore distortion corresponding to an angular change. Fig. 1 gives a schematic representation of this transformation. This geometrical point of view can be expressed in terms of invariants of the strain tensor with an appropriated choice of strain measure. The logarithmic strain measure [37] accounts naturally for the decoupling between volume variation and shape changes because its first invariant is directly linked to V, the relative volume change. It has by now been clearly established in the literature, as in [37-39], for example, that numerical integration of the strain rate tensor D according to an incremental method consistent with the Jauman derivative yields a close approximation of the logarithmic strain, which is commonly used for simulations in the context of finite deformations as in the case of Abaqus. Using this Jaumann derivative in the incremental constitutive equations of s in the viscous and hysteresis contributions (11) and (13), leads to this natural decoupling of the associated (logarithmic) deformation measurement between volume and shape changes. As shown in Eq. (4), the spherical potential is directly linked to V alone. In contrary, as shown by Eq. (7), the Hart-Smith potential is based on volumetric-insensitive invariants. In conclusion, the present HVH model establishes an accurate decoupling between the volumetric and distortional parts. In the subsequent identification, K will be used as an independent coefficient and we will see that the response of the model will not be affected by any of the other parameters in an oedometric test. In addition, the material model will be naturally quasi-incompressible if the internal mechanical power required to change the volume is much larger than that required to change the isochoric distortion. In the HVH model, a sufficiently high value of K will ensure quasi-incompressibility. This can be easily illustrated in the case of a linear elastic model where K  G. In an

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

5

reversible nature of the whole model is therefore given by a combination of the hyperelastic and elasto-hysteresis contributions. It is not difficult, however, to numerically evaluate a proportion of hyperelasticity corresponding to the elasto-hysteresis as a function of the temperature so as to refine the springback of the whole model. This property will be used to account for the decrease in the contact force and its loss in the O-ring test. 2.3. Temperature dependence The temperature dependence first affects the values of the material parameters (i.e. all the material parameters will be temperature dependent), and secondly, the thermal expansion process. The thermal retraction/expansion processes will be taken to be governed by a classical linear behavior expressed by eth = aT • DT • I

(16)

where a T is the coefficient of thermal expansion and DT is the difference in temperature with respect to a reference temperature. This relation gives eth , the thermal part of strain, which is subtracted from the total strain so as to compute the behavior of the material: s = F (e − eth )

(17)

Fig. 1. Schematic representation of the volumetric and isochore parts of a strain path.

uniaxial compression test, for example, the ratio between the strains in the transverse and loading directions will be approximately 0.5, which is the classical theoretical incompressibility value. In the phenomenological approach, the main partitioning usually performed is that between the non-viscous (se + sh ) and viscous (sv ) parts. It is assumed that the viscous part vanishes entirely after a sufficiently long period of time under a constant strain, as previously described by Lion [40]. Fig. 2 gives a typical illustration of this assumption, and shows the method used here to separate the viscous and non-viscous coefficients in an uniaxial relaxation test. ¯ e + sh can be said to be The non-viscous stress s∞ = Se + s the residual stress after a long-term period of relaxation [22]. During the relaxation period, all the viscous parameters can therefore be identified separately, whether the stress s ∞ is hyperelastic or elasto-hysteretic. The separation between reversibility and non-reversibility is another topic of interest, but this cannot be done completely satisfactorily in the case of the present HVH model because the elasto-hysteresis model includes an elastic part (parameter l in Eq. (13)), which can be reversible below a threshold stress value. The

Fig. 2. Diagram of the viscous and non-viscous parts of the stress in a relaxation uniaxial test.

where F denotes the HVH model. The various contributions to the model can have temperature dependent parameters, which can be evaluated by performing linear interpolation over a tabulated set of identified values. These values are obtained by applying the identification process presented below at several temperatures corresponding to the experimental conditions. This method is used directly with the Hart-Smith coefficients C1 , C2 and C3 and the Maxwell coefficients l i , Gi . However, a simple linear interpolation is not suitable for the hysteresis part because of the history dependency during a temperature variation. Special attention is therefore paid to determining whether or not the history of the material recorded at a given temperature may still apply at another very different temperature. In this study, the elastohysteresis model accounts for a very long term state of viscosity which is assumed to be non-reversible. This behavior is specific to a certain state of the elastomer with respect to the temperature. It would not be realistic to assume, for example, that a mechanical history at −30 ◦ C still applies when the temperature returns to 20 ◦ C. In addition, the thermal independence of the hysteresis history plays an important role from the numerical point of view. The parameter Q0 is a stress saturation that cannot be exceeded by the intensity of sh . The lower the temperature, the more the material is assumed to show elasto-plastic like behavior. The value of Q0 will therefore be higher at low temperatures and the hysteresis part sh will be able to reach a high levels. If the temperature subsequently increases again, Q0 will decrease and a huge value of sh might become incompatible with a low Q0 limit. In order to solve this mechanical and numerical problem, the elasto-hysteresis part is split into several smaller parts, each of which is dedicated to a specific temperature range. Interlinked rules on mixtures are used in the following formulas: ⎫ sh = a1 sh,T1 + (1 − a1 )sh,T2 ⎪ ⎪ ⎪ sh,T2 = a2 sh,T2 + (1 − a2 )sh,T3 ⎪ ⎬ sh,T3 = a3 sh,T3 + (1 − a3 )sh,T4 with T1 > T2 > · · · > TN ⎪ ⎪ ⎪ ... ⎪ ⎭ sh,TN = sh,TN

(18)

6

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

where each elasto-hysteresis component sh,Ti has its own parameters l i and Q0,i and is computed from Eq. (13). a i are mixed parameters varying linearly in the [0:1] range in line with the following rules: T < Ti+1 : ai = 0 T−T Ti+1 ≤ T ≤ Ti : ai = T −Ti+1 i

T > Ti : ai = 1

i+1

(19)

In addition, it is most important to note that no sh,Ti is computed if a i = 0. This means that the i-th elasto-hysteresis component does not inherit the previous mechanical history if the temperature is outside of the [Ti : Ti+1 ] range. 3. Experimental database 3.1. Description of the material The material under investigation here was a hydrogenated nitrile butadiene rubber (HNBR) corresponding to the formulation and mechanical database presented in [14]. The present study focused on the mechanical responses of the material, as briefly described in the following sections. However, the range of the glass transition temperature and that the coefficient of linear thermal expansion (CLTE) have to be specified here. On the basis of DSC data, the material has a Tg of about −28 ◦ C. DMA gave a frequency-dependent Tg of −18 ◦ C at 1 Hz, −13 ◦ C at 10 Hz and −7 ◦ C at 95 Hz. The thermal expansion was determined by performing dilatometry on a cubic sample (5 × 5 × 10 mm) in line with the ISO standard 11359-2. In this case, Tg was found to be equal to −25 ◦ C. The value of the expansion coefficient above the Tg was found to be around 175e−6 K−1 and below the Tg , it was found to be around 70e−6 K−1 . 3.2. Mechanical tests The data used here were composed of an oedometric test performed at 23 ◦ C, isothermal relaxation tests under compression performed at various temperatures in the [80:−34] ◦ C range and anisothermal relaxation tests at several squeeze in the [5:20] % range at temperature ranging between room temperature and −60 ◦ C [14]. The relaxation of the reaction force of a O-ring was also studied for validation purposes. This test was carried out at the Cetim for the purpose of this study. The experimental conditions and the set-up used are described in [21]. 3.2.1. Oedometric tests A cylindrical sample, 18 mm high and 25 mm in diameter, was loaded into a containment rig with a piston at a rate of 1 mm/min at

Fig. 3. Force-displacement response in an oedometric test at 23 ◦ C (taking the stiffness of the device into account).

a temperature of 23 ◦ C. To perform these tests, the faces of sample were lubricated with silicone oil to obtain a homogeneous strain. The rig stiffness was subtracted from the raw data. The response curve is shown in Fig. 3. The force evolved slowly at the beginning of the test while the sample was not yet entirely in contact with the rig. The actual oedometric test started when the response curve showed a significant inflection. 3.2.2. Isothermal stress relaxation tests Isothermal stress relaxation tests were performed under uniaxial compression on a cylindrical sample 30 mm high and 18 mm in diameter. A strain of 25% was applied in the direction of the height at a rate of 100%/min. The sample was then subjected to a stress relaxation at a constant strain for 3 h. The test was ended by unloading the force down to 0 N at the same rate. As shown in Fig. 4(a), the stress versus time curves had similar patterns in the range of temperatures studied [80:−17] ◦ C. Only the stress level at the end of the loading phase increased as the temperature decreased, whereas, as can be seen from Fig. 4(b), the response evolved significantly from −25 ◦ C, showing an increasing initial slope during the loading. This evolution was induced by the loss of freedom of the polymer chains when glass transition occurred. The behavior of the material showed a plateau during the loading phase at −34 ◦ C and the stress was ten times greater than that observed at −17 ◦ C. Based on these experimental data, a classical time/temperature equivalence was done with the following relations: sT∗

ref

∗

=

T Tref

exp

sT

(20)

log(t ) = log(t) − log(aT ) where Tref is a reference temperature, sT∗ exp

ref

is the extrapolated stress,

sT is the experimental stress at temperature T, t∗ is the extrapolated time and a T is the shift factor versus the experimental time. This method was used for eight days. The main result obtained with this procedure was the presumably stabilized stress after long time relaxation as shown in the diagram in Fig. 2. Fig. 5 shows the evolution of stress versus the temperature. These stabilized stresses were useful for identifying the Hart-Smith and elasto-hysteresis parts. They correspond to the stabilized, i.e. non-viscous part of the HVH model as a function of the temperature. It can be noted in Fig. 4 that the main part of the relaxation stress was due to viscous effects. Fig. 6 gives the viscous proportion of the total stress during the relaxation period obtained by performing several time/temperature extrapolations. DMA results are also given. In this figure, the Tg values obtained with several methods are recalled. The green curves corresponds to the experimental results i.e. without any extrapolation. The red curve corresponds to an extrapolation performed after an eight days and 6 h of relaxation, at which the nonviscous part is shown on Fig. 5. In this case, the viscous proportion of the total stress was approximately in the [50:90] % range. The other curves were drawn up with various extrapolation time. The maximum values tended to occur at −25 ◦ C, which is consistent with the glass temperature obtained in thermal expansion tests. 3.2.3. Anisothermal stress relaxation tests In these tests, a cylindrical sample 13 mm high with a diameter of 29 mm was compressed diametrically. To study the evolution of the mechanical response with the thermal gradient and the squeeze of an O-ring, we chosen to apply a compression on a “compression set” sample in diametrical direction. The aim of this specific test was to describe the mechanical behavior of a cross-section of O-ring and to measure a significant evolution changes in the reaction force. An initial squeeze was applied and kept constant during a temperature cycle. After a 12-hour stress relaxation period at 20 ◦ C, the temperature was decreased from 20 ◦ C to −60 ◦ C in 10 ◦ C steps with holding

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

7

Fig. 4. Stress versus time curves during 3 h of relaxation under uniaxial compression.

times of 4 h and finally reset at 20 ◦ C a single step. The test ended with a stress relaxation step at 20 ◦ C. This test was performed at several compression levels relative to the initial diameter: 5, 7.5, 10, 12.5, 15, 17.5 and 20%. The temperature was measured as near as possible to the sample. The results are given in Fig. 7 in the case of

Fig. 5. Stabilized stress with an eight days time-temperature equivalence. The continuous curve was obtained by performing a local quadratic interpolation.

Fig. 6. The viscous part of the total stress in uniaxial relaxation tests and the loss factor obtained by DMA. The curves recorded after 15 min and 3 h correspond to the original isothermal relaxation data and the other curves correspond to several data extrapolated on the basis of time/temperature equivalence (d = day, m = month).

three examples at relative compression levels of 10, 15 and 20%. In all these cases, when the temperature decreased, the force decreased, except when the temperature reached the Tg , at which the force remained null and a loss of contact force occurred. However, the reaction force was recovered when heating up to 20 ◦ C was applied. These overall experimental conditions are more representative of an O-ring test than a simple isothermal relaxation test. In this test, the strain field in the sample is heterogeneous. The thermal expansion has a considerable effect on the response of the material in terms of the reaction force and the loss of contact. In addition, this test has a useful practical application in terms of the identification of the HVH model. First, after a short loading phase, the whole test is a relaxation test in which the strain rate D is null during this step. From the modelling point of view, the elasto-hysteresis part is therefore not activated after the end of the loading since it is governed by an incremental material law depending on the total strain rate (see Eq. (13)). Secondly, the temperature change is slow during this test. And apart from the thermal effects such as the expansion, it can be assumed that there will be no additional viscous effects during the temperature cycle. The viscous part of the HVH model is therefore not activated. In conclusion, the evolution of the reaction force is due only to thermal retraction and the hyperelastic part. If the thermal coefficient is known a priori, this test provides a very good means of studying the thermal evolution of the hyperelastic part of the HVH model. In the subsequent identification of material parameters, the loss of force during the temperature steps down to −30 ◦ C will be used. Below −30 ◦ C, the data are not usable because of the total loss of contact force. For example, the relative change of force from 10 to 0 ◦ C

Fig. 7. Relaxation tests with temperature cycle at three compression rates (samples loaded diametrically).

8

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

4. Material identification This section gives an overview of the strategy used to identify the parameters of the material in the HVH model. The following three sections explains each step. An additional section compares the numerical results with those obtained in the various mechanical tests used for the identification procedure.

4.1. Identification procedure

Fig. 8. Relative loss of reaction force at several compression rates. Results obtained on cylindrical samples and in O-ring tests are shown on this graph.

is taken to be the difference between the forces at the times 16 and 20 h. Fig. 8 reports this calculation for the first five temperature steps. This graph gives the results obtained using the same procedure in the O-ring test presented in the following section. The total loss of contact force occurring between −20 and −30 ◦ C corresponds to the temperature step including the glass transition temperature, −25 ◦ C measured using thermal expansion methods. 3.2.4. O-ring relaxation tests under cyclic temperature conditions The validity of the identified model will be tested by simulated these O-ring tests. The seal studied here was an O-ring with an internal diameter of 50.17 mm and a cross-section diameter of 5.33 mm [21], in keeping the ISO 3601-329 (which are the same as the BS-329) standards for O-rings. It was subjected to stress relaxation at a constant squeeze during temperature cycle. A 12 -hour stress relaxation was first applied at 20 ◦ C. The device was then subjected to a cooling phase from 20 to −60◦ C followed by a heating phase from −60 to 20 ◦ C in 10 ◦ C steps during both phases and a holding time of 4 h at each step. The reaction force and the temperature of the vacuum chamber and test rig were monitored. The force and temperature profiles are given in Fig. 9 in the case of a 24% squeezed seal. A total loss of contact occurred between −20 and −30 ◦ C and the force was gradually recovered during the heating. As shown in Fig. 8, the loss of contact which occurred during the temperature steps was similar to the 10% squeezed cylindrical sample.

Fig. 9. Temperature-dependent evolution of the force during the relaxation O-ring test.

Basically, the HVH model involves 17 coefficients: K, C1 , C2 , C3 , l, Q0 , np and 5 × (ti , Gi ). The temperature-dependent evolution of these coefficients will be expressed as a linear interpolation of discrete values. Each of them therefore had to be identified at each of the temperatures of interest, which was quite a difficult task. In the next section, some simplifications will be made in order to reduce the number of independent coefficients to 12 at each temperature, by means of a constant K value and a coupling between the hyperelastic parameters on one side and the hysteretic coefficients on the other side. The minimum experimental temperature was −34 ◦ C in the isothermal tests. In the anisothermal tests, the maximum temperature was equal to 20 ◦ C. No identification could therefore be performed below −34 ◦ C or above 20 ◦ C. In the subsequent simulations, the behavior of the material will not be taken to evolve outside this temperature range of [20:−34] ◦ C. As explained before, the HVH model is therefore clearly separated into its volumetric/distortional parts and viscous/non-viscous parts giving a general guideline for identifying the material parameters. The procedure used for this purpose involved the following steps: 1. Identification of the volumetric part Se (spherical hyperelasticity) From the constitutive equations presented in Section 2.1.2, the bulk modulus K was identified using the oedometric data presented in Section 3.2.1. This evaluation was based on trialand-error tests and therefore does not require a detailed description. 2. Identification of the viscous part sv (the generalized Maxwell model) In this case, using the equations presented in Section 2.1.4, the parameters 5 × Gi and 5 × ti were identified based on the data extrapolated from the isothermal relaxation tests presented in Section 3.2.2. As described in [22], a semi-analytical method followed by an inverse method were used to obtain the generalized Maxwell coefficients. In this step, only the relaxation phase was simulated. The bulk modulus K was known and fictitious Hart-Smith parameters were adopted to obtain the stress value at the end of the relaxation phase (s ∞ in Figs. 2 and 5). ¯ e + sh (Hart-Smith 3. Identification of the non-viscous part s potential and elasto-hysteresis model) During this step, from the equations presented in Sections 2.1.3 and 2.1.5, the parameters C1 , C2 , C3 and l, Q0 , np were identified. The experimental data used for this purpose were based on the anisothermal stress relaxations tests presented in Section 3.2.3 in conjunction with the isothermal tests described in Section 3.2.2. The proportion of the elasto¯e hysteresis part sh relative to deviatoric hyperelastic part s was evaluated in a simplified simulation of the relaxation test under a temperature cycle. The idea here was to perform a simulation involving all the other aspects of the model, such as the spherical part, the viscosity and the thermal expansion. Some coupling between parameters was performed in order to reduce the number of free coefficients.

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

9

Fig. 10 summarizes all these three identification steps. Each step is then described in detail. 4.2. Step 1: the spherical part The bulk modulus K is set at 2700 MPa This value is in a good agreement with the results of the oedometric test, as shown in Fig. 11. Oedometric tests results were available only at a temperature of 23 ◦ C. A sufficiently high bulk modulus ensures quasiincompressibility of the material. Its value can be easily set because the spherical hyperelastic part Se is independent of the other contributions involved. A constant value of K is assumed to have no temperature dependence although it is possible that the compressibility of the present HNBR may vary with the temperature. This assumption seems to be reasonable as long as the material is not confined, which was the case in all the tests of the experimental database, including the O-ring tests. 4.3. Step 2: the viscous part As explained above, five Maxwell branches were used to simulate the experimental stress relaxation responses and the corresponding time/temperature extrapolations. The viscous part sv was separated from the non-viscous part se + sh based on the uniaxial relaxation tests performed at each temperature. In Fig. 2, the maximum stress obtained at the end of the loading was equal to the sum s e + s v + s h whatever the nature of the non-viscous stress s e + s h . By giving s h a null value, fictitious Hart-Smith coefficients can be adopted for viscous identification purposes. Fictitious C1 , C2 and C3 coefficients were obtained using a Newton method to simulate the stabilized stress

Fig. 11. Spherical part of the model vs oedometric tests at 23 ◦ C.

after a long period of relaxation, based on the principle given by the equations Se (t∗ ) + s¯ e (t∗ ) = s∞ , sh = 0

(21)

where s ∞ is the stabilized non-viscous stress obtained from a time/temperature equivalence, as shown in Fig. 5, and t∗ is the maximum extrapolated time. The choice of t∗ should be as long as possible, but a period of eight days was chosen here. In fact, this depends on the duration available in the experimental database and on the possibility of reasonable extrapolating of these data. As mentioned above, only the deviatoric part of the generalized Maxwell model was used. The Poisson ratio mi has no influence and

Fig. 10. Summary of the step by step identification procedure.

10

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

Fig. 12. Time intervals on the five Maxwell branches at any temperature.

therefore does not need to be identified. From Eq. (12), only the shear moduli Gi and the characteristic times ti on each branch i have to be found. As described in [22], the method used here started by using a semi-analytical procedure to determine the characteristic time ti of the Maxwell cells. The shear moduli Gi were then determined using an inverse method involving the simulation of isothermal stress relaxation process with fictitious Hart-Smith parameters. This procedure was repeated at each temperature available based on the experimental isothermal tests at 20, 10, 0.5, −9, −13, −17, −25, −29 and −34 ◦ C. Each characteristic time ti was always computed at the same time intervals ti with t∗ = i ti . Fig. 12 gives an overview of this time distribution. This arbitrary choice was made in order to study the thermal evolution of each branch in a few fixed time domains. The approximate times ti were t1  100 s, t2  1/2 h, t3  2.5 h, t4  3.5 days and t5  4.5 days. The first three branches corresponded to the experimental relaxation time range (3 h). The fourth and fifth branches were dedicated to the extrapolated behavior occurring up to eight days and 6 h. Each characteristic time had to obey the relation 4ti < ti+1 to make each Maxwell branch independent, as assumed in [22]. On the one hand, the validity of the semi-analytical method depends on this inequality. On the other hand, the inverse method of Gi identification is better conditioned if each branch has no major effects on any of the others. The process of Gi identification focused on the maximum stress applied at the end of the loading and the long-term relaxation. The short-term relaxation is less useful to capture because the final purpose of this study was to model the relaxation of an O-ring subjected to constant long-term loading conditions. Figs. 13 and 14 give the thermal evolution of the parameters ti and Gi obtained. The characteristic times were found to be practically

Fig. 13. Thermal evolution of identified Gi parameters.

Fig. 14. Thermal evolution of identified ti parameters.

constant with the temperature. Based on Eq. (12) and the constant characteristic times ti , the evolution of gi was found to be similar to that of Gi and is therefore not reported here. Only the first branch had a characteristic time which decreased slightly with the temperature. However, the shear modulus evolved strongly at low temperatures. In particular, the first and second branches, which were associated with the lowest characteristic times, had the highest Gi values. The kinetic of the relaxation was greatly modified near the glass transition temperature. This fact can be observed again in Fig. 15 by comparing the normalized stress during the relaxation phase of the experimental tests. This trend reflects the occurrence of a proportionally faster decrease in the stress at low temperatures. The choice of fixed times ti in the modelling process gives a phenomenological understanding of the thermal evolution of the experimental relaxation kinetics. 4.4. Step 3: the non-viscous part The non-viscous part was studied using both Hart-Smith potential and the elasto-hysteresis model. The strategy usually adopted to identify this part consists in performing uniaxial test interrupted by relaxation steps [22]. But the experimental database used here did not include any tests of this kind. The results of anisothermal stress relaxation tests were therefore used to complete the identification. As mentioned in Section 3.2.3, anisothermal tests provide a means of studying hyperelasticity independently. More specifically, the main factor responsible for force variations is the temperaturedependent change in hyperelasticity. Nevertheless, the rest of the model is also involved and should be taken into account, especially the thermal expansion. Fortunately, the CLTE, the spherical part and the viscous part have already been determined at this stage in the identification procedure. The final step consists in finding the

Fig. 15. Normalized stress vs time during the relaxation part of several isothermal tests.

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

Fig. 16. Examples of the typical evolution of the Hart-Smith behavior in tensile tests (basic coefficients taken from Eq. (22) and some examples of material parameters from Table 1).

most suitable proportion of hyperelasticity relative to the elastohysteresis so as to able to simulate both the decrease in the force shown in Fig. 8 and the isothermal uniaxial long term stress shown in Fig. 5. The Hart-Smith potential and the elasto-hysteresis model have several parameters which can be used to handle the pattern of the stress-strain behavior. Without performing uniaxial tests interrupted by relaxation steps, it is difficult to exactly find the most suitable setup among the many possibilities. The present study focused on the amount of each stress contribution under a 25% strain in the uniaxial tests and a 10% compression during the anisothermal stress relaxation phase. An arbitrary choice was made by linking up some of the parameters in order to reduce the number of free coefficients. The following rules were adopted for linking the Hart-Smith material parameters together and the elasto-hysteresis material parameters together:

C1 Hart-Smith: C2 = 10 and C3 = Elasto-hysteresis: l = 4 • Q0 ; np = 2

C1 100

(22)

These relations give the uniaxial stress-strain patterns shown in Figs. 16 and 17, in the case of the Hart-Smith and elasto-hysteresis models, respectively, during loading and unloading of 25% of strain. In order to provide some other examples, these figures also show the behavior in the case of some subsequent obtained material parameters. With the conditions proposed by equations in Eq. (22), the proportional stress can be managed by only two parameters, C1 and Q0 . Their values are set using an inverse method. A Newton method

Fig. 17. Examples of the typical evolution of the elasto-hysteresis behavior in tensile tests (basic coefficients taken from Eq. (22) and some examples of material parameters from Table 1).

11

gives values of C1 which match the relative loss of force observed in anisothermal tests. At each change in the value of C1 , a secondary Newton method gives the corresponding Q0 consistent with long term isothermal stress. Interpolating into s ∞ in Fig. 5 simplifies the procedure because there is no exact correspondence between isothermal and anisothermal test temperatures. Lastly, this procedure gives the material parameters C1 and Q0 at the temperatures used in anisothermal tests: 10, 0, −10, −20 and −30 ◦ C. The non-viscous behavior at 20 ◦ C only requires to satisfy the uniaxial s ∞ since the contribution of elasto-hysteresis is assumed to be null (Q0 = 1e−6 MPa). For a similar reason, the behavior of the material at −34 ◦ C is obtained by matching the corresponding uniaxial s ∞ and assuming for the Hart-Smith contribution to be null (C1 = 1e−6 MPa). In order to complete the identification procedure, a quadratic interpolation is performed over these values in order to define additional C1 values at the experimental temperatures used: −29, −25, −17, −13 and −9 ◦ C. This approach gives a slightly better discretization of C1 . From the practical point of view, the FE simulation of the anisothermal test requires three-dimensional modelling, a nonlinear contact approach and time-consuming material integration. In addition, due to the duration of the experimental tests, the computation has to be performed for more than 50 h of physical time. This simulation is not suitable for using an inverse method in a reasonable CPU time. In the identification procedure based on Newton methods, a simplified simulation was therefore performed with a unit hexahedral element. This can be regarded as a Representative Volume Element (RVE) approach to modelling the actual heterogeneous test on a cylindrical sample. However, all the thermal, temporal and mechanical aspects are involved here. The unit element is loaded at a 10% compression rate by contact, the thermal expansion process is activated and the actual physical time is used in the whole HVH model given by Eq. (2). The two parameters obtained in this way are shown in Fig. 18 versus the temperature. The value of C1 decreases strongly when approaching the glass transition temperature, and reaches a null value at low temperatures. The parameter Q0 evolves in the opposite way. The elasto-hysteresis model gives the entire non-viscous behavior observed at temperatures below Tg . 4.5. Final set of material parameters At each temperature, only 8 parameters, K, C1 , Q0 and 5 × Gi were identified using the inverse method. The parameters ti were obtained using a semi-analytical procedure. The other parameters can be computed from Eq. (22) in the case of C2 , C3 , l and from Eq. (12) in that of gi . Table 1 shows all the parameters which were identified. The thermal evolution of the Hart-Smith and Maxwell parameters was

Fig. 18. Thermal evolution of parameters C1 and Q0 versus temperature.

12

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

Table 1 Material parameters in the HVH model in function of temperature. Hyperelasticity T (◦ C) 20 C1 (MPa) 4.98e−1 K (MPa) 2700 C2 = C1 /10, C3 = C1 /100

10 4.57e−1

0 3.89e−1

−10 2.08e−1

−20 2.43e−2

−30 1e−6

−34 1e−6

10 6.79e−2

0 1.47e−1

−10 4.28e−1

−20 7.88e−1

−30 2.14

−34 3.51

Elasto-hysteresis T (◦ C) Q0 (MPa) l = 4Q0 , np = 2

20 1e−06

Generalized Maxwell (5 branches) T (◦ C) G1 (MPa) t1 (s) G2 (MPa) t2 (s) G3 (MPa) t3 (s) G4 (MPa) t4 (s) G5 (MPa) t5 (s)

20 0.86 8.44 0.16 189.01 0.11 2.60e+3 0.08 2.12e+4 0.07 2.26e+5

10 1.01 8.67 0.20 210.54 0.13 2.21e+3 0.10 2.99e+4 0.08 2.84e+5

0.5 1.21 8.30 0.29 232.63 0.16 2.57e+3 0.10 3.13e+4 0.08 2.65e+5

−9 1.95 7.40 0.66 230.37 0.24 2.65e+3 0.22 2.91e+4 0.09 2.78e+5

−13 3.83 7.02 1.09 231.83 0.32 2.40e+3 0.23 2.84e+4 0.14 2.84e+5

−17 7.83 6.18 1.77 206.51 0.40 2.31e+3 0.27 3.08e+4 0.17 3.03e+5

−25 45.38 5.42 3.34 184.75 0.98 2.36e+3 1.05 2.97e+4 0.31 2.19e+5

−29 72.12 4.82 4.01 186.24 1.83 2.76e+3 1.36 2.74e+4 1.06 2.60e+5

−34 83.51 4.73 5.13 164.30 1.59 2.37e+3 1.13 2.82e+4 2.76 2.56e+5

For each branch i: gi = 2Gi ti .

obtained by performing linear interpolation over these values. The bulk modulus K is a constant. The elasto-hysteresis has a special thermal evolution given by Eqs. (18) and (19). The corresponding temperatures Ti are those given in Table 1. The thermal expansion coefficient, obtained from [14], is assumed to be constant and equal to 175e−6 K−1 . The material parameters do not evolve at temperature above 20 ◦ C or below −34 ◦ C in any of the components of the model because of the lack of experimental data available for performing parametric identification. 4.6. Simulations of experimental tests In the following sections, the results of the numerical simulations of the isothermal stress relaxation test and anisothermal stress relaxation test are presented. The material parameters used for this purpose are given in Table 1. 4.6.1. Isothermal stress relaxation The isothermal stress relaxation test was simulated on a single linear axisymmetric element. Comparisons between the experimental and numerical results can be made by looking at stress/time and

stress/strain curves in Figs. 19 and 20, respectively. At all temperatures studied, the HVH model was in good agreement with the relaxation phase and the maximum stress recorded at the end of the loading phase. In addition, the behavior of the material at 20 ◦ C was very accurately identified in all the phases (the loading, relaxation and the unloading phases). This is a necessary prerequisite, since the O-rings are loaded and undergo a long relaxation phase at this temperature. At temperatures under Tg during the loading phase, a discrepancy arose between the data obtained with the HVH model and the experimental results. This difference was greater at the lowest temperatures. But the long term residual stress was modelled accurately. This is an important advantage from the point of view of future O-ring test simulations. 4.6.2. Anisothermal stress relaxation As the result of the inverse process, the simplified simulation of the anisothermal stress relaxation process gave the relative loss of force at 10% compression rate based on the data in Fig. 8. It was worth checking the results in the case of the non-idealized 3-D simulation. The numerical model used here is shown in Fig. 21. Only 1/8 of the sample was modelled using the three symmetries. The loading was applied by friction-less contact with a rigid plate. The results are presented in Fig. 22 in terms of the contact force. The predictions of the

Fig. 19. Numerical results obtained in simulated isothermal stress relaxation tests: stress/time curves (NB: experimental data extrapolated after 3 h).

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

13

Fig. 20. Numerical results obtained in simulated isothermal stress relaxation tests: stress/strain curves.

HVH model were in good agreement with the experimental data. The model data matched the loss of contact occurring between −20 and −30 ◦ C and the recovery of the contact force when the temperature returned to 20 ◦ C.

5. Validation of the model: O-ring test The O-ring test was simulated using the data shown in Section 3.2.4. This test confirmed the validity of the HVH model on a near-industrial case. The numerical model is shown in Fig. 23 with meshes composed of axisymmetric elements. A 24% squeeze was applied to the O-ring in the vertical y-direction by friction-less contact with a coarse mesh consisting of undeformable elements. The X displacements of the extreme top and bottom nodes of the O-ring mesh were fixed. The results obtained with the HVH model during the cooling phase were similar to the experimental data, as can be seen from Fig. 25. The total loss of reaction force was found to occur between

Fig. 21. Numerical model for anisothermal stress relaxation (elastomer: 450 nodes and 76 full quadratic interpolation prismatic elements; plate: 1 rigid hexahedral element).

Fig. 22. Simulation of the anisothermal stress relaxation test in the case of a 10% compression rate.

−20 and −30 ◦ C, as observed experimentally. The stress in the loading direction was negligible at −30 ◦ C as shown in Fig. 24. A significant difference was observed during the heating phase from −60 to 10 ◦ C. However, the total recovery at 20 ◦ C was in good agreement with the experimental test results. 6. Discussion The identification procedure used here on the HVH model involves previously described methods [22] to which several adaptations were made in order to use the present HNBR experimental database. Whatever the experimental results used, the only condition was to be able to extract from them some information about each contribution to the HVH model. Modelling spherical behavior is usually quite straightforward if a purely volumetric test is available. However, the parameter K can vary with the temperature. This would be of some importance if the seal was totally constrained, but this occurs only when the seal undergoes a high in-service pressure, which was not the case in the present study. The viscous Maxwell model and the corresponding identification strategy are highly suitable when the relaxation behavior of material plays an important role in the industrial application under investigation. The uniaxial stress relaxation can always be interpreted in terms of a viscous part and a non-viscous part. Time/temperature extrapolation improves this approach by giving a better definition of the long term stabilized stress. In the present study, the generalized Maxwell model was used because we focused on the stress relaxation. The material evolved considerably during the loading and unloading phases at low temperatures. The great stiffness observed at the beginning of the loading and unloading phases at low temperature was mainly due to the viscosity because the stress was mainly viscous, as shown in Fig. 6. In the present study, the non-viscous part of the material behavior was modelled using a new approach in comparison with previous

Fig. 23. Numerical model of the O-ring test (O-ring: 281 nodes and 128 quadratic axisymmetric triangular elements); test rig: a coarse grid consisting of undeformable axisymmetric elements).

14

J. Troufflard et al. / Materials and Design 156 (2018) 1–15

Fig. 24. Numerical results of the O-ring test: stress in the loading direction. a) at room temperature and b) below glass transition temperature.

studies on the HVH model (see for example [22]). Anisothermal stress relaxation tests have been found to be a useful means of studying hyperelasticity separately. As shown in Fig. 18, the hyperelasticity disappeared completely when the temperature was lower than −25 ◦ C. In a simulated anisothermal test, the predicted loss of contact occurred when the hyperelasticity no longer generated any stress. It can be seen from Fig. 7 that all the tests resulted a loss of contact before −30 ◦ C, when the glass transition occurred. Thus, the evolution of material parameters therefore ensures that the total loss of contact will be accurately described by the model at temperatures ranging between −20 and −30 ◦ C. However, the stress-strain behavior of the Hart-Smith part is quite linear. The identified model will therefore give rise some problems when it comes to matching the sensitivity to the non-linear compression rate shown in Fig. 7. On the other hand, several simplifications have been made by coupling the Hart-Smith coefficients and the elasto-hysteresis coefficients, using the relations in Eq. (22). This procedure constrains the behavior to adopt one particular shape, as shown in Figs. 16 and 17. Accounting for the behavior using C1 and Q0 is the first step towards using the existing experimental database and to obtain overall proportions. The method used to link the Hart-Smith coefficients could be improved by expressing the relative loss of force at all the compression rates used in the anisothermal tests. One of the most interesting results obtained here was the ability of the model to describe the recovery of the contact force when the O-ring was heated. From the numerical point of view, the method based on a combination of rules presented in Section 2.3 helps to overcome the problem due to a risk of exceeding the Q0 threshold in the elasto-hysteresis model. For example, this threshold amount to practically 0 at a temperature of 20 ◦ C. This means that even a numerical artefact occurring during a temperature cycle without any additional loading will be liable to give rise to this problem. From the modelling point of view, the memorization of the thermomechanical history of the material is the main question. The behavior

of the present HNBR elastomer changes considerably at low temperatures. However, its original properties can be restored by heating the material. The elasto-hysteresis model is subdivided into the following independent temperature intervals: [20:10], [10:0], [0:−10], [−10:−20] and [−30:−34] ◦ C. The thermo-mechanical history accumulated at one temperature interval will not have to affect the behavior of the material at another sufficiently different temperature. The use of mixed rules is the first line of investigation which comes to mind, but this would have to be integrated directly into the elasto-hysteresis model. The memorization of the loading history during temperature variations will also have to be first investigated experimentally in order to define the most suitable means of including this process in the model. 7. Conclusions In this study, the hyperelasto-visco-hysteresis model was identified depending on the temperature in the case of a HNBR elastomer used for sealing applications. The aim was to use an existing experimental database and find a set of parameters which can be used to simulate an O-ring test under a temperature cycle including the glass transition temperature. An original strategy was introduced to model the loss of contact force of the seal occurring at low temperatures. Taking the volumetric/distortional and viscous/non-viscous parts separately is an integral part of the modelling framework which served as the basis of the identification procedure and the experimental data required. This study has shown that anisothermal stress relaxation tests give a viscous-free evolution of the contact force. In conjunction with uniaxial tests, this approach can be used to study the thermal evolution of the non-viscous mechanical behavior of the material and that of the model. Lastly, the step by step procedure developed here yielded a set of parameters which are suitable for predicting the contact force in a O-ring test above and below the glass transition temperature. Acknowledgments The authors thank the CITEPH program (http://www.citeph.fr) and Saipem, Schlumberger and TOTAL, France for their sponsorship and their support in this research in the framework of the Citeph-26-2014 project.

References

Fig. 25. Comparison between the numerical and experimental results of the O-ring tests.

[1] A. Douglas, O. Devlen, M. Mitchell, D. Edwin-Scott, M. Neal, Development of a procedure to accurately measure the low-temperature operating limits of elastomeric seals, Seal. Technol. 2016 (9) (2016) 7–12. https://doi.org/10.1016/ S1350-4789(16)30296-3. [2] P. Warren, Low temperature sealing capability of elastomer O-rings, Seal. Technol. 2008 (9) (2008) 7–10. https://doi.org/10.1016/S1350-4789(08)70478-1. [3] DIN ISO 815-2 rubber, vulcanized or thermoplastic - determination of compression set - part 2: at low temperatures 2016.

J. Troufflard et al. / Materials and Design 156 (2018) 1–15 [4] G. Spetz, Review of test methods for determination of low-temperature properties of elastomers, Polym. Test. 9 (1) (1990) 27–37. https://doi.org/10.1016/ 0142-9418(90)90046-G. [5] A. Mostafa, A. Abouel-Kasem, M. Bayoumi, M. El-Sebaie, On the influence of CB loading on the creep and relaxation behavior of SBR and NBR rubber vulcanizates, Mater. Des. 30 (7) (2009) 2721–2725. https://doi.org/10.1016/j.matdes. 2008.09.045. [6] M. Jaunich, A. Kömmling, D. Wolff, Investigations of elastomeric-seals-lowtemperature performance and ageing behaviour, Deformation and Fracture Behaviour of Polymer Materials, Springer Series in Materials Science Springer, Cham, 2017, pp. 431–443. https://doi.org/10.1007/978-3-319-41879-7_30. [7] H.-P. Weise, H. Kowalewsky, R. Wenz, Behaviour of elastomeric seals at low temperature, Vacuum 43 (5) (1992) 555–557. https://doi.org/10.1016/0042207X(92)90076-9. [8] M. Jaunich, W. Stark, D. Wolff, A new method to evaluate the low temperature function of rubber sealing materials, Polym. Test. 29 (7) (2010) 815–823. https://doi.org/10.1016/j.polymertesting.2010.07.006. [9] T. Grelle, D. Wolff, M. Jaunich, Temperature-dependent leak tightness of elastomer seals after partial and rapid release of compression, Polym. Test. 48 (2015) 44–49. https://doi.org/10.1016/j.polymertesting.2015.09.009. [10] S. Ronan, T. Alshuth, S. Jerrams, N. Murphy, Long-term stress relaxation prediction for elastomers using the time-temperature superposition method, Mater. Des. 28 (5) (2007) 1513–1523. https://doi.org/10.1016/j.matdes.2006.03.009. [11] J. Martinez, A. Boukamel, S. Méo, S. Lejeunes, Statistical approach for a hyper-visco-plastic model for filled rubber: experimental characterization and numerical modeling, Eur. J. Mech. A. Solids 30 (6) (2011) 1028–1039. https:// doi.org/10.1016/j.euromechsol.2011.06.013. [12] T. Cui, Y.J. Chao, J.W. Van Zee, Sealing force prediction of elastomeric seal material for PEM fuel cell under temperature cycling, Int. J. Hydrog. Energy 39 (3) (2014) 1430–1438. https://doi.org/10.1016/j.ijhydene.2013.10.086. [13] T. Grelle, D. Wolff, M. Jaunich, Leakage behaviour of elastomer seals under dynamic unloading conditions at low temperatures, Polym. Test. 58 (2017) 219–226. https://doi.org/10.1016/j.polymertesting.2016.12.018. [14] F. Rouillard, P. Heuillet, B. Omnès, Viscoelastic characterization at low temperature on an HNBR compound for sealing applications, Constitutive Models for Rubber VIII, CRC Press. 2013, pp. 591–594. https://doi.org/10.1201/b14964107. [15] A.G. Akulichev, B. Alcock, A.T. Echtermeyer, Compression stress relaxation in carbon black reinforced HNBR at low temperatures, Polym. Test. 63 (Supplement C) (2017) 226–235. https://doi.org/10.1016/j.polymertesting.2017.08. 023. [16] A. Ilseng, B.H. Skallerud, A.H. Clausen, Tension behaviour of HNBR and FKM elastomers for a wide range of temperatures, Polym. Test. 49 (2016) 128–136. https://doi.org/10.1016/j.polymertesting.2015.11.017. [17] Z. Chen, T. Liu, J. Li, The effect of the O-ring on the end face deformation of mechanical seals based on numerical simulation, Tribol. Int. 97 (2016) 278– 287. https://doi.org/10.1016/j.triboint.2016.01.038. [18] X. Li, T. Bai, Z. Li, L. Liu, Influence of the temperature on the hyper-elastic mechanical behavior of carbon black filled natural rubbers, Mech. Mater. 95 (2016) 136–145. https://doi.org/10.1016/j.mechmat.2016.01.010. [19] Y. Li, Y. He, Z. Liu, A viscoelastic constitutive model for shape memory polymers based on multiplicative decompositions of the deformation gradient, Int. J. Plast. 91 (2017) 300–317. https://doi.org/10.1016/j.ijplas.2017.04.004. [20] A.G. Akulichev, A.T. Echtermeyer, B.N.J. Persson, Interfacial leakage of elastomer seals at low temperatures, Int. J. Press. Vessel. Pip. 160 (2018) 14–23. https:// doi.org/10.1016/j.ijpvp.2017.11.014. [21] B. Omnès, P. Heuillet, Leak tightness of elastomeric seal at low temperature: experimental and FEM-simulation, in: I.P. Bohdana Marvalova (Ed.), Constitutive Models for Rubber IX, Prague, Czech Republic, 2015, pp. 609–614. https:// doi.org/10.1201/b18701-106.

15

[22] H. Laurent, G. Rio, A. Vandenbroucke, N. Aït Hocine, Experimental and numerical study on the temperature-dependent behavior of a fluoro-elastomer, Mech. Time-Depend. Mater. 18 (4) (2014) 721–742. https://doi.org/10.1007/s11043014-9247-3. [23] B. Wack, J.M. Terriez, P. Guelin, A hereditary type, discrete memory, constitutiveequation with applications to some geometries, Acta Mech. 30 (1983) 9–37. https://doi.org/10.1007/BF01170438. [24] D. Favier, P. Guélin, A discrete memory constitutive scheme for mild steel type material theory and experiment, Arch. Mech. 37 (3) (1985) 201–219. [25] G. Blès, Bases thermomécaniques de la modélisation du comportement des matériaux tissés et des polymères solides, Université Joseph Fourier - Grenoble I, 2002, http://www.theses.fr/2002GRE10058. (Ph.D. thesis) [26] Herezh++: FEM Software for Large Transformations in Solids, in: G. Rio (Ed.), Certification IDDN-FR-010-0106078-000-R-P-2006-035-20600, Université de Bretagne Sud, dépôt APP (Agence pour la Protection des Programmes). 2006, [27] G. Rio, H. Laurent, G. Blès, Asynchronous interface between a finite element commercial software ABAQUS and an academic research code HEREZH++, Adv. Eng. Softw. 39 (12) (2008) 1010–1022. https://doi.org/10. 1016/j.advengsoft.2008.01.004. [28] P.J. Flory, Thermodynamic relations for high elastic materials, Trans. Faraday Soc. 57 (1961) 829–838. https://doi.org/10.1039/TF9615700829. [29] G. Holzapfel, Nonlinear Solid Mechanics, 2nd ed. ed., John Wiley and Sons, Ltd., Chichester, 2001. [30] J.C. Simo, R.L. Taylor, Quasi-incompressible finite elasticity in principal stretches. continuum basis and numerical algorithms, Comput. Meth. Appl. Mech. Eng. 85 (3) (1991) 273–310. https://doi.org/10.1016/00457825(91)90100-K. [31] L.J. Hart-Smith, Elasticity parameters for finite deformations of rubber-like materials, Z. Angew. Math. Phys. 17 (5) (1966) 608–626. https://doi.org/10. 1007/BF01597242. [32] G. Chagnon, E. Verron, G. Marckmann, L. Gornet, Development of new constitutive equations for the Mullins effect in rubber using the network alteration theory, Int. J. Solids Struct. 43 (22) (2006) 6817–6831. https://doi.org/10.1016/ j.ijsolstr.2006.02.011. [33] P. Guélin, Remarques sur l’hystérésis mécanique, J. Méc. Théor. Appl. 19 (2) (1980) 217–247. [34] G. Blès, W.K. Nowacki, A. Tourabi, Experimental study of the cyclic visco-elastoplastic behaviour of a polyamide fibre strap, Int. J. Solids Struct. 46 (13) (2009) 2693–2705. https://doi.org/10.1016/j.ijsolstr.2009.02.015. [35] B.S. Shariat, Y. Liu, G. Rio, Hystoelastic deformation behaviour of geometrically graded NiTi shape memory alloys, Mater. Des. 50 (2013) 879–885. https://doi. org/10.1016/j.matdes.2013.03.091. [36] H. Laurent, A. Vandenbroucke, G. Rio, N.A. Hocine, A simplified methodology to identify material parameters of a hyperelasto-visco-hysteresis model: application to a fluoro-elastomer, Model. Simul. Mater. Sci. Eng. 19 (8). (2011) 085004. https://doi.org/10.1088/0965-0393/19/8/085004. [37] Y. Jiao, J. Fish, Is an additive decomposition of a rate of deformation and objective stress rates passé? Comput. Methods Appl. Mech. Eng. 327 (2017) 196–225. https://doi.org/10.1016/j.cma.2017.07.021. [38] V. Mora, H. Laurent, G. Rio, Cumulated tensorial strain measure in the logarithmic rotating frame, C.R. Mec. 332 (11) (2004) 921–926. https://doi.org/10. 1016/j.crme.2004.07.005. [39] O.T. Bruhns, H. Xiao, A. Meyers, Self-consistent Eulerian rate type elastoplasticity models based upon the logarithmic stress rate, Int. J. Plast. 15 (1999) 479–520. https://doi.org/10.1016/S0749-6419(99)00003-0. [40] A. Lion, A constitutive model for carbon black filled rubber: experimental investigations and mathematical representation, Contin. Mech. Thermodyn. 8 (3) (1996) 153–169. https://doi.org/10.1007/BF01181853.