Prediction of ground motion due to the collapse of a large-scale cooling tower under strong earthquakes

Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

Contents lists available at ScienceDirect

Soil Dynamics and Earthquake Engineering journal homepage: www.elsevier.com/locate/soildyn

Prediction of ground motion due to the collapse of a large-scale cooling tower under strong earthquakes Feng Lin a, Hongkui Ji a, Yinan Li b, Zhaoxia Zuo b, Xianglin Gu a,n, Yi Li a a b

Department of Building Engineering, Tongji University, 1239 Siping Road, Shanghai 200092, China State Nuclear Electric Power Planning Design & Research Institute, Beijing 100094, China

art ic l e i nf o

a b s t r a c t

Article history: Received 29 December 2013 Received in revised form 30 May 2014 Accepted 2 June 2014 Available online 24 June 2014

The ground motion owing to the collapse of a large-scale cooling tower under strong earthquakes was appropriately predicted using a comprehensive approach. The predicted results can be used for the safety evaluation of nuclear-related facilities adjacent to the cooling tower as well as in the planning of nuclear power plant construction in China. In this study, a cooling tower–soil model was first developed based on a falling weight–soil model, which the authors verified by falling weight tests. Then the collapse process of a cooling tower was simulated, and the collapse-induced ground vibrations were assessed by using the proposed model. Finally, the ground motion, which was a combination of the earthquake-induced ground motion and the collapse-induced ground vibrations, was estimated based on the superposition principle of waves. It was found that the cooling tower may collapse under strong earthquakes with the peak ground accelerations (PGAs) in the range of 0.35–0.45 g in x (EW) and y (NS) directions, respectively. These PGAs are far beyond the PGA range of major earthquakes in the common seismic design in China. The types of the site geologies of towers can significantly affect the collapse-induced ground vibrations. For a typical hard soil consisting of strongly weathered sandy slate, moderate ground vibrations may occur in the considered region. The collapse-induced PGAs were in the range of 0.017–0.046 g for the observed points at distances of 350 m in radial direction. For a rock-like foundation, the collapse-induced radial PGAs may be as high as 0.08 g at distances of 350 m, indicating that the effect of the collapse-induced ground vibrations on the nuclear-related facilities should be seriously assessed in certain scenarios. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Cooling tower Strong earthquake Peak ground acceleration Collapse Ground vibration

1. Introduction Recently, construction of some nuclear power plants (NPPs) has been planned in inland China. Consequently, larger-scale natural draft wet cooling towers will be constructed as a part of the NPPs due to technological requirements. Most of these cooling towers will be more than 200 m high, representing some of the highest planned cooling towers worldwide [1], with typical spacing of about 300 m from the adjacent nuclear islands due to site limitations. This paper addresses the concern related to accidental events, i.e., strong earthquakes far beyond the design level. In these events the huge towers might collapse and cause ground vibrations, which will eventually detrimentally affect the safe operation of the adjacent nuclear-related facilities. However, this potential risk has not been given serious attention according to the current assessment concepts of NPPs [2–4]. It is generally known that nuclear-related structures in NPPs are designed to have overall better seismic performance than

n

Corresponding author. Tel./fax: +86 021 65982928. E-mail address: [email protected] (X. Gu).

http://dx.doi.org/10.1016/j.soildyn.2014.06.001 0267-7261/& 2014 Elsevier Ltd. All rights reserved.

other structures, e.g., cooling towers [5,6]. As a result, it is possible that the nuclear-related structures will not collapse under a sole earthquake, but cooling towers are not as structurally secure. The collapse-induced ground vibrations by super large cooling towers may enhance the earthquake-induced ground motion, resulting in a combination of ground motions, which can cause unexpected damage to nuclear-related structures and facilities. This scenario is not considered in current design concepts. Therefore, to protect nuclear-related structures from this kind of secondary disaster, collapse-induced ground vibrations should be properly assessed. Such an assessment must at least consider the following two critical issues: ● Intensities of the earthquakes under which cooling towers will collapse, and ● proper prediction of the collapse-induced ground vibrations.

For this purpose, a comprehensive approach is presented for the prediction of the combined ground motion caused by the cooling tower collapse as well as the earthquake itself. The obtained results, in terms of the characteristics of ground motion,

44

F. Lin et al. / Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

can be applied as input data for safety evaluation of nuclearrelated facilities and for the planning of NPPs. However, such applications are beyond the scope of this paper. Efforts have been made to study the collapse of cooling towers numerically. Krätzig and Zhuang [7] simulated the collapse behavior of cooling towers induced by wind, aiming primarily at the evaluation of realistic safety factors for design. A detailed collapse process was not involved. Li et al. [8,9] studied the failure of a largescale cooling tower subjected to accidental loads, i.e., large TNT equivalent mass explosions, aircraft impacts and local failure of pile foundations. By using a three-dimensional finite element method (FEM) based model, various collapse processes and collapse modes were obtained and analyzed. However, the ground vibrations were not involved in their works. On the other hand, studies on impactinduced ground vibration have been continuing for quite a long time. Early studies can be traced back as far as 1904 by Lamb [10]. Subsequently, with the rapid development of computational technology, numerical simulation is being applied. Gu and Lee [11] and Pan and Selby [12] studied the ground response to the dynamic compaction (DC) of dry sand and loose soils, respectively, using FEM based models. As a result, the stress wave attenuation and improvement effects were realistically predicted. Recently, Lin et al. [13] predicted the ground vibrations induced by the collapse of a 235 m high cooling tower in events of a sudden removal of 60 columns and an extremely strong wind, respectively. Both were hypothetical events wherein sudden removal of columns acted transiently and strong wind was simplified as static loads. Therefore these loads significantly differ from earthquakes which are cyclic or hysteric, resulting in fundamental differences in the structural responses. Moreover, on-site tests were also performed on ground vibrations excited by various vibration sources, e.g., the collapse of a chimney with a height of 140 m and a weight of 2600 t [14] as well as a falling weight of 25 t used in DC [15]. It should be noted that we assumed that the cooling towers were built under the conditions of rational sitting, design, construction and maintenance. The “strong earthquakes” considered here were far beyond the design level for cooling towers in NPPs and, therefore, were regarded as accidental loads [6,16]. These accidental loads may occur in small probability cases, e.g., the nuclear disaster at Fukushima, Japan, 2011, which was due to both an earthquake and tsunami. Moreover, different accidental loads can cause different collapse modes of cooling towers, and thus lead to different characteristics of ground vibrations [13]. Therefore, this executed study has both academic and engineering significance.

(4) Developing cooling tower–soil models based on the falling weight–soil model. (5) Predicting the ground motion using the cooling tower–soil models. (6) Proposing suggestions for the design of nuclear-related facilities and for plant planning.

3. Falling weight tests and falling weight–soil model The falling weight tests were conducted at a construction site by means of DC [13] which is a well-known method of soil improvement by treating loose soils on site [11,12]. As illustrated in Fig. 2, a heavy weight was lifted up by a crawler crane to a certain height and then freely dropped, impacting the ground surface to make it compact. Meanwhile, the ground was also excited to vibrate. The considered vibration point was L far from the tamping points. The ground vibrations were recorded using acceleration sensors and a recording, storing and data processing system, which utilized acceleration histories of the considered vibration points at various distances in the vertical and radial directions. These ground vibrations were used to verify a FEM-based falling weight–soil model illustrated in Fig. 3. The model was built using commercial software LS-DYNA [17,18] and several techniques were adopted and carefully checked. First, the commonly used viscous boundaries, which are one type of the non-reflecting boundaries, were set in the undersurface and in the four vertical side surfaces of the soil model. By doing this, the wave energy was absorbed when reaching the viscous boundaries and almost no reflections and refractions occurred, as they actually did in the real infinite soils. The effect of the viscous boundaries was checked using a numerical example. As presented in Fig. 4, a slender and rectangular block was

falling weight

L vibration point and sensors

tamping point

recording, storing and data process -ing system

total station to measure the height of falling weight

Fig. 2. Equipment layout of falling weight tests [13].

2. Research methodology The research methodology adopted is illustrated in Fig. 1 and was achieved through the following steps, in which step 6 (with the box filled in gray) is not presented in this paper. (1) Conducting the falling weight tests. (2) Developing a falling weight–soil model and verifying it with the tests. (3) Choosing the earthquake waves as input excitations.

local enlargement Fig. 3. Falling weight–soil model [13].

Fig. 1. Research methodology adopted in this study.

Displacement (mm)

F. Lin et al. / Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

density: 0.00235 g/mm , dynamic shear modulus: 1810 MPa cohesion: 0.03 MPa Poisson ratio: 0.36, internal friction angle: 25° damping ratio: 0.01 shear wave velocity :849 m/s

Displacement (mm)

historical curve of the inputted wave

45

0.10 0.05

0.00 -0.05 -0.10 0.00

0.25

0.50 T(s)

0.25

0.50 T(s)

0.75

1.00

0.10 0.05 0.00 -0.05 -0.10 0.00

model

0.75

1.00

vibration of point A under (b) compressive and (c) shear wave Fig. 4. Numerical example for checking the effect of the non-reflecting boundaries.

105 N·m)

2.5

Energy (

discretized by solid elements (SOLID164) with its material described by the Drucker–Prager model. The right end of the model was set to the viscous boundary and its left end was excited by a shock compressive and shear wave, respectively. The wave then propagated along the block and might reflect at the right end. The vibration of point A near the left end was recorded in the form of historical displacements. It can be seen that the first peak is clear, which was caused by input wave; while the second peak is almost unrecognizable and is actually only about 5.5% of the first peak, which was caused by reflection wave. These results indicated an appropriate effect of viscous boundary both for the compressive and shear wave. Second, contact and collision actions occurred between the falling weight and soils. To describe these actions, the widely used penalty function method was adopted to implement the contact– collision algorithm. This was implemented by using the keyword “nCONTACT_AUTOMATIC_SINGLE_SURFACE” with its parameters “static coefficient of friction (FD)” and “dynamic coefficient of friction (FS)” both set to be 0.65, with the “exponential decay coefficient (DC)” assumed to be 0 (default). This method was verified by Hou et al. [19]. Third, one center integral point in solid elements was applied, resulting in computational efficiency as well as undesired zero energy modes (hourglass modes). To control hourglass, viscous damping and small elastic stiffness were applied by using the keyword “nHOURGLASS.” The relevant parameters “Hourglass control type (IQH)” and “Bulk viscosity type (IBQ)” were selected as “the Flanagan–Belytschko stiffness form” and “Bulk viscosity type,” respectively. The parameters “Hourglass coefficient (QM),” “Quadratic bulk viscosity coefficient (Q1)” and “Linear bulk viscosity coefficient (Q2)” were set to be 0.10, 1.5 and 0.06 suggested by [18], respectively. Generally, the hourglass is appropriately controlled when the hourglass energy is less than 10% of the total internal energy. In this simulation, the historical development of the hourglass energy and the total internal energy of the soil due to the vibrations are presented in Fig. 5 where it can be seen that the hourglass is well controlled. Fourth, an approximately constant damping (frequency-independent) over a range of frequencies was used for the soil damping. This was implemented by using the keyword “nDAMPING_FREQUENCY_RANGE” with its parameters “damping ratio (DR)” fixed to be 5% based on engineering experience, “lowest frequency (FLOW)” and “highest frequency (FHIGH)” set to be 5 Hz and 15 Hz which were verified by the simulated results, respectively. Finally, satisfactory agreements were achieved both in the acceleration amplitudes and vibration durations [13]. As a result, the falling weight–soil model was well verified.

2.0 1.5

internal energy

1.0

hourglass energy

0.5 0.0 0.0

0.5

1.0 T (s)

1.5

2.0

Fig. 5. Numerical check of the hourglass control for falling weight test.

4. Cooling tower–soil model A FEM-based cooling tower–soil model was developed from the falling weight–soil model for the prediction of ground motions caused by the collapse of a cooling tower under strong earthquakes. The model was built using the commercial finite element program ANSYS/LS-DYNA. 4.1. Modeling of the cooling tower The planned cooling tower was a thin-walled structure constructed using reinforced concrete as shown in Fig. 6. The total height of the cooling tower was 215 m including the height of 46 columns being 19.5 m. The bottom and the top diameter of the shell body were 155 m and 103 m, respectively. The thickness of the shell body varied continuously ranging from 1.8 m (at the bottom) to 0.27 m (at the throat) and further to 0.31 m (at the top). The concrete compressive and tensile strengths were designed to be 29.6 MPa and 2.51 MPa, respectively. The yield strength of the reinforcing steel bars was 400 MPa. A three-dimensional FEM-based model was built as shown in Fig. 7. The shell body was divided into 152 vertical “levels” modeled by four-node shell elements with both bending and membrane capabilities. In this way, the continuously varied shell thickness was properly modeled by changing the thickness of shell elements at each “level.” All the shell elements were divided into 15 “layers” along the direction of thickness. Each layer could be assigned to a specific material, i.e., concrete, meridional steel reinforcements or circumferential steel reinforcements. The fraction of reinforcements in meridional and circumferential directions varied and was specified by design documents. For the modeling of the columns, hexahedral elements (SOLID164) and beam elements (BEAM161) were adopted for the concrete and the

46

F. Lin et al. / Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

103m

throat shell body

215m

ring-beam 155m Fig. 6. Profile of the cooling tower.

circumferential reinforcements

meridional reinforcements

concrete shell element

with [22]. A plastic hardening kinematic model (“nMAT_PLASTIC_KINEMATIC”) was used for the steel bars in columns with consideration of strain rate effect. The parameters “Mass density,” “Young's modulus,” “Poisson's ratio” and “Effective plastic strain” were set to 7800 kg/m3, 2
105 N/mm2, 0.3 and 0.2, respectively. The parameters “Tangent modulus,” “Hardening parameter (BETA)” and “Strain rate parameter for Cowper Symonds strain rate model (SRC)” were empirically adopted as 2000 MPa, 40.4 and 5. For the shell elements, the commonly used material model with the keyword nMAT_CONCRETE_EC2 was applied for the concrete and reinforcing steel bars. Material equations governing the behavior were taken from Eurocode 2 Part 1.2 [23]. Tension softening was governed by the parameter “ECUTEN (strain to fully open a crack)” followed LS-DYNA's suggestion. Strain rate effect was also considered following the suggestions in [22]. Fig. 9 illustrates the stress–strain relationship for concrete and reinforcing steel bars in the shell, respectively, both subject to loading, unloading and reloading in uniaxial compression and tension. Mesh-dependencies were numerically checked, and it was found that the adopted mesh sizes ( E1 m
1 m) were fine enough to describe the dynamic behavior of the cooling tower with consideration of computational efficiency. 4.2. Modeling of the soils

reinforcing steel bars elements

concrete elements

Fig. 7. The FEM-based model of the cooling tower.

solid elements (with embedded beam elements)

shell elements

(D) A

B

nodal fiber C

Fig. 8. Connection of shell elements to solid elements.

reinforcing steel bars, respectively, without consideration of the slip behavior between them. Elements of different types were appropriately connected in the model. Solid and beam elements in the columns were connected at the common nodes. However, the shell elements could not directly connect to the solid and beam elements, because the shell and solid elements had inconsistent nodal degrees of freedom and the moments could not appropriately be transferred. To connect them in a realistic way, we used a specific technique named with the keyword “CONSTRAINED_SHELL_TO_SOLID.” As illustrated in Fig. 8, the nodes A, B and C in a solid element were first constrained as a group, lying along the nodal fiber. Then the node D in the adjacent shell element was tied to the nodes A, B and C. By doing this, the three rotational displacements in the shell–column connections were continuous and compatible, resulting in appropriate transferring of the moments. This function was also used by Hufenbach et al. to solve a similar problem [20]. The Johnson–Holmquist concrete model [21] was adopted for modeling the concrete in the columns. This is a 3D plastic damage model which may describe complex behavior of concrete subjected to large strains, high strain rates and high pressures as in situations related to impacts or explosions. The model parameter “Strain rate coefficient” was adopted as 0.007 in accordance

Strongly weathered sandy slate was adopted in this study to represent the usual site geology of cooling towers in NNPs. The soil properties were excerpted from a practical project and are presented in Table 1. Fig. 10 illustrates the FEM-based model of the soil. The dimension was adopted as 1000 m
1000 m in plane with a thickness of 35 m, which was adequately accurate according to the theory of R-wave propagation and the study results presented by Lou et al. [24]. The viscous boundaries were set in the undersurface and in the four vertical side surfaces of the soil model. Eight-node isoparametric elements (SOLID164) were used to mesh the soil model. It was assumed that the mechanical behaviors of the layered soils were ideal elasto-plastic, which could be described using the Drucker–Prager model [25]. Gu and Lee [11] applied this model and successfully predicted the ground response to the DC of dry sand. An appropriate mesh size of the soil model depended primarily on the accuracy and numerical efficiency. In general, based on the wave propagation theory, the maximum mesh size, le, used in a FEM-based dynamic analysis should fit Eq. (1) [26]:     1 1 1 1 v  λT ¼  ð1Þ le r 12 6 12 6 f T where λT is the wave length corresponding to the dominant wave frequency, fT; v denotes the propagation velocity of the wave under consideration. For the impact-induced ground vibrations, the Rayleigh wave (R-wave) was the predominant component of the surface waves and played a predominant role [27]. Basically, the dominant wave frequency, fT, can be obtained from the Fourier acceleration amplitude spectra of the vibration points presented in Section 5.2 after the numerical computations. The propagation velocity of the R-wave can be derived according to the soil properties in Table 1. Eventually, the mesh size, le, of 5 m
5 m
5 m was adopted and verified by trial computations. Contact and collision actions occurred between tower fragments, as well as between the fragments and the surface soil during the collapse process. To describe these actions, the penalty function method was adopted. In addition, hourglass was well controlled, as the after-mentioned results revealed that the hourglass energy was much less than the total internal energy of the soil due to the vibrations (Fig. 11). Finally, parameter “lowest frequency (FLOW)” and “highest frequency (FHIGH)” were set to

40

600

30

400

Stress (N/mm2)

Stress (N/mm2)

F. Lin et al. / Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

20 10 0 -10 -0.01

0.01 Strain

0.03

47

200 0 -200 -400 -600 -0.004

-0.002

0 0.002 Strain

0.004

Fig. 9. The stress–strain relationship for (a) concrete and (b) reinforcing steel bars in the shell subject to loading, unloading and reloading in uniaxial compression and tension. Table 1 Soil properties for Strongly weathered sandy slate. Type

Density (kg/m3)

Dynamic shear modulus (MPa)

Poisson ratio

Cohesion (MPa)

Internal friction angle (deg)

Damping ratio

Shear wave velocity (m/s)

Strongly weathered sandy slate

2350

1810.0

0.36

0.03

25

0.01

849.0

Fig. 10. The FEM-based model of the soil.

Energy (

108 N·m)

2.5 2.0 1.5

internal energy

1.0

hourglass energy

0.5 0.0 10.0

20.0

30.0

40.0

T (s) Fig. 11. Numerical check of the hourglass control for the cooling tower under earthquake RG1.60 and with strongly weathered sandy slate.

15 Hz and 25 Hz which were also verified by the after-mentioned results, respectively. For modeling the cooling tower and soils, many parameters were used including seven free parameters in the authors' viewpoint. The information for these parameters is presented in Table 2. 4.3. Input of earthquake waves As is well known, appropriate earthquake input is critical for a rational time-history analysis of structures. Generally, the peak ground acceleration (PGA), duration and spectrum characteristics are the essential factors for the description and selection of earthquake waves. In this study six earthquake waves were chosen with their characteristics listed in Table 3. As an example, Fig. 12 presents the acceleration response spectrum (ARS) of the Chi-Chi earthquake (TCU073) and the standard ARS specified in Chinese

code [5]. One can see that the former curve is close to the latter. Similarity in curves was also observed for other earthquake waves. A special strategy was executed in order to save computational time when personal computers (PCs) were used, which was the practice through this study. Obviously it was time consuming when a soil model was activated and incorporated into the computation before the tower fragments impacted the ground. To improve computational efficiency, in our model the columns actually did not connect to the soil surface but were separated from it with a spacing set as 2.0 m. The loads were applied in the form of ground acceleration histories in x (EW), y (NS) and z (vertical) directions, which were defined by the adopted earthquake waves. The loads acted on all the bottom nodes of the elements on column foots, but not the nodes of the soil surface. As a result, it was possible for the separated soil model to be activated and then to take part in the calculation until the tower began to collapse, resulting in a significant reduction of numerical efforts. However, the shortcoming of this strategy was that the interaction between the cooling tower and the soil was ignored. This was because the earlier falling fragments impacted the ground, and then induced ground vibrations which also affected the collapse progress of the residual part of the cooling tower. In spite of this, this effect was believed to be insignificant on the ground motion. As the simulation results in this study indicated, on the one hand, the collapse mode did not change which resulted in similar collapse-induced ground vibrations. On the other hand, the collapse-induced ground vibrations happened always after the occurrence of the PGAs of earthquakes, and its amplitudes were far less than the PGAs of the earthquakes. This meant that the combined ground motions were expected to show no significant differences whether the soil model and cooling tower model were connected or not. An iterative procedure was used to determine the intensities of the earthquakes under which the cooling tower would collapse. The PGAs of the earthquake waves in Table 3 were first normalized, i.e., set to 0.05 g in the two horizontal (x and y) directions and 0.033 g in the vertical (z) direction, respectively, with the ratio 0.05:0.033¼1:0.67 in conformance with Chinese code [28]. Then we performed the iterative computation. If the tower did not collapse, the PGAs of the earthquake waves were added with a fixed step increment of 0.05 g in x and y directions and 0.033 g in z direction until the tower finally collapsed. Note that in such a procedure the ratio of the PGAs in x and y directions actually did

48

F. Lin et al. / Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

Table 2 Information for the parameters in the cooling tower–soil model. Description

Keyword

Concrete in columns

n

Reinforcing steel bars in columns

Parameter

MAT_JOHNSON_HOLMQUIST_CONCRETE Mass density (RO) ρ

n

MAT_PLASTIC_KINEMATIC

Parameter value

2400 kg/m3

Reinforcing steel bars in shell

n

n

MAT_CONCRETE_EC2

MAT_CONCRETE_EC2

Location

No

Technical manual Technical manual Engineering experience Design documents Design documents

Section 4.1

Technical manual Chinese code Chinese code

Section 4.1

Shear modulus (G) G

1.486E þ 04 N/mm

Strain rate coefficient (C)

0.007

Yes

Quasi-static uniaxial compressive strength (FC) Maximum tensile hydrostatic pressure (T)

29.6 MPa

No

2.51 MPa

No

Mass density (RO) ρ

7800 kg/m3

Young' s modulus (E) E 200000 N/mm Poisson' s ratio of reinforcement 0.3 (PR) v Yield stress (SIGY) fy 400 MPa

Concrete in shell

2

Free Source parameter

No

No 2

No No No

Tangent modulus (ETAN) Etan

2000 N/mm2

Yes

Effective plastic strain for eroding elements (FS) Hardening parameter (BETA)

0.2

No

40.4

Yes

Strain rate parameter for Cowper Symonds strain rate model (SRC)

5

Yes

Mass density (RO) ρ

2400 kg/m3

No

Compressive strength of concrete (FC) fc Tensile stress to cause cracking (FT) ft Tension stiffening (ECUTEN)

29.6 MPa

No

2.51 MPa

No

0 (default)

Yes

200,000 N/mm2 0.3

Young' s modulus E Poisson' s ratio of reinforcement (PRREINF) v Ultimate stress of reinforcement (SUREINF) fd Fraction of reinforcement (xaxis) (FRACRX) Fraction of reinforcement (yaxis) (FRACRY)

Design documents Engineering experience Chinese code Engineering experience Engineering experience

Technical manual Design documents Design documents Suggested by LS-DYNA

Section 4.1

No No

Chinese code Chinese code

Section 4.1

400 MPa

No

Specified by design documents and varied along x-axis Specified by design documents and varied along y-axis

No

Design documents Design documents

No

Design documents

Soil

n

MAT_DRUCKER_PRAGER

Mass density (RO) ρ Elastic shear modulus (GMOD) G Poisson ratio (RNU) υ Angle of friction (PHI) Φ Cohesion value (CVAL) c

Table1 of revised version

No

Geological report available

Sections 3 and 4.2

Soil damping

n

DAMPING_FREQUENCY_RANGE

The damping ratio (DR) ξ

1% for stone, 5% for weak soil 15 Hz 25 Hz

Yes

Engineering experience Verified by the numerical results

Sections 3, 4.2 and 6.2.

Model Code 2010 Technical manual Educated guess in [19]

Sections 3 and 4.2

Lowest frequency (FLOW) flow Highest frequency (FHIGH) fhigh

Penalty function method for contact action

n

CONTACT_AUTOMATIC_SINGLE_SURFACE Static coefficient of friction (FD) 0.65 fs Dynamic coefficient of friction 0.65 (FS) fd exponential decay coefficient 0 (DC)

No No

No No Yes

F. Lin et al. / Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

49

Table 2 (continued ) Description

Keyword

Parameter

Parameter value

Free Source parameter

Hourglass control

n

Hourglass coefficient (QM) Quadratic bulk viscosity coefficient (Q1) Linear bulk viscosity coefficient (Q2)

0.1 1.5

No No

0.06

No

HOURGLASS

Location

Suggested by LS- Sections DYNA 3 and 4.2

Note: Free parameter is a variable used in a mathematical model and for scientific modeling. The values of free parameters used in this study were provided by engineering experience or literature.

Table 3 Characteristics of the earthquake waves used in this study. Earthquake name

RG1.60 Chi-Chi Hector Mine Landers Tabas Chi-Chi

Station name

Peak ground accelerations (g)

Artificial wave TCU073 Twentynine Palms Silent Valley Poppet Flat Tabas TCU039

Duration (s)

x direction

y direction

z direction

0.15000 0.03362 0.06581 0.04988 0.83580 0.05478

0.15000 0.02055 0.06706 0.04105 0.85170 0.07093

0.10000 0.01679 0.04307 0.03804 0.68850 0.02405

25.00 64.00 59.96 54.96 32.84 64.00

Note: (1) g denotes the acceleration of gravity. (2) x, y and z direction was set to be the east-west (EW), north-south (NS) and vertical (V) direction, respectively, as shown in Fig. 10.

2.00

2.00 z-direction

1.50

y-direction

1.00

standard (horizontal)

1.50

Acc (m/s2)

Acc (m/s2)

x-direction

0.50

standard (vertical)

1.00 0.50 0.00

0.00 0

1

2

3

0

4

1

2

T (s)

3

4

T (s)

Fig. 12. The acceleration response spectrum (ARS) of the Chi-Chi earthquake (TCU073) and the standard ARS specified in Chinese code [5].

not affect the model behavior. This was because the different amplitude ratios only affected the direction of the resultant wave in this iterative procedure. Fortunately, both the cooling tower and the supported soil are central symmetric, which means that the dynamic behavior of the cooling tower is independent of the direction of the input earthquake wave. On the other hand, spectra characteristics of an earthquake wave in the x and y directions are generally similar (Fig. 12).

5. Ground motion 5.1. Collapse process of the cooling tower Table 4 lists the PGAs of the six selected earthquakes under which the cooling tower collapsed. The collapse modes of the tower under all of the earthquakes were similar. As an example, Fig. 13 shows the collapse process of the tower under the Chi-Chi earthquake (TCU073). For this earthquake, relatively high PGAs began to occur after 10 s and continued on for another 10 s. Consequently, local failure of the cooling tower first appeared on the column tops after approximately 14 s. The failure happened in this region due to concrete crush. This was because the stresses concentrated on the column tops where the dimension was small compared to the ring beam and the column feet. Then failure

Table 4 Peak ground accelerations (PGAs) of the ground motion. Name

PGAs of earthquakes under which the cooling tower collapsed (g)

PGAs of the combined ground motion after superposition (g)

x and y direction z direction x and y direction z direction RG1.60 Chi-Chi (TCU073) Hector Mine Landers Tabas Chi-Chi (TCU039)

0.35 0.35 0.35 0.45 0.35 0.45

0.233 0.233 0.233 0.300 0.233 0.300

0.35 0.35 0.35 0.45 0.35 0.45

0.233 0.233 0.233 0.300 0.233 0.300

extended quickly to the adjacent column tops. The cooling tower began to incline at approximately 21 after about 23 s due to loss of balance of the columns. After that, the huge structure retained such an inclined posture and began to collapse, but did not roll over partially because of its small height–width ratio. The collapse duration lasted for about eight seconds, which was similar to the results of other earthquake cases. However, the moments when the cooling tower began to collapse were different for the six earthquake cases owing to different characteristics of individual earthquake waves. The collapse-induced ground vibrations could

50

F. Lin et al. / Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

weak parts

22nd s

24th s

26th s

27th s

weak part (21.2nd s)

0.00 -0.60 -1.20 20.00

25.00

30.00

1.20 0.60 0.00 -0.60 -1.20 20.00

35.00

L=250m, t-direction

25.00

L=350m, r-direction

0.00 -0.60 -1.20 20.00

25.00

30.00

Acc (m/s2)

Acc (m/s2)

0.60

1.20 0.60

L=450m, r-direction

-0.60 25.00

0.00 -0.60

25.00

0.00

25.00

30.00

1.20 0.60

30.00

35.00

1.20 0.60

0.00

25.00

0.00

T (s)

30.00

T (s)

30.00

35.00

T (s)

L=450m, t-direction

25.00

35.00

-0.60 -1.20 20.00

35.00

-0.60 -1.20 20.00

30.00

L=350m, z-direction

T (s)

0.00

-1.20 20.00

L=250m, z-direction

T (s)

-0.60 -1.20 20.00

35.00

Acc (m/s2)

Acc (m/s2)

0.60

0.60

-1.20 20.00

35.00

L=350m, t-direction

T (s) 1.20

1.20

T (s)

T (s) 1.20

30.00

Acc (m/s2)

L=250m, r-direction

Acc (m/s2)

0.60

35.00

Acc (m/s2)

1.20

Acc (m/s2)

Acc (m/s2)

Fig. 13. Collapse process of the cooling tower under the Chi-Chi earthquake (TCU073).

1.20 0.60

L=450m, z-direction

0.00 -0.60 -1.20 20.00

25.00

30.00

35.00

T (s)

Fig. 14. Collapse-induced acceleration histories of vibration points at various distances under the Chi-Chi earthquake (TCU073).

be detected after about 23 s when the shell body impacted the ground in this example. The horizontal PGAs in Table 4 lay in the range of 0.35–045 g and are far beyond the major earthquakes adopted in the common seismic design specified in China [28]. For example, the usual horizontal PGAs of the major earthquakes in x and y directions are 0.10 g for cooling towers. These results of PGAs can be regarded as reasonable in comparison with the results for the reinforced concrete cooling tower investigated in [29]. Moreover, the collapse process also revealed that the weak parts (shown in Fig. 13) of this cooling tower under earthquakes were the column tops connected to the shell body. The weak parts initially failed and further induced the collapse of the tower.

5.2. Collapse-induced ground vibrations The collapse-induced ground vibrations were obtained in the form of acceleration histories of vibration points at various distances. As an example, the results of the Chi-Chi earthquake (TCU073) are illustrated in Fig. 14 in the r (radial), t (tangent) and z directions for the points at distances (L) of 250 m, 350 m and 450 m. These points were particularly chosen because they vibrated with the maximum PGAs at the distances, as signed in Fig.10. This region with distances 250–450 m was of interest because the nuclear island would most likely be located here. For clarity, the collapse-induced PGAs and vibration durations for the considered earthquakes are listed in Table 5. The derived

Fourier acceleration amplitude spectra of the vibration points at various distances are presented in Fig. 15. It was found that: (1) The observed ground region may experience moderate vibrations in comparison with the earthquakes under which the cooling tower just collapsed. For example, the collapseinduced PGAs of the six earthquakes in the r direction were in the range of 0.017 g to 0.046 g at the distance (L) of 350 m. The collapse-induced PGAs in r, t and z directions for different earthquakes did not show significant variation because the collapse mode of each cooling tower was almost identical. (2) The vibrations attenuated in directions with increase in distance. The attenuation rate showed no statistical rules in directions. Wave propagations are direction oriented and exhibit different attenuation rates depending on the frequency of the vibrations, the type of soils and the radiation damping, etc. [30]. Based on this knowledge, the authors believe that the observed characteristics of the vibration attenuation are casedependent and cannot be considered as general conclusions. (3) The vibration durations of the six earthquake cases were similar and lasted for about 7.1–8.0 s. The dominant frequency band of the vibration points at various distances lay mainly in the range of 15–25 Hz. These two observations mainly contributed to the similar collapse mode of the cooling tower under different earthquakes. To clearly describe the results, a contour map of the collapseinduced PGAs distribution under the Chi-Chi earthquake (TCU073)

F. Lin et al. / Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

51

Table 5 Collapse-induced PGAs and vibration durations. PGAs (g)

Duration (s)

L ¼250 m t

z direction

r

t

z direction

r

t

z direction

0.108 0.054 0.048 0.099 0.057 0.022

0.050 0.012 0.028 0.022 0.025 0.018

0.149 0.113 0.078 0.148 0.053 0.037

0.037 0.046 0.030 0.046 0.039 0.017

0.015 0.011 0.021 0.010 0.012 0.012

0.037 0.061 0.047 0.049 0.031 0.035

0.013 0.014 0.026 0.016 0.011 0.018

0.013 0.012 0.017 0.008 0.011 0.011

0.015 0.033 0.028 0.023 0.012 0.019

L=250m, r-direction

0.01

0.02

0.00

L=250m, t-direction

0.01

0.00 5

10

15

20

25

5

f (Hz)

L=350m, r-direction

0.01

0.02

0.00

10

15

20

5

10

15

0

20

0.01

25

0.02

0.00

5

10

15

20

0.02

20

25

L=350m, z-direction

0

5

10

15

20

25

f (Hz)

L=450m, t-direction

0.01

25

20

0.01

25

0.00

10 15 f (Hz)

15

0.00 0

Fourier Amplitude(m/s)

0.01

5

10

f (Hz)

L=450m, r-direction

0

5

f (Hz)

L=350m, t-direction

f (Hz) 0.02

0.01

25

0.00

0

L=250m, z-direction

f (Hz) Fourier Amplitude(m/s)

Fourier Amplitude(m/s)

0.02

8.0 7.8 7.7 7.3 7.8 7.1

0.00 0

Fourier Amplitude(m/s)

0

Fourier Amplitude(m/s)

0.02

Fourier Amplitude(m/s)

0.02

Fourier Amplitude(m/s)

L ¼450 m

r

Fourier Amplitude(m/s)

RG1.60 Chi-Chi (TCU073) Hector Mine Landers Tabas Chi-Chi (TCU039)

L ¼350 m

0.02 Fourier Amplitude(m/s)

Earthquake

L=450m, z-direction

0.01

0.00 0

5

10 15 f (Hz)

20

25

0

5

10 15 f (Hz)

20

25

Fig. 15. Collapse-induced Fourier acceleration amplitude spectra of vibration points at various distances under the Chi-Chi earthquake(TCU073).

The underground vibrations were also investigated for the safety evaluation of the facilities embedded in the subsoil. Table 6 presents the collapse-induced acceleration amplitudes of the underground vibrations, with respect to depth, H, in the case of the Chi-Chi earthquake (TCU073). One can see that the vibrations attenuated with the increase in depth, indicating that increasing the burial depth contributes to the significant reduction of vibrations.

500 0.05g 0.10g

Distance (m)

250

0.20g

0

5.3. Superposition of ground motions

-250

>0.2g 0.1g 0.2g 0.05g 0.1g <0.05g

-500 -500

-250

0

250

500

Distance (m) Fig. 16. The contour map of collapse-induced PGAs distribution under the Chi-Chi earthquake (TCU073).

is shown in Fig. 16 as an example. This map is convenient for the safety evaluation of nuclear-related facilities and the plant planning of NPPs.

Actually, the ground motion was the combination of the collapse-induced ground vibrations and the earthquake-induced ground motion. To predict the combined ground motion, the widely adopted superposition principle in earthquake engineering was used here. Fig. 17 shows an example of the superposition process of ground vibrations/motion in time order at a distance of 350 m in the z direction, as in the case of the Chi-Chi earthquake (TCU073). One can see that the curve shape shows almost no difference and the PGAs are identical before and after the superposition. Similar results were also observed in other directions and for other earthquake cases as listed in Table 4. This was because, for the collapse-induced vibrations and the earthquake-induced

52

F. Lin et al. / Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

Table 6 Collapse-induced maximum acceleration amplitudes of the underground vibrations for the case of the Chi-Chi earthquake (TCU073). Depth (m)

Maximum acceleration amplitudes (g) L ¼ 250 m

0 10 20

L ¼ 350 m

r

t

z direction

r

t

z direction

r

t

z direction

0.054 0.025 0.019

0.012 0.013 0.011

0.113 0.083 0.053

0.046 0.026 0.017

0.011 0.012 0.009

0.061 0.036 0.014

0.014 0.011 0.004

0.012 0.007 0.007

0.033 0.015 0.011

4.00 z-Acc (m/s2)

L ¼ 450 m

earthquake - induced ground motion

2.00 0.00

-2.00 -4.00 0

10

20

z-Acc (m/s2)

2.00

30 40 T (s)

50

60

collapse-induced ground vibration

1.00 0.00

-1.00 -2.00 0

20

superposition region

4.00 z-Acc (m/s2)

10

2.00

30 40 T (s)

50

60

combined ground motion

0.00

Fig. 18. Collapse mode with “collapse in fragments” for another cooling tower under the Tabas earthquake [31].

collapse modes and the site geology on the collapse-induced vibrations.

-2.00 -4.00 0

10

20

30 40 T (s)

50

60

Fig. 17. The superposition process of ground vibration/motion in time order at distance of 350 m in the z direction under the Chi-Chi earthquake (TCU073).

ground motion, the PGAs of the former were relatively small, i.e., only 17% of those of the latter. Moreover, their corresponding PGAs did not appear simultaneously, i.e., the PGAs of collapse-induced vibrations appeared always after the PGAs of earthquake-induced ground motion in these study cases. It was meaningful to estimate the upper limit value of the combined PGAs, when the actual randomness of earthquakes was considered. It is assumed that: (1) The linear superposition principle is valid for particles vibrating in soils. (2) The PGAs of the collapse-induced vibrations and the earthquakeinduced ground motion appear simultaneously. Based on the obtained results, for the shear wave velocity of the soil being about 849 m/s and the observed region with a distance of 350 m from the cooling tower, the upper limit value of the combined PGAs could be proposed as 1.17 times that of the earthquake-induced PGAs in the most unfavorable conditions.

6.1. Influence of the collapse modes The observed collapse mode in Section 5.1 belonged to “collapse in integrity” [13] or “pancake-type collapse” [31]. As described above, the local failure originally appeared on the column tops and extended to the adjacent ones. Then the shell body began to incline and fall down primarily under dead weight after loss of balance. During this process, the structure impacted the ground and finally disintegrated into fragments. As demonstrated in [13], two collapse modes were observed in different load cases, i.e. “collapse in integrity” in the event of sudden removal of 60 columns, or “collapse in fragments” in the load case of an extremely strong wind in which only part of the shell body disintegrated. It was noted that different collapse modes could lead to significantly different ground vibrations. Even for the same load case of earthquakes, unlike in this study, “collapse in fragments” was also observed for another cooling tower reported in [32] and is presented in Fig. 18, most likely owing to a different design strategy for that cooling tower in that study case. Conceptually, in that study, the ground vibrations were believed to be different from the results in this study, although a direct comparison between the two results could not be performed due to totally different engineering factors, e.g., tower weight and site geology. 6.2. Influence of site geology

6. Discussion As is generally known, the ground vibrations depend on the characteristics of the vibration source and the propagation medium [30], i.e., the site geology. The characteristics of the vibration source were affected critically by the collapse mode of the cooling tower. The following sections will cover the influences of the

To demonstrate the influence of site geology on the collapseinduced vibrations, two additional and common types of site geology for cooling towers were assumed, i.e., homogeneous pebbly clay and quartz sandstone, which each filled the whole soil model, respectively. Their soil properties are presented in Table 7. A comparison on the collapse-induced PGAs for different site geologies at various distances in the r, t and z directions is

F. Lin et al. / Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

53

Table 7 Soil properties for pebbly clay and quartz sandstone. Type

Density (kg/m3)

Dynamic shear modulus (MPa)

Poisson ratio

Cohesion (MPa)

Internal friction angle (deg)

Damping ratio

Shear wave velocity (m/s)

Pebbly clay Quartz sandstone

1900 2580

45.0 9500.0

0.40 0.30

0.04 9.98

12.00 42.04

0.05 0.01

140.0 1929.0

Table 8 Collapse-induced PGAs for different site geologies at various distances. Site geology

PGAs (g)

Duration (s)

L ¼ 250 m

Pebbly clay Strongly weathered sandy slate Quartz sandstone

L ¼ 350 m

L ¼450 m

r-

t-

z-direction

r-

t-

z-direction

r-

t-

z-direction

0.017 0.054 0.172

0.019 0.012 0.054

0.018 0.113 0.421

0.011 0.046 0.080

0.014 0.011 0.033

0.015 0.061 0.161

0.008 0.014 0.031

0.009 0.012 0.017

0.011 0.033 0.036

presented in Table 8. It was found that the ground vibrations at the quartz sandstone site became dramatically more intensive (e.g., about 2.6–11 times for L ¼350 m) and the vibration duration was slightly shorter compared to that of the pebbly clay site. This finding was different from the common knowledge in earthquake engineering where rock-like foundations help to reduce the ground motion. The essential reason is that the vibration source and the concerned wave propagation in a collapse scenario are on the ground surface, not in the underground where earthquake happens. Besides, the proposed upper limit value of the combined PGAs in section 5.3 should also be reevaluated when the soil is different.

7. Conclusions This study presents a comprehensive approach for the prediction of ground motion due to the collapse of a large-scale cooling tower under strong earthquakes. The ground motion was a combination of the collapse-induced ground vibrations and the earthquake-induced ground motion. The developed cooling tower–soil model was based on the falling weight–soil model verified by the falling weight tests. The collapse-induced ground vibrations were predicted in the form of acceleration histories of points at various distances. The following conclusions can be drawn from the study: (1) The cooling tower may collapse under strong earthquakes far beyond the seismic design level, with horizontal PGAs ranging from 0.35 g to 0.45 g. in the x (EW) and y (NS) directions, respectively. The weak parts of the cooling towers under earthquakes were the column tops in this study. (2) Moderate ground vibrations caused by the collapse of the tower may occur in the considered region. The collapseinduced radial PGAs of the selected six earthquakes were in the range of 0.011–0.080 g at a distance of 350 m for shear wave velocities of the soils in the range of 140–1929 m/s. The collapse-induced vibrations attenuated with increase in distance. For the combined ground motion, the estimated upper limit value of PGAs can be proposed as 1.17 times of the earthquake-induced PGAs in this study case. (3) The collapse modes and site geologies can significantly affect the collapse-induced ground vibrations.

8.7 7.8 6.8

In the current design concepts of the NPPs, the considered accidental loads do not include the collapse-induced ground vibrations due to earthquakes. However, this study demonstrates that in certain scenarios, especially those concerning large-scale cooling tower and rock-like foundations, the effect of collapseinduced ground vibrations on nuclear-related facilities should be seriously assessed.

Acknowledgments This research was sponsored by the National High-Tech Development Plan (863 Program) under Grant no. 2012AA050903. The authors would like to extend their sincere gratitude to the Ministry of Science and Technology, China for their supports. References [1] Busch D, Harte R, Krätzig WB, Montag U. New natural draft cooling tower of 200 m of height. Eng Struct 2002;24(12):1509–21. [2] International Atomic Energy Agency (IAEA), Seismic design and qualification for nuclear power plants: safety guide, Safety Standards Series No. NS-G-1.6. Vienna; 2003. [3] National Nuclear Safety Administration, General safety principles for design of nuclear power plants, HAD102/01. Beijing: China Legal Publishing House; 1989 (in Chinese). [4] American Society of Civil Engineers, Seismic analysis of safety-related nuclear structures and commentary, ASCE 4-98. 2000. [5] Standardization Administration of the People’s Republic of China, Code for seismic design of nuclear power plants, GB 50267-97. Beijing: China Planning Press; 1998 (in Chinese). [6] VGB-Technical Committee, Structural design of cooling towers. VGB PowerTech e.V., Essen, Germany; 2005. [7] Krätzig WB, Zhuang Y. Collapse simulation of reinforced concrete natural draught cooling tower. Eng Struct 1992;14(5):291–9. [8] Li Y, Lin F, Gu XL, Lu XQ. Numerical research of a super-large cooling tower subjected to accidental loads. Nucl Eng Des 2014;269:184–91. [9] Li Y, Huang SK, He JT, Gu XL, Shi FM, Lin F Numerical simulation on collapse of a super-large cooling tower caused by local failure of pile foundation. In: Proceedings of 22nd International Conference on Structural Mechanics in Reactor Technology (SMiRT22). San Francisco, California, USA, 2013, Division III, p. 962–71. [10] Lamb H. On the propagation of tremors over the surface of an elastic solid. Philos Trans R Soc Lond 1904;203:359–71. [11] Gu Q, Lee FH. Ground response to dynamic compaction of dry sand. Geotechnique 2002;52(7):481–93. [12] Pan JL, Selby AR. Simulation of dynamic compaction of loose granular soils. Adv Eng. Softw. 2002;33(7–10):631–40. [13] Lin F, Li Y, Gu XL. Prediction of ground vibration due to the collapse of a 235 m high cooling tower under accidental loads. Nucl Eng Des 2013;258:89–101.

54

F. Lin et al. / Soil Dynamics and Earthquake Engineering 65 (2014) 43–54

[14] Stangenberg F, Zinn R Bodenerschütterungen beim sprengtechnischen Abbruch von 140 m hohen Stahlbetonschornsteinen. In: Proceedings of the finite elemente in der baupraxis. Modellierung, Berechnung und Konstruktion. Beiträge zur Tagung FEM’98, Darmstadt, Germany, 1998, p. 471–80 (in German). [15] Kwang JH, Tu TY. Ground vibration due to dynamic compaction. Soil Dyn Earthq Eng 2006;26:337–46. [16] Ministry of Housing and Urban-Rural Construction of the People’s Republic of China. Specification for design of reinforced concrete shell structures, JGJ222012. Beijing: China Architecture & Building Press; 2012 (in Chinese). [17] Hallquist JO. LS-DYNA theory manual. Livermore, California, USA: Livermore Software Technology Corporation; 2006. [18] Hallquist JO. LS-DYNA keyword user’s manual. Livermore, California, USA: Livermore Software Technology Corporation; 2012. [19] HOU J, LIN F, GU XL. Impulse model of impact between concrete blocks. J Vib Shock 2007;26(10):1–6 (in Chinese). [20] Hufenbach W, Gude M, Ebert C, Zscheyge M, Hornig A. Strain rate dependent low velocity impact response of layerwise 3D-reinforced composite structures. Int J Impact Eng 2011;38(5):358–68. [21] Holmquist TJ, Johnson GR, Cook WH A computational constitutive model for concrete subjected to large strains, high strain rates, and high pressures. In: Proceedings of the 14th International Symposium on Ballistics; 1993, p. 591–600. [22] Commitee Euro-International du Beton, CEB-FIP Model Code 1990. Thomas Telford; 1993.

[23] European Committee for Standardization, Design of concrete structures-part 1: general rules and rules for buildings, EN 1992-1-1. Brussels; 2004. [24] Lou ML, Pan DG, Fan LC. Effect of vertical artificial boundary on seismic response of soil layer. J Tongji Univ (Nat Sci) 2003;31(7):757–61. [25] Drucker DC, Prager W. Soil mechanics and plastic analysis for limit design. Q Appl Math 1952;10(2):157–65. [26] Kausel E, Manolis G. Wave motion in earthquake engineering. Boston: WIT Press; 2000. [27] Bolt BA. Earthquake. New York: Freeman; 1988. [28] Standardization Administration of the People’s Republic of China, Code for seismic design of buildings, GB 50011-2010. Beijing: China Architecture & Building Press; 2010 (in Chinese). [29] Sabouri-Ghomi S, Nik FA, Roufegarinejad A, Bradford MA. Numerical study of the nonlinear dynamic behaviour of reinforced concrete cooling towers under earthquake excitation. Adv Struct Eng 2006;9(3):433–42. [30] Athanasopoulos GA, Pelekisa PC, Anagnostopoulos GA. Effect of soil stiffness in the attenuation of Rayleigh-wave motions from field measurements. Soil Dyn Earthq Eng 2000;19:277–88. [31] Starossek U. Progressive collapse of structures. London: Thomas Telford Limited; 2009. [32] Gu XL, Li Y, Ji HK. Study on the appropriate distance of common countercurrent cooling towers with the area of water drenching of 18000 m2 and 24000 m2 to nuclear island. Shanghai: Tongji University; 2013 (in Chinese).