An Integrated Risk Assessment Methodology of In-Service Hydrogen Storage Tanks Based on Connection Coefficient Algorithms and Quintuple Subtraction Set Pair Potential

Abstract

At present, there have been a number of hydrogen storage tank explosions in hydrogen filling stations, causing casualties and property losses, and having a bad social impact. This has made people realize that the risk assessment and preventive maintenance of hydrogen storage tanks are crucial. Therefore, this paper innovatively proposes a comprehensive risk assessment model based on connection coefficient algorithms and quintuple subtractive set pair potential. First of all, the constructed index system contains five aspects of corrosion factors, material factors, environmental factors, institutional factors and human factors. Secondly, a combined weighting analysis method based on FAHP and CRITIC is proposed to determine the weight of each indicator. The basic indicators influencing hydrogen storage tanks are analyzed via the quintuple subtraction set pair potential and full partial connection coefficient. Finally, the risk level and development trend of hydrogen storage tanks in hydrogen filling stations are determined by a combination of the three-category connection coefficient algorithms and the risk level eigenvalue method. The results of our case analysis show that the proposed risk assessment model can identify the main weak indicators affecting the safety of hydrogen storage tanks, including installation quality, misoperation and material quality. At the same time, it is found that the risk of high-pressure hydrogen storage tanks is at the basic safety level, and the development trend of safety conditions holds a critical value. The evaluation results can help establish targeted countermeasures for the prevention and maintenance of hydrogen storage tanks.

1. Introduction

According to statistical data of H₂ stations, the number and scale of global hydrogen filling stations will continue to grow, and it is expected that the number of global hydrogen filling stations will exceed 2000 by 2026, as shown in Figure 1 [1]. At the same time, China has also released a series of hydrogen filling station policies in recent years [2]. In particular, in March 2022, the Development and Reform Commission had clearly proposed that a large number of hydrogen filling stations will be deployed before 2025 [3].

The rapid construction and operation of hydrogen filling stations will inevitably produce a series of problems, the most important of which is safety risk. Gaseous hydrogen filling stations are built in large quantities in China. They usually consist of seven systems, among which high-pressure key facilities such as the pressurization system, hydrogen storage system and hydrogen filling system have complicated structures [4,5].

According to the statistics of 48 accident cases in hydrogen filling stations and joint construction stations, the accident causes include operation errors or improper maintenance, sealing problems, design problems, etc., as shown in Figure 2 [6].

The safety of high-pressure hydrogen storage tanks cannot be ignored. If leakage occurs, the amount of leakage will be large, which will lead to more serious consequences than with other equipment [7,8]. In the last five years, the U.S., Norway and South Korea have successively experienced a number of hydrogen explosion accidents, most of which were caused by the failure of high-pressure hydrogen storage devices; the details of the accidents are shown in

Table 1. Typical hydrogen explosion accidents.

Date	Accidents	Cause of Accidents
10 June 2019	Norwegian joint venture hydrogen filling station explosion.	A mistake in the installation of a special structure in the end part of the high-pressure hydrogen storage tank led to a leak of hydrogen gas. It went from a minor leak to a concentrated rapid leak, which in turn led to a fire and explosion in a relatively small space.
1 June 2019	Hydrogen tanker fire in Silicon Valley, USA	As a result of the fire, there was a gas leak from about 5 to 10 hydrogen tanks, which in turn led to an explosion.
23 May 2019	Hydrogen fuel storage tank explosion in South Korea	In the course of performing a hydrogen test by water electrolysis, an operational error caused a massive explosion in a hydrogen tank with a capacity of 400 L.

[9,10,11]. It can be seen that high-pressure hydrogen storage device accidents occur frequently, so it is urgent to conduct hazard identification and risk assessment for high-pressure hydrogen storage devices so as to provide safe technical support for the hydrogen storage device.

At present, more and more scholars are focusing on the safety of hydrogen filling stations. Li et al. studied the hazard impact distance of major accidents in hydrogen filling stations [12]. Various types of accident consequences were assumed, including flash fire, jet fire, and cloud explosions in confined spaces. Zhao Wenqing conducted a comparative study of hydrogen filling stations in Asia, Europe, North America and other regions, and developed a risk-based algorithm to determine safe spacing, providing a basis for the scientific and reasonable formulation and revision of the safety spacing [13]. Kikukawa et al. conducted a risk assessment of a 70 MPa hydrogen filling station in Japan, identifying 721 accident scenarios and calculating the consequences such as explosion pressure and jet fire, and the results show that a safety spacing of 6 m is sufficient [14]. In addition, the international standard for risk management in hydrogen storage facilities is ISO 19880 “Risk Assessment and Management of Hydrogen Energy Infrastructure” [15]. ISO 19880 specifies safety recommendations for the risk management of hydrogen refueling stations, mitigation measures to enhance system safety, safe distances for hydrogen exhaust ports and protection measures against non-hydrogen hazards. It defines minimum design, installation, commissioning, operation, inspection and maintenance requirements for the safety of hydrogen refueling stations. It can be seen that this standard is a relevant requirement from the design stage, so it is difficult to grasp the risk status of key facilities in the hydrogen filling station during operation, and it is difficult to put forward targeted risk management and control strategies.

Most of the above studies focus on safe design during the building period, but safety evaluation during the operation period is equally important [16,17]. The development of risk assessment technology for key facilities of hydrogen filling stations in China is still at a preliminary stage, and various risk analysis methods are not systematic enough. To this end, it is necessary to further develop a comprehensive risk assessment methodology and conduct detailed assessments of key facilities, so as to prevent the occurrence of uncontrolled risk factors.

In this study, by taking advantage of the connection coefficient algorithms and quintuple subtraction set pair potential, a comprehensive risk assessment methodology is proposed. Detailed procedures of the proposed methodology and related algorithms are shown in Section 2. In Section 3, a practical application of the proposed integrated risk evaluation methodology for hydrogen filling station is presented. Section 4 presents the conclusions and recommendations of the study.

2. Procedures

The structure of the proposed methodology is shown in Figure 3. The structure includes the construction of the failure risk indicator system, the determination of indicator weights based on the FAHP and CRITIC methods, the analysis of vulnerability and development trends of key risk indicators, and a final risk level eigenvalue judgement.

Among them, the failure indicator system is based on the comprehensive analysis of multi-source information factors, hazard identification and indicator screening. The FAHP–CRITIC combined weighting method can comprehensively consider the importance of expert experience and data factors. The quintuple subtraction set pair potential can diagnose the vulnerability indicators in basic events. The full partial connection coefficient can be used to identify trends in each indicator.

The overall risk status of high-pressure hydrogen storage tanks is determined by comprehensively analyzing the risk level eigenvalues and trends of the three types of connection coefficient algorithms for high-pressure hydrogen storage tanks.

2.1. Determination of Index Weight

2.1.1. FAHP

FAHP is a subjective weighting method based on the traditional AHP and fuzzy mathematics theory [18,19,20]. Different from the traditional AHP method, FAHP needs to build a fuzzy complementary judgment matrix, obtain a fuzzy consistent matrix through mathematical transformation, and obtain the indicator weight using the row sum normalization method. The n indicator factors $u_1$, $u_2$, $u_3$,...,$u_n$ with the same affiliation at the same level are compared, and the importance scale is expressed by 1~9. The judgment matrix established is:

$$X_{n\times n}={\left(x_{ij}\right)}_{n\times n}$$

(1)

where $x_{ij}$ is the AHP importance scale of $u_i$ and $u_j$ relative to the upper layer factor, and $x_{ii}=0.5$, $x_{ij}=1-x_{ji}$.

The elements in the fuzzy complementary judgment matrix are summed by rows, as follows:

$$y_i={\sum }_{k=1}^n{x}_{ik}$$

(2)

The fuzzy consistent matrix $Y={\left(y_{ij}\right)}_{n\times n}$ is obtained by mathematical transformation, and the calculation formula is as follows:

$$y_{ij}=\frac{y_i-y_j}{2\left(n-1\right)}+0.5$$

(3)

The elements in the fuzzy consistent matrix are summed by rows and normalized, and the index weight is as follows:

$$W_i=\frac{{\sum }_{j=1}^n{y}_{ij}}{{\sum }_{i=1}^n{\sum }_{j=1}^n{y}_{ij}}=\frac{{\sum }_{j=1}^n{x}_{ij}+\frac{n}{2}-1}{n\left(n-1\right)}$$

(4)

where $i,j=\mathrm{1,2},...,n$.

The n-order eigenvalue matrix $W^{\ast }$ is constructed using the weight vector, and the consistency indicator of the fuzzy judgment matrix and its eigenvalue matrix is used to test its consistency. The conditions to satisfy the consistency test are as follows:

$$W^{\ast }=\frac{W_i}{W_i+W_j}$$

(5)

$$I\left(A,W^{\ast }\right)=\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n\left|x_{ij}+W_{ij}^{\ast }-1\right|}{n^2}\le \alpha$$

(6)

where, when $I(A,W)\le 1$, the judgment matrix has consistency, or else it needs to be reconstructed.

2.1.2. FAHP

CRITIC is an objective weighting method that uses the objective properties of the data for scientific evaluation [21,22]. It is a better method than the entropy weight method [23,24]. At present, the CRITIC method is widely used in risk assessment, and the method has certain advantages over the objective weighting method. The specific calculation steps are as follows.

Step 1: Contrast strength.

The contrast strength is expressed as the difference between samples within an indicator. The greater the contrast strength, the greater the corresponding weights. The contrast strength formula of the indicator is:

$$S_j=\sqrt{\frac{{\sum }_{i=1}^p{\left(x_{ij}-\frac{1}{n}{\sum }_{i=1}^n{x}_{ij}\right)}^2}{n-1}}$$

(7)

where $x_{ij}$ is data after standardization, $S_j$ represents the contrast strength of the jth index, n is the number of samples, and p is total number of indicators.

Step 2: Relevance.

The stronger the relevance between an indicator and other indicators, the smaller the conflict and the more the information is repeated. The formula for the relevance coefficient is:

$$R_j={\sum }_{i=1}^p\left(1-r_{ij}\right)$$

(8)

where $R_j$ is the relevance coefficient between the jth indicator and other indicators, and $r_{ij}$ represents the relevance coefficient between the ith indicator and the jth indicator.

Step 3: Objective weight determination.

The objective weight formula is:

$$W_{2j}=\frac{S_j\times R_j}{{\sum }_{j=1}^p{S}_j\times R_j}$$

(9)

where $W_{2j}$ is the objective weight of the jth indicator.

2.1.3. Combination Weighting Method Based on FAHP and CRITIC

The AHP and CRITIC are both widely used weight analysis methods [25]. Among them, AHP relies too much on expert experience to score, and so it is more subjective. CRITIC is free from the constraint of subjectivity, but ignores the importance of indicator factors. In order to reasonably determine the indicator weights of hydrogen storage tank failure risk evaluation, this paper invokes combination assignment thinking to determine the index weights. It combines the advantages of a subjective weight method and objective weight method, taking into account both expert experience in indicator weights and actual objective data. Compared with traditional AHP weighting methods, it can better reflect the true weight values of various indicators.

To this end, it is necessary to determine the ratio of subjective weight and objective weight, so that the comprehensive weight can better reflect the importance difference between multiple evaluation indicators. $W_1$ is the weight obtained according to the FAHP, and $W_2$ is the weight obtained according to the CRITIC. The linear combination W of the two weighting methods is:

$$W=\alpha W_1+\beta W_2$$

(10)

In the above formula, $\alpha$ and $\beta$ represent the weight distribution coefficients. In order to find the best comprehensive weight vector, it is necessary to find the optimal weight distribution coefficient to minimize the standard deviation between W and the two groups of weight vectors, namely:

$$x_i={\alpha }_i{w}_{1i}+{\beta }_i{w}_{2i}$$

(11)

$$\delta =\sqrt{\frac{{\sum }_{i=1}^j\left(x_i-{\sum }_{i=1}^j\frac{x_i}{j}\right)}{j}}$$

(12)

where $j=\mathit{max}\left(i\right)$, $0\le \alpha \le 1$, $0\le \beta \le 1$, $\alpha +\beta =1$.

2.2. Connection Coefficient Algorithms Based on Set Pair Analysis Method

The set pair analysis method was first constructed by Zhao Keqin [26,27]. The general form of the evaluation sample connection coefficient, the evaluation indicator value connection coefficient and the average connection coefficient is:

$$u=a+bI+cJ$$

(13)

where a is the same parameter, b is the difference parameter, and c is the opposing parameter, a, b, c∈[0, 1], and a + b + c = 1. I is the difference parameter coefficient, the value range of which is [−1, 1]; sometimes it only serves as a difference marker. J is the opposite parameter coefficient, which is defined as −1 and sometimes only serves as a sign of opposites.

2.2.1. Evaluation Sample Connection Coefficient

The evaluation sample connection coefficient is obtained according to the index number of each index value for the evaluation sample in the k-level evaluation grade [28,29]. The calculation formula of the evaluation sample connection coefficient is as follows:

$$u_{1i}=v_{1i0}+v_{1i1}{I}_1+v_{1i2}{I}_2+v_{1i3}{I}_3+v_{1i4}J$$

(14)

$$v_{1i0}={\displaystyle \sum _{j=1}^{n_1}{w}_j}, v_{1i1}={\displaystyle \sum _{j=n_1+1}^{n_1+n_2}{w}_j}, v_{1i2}={\displaystyle \sum _{j=n_1+n_2+1}^{n_1+n_2+n_3}{w}_j}, v_{1i3}={\displaystyle \sum _{j=n_1+n_2+n_3+1}^{n_1+n_2+n_3+n_4}{w}_j}, v_{1i4}={\displaystyle \sum _{j=n_1+n_2+n_3+n_4+1}^{n_1+n_2+n_3+n_4+n_5}{w}_j}$$

(15)

where $u_{1i}$ is evaluation sample connection coefficient of sample i. I₁, I₂ and I₃ are difference parameter coefficients. J is the opposite parameter coefficient. $v_{1ik}$ is the connection coefficient component of $u_{1ik}$. $w_j$ is the weight of the jth indicator. In addition, n₁, n₂, n₃, n₄, and n₅ are the number of indexes in which sample i falls into each evaluation level, respectively.

2.2.2. Connection Coefficient of the Evaluation Indicator Value

Based on the proximity of the sample value $x_{ij}$ of indicator j in the evaluation sample i to the level of the evaluation criteria, the connectivity coefficient for the value of the evaluation indicator can be calculated [30].

When the positive index $s_{0j}<x_{ij}\le s_{1j}$ or the reverse index $s_{0j}>x_{ij}\ge s_{1j}$,

$$\left\{\begin{array}{c}{u}_{2ij0}=1\\ u_{2ij1}=1-2\frac{\left(s_{1j}-x_{ij}\right)}{\left(s_{1j}-s_{0j}\right)}\\ u_{2ij2}=-1\\ u_{2ij3}=-1\\ u_{2ij4}=-1\end{array}\right.$$

(16)

When the positive index $s_{1j}<x_{ij}\le s_{2j}$ or the reverse index $s_{1j}>x_{ij}\ge s_{2j}$,

$$\left\{\begin{array}{c}{u}_{2ij0}=1-2\frac{\left(x_{ij}-s_{1j}\right)}{\left(s_{2j}-s_{1j}\right)}\\ u_{2ij1}=1\\ u_{2ij2}=1-2\frac{\left(s_{2j}-x_{ij}\right)}{\left(s_{2j}-s_{1j}\right)}\\ u_{2ij3}=-1\\ u_{2ij4}=-1\end{array}\right.$$

(17)

When the positive index $s_{2j}<x_{ij}\le s_{3j}$ or the reverse index $s_{2j}>x_{ij}\ge s_{3j}$,

$$\left\{\begin{array}{c}{u}_{2ij0}=-1\\ u_{2ij1}=1-2\frac{\left(x_{ij}-s_{2j}\right)}{\left(s_{3j}-s_{2j}\right)}\\ u_{2ij2}=1\\ u_{2ij3}=1-2\frac{\left(s_{3j}-x_{ij}\right)}{\left(s_{3j}-s_{2j}\right)}\\ u_{2ij4}=-1\end{array}\right.$$

(18)

When the positive index $s_{3j}<x_{ij}\le s_{4j}$ or the reverse index $s_{3j}>x_{ij}\ge s_{4j}$,

$$\left\{\begin{array}{c}{u}_{2ij0}=-1\\ u_{2ij1}=-1\\ u_{2ij2}=1-2\frac{\left(x_{ij}-s_{3j}\right)}{\left(s_{4j}-s_{3j}\right)}\\ u_{2ij3}=1\\ u_{2ij4}=1-2\frac{\left(s_{4j}-x_{ij}\right)}{\left(s_{4j}-s_{3j}\right)}\end{array}\right.$$

(19)

When the positive index $s_{4j}<x_{ij}\le s_{5j}$ or the reverse index $s_{4j}>x_{ij}\ge s_{5j}$,

$$\left\{\begin{array}{c}{u}_{2ij0}=-1\\ u_{2ij1}=-1\\ u_{2ij2}=-1\\ u_{2ij3}=1-2\frac{\left(x_{ij}-s_{4j}\right)}{\left(s_{5j}-s_{4j}\right)}\\ u_{2ij4}=1\end{array}\right.$$

(20)

In the above formula, the positive index indicates that the value of the index increases with the increase in the evaluation criterion level k and the negative index indicates that the value of the index decreases with the increase in the evaluation criterion level k. For positive indicators, $(s_{0j},s_{1j}]$, $(s_{1j},s_{2j}]$,$(s_{2j},s_{3j}]$, $(s_{3j},s_{4j}]$, and $(s_{4j},s_{5j}]$ are thresholds of levels 0, 1, 2, 3, and 4, respectively.

The connection coefficient can be used as a means of evaluating the proximity between the sample value $x_{ij}$ and the rating of the evaluation criteria $s_{kj}$. The variable fuzzy relationship is a relative difference function and the corresponding relative membership $v_{2ijk}^{\ast }$ is:

$$v_{2ijk}^{\ast }=0.5+0.5u_{2ijk}$$

(21)

After normalization, we can derive the connection coefficient component $v_{2ijk}$ of the evaluation index value, as shown in the following formula:

$$v_{2ijk}=v_{2ijk}^{\ast }/{\sum }_{k=0}^4{v}_{2ijk}^{\ast }$$

(22)

Then, the connection coefficient of the evaluation index value is:

$$u_{2ij}=v_{2ij0}+v_{2ij1}{I}_1+v_{2ij2}{I}_2+v_{2ij3}{I}_3+v_{2ij4}J$$

(23)

Therefore, the connection coefficient of the evaluation index value for hydrogen storage tank assessment sample i can be obtained, as follows:

$$\begin{array}{l}{u}_{2i}=v_{2i0}+v_{2i1}{I}_1+v_{2i2}{I}_2+v_{2i3}{I}_3+v_{2i4}J=\\ {\displaystyle \sum _{j=1}^{n_j}{w}_j{v}_{2ij0}+}{\displaystyle \sum _{j=1}^{n_j}{w}_j{v}_{2ij1}{I}_1+}{\displaystyle \sum _{j=1}^{n_j}{w}_j{v}_{2ij2}{I}_2+{\displaystyle \sum _{j=1}^{n_j}{w}_j{v}_{2ij3}{I}_3+{\displaystyle \sum _{j=1}^{n_j}{w}_j{v}_{2ij4}J}}}\end{array}$$

(24)

2.2.3. Average Connection Coefficient of the Evaluation Sample

The average connection coefficient of the evaluation sample is:

$$v_{ik}={\left(v_{1ik}{v}_{2ik}\right)}^{0.5}/{\sum }_{k=0}^4{\left(v_{1ik}{v}_{2ik}\right)}^{0.5}$$

(25)

$$u_i=v_{i0}+v_{i1}{I}_1+v_{i2}{I}_2+v_{i3}{I}_3+v_{i4}J$$

(26)

The evaluation sample connection coefficient and the evaluation indicator value connection coefficient reflect the degree of correlation between the evaluation sample and the level of the evaluation criteria at the macro level and micro level, respectively, and the connection coefficients obtained by the two types are consistent in terms of the overall trend and distribution [31].

The indicator value connection coefficient is obtained by fuzzy matching using the continuous proximity information between the indicator value of the evaluation sample and the level of the evaluation standard, and its distribution is more uniform than the connection coefficient of the index number based on the discrete value. In order to further optimize the reasonable distribution of connection coefficients and deeply analyze the correlation degree between evaluation samples and evaluation standard grades, the principle of minimum relative entropy is applied to couple the two kinds of connection coefficients.

2.2.4. Full Partial Connection Coefficient

The concept of a partial connectivity coefficient is based on set pair analysis [32,33]. The full partial connection coefficient $p_f\left(u\right)$ essentially describes the relatively determined development trend represented by the connection coefficient at the micro level, so it can be defined as:

$$p_f\left(u\right)=\partial u^{+}-\partial u^{-}=\left[\frac{a}{a+b}+\frac{b}{b+c}{I}_1\right]-\left[\frac{b}{a+b}+\frac{c}{b+c}{I}_2\right]$$

(27)

where $I_1$ and $I_2$ are the difference coefficients of partially positive and partially negative connection coefficients, respectively. The value interval is [−1, 1], which can be generally determined by the ratio value method.

When the full partial connection coefficient $p_f\left(u\right)$ is greater than zero, equal to zero or less than zero, the trend of the micro-level development expressed in terms of the connection coefficient can be diagnosed as positive development trend, critical trend or negative development trend, respectively. The indicator of negative development trend is the main factor causing the weakness of hydrogen storage tank.

2.3. Risk Level Eigenvalue Method

The risk level eigenvalue method can be used to calculate the risk level of the evaluation sample’s connection coefficient, the evaluation indicator value’s connection coefficient and the average connection coefficient [34]. The calculation formula of risk level evaluation based on the level eigenvalue method is as follows:

$$h_1\left(i\right)={\sum }_{k=1}^5{v}_{1ik}k$$

(28)

$$h_2\left(i\right)={\sum }_{k=1}^5{v}_{2ik}k$$

(29)

$$h\left(i\right)={\sum }_{k=1}^5{v}_{ik}k$$

(30)

2.4. Set Pair Potential

2.4.1. Division Set Pair Potential

The division set pair potential is calculated as follows:

$$s_{f1}\left(u\right)=a/c$$

(31)

When c is small, the division set pair potentials tend to be unstable [35]. For example, the connection coefficients 6 + 0.06I + 0.2J and 6 + 0.06I + 0.02J have little difference, but their division set pair potential $s_{f_1}\left(u\right)$ values are 30 and 300, respectively, which contain very obvious differences [36]. Moreover, the division set pair potential undergoes a large range of changes, so it is not easy to classify. Therefore, it is necessary to carefully analyze and apply the division set pair potential.

2.4.2. Quintuple Subtraction Set Pair Potential

The calculation formula for quintuple subtraction set pair potential is as follows:

$$s_{f2}\left(u\right)=\left(a-c\right)\left(1+b_1+b_2+b_3\right)+0.5\left(b_1-b_3\right)\left(b_1+b_2+b_3\right)$$

(32)

The above formula is obtained from the set pair potential formula of ternary subtraction through empirical correction via practical applications [37,38]. Obviously, $s_{f2}\left(u\right)\in [-1,1]$. The numeric size between the subtraction set pair potentials can be compared directly based on the $s_{f2}\left(u\right)$, which greatly simplifies the tedious work of determining the order of the set pair potentials.

The subtraction set pair potential $s_f\left(u\right)$ is divided into five potential levels, among which the reverse potential $s_{f2}\left(u\right)\in [-1.0, -0.6]$, partial reverse potential $s_{f2}\left(u\right)\in [-0.6, -0.2]$, equipotential $s_{f2}\left(u\right)\in [-0.2, 0.2]$, partial same potential $s_{f2}\left(u\right)\in [0.2, 0.6]$, and same potential $s_{f2}\left(u\right)\in [0.6, 1]$. Similarly, when $s_{f2}\left(u\right)\in [-1.0, -0.2]$, the indicator is the main factor leading to the poor grade of the evaluation object, and so it can be identified as a vulnerability indicator. Therefore, the vulnerability indicator that can be identified for hydrogen storage tank failure should be improved.

3. Application of the Proposed Methodology in the Hydrogen Filling Station

In order to verify the practicability and accuracy of the proposed risk assessment methodology, this paper takes a hydrogen filling station in China as an example for in-depth analysis. The hydrogen filling station covers an area of 880 square meters, and the maximum storage capacity in the station is 800 kg. The distribution of the equipment in the hydrogen filling station is shown in Figure 4.

Hydrogen is transported into the station from the hydrogen source end by a hydrogen bundle truck, compressed to 45 MPa (high) pressure by a hydrogen compressor and stored in a 45 MPa hydrogen storage tank. A two-gun hydrogen filling machine at the filling island can fill hydrogen into cars and buses with a filling pressure of 35 MPa; the specific equipment details are shown in

Table 2. Specific equipment details.

Equipment Name	Quantity	Remarks
Hydrogen Tubular Vehicle	8	The volume is 2.3 m3, and the hydrogen pressure can not exceed 20 MPa.
45 MPa hydrogen storage tank	9	It consists of 9 interconnected cylindrical pressure vessels, each of which has a volume of 0.77 m3.
Compressor	2	There are two compressors in total, one of which is a backup.
35 MPa hydrogen double-shot hydrogenation machine	1	The maximum pressure to refill the car is 35 MPa.

. Since the hydrogen filling station has been in safe operation for 16 years, much of its equipment may be aging, so it is necessary to carry out risk assessment on some key equipment.

3.1. Construction of the Failure Risk Index System for the Hydrogen Storage Tank in a Hydrogen Filling Station

Due to the special geographical location, natural environment and special nature of the hydrogen filling station, the index system that affects the safe operation of the hydrogen storage tank should be established first. The indicators were screened according to the principles of specificity, representativeness and practicability of the objectives. Combined with relevant accident cases, experts experience and a large number of studies, the indicators were further screened. Based on multi-source data, a safe operation index system of a high-pressure hydrogen storage tank in the hydrogen filling station has been established, with five aspects: corrosion factor, material factor, environmental factor, system factor and human factor, as shown in Figure 5.

3.2. Calculation of Index Weights Based on FAHP and CRITIC Method

Step 1: Determination of subjective weight based on FAHP.

The secondary indicators fuzzy complementary judgment matrix is constructed as follows:

$$A=\left[\begin{array}{l}0.5 0.1 0.4 0.3 0.2\\ 0.9 0.5 0.7 0.7 0.6\\ 0.6 0.3 0.5 0.4 0.3\\ 0.7 0.3 0.6 0.5 0.4\\ 0.8 0.4 0.7 0.6 0.5\end{array}\right]$$

The weight of the fuzzy judgment matrix A is obtained through row summation, transformation and normalization. The weight of the fuzzy judgment matrix is as follows:

$$W_1={\left(0.15, 0.245, 0.18, 0.2, 0.225\right)}^T$$

In order to test whether the weight of the fuzzy judgment matrix is reliable, a consistency test of the fuzzy judgment matrix is required. The n-order eigenvalue matrix $W^{\ast }$ of the fuzzy judgment matrix is constructed using the weight vector, as follows:

$$W^{\ast }=\left[\begin{array}{l}0.5 0 0.38 0.45 0.43 0.40\\ 0.62 0.50 0.58 0.55 0.52\\ 0.55 0.42 0.50 0.47 0.44\\ 0.57 0.45 0.53 0.5 0 0.47\\ 0.60 0.48 0.56 0.53 0.50\end{array}\right]$$

After the consistency check, $I\left(A,W\right)=0.10\le 0.1$. This shows that the weight value meets the consistency requirement of the fuzzy judgment matrix.

Step 2: Determination of objective weight based on CRITIC

Five experts in hydrogen energy-related fields were invited to score the five secondary indicators M₁~M₅ respectively, and MATLAB R2018b software was used to calculate the standard deviation corresponding to the contrast strength and the correlation coefficient corresponding to the correlation, so as to obtain the final objective weight. The calculation results are shown in

Table 3. Scoring of secondary indicators by 5 experts.

Secondary Indicators	Expert A	Expert B	Expert C	Expert D	Expert E	Standard Deviation	Correlation Coefficient	Index Weight W2
M1	7.8	8.6	7.9	8.4	7.3	0.46	5.84	0.38
M2	9.2	9.5	9.4	9.8	9.6	0.20	3.65	0.10
M3	7.2	7.9	8.3	7.1	8.4	0.54	3.15	0.24
M4	8.5	8.4	8.9	8.1	9.2	0.39	3.5	0.20
M5	9.1	9.0	9.2	9.2	9.5	0.17	3.24	0.08

.

Step 3: Determination of combination weight.

The combination weighting method can better reflect the importance difference between multiple evaluation indicators. According to the calculation methods of Formula (10)–(12) and W₁ and W₂ obtained in the first two steps, the following is obtained:

$$\begin{array}{l}W=\alpha W_1+\beta W_2\\ =\left(-0.23\alpha +0.38 0.145\alpha +0.10 -0.06\alpha +0.24 0.2 0.145\alpha +0.08\right)\end{array}$$

In order to obtain the minimum standard deviation of W, this paper uses the enumeration method to solve it, and obtains the minimum standard deviation when the subjective weight proportion is 0.76, as shown in Figure 6.

Therefore, the comprehensive weight W is:

$$\begin{array}{l}W=0.76W_1+0.24W_2\\ =\left(0.21 0.21 0.19 0.20 0.19\right)\end{array}$$

(33)

The weight calculation results show that for the second-level indicators M₁~M₅, M₁ and M₂ have the largest weights, and the weights of other indicators in descending order are M₄, M₃ and M₅. Similarly, according to the above method of calculating combination weights, the calculation results of basic event weights are shown in

Table 4. Scoring of third indicators by 5 experts.

Basic Event	Subjective Weight W1	Objective Weight W2	Relative Comprehensive Weight W3	Comprehensive Weight W
X1	0.4	0.31	0.4	0.15
X2	0.6	0.69	0.6	0.23
X3	0.42	0.38	0.41	0.04
X4	0.25	0.55	0.31	0.03
X5	0.33	0.07	0.28	0.03
X6	0.43	0.29	0.36	0.09
X7	0.35	0.46	0.41	0.10
X8	0.22	0.24	0.23	0.05
X9	0.3	0.19	0.29	0.06
X10	0.24	0.17	0.23	0.05
X11	0.23	0.48	0.26	0.05
X12	0.23	0.16	0.22	0.04
X13	0.47	0.11	0.38	0.03
X14	0.28	0.61	0.36	0.03
X15	0.25	0.28	0.26	0.02

.

3.3. Determination of Vulnerability and Development Trend of Key Risk Indicators

According to the actual situation of the site and expert experience, the evaluation indicators can be divided into five levels, namely, safe, relatively safe, basically safe, relatively dangerous and dangerous. The specific division scope of each grade is shown in

Table 5. Risk level classification.

Risk Level	Safe	Relatively Safe	Basically Safe	Relatively Dangerous	Dangerous
Risk value	1	2	3	4	5

.

Five teams of professors in the fields of corrosion, materials, environment, business management and human–computer interaction were invited to score the relevant indicators. We substituted the expert evaluation data of the index into Formulae (16)–(24) to obtain the connection coefficient of the evaluation indicator value, then calculated the full partial connection coefficient from Formula (27), and finally calculated the subtraction set pair potential according to Formula (35). The specific calculation results are shown in

Table 6. Connection coefficient of each assessment sample for basic events of hydrogen storage tank failure risk.

Basic Events	Evaluation Indicator Value Connection Coefficient	Quintuple Subtraction Set Pair Potential	Full Partial Connection Coefficient
X1	0 + 0$I_1$ + 0$I_2$ + 0.43$I_3$ + 0.57J	−0.91	−0.93
X2	0 + 0.35$I_1$ + 0.5$I_2$ + 0.15$I_3$ + 0J	0.1	0
X3	0 + 0$I_1$ + 0$I_2$ + 0.39$I_3$ + 0.61J	−0.92	−0.99
X4	0 + 0.01$I_1$ + 0.5$I_2$ + 0.49$I_3$ + 0J	−0.24	0
X5	0 + 0$I_1$ + 0$I_2$ + 0.29$I_3$ + 0.71J	−0.96	−1.13
X6	0 + 0$I_1$ + 0$I_2$ + 0.45$I_3$ + 0.55J	−0.90	−0.98
X7	0 + 0.03$I_1$ + 0.5$I_2$ + 0.47$I_3$ + 0J	−0.22	0
X8	0 + 0$I_1$ + 0.19$I_2$ + 0.5$I_3$ + 0.31J	−0.70	−0.55
X9	0 + 0$I_1$ + 0$I_2$ + 0.4$I_3$ + 0.6J	−0.92	−0.98
X10	0 + 0$I_1$ + 0.31$I_2$ + 0.5$I_3$ + 0.19J	−0.55	−0.35
X11	0 + 0$I_1$ + 0.43$I_2$ + 0.5$I_3$ + 0.07J	−0.37	−0.14
X12	0 + 0$I_1$ + 0.02$I_2$ + 0.5$I_3$ + 0.48J	−0.86	−0.80
X13	0 + 0$I_1$ + 0$I_2$ + 0.37$I_3$ + 0.63J	−0.93	−1.02
X14	0 + 0.03$I_1$ + 0.5$I_2$ + 0.47$I_3$ + 0J	−0.22	0
X15	0 + 0.22$I_1$ + 0.5$I_2$ + 0.28$I_3$ + 0J	−0.03	0

. As shown in Table 6, the basic events X1, X3, X5, X6, X8, X9, X12 and X13 are in the opposite direction, which gives the vulnerability index of the high-pressure hydrogen storage tank in the hydrogen filling station. This is the main prevention and control object related to safety management to ensure that the risk level is controlled within an acceptable level.

At the same time, it can be found that the full partial connection coefficient and the subtraction set pair potential of the basic event have a high degree of similarity in Figure 7. The main reason is that both can describe the uncertain state and development trend of indicators to a certain extent. Among them, the subtraction set pair potential focuses on identifying the uncertain state of each basic event, thereby identifying the vulnerable indicators. The full partial connection coefficient focuses on describing the development trend of each basic event.

Therefore, based on the identification of X1, X3, X5, X6, X8, X9, X12, and X13 vulnerability indicators, it can be found that the full partial correlation coefficients of vulnerability indicators are all less than 0, indicating a negative development trend. The sequence of negative development trends can be ranked as X5 > X13 > X3 > X6 = X9 > X1 > X12 > X8. Priority control should be given to installation quality, misoperation, material quality, natural disasters, internal corrosion, emergency planning, and air humidity and temperature, in order of priority.

3.4. Determination of Risk Level and Development Trend of High-Pressure Hydrogen Storage Tank

Based on the data in

Table 7. Assessment results of failure risk for a high-pressure hydrogen storage tank.

Connection Coefficient	Connection Coefficient Value	Evaluation Grade	Subtractive Set Pair Potential	Full Partial Connection Coefficient
Evaluation sample connection coefficient	0.4 + 0.14I1 + 0.46I2 + 0I3 + 0J	2.06	0.68	0.69
Connection coefficient of the evaluation indicator value	0 + 0.09I1 + 0.25I2 + 0.38I3 + 0.28J	3.85	−0.59	−0.50
Average connection coefficient	0 + 0.25I1 + 0.75I2 + 0I3 + 0J	2.75	0.13	0

and the index weights in Table 4, the final connection coefficient of the evaluation index number is obtained using Formulae (14) and (15). Then, the final connection of the evaluation index value can be obtained from Formulae (16)–(24), and we can obtain the final average connection coefficient of the evaluation sample from Formulae (25) and (26). Finally, we can determine the final risk level and development trend of a high-pressure hydrogen storage tank.

According to the comprehensive analysis of the three connection coefficients, the risk level of the high-pressure hydrogen storage tank is close to Level 3, which means it is basically safe. The subtraction set pair potential approaches 0.1 and is in equilibrium, indicating that the high-pressure hydrogen storage tank has a certain safety bearing capacity. The full partial connection coefficient approaches 0, indicating that the safety situation of high-pressure hydrogen storage tanks is undergoing a critical development trend.

Based on the above analysis results, it can be seen that the high-pressure hydrogen storage tank of the hydrogen filling station is at a medium risk level, and the risk development trend is in a critical state, as shown in Table 7. It is necessary to strengthen control over weak indicators X1, X3, X5, X6, X8, X9, X12, and X13. After our investigation, we find the risk status of high-pressure hydrogen storage tanks is consistent with the actual situation on site, which confirms the effectiveness and practicality of the proposed method.

4. Conclusions

In order to establish a highly applicable risk assessment method for high-pressure hydrogen storage tanks in hydrogen filling stations, this paper proposes a comprehensive risk assessment methodology based on connection coefficient algorithms and quintuple subtraction set pair potential. The methodology mainly includes the FAHP–CRITIC combined weighting method, quintuple subtraction set pair potential, the full partial coefficient method and the risk level eigenvalue method.

The FAHP–CRITIC combined weighting method is a good method that can comprehensively consider the subjective experience and objective data, so it can effectively avoid the uncertainty caused by a single-type weighting method. In order to effectively deal with the uncertainty between evaluation indicators and evaluation criteria, a quintuple subtraction set pair potential algorithm is proposed, and it is used to identify the uncertainty state of each indicator of a high-pressure hydrogen storage tank. At the same time, the full partial connection coefficient can identify the development trend of each indicator.

The usability and validity of the methodology proposed in this study were verified through a case study on the risk encountered by high-pressure hydrogen storage tanks. The results show that the risk level of a high-pressure hydrogen storage tank in the hydrogen filling station is grade 3, the subtraction set pair potential approaches 0.1, and the full partial connection coefficient approaches 0. By comparing the result with the actual situation, we found that the result was consistent with the actual situation on site.

Compared with traditional risk assessment methods, the proposed risk assessment method not only identifies the risk level and weak indicators of the target, but also determines the development trend of each weak indicator. The risk assessment results are more accurate and specific, which can provide direction to guide site managers in risk control.