The Ignalina Nuclear Power Plant (NPP) contains two RBMK-1500 reactors. ... of the RBMK reactor design series (actually the only two of this type that were ...
Proceedings of the 8th International Conference on Probabilistic Safety Assessment and Management May 14-18, 2006, New Orleans, Louisiana, USA
PSAM-0471 PROBABILISTIC MODELING OF AIRCRAFT CRASH AND IMPACT ON IGNALINA NPP CONSIDERING UNCERTAINTY Robertas Alzbutas/Lithuanian Energy Institute
Gintautas Dundulis/Lithuanian Energy Institute
Juozas Augutis/Lithuanian Energy Institute
Eugenijus Ušpuras/Lithuanian Energy Institute
ABSTRACT The objective of considered work is to present the approach and analysis of the probabilistic modeling of the aircraft crash as well as results of the probabilistic modeling of aircraft crash impact on the main building of Ignalina nuclear power plant. In order to estimate the probability of this critical event occurrence near the Ignalina nuclear power plant in Lithuania, the model of aircraft crash frequency is developed and analyzed. The real statistical data are used. The impact of aircraft crash on the building of the plant is considered. The simulation and probabilistic estimation of this event occurrence, the uncertainty and sensitivity analysis and the investigated modeling approach as well as used software application possibilities are discussed. And finally, the probabilistic analysis of the most common aircraft used in Lithuania, i.e. two engine commercial airliner crash impact on the plant is considered. Key words: probabilistic modeling, aircraft crash frequency, aircraft crash impact, uncertainty and sensitivity analysis. INTRODUCTION The Ignalina Nuclear Power Plant (NPP) contains two RBMK-1500 reactors. This is the most advanced version of the RBMK reactor design series (actually the only two of this type that were built). "RBMK" is a Russian acronym for "Channelized Large Power Reactor". Compared to the Chernobyl NPP, the Ignalina NPP is more powerful (1500 MW versus 1000 MW), and is provided with an improved ALS (Accident Localization System). In the most other respects, the plants are quite similar to their predecessors. They have two cooling loops, a direct cycle, fuel clusters are loaded into individual channels rather than a single pressure vessel, the neutron spectrum is thermalized by a massive graphite moderator block. The plant can be refueled on line and uses slightly enriched nuclear fuel. The Ignalina NPP as a lot of others NPP’s made significant safety improvements after Chernobyl accident. Efforts to perform additional deterministic and probabilistic analysis and improve safety of NPP’s were also accelerated in Lithuania. The results of Probabilistic Safety Assessment (PSA) of the Ignalina NPP clearly indicated the need of continued attention for the risk analysis and safety improvements. More advanced models were created and new approaches for analysis were used. The sensitivity and uncertainty analysis related to aircraft crash was considered and related to the modeling of aircraft crash. The modeling of aircraft crash occurrence and impact on Ignalina NPP is significant in order to evaluate risk and prevent potential failure of the systems important to safety. In order to estimate aircraft crashes influence to safety of Ignalina NPP a general concept for analyzing this event is developed and the new simulation approach of event occurrence is proposed. The probabilities of such event occurrence and impact on the plant are determined using the analyzed statistical data, mathematical models and probabilistic simulations.
Copyright © 2006 by ASME
DATA FOR PROBABILISTIC MODEL OF THE AIRCRAFT CRASH Within the analysis of an aircraft or other flying device (helicopter, balloon etc.) crash, as a rule, unintended aircraft accidents are examined, but terrorist actions or other not ordinary human activities are excluded. However the general probabilistic analysis of the aircraft crash impact on the plant, including the analysis of consequences in case of the terrorist actions, was performed separately and results are presented in section 5. Statistically, the analyzed frequency of aircraft crashes depends on the intensity of flights near the target object, the technical condition of the aircrafts, the experience of crew, the navigational elevated aids (radio beacons), the meteorological conditions and other factors. The initial data for estimating the aircraft crash probability includes: • Distance of the NPP from civil or military airports; • Arrangement of air transport corridors and their distance to the NPP; • Intensity of flights in air transport corridors; • Distribution of flying aircrafts by types; • Generalized world statistics for aircraft crashes by weight and types; • Statistical data on serious incidents at the state. An aircraft crash into the NPP can cause severe damages to the plant. Even if the catastrophe takes place near the NPP it is usually followed by explosion and fire and can significantly influence safe operation of the plant. There is no exact deterministic way to model a problem of this type and some of the factors cannot be easily quantified . The frequency of aircraft crashes on studied objects depends on the amount of aircraft flying close to these objects. More difficult it is to evaluate crash factors such as aircraft technical state, crew experience, reliability of navigation techniques, meteorological, and other conditions. Because of high uncertainty of the accident initiating factors, aircraft crash modeling usually is based or leads to the conservative assumptions or the uncertainty analysis of the results . At first, the simple initial probabilistic model of the aircraft crash was investigated in order to analyze data and estimate probability of this critical event occurrence near the Ignalina NPP. As there are no airports in the surroundings of the Ignalina NPP, the aircraft strike probability per year can be estimated as: P = Pl⋅Nc⋅A⋅F,
Where P is aircraft strike probability per year, Pl is aircraft strike frequency per flight kilometer, Nc is flight number per year, A is target area, F is function of deviation per flight kilometer from initial flight route during an accident. As all flights below 6000 meters are prohibited in the area of 10 km radius around the Ignalina NPP, the main aircraft accident sources were considered to be the closest international air corridors and . According to the data of Lithuanian flight control centre flight intensity in the Lithuanian corridor was approximately 20000 flights per year. As flight intensity data in the Belarussian-Latvian corridor were not exactly known the maximum number of flights was conservatively estimated to be 50000 per year. To maintain conservative approach air corridors of 50 km width were modeled to pass in 10 kilometers distance from the Ignalina NPP. Table 1. Relatively new and old plane crash frequencies Aviation Catastrophe Frequency / year For 1 flying hour For 1 flying km
Relatively new planes Weight Weight up to 5700 kg over 5700 kg 2.1·10-5 9.0·10-7 8.4·10-8 1.2·10-9
Relatively old planes Weight Weight up to 5700 kg over 5700 kg 2.5·10-5 1.0·10-6 1.0·10-7 1.3·10-9
Aircraft strike frequency per flight kilometer was obtained from Department of Civil Aviation and is shown in Table 1. Flights are carried out by both relatively new and relatively old (soviet type) aircrafts, which are different in weight. Therefore, catastrophe intensity was taken separately for different aircrafts with weight less than W and with weight W and higher (used limiting weigh W = 5700 kg).
Copyright © 2006 by ASME
PROBABILISTIC MODEL OF THE AIRCRAFT CRASH For aircraft crash probability estimation, when flying route distance from the object territory is s, the following mathematical model was improved:
P ( s ) = Pl ⋅ N c ⋅ A ⋅
g − gs ⋅e , 2
Were g - a constant dependent on type of planes and describing likelihood of close crash (for transport planes g = 1.00, for military planes g = 0.63 and for passenger planes g = 0.23). This means that probability of the plane to deviate 10 kilometers aside its initial trajectory during the accident is 10 times less than to crash on it. Considering that 10 km radius zone around Ignalina NPP at height lower than 5950 m are non-flying zone, conservatively it was assumed that all flights from 50 km radius zone, take place at the periphery of 10 km radius zone. As the civil planes compose major number of flights, value of coefficient g was taken 0.23, which at the same time is the most conservative from the above mentioned. This value was also used in some references, e.g. 0.
Route of flight X, km 50
x y s
Figure 1. The model of aircraft crash calculation on territory of NPP with radius r Plane crash probability onto r radius zone Z (Figure 1) with condition if plane loses control and start to fall in distance s from the zone centre is expressed by:
Pr ( s ) = ∫∫ p ( s ) dsd ϕ ≅ Z
r2 ⋅ g2 ⋅e 2
π g − gs = ⋅ e , that g 2
− g ⋅s
s ⋅ p(s)dsd ϕ = 1
If plane flies 50-km radius route that makes tangent for 10-km radius zone, then the distance to the analyzed zone centre is expressed by:
s = x 2 + y 2 , y = 10 and x∈(-50, 50)
From the flying route common plane crash onto r radius zone is:
N c Pl r 2 g 2 2
x 2 + 100
Using this formula we can evaluate crash probabilities on different radius territory that corresponds to the area of Ignalina NPP buildings or reactor.
Copyright © 2006 by ASME
In order to determine the dependencies of results on the radius r the calculations with different r values and plane crash frequencies (Figure 2) were performed. 8.0E-06
P/year New type
Figure 2. Light (a) and heavy (b) plane crash frequencies
Obtained estimations of heavy planes crash probabilities are lower than those obtained by probabilistic analysis for most NPPs in USA, where the number of flights is higher. SIMULATION AND UNCERTAINTY ANALYSIS APPROACH The approach suggested for uncertainty and sensitivity analysis, and illustrated in this paper by an application of SUSA (Software System for Uncertainty and Sensitivity Analyses), is based on well-established concepts and tools from probability calculus and statistics. It requires identification of the potentially important contributors to the uncertainty of the model results and the quantification of the respective state of knowledge by subjective probability distributions . Such a distribution expresses how well an uncertain parameter of the model application is known. The sensitivity analysis is used to identify the main contributors to the possible variations of results. Sensitivity analysis is performed in connection with uncertainty analysis in order to see the combined influence of all the potentially important uncertainties on the result. Sampling and Uncertainty Measures The quantitative uncertainty analysis results can be expressed as quantiles (as example 5% and 95%) of the result distribution. They could be obtained easily if result distribution was known. In practice these quantiles are estimated using parameters subjective probability distributions and Monte Carlo simulations. The quantitative uncertainties of model parameters can be expressed by the parameters distribution with mean and standard deviation values. Standard deviation in normal distribution case can be assumed as a value, which is three times less than the interval between maximum and minimum values. If distribution is untruncated, then the probability, that parameter value belongs to this interval, is 0.866. Otherwise, if distribution is truncated, then the probability is equal to one. The simulations of aircraft crash are performed for each sample set (all random parameters varied simultaneously). For this purpose, the special software was created and integrated into the SUSA, so that the desired quantiles can be estimated from this sample by standard statistical techniques. In addition, the possible impact of the sampling error is considered. Usually this can be done by the computing (u, v) statistical tolerance limits . Where v is the confident level that maximum model result will not be exceeded with the probability u (or u% quantile, which reflects the amount of combined influence of all quantified uncertainties) of the corresponding output distribution, which is to be compared to the acceptance criterion. According to the classical statistical approach the confidence statement quantifies the possible influence of the fact that only a limited (frequently small) number of model runs have been performed. For example, according Wilks’ formula , 93 runs are sufficient to have two sided (0.95, 0.95) statistical tolerance limits. The required number n1 of runs for one-sided tolerance limits and correspondingly the number n2 for two-sided statistical tolerance intervals can be expressed as following: n1 ≥ ln(1 − v) / ln(u ) ; n2 ≥ (ln(1 − v) − ln((n2 / u ) + 1 − n2 )) / ln(u ) .
The minimum number of model runs needed for these limits is independent of the number of uncertain quantities taken into account and depends only on the two probabilities u and v given above. The amount of runs is a result from nonparametric statistics. Its advantage is that this amount is completely independent to the number of uncertainties taken into account and does not assume any particular type of underlying distribution. The distribution event does not need to be continuous .
Copyright © 2006 by ASME
Sensitivity Measures Results from applications of models are subject to uncertainty. Usually an uncertainties analysis can provide a statement about the combined influence of potentially important uncertainties on the results. Often more important, it provides quantitative sensitivity statements that rank the uncertainties with respect to their contribution to the model output uncertainty . In order to rank uncertainties according to their contribution to model output uncertainty, standardized regression coefficients (SRCs) are chosen here from the many measures available. They are capable of indicating the direction of the contribution (negative means inverse proportion). Additionally, the correlation ratios are computed. The correlation ratio is the square root of the quotient of the variance of the conditional mean value of the model result (conditioned on the uncertain parameter) divided by the total variance of the model result due to all uncertainties taken into account. So it serves as a measure, how one model uncertainty was quantified through a set of alternative model formulations. The correlation ratio quantifies degrees of parameters and results relationship. A standardized regression coefficient is supposed to tell by how many standard deviations the model result will change if the uncertain parameter is changed by one standard deviation. How well this is achieved in practice depends on the degree of linearity between the model result and the uncertain parameter. In case the numbers of uncertainties are large and sample size is small, spurious correlations can play a non-negligible role. The effect of spurious correlations on sensitivity measures may be avoided if the estimates of correlation coefficients and standardized regression coefficient are compared . The options described above and chosen for this illustration of the approach are part of the software system SUSA . Aircraft Crash Sampling A part of the state of knowledge quantification is presented in Table 2. The main parameters, which may impact the calculation uncertainty, were divided into two main groups related to the initial conditions and to the model parameters. Table 2. Selection of parameters, which may impact the uncertainty of calculation results #
Initial conditions 1 Pl, aircraft impact frequency per 1 km 2 Nc, flight number per 1 year Model parameter 3 r, aircraft impact radius, km 4 g, aircraft impact deviation coefficient
Ranges Min. Max.
Reference (mean m)
Normal, truncated Normal, truncated
Normal, truncated Normal, truncated
As some initial conditions and plane crash model parameters are not exact or have different values for different plane manufacturing, there was performed uncertainty analysis. Referring Wilks formula and uncertainty methodology investigated in this paper, there were performed 100 numerical tests simulating different plane crash variations. As quantitative uncertainty measures two-sided statistical tolerance limits (with given probability u = 0.95 and confidence v = 0.95) were chosen for each of the model results of interest. These limits contain at least 95% of the combined influence of all quantified uncertainties at a classical statistical confidence level of at least 95% (resulting v = 0.9629). Two sided tolerance limits formed by sample extremes are 1.4·10-9 and 5.9·10-8. The empirical distribution function is presented together with fitted Gamma distribution, which can be used as quite good approximation of model output.
Copyright © 2006 by ASME
Table 3. Main statistical characteristics of aircraft impact model output sample Result with (0.95, 0.95) tolerance limits P aircraft impact probability per year
Standard Confidence limits deviation (5%, 95%) 5.9·10-8 1.7·10-8 1.2·10-8 Lower: 3.4·10-9 Upper: 3.8·10-8
The minimum, maximum, mean, standard deviation and other characteristics of the model output sample are presented in Table 3. The lower and upper confidence limits are 5% and 95% quantiles derived from suitable observations by classical statistics. The 95% quantile, for instance says that the appropriate result value is below this quantile value with the 95% degree of belief. Sensitivity and Uncertainty Analysis Sensitivity estimations (Figure 3) evaluated from statistic analysis describe how plane crash initial conditions and model parameters (Table 2) influence the results. From the sensitivity results it is seen that the parameter Pl has the highest influence on the results. 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
-0.35 Parameter index 0
Figure 3. Standard regression coefficients (a), empirical correlation coefficients (b)
The sensitivity factors (see Fig. 3), which are a by-product of the probabilistic analysis, tell the analyst which sources of uncertainty are contributing most to the uncertainty in the predicted performance. It is possible that sensitivity measure SRC explain too small a fraction of the variability of the model output values if coefficient of determination (the correlation coefficient between the computer model output and the corresponding output from the multiple regression model) R2 < 0.5, for instance. However, for analyzed sample the coefficient of determination R2 = 0.7. Thus, the analyst can determine which variables (model parameters) should be regulated and controlled better in order to decrease unfavorable event occurrence. Alternatively, the analyst can determine which tolerances could be relaxed without substantially affecting system risk. Summarizing, the maximal and minimal values of air crash per year probability are equal to 5.9E-08 and 1.4E09, respectively. The average value of the aircraft crash per year probability is equal to 1.7E-08, the standard deviation is equal to 1.2E-08. The Gamma distribution is used to represent the obtained results. The main results of modeling are presented in the following Figure 4. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.E+00
Probability per year 1.E-08
Figure 4. Distribution function of modeling results
Copyright © 2006 by ASME
PROBABILISTIC CONSEQUENCE ANALYSIS After the terrorist attacks in New York and Washington D. C. using civil airplanes, the evaluation of civil airplane crashes has become one of the main external event considered for current NPP’s as well as for new NPPs’ . The interceptions of many terrorists’ communications reveal that the use of commandeered commercial aircraft is still a major part of their plans for destruction. Aircraft crash or other flying objects in the territory of Ignalina NPP represents a very large danger to the plant, including the reactor. Aircraft traveling at high velocity has a huge destructive potential. The aircraft crash may destroy the roof and walls of buildings, pipelines, electric motors, cases of power supplies, power cables of electricity transmission and other elements and systems, important for safety. Therefore, the analysis of aircraft crash impact on the plant is very important from safety point of view and was performed for Ignalina NPP. The probabilistic analysis of the aircraft (the two engine commercial airliner of the most common aircraft used in Lithuania) crash impact on the plant building consequences was carried out. The software systems ProFES  and NEPTUNE  are used for probabilistic analysis of failure of the walls of the Ignalina NPP building. ProFES is a probabilistic analysis system that allows designers to perform probabilistic finite element analysis (FEA) in a 3D environment that is similar to modern deterministic FEA. In this study ProFES is coupled to the deterministic finite element code NEPTUNE. By default the NEPTUNE input can not be imported into the ProFES code. Thus, a special ProFES/NEPTUNE translator  was developed to facilitate this effort. The coupled code approach simplifies the modeling of uncertainties and evaluates the effect of these uncertainties on model and system reliability. Taking into consideration that the probabilistic analysis takes a lot of time, a simplified Finite Element (FE) model of the Ignalina NPP building was used for this analysis. The impacted wall and the adjacent walls and ceilings are included in the FE model of the Ignalina NPP building. This FE model was created using the ALGOR preprocessor. The building FE model prepared for probabilistic analysis is presented in Figure 5.
a) b) Figure 5. The finite element model for analysis of the Ignalina NPP building: a – FE model view in the ALGOR window, b - FE model view in the ProFES window
The aim of the uncertainty analysis is to identify and quantify all potential important uncertainty parameters. Ranges and subjective probability distribution describe the state of knowledge about all uncertain parameters. In probabilistic analysis, uncertainties in numerical values are modeled as random variables. The following mechanical properties and geometrical parameters important to strength of structures were used as random variables: • Mechanical properties: 1. Concrete – Young’s modulus, points of the stress-strain curve of the impacted and neighboring walls; 2. Reinforcement bar – points of the stress-strain curve of the impacted wall and neighboring walls. • Geometry data: 3. Concrete wall – thickness of the impacted wall and neighboring walls; 4. Reinforced concrete – rebar area of the impacted wall and neighboring walls.
Copyright © 2006 by ASME
The testing of material properties of reinforced concrete were carried out at the Ignalina NPP. The coefficient of variables was adopted according to the data of tested samples. The logarithmic normal distribution for mechanical properties and geometry parameters was used in this analysis. The peak points of load curves are modeled as random variables. The normal distribution for load points was used in this analysis. During an airplane crash, the dynamic impact loading is uncertain. Therefore it is important to estimate the dependence of failure probability of the building due to the uncertainty in the loading. The RS/MCS method was used for the determination of such a relation expressed by the probability-loading function. The RS/MCS method fits a response surface about the mean of the random variable and then computes statistics of response variables and probabilities for each limit state using MCS. First the RS method was used to determine the probability function for probabilistic analysis of failure of the walls of Ignalina NPP building. Then the MC method was used to study the probability of failure of the building walls. In this analysis the probability function was used as internal function in ProFES to determine the failure probability. The results of the probabilistic consequence analysis are presented in Figure 6. 1 0.9 0.8
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
Figure 6. Probability of failure in the impacted concrete wall due to the airplane crash
According to this result the relation between probability of wall failure and applied load was defined. The failure (concrete element reaches ultimate strength for tension) of the impacted wall due to resultant force in the impacted wall occurs at resultant force approximately equal to 10 MN, and the failure of the impacted wall occurs with the probability of 1 at the resultant force approximately equal to 80 MN (Figure 6). The maximum impact force due to the aircraft (the two engine commercial airliner of the most common aircraft used in Lithuania) used in the deterministic transient analysis is 58 MN . Thus, it is seen that the failure (concrete element reaches ultimate strength for tension) of the impacted wall due to resultant force in the impacted wall occurs with the probability approximately equal to 0.9. However, there is estimated that no rebar failure should occur. SUMMARY Calculations show that probability of analyzed aircraft crash occurrence is negligible, however the event occurrence and propagation can be sufficiently uncertain. Estimate of air crash per year probability, with confidence levels that correspond to 5% and 95% quantile, is in the interval between 3.4E-09 and 3.8E-08. This event occurrence is unlikely probable, however similar events are analyzed due to high negative consequence of impact. The results of the probabilistic analyses of failure of the impacted wall show that the aircraft (the two engine commercial airliner of the most common aircraft used in Lithuania) in case of crash will damage the impacted wall of the Ignalina NPP building with the probability of 0.9. However, this damage is related only with wall cracks (concrete element reaches ultimate strength for tension), and no penetration due to the aircraft crash should occur. Even in case when the event analysis show rather limited danger, the sensitivity analysis can identify the highest influence factors. The appropriate control of them is very significant for risk related requirements keeping and risk-based decisions making.
Copyright © 2006 by ASME
ACKNOWLEDGMENTS The research described in this publication was supported in part by the Lithuanian State Science and Studies Foundation. The authors would like to acknowledge the support and access to the newest NEPTUNE code provided by the US DOE and ANL. The authors also want to extend thanks to the administration and technical staff at the Ignalina NPP, for providing information regarding operational procedures and operational data. REFERENCES  Kobayashi, T., 1988, “Probabilistic Analysis of an Aircraft Crash to a Nuclear Power Plant”, Nuclear Engineering and Design, 110, pp. 207-211.
 NUREG/CR-4550, 1986, “Analysis of Core Damage Frequency, Surry Power Station, Unit 1, External Events”, Sandia National Laboratories, SAND86-2084.  Hofer, E., 1999, “Sensitivity analysis in the context of uncertainty analysis for computationally intensive models”, Computer Physics Communications, 117, Elsevier Science, pp. 21-34.  Sachs, L., 1984, Applied Statistics, Second edition, Springer Verlang.  Krzykacz, B., Hofer, E., and Kloos, M., 1994, “A software System for Uncertainty and Sensitivity Analysis of Results from Computer Models”, Proc. Int. Conf. PSAM-II, 2, Session 063, San Diego, CA, pp. 20-25.  Alzbutas R., Augutis J., Maioli A., Finnicum D.J., Carelli M.D., Petrovic B., Kling C.L., Kumagai Y., 2005, “External events analysis and probabilistic risk assessment application for iris plant design”, 13th international conference on nuclear engineering: ICONE 13, Beijing, China, May 16-20, ISBN 7-5022-3400-4, pp. 1-8.  Cesare, M. A., and Sues, R.H., 1999, “PROFES Probabilistic Finite Element System – Bringing Probabilistic Mechanics to the Desktop”, American Institute of Aeronautics and Astronautics, AIAA 99-1607, pp. 111 p.  Kulak, R. F., and Fiala, C., 1988, “NEPTUNE, A System of Finite Element Programs for Three Dimensional Nonlinear Analysis”, Nuclear Engineering and Design, Vol. 106, pp. 47-68.  Kulak, R. F., and Marchertas, P. V., 2003, “Development of a Finite Element Based Probabilistic Tool”, Transactions 17th International Conference on Structural Mechanics in Reactor Technology, Prague, Czech Republic, August 12-17, Paper # B01-2.  Dundulis G., Ušpuras E., Kulak R., 2005, “Structural reliability of an Ignalina NPP building subjected to an airplane crash”, Proceedings of the structural engineering convection SEC 2005, Bangalore, India, December 14 -16, pp. 1-8.
Copyright © 2006 by ASME