Finite Element Analysis on the Temperature- Dependent Burst Behavior of Domed 316L Austenitic Stainless Steel Rupture Disc

Abstract

As a safety device, a rupture disc instantly bursts as a nonreclosing pressure relief component to minimize the explosion risk once the internal pressure of vessels or pipes exceeds a critical level. In this study, the influence of temperature on the ultimate burst pressure of domed rupture discs made of 316L austenitic stainless steel was experimentally investigated and assessed with finite element analysis. Experimental results showed that the ultimate burst pressure gradually reduced from 6.88 MPa to 5.24 MPa with increasing temperature from 300 K to 573 K, which are consistent with the predicted instability pressures acquired by nonlinear buckling analysis using ABAQUS software. Additionally, it was found that a gradual transition from opening ductile mode to cleavage mode happened with increasing temperature due to more cross slips occurring under serious plastic deformation. The equivalent stress, equivalent strain and strain hardening rates acquired by static analysis were effective at rationalizing the temperature-dependent fracture behavior of the domed rupture discs.

1. Introduction

Rupture discs, also known as bursting discs, are widely used in the nuclear power, fire protection and petrochemical industries owing to the advantages provided by their simple structure, high sensitivity and strong venting ability [1,2,3,4,5]. Hydrogen has especially been regarded as a promising alternative energy source; however, safety issues associated with high-pressure hydrogen hinders its application [6]. Several studies related to the spontaneous ignition of high-pressure hydrogen during the utilization of rupture discs have been conducted [7,8,9]. Presently, a rupture disc can be roughly classified as a domed rupture disc, a reverse domed rupture disc and a flat rupture disc. The domed and flat rupture discs undergo tensile failure of materials, while the reverse domed rupture disc undergoes structure instability failure [10,11,12].

Schrank [12] investigated the instability behavior of reverse discs by putting them through vacuum cycles to simulate flushing the facility and pressure cycles as emergencies with a sudden rise in pressure, and they found that the tested discs showed no sign of degradation after vacuum and pressure cycling compared to the unstressed discs. As rupture discs are a pressure relief component, one of the critical issues is to accurately predict their ultimate burst pressure (P_b). Lake and Inglis [13] suggested to estimate the burst pressure using

$$P_b=2.6{\sigma }_s\times \frac{S}{a}$$

(1)

based on a constant volume assumption, where ${\sigma }_s$ is the ultimate tensile strength measured from a conventional uniaxial tensile test, and s and a are the thickness and interdiamater of the domed rupture disc, respectively, as shown in Figure 1. Kanazawa et al. [14] proposed another equation

$$P_b=4{\sigma }_s\times {(\frac{2}{3})}^n\times \frac{S}{a}$$

(2)

based on the tensile instability conditions [15], where n is the strain hardening coefficient of the material.

The effect of environmental temperature can, to some extent, be indirectly identified from the tensile strength. Actually, the ultimate burst pressure of a domed rupture disc is closely linked to its relief caliber, thickness and fillet radius with regard to a particular material. Much effort has been made using commercial simulation software [16,17,18,19,20], for example, Jeong et al. [20] carried out a structural analysis of domed stainless-steel rupture discs from the viewpoint of materials failure, and the relationship amongst burst pressure, thickness and superficial groove depth were established; however, no experimental evidence was provided.

In the present study, we used austenitic stainless steel 316L as a model material, to predict its ultimate burst pressure. The equivalent stress, equivalent strain and strain hardening rates were also acquired using the finite element method (ABAQUS/CAE 6.14-4, Dassault Systèmes Simulia Corp., Providence, RI, USA) to better understand the temperature-dependent burst behavior of domed rupture discs.

2. Experimental Procedure

In the present study, a 316L austenitic stainless sheet steel with a thickness of 0.064 mm was adopted to fabricate the domed rupture discs. The chemical composition of the 316L austenitic stainless steel is listed in

Table 1. The chemical composition of the studied steel (wt.%).

C	Si	Mn	P	S	Ni	Mo	Cr	Fe
≤0.03	≤1.00	≤2.00	≤0.035	≤0.03	10.00–14.00	2.00–3.00	16.00–18.00	bal.

. The rupture discs were fabricated by press forming within a die with the dimensions shown in Figure 1. All the rupture discs were annealed at 473 K for 2 h to relieve stress. The bursting tests were conducted in a self-designed rupture set-up, as shown in Figure 2, at 300 K, 373 K, 473 K and 573 K, respectively, where the temperature was controlled by a K-type thermocouple welded between holders with a resolution of ±2 K, and the pressure was applied using compressed air. All the tests for each temperature were performed on at least five samples. The experimental process was as follows: (1) the rupture disc was fixed between the top and bottom holders; (2) the holders were heated to the preset temperature and held for 30 min; and (3) the pressure was applied by filling with compressed air until the fracture happened. The fracture surfaces were taken from the discs and directly examined using a scanning electron microscope (SEM; JSM-6700F, JEOL, Akishima Tokyo, Japan). The finite element analysis procedures are provided in the next section.

3. Results and Discussion

The experimental results showed that all the rupture discs burst at the dome, as shown in Figure 3, and the average ultimate burst pressure was 6.88 MPa at 300 K, gradually decreasing to 6.20 MPa, 5.59 MPa and 5.24 MPa when the temperature rose to 373 K, 473 K and 573 K, respectively. The fracture morphologies are shown in Figure 4. In the case of 300 K, one could observe some parabolic (elongated) dimples and cleavage features besides equiaxed dimples near the inner surface, as marked by the yellow arrows in Figure 4a. Such features were different from the full equiaxed dimples obtained during the conventional uniaxial tensile tests [21,22] and the small punch tests (SPT) [23,24]. It could be inferred that the severe plastic deformation was first initiated by the rapid tearing stress from the outside surface, and then the final fracture changed to an opening mode caused by the direct stress. With rising temperatures, as shown in Figure 4b–d, the proportion of equiaxed/elongated dimples gradually decreased and more cleavage features appeared. The cleavage mode dominated the fracture at 573 K.

In the present study, finite element analyses, including nonlinear buckling analysis and static analysis during successive loading, were carried out using ABAQUS software. A nonlinear buckling analysis was performed, as the linear elastic constitutive relationship was not applicable for local buckling of the rupture disc. Static analysis was conducted to acquire the equivalent stress, equivalent strain and strain hardening rates to explain the fracture behavior. A Nlgeom option was applied in both the nonlinear buckling analysis and static analysis. The physical properties of the 316L austenitic stainless steel are listed in

Table 2. The physical properties of the 316L austenitic stainless steel [25,26].

Temperature	300 K	373 K	473 K	573 K
Ultimate tensile strength (MPa)	632	590	492	485
Break elongation (%)	47	36	27	26
Young’s modulus (MPa)	195,000	189,000	183,000	176,000
Poisson ratio	0.3	0.3	0.3	0.3

[25,26].

A nonlinear buckling analysis can be used to simulate structural responses after loss of balance, and generally, it is conducted as follows: (1) to analyze eigenvalue buckling to calculate the eigenvalues and obtain buckling modes; (2) to introduce an appropriate number of eigenmodes into the perfect geometry as an initial imperfection; and (3) to obtain a reaction force (i.e., the force acting on the normal direction of the dome surface) versus the arc length curves using Riks method which was used to calculate the post-buckling strength and post-buckling mode [27,28,29,30]. A rigid body hemisphere with a radius of 12 mm was applied to provide pressure, and the flat clamped by holders was deformable, with a 360° revolution and preset thickness of 0.064 mm, as shown in Figure 5. The displacement of the hemisphere was set to be 5.4 mm to get a dome of the flat. An encastre option (neither displacement nor rotation) was applied on the edge part, and a displacement without rotation option was applied on the dome part, as shown in Figure 6a. A tet option (element shape) was chosen for the irregular and complex areas. A quadratic was adopted in the geometric order, and other options were subjected to defaults. In addition, we set the maximum load proportionality factor (the ratio of calculated load to applied load) in the stopping criteria (instability criteria) option as one and the DOF value (degree of freedom) as three (the three-dimension mode) on the top of the dome (the node region option). Meanwhile, the Riks analysis would automatically stop once the vault region was fully plastic by presetting the “other” option in the ABAQUS software. C3D10 was a 10-node quadratic tetrahedron for noncontact cases, as shown in Figure 6b. There are 57,153 elements in total and two elements along the wall thickness.

A series of reaction force (F) distributions during successive loading at different temperatures could be obtained. Two representative cloud maps at 300 K were taken as examples, as shown in Figure 7, and the corresponding plastic deformation distributions are shown in Figure 8. One can notice that the maximum reaction force/plastic deformation appeared at the top of the dome, which agreed well with the experimental results, as shown in Figure 3. By extracting the maximum reaction force and the corresponding arc length [31,32], the relationship curves between the reaction force and arc length were plotted in Figure 9a, and then the relationship between the pressure and arc length were obtained, as shown in Figure 9b. According to Figure 9b, the instability pressures were determined to be 7.20 MPa, 6.50 MPa, 5.97 MPa and 5.60 MPa at 300 K, 373 K, 473 K and 573 K, respectively, which were close to the experimental results; here, the instability pressure was defined to be the maximum reaction force divided by the original dome area. Meanwhile, the ultimate burst pressures estimated from Equations (1) [13] and (2) [14] were integrated into Figure 10 for comparison.

According to Figure 10, all the predicted data from the Lake model (Equation (1)) seem to be less than the experimental results, while the data from the Kanazawa model (Equation (2)) are slightly higher. The deviation is explained as follows. Firstly, the Lake model [13] was derived from

$$\sigma =\frac{P\times a}{S}\frac{{\left(1+4{\phi }^2\right)}^2}{16\phi }$$

(3)

where σ, p, S, φ and a are equivalent stress, pressure, thickness, relative arch height (vertical displacement divided by interdiameter of the disc) and relief caliber, respectively. The equation was proposed with an assumption that the rupture disc deforms uniformly; however, it is inevitable that the vault of the dome thins the most [33]. Secondly, Kanazawa model was obtained by modifying the equations of

$$\frac{P\times a}{S\times A}=\frac{2^{2+n}{\left[\ln \left(1+{\phi }^2\right)\right]}^n\times \phi }{{\left(1+{\phi }^2\right)}^3}$$

(4)

$$n=(5-\frac{1}{{\phi }^2})\times \frac{\ln \left(1+{\phi }^2\right)}{2}$$

(5)

where A is the hardness coefficient (kg/mm²) and n is the strain hardening coefficient. Although nonuniform deformation is considered in this model by introducing an arc of a smaller radius, the φ value becomes higher and thus results in an over-predicted value.

Comparatively speaking, the predicted values of the present work using the finite element method are closer to the experimental results, indicating the reasonability of the simulation model described above. P_total = P₀ + λ(P_ref − P₀) was introduced to the ABAQUS software to evaluate the buckling load, where P₀ is a constant value (dead load), P_ref is the reference load, and λ is the load proportionality factor. When conducting an ABAQUS/Riks analysis, P_ref is altered automatically by the ABAQUS program to obtain the corresponding P_total during the analysis process until P_total equals P_ref, representing that the structure reached stable (λ = 1). λ changes to maintain P_total equals P_ref if P_ref continues to increase when the structure has lost stability, and the nonlinear bucking strength is expressed as the product of λ and P₀ [27,32]_.

The equivalent stress (σ_T) and equivalent strain (ε_T) can be further determined by static analysis with the same finite element model. The experimental data of the ultimate burst pressure were applied. The equivalent stress and equivalent strain distributions of the domed rupture disc at 300 K are shown in Figure 11 as an example, and, correspondingly, the relationship between the equivalent stress and equivalent strain can be obtained by extracting data, as shown in Figure 12a. As the temperature rises, the equivalent stress in the dome part gradually decreases. Figure 12b presents the relationship between the strain hardening rate and the equivalent strain by differentiating the curves of Figure 12a. At the low strain stage before yielding, the deformation mode of the 316L austenitic stainless steel is planar slip with a high strain hardening rate [34,35]. As a face-centered cubic structured material, the deformation mode of the 316L austenitic stainless steel gradually changes from planar slip to cross slip with increasing strain, causing saturated dislocation tangles, which may exhaust the strain hardening capability [36]. On the other hand, the thermally activated process could intensify at a higher temperature, as reflected by the enlargement of the cleavage area in Figure 4; that is, the strain hardening capability could be more easily consumed, resulting in a faster burst fracture [37].

4. Conclusions

The effect of temperature on the ultimate burst pressure and fracture behavior of the domed 316L austenitic stainless steel rupture discs was investigated. Finite element analysis based on nonlinear buckling analysis and static analysis during successive loading was proven to be an effective way to predict the burst pressure. The following conclusions were obtained:

(1) Experimental results showed that the ultimate burst pressure decreased from 6.88 MPa to 5.24 MPa with increasing temperature from 300 K to 573 K. Instability pressures acquired by nonlinear buckling analysis under different temperatures agreed well with the experimental results.

(2) The fracture morphologies showed a gradual transition from an opening ductile mode to a cleavage mode with increasing temperature due to more cross slips occurring under serious plastic deformation. The equivalent stress, equivalent strain and strain hardening rates acquired by static analysis can be used to effectively rationalize the temperature-dependent fracture behavior of domed rupture discs.