Typhoon-Induced Failure Process and Collapse Mechanism of Super-Large Cooling Tower Based on WRF-CFD-LS/DYNA Nesting Technology

: There have been several cases of large cooling towers being damaged by wind in history. A typhoon has the characteristics of a strong wind ﬁeld energy and large shear wind speed. This paper simulates the entire collapse process of large hyperbolic cooling towers by the action of typhoons and reﬁnes the typhoon-induced failure mechanism for cooling towers. Firstly, based on WRF-CFD wind ﬁeld downscaling technology, a ﬁne simulation of the near-ground multiscale wind ﬁeld produced by China’s strongest typhoon “Typhoon Rammasun” is performed to extract effective three-dimensional (3D) typhoon load input parameters. Then, by loading the obtained 3D wind load on the ﬁnite element model, a pseudo-dynamic analysis of the world’s tallest cooling tower “Luan Cooling Tower” is performed based on LS-DYNA explicit dynamic analysis, and the typhoon-induced collapse process is simulated. Finally, the stress distribution and distortions of the tower and the response time history of key units are compared and analyzed to determine the collapse mechanism. The process of collapse begins with large deformation of the windward surface of the tower throat, which shows folds in the range of 62 ◦ on both sides. Eventually, collapse occurs due to uncoordinated deformation. The collapse mechanism can be divided into a bending arch mechanism and a suspension wire mechanism.


Introduction
There are several cases [1][2][3][4] of cooling tower damage caused by wind, suggesting that their structural design needs greater emphasis on the temporal-spatial distribution of wind loads. With the implementation of the "Turn on the big generator set and shut down the small generator set" policy of the Chinese thermal power industry, the trend of new cooling towers becoming bigger and bigger has made them more sensitive to wind. Significantly, the structural characteristics of tall and thin shells increase the wind damage problem [5], especially in southeastern coastal areas of China where typhoons are frequent. China's current wind resistance design codes [6] for cooling towers are based on the normal wind. Typhoons involve strong wind field energy, high-shear wind speeds, and complex near-surface wind profiles [7]. If these are not fully considered, the aerodynamic and structural stability of existing cooling tower systems may be potentially unsafe; therefore, it is necessary to study the mechanisms of super-large cooling tower collapse due to typhoons.
Existing research [8][9][10] on the typhoon resistance of buildings is mostly based on limited measured data. The typhoon engineering models commonly used in structural wind engineering are somewhat simplified and have difficulty in truly reflecting the typhoon wind field characteristics of strong wind field energy, high shear wind speeds, and disordered wind fields. Due to the huge range of typhoons, its minimum simulation The ring plate foundation is a cast-in-situ reinforced concrete structure with a width of 10

Typhoon
In this paper, the Super Typhoon Rammasun is considered, which occurred in the Pacific Ocean in 2014. It was the strongest typhoon (140 Kt) to land in China since the founding of the People's Republic of China. Based on the final analysis data (FNL) of the Global Forecasting System (GFS), it landed in Wenchang, Hainan, China at 15:00 on 18 July 2014. At the time of landing, its maximum wind speed was above 60 m/s (Category 17 typhoon) and recorded a minimum pressure of 899 hPa.

WRF-CFD Wind Field Downscaling Technology
Based on the WRF-CFD technology, a downscaled numerical simulation of the typhoon-induced wind field was performed to extract the three-dimensional distribution of wind loads on the surfaces of the large cooling tower tube. This provided wind load input parameters for the subsequent simulation of the entire collapse process. Figure 2 shows the flowchart of WRF-CFD wind field downscaling technology, where "*" represents the omitted character.

WRF Mode Introduction
According to the meteorological scale, typhoons are mesoscale weather systems and their mesh resolution is usually on the order of kilometers. The WRF model is currently one of the most advanced mesoscale weather research and forecast systems. This paper adopts the WRF model, which is based on the non-static equilibrium Euler equation model, to comprehensively consider the impact of physical schemes such as pavement

Typhoon
In this paper, the Super Typhoon Rammasun is considered, which occurred in the Pacific Ocean in 2014. It was the strongest typhoon (140 Kt) to land in China since the founding of the People's Republic of China. Based on the final analysis data (FNL) of the Global Forecasting System (GFS), it landed in Wenchang, Hainan, China at 15:00 on 18 July 2014. At the time of landing, its maximum wind speed was above 60 m/s (Category 17 typhoon) and recorded a minimum pressure of 899 hPa.

WRF-CFD Wind Field Downscaling Technology
Based on the WRF-CFD technology, a downscaled numerical simulation of the typhooninduced wind field was performed to extract the three-dimensional distribution of wind loads on the surfaces of the large cooling tower tube. This provided wind load input parameters for the subsequent simulation of the entire collapse process. Figure 2 shows the flowchart of WRF-CFD wind field downscaling technology, where "*" represents the omitted character.
The ring plate foundation is a cast-in-situ reinforced concrete structure with a width of 10.5 m and a height of 2.2 m.

Typhoon
In this paper, the Super Typhoon Rammasun is considered, which occurred in the Pacific Ocean in 2014. It was the strongest typhoon (140 Kt) to land in China since the founding of the People's Republic of China. Based on the final analysis data (FNL) of the Global Forecasting System (GFS), it landed in Wenchang, Hainan, China at 15:00 on 18 July 2014. At the time of landing, its maximum wind speed was above 60 m/s (Category 17 typhoon) and recorded a minimum pressure of 899 hPa.

WRF-CFD Wind Field Downscaling Technology
Based on the WRF-CFD technology, a downscaled numerical simulation of the typhoon-induced wind field was performed to extract the three-dimensional distribution of wind loads on the surfaces of the large cooling tower tube. This provided wind load input parameters for the subsequent simulation of the entire collapse process. Figure 2 shows the flowchart of WRF-CFD wind field downscaling technology, where "*" represents the omitted character.

WRF Mode Introduction
According to the meteorological scale, typhoons are mesoscale weather systems and their mesh resolution is usually on the order of kilometers. The WRF model is currently one of the most advanced mesoscale weather research and forecast systems. This paper adopts the WRF model, which is based on the non-static equilibrium Euler equation model, to comprehensively consider the impact of physical schemes such as pavement

WRF Mode Introduction
According to the meteorological scale, typhoons are mesoscale weather systems and their mesh resolution is usually on the order of kilometers. The WRF model is currently one of the most advanced mesoscale weather research and forecast systems. This paper adopts the WRF model, which is based on the non-static equilibrium Euler equation model, to comprehensively consider the impact of physical schemes such as pavement processes, microphysical processes, and wave radiation. By simulating the characteristics of airflow, air pressure, and wind field of mesoscale typhoon weather, the inflow boundary conditions of subsequent CFD models were extracted.

Physical Scheme Selection and Parameter Settings
In the WRF model, the Arakawa C mesh is used in the horizontal direction and the terrain follows the mass coordinate system in the vertical direction. Through the use of triple-nested mesh technology, the WRF simulation accuracy is 1.5 km. The mesh areas of each layer are referred to as D01, D02, and D03, respectively. Area D03 is the center of the simulation area and is set as the Guangdong region (110.3 • N, 20.3 • E). The outermost coarse mesh spacing is 13.5 km and the number of nodes is 211 × 211. The second mesh spacing is 4.5 km with 217 × 217 nodes. The innermost mesh spacing is 1.5 km with 241 × 241 nodes. The map projection uses the Lambert isometric projection scheme. Based on the Hybrid Vertical Coordinate, the computational domain is divided into 30 layers along the height direction.
During the numerical simulation of the typhoon, after several screening tests, the selection was the Mellor-Yamada-Janjic' (MYJ) boundary layer scheme and the Kain-Fritsch cumulus convection parameterization scheme. To analyze the daily changes in the boundary layer, Typhoon Rammasun was simulated for 48 h. For precise numerical simulation, the first 7 h was selected as the mode's spin-up time. The output frequency of the simulation result was once per hour. The detailed parameters of the WRF model are shown in Table 1. For details on each scheme, please refer to the WRF user manual [25].  Figure 3 shows the simulation area of the WRF model and a map of Typhoon Rammasun's path. The blue line is the best path provided by the Tropical Cyclone Data Center of the China Meteorological Administration. The red line is the WRF-simulated path. From the chart, it can be seen that Typhoon Rammasun, as a whole, moved from the southeast coast to the northwest one. There is a certain deviation between the simulated and observed paths. In this paper, the classic turbulent closure scheme is used, which considers local turbulent transport. Although the typhoon wind speed can be simulated well, there is a certain deviation in the wind direction. This study focuses on the extraction of near-surface wind profiles under the action of Typhoon Rammasun. In general, the simulation of Typhoon Rammasun by the WRF model can effectively provide the required wind profile for the CFD inlet boundary.   Figure 4 shows the comparison curve between the maximum stable wind speed near the typhoon center simulated by the WRF model and the data [26] measured by the Japan Meteorological Agency (JMA). The comparison curve shows that the trend of the maximum wind speed near the typhoon center over time is consistent with the corresponding measured data, which both increase first and then decrease. The initial and mid-to-late simulation data are in good agreement with the measured data, while the mid-term simulation data exhibit a slight overestimation with respect to the measured data. The WRF simulation data concern the maximum wind speed in all areas near the typhoon center, while the JMA measured data concern the maximum wind speed at the location of the wind measurement mast.  The effectiveness of the typhoon simulation is verified by analyzing the wind field characteristics of the typhoon. Figure 5 shows the WRF model simulation results for the wind field sea-level pressure, wind speed, equivalent potential temperature, and absolute vorticity when Typhoon Rammasun was logged on 18 July 2014. The analysis shows that the typhoon center has a low-speed and low-pressure distribution. The wind speed and pressure around the center gradually increase and the high-level airflow with higher potential temperature flows into the boundary layer through the coiling effect and then mixes and warms the low-level airflow. The strong wind field energy of the typhoon causes the absolute vorticity of the surrounding wind field to weaken outward from the typhoon center. The wind pressure in the central area of the typhoon simulation is 930-970 hPa and the wind speed is 30-34 m/s, which is consistent with the information recorded by the Central Meteorological Typhoon Network.  Figure 4 shows the comparison curve between the maximum stable wind speed near the typhoon center simulated by the WRF model and the data [26] measured by the Japan Meteorological Agency (JMA). The comparison curve shows that the trend of the maximum wind speed near the typhoon center over time is consistent with the corresponding measured data, which both increase first and then decrease. The initial and mid-to-late simulation data are in good agreement with the measured data, while the mid-term simulation data exhibit a slight overestimation with respect to the measured data. The WRF simulation data concern the maximum wind speed in all areas near the typhoon center, while the JMA measured data concern the maximum wind speed at the location of the wind measurement mast.   Figure 4 shows the comparison curve between the maximum stable wind speed near the typhoon center simulated by the WRF model and the data [26] measured by the Japan Meteorological Agency (JMA). The comparison curve shows that the trend of the maximum wind speed near the typhoon center over time is consistent with the corresponding measured data, which both increase first and then decrease. The initial and mid-to-late simulation data are in good agreement with the measured data, while the mid-term simulation data exhibit a slight overestimation with respect to the measured data. The WRF simulation data concern the maximum wind speed in all areas near the typhoon center, while the JMA measured data concern the maximum wind speed at the location of the wind measurement mast.  The effectiveness of the typhoon simulation is verified by analyzing the wind field characteristics of the typhoon. Figure 5 shows the WRF model simulation results for the wind field sea-level pressure, wind speed, equivalent potential temperature, and absolute vorticity when Typhoon Rammasun was logged on 18 July 2014. The analysis shows that the typhoon center has a low-speed and low-pressure distribution. The wind speed and pressure around the center gradually increase and the high-level airflow with higher potential temperature flows into the boundary layer through the coiling effect and then mixes and warms the low-level airflow. The strong wind field energy of the typhoon causes the absolute vorticity of the surrounding wind field to weaken outward from the typhoon center. The wind pressure in the central area of the typhoon simulation is 930-970 hPa and the wind speed is 30-34 m/s, which is consistent with the information recorded by the Central Meteorological Typhoon Network. The effectiveness of the typhoon simulation is verified by analyzing the wind field characteristics of the typhoon. Figure 5 shows the WRF model simulation results for the wind field sea-level pressure, wind speed, equivalent potential temperature, and absolute vorticity when Typhoon Rammasun was logged on 18 July 2014. The analysis shows that the typhoon center has a low-speed and low-pressure distribution. The wind speed and pressure around the center gradually increase and the high-level airflow with higher potential temperature flows into the boundary layer through the coiling effect and then mixes and warms the low-level airflow. The strong wind field energy of the typhoon causes the absolute vorticity of the surrounding wind field to weaken outward from the typhoon center. The wind pressure in the central area of the typhoon simulation is 930-970 hPa and the wind speed is 30-34 m/s, which is consistent with the information recorded by the Central Meteorological Typhoon Network. Appl. Sci. 2022, 12, x FOR PEER REVIEW 6 of 23  Figure 6 shows some WRF simulation results for cloud maps of the velocity streamline of Typhoon Rammasun during different landing stages. As can be seen from the figure, before the typhoon landed, there was a complete vortex structure in the typhoon eye area, cloud wall area, and spiral rain belt area. After the typhoon passed and advanced to the northwest, the surrounding airflow continued to converge toward the center of the typhoon. The sea and land characteristics of Hainan Island separated the vortex structure and the two regions showed a trend of chaotic wind speeds and irregular weakening. From the simulation results of the wind field, it can be seen that the WRF model can effectively simulate the landing and weakening of Typhoon Rammasun. It also has a certain ability to simulate the movement of weak low-pressure circulation on land after landing.  Figure 6 shows some WRF simulation results for cloud maps of the velocity streamline of Typhoon Rammasun during different landing stages. As can be seen from the figure, before the typhoon landed, there was a complete vortex structure in the typhoon eye area, cloud wall area, and spiral rain belt area. After the typhoon passed and advanced to the northwest, the surrounding airflow continued to converge toward the center of the typhoon. The sea and land characteristics of Hainan Island separated the vortex structure and the two regions showed a trend of chaotic wind speeds and irregular weakening. From the simulation results of the wind field, it can be seen that the WRF model can effectively simulate the landing and weakening of Typhoon Rammasun. It also has a certain ability to simulate the movement of weak low-pressure circulation on land after landing.   Figure 6 shows some WRF simulation results for cloud maps of the velocity streamline of Typhoon Rammasun during different landing stages. As can be seen from the figure, before the typhoon landed, there was a complete vortex structure in the typhoon eye area, cloud wall area, and spiral rain belt area. After the typhoon passed and advanced to the northwest, the surrounding airflow continued to converge toward the center of the typhoon. The sea and land characteristics of Hainan Island separated the vortex structure and the two regions showed a trend of chaotic wind speeds and irregular weakening. From the simulation results of the wind field, it can be seen that the WRF model can effectively simulate the landing and weakening of Typhoon Rammasun. It also has a certain ability to simulate the movement of weak low-pressure circulation on land after landing.   Figure 7a shows the wind speed profiles near the center of the typhoon at different times. It can be seen that when the typhoon landed on 18 July, as a result of the influence of ground roughness, the closer to the ground the greater the change in wind speed. After some hours, the near-surface wind speed decreased by 33%. To compare the difference between the wind field of the normal wind and the typhoon, the speed profile of class A landform (GBT 50102-2014 2014) is studied with the same method, setting the wind speed of the normal wind at 300 m equivalent to that of the typhoon wind field. Based on this, the inflow profile boundary conditions for the typhoon and the normal wind of the subsequent small-scale CFD model were extracted. Figure 6b shows a comparison diagram of the average wind speed profile between the typhoon and the normal wind, which was obtained by the nonlinear Least-Squares Method (LSM). The near-surface wind speed of the typhoon increases with altitude faster than that of normal wind. In the expression of the standard wind speed profile, the ground roughness index of normal wind is 0.12, the fitted value of the typhoon ground roughness index is 0.087, while the basic wind speed at the reference height of the normal wind and the typhoon was 14.6 m/s and 16.59 m/s, respectively.

CFD Model Introduction
At present, high-rise buildings, large-span stadiums, and super-long bridges belong to small-scale structural systems whose overall dimensions are only in the order of 100 m. To accurately predict vortex shedding and aerodynamic loads on a structure's surface, the boundary layer mesh size is strictly limited in the order of centimeters. This paper adopts the CFD model to obtain equivalent static wind loads and wind field pulsation character-  Figure 7a shows the wind speed profiles near the center of the typhoon at different times. It can be seen that when the typhoon landed on 18 July, as a result of the influence of ground roughness, the closer to the ground the greater the change in wind speed. After some hours, the near-surface wind speed decreased by 33%. To compare the difference between the wind field of the normal wind and the typhoon, the speed profile of class A landform (GBT 50102-2014 2014) is studied with the same method, setting the wind speed of the normal wind at 300 m equivalent to that of the typhoon wind field. Based on this, the inflow profile boundary conditions for the typhoon and the normal wind of the subsequent small-scale CFD model were extracted. Figure 6b shows a comparison diagram of the average wind speed profile between the typhoon and the normal wind, which was obtained by the nonlinear Least-Squares Method (LSM). The near-surface wind speed of the typhoon increases with altitude faster than that of normal wind. In the expression of the standard wind speed profile, the ground roughness index of normal wind is 0.12, the fitted value of the typhoon ground roughness index is 0.087, while the basic wind speed at the reference height of the normal wind and the typhoon was 14.6 m/s and 16.59 m/s, respectively.  Figure 7a shows the wind speed profiles near the center of the typhoon at different times. It can be seen that when the typhoon landed on 18 July, as a result of the influence of ground roughness, the closer to the ground the greater the change in wind speed. After some hours, the near-surface wind speed decreased by 33%. To compare the difference between the wind field of the normal wind and the typhoon, the speed profile of class A landform (GBT 50102-2014 2014) is studied with the same method, setting the wind speed of the normal wind at 300 m equivalent to that of the typhoon wind field. Based on this, the inflow profile boundary conditions for the typhoon and the normal wind of the subsequent small-scale CFD model were extracted. Figure 6b shows a comparison diagram of the average wind speed profile between the typhoon and the normal wind, which was obtained by the nonlinear Least-Squares Method (LSM). The near-surface wind speed of the typhoon increases with altitude faster than that of normal wind. In the expression of the standard wind speed profile, the ground roughness index of normal wind is 0.12, the fitted value of the typhoon ground roughness index is 0.087, while the basic wind speed at the reference height of the normal wind and the typhoon was 14.6 m/s and 16.59 m/s, respectively.

CFD Model Introduction
At present, high-rise buildings, large-span stadiums, and super-long bridges belong to small-scale structural systems whose overall dimensions are only in the order of 100 m. To accurately predict vortex shedding and aerodynamic loads on a structure's surface, the boundary layer mesh size is strictly limited in the order of centimeters. This paper adopts the CFD model to obtain equivalent static wind loads and wind field pulsation character-

CFD Model Introduction
At present, high-rise buildings, large-span stadiums, and super-long bridges belong to small-scale structural systems whose overall dimensions are only in the order of 100 m. To accurately predict vortex shedding and aerodynamic loads on a structure's surface, the boundary layer mesh size is strictly limited in the order of centimeters. This paper adopts the CFD model to obtain equivalent static wind loads and wind field pulsation characteristics of the cooling tower. Based on user-defined functions (UDF), the typhoon wind profile obtained by the WRF model is integrated into the Fluent velocity inlet to implement WRF-CFD wind field downscaling. On this basis, fine 3D typhoon load input parameters for subsequent explicit dynamic analysis are extracted.

Mesh Division and Parameter Settings
To ensure that the wake can fully develop, the calculation domain size is set as x = 16 × D, y = 20 × D, and z = 4 × H, where x is the crosswind direction, y is the downwind direction, z is the vertical direction, D is the diameter of the bottom of the tower, and H is the height of the tower. The maximum blockage rate of the tower shape does not exceed 5%, which meets the requirements of the design code. To facilitate calculation efficiency and accuracy while considering the complex shape of the cooling tower and the more complex internal wind field near the structure, the entire calculation domain is divided by a hybrid network. The discrete form of the mesh divides it into external regions and internal encrypted regions. The external regions are divided by high-quality hexahedral structured mesh, while the internal encrypted regions are divided by tetrahedral unstructured mesh.
Before determining the final mesh scheme, mesh independence was verified. Table 2 shows the mesh quality and pressure coefficient of the windward surface of the tower throat area under different mesh schemes. As can be seen from the table, by increasing the total number of cells, the mesh quality, mesh skew, and windward pressure coefficient gradually converge. There is no significant difference in the mesh quality and calculation results when using 12.8 million or 18.8 million cells. Considering the efficiency and accuracy comprehensively, this paper selects the total number of 12.8 million cells as the final mesh scheme. In the internal encryption area, the minimum mesh size of the structural model is 1m and that of other calculation domains is 15m. The wall surface y + value is from 29.8 to 35.2, which can ensure that the underlying mesh is located in the logarithmic law layer. In conclusion, the quality and quantity of the mesh meet the calculation requirements. Due to the full-scale model and standard air materials, the Reynolds number of the numerical simulation is on the order of 10 8 , which is similar to the Reynolds number of the actual cooling tower structure at the design wind speed. The boundary conditions of the calculation domain include velocity inlets, pressure outlets, free-slip walls (sides and tops), and non-slip walls (ground and tower walls). Figure 8 shows the CFD mesh division, measurement point layout, and boundary conditions. The outlet boundary conditions are defined as pressure outlets, where the relative pressure is set to 0. Non-slip walls are used on the ground and the tower surface of the calculation domain and free-slip walls are used on both the sides and top surface of the calculation domain. The wind field velocity inlet was defined by a user-defined function (UDF) of the above-mentioned mesoscale WRF model typhoon field simulation results. The turbulence intensity profile [27] adopts the distribution form of the Chinese design code A-type landform. The above velocity entry is shown as follows: where z is the height above the ground; α is the ground roughness index; U 0 is the basic wind speed; I 10 = 0.12 is the nominal turbulence intensity at a height of Z 0 = 10 m in the area; k is turbulent kinetic energy; ε is the turbulent dissipation rate; C µ = 0.09 is an empirical constant; l is the scale of turbulence.
where z is the height above the ground;  is the ground roughness index; 0 U is the basic wind speed; The numerical calculation uses a 3D single-precision pressure-based solver. The air model uses an incompressible air wind field, the turbulence model uses the RNG k-ε turbulence model, and the pressure-speed coupling equations are solved using the SIMPLEC algorithm. The gradient discretization uses the least-squares cell-based option. The pressure dispersion adopts the standard format, the dynamic dispersion adopts the secondorder upwind style, and the calculation residual of the governing equation is set to 1 × 10 −6 . For the surface roughness characteristic simulation, the surface roughness parament is set by making the simulation results of the average wind pressure distribution on the surface of the cooling tower correspond to the existing measured results (GBT 50102-2014 2014).

Validity and Aerodynamic Characteristics
To verify the validity of CFD, the pulsating wind field of the design code (GBT 50102-2014 2014) type A landform was numerically simulated based on the small-scale CFD model. The wind pressure coefficient is defined as follows: The numerical calculation uses a 3D single-precision pressure-based solver. The air model uses an incompressible air wind field, the turbulence model uses the RNG k-ε turbulence model, and the pressure-speed coupling equations are solved using the SIM-PLEC algorithm. The gradient discretization uses the least-squares cell-based option. The pressure dispersion adopts the standard format, the dynamic dispersion adopts the secondorder upwind style, and the calculation residual of the governing equation is set to 1 × 10 −6 . For the surface roughness characteristic simulation, the surface roughness parament is set by making the simulation results of the average wind pressure distribution on the surface of the cooling tower correspond to the existing measured results (GBT 50102-2014 2014).

Validity and Aerodynamic Characteristics
To verify the validity of CFD, the pulsating wind field of the design code (GBT 50102-2014 2014) type A landform was numerically simulated based on the small-scale CFD model. The wind pressure coefficient is defined as follows: where µ i,θ is the wind pressure coefficient; C Pi,θ is the dimensionless wind pressure; Z i is the height of measuring point i; h is the height of reference point; P i is the pressure of measuring point i; P 0 and P ∞ are the total pressure and static pressure at the reference height, respectively. Figure 9 shows the numerical analysis results of the wind pressure coefficient of the typical weak throat layer (layers 8, 9, and 10) of the type A wind field and compares it with the results of existing related research [6,28,29]. The analysis shows that the numerically predicted wind pressure coefficients in the upwind and leeward areas are in good agreement, while in the crosswind area they are slightly larger than those of the norm; however, they are more consistent with the results of the wind tunnel test.
i Z is the height of measuring point i; h is the height of reference point; i P is the pressure of measuring point i; 0 P and P  are the total pressure and static pressure at the reference height, respectively. Figure 9 shows the numerical analysis results of the wind pressure coefficient of the typical weak throat layer (layers 8, 9, and 10) of the type A wind field and compares it with the results of existing related research [6,28,29]. The analysis shows that the numerically predicted wind pressure coefficients in the upwind and leeward areas are in good agreement, while in the crosswind area they are slightly larger than those of the norm; however, they are more consistent with the results of the wind tunnel test. To compare and analyze the aerodynamic characteristics of the typhoon wind field, Figure 10 shows the average wind pressure coefficient distribution of each layer and Figure 11 shows a schematic comparison of the outer surface pressure coefficient of a cooling tower between typhoon and type A normal wind. Based on this, the accurate 3D typhoon wind load distribution parameters that can be used as input for subsequent analysis of the typhoon-induced collapse of the entire cooling tower were extracted. As can be seen from the figure, the maximum wind pressure coefficient of each layer of the cooling tower appears on the 0° meridian; furthermore, this maximum is greater at the top of the tower (Layer 1) and becomes smaller towards the bottom of the tower (Layer 12). The minimum wind pressure coefficient of each layer of the cooling tower appears at 72° and 288° meridian. Compared with the normal wind, the gradient of the average pressure coefficient of the windward area of the tower is more significant under the action of the typhoon. The average pressure coefficients of the windward and crosswind surfaces in the typhoon case are increased significantly compared with the normal wind case. In particular, the maximum increase of the windward side was 46.8% and the maximum increase of the crosswind surface was 34.6%, which all occurred at the bottom of the tower. Results are as expected due to wind profile change with higher wind load and a lower height. To compare and analyze the aerodynamic characteristics of the typhoon wind field, Figure 10 shows the average wind pressure coefficient distribution of each layer and Figure 11 shows a schematic comparison of the outer surface pressure coefficient of a cooling tower between typhoon and type A normal wind. Based on this, the accurate 3D typhoon wind load distribution parameters that can be used as input for subsequent analysis of the typhoon-induced collapse of the entire cooling tower were extracted. As can be seen from the figure, the maximum wind pressure coefficient of each layer of the cooling tower appears on the 0 • meridian; furthermore, this maximum is greater at the top of the tower (Layer 1) and becomes smaller towards the bottom of the tower (Layer 12). The minimum wind pressure coefficient of each layer of the cooling tower appears at 72 • and 288 • meridian. Compared with the normal wind, the gradient of the average pressure coefficient of the windward area of the tower is more significant under the action of the typhoon. The average pressure coefficients of the windward and crosswind surfaces in the typhoon case are increased significantly compared with the normal wind case. In particular, the maximum increase of the windward side was 46.8% and the maximum increase of the crosswind surface was 34.6%, which all occurred at the bottom of the tower. Results are as expected due to wind profile change with higher wind load and a lower height. Appl. Sci. 2022, 12, x FOR PEER REVIEW 11 of 23

. Finite Element Model
A 3D finite element model of the super-large cooling tower is established based on the macro-model modeling method. The shell of the tower uses Element Shell163, which is divided into 128 layers along the height direction and 256 units along the ring direction. The X-type pillars use Element Beam161, which is divided into two beam elements by cross position. To achieve displacement coordination, the shell element of the upper tower and the beam element of the lower column are coupled by the rigid domain. The ground uses a rigid wall plane. Figure 12 shows a schematic diagram of the finite element model of the cooling tower.

Finite Element Model
A 3D finite element model of the super-large cooling tower is established based on the macro-model modeling method. The shell of the tower uses Element Shell163, which is divided into 128 layers along the height direction and 256 units along the ring direction. The X-type pillars use Element Beam161, which is divided into two beam elements by cross position. To achieve displacement coordination, the shell element of the upper tower and the beam element of the lower column are coupled by the rigid domain. The ground uses a rigid wall plane. Figure 12 shows a schematic diagram of the finite element model of the cooling tower. The wall of the super-large cooling tower is thin and its longitudinal and circumferential reinforcement ratios are large. The strong tensile strength of the rebar makes the tower wall form a well-integrated unit, which has obvious characteristics of plastic followup. Madurapperuma [30] found that the concrete and rebar can be considered comprehensively, and the failure behavior of the component during the collapse process is simulated by the unit failure criterion. The plastic follow-up model (PLAW) with the failure mechanism that is adopted in this paper considers the material characteristics of the tower as a whole. Table 3 summarizes the parameters of the material model. Based on the research results of the ultimate bearing capacity of the cooling tower [24], we conducted the material model debugging and determined the material model parameters for this paper. The block Lanczos method was used to analyze the dynamic characteristics of the 3D finite element model of the super-large cooling tower. Figure 13 shows a curve of the 200order natural frequency of the cooling tower as a function of the mode order. According to this graph, the odd and even order frequencies of the cooling tower are the same and have symmetry. The natural frequency distribution is very dense, the fundamental frequency is only 0.68 Hz, and the first 50 order mode frequencies are mainly concentrated The wall of the super-large cooling tower is thin and its longitudinal and circumferential reinforcement ratios are large. The strong tensile strength of the rebar makes the tower wall form a well-integrated unit, which has obvious characteristics of plastic follow-up. Madurapperuma [30] found that the concrete and rebar can be considered comprehensively, and the failure behavior of the component during the collapse process is simulated by the unit failure criterion. The plastic follow-up model (PLAW) with the failure mechanism that is adopted in this paper considers the material characteristics of the tower as a whole. Table 3 summarizes the parameters of the material model. Based on the research results of the ultimate bearing capacity of the cooling tower [24], we conducted the material model debugging and determined the material model parameters for this paper. The block Lanczos method was used to analyze the dynamic characteristics of the 3D finite element model of the super-large cooling tower. Figure 13 shows a curve of the 200-order natural frequency of the cooling tower as a function of the mode order. According to this graph, the odd and even order frequencies of the cooling tower are the same and have symmetry. The natural frequency distribution is very dense, the fundamental frequency is only 0.68 Hz, and the first 50 order mode frequencies are mainly concentrated between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].
Appl. Sci. 2022, 12, x FOR PEER REVIEW 13 of 23 between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].  Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the low-order mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small. This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32].
Structural failure in large-scale cooling towers under wind loads is similar to that of material failure in buckling mode and the main effect of wind-induced collapse load is the average component [33], so the wind load can be regarded as a pseudo-static load. There is a response spectrum for step force with finite rise time [34]. When the ratio of rising  Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the loworder mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small. between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].  Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the low-order mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small. This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32].
Appl. Sci. 2022, 12, x FOR PEER REVIEW 13 of 23 between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].  Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the low-order mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small.

Contact Model and Loading Method
This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32].
Appl. Sci. 2022, 12, x FOR PEER REVIEW 13 of 23 between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].   Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the low-order mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small.

Contact Model and Loading Method
This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32].
Appl. Sci. 2022, 12, x FOR PEER REVIEW 13 of 23 between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].   Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the low-order mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small.

Contact Model and Loading Method
This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32].
Appl. Sci. 2022, 12, x FOR PEER REVIEW 13 of 23 between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].   Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the low-order mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small.

Contact Model and Loading Method
This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32]. between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].   Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the low-order mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small.

Contact Model and Loading Method
This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32]. between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].   Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the low-order mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small.

Contact Model and Loading Method
This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32]. between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].  Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the low-order mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small.

Contact Model and Loading Method
This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32]. between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].  Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the low-order mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small.

Contact Model and Loading Method
This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32]. between 0.6 and 1.9 Hz. To consider the influence of fluctuating wind pressure, the wind vibration coefficient used in the subsequent explicit dynamic analysis is set to 1.9 [31].  Table 4 shows the typical mode shapes of the cooling tower structure. From the analysis of the graphs, it is found that the vibration modes of the tower are greatly different in the hoop and meridional directions; in particular, with the increase of the order, the number of the hoop and meridional harmonics increases significantly, while in the low-order mode, the meridional direction shows 1 to 3 harmonics and the hoop direction shows 6 to 12 harmonics. Among them, in the first few modes, the tower top vibration is the main type. In the bottom tower vibration mode, there are many ring-shaped harmonics and the tower top amplitude is small.

Contact Model and Loading Method
This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32].

Contact Model and Loading Method
This paper adopted the contact model (*Contact Automatic Single Surface) and the friction coefficient of the material is uniformly set to 0.25. The single-point integration algorithm is used in LS-DYNA to improve calculation efficiency. To weaken the unfavorable zero-energy mode (hourglass mode), the energy consumed by the hourglass mode is controlled by the viscous damping method and the small stiffness method. The specific scheme settings refer to the existing cooling tower collapse simulation results [17,32].
Structural failure in large-scale cooling towers under wind loads is similar to that of material failure in buckling mode and the main effect of wind-induced collapse load is the average component [33], so the wind load can be regarded as a pseudo-static load. There is a response spectrum for step force with finite rise time [34]. When the ratio of rising time to the natural period (tr/Tn) is greater than three, the dynamic coefficient of this step force at a relatively long rise time tends to one, which means that it can be considered a static load.
To reduce the dynamic effect during the loading process, after debugging based on the error rate, the average wind pressure is linearly increased to the specified value within 15 × T1 (T1 being the first-order natural vibration period) to allow for a time history analysis. This is loaded to the corresponding loading point of the tower in the form of point load.
To avoid additional vibration in the tower due to the gravity load, the total load is applied in two stages: First, the gravity load is applied and a high amount of damping is applied to the model at the same time. Iterative calculations are performed until the cooling tower stabilizes under the action of gravity. Then, the "Full Restart" technology is used to redefine the equivalent damping ratio of the cooling tower to 2% [35] and to perform iterative calculations under the effect of wind load.

Validation
To verify the effectiveness of the numerical simulation of the whole process of windinduced collapse of the cooling tower based on CFD and LS-DYNA technology, a wind tunnel test [36] [36].
The buckling wind speed of the cooling tower in the tunnel test was 31 m/s. Based on the flow similarity theory, the actual full-scale collapse wind speed was 88.48 m/s, which has a 9.6% difference from the numerical simulation results in this paper. Additionally, a stability analysis was performed, where the E is the elastic modulus of the cooling tower, q cr is the critical wind pressure, h is the wall thickness of the shell, and r 0 is the throat radius. The numerical simulation results in this article are converted by the same method and compared with the tunnel test results. Figure 14 shows the comparison between the stability analysis results of numerical simulation against those of the tunnel test. The difference between the numerical simulation and the nonlinear stability results is only 0.15%. This is consistent with the basic assumptions of the explicit dynamic analysis (EDA) adopted in this paper, which shows that the numerical simulation method presented in this paper is effective.
time to the natural period (tr/Tn) is greater than three, the dynamic coefficient of this step force at a relatively long rise time tends to one, which means that it can be considered a static load. To reduce the dynamic effect during the loading process, after debugging based on the error rate, the average wind pressure is linearly increased to the specified value within 15 × T1 (T1 being the first-order natural vibration period) to allow for a time history analysis. This is loaded to the corresponding loading point of the tower in the form of point load.
To avoid additional vibration in the tower due to the gravity load, the total load is applied in two stages: First, the gravity load is applied and a high amount of damping is applied to the model at the same time. Iterative calculations are performed until the cooling tower stabilizes under the action of gravity. Then, the "Full Restart" technology is used to redefine the equivalent damping ratio of the cooling tower to 2% [35] and to perform iterative calculations under the effect of wind load. The buckling wind speed of the cooling tower in the tunnel test was 31 m/s. Based on the flow similarity theory, the actual full-scale collapse wind speed was 88.48 m/s, which has a 9.6% difference from the numerical simulation results in this paper. Additionally, a stability analysis was performed, where the E is the elastic modulus of the cooling tower, cr q is the critical wind pressure, h is the wall thickness of the shell, and 0 r is the throat radius. The numerical simulation results in this article are converted by the same method and compared with the tunnel test results. Figure 14 shows the comparison between the stability analysis results of numerical simulation against those of the tunnel test. The difference between the numerical simulation and the nonlinear stability results is only 0.15%. This is consistent with the basic assumptions of the explicit dynamic analysis (EDA) adopted in this paper, which shows that the numerical simulation method presented in this paper is effective.  Figure 14. Comparison of stability analysis results. Figure 14. Comparison of stability analysis results. Figure 15 shows comparative schematic diagrams of the above experiment and the present numerical method for the wind-induced instability collapse of the cooling tower. From the comparison of the figures, it can be seen that the diagrams of both methods have the same collapse trend from the ellipse of the throat to the appearance of folds and finally the tearing and destruction of the rigid ring. The specific distortion shape of the two is only slightly different. Figure 15 shows comparative schematic diagrams of the above experiment and the present numerical method for the wind-induced instability collapse of the cooling tower. From the comparison of the figures, it can be seen that the diagrams of both methods have the same collapse trend from the ellipse of the throat to the appearance of folds and finally the tearing and destruction of the rigid ring. The specific distortion shape of the two is only slightly different.

Weakness Analysis
The basic design wind speed of Shanxi Lu'an Tower is 28.3 m/s. Finite element analysis of the cooling tower under two wind load distribution modes (type A wind field and typhoon wind field) was performed to extract the structural weak points during the collapse of the cooling tower. Figure 16 shows the comparative distribution of the stress ratios (a dimensionless ratio of the first principal stress to the strength of concrete) of the large cooling tower under the action of the two wind fields. It is found from the figure that when the cooling tower is in a stable state, the three-dimensional meridional stress is symmetrical around the meridian axis of 0°. The meridional stress increases from the top of the tower to the throat and the meridional stress from the throat to the bottom of the tower is similar. When the cooling tower is stable under wind load, the hoop stress in the throat is symmetrical around the meridian at 0° and the maximum first principal stress occurs not in the middle but at ±62° in the meridional direction, thereby forming a folding zone around the throat. Compared with the normal wind, the overall stress ratio of the typhoon is similar, but the distribution range is significantly expanded and derived to the leeward area. In summary, in the analysis of the whole process of wind-induced collapse of the cooling tower, attention should be paid to the section of the concrete tower that is located at the intersection of height 40 m-120 m and meridians ±62°. Additionally, the fold area between the highand low-stress ratios, where the concrete tower is prone to shear cracking, must be considered comprehensively.

Weakness Analysis
The basic design wind speed of Shanxi Lu'an Tower is 28.3 m/s. Finite element analysis of the cooling tower under two wind load distribution modes (type A wind field and typhoon wind field) was performed to extract the structural weak points during the collapse of the cooling tower. Figure 16 shows the comparative distribution of the stress ratios (a dimensionless ratio of the first principal stress to the strength of concrete) of the large cooling tower under the action of the two wind fields. It is found from the figure that when the cooling tower is in a stable state, the three-dimensional meridional stress is symmetrical around the meridian axis of 0 • . The meridional stress increases from the top of the tower to the throat and the meridional stress from the throat to the bottom of the tower is similar. When the cooling tower is stable under wind load, the hoop stress in the throat is symmetrical around the meridian at 0 • and the maximum first principal stress occurs not in the middle but at ±62 • in the meridional direction, thereby forming a folding zone around the throat. Compared with the normal wind, the overall stress ratio of the typhoon is similar, but the distribution range is significantly expanded and derived to the leeward area. In summary, in the analysis of the whole process of wind-induced collapse of the cooling tower, attention should be paid to the section of the concrete tower that is located at the intersection of height 40-120 m and meridians ±62 • . Additionally, the fold area between the high-and low-stress ratios, where the concrete tower is prone to shear cracking, must be considered comprehensively. Appl. Sci. 2022, 12 Figure 17 shows the comparative distribution of the horizontal displacement of large cooling towers under the action of the two wind fields. It can be found from the Figure 17 that when the cooling tower is in a stable state, the horizontal displacement is symmetrical around the 0° meridian axis and the maximum horizontal displacement appears 10-20 m below the throat. The horizontal displacement keeps increasing from the top of the tower to the throat. At the meridians of ±62° on both sides of the throat, the tower of the cooling tower protrudes, thereby forming a weak fold and tear zone. The horizontal displacement under the effect of the typhoon is similar to that under the effect of normal wind, but the distribution range is expanded and an overall horizontal displacement occurs. In summary, the position 15 m below the throat is marked as the weak position of the entire collapse process in this paper.  Figure 17 shows the comparative distribution of the horizontal displacement of large cooling towers under the action of the two wind fields. It can be found from the Figure 17 that when the cooling tower is in a stable state, the horizontal displacement is symmetrical around the 0 • meridian axis and the maximum horizontal displacement appears 10-20 m below the throat. The horizontal displacement keeps increasing from the top of the tower to the throat. At the meridians of ±62 • on both sides of the throat, the tower of the cooling tower protrudes, thereby forming a weak fold and tear zone. The horizontal displacement under the effect of the typhoon is similar to that under the effect of normal wind, but the distribution range is expanded and an overall horizontal displacement occurs. In summary, the position 15 m below the throat is marked as the weak position of the entire collapse process in this paper.

Analysis of Collapsed Form
Based on the incremental dynamic analysis (IDA) method, the structure is analyzed non-linearly with the basic wind speed starting at 25 m/s for stepwise loading. The typical weakness (Section 4.2.1) is used as the collapse control member. The loading wind speed step is 5 m/s. When the wind speed increases to the tensile failure of the concrete (tensile strength ≥ 1.71 MPa), the concrete cracks in the local area and the rebar are stretched. As the wind speed increases further, the compression zone of the tower is close to the limit stress state (compressive strength ≥ 19.1 MPa) and the wind-induced response of the cooling tower increases significantly. Therefore, the rapid increase in the maximum displacement of the cooling tower is used as the basis for judging its collapse. Figure 18 shows the maximum displacement of the cooling tower and throat under the typhoon wind field. It can be found from this figure that there is an obvious critical wind speed (80 m/s) in the displacement wind speed curve. Before reaching this critical wind speed, the displacement and wind speed are linear (0-60 m/s) and some nonlinear characteristics begin to appear (60-80 m/s). The cooling tower undergoes a small deformation and is in the approximate elastic stage. After this critical wind speed, the displacement suddenly increases, and the cooling tower undergoes a large deformation and enters the stage of plastic continuous collapse. Research shows that the critical wind pressure of the numerical simulation result is 18.62 kPa. The critical wind pressure of the design code (GBT 50102-2014) when considering the overall stability is 17.67 kPa. The difference between the two is 5.10% and the code is relatively conservative compared with the numerical simulation results. In summary, under the coupling of wind load and self-weight load, the critical wind speed for the collapse of this large cooling tower is 80 m/s.

Analysis of Collapsed Form
Based on the incremental dynamic analysis (IDA) method, the structure is analyzed non-linearly with the basic wind speed starting at 25 m/s for stepwise loading. The typical weakness (Section 4.2.1) is used as the collapse control member. The loading wind speed step is 5 m/s. When the wind speed increases to the tensile failure of the concrete (tensile strength ≥ 1.71 MPa), the concrete cracks in the local area and the rebar are stretched. As the wind speed increases further, the compression zone of the tower is close to the limit stress state (compressive strength ≥ 19.1 MPa) and the wind-induced response of the cooling tower increases significantly. Therefore, the rapid increase in the maximum displacement of the cooling tower is used as the basis for judging its collapse. Figure 18 shows the maximum displacement of the cooling tower and throat under the typhoon wind field. It can be found from this figure that there is an obvious critical wind speed (80 m/s) in the displacement wind speed curve. Before reaching this critical wind speed, the displacement and wind speed are linear (0-60 m/s) and some nonlinear characteristics begin to appear (60-80 m/s). The cooling tower undergoes a small deformation and is in the approximate elastic stage. After this critical wind speed, the displacement suddenly increases, and the cooling tower undergoes a large deformation and enters the stage of plastic continuous collapse. Research shows that the critical wind pressure of the numerical simulation result is 18.62 kPa. The critical wind pressure of the design code (GBT 50102-2014) when considering the overall stability is 17.67 kPa. The difference between the two is 5.10% and the code is relatively conservative compared with the numerical simulation results. In summary, under the coupling of wind load and self-weight load, the critical wind speed for the collapse of this large cooling tower is 80 m/s.

Analysis of Collapsed Form
Based on the incremental dynamic analysis (IDA) method, the structure is analyzed non-linearly with the basic wind speed starting at 25 m/s for stepwise loading. The typical weakness (Section 4.2.1) is used as the collapse control member. The loading wind speed step is 5 m/s. When the wind speed increases to the tensile failure of the concrete (tensile strength ≥ 1.71 MPa), the concrete cracks in the local area and the rebar are stretched. As the wind speed increases further, the compression zone of the tower is close to the limit stress state (compressive strength ≥ 19.1 MPa) and the wind-induced response of the cooling tower increases significantly. Therefore, the rapid increase in the maximum displacement of the cooling tower is used as the basis for judging its collapse. Figure 18 shows the maximum displacement of the cooling tower and throat under the typhoon wind field. It can be found from this figure that there is an obvious critical wind speed (80 m/s) in the displacement wind speed curve. Before reaching this critical wind speed, the displacement and wind speed are linear (0-60 m/s) and some nonlinear characteristics begin to appear (60-80 m/s). The cooling tower undergoes a small deformation and is in the approximate elastic stage. After this critical wind speed, the displacement suddenly increases, and the cooling tower undergoes a large deformation and enters the stage of plastic continuous collapse. Research shows that the critical wind pressure of the numerical simulation result is 18.62 kPa. The critical wind pressure of the design code (GBT 50102-2014) when considering the overall stability is 17.67 kPa. The difference between the two is 5.10% and the code is relatively conservative compared with the numerical simulation results. In summary, under the coupling of wind load and self-weight load, the critical wind speed for the collapse of this large cooling tower is 80 m/s.   Figure 19 shows the displacement schematic diagram and the first principal stress schematic cloud diagram of the numerical simulation of the collapse of the large cooling tower at the critical wind speed of 80 m/s. From the analysis of this figure, it can be seen that the whole process of wind-induced cooling tower collapse can be roughly divided into four stages. First, the throat is depressed and becomes ovoid when viewed from the top. At this time, the first principal stress in the cooling tower extends from the throat to the bottom on the windward side.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 18 of 23 Figure 19 shows the displacement schematic diagram and the first principal stress schematic cloud diagram of the numerical simulation of the collapse of the large cooling tower at the critical wind speed of 80 m/s. From the analysis of this figure, it can be seen that the whole process of wind-induced cooling tower collapse can be roughly divided into four stages. First, the throat is depressed and becomes ovoid when viewed from the top. At this time, the first principal stress in the cooling tower extends from the throat to the bottom on the windward side. Then, the depression continues to develop in the throat and the rigid ring on the top is depressed more slowly. The throat is dragged and two folds are formed obliquely upward from the windward point of the cooling tower throat. At the same time, the first principal stress also develops upward with the pleats transversely. Then, as the rigid ring on the top of the cooling tower continues to sag, the strain reaches the failure point. During the collapse, the rigid ring collapses inward from the windward point and, at the same time, the shear tear at the crease of the throat breaks down. The first main stress develops from the throat to the rigid ring in the meridional direction. Finally, the top rigid ring fragment impacts the back of the cooling tower. Under the combined action of gravity and wind loading, the shell fails and continuous collapse occurs. Compared with the effect of the normal wind field (Wu, 2020), it is found that there are some differences in the weak position and the degree of continuous collapse, but there is a similar failure trend under the effect of typhoon and the conventional wind field. Generally speaking, there is a similar wind-induced failure mechanism. Then, the depression continues to develop in the throat and the rigid ring on the top is depressed more slowly. The throat is dragged and two folds are formed obliquely upward from the windward point of the cooling tower throat. At the same time, the first principal stress also develops upward with the pleats transversely. Then, as the rigid ring on the top of the cooling tower continues to sag, the strain reaches the failure point. During the collapse, the rigid ring collapses inward from the windward point and, at the same time, the shear tear at the crease of the throat breaks down. The first main stress develops from the throat to the rigid ring in the meridional direction. Finally, the top rigid ring fragment impacts the back of the cooling tower. Under the combined action of gravity and wind loading, the shell fails and continuous collapse occurs. Compared with the effect of the normal wind field (Wu, 2020), it is found that there are some differences in the weak position and the degree of continuous collapse, but there is a similar failure trend under the effect of typhoon and the conventional wind field. Generally speaking, there is a similar wind-induced failure mechanism.

Time-History Analysis of Collapse
To analyze the response coupling characteristics of the typical weakness (Section 4.2.1) in the entire collapse process, Figure 20 shows the displacement time history of the typical transient maximum displacement response node (key weak node). Among them, Node16385(A), Node18433(B), Node20481(C), Node22529(D), and Node23809(E) are the maximum displacement response nodes at transient t = 0.2 s, t = 0.4 s, t = 0.4 s, t = 1.4 s, and t = 3.4 s. The comparative study found that the maximum displacement position of the cooling tower first appears at 167.90 m (Node23809(E)). With the progression of collapse, it moves down to 124.48 m (Node16385(A)). Then, the maximum displacement position moves quickly to the tower's top rigid ring area after the rigid ring is broken. In the deformation section of the throat (0 s-2.2 s), the maximum displacement occurs in the throat area. The displacement of each key weak node increases with the degree of collapse. The displacement growth rate of the throat area is higher than that of the rigid ring area. In the transitional period (2.2 s-4.4 s), the displacement growth rate of each key weak node slowed down and the maximum displacement shifts from the throat to the windward point of the rigid ring, while the displacement growth rate of the rigid ring area exceeds the throat area. To analyze the response coupling characteristics of the typical weakness (Section 4.2.1) in the entire collapse process, Figure 20 shows the displacement time history of the typical transient maximum displacement response node (key weak node). Among them, Node16385(A), Node18433(B), Node20481(C), Node22529(D), and Node23809(E) are the maximum displacement response nodes at transient t = 0.2 s, t = 0.4 s, t = 0.4 s, t = 1.4 s, and t = 3.4 s. The comparative study found that the maximum displacement position of the cooling tower first appears at 167.90 m (Node23809(E)). With the progression of collapse, it moves down to 124.48 m (Node16385(A)). Then, the maximum displacement position moves quickly to the tower's top rigid ring area after the rigid ring is broken. In the deformation section of the throat (0 s-2.2 s), the maximum displacement occurs in the throat area. The displacement of each key weak node increases with the degree of collapse. The displacement growth rate of the throat area is higher than that of the rigid ring area. In the transitional period (2.2 s-4.4 s), the displacement growth rate of each key weak node slowed down and the maximum displacement shifts from the throat to the windward point of the rigid ring, while the displacement growth rate of the rigid ring area exceeds the throat area.  Figure 21 shows the time history of the strain energy per unit volume during the entire collapse process. From the perspective of structural energy dissipation, it is found that during the initial small deformation elastic stage (0-2.2 s), the unit volume strain energy is close to 0. There is no way to bear the extreme wind loads by energy consumption. In the plastic development phase of the tower (2.2-8.6 s), the unitary volume of strain energy begins to accumulate at Node A, which has a certain energy consumption effect. After complete collapse (≥8.6 s), each unit bears a certain amount of energy dissipation, but due to the release of gravitational potential energy the ring is completely collapsed. Comprehensive analysis indicates that increasing the duration of the plastic development phase can improve the ability to resist continuous collapse.  Figure 21 shows the time history of the strain energy per unit volume during the entire collapse process. From the perspective of structural energy dissipation, it is found that during the initial small deformation elastic stage (0-2.2 s), the unit volume strain energy is close to 0. There is no way to bear the extreme wind loads by energy consumption. In the plastic development phase of the tower (2.2-8.6 s), the unitary volume of strain energy begins to accumulate at Node A, which has a certain energy consumption effect. After complete collapse (≥8.6 s), each unit bears a certain amount of energy dissipation, but due to the release of gravitational potential energy the ring is completely collapsed. Comprehensive analysis indicates that increasing the duration of the plastic development phase can improve the ability to resist continuous collapse. Appl. Sci. 2022, 12, x FOR PEER REVIEW 20 of 23

Analysis of Force Mechanism
The tower of a large cooling tower is a typical hyperbolic thin-shell system. The collapse process can be summarized as undergoing two stages: a small deformation elastic stage and a large deformation plastic stage. Figure 22 shows the schematic diagram of each stage of the cooling tower collapse and the corresponding mechanical model diagram. The deformation before instability is small and in the elastic stage, and the effects of its geometric nonlinearity and material nonlinearity can be ignored. Under small deformation, large-scale cooling towers mainly resist structural collapse by the bending bearing capacity of their thin shells. The collapse of the tower occurs suddenly. First, local damage occurs near the throat on the windward side, then the fragmentation range continues to expand along the height and around the ring. Under large deformation, a plastic hinge is generated at the windward point of the throat and the tower enters plastic deformation, thus exerting the catenary mechanism.

Analysis of Force Mechanism
The tower of a large cooling tower is a typical hyperbolic thin-shell system. The collapse process can be summarized as undergoing two stages: a small deformation elastic stage and a large deformation plastic stage. Figure 22 shows the schematic diagram of each stage of the cooling tower collapse and the corresponding mechanical model diagram. The deformation before instability is small and in the elastic stage, and the effects of its geometric nonlinearity and material nonlinearity can be ignored. Under small deformation, large-scale cooling towers mainly resist structural collapse by the bending bearing capacity of their thin shells. The collapse of the tower occurs suddenly. First, local damage occurs near the throat on the windward side, then the fragmentation range continues to expand along the height and around the ring. Under large deformation, a plastic hinge is generated at the windward point of the throat and the tower enters plastic deformation, thus exerting the catenary mechanism.

Analysis of Force Mechanism
The tower of a large cooling tower is a typical hyperbolic thin-shell system. The collapse process can be summarized as undergoing two stages: a small deformation elastic stage and a large deformation plastic stage. Figure 22 shows the schematic diagram of each stage of the cooling tower collapse and the corresponding mechanical model diagram. The deformation before instability is small and in the elastic stage, and the effects of its geometric nonlinearity and material nonlinearity can be ignored. Under small deformation, large-scale cooling towers mainly resist structural collapse by the bending bearing capacity of their thin shells. The collapse of the tower occurs suddenly. First, local damage occurs near the throat on the windward side, then the fragmentation range continues to expand along the height and around the ring. Under large deformation, a plastic hinge is generated at the windward point of the throat and the tower enters plastic deformation, thus exerting the catenary mechanism.

Conclusions
This paper proposes a method for numerically simulating the entire process of typhoon-induced collapse in a large-scale cooling tower. It is based on the WRF-CFD-LS/DYNA nesting technology and numerically reproduces the entire process of windinduced collapse of the highest super-large cooling towers in the world (220 m). Based on this, the force characteristics and mechanisms of wind-induced collapse in such super-large cooling towers are studied comparatively, providing a theoretical basis for the windresistant design of large-scale cooling towers based on structural robustness. The research shows that the numerical method described in this paper can effectively simulate the whole process of typhoon-induced collapse in large-scale cooling towers. The main research conclusions are as follows: (1) The Rammasun typhoon wind speed profile obtained based on the WRF model is considered. The fitted value of the typhoon ground roughness index is 0.087, which is 27.5% less than the normal wind speed profile index of Type A, and the basic wind speed of Typhon "Rammasun" is 16.59 m/s, which is increased by 13.6% compared with Type A wind. (2) The pressure at the typical area (15 m below the throat) under typhoon loading is significantly greater than that under the normal wind. The maximum increase in average pressure coefficients at the windward side and side wind surfaces of the cooling tower is 46.8% and 34.6%, respectively. There is a larger amplitude in the 3D aerodynamic distribution of the tower and more severe airflow vortex shedding and turbulence under typhoons. (3) The weak points in the process of wind-induced collapse are located in the high-stress ratio area (40-120 m above the ±62 • meridians of the tower) and the fold area, which is prone to shear failure. The maximum displacement of the cooling tower first appears above the throat. With the progression of collapse, it moves down to the throat and, finally, quickly moves to the tower's top rigid ring area after the rigid ring breaks. The typhoon-induced collapse process begins with a large deformation of the throat, which shows folds within the ±62 • meridians on both sides, eventually leading to collapse due to uncoordinated deformation. (4) The wind-induced collapse mechanism can be divided into two types: an arch mechanism and a suspension line mechanism. In the elastic phase of low deformation, resistance to structural collapse mainly comes from the bending bearing capacity of the arch-shaped thin shell. After entering the plastic phase of large deformation, a plastic hinge is formed at the windward point of the tower top and the structure's continuous collapse is resisted by the axial tension of the tower. Data Availability Statement: Some or all data, models, or codes that support the findings of this study are available from the corresponding author upon reasonable request.