Structural Modeling and Dynamic Analysis of a Nuclear Reactor Building

Abstract

Increasing energy demand urge the developing countries to consider different types of energy sources. Owing the fact that the energy production capacity of renewable energy sources is lower than a nuclear power plant, developed countries like US, France, Japan, Russia and China lead to construct nuclear power plants. These countries compensate 80% of their energy need from nuclear power plants. Further, they periodically conduct tests in order to assess the safety of the existing nuclear power plants by applying impact type loads to the structures. In this study, a sample third-generation nuclear reactor building has been selected to assess its seismic behavior and to observe the crack propagations of the prestressed outer containment. First, a 3D model has been set up using ABAQUS finite element program. Afterwards, modal analysis is conducted to determine the mode shapes. Nonlinear dynamic time history analyses are then followed using an artificial strong ground motion which is compatible with the mean design spectrum of the previously selected ground motions that are scaled to Eurocode 8 Soil type B design spectrum. Results of the conducted nonlinear dynamic analyses are considered in terms of stress distributions and crack propagations.

1. Introduction

Large growth of energy consumption is an essential issue that shall be considered, especially for developing countries. Nuclear power industry seems to take a key role over the continuously increasing energy demand. First known nuclear power station for the generations of electricity was built for commercial use in Russia in 1954, with an output of 5 MW(e). Afterwards, USA constructed a pressurized light water reactor (PWR) with 60 MW(e) which first produced electricity. Many countries over the world followed this first attempt and a rapid growth of nuclear energy production started. Though renewable energy seems to take role in the approaching era, it is estimated that the percentage of the generated electricity from these sources will not be deemed to nuclear power plants. It can be forecasted that up to 2050, nuclear energy and nuclear power plants (NPP) will play a major role in producing electricity worldwide.

Nowadays, supporters of nuclear power admit that this is indeed a technology which is more expensive than its counterparts. Throughout its development, starting from the 1980s, the sector improved the design of reactors. Furthermore, the scientists working in this era aimed to overcome human error or equipment malfunctions. These resulted more robust, fuel efficient and advanced reactor designs where the cost of the construction is also reduced.

Next generation nuclear reactors may be classified into two broad categories: evolutionary and revolutionary [1]. While the Generation III and III+ types belong to the former, Generation IV reactors fell into the latter category. Crimello [2] reported that the design of Generation III type nuclear reactors based on minimizing the risk and thus increase safety. Advanced Boiling and Pressurized Water reactors and Enhanced CANDU 6 are all examples to this type. Passive safety features are incorporated into Generation III type reactors to reach the revolutionary reactors Generation III+. Some examples to these reactors may be: AP1000, VVER 1200 and APR 1400.

Generally, an NPP is composed of five principal buildings: the nuclear island, the annex building, the turbine building, the diesel generator building and radwaste buildings. The nuclear island consists of the containment building, the shield building and the auxiliary building. The containment is one of the most important components of an NPP because it serves as the final barrier under postulated accident conditions.

In the last decade, there is a keen interest on modeling the reactor buildings by using simplified 2D lumped-mass stick model [3, 4, 5, 6, 7]. The results of the analysis may not be considered as accurate as 3D finite element modeling, although the latter one is more time consuming. In the 3D approach, a containment building is generally modeled using either by shell or 3D brick elements [8, 9, 10]. In a recent approach to save time, 3D lumped mass stick models were developed and used in the literature to represent the seismic behavior of containment building [11]. This approach is mostly preferred when the coupling effects between the containment building and auxiliary buildings are considered. Improvements in computational efficiency increased the use of 3D finite element models which are capable of defining high levels of structural detail.

While a possible damage in an infra-structure can be repaired or retrofitted within mean time, a possible damage in an NPP may cause catastrophic damage. The massive damage in an NPP may be observed when an internal accident happens such as loss of coolant accident (LOCA) or when an external event (airplane crash, earthquake, explosions etc.) happens. To overcome these deficiencies the majority of the latest NPP’s are constructed with double containment. Further, in advance, the containment structure constitutes an ultimate barrier against the dissemination of fissile products towards the general public [12].

Kwak and Kim [13] highlighted that most of the recent containments are composed of a dome, a wall and a foundation which are laterally prestressed. The International Federation for Structural Concrete (fib) [14] reported that the response of these shell-type concrete structures due to external events should be experimentally and analytically studied to evaluate their safety. There are numerous experiments in the literature which can be regarded as representative experimental studies. Sandia National Laboratories [15] brought foreword the cost and time-consuming features of the conducted experiments and emphasized that these costly experiments often do not precisely simulate the loading and support conditions of the actual structure.

In this study, seismic behavior of the prestressed outer containment of a third-generation nuclear reactor building is evaluated. Since Turkish Building Earthquake Code (TBEC) [16] does not cover the earthquake resistant design of nuclear structures, international standards such as ASCE 4-16 [17] and ASCE 43-05 [18] have been adopted for the analysis. First, a 3D model of the outer containment vessel has been set up using ABAQUS finite element program. Afterwards, modal analysis is conducted to determine the mode shapes and followed by nonlinear dynamic time history analysis using previously selected ground motions to determine the stress distribution and crack propagation.

2. Details of the 3D analytical model of case-studied nuclear reactor building

Both the Sandia National Laboratories and United States Nuclear Regulatory Commission (USNRC) conducted test programs concerning about the seismic performance of prestressed containments under severe accident conditions. The experiment model was 1:4 prestressed containment and it was conducted in two phases. Determination of the ultimate load capacity and seismic response were carried in the first phase, whereas beyond design basis response was the aim of second phase. More details about this experiment may be found in Hessheimer et al. [19, 20].

A sectional view of the considered test building is illustrated in Figure 1. The structure has an overall height of 16.4 m, the inner radius of the dome is 5.37 m and the cylinder wall has a thickness of 0.325 m. The bottom part of the structure is fixed, gallery and buttresses are incorporated into the system. Prestressing force is applied by using horizontal and vertical tendons which are anchored at the buttresses and gallery, respectively. Around the gallery and buttresses steel rebars are used to enhance the capacity of the double-shell containment building.

Figure 1.

Sectional view of the considered test building.

The three-dimensional model of the test structure has been modeled using shell elements. Smeared steel layers have been used to simulate both the steel reinforcements and prestressing tendons. The created three-dimensional model has been illustrated in Figure 2. To satisfy the experiment conditions the containment is assumed as fix, thus boundary conditions are assigned to the structure at the bottom. In literature two alternatives are valid to define the prestressing force: assigning an initial stress value and altering with the temperature of the structure. Here the former methodology is applied to properly simulate the prestressing force. The abrupt changes around the openings are neglected. Incremental load steps are used to define the internal pressure values.

Figure 2.

Illustration of the three-dimensional model of the reactor building.

Yielding of reinforcing steel and prestressed tendons has been selected as the criteria of failure in any location in the containment. The total mass of the structure has been calculated as 2,956,294 kg and the results of the eigenvalue analysis are summarized in Table 1.

Mode No	Eigenvalue	Frequency
		(rad/time)	(cycles/time)
1	0.36315	0.60262	9.5910E-02
2	0.52506	0.72461	0.11533
3	4.6031	2.1455	0.34146
4	6.6812	2.5848	0.41138
5	12.831	3.5820	0.57009
6	18.840	4.3405	0.69081
7	31.426	5.6059	0.89220
8	40.622	6.3735	1.0144
9	43.410	6.5887	1.0486
10	45.345	6.7339	1.0717

Table 1.

Eigenvalue output.

3. Constitutive material models

3.1 Constitutive model for concrete

The constitutive relation of concrete material is simulated by concrete damaged plasticity model (CDP) of ABAQUS, which is proposed by Lubliner et al. [21] and by Lee and Fenves [22]. This model is preferred by the scientists owing the fact that the model is capable of defining the concrete behavior under cyclic or dynamic loading properly. Further, the material model aims to properly simulate the effects of residual damage. Both the tensile and compressive damages are taken into account in the CDP model. Main contents of the CDP model have been discussed and summarized in the followings.

a. Decomposition of the strain rate.

This is done by assuming an additional strain rate decomposition for the rate independent model which is defined by Eq. (1).

$$\acute{\epsilon }={\acute{\epsilon }}^{\mathit{el}}+{\acute{\epsilon }}^{\mathit{pl}}$$E1

where $\acute{\epsilon }$is the total strain rate, ${\acute{\epsilon }}^{\mathit{el}}$is the elastic part of the strain rate, and ${\acute{\epsilon }}^{\mathit{pl}}$is the plastic part of the strain rate.

b. Stress–strain relations.

This is dominated by scalar damaged elasticity value which can easily be calculated using Eq. (2).

$$\sigma =\left(1-d\right)D_0^{\mathit{el}}:\left(\acute{\epsilon }-{\acute{\epsilon }}^{\mathit{pl}}\right)=D^{\mathit{el}}:\left(\acute{\epsilon }-{\acute{\epsilon }}^{\mathit{pl}}\right)$$E2

where $D_0^{\mathit{el}}$; $D^{\mathit{el}}=\left(1-d\right)D_0^{\mathit{el}}$; and $d$, show the initial stiffness, the degraded elastic stiffness and the scalar stiffness degradation variable, respectively. The scalar stiffness degradation may be associated with the failure mechanism of the concrete which has a range $0\le d\le 1$. Thus, the effective stress is calculated using Eq. (3), and the Cauchy stress by Eq. (4).

$${\overline{\sigma }}^{\mathit{def}}=D_0^{\mathit{el}}:\left(\dot{\epsilon }-{\dot{\epsilon }}^{\mathit{pl}}\right)$$E3

$$\sigma =\left(1-d\right)\acute{\sigma }$$E4

$\left(1-d\right)$in Eq. (4) is the factor that represents the effective load carrying area. When $d=0$, reduction in the elastic stiffness is not expected and it is the condition where $\acute{\sigma }=\sigma$. When $d\ne 0$, damage is possible and use of the effective stress value will give more robust results. Hence, in problems dealing with plasticity $\acute{\sigma }$value is preferred.

c. Hardening variables.

It is also possible to relate $d$value by a set of hardening variables, ${\stackrel{\sim }{\epsilon }}^{\mathit{pl}}$, and $\overline{\sigma }$; such, $d=d\left(\overline{\sigma },{\stackrel{\sim }{\epsilon }}^{\mathit{pl}}\right)$.${\stackrel{\sim }{\epsilon }}_t^{\mathit{pl}}$and ${\stackrel{\sim }{\epsilon }}_c^{\mathit{pl}}$are referred to define damage states in tension and compression, respectively. The governing equation of the hardening variables is given in Eq. (5).

$${\stackrel{\sim }{\epsilon }}^{\mathit{pl}}=\left[\begin{array}{c}{\stackrel{\sim }{\epsilon }}_t^{\mathit{pl}}\\ {\stackrel{\sim }{\epsilon }}_c^{\mathit{pl}}\end{array}\right];{\dot{\stackrel{\sim }{\epsilon }}}^{\mathit{pl}}=h\left(\overline{\sigma },{\stackrel{\sim }{\epsilon }}^{\mathit{pl}}\right).{\dot{\epsilon }}^{\mathit{pl}}$$E5

These variables control the yield surface, the degradation of the elastic stiffness and the dissipated fracture energy.

d. Yield function.

The yield function, $F\left(\overline{\sigma },{\stackrel{\sim }{\epsilon }}^{\mathit{pl}}\right)$, represents a surface in effective stress space, which determines the states of failure or damage. For the inviscid plastic-damage model, $F\left(\overline{\sigma },{\stackrel{\sim }{\epsilon }}^{\mathit{pl}}\right)\le 0$.$F$can also be defined in terms of $\overline{\sigma }$, as in Eq. (6).

$$F\left(\overline{\sigma },{\stackrel{\sim }{\epsilon }}^{\mathit{pl}}\right)=\frac{1}{1-\alpha }\left(\overline{q}-3\alpha \overline{p}+\beta ({\stackrel{\sim }{\epsilon }}^{\mathit{pl}}\right)⟨{\widehat{\overline{\sigma }}}_{\mathit{max}}⟩-\gamma ⟨-{\widehat{\overline{\sigma }}}_{\mathit{max}}⟩)-{\overline{\sigma }}_c\left({\stackrel{\sim }{\epsilon }}_c^{\mathit{pl}}\right)\le 0$$E6

$\alpha$and $\gamma$in Eq. (6) are dimensionless constants. $\overline{p}$and $\overline{q}$are the effective hydrostatic pressure and the Von Mises equivalent effective stress values, respectively. Further these parameters can be calculated by Eqs. (7) and (8). The deviatoric stress may then be expressed using Eq. (9). Here ${\widehat{\overline{\sigma }}}_{\mathit{max}}$refers to the maximum eigenvalue of the effective stress value.

$$\overline{p}=-\frac{1}{3}\overline{\sigma }:I$$E7

$$\overline{q}=\sqrt{\frac{3}{2}\overline{S}:\overline{S}}$$E8

$$\overline{S}=\overline{p}I+\overline{\sigma }$$E9

When ${\widehat{\overline{\sigma }}}_{\mathit{max}}=0$is the case, Eq. (6) reduces to the Drucker-Prager yield statement. Then the $\alpha$coefficient defined in (6) may be determined using the initial equibiaxial (${\sigma }_{b0})$and uniaxial compressive yield stress (${\sigma }_{c0}$) as presented in Eq. (10). Experimental results revealed the fact that $0.08\le \alpha \le 0.12$[21].

$$\alpha =\frac{{\sigma }_{b0}-{\sigma }_{c0}}{2{\sigma }_{b0}-{\sigma }_{c0}}$$E10

The $\beta \left({\stackrel{\sim }{\epsilon }}^{\mathit{pl}}\right)$is then defined using these $\alpha$parameters as given in Eq. (11). Here ${\overline{\sigma }}_t$and ${\overline{\sigma }}_c$are the effective tensile and compressive cohesion stresses, respectively.

$$\beta \left({\stackrel{\sim }{\epsilon }}^{\mathit{pl}}\right)=\frac{{\overline{\sigma }}_c\left({\stackrel{\sim }{\epsilon }}_c^{\mathit{pl}}\right)}{{\overline{\sigma }}_t\left({\stackrel{\sim }{\epsilon }}_t^{\mathit{pl}}\right)}\left(1-\alpha \right)-\left(1+\alpha \right)$$E11

e. Flow rule.

Defining a flow potential value $G$, plastic flow rule may be expressed using Eq. (12). Here $\dot{\lambda }$is a plastic multiplier which always has a positive value. Following the Kuhn-Tucker relations, which are $\dot{\lambda }F=0;\dot{\lambda }\ge 0;F\le 0$, the $G$parameter is defined in the $\overline{\sigma }$space.

$${\dot{\epsilon }}^{\mathit{pl}}=\dot{\lambda }\frac{\partial G\left(\overline{\sigma }\right)}{\partial \overline{\sigma }}$$E12

Using the Drucker-Prager hyperbolic function, the $G$parameter may be defined using Eq. (13). Here $\psi$is the dilation angle at high confining pressure; ${\sigma }_{t0}$is the uniaxial tensile stress at failure; and $ϵ$refers to an eccentricity parameter.

$$G=\sqrt{{\left(ϵ{\sigma }_{t0}\mathit{\text{tanψ}}\right)}^2+{\overline{q}}^2}-\overline{p}\mathit{\text{tanψ}}$$E13

f. Damage and stiffness degradation.

This section will be considered for both uniaxial and cyclic loading protocols.

• Uniaxial loading

The assumption relies on the fact that the uniaxial stress–strain curves can be represented in terms of stress and plastic strain values using Eq. (14).

$${\sigma }_t={\sigma }_t\left({\stackrel{\sim }{\epsilon }}_t^{\mathit{pl}},{\dot{\stackrel{\sim }{\epsilon }}}_t^{\mathit{pl}},\theta ,f_i\right),{\sigma }_c={\sigma }_c\left({\stackrel{\sim }{\epsilon }}_c^{\mathit{pl}},{\dot{\stackrel{\sim }{\epsilon }}}_c^{\mathit{pl}},\theta ,f_i\right)$$E14

In Eq. (14) the subscripts $t$and $c$refer to tension and compression, respectively; ${\dot{\stackrel{\sim }{\epsilon }}}_t^{\mathit{pl}}$and ${\dot{\stackrel{\sim }{\epsilon }}}_c^{\mathit{pl}}$are the equivalent plastic strain rates, ${\stackrel{\sim }{\epsilon }}_t^{\mathit{pl}}={\int }_0^t{\dot{\stackrel{\sim }{\epsilon }}}_t^{\mathit{pl}}\mathit{dt}$and ${\stackrel{\sim }{\epsilon }}_c^{\mathit{pl}}={\int }_0^t{\dot{\stackrel{\sim }{\epsilon }}}_c^{\mathit{pl}}\mathit{dt}$are the equivalent plastic strains, $\theta$is the temperature, and $f_i$, ($i=1,2,..$) are other predefined field variables. Under uniaxial loading conditions the effective plastic strain rates are given as: ${\dot{\stackrel{\sim }{\epsilon }}}_t^{\mathit{pl}}={\dot{\stackrel{\sim }{\epsilon }}}_{11}^{\mathit{pl}}$, in uniaxial tension and ${\dot{\stackrel{\sim }{\epsilon }}}_c^{\mathit{pl}}=-{\dot{\stackrel{\sim }{\epsilon }}}_{11}^{\mathit{pl}}$, in uniaxial compression (ABAQUS 2020). The response of concrete to uniaxial loading is as illustrated in Figure 3. The left illustration is for when the system is under effect of compression only and the right one corresponds to tension case.

Figure 3.

Response of concrete to uniaxial loading.

To properly evaluate the degraded response of concrete two independent damage variables which are assumed as functions of the strain values are introduced: $d_t$and $d_c$. These variables can be calculated following Eqs. (15) and (16). If there is no damage in the material then these variables are assumed to be equal to zero, on the contrary these variables are taken as one when the material is fully damaged.

$$d_t=d_t\left({\stackrel{\sim }{\epsilon }}_t^{\mathit{pl}},\theta ,f_i\right),\left(0\le d_t\le 1\right)$$E15

$$d_c=d_c\left({\stackrel{\sim }{\epsilon }}_c^{\mathit{pl}},\theta ,f_i\right),\left(0\le d_c\le 1\right)$$E16

Defining $E_O$as the initial modulus of the material, the stress values may be determined using Eq. (17).

$${\sigma }_t=\left(1-d\right)E_O\left({\epsilon }_t-{\stackrel{\sim }{\epsilon }}_t^{\mathit{pl}}\right),{\sigma }_c=\left(1-d\right)E_O\left({\epsilon }_c-{\stackrel{\sim }{\epsilon }}_c^{\mathit{pl}}\right)$$E17

When tension loading protocol is the case, cracks propagate in a direction transverse to the stresses. Further, cracks reduce the load bearing capacity of the material and this increases the $\overline{\sigma }$value. On the contrary, when compressive loading is the case, cracks propagate in a direction parallel to the stresses. This minimize the effect of cracks. Naming ${\overline{\sigma }}_t$and ${\overline{\sigma }}_c$as the effective uniaxial cohesion stresses values at tension and compression, respectively, Eqs. (18) and (19) may be applicable to determine the size of the yield surface.

$${\overline{\sigma }}_t=\frac{{\sigma }_t}{\left(1-d_t\right)}=E_0\left({\epsilon }_t-{\stackrel{\sim }{\epsilon }}_t^{\mathit{pl}}\right),$$E18

$${\overline{\sigma }}_c=\frac{{\sigma }_c}{\left(1-d_c\right)}=E_0\left({\epsilon }_c-{\stackrel{\sim }{\epsilon }}_c^{\mathit{pl}}\right)$$E19

• Cyclic loading

The interaction of the cracks makes this type of loading protocol more complex to understand the behavior. Previous experiments revealed the fact that there is indeed some recovery of the stiffness and this is named as “the stiffness recovery effect (SRE)”. This is a significant behavior of concrete under cyclic loading, especially when the load changes from tension to compression. To correlate $d$and $E$values, Eq. (20) has been proposed for the CDP model.

$$E=\left(1-d\right)E_0$$E20

It is previously shown that $d$may be related with $d_t$and $d_c$values. When cyclic loading protocol is applied ABAQUS follows the relation given in Eq. (21). In Eq. (21)$s_t$and $s_c$are stress functions which includes the SRE associated with stress reversals.

$$\left(1-d\right)=\left(1-s_t{d}_c\right)\left(1-s_c{d}_t\right),0\le s_t,s_c\ge 1$$E21

Introducing ${\omega }_t$and ${\omega }_c$material weight factors responsible from the tensile and compressive stiffness recovery, Eqs. (22) and (23) is proposed. The constant $r$in Eqs. (22) and (23) is determined from Eq. (24).

$$s_t=1-{\omega }_{t}r\ast \left({\overline{\sigma }}_{11}\right);0\le {\omega }_t\le 1$$E22

$$s_c=1-{\omega }_c\left(1-r\ast \left({\overline{\sigma }}_{11}\right)\right);0\le {\omega }_c\le 1$$E23

$$r\ast \left({\overline{\sigma }}_{11}\right)=H\left({\overline{\sigma }}_{11}\right)=\left\{\begin{array}{c}1\mathit{\text{if}} {\overline{\sigma }}_{11}>0\\ 0\mathit{\text{if}} {\overline{\sigma }}_{11}<0\end{array}\right.$$E24

This observation in concrete may be accepted as a proof of SRE when cracks close as the load goes from tension to compression. When vice versa is valid, then the tensile stiffness is not recovered. This is ensured in ABAQUS by assuming ${\omega }_t=0$and ${\omega }_c=1$. Cyclic load behavior is illustrated in Figure 4.

Figure 4.

Cyclic loading protocol for${\omega }_t=0$and${\omega }_c=1$.

3.2 Constitutive model for steel

Isotropic hardening model is used to simulate the elastic–plastic steel behavior. As well as the rebars, prestressed tendons are defined using steel material. This indicates that the yield surface changes in all directions. To define the yield function $f\left(\sigma \right)$properly, equivalent stress ${\sigma }^0$is introduced first. This stress value is indeed a function of strain and temperature, as presented in Eq. (25). Here $\theta$is the considered temperature and ${\epsilon }^{\mathit{pl}}$gives the value of the equivalent plastic strain, given in Eq. (26).

$$f\left(\sigma \right)={\sigma }^0\left({\epsilon }^{\mathit{pl}},\theta \right)$$E25

$${\sigma }^0{\acute{\epsilon }}^{\mathit{pl}}=\sigma :{\acute{\epsilon }}^{\mathit{pl}}$$E26

Rice [23] brought forward the use of this model especially when the Bauschinger effect is forceful. The strain rate decomposition may be calculated from Eq. (27). When Eq. (27) is integrated through the contour of the yield surface Eq. (28) is obtained.

$$\mathit{d\epsilon }=d{\epsilon }^{\mathit{el}}+d{\epsilon }^{\mathit{pl}}$$E27

$$\epsilon ={\epsilon }^{\mathit{el}}+{\epsilon }^{\mathit{pl}}$$E28

The elasticity can be then written in terms of two temperature-dependent material parameters. These parameters are generally selected as the bulk modulus, $K$, and the shear modulus, $G$, which are given Eqs. (29) and (30).

$$K=\frac{E}{3\left(1-2\upsilon \right)}$$E29

$$G=\frac{E}{2\left(1+\upsilon \right)}$$E30

The constitutive material parameters used in this study for concrete and steel are as given in Tables 2 and 3, respectively.

Parameters	Values
Elastic modulus (MPa)	3.6 × 10^4
Ultimate tensile strength (MPa)	2.85
Ultimate compressive strength (MPa)	38.5
Poisson ratio	0.2
Density (kg/m^3)	2400
Thermal expansion coefficient	1 × 10^−5
Viscosity parameter	0.005
Dilation angle	36.31
Eccentricity	0.1
${\sigma }_{b0}/{\sigma }_{\mathit{co}}$	1.16
$K_C$	0.667

Table 2.

Parameters of the concrete material.

Parameters	Rebars	Tendons
Elastic modulus (MPa)	2 × 10^5	1.95 × 10^5
Yield strength (MPa)	486.55	1860
Yield strain (m/m)	0.002613	0.008745
Poisson ratio	0.3	0.3
Density (kg/m^3)	7850	7850
Thermal expansion coefficient	1 × 10^−5	1 × 10^−5

Table 3.

Parameters of the steel material.

4. Seismic performance assessment of the selected reactor containment building

Two basic regulations for NPP’s may be named as the 10 CFR 50 and 10 CFR 100 which are published by the Nuclear Regulatory Commission (USNRC) [24, 25]. The former one covers the licensing issues and the latter one explains the steps of evaluating a license. In the latter regulation, two basic earthquakes are defined when performing a seismic hazard assessment: safe shutdown earthquake (SSE) and operating basis earthquake (OBE). While SSE is considered to be the maximum earthquake which could occur at the investigation site, OBE is defined as the earthquake which could be expected to occur at the site during the lifetime of the plant. Further in the regulation, it is stated that the maximum vibratory ground acceleration of the OBE must be at least 33% of the maximum vibratory ground acceleration of the SSE [24, 25].

Nowadays, the methodology for risk analysis involves the use of component fragility curves developed using ground-motion parameters. To obtain the fragility curves the capacity and the demand parameters should be evaluated at first instance. It is a well-known fact that failure of both structural and nonstructural components of an NPP are much more involved with the structural response than the ground parameters.

The FEMA P-58-1 Guidelines [26] provide a basis to improve the risk assessment procedure for NPPs. The guideline develops next-generation tools and new approaches for performance assessment of buildings, with a focus on measuring performance in terms of direct economic loss, casualties and downtime. The FEMA P-58 Guidelines [26] present procedures for performance assessment using a probabilistic framework, which provides a robust methodology to integrate hazard curves, component fragility curves and consequence functions and to capture the dispersions in each of these elements for evaluating the performance of a building. Importantly, the fragility curves used in the analysis are defined in terms of structural response parameters.

ASCE 43-05 [18] provides seismic design criteria for structures, systems, and components (SSC) that are used in nuclear facilities. This is executed by using a graded approach where the design criteria are proportional with the relative importance to safety, magnitude of the seismic event and other factors. To achieve this 20 Seismic Design Bases (SDB) are defined, where for each Seismic Design Category (SDC) probabilistic target performance goals are used. This graded approach is summarized in Table 4. Further in the same guide for each Limit State (LS) an acceptable level of structural response goal is defined. While International Building Code (IBC) [27] may be followed for SDBs defined by SDC1 and 2; for SDBs defined by SDC 3, 4, and 5, ASCE standard should be followed as illustrated in Figure 5. Further in the mentioned code, the Design Basis Earthquake (DBE) is defined by Uniform Hazard Response Spectra (UHRS).

SDC	Limit state
	A; Large permanent distortion (short of collapse)	B; Moderate permanent distortion	C; Limited permanent distortion	D; Essentially elastic
1	SDB-1A; ASCE7	SDB-1B; ASCE7	SDB-1C; ASCE7	SDB-1D; NA
2	SDB-2A; ASCE7	SDB-2B; ASCE7	SDB-2C; NA	SDB-2D; NA
3	SDB-3A; ASCE 43-05	SDB-3B; ASCE 43-05	SDB-3C; ASCE 43-05	SDB-3D; ASCE 43-05
4	SDB-4A; ASCE 43-05	SDB-4B; ASCE 43-05	SDB-4C; ASCE 43-05	SDB-4D; ASCE 43-05
5	SDB-5A; ASCE 43-05	SDB-5B; ASCE 43-05	SDB-5C; ASCE 43-05	SDB-5D; ASCE 43-05

Table 4.

Graded approach defined in ASCE 43-05 [18].

Figure 5.

Seismic design procedure defined in ASCE 43-05 [18] for SDC 3, 4, and 5.

For DBE, Limit State A and D are defined as the intensity of the high and low structural damage, respectively. In other words, Limit State D is where complete elastic behavior is dominant. Damage regions between A and D are called the intermediate levels. The deformation limits associated with each Limit State are described in Table 5.

Limit state	Structural deformation limit
A	Large permanent distortion, short of collapse; Significant damage
B	Moderate permanent distortion; Generally repairable damage
C	Limited permanent distortion; Minimal damage
D	Essentially elastic behavior; No damage

Table 5.

Structural deformation limits for limit state defined in ASCE 43-05 [18].

To properly determine the seismic demand, ASCE 43-05 [18] allows use of linear equivalent static analysis, linear dynamic analysis, complex frequency response methods, or nonlinear analysis.Moreover, when the fundamental mode is dominant then a single step pushover methodology is allowed while conducting nonlinear analysis. These procedures shall be conducted by following the criteria given in FEMA-356 [28] or ATC-40 [29]. This should be clarified that nonlinear analysis is only permitted when beyond design earthquake is considered for low-rise regular NPPs, where higher mode effects are negligible. Otherwise, elastic models are preferred in the seismic design of safety-related structures.

5. Results of the conducted analysis

The 3D finite element model has been subjected to gravity, wind and snow loads in accordance with Eurocode 1 [30], Section 4. First, to determine the fundamental periods and vibration modes of the considered structure, eigenvalue analyses are performed and eigenvectors have been determined. Figure 6 shows the results of the modal analysis considering first six modes.

Figure 6.

Calculated mode shapes and periods.

As it can be inferred from Figure 6, first three modal results reveal the fact that the dome of the reactor building is the most vulnerable section. The displacement value of the dome is calculated approximately 1 mm for the fundamental mode. When considered the higher modes, stress concentrations around operational gaps are observed.

After determining the mode shapes of the outer containment vessel, nonlinear dynamic time history analyses are conducted to assess the seismic behavior of the structure and the stress distribution. For this purpose, PEER database has been used and seven earthquake excitations are first selected in accordance with ATC 58 [31] and scaled to fit Eurocode 8 [32] Soil type B design spectrum. Figure 7 compares Eurocode 8 [32] design spectrum with the spectra of the selected ground motions. The details of the selected strong ground motions are presented in Table 6. Following, an artificial time history record which is compatible with the mean spectrum has been generated and used in the nonlinear dynamic analyses. Figure 8 represents the artificial ground motion.

Figure 7.

Comparison of Eurocode 8 design spectrum with the selected ground motions spectra.

No	Earthquake	Station	Mw	Mechanism	Vs30 (m/s)	Rjb (km)	Rrup (km)
1	Imperial Valley (1979)	Cerro Prieto	6.53	Strike slip	471.53	15.19	—
2	Irpinia Italy (1980)	Rionero In Vulture	6.2	Normal	574.88	22.68	22.69
3	Loma Prieta (1989)	Coyote Lake Dam - Southwest Abutment	6.93	Reverse Oblique	561.43	19.97	20.34
4	Chi-Chi Taiwan (1999)	CHY010	7.63	Reverse Oblique	538.69	19.93	19.96
5	Chi-Chi Taiwan (1999)	TCU075	6.3	Reverse	573.02	24.34	26.31
6	Landers (1992)	North Palm Springs Fire Sta #36	7.28	Strike slip	367.84	26.95	26.95
7	Chuetsu-oki Japan (2007)	Matsushiro Tokamachi	6.8	Reverse	640.14	18.16	25.03

Table 6.

Details of the selected strong ground motions.

Figure 8.

Artificial ground motion.

Modified Newton–Raphson iterative procedure has been followed during the dynamic analysis. ABAQUS offers several convergence norms. Besides stopping the iteration in case of convergence, the iteration process is also stopped if a specified maximum number of iterations has been reached or if the iteration obviously leads to divergence. The detection of divergence is based on the same norms as the detection of convergence. A preferred way to check the convergence is the energy norm.

As the consequence of dynamic analysis total displacement, total deformation and total stress concentration values are illustrated in Figure 9. It should also be highlighted that results of the modal analysis are in good correlation with the dynamic results and the most vulnerable section is assessed as the dome of the structure. In Figure 9(a), total displacement graph has been presented. The maximum displacement value has been calculated as 7.625 mm at the roof. The deformation graph, which is represented in Figure 9(b), is compatible with the stress concentration that is represented in Figure 9(c). Total displacement and stress concentration graphs of the prestressed tendons are given in Figure 10(a) and (b), respectively.

Figure 9.

Nonlinear dynamic time history results of the outer containment: (a) total displacement, (b) total deformation and (c) total stress concentration.

Figure 10.

Nonlinear dynamic time history results of the prestressed tendons: (a) total displacement and (b) total stress concentration.

6. Conclusions

In this study, seismic behavior of the outer containment of a sample reactor building has been analyzed. First, details of the 3D model have been presented, and material models are described in detail. Following that, results of the eigenvalue analysis are discussed. Then, nonlinear time history analyses are conducted using the artificial record which is compatible with the mean design spectrum of the selected ground motions that are previously scaled to Eurocode 8 Soil Type B design spectrum.

The primary outcome of the study is that the most critical section of a reactor building may be assessed as the dome of the structure. More care should be given while designing the dome. In the recent literature, to overcome the deficiencies of prestressed concrete, use of fiber reinforced concrete is proposed due to the fact that these fibers in plain concrete directly effects both the ultimate capacity and post-cracking behavior of a conventional prestressed containment. It may also be interpreted from these ongoing studies that the structural performance of components of an NPP improves when these fibers are used with plain concrete.

It is believed that the study should be developed to consider inelastic behavior of soil and soil-structure interaction effects should be taken into account in the dynamic analysis.

Acknowledgments

This paper is generated from the PhD study of the corresponding author. The author is grateful for the endless support of her supervisor.

Conflict of interest

The authors declare no conflict of interest.