Early Warning of Internal Leakage in Heat Exchanger Network Based on Dynamic Mechanism Model and Long Short-Term Memory Method

Abstract

In the process of butadiene rubber production, internal leakage occurs in heat exchangers due to excessive pressure difference. It leads to the considerable flow of organic matters into the circulating water system. Since these organic matters are volatile and prone to explode in the cold water tower, internal leakage is potentially dangerous for the enterprise. To prevent this phenomenon, a novel intelligent early warning and risk assessment method (DYN-EW-QRA) is proposed in this paper by combining dynamic simulations (DYN), long short-term memory (LSTM), and quantitative risk assessment (QRA). First, an original internal leakage mechanism model of a heat exchanger network is designed and simulated by DYN to obtain datasets. Second, the potential relationships between variables that have a direct impact on the hazards of the accident are deeply learned by LSTM to predict the internal leakage trends. Finally, the QRA method is used to analyze the range and destructive power of potential hazards. The results show that DYN-EW-QRA method has excellent performance.

1. Introduction

The production of butadiene rubber (BR) is an orientated polymerization process with butadiene as the monomer [1]. BR ranks second in production among all kinds of rubber product in the world, only behind styrene butadiene rubber. It has the advantages of good elasticity, strong wear resistance, and low temperature resistance. Due to its excellent properties, BR has been extensively applied in the production of impact-resistant plastics, tires, tapes, hoses, rubber shoes, and other rubber products [2]. Therefore, it is of great significance to ensure the safe and stable production of BR. In practice, the internal leakage of the heat exchanger network has always existed in the BR industry due to the deterioration of equipment. It is estimated that 10% of all corrosion damage in industrial systems results in leakage. However, internal leakage is usually hard to detect, e.g., 6.5 tons of liquid chlorine leaked unknowingly from a heat exchanger at Honeywell’s refrigerant production facility, leading to an emergency evacuation of the entire facility [3]. Pressure difference is generally one of the main causes of internal leakage [4]. The ionic components in the leakage flow may cause reduction reactions between water and the cathodic metallization and produce hydrogen gas and hydroxide ions [5]. Consequently, internal leakage caused by pressure difference and corrosion can lead to adverse safety accidents and environmental impacts. Early warning and risk assessment are absolutely necessary for avoiding such leakage accidents.

In some emergency evacuation cases, the improper dissemination of evacuation warnings is the main cause of casualties [6]. The cause analysis of major accidents shows that some accidents could have been avoided if appropriate preventive measures were put in place [7]. Nevertheless, many early warning studies are difficult to provide a reliable description of emergency operations and further related preventive measures [8]. In addition, the results of risk analysis are not comprehensive, so the accuracy and pertinence of monitoring and early warning are difficult to guarantee. Therefore, early warning methods need to be improved. Deep learning is considered as one of the most promising artificial intelligence technologies. It has attracted more and more research attention in recent years as early warning and data prediction using deep learning have great implication in the safety analysis of chemical processes. Based on deep learning, a generative adversarial networks-spearman’s rank correlation coefficient-deep belief networks (GAN-SRCC-DBN) method was proposed for the identification of abnormal conditions in chemical processes [9]. However, most methods are limited to statistical methods and unable to capture the dynamic changes in the data deeply. Bayesian network has been widely used in the field of probabilistic safety analysis, but it ignores the time variable and is incapable of reflecting the evolution of the system variables over time [10]. Another drawback of traditional neural networks is their inability to interpret time series based on specific context. Early attempts to solve this problem were Recurrent Neural Networks (RNN) based on the feedback method [11]. However, it frequently leads to a gradient disappearance problem in certain applications. Long short-term memory (LSTM) is capable of alleviating these problems as a variant of RNN [12]. LSTM was first proposed by Sepp Hochreiter and Jurgen Schmidhuber [13]. It can maintain memory in their cellular state to figure out sequence and time problems without losing gradients that affect their performance. Its greatest strength consists in that it can learn the time variable introduced by dynamic simulation platforms to predict changes of risk variables over time. Hence, LSTM can be utilized to study the underlying mechanism in the model; then, it can be used to obtain quantitative results in disturbed conditions.

Dynamic simulation plays a significant role in process control strategy design, validation, and safety analysis. It introduces a time variable that reflects transient changes into the process when the system is disturbed. Therefore, it is widely used in chemical processes as the major tool for system safety analysis. For example, a novel DYN-SIL method combining dynamic simulation with the safety integrity level (SIL) was put forward to analyze the risk of the fluid catalytic cracking (FCC) fractionating system [14]. Risk analysis and risk matrices combining dynamic simulation were presented to determine consequences and quantify severity [15]. A semi-batch production of nitrile from fatty acid and ammonia was simulated using dynamic simulation to reduce batch cycle time [16]. Dynamic simulation was used to formulate a suitable model for evaluating the overall dynamic characteristics of complex thermal recovery steam cycles [17]. In general, dynamic simulation can be used to investigate the impact of variable fluctuations on safety response time when the variables reach the critical limit. So, the internal leakage within the heat exchanger network is studied in this paper utilizing dynamic simulation.

In recent years, a large number of safety accidents occurred in the chemical industry, resulting in a large number of deaths. Risk assessment is an essential requirement to ensure the safe design and operation of processes [18]. Along with the increasing public appeal and strict government management, substantial sophisticated safety assessment methods and their improved methods have been widely employed in the past decade [19]. Especially, SIL and hazard and operability analysis (HAZOP) are two of the most commonly used methods. However, those traditional methods cannot measure risk with high credibility [20]. In order to address this problem, quantitative risk assessment (QRA) has attracted extensive attention from experts and scholars. QRA has been widely applied in various industries, including the nuclear industry, oil and gas industry, and steel industry [21]. For example, a new QRA method integrated with dynamic simulation and accident simulation was proposed to discover inherent risks that are undetectable by conventional risk analysis methods based on steady-state conditions [22,23]. A risk-based accident model was used to conduct QRA and find the potential explosion risk for leakage failure of a submarine pipeline [24,25]. Dynamic variables were added to QRA to improve the operational performance analysis of chemical processes [26]. Therefore, in this paper, QRA was employed to conduct the safety assessment of the heat exchanger network and obtain the quantitative ranges and destructive power of potential hazards caused by the internal leakage.

This paper realizes the early warning and risk assessment by combining DYN, LSTM, and QRA with regard to internal leakage of the heat exchanger network in the butadiene rubber polymerization unit. First, an internal leakage mechanism model of heat exchanger is built utilizing dynamic simulation. Based on this model, LSTM is trained to predict risk variables for constructing an early warning model. Finally, the ranges and destructive power of potential hazards are obtained using the QRA approach followed by the quantitative safety measures recommended. The remainder of this paper is arranged as follows. Section 2 describes the BR process and dynamic simulations. Section 3, Section 4 and Section 5 illustrates the heat exchanger network simulation, the early warning modeling, and risk assessment processes, respectively. Some important results are given in the final conclusion section.

2. Description of BR Process and DYN-EW-QRA Method

The BR production process consists of nine units: metering, polymerization, rubber tank, coagulation, rubber washing and drying, briquetting, packaging, solvent recovery, and waste gas treatment. Among them, the polymerization unit is the major research object. The glue generated by the polymerization reaction leaves from the top of the first reactor; then, it enters the bottom of the second reactor and other reactors in series. The jacket of each polymerization reactor is filled with frozen brine at −3 to −7 °C to adjust the reaction temperature. The schematic diagram of the polymerization unit in the BR process is shown in Figure 1.

The internal leakage of the heat exchanger network in the feed section of the polymerization unit is studied in this internal leakage work. The main causes of heat exchanger internal leakage are the pressure difference and the corrosion phenomenon. In this heat exchanger network, a large pressure difference between the matter flow and the circulating water exists all the time. The circulating water system uses salt water, while the matter flow contains toxic, explosive, and flammable dangerous substances. Once this matter enters the circulating water, organic matters will become the nutrient of the microorganism, accelerating the reproduction of microorganisms and further aggravating the corrosion of the heat exchanger network. As a result, heat exchangers are susceptible to corrosion. Therefore, according to distributed control system (DCS) data, the heat exchanger with the largest pressure difference is selected for leakage case simulation.

Due to the limitation of software functions, it is impossible to simulate the internal leakage of a single heat exchanger in detail and obtain the dynamic datasets of variables directly. As a consequence, we propose an internal leakage mechanism model, in which the leakage point of each heat exchanger is assumed to be one, and the internal leakage of each heat exchanger is simulated with two heat exchangers by adding mixers and splitters. The quantity of internal leakage is represented by the flow rate from the splitters, which is marked with the red line, as shown in Figure 2.

During BR production, the refined butadiene from the butadiene workshop is preheated, and it finally leaves the polymerization unit from the bottom of the last kettle. In Figure 2, in order to simulate the matter leakage from the high-pressure side to the low-pressure side, the flow rate is changed by adjusting the parameters of splitters. In dynamic simulation, the flow can be adjusted by changing the valve opening. By adjusting the leakage rate from 0 to 20%, the corresponding temperature, pressure, flow rate, and other process data can be obtained.

For the BR process, the DYN-EW-QRA method consists of three parts: (1) internal leakage simulation; (2) early warning mechanism model construction; and (3) risk assessment. Its framework is shown in Figure 3.

The steady-state simulation is carried out to build a mechanism model based on the industrial parameters of the process. Then, the dynamic simulation is conducted based on the steady-state simulation, in which a control scheme is established to ensure the normal operation. Since the dynamic simulation takes into account the cumulative terms in the balance equation, the pressure magnitude of each module such as valves and mixers with retention capability will be adjusted. It also simulates leakage by adding disturbance to obtain time series process data. The VBA and ASM language (ASM dedicated to Aspen software) are used to simulate different internal leakage disturbance scenarios.

In the early warning stage, the normalized dynamic data are divided into a training set and a test set to carry out network prediction. LSTM is used to deeply analyze the influence of disturbance variables on the process. In order to obtain the optimal prediction results, the network hyper parameters need to be adjusted by the orthogonal parameter method. In the safety analysis stage, the range and destructive power of vapor cloud explosion hazards such as jet fire, flash fire, and explosive overpressure are calculated through the risk model with environmental parameters and variables.

3. Heat Exchanger Network Simulation

3.1. Dynamic Simulation

Steady-state simulation is the basis for the construction of the dynamic mechanism model, which guarantees the normal operation of the process and provides data for the dynamic simulation. In this section, controllers are added on the basis of steady-state simulation to facilitate the simulation of internal leakage accurately. The parameter disturbance is introduced to observe the variation of variables over time. Dynamic simulation provides the training dataset for the LSTM and also data proof for the prediction accuracy of the LSTM. The information of all controllers is listed in

Table 1. Controller parameters.

Controller Name	Tag Number	Controlled Variable	Manipulated Variable	Initial Value (kmol/h)
Flow controller	FC1	B1 cold stream feed flow	V1 opening	18,143.7
Flow controller	FC2	Feed flow	V2 opening	10,432.6
Flow controller	FC3	B4 hot stream feed flow	V3 opening	45,359.2

. The dynamic simulation flow chart is shown in Figure 4.

In Figure 4, to facilitate dynamic datasets, three independent controllers (FC1, FC2, FC3) are added to streams S1, S9, and S25, respectively. These controllers play a pivotal role in simulating an internal leakage state when the data deviates. The simulation results of main streams are listed in

Table 2. Simulation results of the main streams.

Stream	Temperature (°C)			Pressure (kPa)			Mass Flows (tons/h)
	Simulation Value	Reality Value	Relative Error %	Simulation Value	Reality Value	Relative Error %	Simulation Value	Reality Value	Relative Error %
S1	20	20	0.00	218	218	0.00	20	20	0.00
S9	41	41	0.00	520	520	0.00	11.5	11.5	0.00
S8	27.2	27	0.74	218	218	0.00	20	20	0.00
S18	21.1	21	0.48	519.6	520	−0.08	11.5	11.5	0.00
S24	78	78	0.00	519.3	520	−0.13	11.5	11.5	0.00
S25	80	80	0.00	218	218	0.00	50	50	0.00
S32	54.1	54	0.18	218	213	0.00	50	50	0.00

and

Table 3. Simulation results of heat exchangers.

HeatX	Inlet Temperature (°C)		Outlet Temperature (°C)		Heat Duty(kW)
	Hot	Cold	Hot	Cold
B1	30	20	21	24	66.033
B2	41	23	30	28	84.062
B3	62	21	54	47	411.543
B4	80	47	62	78	973.238

.

It can be seen from Table 2 that the simulated values of temperature, pressure, and mass flow for each stream have low relative error with the actual industrial data provided by the chemical enterprise. Table 3 illustrates the temperature at the inlet and outlet point of each heat exchanger along with heat duty under steady conditions.

3.2. Leakage Scenarios Simulation

To simulate the leakage scenarios, the Fsplit module is added to control the amount of butadiene entering the circulating water. In the corresponding pipeline, the dynamic changes of the temperature, pressure, flow rate of the matter, and circulating water can be obtained after the disturbance is added. By the exhaustive method, five different leakage scenarios are simulated as listed in

Table 4. Five scenarios list.

Cases Number	Description
(1)	B1 only leaks
(2)	B2 only leaks
(3)	B1 leaks first B2 then leaks
(4)	B2 leaks first B1 then leaks
(5)	B1and B2 leaks at the same time

.

For the sake of brevity, the third leaking scenario is selected to show the changes in flow, temperature, and pressure caused by internal leakage as shown in Figure 5, Figure 6 and Figure 7. In Figure 5a, the opening of the control valve V4 changes from 0 to 20% in 3 to 4 h, which means that the leakage disturbance reaches 20% of the total matter. Then, the maximum internal leakage flow of S15 reaches 584.83 kg/h and keeps stable thereafter. Similarly, the 20% internal leakage disturbance of preheating heat exchanger (B2) is added at 4 h and ends at 5 h in Figure 5b. The maximum internal leakage flow of S34 reaches 120.44 kg/h and keeps stable later. These results prove that the controllers can restore the process to a stable state when confronted with disturbance from the perspective of flow change.

In Figure 6a, when the opening of valve V4 is gradually increasing at 3 h through 4 h, part of the hot matter flows into the circulating water at lower temperature. Since the circulating water from S1 receives less heat exchange in this case, it enters B2 at a lower temperature. Therefore, compared with the stable state in the first three hours, the temperature of circulating water S8 flowing out of B2 shows a short drop at 3 h. Then, with the increase of flow rate of S15, the temperature of S8 gradually increased, and it finally stabilizes at 32.9 °C. In Figure 6b, due to the increase of flow rate of S15, the hot flow rate decreases while the cold flow rate does not change, so the temperature of S18 also drops briefly and then becomes stable at 23.5 °C again under the regulation of controllers. Similarly, when the internal leakage of B3 and B4 occurs at 4 h through 5 h, the temperature of S24 rises continuously and then keeps stable at 79.2 °C, as shown in Figure 6c. In Figure 6d, when the opening of valve V5 changes from 0 to 20% in 4 to 5 h, the internal leakage flow represented by S34 increases, and the cold flow from S18 hotter than the original temperature enters B3, leading to a gradual rise of the temperature of S32. Due to the hysteresis of FC1 and FC2, under the action of FC3, the temperature of S32 drops briefly and then stabilizes at 56.3 °C within 8 h.

In Figure 7, when disturbance is added at 3 h, the internal leakage disturbance is minor. Since the matter flow is increased under the control of FC2, the pressure rises slightly. At 5 h, the internal leakage disturbance of B1 and B2 reaches 20%, so the inside pressure reaches the peak of 4.8014 bar. Subsequently, as a large amount of matter flows into the circulating water, the pressure inside the tube drops continuously and the maximum pressure reaches 4.8014 bar; then, it remains stable at 4.7980 bar under the control of FC2 ultimately.

4. Early Warning Model Construction

4.1. Model Construction and Optimization of LSTM

The dynamic dataset is obtained by simulating five internal leakage scenarios. For LSTM modeling, we simulate multiple disturbances as the training set and a single disturbance as a test set under the same conditions to verify the predictive result of the network. Furthermore, each set of 2000 continuous time samples is fed into the LSTM input gate in batches, constituting 22 sets in total. The training set and test set include 2000 and 4000 samples, respectively.

The difference between LSTM and other neural networks is the temporal link between LSTM units in the hidden layer. A typical LSTM cell usually consists of an input gate, an output gate, a forget gate, a unit input, a cell state, and peepholes. The cell state of the LSTM changes over time to preserve long-term memory. LSTM uses storage cells instead of simple neurons to solve the problem of gradient explosion and gradient vanishing. Figure 8 presents the simplified LSTM units [27]. It consists of several layers and point-by-point operations that provide information about the state of the LSTM unit and maintain long-term memory through the network and input.

Through the introduction of the memory state and three gates, LSTM is able to find what historical data is worth remembering and what data should be forgotten. In general, LSTM can represent or simulate human behavior, logical development, and the cognitive process of neural tissue more realistically.

To further enhance the accuracy of network prediction, the number of network layer nodes, activation functions, and batch size are optimized by an orthogonal experiment. The primary and secondary relationship can be distinguished from many influential factors to quickly affirm the optimal network parameters. Then, the required number of the experiment is reduced, and the efficiency is intensified. The experimental schemes of the orthogonal experiment are listed in

Table 5. The schemes of the orthogonal experiment.

Scheme Number	Horizontal Combination	Optimization Parameter
		Number of Network Layer Nodes	Activation Function	Batch Size
A	A1 B1 C1	5	sigmoid	50
B	A1 B2 C2	5	relu	100
C	A1 B3 C3	5	tanh	200
D	A2 B3 C1	10	tanh	50
E	A2 B1 C2	10	sigmoid	100
F	A2 B2 C3	10	relu	200
G	A3 B2 C1	30	relu	50
H	A3 B3 C2	30	tanh	100
I	A3 B1 C3	30	sigmoid	200

.

The activation function is the mechanism by which neurons process and transmit information through the neural network. The activation function helps determine whether we need to activate neurons and the intensity of those neurons signals. Different activation functions may lead to differences in network classification or prediction performance, so we compare three common activation functions including Sigmoid, Relu, and Tanh for training neural networks. For the number of layers in the network, the larger it is, the stronger the learning and expression abilities of the network to data. However, a larger layer node number does not definitely mean a better performance, because on the one hand, it will result in a training cost, and on the other hand, it will result in unnecessary fitting brought by a more complex network structure. So, when setting up the network model, we need to choose an appropriate size of the dataset and network layer node number. Batch Size refers to the number of samples selected once training. Its size affects the optimization degree and speed of the model. A large Batch Size means accurate gradients. However, it may cause enormously different gradients after a certain level, therefore making it difficult to use a global learning rate. An appropriate Batch Size can reduce the number of iterations for network training, improve the training speed, and make the gradient descent direction more accurate. In order to improve the accuracy of the internal leakage prediction, nine calibration schemes of network parameters were constructed through orthogonal tests. The ultimate purpose is to make the predicted result slightly greater than the actual leakage quantity for an effective subsequent risk assessment as well as prediction accuracy.

In Table 5, the normal multinomial regressive experimental design of three factors and three levels are adopted. The optimization parameters are shown as follows:

• The number of network layer nodes are 10, 30, and 50.
• Activation functions are sigmoid, relu, and tanh.
• Batch sizes are respectively 50, 100, and 200.

The optimized prediction effect is intuitively shown in Figure 9. It can be seen that the predictive effect of Figure 9a,d,g is better than others. However, the predicted value in Figure 9d is smaller than the observed value overall. So, it may underestimate the harm of risk consequences, resulting in inadequate protection measures preparation and further unnecessary losses. Conversely, the prediction accuracy of Figure 9a is more excellent than that of Figure 9g, so the orthogonal experiment of scheme (a) is the most effective. Therefore, five layer nodes, sigmoid activation function, and a Batch Size of 50 are selected as optimal LSTM network parameters in this paper.

4.2. Monitoring Process

On the basis of the analysis in Section 4.1, adam is select as the optimizer. The network adopts the sigmoid activation function, and its recurrent core adopts the hard-sigmoid function for prediction with the mean absolute error as the loss function. The dropout layer is added to prevent the training process from over fitting.

An early warning analysis model is the core of early warning implementation. To improve the timeliness and accuracy, the early warning indicator is selected among data-driven statistic methods. Principal Component Analysis (PCA) is a data-driven process monitoring method that is widely applied in engineering practice. The assumption that the data obeys Gaussian distribution causes its monitoring effect to be poor and the extraction of nonlinear features to be unrobust. In a consequence, Johnson Gaussian transformation is performed on process data to make it obey Gaussian distribution. The implementation steps are similar to PCA. The advantage of this combination is that it well captures the non-Gaussian changes of the process more robustly. The process is shown in Figure 10.

The monitoring steps are divided into two stages. Stage 1 is offline modeling with following steps: (1) standardize the normal historical process data to make the data have zero mean and unit variance; (2) perform Johnson transformation on the standardized data to make them obey normal distribution; (3) extract the feature vectors and construct the feature space; (4) project the normal data into the feature space and calculate the thresholds of T². Stage 2 is online monitoring with following steps: (1) standardize the online process data to make it have zero mean and unit variance; (2) perform Johnson transformation on the standardized data to make the data obey normal distribution; (3) project the data into the feature space and calculate the T² statistics; (4) judge whether the online data of T² exceed its threshold value. If the threshold is exceeded, go to step 5; otherwise, go to step 6; (5) internal leakage early warning; (6) go back to step 1 and continue to monitor the next set of data.

As an early warning indicator, the T² statistic represents the distance between the observed value in the characteristic subspace and the data center. Its calculation formula is given as Equation (1):

$${\mathit{T}}^2={\mathit{S}}_{\mathit{i}}^{\mathit{T}}{\Lambda }_{\mathit{a}}^{-1}{\mathit{S}}_{\mathit{i}}$$

(1)

where Λa represents the diagonal matrix of the diagonal elements corresponding to the previous a eigenvectors. The threshold formula of T² is as shown in Equation (2):

$${\mathit{T}}_{\mathit{c},\alpha }^2=\frac{\mathit{l}({\mathit{n}}^2-1)}{\mathit{n}\left(\mathit{n}-1\right)}{\mathit{F}}_{\alpha }\left(\mathit{l},\mathit{n}-1\right)$$

(2)

where Fα(l, n-l) represents the F distribution with the confidence level of α, the degree of freedom of l, and n − 1, n is the number of samples in the training set, and l is the number of selected eigenvectors.

The predicted results of Cases (1)–(4) listed in Table 4 are shown in Figure 11. Shown as the dotted red line in Figure 11, the indicator threshold value four cases are 3.54, 3.36, and 2.56 kg/h, respectively. Prediction curves of Cases (1), (3), and (4) are tightly coincident with observation line. It can be seen that the indicator can early warn of internal leakage when the leakage flow is relatively small. It should be noted that the amount of test data in Case (2) is only 1500. On account of the bad effect compared with the other three cases, we retrain the network with the first 500 data of the test set as the training set.

5. Risk Assessment

Butadiene and hexyl hydride are flammable and can form an explosive mixture when mixed with air. Once they are exposed to high temperature, spark, or oxidants, combustion explosions are bound to happen. Furthermore, butadiene and hexyl hydride vapor is heavier than air, so they can diffuse far away from the surface of the earth, leading to a rekindling of fire. Therefore, it is necessary to conduct safety analysis with regard to jet fire, flash fire, and vapor cloud explosion. For Cases (1)–(4) caused by internal leakage, the quantity of leakage is roughly divided into three categories. Parameters related to leakage flow are listed in

Table 6. Parameters of leaked flow.

Case	Leak Flow (kg/h)	Leakage Type	Atmospheric Stability	Wind Speed (m/s)
1	120.4	Internal leakage	D	5
2	584.8
3	705.2

.

5.1. Jet Fire

According to the different flow pattern of flammable matter from the pipeline, there are usually three types of jet form: floating plume, floating jet, and pure jet. Although the diffusion modes of three jets are different, their damage models are the same. The Hawthorne flame length can be calculated as:

$$\mathit{L}=(5.3\mathit{d}/{\mathit{C}}_{\mathit{t}}){\{({\mathit{T}}_{\mathit{f}}/\alpha {\mathit{T}}_{\mathit{n}})[{\mathit{C}}_{\mathit{t}}+(1-{\mathit{C}}_{\mathit{t}})\mathit{M}\mathit{s}/\mathit{M}\mathit{a}]\}}^{1/2}$$

(3)

where L is the flame length (m), d is the injection diameter (m); C_t is the molar concentration of fuel air mixture (mol/L); T_f is the flame temperature (K); α is the ratio of the number of moles of reactants in the mixture to the number of moles of combustion products; T_n is the ambient temperature (K); M_S is the molecular weight of fuel gas (g/mol); and Ma is the air quality (g/mol).

The flame energy of a jet fire radiates from its surface to the surroundings. If its surface and side area can be approximately cylindrical, and the maximum width of 2Y_smax are used as the diameter of the bottom surface, the side area is 2πY_smaxL. The heat flux on the flame surface can be calculated by the following formula:

$${\mathit{q}}_0={\mathit{m}}_{\mathit{f}}\mathit{H}\mathit{f}/(2\pi {\mathit{Y}}_{\mathit{s}\max }\mathit{L})$$

(4)

where f is the thermal radiation coefficient, which can be 0.15; M_f is the combustion rate (kg/s). The following formula is used to solve the leakage model:

$${\mathit{m}}_{\mathit{f}}={\mathit{c}}_0{\mathit{A}}_{\mathit{h}}{[\mathit{P}{\rho }_0\gamma {(\frac{2}{\gamma +1})}^{\frac{\gamma +1}{\gamma -1}}]}^{1/2}$$

(5)

where M_f is the gas mass leakage rate (kg/s); C₀ is the leakage coefficient; A_h is the area of the leakage hole (m²); P is absolute pressure (Pa); ρ₀ is the gas density in the tank (kg/m³); and γ is the adiabatic coefficient of the gas.

In Figure 12, jet fires caused by three different leakage scenarios are simulated based on leakage flowrates of 120.4, 584.8, and 705.2 kg/h, respectively. The hazard radius of jet fire is proportional to the leakage amount, which is 3.5, 5.5, and 9 m. As the distance increases, the heat flux gradually decreases. The heat flux closest to the center of the flame is 37.5 kW/h, the outermost is 4 kW/h, and the middle is 12.5 kW/h. Obviously, different heat fluxes have different destructive forces.

Table 7. The hazardous level of radiant heat flux.

Radiant Heat Flux (kW/m2)	Damage to Equipment	Damage to People
37.5	All operating equipment is damaged	10 s, 1% death 1 min, 100% death
25	Minimum energy for wood burning under flameless, long-term radiation	10 s, major casualties 1 min, 100% death
12.5	Minimum energy for wood burning and plastic melting in the presence of flame	10 s, first degree burn 1 min, 1% death
4.0		20 s or more, pain

shows that people will suffer from minor burns inside the blue line and outside the green line. Wood will burn and people will get first-degree burns inside the green line and outside the red line. All operating equipment in the red line will be damaged, and people will die within 1 min.

5.2. Flash Fire

Flash fire is a non-explosive combustion process that occurs when combustible gas or vapor leaks into the air and is ignited after mixing with air. There is no self-accelerating condition in the flame propagation process; that is, the leakage source does not produce turbulence, and there is no local constraint condition in the space. The consequences are mainly from thermal radiation and direct contact with flames, causing large-scale fires to burn property and do harm to people. The hazard range of flash fire under different leakage quantity of 120.4, 584.8, and 705.2 kg/h is simulated and shown in Figure 13.

Combining the conservation relationship of mass, momentum, and energy with some empirical facts, the relationship of flame width W with time is given:

$$\mathit{W}=2{[{\mathit{R}}^2-{(\mathit{R}-\mathit{S}\mathit{t})}^2]}^{1/2}$$

(6)

where R is half width.

In Figure 13, the blue line represents the lower limit of flash concentration, and the red line represents the upper limit of flash concentration. The radius of the gas concentration exceeding the upper limit of the flash fire is proportional to the amount of leakage, which is 2, 3, and 5 m respectively. The area between the upper and lower limits of the flash fire concentration can also cause a certain degree of harm. Gas concentrations beyond lower limits of the flash fire do not cause danger.

5.3. Vapor Cloud Explosion Overpressure

After leakage occurs, the following four conditions should be met for the occurrence of the vapor cloud explosion with destructive overpressure: (1) The leakage material must be flammable with appropriate temperature and pressure conditions. (2) A sufficiently large cloud must be formed in the diffusion stage prior to ignition. (3) Only when a sufficient number of clouds lie within the explosive limit of the material can significant overpressure be generated. (4) Turbulence exists in confined space propagation or initial ignition energy exists due to source excitation. The TNT equivalent method is used to calculate the mathematical model of vaper cloud explosion consequence, which is a method to connect the explosive energy generated by accidental explosion with the quantitative TNT. It can be determined by the following formulas:

$$\mathit{F}=1-\exp [\frac{-{\mathit{C}}_{\mathit{p}}\Delta \mathit{T}}{\mathit{L}}]$$

(7)

$${\mathit{W}}_{\mathit{f}}=2\times \mathit{F}\times \mathit{W}$$

(8)

$${\mathit{W}}_{\mathit{T}\mathit{N}\mathit{T}}={\alpha }_{\mathit{e}}{\mathit{W}}_{\mathit{f}}{\mathit{H}}_{\mathit{f}}/{\mathit{H}}_{\mathit{T}\mathit{N}\mathit{T}}$$

(9)

$$\overline{\mathit{R}}=\mathit{R}/{\mathrm{W}}_{\mathrm{TNT}}^{1/3}$$

(10)

where F is the evaporation coefficient; C_p is the average specific heat of fuel (kJ/(kg∙K)); ΔT is the temperature difference between the temperature in the container under ambient pressure and the temperature of the boiling point (K); L is the latent heat of vaporization (kJ/kg); W_f is cloud fuel mass (kg); W is the amount of fuel leaked (kg); W_TNT is equivalent to TNT (kg); H_f is the combustion heat of fuel (MJ/kg); H_TNT is the detonation heat of TNT (MJ/kg); α_e is the equivalent coefficient of TNT, R is the actual distance from the explosion point (m), and R is the dimensionless distance.

Overpressure refers to the pressure difference between the explosion wave pressure and atmospheric pressure. It is divided into positive pressure and negative pressure. Positive pressure is harmful to the human body and buildings, so only positive pressure is considered here. The relationship between overpressure and destructive force is shown in

Table 8. The harm of overpressure to people and buildings.

Overpressure (P/bar)	Damage to Buildings	Damage to Body
>0.75	House collapse	Instant death
0.45–0.75	Destroy all doors and windows	Serious injury
0.25–0.45	Destroy most doors and windows	Moderate injury
0.1–0.25	Destroy part of doors and windows	Slight injury
<0.01	Basically no damage	Safety

. When a vapor cloud explosion occurs, it will cause local turbulence and accelerate the propagation speed. The relationship between overpressure and distance is shown in Figure 14.

The maximum leakage is selected as 705.2 kg/h, and the leakage time is one hour. Combining Figure 14 and Table 8, the overpressure range between the red lines can be obtained as 0.25 to 0.75 bar. The overpressure is inversely proportional to distance in general. Therefore, the overpressure is greater than 0.75 bar within a distance of 0–12 m, which will cause severe damage to people and property. The overpressure decreases gradually in the range of 12–23 m, but it is still destructive. Over 23 m, the overpressure damage tends to be in a safe range.

5.4. Protective Measures

According to the simulated risk consequences of different disturbances, the reasonable protection measures are proposed for the internal leakage. To provide reliable security protection information for the chemical industry, the safety analysis report of the internal leakage of DYN-EW-QRA is listed in

Table 9. The safety analysis report of DYN-EW-QRA.

Accident	Internal Leakage
Consequence	a. Human death within 12 m, wall and roof collapsed within 23 m. b. Combustion within 9 m, equipment structure damage, metal plate distortion. c. Heat flux in the hazard range of 9 m. d. 0 to 23 m exploded with open fire
Recommended measures	a. Add fire-resistant coating for the important equipment within 30 m from the leak point. b. Install the electrostatic discharge devices within 30 m. c. Install the flammable and toxic gas detector at the flange interface. d. Add necessary interlocking and alarm devices in the process

. Table 9 shows that DYN-EW-QRA provides a detailed analysis for the consequences of internal leakage, including casualties, damage to houses, equipment, and explosion limits in meter. DYN-EW-QRA also shows dynamic hazard consequences based on the deviation of process parameters. Through comprehensive analysis of risk consequences and frequency, the recommended measures given by DYN-EW-QRA are highly quantitative and reliable.

6. Conclusions

To realize timely and quantitative early warning and risk analysis of internal leakage, a novel intelligent early warning and risk assessment strategy is proposed based on the internal leakage mechanism model. The simulated values of temperature, pressure, and mass flow for each stream are calculated by steady-state simulation with results very close to their real industrial values. The early warning analysis model is established by LSTM, in which PCA is selected for the early warning indicator. The QRA approach is applied to conduct risk assessment, whose results show that the hazard radius of jet fire is proportional to the leakage amount, which is 3.5, 5.5, and 9 m. The radius of the gas concentration exceeding the upper limit of the flash fire is proportional to the amount of leakage, which is 2, 3, and 5 m. The overpressure will cause severe damage to people and property within 23 m. The vapor cloud of organic matters caused by internal leakage will explode if it encounters open fires between 15 and 30m. According to the analysis of dynamic results, the important equipment within 30 m from the explosion center should be coated with fire-retardant coating. Simultaneously, the electrostatic discharge devices should be installed within 30 m of the factory equipment group. At the appropriate location in the factory, the detectors and alarm devices need to be increased.

The dynamically updating risk realized in this paper indicates more effective and real-time warning information. However, due to the limitation of the process, the applicability and reliability of the method need to be verified in a larger process in the future.