Under the service conditions of ultra-high-pressure hydrogen storage vessels, cyclic loading—particularly concentrated at threaded connections—is a common occurrence. The non-uniform load distribution along the threaded engagement, coupled with pronounced stress concentration at the thread root, renders this region highly susceptible to fatigue crack initiation. These cracks typically nucleate at the thread root and propagate along the radial-circumferential plane, as illustrated in Figure 1.

Although the internal high-pressure hydrogen does not directly act on the crack surface, the severe stress concentration at the thread root generates a steep longitudinal stress gradient, which significantly accelerates crack propagation. In this study, to analyze the crack growth behavior in high-stress regions, the following assumptions are adopted: the crack propagates within a plane perpendicular to the direction of the maximum tensile stress, and the initial crack configuration is idealized as a semi-elliptical surface crack with an aspect ratio (a/ℓ) of 1/3. Variations in this aspect ratio during crack growth are considered to reflect the evolving crack morphology.

The initiation and propagation of cracks are predominantly governed by the combined effects of uneven load distribution and stress concentration at the thread root. Despite the absence of direct exposure of the crack surface to the high-pressure hydrogen, the localized stress concentration is sufficient to generate a sharp longitudinal stress gradient, promoting crack growth along the radial-circumferential plane. Accurate prediction of crack propagation behavior therefore requires precise determination of the crack initiation location and the propagation path, which typically aligns with the plane normal to the maximum tensile stress. Moreover, the continuous evolution of the crack aspect ratio during growth must be incorporated into the analysis to capture the dynamic characteristics of crack propagation under cyclic loading conditions.

The analysis process begins by setting initial conditions and key parameters, with the critical stress intensity factor K_c and hydrogen-induced crack threshold representing the material’s fracture toughness and resistance to hydrogen embrittlement, respectively. The analysis starts by assuming an initial crack size and an aspect ratio of 1/3, based on the geometry of common cracks in engineering.

Loading conditions are set to simulate the stress state of the material in the actual usage environment. Then, using crack geometry and stress distribution analysis, the stress intensity factor is calculated as long as the maximum stress intensity factor K_max does not exceed the threshold values of K_c or K_H, and the crack’s aspect ratio remains below 0.8. Based on this condition, the stress intensity factor range Δ_K is calculated Δ_K is negative, the crack is assumed to be closed under compressive loading, and its growth rate da/d_N is assumed to be 0 mm/cycle. In this case, crack propagation is not expected.

For positive values of Δ_K, the crack propagation rate is calculated using an equation that incorporates material constant C, a function f(R_k) of the stress ratio R_k, and the crack growth exponent m. This formula considers material properties and loading conditions, and is key to predicting crack propagation. By calculating the crack propagation rate, the growth rate of the crack is determined, which provides the basis for subsequent cycle number calculations.

Once the crack growth rate is determined, the next step is to calculate the crack growth increment Δ_(a) and the corresponding cycle number increment Δ_N. By iterating calculations of crack size and cycle number, the crack growth and cycle numbers are accumulated, and the complete crack growth trajectory with time is obtained. When the crack reaches its critical size, it is considered that the material can no longer withstand the load. The total failure cycle number N_d is the minimum number of cycles required for the crack to grow from the initial crack size to the critical size. This result can also be defined by another standard, which is the number of cycles required for the crack to grow to the allowed final depth. The smaller of these two values is taken as the final fatigue life prediction result.

This crack growth-based fatigue life analysis method, by combining material fracture toughness, stress distribution, and crack geometry characteristics, can accurately predict fatigue life. Additionally, the introduction of the hydrogen-induced crack threshold K_IH is considered, especially under extreme conditions such as ultra-high pressure hydrogen storage, making the model more reliable and scientific for predicting fatigue life in special environments. Ultimately, this method not only predicts the material’s service life but also provides a theoretical basis and data support for design optimization and safety evaluation.

Based on technical data from an experimental sample of an ultra-high-pressure vessel as shown in Table 1, it is conservatively assumed that no cracks deeper than 0.8 mm exist due to manufacturing requirements and non-destructive testing capabilities. Therefore, the maximum stress location on the vessel’s thread is assumed to have an initial crack depth of 0.8 mm. The crack type is assumed to be radial and circumferential cracks along the inner thread root. This conservative assumption is made to ensure a safe margin in crack impact assessment.

The ultra-high-pressure hydrogen storage vessel is subjected to a pressure cycling profile over its 20-year designed service life. The maximum pressure reaches 134.6 MPa, which corresponds to the hydrostatic pressure test conducted biennially. It should be noted that this test pressure does not represent the vessel’s normal operating conditions but is part of a periodic structural integrity assessment, simulating extreme pressure scenarios the vessel may encounter. Figure 2 illustrates this pressure cycling profile.

During regular operation, the vessel undergoes approximately 27,740 pressure fluctuation cycles between 71 MPa and 89.6 MPa, reflecting typical variations during hydrogen charging and discharging. Additionally, about 40 full pressure release cycles are considered, where the pressure rises from 0 MPa to 89.6 MPa, simulating complete filling and emptying events, 10 cycles from 0 to 134.6 MPa, which are conducted as part of periodic hydrostatic requalification tests.

This pressure cycling profile is crucial for assessing the vessel’s fatigue behavior and long-term structural integrity. The detailed load history supports quantitative analysis of the vessel’s durability and aids in validating design assumptions, predicting fatigue life, and developing safety guidelines and inspection protocols for ultra-high-pressure hydrogen storage systems.

Figure 3 illustrates the end-cap sealing mechanism and threaded connection system of the ultra-high-pressure hydrogen storage vessel. This structure plays a vital role in connecting the end-cap to the vessel body, providing sealing performance, and ensuring mechanical locking under high-pressure conditions.

As shown in Figure 3a, the end-cap is connected to the vessel through a threaded system with precisely defined geometric dimensions. A step reduction in diameter is designed at the connection region to enhance the load-bearing capacity and improve structural safety under cyclic pressure loads.

Figure 3b presents the detailed thread specification, identified as 4-8UN-2A/2B, indicating unified standard threads with specific pitch and tolerance grades. The accurate design and machining of the thread ensure the reliability of the connection and the sealing performance required for ultra-high-pressure operation.

To investigate the axial stress distribution of the threaded connection under cyclic loading, a local axisymmetric finite element model of the end-cap, vessel neck, and thread was developed, as shown in Figure 4. The model adopts the center of the end-cap’s top surface as the origin, with the Y-axis pointing upward, the Z-axis extending outward from the paper, and the X-axis oriented horizontally to the upper left. The 8-node axisymmetric solid element (PLANE183) is used for the structural components, while contact interactions at the thread and sealing surfaces are simulated using CONTA172 and TARGE169 elements. (2) Finite Element Model.

Boundary conditions are applied as shown in Figure 5. Axial displacement constraints are imposed at the bottom of the internal thread, and symmetry conditions are applied along the axial symmetry plane of the thread plug. Internal pressures of 71 MPa, 89.6 MPa, and 134.6 MPa are applied to the lower surface of the thread plug and the inner surface of the vessel neck, corresponding to different service and test conditions. This modeling approach enables an accurate evaluation of the stress distribution and structural performance of the threaded connection under high-pressure cyclic loading.

A local axisymmetric finite element model of the end cap thread connection was developed based on Saint-Venant’s principle and the structure’s axisymmetric characteristics. The model uses PLANE183 elements and considers geometric nonlinearity to capture the stress response under cyclic loading accurately. Axial displacement constraints were applied at the bottom of the internal thread, and symmetry constraints were imposed along the axisymmetric plane. Internal pressures of 71 MPa, 89.6 MPa, and 134.6 MPa were applied to simulate different working conditions.

As shown in Figure 6, the load distribution and boundary conditions of the model reflect the actual stress state of the thread connection under ultra-high pressure. Figure 7 shows the stress and strain contours, where stress concentrations are mainly located at the thread crest and root potential initiation sites for fatigue cracks. The strain distribution also reveals localized plastic deformation, indicating areas susceptible to damage under high cyclic loads.

The longitudinal stress distribution of the threaded connection was analyzed under internal pressures of 71 MPa, 89.6 MPa, and 134.6 MPa. As shown in Figure 8, at 71 MPa, stress is mainly concentrated at the thread crest and root, which are potential fatigue crack initiation sites.

With increasing pressure to 89.6 MPa (Figure 9), the overall stress level rises, and stress concentration becomes more pronounced in critical areas. At the maximum pressure of 134.6 MPa (Figure 10), stress near the thread root reaches its peak, indicating a significantly elevated risk of fatigue failure.

The radial stress variation along path A-A, from the inner to outer diameter, is shown in Figure 11 for three different internal pressure conditions. The stress is most concentrated near the inner thread diameter, with a marked increase as pressure rises. A sharp decrease in stress occurs at the transition to the outer diameter, indicating a more uniform stress distribution outside the thread region.

Figure 12 displays the stress distribution contour in the thread region of the ultra-high-pressure hydrogen storage vessel. Stress levels are represented by a color gradient, from blue (low stress) to red (high stress). The red areas highlight regions of maximum stress concentration, primarily caused by the thread geometry. These areas are critical for the thread’s fatigue performance and integrity, being potential sites for fatigue crack initiation.

The contour of plastic strain distribution in the thread region is shown in Figure 13. Using color coding, the diagram represents the magnitude of plastic strain, with blue indicating the lowest strain values and transitioning to green, yellow, and red as the strain increases. Areas of concentrated strain suggest that the material may have reached or exceeded its yield limit, resulting in plastic deformation.

The relationship between the axial stress s and the distance x from the inner surface is expressed as a cubic polynomial for each pressure condition, as follows (see Equation 1):

Using the least squares method, the cubic polynomial coefficients for each region under different pressures are determined, as shown in Table 2.

Based on ASME VIII-3-2019 Appendix D, Formula D-403, the stress intensity factor K₁ corresponding to the crack depth a is calculated for 71 MPa, 89.6 MPa, and 134.6 MPa (see Equation 2):

Where ( 1 − ( F S y ) 2 ) 1 / 6 —represents the plastic zone correction factor.

When used for crack propagation, it is set to 1. The crack propagation formula is) (see Equation 3):

The formula is transformed into (see Equation 4):

F1 is expressed as (see Equation 5):

F2 is expressed as (see Equation 6):

F3 is expressed as (see Equation 7):

F4 is expressed as (see Equation 8):

When considering the multi-region stress distribution model, the stress intensity factor K₁ should theoretically remain continuous at the boundaries between regions. However, due to differences in the cubic polynomial coefficients A₀, A₁, A₂, and A₃ between regions, the calculated K₁ values To address this, the difference in 𝐾₁ values at the boundaries between Region 1 and Region 2 and between Region 2 and Region 3, denoted as ΔK_1,2 and ΔK_2,3, calculated as shown in Table 3. This correction ensures the continuity of K₁ across the entire analysis range, allowing for more accurate crack propagation predictions.

In the K₁ calculation values for Region 2 and Region 3, in order to ensure continuity across regions and accurately predict crack propagation, we apply corrections to these values. The specific correction process involves adding ΔK_1,2 to the ΔK_1,2 calculation value in Region 2, and adding to the ΔK_1,2 + ΔK_2,3 calculation value in Region 3. This adjustment ensures the continuity of K₁ at the region boundaries, as presented in Figure 14, and allows us to obtain a corrected K₁ value, providing a more accurate calculation foundation for crack propagation analysis.

Based on the pressure cycle characteristics shown in Figure 3, the stress intensity factor at the crack tip is quantitatively analyzed under different cyclic pressures. Specifically, for each pressure cycle, the stress intensity factors at the deepest crack point and the crack surface under minimum cycle pressure (P_min) and maximum cycle pressure (P_max) are calculated, and the stress intensity factors K_1min and K_1max are derived. Then, the crack growth rates are computed using the following formulas (see Equation 9):

The stress ratio R is defined as (see Equation 10):

The crack growth rate at the deepest crack point is given by (see Equation 11):

According to KD-430 regulations, the crack growth at the current crack depth is (see Equation 12):

Using the expression for f(R_k) as defined in Code Case 2938-1, we calculate the crack growth rate at the crack surface f(R_k) (see Equation 13):

f(Rk) is expressed as (see Equation 14):

The crack growth rate at the crack surface is (see Equation 15):

From the initial crack depth a₀, the crack growth is iteratively calculated based on the loading order and pressure amplitude shown in Figure 3 for each cyclic pressure. After each cycle, the crack depth is updated as a = a + Δa, and this iteration continues until the maximum stress intensity factor K_1max at the crack tip reaches the critical hydrogen-induced cracking stress intensity factor threshold K_IH. After 100,000 cycles, the crack depth increases from the initial 0.8 mm to 1.3028925 mm. Under the highest working pressure of 89.6 MPa, the stress intensity factor at the crack tip is 15.478 MPa · m at the experimental pressure of 134.6 MPa, it is 23.217 MPa · m , both significantly lower than the critical threshold of K_1H = 40 MPa · m , indicating that the crack growth rate and final size meet the design cycle life requirements.

The design cycle life is defined as the minimum of the following two conditions:

The crack propagation calculation results show that expanding the crack from the initial depth of 0.8 mm to the critical depth of 2.874 mm requires N_d = Min(Nc/2, Np) = Min(122,638/2, 64,361) = 613,19 cycles, while expanding it to the allowed final depth of 1.319 mm requires.