To better evaluate the trend of the horizontal acceleration amplitude, the acceleration amplification factor (AAF) is obtained by dividing the peak horizontal acceleration at each measuring point by the peak horizontal acceleration of input motion. Therefore, the AAF at the base of the model (at soil depth of 1.6 m in Figure 3a) is equal to unity (1). This definition also applies to the following parts. Figure 8 shows the AAF of the ground along the depth under different RVH. The horizontal acceleration amplitude of the soil tends to increase from the bottom to the ground surface. Under the four conditions the AAFs on the ground surface are about 2.3–3.1 and the soil presents significant amplification. This is reasonable since the amplitude of all the input motions are relatively low.

**Figure 8.** Acceleration amplification factor of ground under different RVH.

In terms of the acceleration levels in this study, with increasing RVH, soil horizontal acceleration increases and the increase in the amplitude grows significantly with increasing RVH on the whole. It demonstrates that soil horizontal acceleration is more sensitive to higher RVH or higher vertical acceleration.

#### *3.2. Acceleration Di*ff*erence between the Side Wall and Adjacent Soil*

To investigate the influence of the RVH on the acceleration response difference between the side wall and adjacent soil, Figure 9 illustrates the comparison of the acceleration time histories between side wall and adjacent soil at three different depths, which are ceiling slab (Figure 9a), middle slab (Figure 9b), and base slab (Figure 9c), respectively.

**Figure 9.** *Cont.*

**Figure 9.** *Cont.*

**Figure 9.** Comparison of the acceleration time histories between side wall and adjacent soil. at three different depths: (**a**) comparison at depth of ceiling slab; (**b**) comparison at depth of middle slab; (**c**) comparison at depth of base slab.

As can be seen from Figure 9, the phases of the side wall and adjacent soil are always consistent at three depths. The acceleration difference between the side wall and the adjacent soil mainly lies on the amplitude. For the side wall and the adjacent soil, the amplitude of them is almost the same at the depth of the middle slab, where the amplitude differences are about 0–0.025 g at the peak-strain moment. The amplitude of them has minor differences at depth of base slab, where the amplitude differences are about 0.023–0.035 g at the peak-strain moment. It is worth mentioning that some larger differences between the side wall and adjacent soil can be observed at the depth of the ceiling slab, where the amplitude differences reach about 0.045–0.084 g at the peak-strain moment. Under RVH = 1.5 and horizontal input acceleration of 0.1 g, the aforementioned amplitude difference even reaches 0.084 g. It reveals larger differences in the horizontal acceleration between the side wall and the adjacent soil at the depth of the ceiling slab. The differences can be attributed to that the ceiling slab lies too close to the ground surface and there is a significant amplification of the acceleration responses since fewer soil constraints exist there. It reveals that for the underground structures with zero or near-zero buried depth, such as atrium-style metro stations, special attention should be paid to the acceleration response difference between the structure and adjacent soil.

In order to investigate the differences of the horizontal acceleration amplitude between the side wall and the adjacent soil with increasing RVH, the results at three different depths are depicted in Figure 10. It is found that at any depth, the acceleration differences between side wall and adjacent soil appear a linear increase in the amplitude as RVH increases from 0 to 1.5. With regard to the amplitude, the acceleration differences at the depth of the ceiling slab are larger than those at the depth of the base slab. The acceleration differences at the depth of the middle slab are the smallest ones among the three depths. With regard to the increasing rates, the increasing rate at the depth of the ceiling slab is about 2.56%, which is larger than that of 1.65% at the depth of the middle slab. The smallest one is 0.82% at the depth of the base slab.

In conclusion, for an atrium-style metro station, the differences in the horizontal acceleration amplitude between the structure and the adjacent soil rise with increasing RVH, which are different at different depths. The most significant differences occur at the depth of the ceiling slab. The RVH has a significant influence on the dynamic soil-structure interaction, especially for higher RVH.

**Figure 10.** The difference of horizontal acceleration amplitude with increasing RVH between side wall and adjacent soil at three different depths.

#### *3.3. Influence of RVH on Dynamic Soil Normal Stress*

Dynamic soil normal stress (DSNS) is defined as the difference between the total stress and initial static stress along the side wall of the model station. As an example, Figure 11 compares the two time histories of the DSNS between RVH = 0 and RVH = 1. It is found that the time histories of the DSNS resemble that of the input motion. The essence of this phenomenon is also revealed by Kramer, who points out that the input seismic excitation in the form of an acceleration can be converted into a stress wave [35]. With different RVH their phases are basically consistent at all measuring points, though the amplitudes of the DSNS behave differently at different depths of the side wall. For example, the amplitudes are almost the same near the bottom of the side wall (e.g., PR1 and PR2) while the amplitudes under RVH = 1 are always larger than those under RVH = 0 at other depths of the side wall (e.g., PR3~PR11).

**Figure 11.** *Cont.*

**Figure 11.** *Cont.*

**Figure 11.** Comparison of time histories of dynamic soil normal stress (DSNS) along right side wall between RVH = 0 and RVH = 1.

To figure out the influence of the RVH on the DSNS along the side wall, Figure 12 displays the peak DSNS along the left and right side walls under different RVH.

**Figure 12.** Distribution of peak DSNS under different RVH along (**a**) left side wall, and (**b**) right side wall.

As can be seen from Figure 12a, under the input motion of Loma Prieta, peak DSNS along the left side wall follows an approximate L–shaped distribution. With increasing RVH, on the whole, there is no significant change in the amplitude of the peak DSNS. Compared to other RVH, a relatively large increase at the depth of the platform floor can be observed when RVH equals to 1.5.

With respect to the peak DSNS along the right wall, as shown in Figure 12b, at the depth between the top of the right side wall and the middle of the station hall floor, they all follow an approximately linear distribution. At the depth between the middle of the station hall floor and the bottom of the right side wall, they all follow an approximately –shaped distribution. It is worth mentioning that the maximum stress occurs at the bottom of the side wall basically. With increasing RVH, peak DSNS along the right wall increase on the whole, except for the bottom of the side wall (base slab level). Compared to the results under only horizontal input motion (RVH = 0), the maximum increase in peak DSNS under RVH = 0.5, 1, and 1.5 reach about 148%, 159%, and 199%, respectively. For the bottom of the side wall (base slab level), with increasing RVH, the peak DSNS has no significant change, which is different

from the situation at other position of the side wall and has also been concluded above from Figure 11. It might be attributed to the sharp transition in the corner and a possible stress concentration there.

When comparing the distribution of the peak DSNS between the left and right side walls, some obvious differences between them can be found. A similar finding was also drawn from a numerical study of a shallow buried rectangular underground structure under earthquake loading with both the horizontal and vertical components [36]. The differences in stress distribution of the peak DSNS between the left and right side walls can be explained by the shadow effect, which is caused by the entrapment of the waves between the ground and the underground structure [37]. With increasing RVH, the above-mentioned differences between the left and right side walls also increase. From the perspective of dynamic soil normal stress, the experimental results have proven that RVH has significant influence on the dynamic soil-structure interaction.

#### *3.4. Influence of RVH on Structural Dynamic Strain*

The dynamic strain is defined as the difference between the total strain during shaking and initial static strain before shaking. As an example, Figure 13 compares the time histories of the dynamic strain between RVH = 0 and RVH = 1 from one side of every cross section (Figure 3d). Just as the characteristics of the DSNS mentioned earlier, the time histories of the dynamic strain also resemble that of the input motion. In terms of the two different RVH, their phases are extremely consistent for all the measuring points and their amplitudes have minor differences for most of the measuring points.

Figure 14 shows the variation of the peak dynamic tensile strain (DTS) with different RVH for different measuring points of the model station. To easily observe the rate of change in the dynamic strain over different RVH, the value of the y–coordinate in Figure 14 represents the ratio of the peak DTS under RVH = i (i = 0, 0.5, 1, and 1.5) to that under RVH = 0. Since two strain gauges are arranged for every cross section (Figure 3d), the strain results on one side are displayed in Figure 14a and the opposite ones are shown in Figure 14b, where the layouts of corresponding strain gauges are plotted again. As can be seen from Figure 14a, with increasing RVH, the DTS will change in an undulating fashion (decrease, increase, and decrease again), which is different from the response characteristics of the soil, in that the soil acceleration at any depth increases rapidly with increasing RVH. It means the influence of the RVH on the structural DTS is rather complex.

**Figure 13.** *Cont.*

**Figure 13.** Comparison of time histories of dynamic strain on one side of structural components between RVH = 0 and RVH = 1.

**Figure 14.** Peak dynamic tensile strain (DTS) of model station under different RVH: (**a**) strains on one side; and (**b**) strains on another side.

The soil and structure behave differently in the seismic responses with increasing RVH. It can also be seen from Figure 14a that the DTS with bidirectional input motion may be less than those with only horizontal input motion (the value of the y–coordinate is less than 1). This conclusion is verified by the numerical simulation conducted on a prototype atrium-style metro station [31], and can also be found in [38]. This phenomenon can be explained as follows: the presence of a cavity can cause, under certain conditions, intense and selective de-amplification of the free-filed motion, which is referred to as "shadow zone" [37]. When comparing the results of the DTS in Figure 14a,b, for every cross section, the DTS on one side is not equal to the opposite one. It reveals that the stresses in the cross section of every structural element are uneven under bidirectional input motion.

#### **4. Conclusions**

An atrium-style metro station is a favorite choice owing to its excellence in providing a larger and clearer space for the passengers and commercial clients at the station hall floor. However, concern for the seismic safety of such a station is growing since its hall slabs are replaced with the flat-beams. Tests are needed for the evaluation of its seismic safety. A series of 1 g shaking table tests were conducted to investigate the seismic responses of an atrium-style metro station.

In this study, extensive experimental results were presented to investigate the influence of the RVH (ratio of vertical to horizontal peak ground acceleration) on the seismic responses of both the soil and the underground metro station. Results show that the acceleration amplification factors (AAFs) on the ground surface are about 2.3–3.1, and that the soil presents significant amplification under several considered conditions. Soil horizontal acceleration is more sensitive to the higher RVH or higher vertical acceleration. Significant differences in horizontal acceleration amplitude between the side wall and the adjacent soil at the depth of the ceiling slab are found, and attention should be paid to these kind of differences for the underground structures with zero or near-zero buried depth, such as atrium-style metro stations. With increasing RVH, the differences in the horizontal acceleration amplitude between the side wall and the adjacent soil also rise. The RVH has a significant influence on the dynamic interaction, especially for higher RVH. Under the input motion of Loma Prieta, the distribution of the peak DSNS (dynamic soil normal stress) along the left and right side walls represents significant difference. The peak DSNS along the left side wall follows an approximate L-shaped distribution, while it is different for the right side wall. With increasing RVH, there is no significant change in the amplitude of the peak DSNS along the left side wall on the whole, while it increases a lot (the maximum increase reach about 148%, 159%, and 199% under RVH = 0.5, 1, and 1.5 when comparing with that under RVH = 0) for the right side wall. The soil and structure behave differently in the seismic responses with increasing RVH. The stresses in the cross-section of every structural element are uneven under horizontal–vertical earthquake excitation.

In conclusion, from the perspective of the soil acceleration, structural acceleration, and dynamic soil normal stresses along the side wall of station, experimental results have proven that ratio of vertical to horizontal peak ground acceleration (RVH) has significant influence on the dynamic soil-structure interaction. It is believed that under extreme earthquake loading, such as near fault zones, RVH is a parameter of paramount importance and should be accounted for in the seismic analyses and seismic performance assessments of the underground structures, especially for those with zero or near-zero buried depth, such as atrium-style metro stations.

It is worth mentioning that further investigation of the effect of the vertical earthquake motion is essential through a series of numerical analyses using a minimum number of the earthquake records specified by a certain code. In this way, a cross-check of the findings from the multiple records could be expected to accomplished and this is a work in progress.

**Author Contributions:** Test scheme, Z.Z., Y.Y., H.Z., and H.Y.; experimental works, Z.Z.; discussion and analysis, E.B. and Y.Y.; writing—original draft preparation, Z.Z.; writing—review and editing, E.B. and Y.Y.; funding acquisition, Y.Y. and H.Z.

**Funding:** This research was funded by the National Key Research and Development Plan of China (Grants Nos. 2018YFC1504305, 2018YFC0809600, 2018YFC0809602, and 2017YFC1500703), the State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University (Grant No. SLDRCE15-02), and the China Scholarship Council (Grant No.<PHONE_NUMBER>43).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

## **An Analytical Framework for the Investigation of Tropical Cyclone Wind Characteristics over Di**ff**erent Measurement Conditions**

### **Lixiao Li 1,2, Yizhuo Zhou 1, Haifeng Wang 3, Haijun Zhou 1,2, Xuhui He 4 and Teng Wu 3,\***

Received: 1 November 2019; Accepted: 6 December 2019; Published: 9 December 2019

**Abstract:** Wind characteristics (e.g., mean wind speed, gust factor, turbulence intensity and integral scale, etc.) are quite scattered in different measurement conditions, especially during typhoon and/or hurricane processes, which results in the structural engineer ambiguously determining the wind parameters in wind-resistant design of buildings and structures in cyclone-prone regions. In tropical cyclones (including typhoons and hurricanes), the inconsistent wind characteristics may be in part ascribed to the complex flow structure with the coexistence of both mechanical and convective turbulence in the boundary layer of tropical cyclones. Another significant contribution to the scattered wind characteristics is due to various measurement conditions (e.g., terrain exposure and height) and data processing schemes (e.g., averaging time). The removal of the inconsistency in the field-measurement system may offer a more rational comparison of measured wind data from various observation platforms, and hence facilitates a better identification scheme of the wind characteristics to guide the urban planning design and wind-resistant design of buildings and structures. In this study, an analytical framework was firstly proposed to eliminate the potential observation-related effects in wind characteristics and then the wind characteristics of seven field measured tropical cyclones (four typhoons and three hurricanes) were comparatively investigated. Specifically, field measurements of wind characteristics were converted to a standard reference station with a roughness length of 0.03 m, observation duration of 10 min for mean wind and averaging time of 3 s for gusty wind at a 10 m height. The differences of the measured wind characteristics between the typhoons and hurricanes were highlighted. The standardized turbulent wind characteristics under the analytical framework for typhoons and hurricanes were compared with the corresponding recommendations in standard of American Society of Civil Engineers (ASCE 7-10) and Architectural Institute of Japan Recommendations for Loads on Buildings (AIJ-RLB-2004).

**Keywords:** wind characteristics; boundary layer; typhoon; hurricane; field measurement

#### **1. Introduction**

Wind characteristics (e.g., mean wind speed, gust factor, turbulence intensity and integral scale, etc.) are the critical factors for wind-resistant design of the wind-sensitive infrastructures and urban planning. Resisting wind effects and reducing wind-induced damage in tropical cyclones is the challenge for the wind sensitive buildings and structures in cyclone-prone regions, as these regions are normally in the economically developed areas with crowded populations and large-scale landmark buildings and structures. Therefore, a rational analytical framework for investigating wind characteristics in tropical cyclones is essential to understand the nature of winds, calibrate codes of practice for wind-resistant design of the large-scale structures and enhance wind tunnel simulations and numerical modeling [1]. Oncoming winds of buildings and bridges are usually simplified as the steady flow part featured by mean wind speed and corresponding vertical profile, and the fluctuating flow part characterized by turbulence intensity, integral scale, gust factor, peak factor, probability distribution, and power spectrum.

Tropical cyclones are characterized by the asymmetric helical flow structure and complex turbulence driven mechanism (both convective and mechanical turbulence). The spatial distribution of the flow structure also varies significantly in the footprint of tropical cyclones. Due to the limited measurements in the lower boundary layer of tropical cyclones, a basic premise of the existing codes and standards is that the turbulent wind characteristics in tropical cyclones are similar to those observed in the boundary layer winds of extratropical storms. However, it is well known that the downward transport of convective cells generated at higher levels together with the boundary layer rolls could modulate the wind structure and turbulence in the lower tropical cyclone boundary layer. These thermodynamics-related activities may lead to the turbulent wind characteristics of the hurricanes/typhoons different from those of the extratropical winds [2–4].

A direct and reliable approach to examine the turbulent wind characteristics is based on the field observations in the paths of landfalling tropical cyclones. Thus, a number of field measurement programs were initiated in the tropical cyclone-prone regions to monitor the hurricane/typhoon winds [5–14]. The field-measured wind characteristics from different observation stations for various tropical cyclones are quite scattered and hard-to-reach unified conclusions to guide the wind-resistance design of buildings and structures in the cyclone-prone regions. The inconsistent wind characteristics of tropical cyclones may be attributed in part to the complexity of turbulence driven mechanisms, e.g., shear (namely roll and streak structures near the surface), convection, rotation, blocking and sheltering effects at the boundary layer, and also the interactive motions of multi-scale eddies in the flow fields of tropical cyclones [15,16]. On the other hand, the underlying surface and the employed schemes to obtain the turbulence parameters may also significantly influence on the variability of wind characteristics. Since the tracks of tropical cyclones are random, most of the field observations were conducted by installing anemometers and accelerometers on structures or observation towers, which were built in the regions frequently attacked by tropical cyclones. The measured turbulent wind characteristics from these observation stations are quite different from one another because of the underlying surrounding terrain conditions and the lack of well-established guidelines for an appropriate documentation of the near surface wind filed in tropical cyclones. In the China wind codes, wind characteristics are specified over standard terrain with roughness length of 0.03 m, averaging time of 10 min for mean wind and duration time of 3 s for gusty wind at 10 m height. Accordingly, it is essential to convert the turbulence characteristics obtained from various stations to the standard terrain and investigate the wind nature in a unified analytical framework. The "standardized" wind characteristics due to their universality could be useful in instructing the structural design in cyclone-prone regions.

This study first presented an analytical framework in which the mean wind speed, turbulence intensity, integral scale, gust factor, and peak factor measured at various terrains, heights and averaging times were properly standardized. Then, the typhoon and hurricane wind data analyzed here were briefly described. Finally, field-observed turbulent wind characteristics of four typhoons and three hurricanes were converted to the standard condition and comparatively investigated. The standardized results were also compared with the corresponding recommendations in ASCE7-10 [17] and AIJ-RLB-2004 [18]. The difference of the wind characteristics in hurricanes and typhoons were also highlighted.

#### **2. Analytical Framework**

In this section, an analytical framework will be proposed in which the turbulent characteristics measured in various terrain conditions, heights and averaging times could be converted to a standard station. The standardization of the wind characteristics is based on the assumption of the equilibrium boundary-layer theory [19]. The atmosphere stratification in the boundary layer of tropical cyclones is assumed to be neutrally stable, which implies that the turbulent structure within this region is driven mainly by the local surface roughness effects [20]. In the non-equilibrium boundary-layer, this analytical framework may need further investigations subjecting to specific terrain conditions.

#### *2.1. Mean Wind Speed*

To analyze the wind characteristics, an essential step is to convert the wind speeds measured at different station conditions (i.e., various exposures, heights, and averaging times) to the standard condition. The standardization of the mean wind speed in this study follows the three steps: (1) determine the exposure type of the observation station; (2) calculate the gradient wind speeds over the observation exposure; and (3) calculate the mean wind speed at the reference height (10 m high) over standard exposure (open flat terrain) by assuming the gradient wind speeds are equal at the gradient height over different exposures.

#### 2.1.1. Logarithmic Law Wind Profile

Normally, the observations by Global Position System (GPS) dropsonde and Doppler radar show that the variation of mean wind speed with height follows the logarithmic law in the lower part of tropical cyclone boundary layer [21–25]. Thus, the logarithmic law can be used to describe lower boundary layer and the outer-vortex regions of a tropical cyclone:

$$\mathcal{U}L\_{\mathfrak{s}}(z\_{\mathfrak{s}}) = \frac{u\_{\mathfrak{s}\_{\mathfrak{s}}}}{k} \ln \left(\frac{z\_{\mathfrak{s}}}{z\_{0\mathfrak{s}}}\right) \tag{1}$$

where *Us*(*zs*) represents the mean wind speed at height *zs* over the standard exposure. Specifically, the standard exposure in this study corresponds to the roughness length *z*0*s* = 0.03 m, the reference height is 10 m, and the time scale for the average value is 10 min. *u*∗*s* denotes the friction velocity over the standard exposure, and *k* ≈ 0.40 is the von Kármán constant.

According to Equation (1), the key procedure to standardize the mean wind speed is to determine the relationship of the friction velocities in various terrains. As all anemometers used in this study are set between 10–60 m height, it is reasonable to assume that the friction velocity in the lower tropical cyclone boundary layer is a height-independent constant [25–27]. Based on the assumption of local equilibrium conditions, the transition model in Engineering Sciences Data Unit (ESDU) [28], which has been applied to convert the 3 s peak speed over open-terrain and the 1-min mean wind speed above open water in hurricane by Simiu et al. [29], is employed:

$$\frac{\mu\_{\ast\_{S}}}{\mu\_{\ast\_{M}}} = \frac{\ln\left(\frac{10^5}{z\_{0\_s}}\right)}{\ln\left(\frac{10^5}{z\_{0\_{0}}}\right)},\tag{2}$$

where *u*∗*m* is the friction velocity over the field measured exposure with roughness length of *z*0*m*. Then the relation of mean wind speeds with different terrains can be accordingly expressed as:

$$\frac{dL\_s(z\_s)}{dL\_m(z\_m)} = \frac{\ln\left(\frac{10^5}{z\_{0s}}\right)}{\ln\left(\frac{10^5}{z\_{0m}}\right)\ln\left(\frac{z\_m}{z\_{0m}}\right)}\tag{3}$$

where *Um*(*zm*) is the mean wind speed measured at experiment station with height *zm* and roughness *z*0*m*. In this model, the gradient balance assumption, which has been demonstrated to be valid at a sufficiently high altitude [30], is adopted. The super-gradient flows were observed in the boundary layer of some tropical cyclones, however, it is not systematic in tropical cyclones, especially in overland conditions [31,32]. The case of super-gradient flows in tropical cyclones will discussed in next section.

The relation among mean wind speeds of various averaging times can be expressed as [33]:

$$\mathcal{U}\_{\pi}(z) = \mathcal{U}\_{3600}(z) \left[ 1 + \frac{\beta^{0.5} c(\pi)}{2.5 \ln\left(\frac{z}{z\_0}\right)} \right] \tag{4}$$

where τ denotes the averaging time; *U*τ(*z*) and *U*3600(*z*) are respectively τ-s mean and 1-h mean wind speeds; β represents the ratio of the fluctuating wind speed variance to the square of friction velocity; *c*(τ) is an averaging time-related parameter that determined by statistical characteristics of wind speed measurements.

#### 2.1.2. Super-Gradient Wind Profile

The field measurements show the existence of super-gradient wind over ocean surface and the sea land transition regions in tropical cyclones, and the variation of mean wind with height following a logarithmic-quadratic profile [34]. Based on the field measurements in hurricanes over land and ocean surface, Snaiki and Wu [35] proposed a semi-empirical model to depict the mean wind profile. As the empirical model is convenient and accurate, it is adopted here to convert the mean wind speed in landfalling typhoons. The power law-based wind profile is used as follows:

$$\mathcal{U}I\_s(z\_s) = \mathcal{U}I\_{10\_s} \left[ \left(\frac{z\_s}{10}\right)^{\alpha\_s} + \eta\_1 \sin\left(\frac{z\_s}{\delta\_s}\right) \exp\left(-\frac{z\_s}{\delta\_s}\right) \right] \tag{5}$$

where *Us*(*zs*) and *U*10*s* are the mean wind speed at height *zs* and 10 m over the standard exposure; α*s* is the power law exponent over the standard exposure; δ*s* is the height of the wind maximum over the standard exposure; η1*s* is derived to be:

$$\eta\_{1s} = \frac{\left(\frac{\delta\_s}{10}\right)^{\alpha\_s} \alpha\_s \varepsilon}{\sin 1 - \cos 1'} $$

Analogously, the field measured mean wind over the experiment exposure is:

$$\mathcal{U}l\_{m}(z\_{m}) = \mathcal{U}l\_{10\_{m}} \left[ \left(\frac{z\_{m}}{10}\right)^{a\_{m}} + \eta\_{1m} \sin\left(\frac{z\_{m}}{\delta\_{m}}\right) \exp\left(-\frac{z\_{m}}{\delta\_{m}}\right) \right] \tag{6}$$

where *Um*(*zm*) and *U*10*m* are the mean wind speed at height *zm* and 10 m over the measured exposure; α*m* is the power law exponent over the measured exposure; δ*m* is the height of the wind maximum over the measured exposure; η1*s* is:

$$
\eta\_{1m} = \frac{\left(\frac{\delta\_m}{10}\right)^{\alpha\_m} \alpha\_m e}{\sin 1 - \cos 1}
$$

By adopting the assumption that the wind speeds at the wind maximum height ( δ*s* and δ*m*) are equal, the following expression can be deduced:

$$\frac{dL\_s(z)}{dL\_m(z)} = \frac{\left(\frac{\delta\_m}{10}\right)^{a\_m} \left[ \left(\frac{z\_s}{10}\right)^{a\_s} + \eta\_{1s} \sin\left(\frac{z\_s}{\delta\_s}\right) \exp\left(-\frac{z\_s}{\delta\_s}\right) \right]}{\left(\frac{\delta\_s}{10}\right)^{a\_s} \left[ \left(\frac{z\_m}{10}\right)^{a\_m} + \eta\_{1m} \sin\left(\frac{z\_m}{\delta\_m}\right) \exp\left(-\frac{z\_m}{\delta\_m}\right) \right]},\tag{7}$$

Equation (7) could be used to convert the field wind speeds to the standard exposures in the tropical cyclones with super-gradient flow. Actually when the wind speed measured in the lower regions following the logarithmic law, the Equation (7) will merge into the logarithmic law or power law.

#### *2.2. Turbulence Intensity*

It is conventional to treat the turbulence ratio (the ratio of the standard deviation of longitudinal wind velocity component σ*u* to the friction velocity *u*∗) as terrain-independent in the equilibrium boundary layer [36,37]. On the other hand, Harris and Deaves [38] proposed an empirical model to consider the variation of turbulence ratio with height as:

$$\frac{\sigma\_{\rm II}}{\mu\_{\star}} = 2.63\eta \left[ 0.538 + 0.090 \ln \left( \frac{z}{z\_0} \right) \right]^{\eta^{16}} \text{ .} \tag{8}$$

where η = 1 − *z*/*h*; *h* = *u*∗/(6 *f*); *f* = 1.458 × 10−4 sin φ is the Coriolis parameter; and φ denotes the latitude of the observation site. Due to the fact that the derivation of Equation (8) was partly based on non-equilibrium-condition data, the estimation of maximum turbulence ratio, [σ*u*/*u*∗]*max* obtained from Equation (8) is dependent on the terrain roughness length, which is in contradiction to the equilibrium assumption. To correct this issue, the empirical variation of σ*u*/*u*∗ with respect to terrain roughness is introduced in ESDU [39] to obtain approximately a constant [σ*u*/*u*∗]*max* for various terrain roughness lengths:

$$
\left[\frac{\sigma\_u}{\mu\_\*}\right](z\_0) = 1 + 0.156 \ln\left(\frac{\mu\_\*}{fz\_0}\right) \tag{9}
$$

Since the field measurements give [σ*u*/*u*∗]*max* = 2.85 for the terrain with a roughness length of 0.03 m, Equation (8) can be corrected by factoring 2.85/ / 1 + 0.156 ln[*u*∗/(*f z*0)]0 , resulting in an improved model to calculate turbulence ratio as in ESDU [39]:

$$\frac{\sigma\_{\rm II}}{\mu\_{\rm \ast}} = \frac{7.496 \eta \left[0.538 + 0.090 \ln\left(\frac{z}{z\_0}\right)\right]^{\eta^{16}}}{1 + 0.156 \ln\left(\frac{\mu\_{\rm \ast}}{f^z z\_0}\right)}. \tag{10}$$

As expected, Equation (9) gives a maximum turbulence ratio [σ*u*/*u*∗]*max* of approximately 2.85 for various roughness lengths. However, field measurements show that turbulence ratios in tropical cyclones are usually greater than the values in extratropical storms [12,24,40], making the selection of [σ*u*/*u*∗]*max* = 2.85 inapplicable to tropical cyclones. On the other hand, a height-independent relation between turbulence ratio σ*u*/*u*∗ and underlying surface roughness length *z*0 was proposed by Li et al. [13] based on the analysis of field measurements in typhoons. In this study, this height-independent relation proposed by Li et al. [13] is adopted as:

$$
\left[\frac{\sigma\_u}{\mu\_\*}\right]\_{\text{max}} = 2.72 - 0.25 \log z\_0. \tag{11}
$$

As a result, the turbulence ratio in tropical cyclone will be corrected by multiplying Equation (8) by the following factor:

$$\frac{2.72 - 0.25 \log z\_0}{1 + 0.156 \ln[\mu\_\*/(f z\_0)]}. \tag{12}$$

Then the turbulence ratio could be expressed as:

$$\frac{\sigma\_{\rm II}}{\mu\_{\star}} = \frac{2.63\eta \left[0.538 + 0.009\ln\left(\frac{z}{z\_0}\right)\right]^{\eta^{16}} \left[2.72 - 0.25\log(z\_0)\right]}{1 + 0.156\ln\left(\frac{\mu\_{\rm I}}{f\bar{z}\_0}\right)}.\tag{13}$$

For a standard terrain condition (*z* = 10 m; *z*0 = 0.03 m; assuming *u*∗ = 1 m/s and φ = 25◦), the turbulence ratios estimated by Equations (9) and (11) are 2.55 and 2.78, respectively. In this

study, Equation (13) is utilized to convert the measured turbulence ratio to the standard condition. Accordingly, the longitudinal turbulence intensity in a standard exposure can be calculated as:

$$[I\_u]\_s = \frac{2.63\eta\_s [0.538 + 0.090 \ln\left(\frac{z}{z\_0}\right)]^{\eta\_s^{1.6}} [2.72 - 0.25 \ln(z\_{0t})] [u\_\*]\_s}{\left[1 + 0.156 \ln\left(\frac{u\_\*}{f^z z\_0}\right)\right] [L]\_s}. \tag{14}$$

#### *2.3. Integral Scale*

The approach of integrating correlation function by invoking Taylor's hypothesis is frequently used to estimate the integral scale as it has a clear physical meaning [24,33]:

$$L\_u^x = \frac{\mathcal{U}}{\sigma\_u^2} \int\_0^{\mathcal{R}\_{\text{int}} = 0.05r\_u} R\_{uu}(\tau) \,\text{d}\tau,\tag{15}$$

where *Ruu* is the autocorrelation function of the longitudinal fluctuating component.

The integration of autocorrelation function, however, usually overestimates the value of integral scale and will result in a deviation of inertial sub-range in the estimated von Kármán-type spectrum [41] compared to that in the field-measured spectra. To improve the accuracy, Harris and Deaves [38] suggested the following model to estimate the longitudinal integral scale:

$$L\_u^x = \frac{A^{\frac{3}{2}} \left(\frac{\sigma\_u}{u\_\*}\right)^3 z}{2.5 K\_z^{\frac{3}{2}} \left(1 - \frac{z}{h}\right)^2 \left(1 + 5.75 \frac{z}{h}\right)},\tag{16}$$

where *z* is the height from the ground.

$$A = 0.115 \left( 1 + 0.315 \eta^6 \right)^{\frac{2}{3}} \tag{17}$$

and

$$K\_z = 0.19 - (0.19 - K\_0) \exp\left[ -B \left( \frac{z}{h} \right)^N \right] \tag{18}$$

in which

$$K\_0 = \frac{0.39}{R\_\phi^{0.11'}}\tag{19}$$

$$B = 24R\_o^{0.155} \text{\AA} \tag{20}$$

$$N = 1.24 R\_{\vartheta} \stackrel{0.008}{\text{ \\$}},\tag{21}$$

$$R\_{\theta} = \frac{u\_{\ast}}{f z\_{0}}.\tag{22}$$

The longitudinal integral scale over standard exposure [*Lx u*]*s* can be estimated according to Equation (16) by introducing the corresponding values of [*u*∗]*s* and [*z*0]*s*.

#### *2.4. Peak Factor*

Peak factor *gu* is defined as the ratio of maximum wind speed fluctuation in a duration τ to the standard deviation of the fluctuating wind speed within an observation period of *T:*

$$g\_{\boldsymbol{u}}(\boldsymbol{\tau},T) = \frac{\max[\boldsymbol{u}(\boldsymbol{\tau},T)]}{\sigma\_{\boldsymbol{u}}(\boldsymbol{\tau},T)} \frac{\sigma\_{\boldsymbol{u}}(\boldsymbol{\tau},T)}{\sigma\_{\boldsymbol{u}}}.\tag{23}$$

For a stationary stochastic process following Gaussian distribution, the peak factor with τ → 0 and *T* ≥ 3600 *s* could be calculated as [42]:

$$g\_{\mu}(\tau, T) = \sqrt{2\ln[\nu(\tau, T)T]} + \frac{0.5772}{\sqrt{2\ln[\nu(\tau, T)T]}},\tag{24}$$

where υ(τ, *T*) is the zero up-crossing rate. It can be calculated by [42,43]:

$$\nu^2(\tau, T) = \frac{\int\_0^\infty n^2 S\_\mathrm{u}(n) \chi^2(n, \tau, T) \mathrm{d}n}{\int\_0^\infty S\_\mathrm{u}(n) \chi^2(n, \tau, T) \mathrm{d}n},\tag{25}$$

in which *Su*(*n*) represents the wind velocity spectrum; *n* denotes the frequency in Hertz, and χ2(*n*, τ, *T*) is a filter function used to consider the influence of sampling frequency, averaging time and response characteristics of the anemometer. In this study, the von Kármán-type spectrum is employed as:

$$\frac{nS\_{\rm II}(n)}{\sigma\_{\rm u}^{2}} = \frac{4\left(\frac{nL\_{\rm u}^{x}}{\ell I}\right)}{\left[1 + 70.8\left(\frac{nL\_{\rm u}^{x}}{\ell I}\right)^{2}\right]^{\frac{5}{6}}},\tag{26}$$

The filter function is chosen as following for sonic anemometers [44]:

$$\chi^2(n,\tau,T) = \left[\frac{\sin(\pi n \tau)}{\pi n \tau}\right]^2 - \left[\frac{\sin(\pi n T)}{\pi n T}\right]^2,\tag{27}$$

For propeller anemometers, the following filter function, which takes the mechanical features of propeller anemometers into consideration, is adopted [43]:

$$\chi^2(n,\tau,T) = \left\{ \left[\frac{\sin(\pi n \tau)}{\pi n \tau}\right]^2 - \left[\frac{\sin(\pi n T)}{\pi n \tau}\right]^2 \right\} \frac{1}{1 + \left(\frac{2\pi n \lambda}{lI}\right)^2},\tag{28}$$

where λ is the distance constant of the propeller anemometer.

Equation (24) is valid for calculating the average of instantaneous peak factor ( τ → 0) from a long enough wind speed record (e.g., *T* ≥ 1 h). With a finite averaging time, τ, and a finite observation period, *T*, the estimation of standard deviation in Equation (26) could be biased since the measured spectrum is truncated in both high-frequency and low-frequency regions, and might eventually lead to an inaccurate estimation of the peak factor. In the case that these conditions are not satisfied, the following relation is necessary to be introduced to consider the effects of the variance reduction due to the truncation of the velocity spectrum:

$$\frac{\sigma\_{\rm u}(\tau, T)}{\sigma\_{\rm u}} = \frac{\int\_{0}^{\infty} S\_{\rm u}(n) \chi^{2}(n, \tau, T) \mathrm{d}n}{\int\_{0}^{\infty} S\_{\rm u}(n) \mathrm{d}n},\tag{29}$$

The 3-s peak factor, [*gu*] *s*, in time scale [*T*]*s* = 600 s in the standard terrain can be estimated by introducing *Us* and *Lx us*, which were respectively calculated through Equations (3) or (7) and (16).

#### *2.5. Gust Factor*

Gust factor, *Gu*(τ, *T*), herein is defined as the ratio of gust speed with gust duration τ to the mean wind speed *U*(*T*) with an observation period of *T*:

$$G\_{\rm u}(\tau, T) = 1 + \frac{\max[\boldsymbol{u}(\tau, T)]}{\sigma\_{\rm u}(\tau, T)} \frac{\sigma\_{\rm u}(\tau, T)}{\sigma\_{\rm u}(\Delta t, T)} \frac{\sigma\_{\rm u}(\Delta t, T)}{\mathcal{U}(T)},\tag{30}$$

where Δ*t* is the sampling interval.

Substituting the peak factor and turbulence intensity into the corresponding terms of Equation (30), the gust factor can be re-expressed as:

$$G\_{\rm u}(\tau, T) = 1 + g\_{\rm u}(\tau, T) I\_{\rm u} \frac{\sigma\_{\rm u}(\tau, T)}{\sigma\_{\rm u}(\Delta t, T)},\tag{31}$$

where σ*u*(τ, *T*)/σ*u*(Δ*t*, *T*) can be calculated by:

$$\frac{\sigma\_{\rm u}(\tau, T)}{\sigma\_{\rm u}(\Delta t, T)} = \frac{\int\_0^\infty S\_{\rm u}(n) \chi^2(n, \tau, T) \mathrm{d}n}{\int\_0^\infty S\_{\rm u}(n) \chi^2(n, \Delta t, T) \mathrm{d}n},\tag{32}$$

As a result, the gust factor in the standard exposure, [*Gu*(τ, *T*)]*s*, can be estimated by:

$$[G\_{\rm u}(\tau, T)]\_{\rm s} = 1 + [g\_{\rm u}(\tau, T)]\_{\rm s} [I\_{\rm u}]\_{\rm s} \left[\frac{\sigma\_{\rm u}(\tau, T)}{\sigma\_{\rm u}(\Delta t, T)}\right]\_{\rm s} \tag{33}$$

#### **3. Data Sources**

#### *3.1. Tropical Cyclones and Instruments*

In this study, the data of four typhoons (0601 typhoon Chanchu, 0606 typhoon Prapiroon, 0812 typhoon Nuri, and 0814 typhoon Hagupit) and three hurricanes (0504 hurricane Katrina, 0510 hurricane Rita, and 0512 hurricane Wilma) were comparatively analyzed. The detailed descriptions of the observation site exposures and the observation tower configurations for the four typhoons and three hurricanes were presented in Li et al. [16] and Masters et al. [12], respectively. The GPS coordinates of the observation stations were listed in Table 1. As the latitudes of all observation stations are around 25◦, the latitude of the standard condition is set to be 25◦ for the convenient of calculation.

**Table 1.** The GPS coordinates of observation towers. Reproduced with permission from [16], Elsevier, 2019.

It is noted that the distance constant λ of the propeller anemometer is an important factor to calculate the peak factor and gust factor as this type of anemometer mechanically filters the amplitudes of gusts with wavelengths less than 2πλ due to the mechanical limitations [10]. In this study, the propeller anemometers of models R.M. Young 05103L and R.M. Young 27106R were respectively used to measure typhoons and hurricanes. The specifications of these two propeller anemometers are listed in Table 2. Based on the parameters in Table 2, the data measured by propeller anemometers were corrected according to Equation (24). In addition to propeller anemometers, the sonic anemometers were also utilized in the field measurement of typhoons. Specifically, two 3-D ultrasonic anemometers (WindMaster™ Pro., Gill Instruments Ltd., Lymington, UK) were installed on tower RBT and one 3-D ultrasonic anemometer (HD2003, Delta Ohm Srl, Selvazzano Dentro, Italy) were setup on tower OT for the measurements of Typhoon Chanchu; one HD2003 anemometer were installed on tower BT to acquire data from Typhoon Prapiroon, and towers DIT and ST were equipped with Gill WindMasterTM Pro. anemometers. The specifications of the utilized sonic anemometers are also listed in Table 2. For the data obtained from sonic anemometers, the filter function presented by Equation (23) were used for the correction.

**Table 2.** Specifications of anemometers.

#### *3.2. Data Quality Control and Data Source*

Tropical cyclones are characterized by strong winds accompanied by torrential rain, ocean waves, and storm surge. The representative wind records are usually located in the eyewall regions of tropical cyclones. In the field measurements of tropical cyclones, however, the anemometers in heavy rain bands usually present some spikes and errors. Therefore, the data quality-control procedure is a necessary step before the analysis of the wind characteristics. In this study, the data quality-control schemes and the criteria for the selection of samples were referred to Li et al. [13]. Specifically, the spikes and errors in the data were first identified and replaced by the five-point weighted averages. Then, the reverse arrangement test [45] and run test [46] with a 95% significance level were employed to test the stationarity of the recorded winds. The datasets that failed to pass both two types of stationary tests were removed from the analysis. The stationarity test ensured that the analyzed data could satisfy the local equilibrium boundary-layer assumption, where the friction velocity is independent of the location in the along-wind direction and the Reynolds number [37].

Table 3 briefly summarized the datasets utilized in this study, together with the observation heights, types of the anemometers used, number of the runs in each group (sample size) and average of the 10 min mean wind speeds. The detailed analysis of the turbulent wind characteristics of the original datasets both in typhoons and hurricanes were presented in Li et al. [16].

#### **4. Results and Discussions**

The selected datasets from the four typhoons and the three hurricanes summarized in the preceding section were investigated in the analytical framework presented in Section 2. Specifically, both the datasets in typhoons and hurricanes were converted to the standard exposure with a roughness length of 0.03 m at 10 m height and an observation time scale of 10 min. The latitude of the standard terrain is assumed to be 25◦.

#### *4.1. Turbulence Intensity*

The turbulence intensities were extracted based on the analytical framework and shown in Figures 1a and 2a respectively for investigated typhoons and hurricanes. The corresponding probability density functions (PDFs) are shown in Figures 1b and 2b. For typhoons, the average value of longitudinal turbulence intensities is 0.1952 and the standard deviation is 0.0032. For hurricanes, the average value of longitudinal turbulence intensities is 0.1906 and the standard deviation is 0.0022. The longitudinal turbulence intensity of these four typhoons presents slightly higher values in terms of both mean and standard deviation compared to those three hurricanes. This observation can be further demonstrated by comparing Figures 1b and 2b, where the probability distribution of the longitudinal turbulence intensities in both typhoons and hurricanes follows the normal distribution quite well. One can easily conclude that in these four typhoons the turbulence intensity has a larger value than that in those three hurricanes under the same probability of exceedance.

In ASCE7-10, the longitudinal turbulence intensity is given by:

$$I\_{\rm II} = c \left(\frac{10}{z}\right)^{\frac{1}{b}},\tag{34}$$

where *c* equals to 0.30, 0.20, and 0.15 for category B, C (corresponding to the standard terrain in this study), and D exposures, respectively. In AIJ-RLB-2004 code, the longitudinal turbulence intensity over flat terrain categories is given by:

$$I\_{\rm li} = 0.1 \left(\frac{z}{z\_G}\right)^{-\alpha - 0.05},\tag{35}$$

in which α and *zG* are parameters reflecting the category of exposures. In category II exposure, which is the closest to the standard terrain in this study, α and *zG* are respectively 0.15 and 350. The longitudinal turbulence intensities obtained from ASCE7-10 and AIJ-RLB-2004 are 0.2000 and 0.2036, respectively, for the standard exposure. As depicted in Figures 1a and 2a, both ASCE7-10 and AIJ-RLB-2004 present a slightly higher estimation of longitudinal turbulence intensity for hurricanes and typhoons. Generally, the estimation of ASCE7-10 is relatively better compared to that of AIJ-RLB-2004.

**Figure 1.** Turbulence intensities and their probability density functions (PDF) for typhoons.

**Figure 2.** Turbulence intensities and their PDF for hurricanes.

#### *4.2. Integral Scale*

The integral scales and the corresponding probability distributions obtained are presented in Figures 3 and 4 for typhoons and hurricanes, respectively. The longitudinal integral scale in typhoons has an average value of 146.4482 and a standard deviation of 9.2143. The length scales extracted from hurricane measurements have a slightly higher average value of 157.5796 and a significantly lower standard deviation of 3.8106 compared to those of typhoons. As shown in Figures 3a and 4a, the range of integral scales extracted from typhoon measurements is significantly larger than that of hurricanes, which can be better illustrated by the probability distributions of longitudinal integral scales in typhoons and hurricanes. The probability distribution of typhoons follows the Weibull distribution with scale parameter of 150.368 and shape parameter of 20.0454, while the probability distribution of longitudinal integral scales for hurricanes follows the generalized extreme value distribution with scale parameter of −0.3259, shape parameter of 3.8960 and location factor of 156.311. Compared with Figure 4b, where the observed hurricane length scale shows a narrow distribution, the probability distribution of the observed typhoon integral scale of Figure 3b has a significantly wider range. This phenomenon may be in part attributed to the exposures of the observation station for the original datasets. The observation station in those three hurricanes are located in homogeneous open flat terrain as stated in Masters et al. [12]. However, in the observation of those four typhoons, the exposures of the measured stations are a little bit inhomogeneous. Another influence could be ascribed to the differences of turbulent structures of those typhoons and hurricanes. As noted in Li et al. [16], at the same roughness regime, the field measured integral scales in these four typhoons were greater than that in those three hurricanes. The different distributions of the observed hurricane and typhoon length scales might indicate that energy-containing eddies in the observed typhoons have various representative length-scales while those of the observed hurricanes are concentrated around the mean value. The multiple-scale eddy interactions in typhoons and hurricanes need further investigations before fully understanding the observed difference.

In ASCE7-10, the longitudinal integral scale is computed by:

$$L\_u^{\mathbf{x}} = l \left(\frac{z}{10}\right)^{\overline{7}},\tag{36}$$

where *l* and ε are respectively 152.4 m and 0.2 for category C exposure (standard exposure in this study). In AIJ-RLB-2004 code, the turbulence integral scale is defined independently of the terrain categories and is given by:

$$L\_u^x = \begin{cases} 100 \left(\frac{z}{30}\right)^{0.5} & 30 \text{ m} < z < z\_G\\ 100 & z \le 30 \text{ m} \end{cases} \tag{37}$$

where *zG* equals to 350 m for category II exposure, corresponding to the standard exposure in this study. Accordingly, the longitudinal integral scales obtained from ASCE7-10 and AIJ-RLB-2004 are respectively, 152.4 m and 100 m, for the standard exposure. It is noted that ASCE7-10 presents a reasonable estimation of the longitudinal integral scales for both typhoons and hurricanes. However, AIJ-RLB-2004 underestimates the longitudinal integral scales for both typhoons and hurricanes, suggesting that the usage of AIJ-RLB-2004 may lead to an inaccurate estimation of the power spectrum.

**Figure 3.** Longitudinal integral scales and their PDF for typhoons.

**Figure 4.** Longitudinal integral scales and their PDF for hurricanes.

#### *4.3. Peak Factor*

The peak factor is usually utilized for the estimation of gust factor, which plays an important role in determining the wind load on structures [43]. The estimated peak factors of typhoons and hurricanes are respectively presented in Figures 5 and 6, together with the corresponding probability distributions. For typhoons, the average value of the peak factor is 2.5211 and the standard deviation is 0.0198. The fitted PDF is shown in Figure 5b. For hurricanes, the average value of peak factor is 2.5123, slightly smaller than that of typhoons, and the standard deviation of peak factor is 0.0102, significantly smaller than that of typhoons. The probability distribution of peak factors in hurricanes follows t location-scale distribution with location parameter of 2.5174, scale parameter of 0.0022 and shape parameter of 0.9845.

Neither the expression nor the value of peak factor is explicitly prescribed in ASCE7-10. By matching the gust factor over open terrain ( *Gu* = 1.53) and the turbulence intensity *Iu* of

0.2 in Equation (34), a peak factor of 2.65 could be obtained. It should be noted that this calculation is based on an averaging time of 1 h. For a duration of *tg* , the gust wind speed can be expressed as:

$$
\mathcal{U}\{t\_{\mathcal{G}'}T\} = \overline{\mathcal{U}}(T) + \mathcal{g}\_{\mathfrak{u}}(t\_{\mathcal{G}'}T)\sigma\_{\mathfrak{u}}.\tag{38}
$$

Suppose the turbulent wind fluctuations follow the Gaussian distribution, the peak factor will be associated with the exceedance probability of the standard normal distribution. The probability of exceedance of wind gust with a duration of *tg* within an observation period of *T* could be calculated as [47,48]:

$$P\left[\mathcal{U} > \mathcal{U}\left(t\_{\mathcal{S}'}\right)\right] = \frac{t\_{\mathcal{S}}}{T}.\tag{39}$$

Thus, the peak factor should satisfy:

$$\log \left( t\_{\mathcal{S}'} \, T \right) = \Phi^{-1} \left( 1 - \frac{t\_{\mathcal{S}}}{T} \right) . \tag{40}$$

With the gust duration of 3 s and the averaging time of 10 min, the peak factor is around 2.575, which is slightly higher than the measured values in typhoons and hurricanes. AIJ-RLB-2004 carries out a performance-based wind resistant design procedure. Accordingly, the peak factor is included in the required performance of wind load level and the return period of wind speed. Hence, the comparison of the measurement results with AIJ-RLB-2004 is not discussed here.

**Figure 5.** Peak factors and their PDF for typhoons.

**Figure 6.** Peak factors and their PDF for hurricanes.

#### *4.4. Gust Factor*

The gust factor in steady wind conditions depends on several wind characteristics, such as the intensity and integral scale, hence, it is a basic representation of the dynamic properties of wind loads [49]. Figure 7 depicts the gust factors and the corresponding probability distribution of typhoon winds obtained based on the unified analysis framework. For typhoons, the average value of gust factors is 1.4919 and the standard deviation is 0.0069. The gust factors of hurricanes are presented in Figure 8a, where the gust factors have a mean of 1.4787 and a standard deviation of 0.0071. The mean of gust factors for typhoons are slightly higher than that for hurricanes, while the standard deviations of gust factors for typhoons and hurricanes are almost identical with similar probability distribution shapes. The probability distribution of gust factors for typhoons follows the extreme value distribution with a location parameter of 1.4948 and a scale parameter of 0.0051, while the probability distribution of gust factors for hurricanes follows the generalized extreme value distribution with shape parameter of −0.5535, location parameter of 1.4771, and scale parameter of 0.0076.

In ASCE7-10, the calculation of gust factor is referenced to the gust factor curve proposed by Durst [50]. The averaging time of the mean wind speed is 1 h in the Durst gust factor curve, while the gust factor is calculated based on an averaging time of 10 min in this study. Therefore, the conversion scheme for the gust factors with different averaging times presented in Vickery and Skerlj [20] was utilized here. Gust factor with a duration of 3 s and an averaging time of 10 min can be expressed as:

$$G\_{\mu} = 1 + (SLI)[SD(600, 3)],\tag{41}$$

where *SU* is the value of the standard normal deviation associated with the exceedance probability of 0.5% and equals to 2.575. The *SD*(600, 3) could be estimated by the following formula:

$$SD(600, \,\,3) = \left[SD^2(3600, \,\,3) - SD^2(3600, \,\,600)\right]^{1/2},\tag{42}$$

where *SD*(3600, 3) and *SD*(3600, 600) can be interpolated as indicated in Vickery and Skerlj [20] and equal to 0.1617 and 0.0650, respectively. The gust factor with duration of 3 s and averaging time of 10 min based on Equations (37) and (38) is around 1.3814, which indicates that the gust factors for both typhoons and hurricanes are greater than those for extratropical storms.

**Figure 7.** Gust factors for typhoons.

**Figure 8.** Gust factors for hurricanes.

#### **5. Concluding Remarks**

An analytical framework was introduced in this study to standardize the turbulent wind characteristics, namely turbulence intensity, integral scale, peak factor and gust factor for various terrain conditions, heights, and averaging times in tropical cyclones. This analytical framework is based on the equilibrium boundary-layer theory and the assumption that the lower tropical cyclone boundary layer is neutrally stable. Field-measured data of the four typhoons and three hurricanes were standardized to the reference exposure with roughness length of 0.03 m, height of 10 m, and averaging time of 10 min, and then utilized for the extraction of wind characteristics under standard exposure. The differences of obtained wind characteristics between typhoons and hurricanes were highlighted, which may be attributed to the basins and latitudes of the genesis of hurricanes and typhoons, the influence of local topography and sea-land transition zone and the differences in turbulent flow structures of typhoons and hurricanes that need further investigations. More specifically, the wind characteristics of these observed typhoons typically present larger values compared to those of observed hurricanes, except for the turbulence integral scale. The turbulence integral lengths of typhoons have a wider distribution compared with those of hurricanes. The obtained turbulent wind characteristics based on the unified analysis framework were comparatively investigated together with the recommendations in ASCE7-10 and AIJ-RLB-2004. The difference between the standardized turbulent characteristics and the corresponding suggested values in the standards (ASCE7-10 and AIJ-RLB-2004) indicates that the tropical cyclone-induced wind loads need be taken into consideration in standards for tropical cyclone-prone regions. It is noted that the ASCE7-10 presents good estimations of the longitudinal turbulence intensity and integral scale for both typhoons and hurricanes, while the peak factor was slightly overestimated and the gust factor was underestimated. The AIJ-RLB-2004 makes a slightly higher estimation of the longitudinal turbulence intensity and a lower estimation of the longitudinal integral scale for both typhoons and hurricanes. The potential reason may be ascribed to the limitation of datasets which used to specify the wind characteristics, although it includes both tropical and extratropical winds. As noted in the AIJ-RLB-2004, the integral scale was treated to be terrain independent. However, the scales of wind eddies are strongly affected by the local roughness.

**Author Contributions:** Conceptualization, L.L. and T.W.; methodology, L.L. and H.Z.; software, Y.Z. and H.W.; validation, T.W. and H.W.; formal analysis, Y.Z.; investigation, Y.Z. and H.W.; resources, X.H.; data curation, Y.Z.; writing—original draft preparation, L.L. and T.W.; writing—review and editing, L.L., X.H. and T.W.; visualization, Y.Z.; supervision, L.L., H.Z., X.H. and T.W.; project administration, L.L.; funding acquisition, L.L. and T.W.

**Funding:** This research was funded by National Natural Science Foundation of China (Grant No. 51778373), the Knowledge Innovation Project of Shenzhen (Grant No. JCYJ20170302143625006), Natural Science Foundation of SZU (Grant no. 082017), Natural Science Foundation Grant # CMMI 15-37431 and National Key Research and Development Program of China (2017YFB1201204).

**Acknowledgments:** The authors gratefully acknowledge K. Gurley of the University of Florida for providing the hurricane data and Lili Song of China Meteorological Administration for providing the typhoon data for this study.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Full-Scale Train Derailment Testing and Analysis of Post-Derailment Behavior of Casting Bogie**

#### **Hyun-Ung Bae 1, Jiho Moon 2, Seung-Jae Lim 3, Jong-Chan Park 4 and Nam-Hyoung Lim 4,\***

Received: 10 November 2019; Accepted: 18 December 2019; Published: 19 December 2019

**Abstract:** In this study, a full-scale train bogie derailment test was conducted. For this, test methodologies to describe the wheel-climbing derailment of the train bogie and to obtain accurate test data were proposed. The derailment test was performed with the casting bogie for a freight train and a Rheda 2000 concrete track. Two different derailment velocities (28.08 km/h and 55.05 km/h) were considered. From the test, it was found that humps in the concrete track affected the post-derailment behavior of the bogie when the derailment velocity was 28.08 km/h. For a higher derailment velocity (55.05 km/h), significant lateral movement of the derailed bogie was observed. This lateral movement was first controlled by wheel–rail contact, followed by contact with the containment wall. Finally, the train was returned to the track center.

**Keywords:** train derailment; derailment containment provisions; collision testing; post-derailment behavior

#### **1. Introduction**

In Korea, there were 33 train accident cases in total during the five-year period 2012–2016. Among these, derailment accidents accounted for 78.8% (26 cases) [1], and derailment occurred more often than other types of train accidents. Derailment accidents can cause catastrophic damage to a community. It is hard to prevent 100% of derailment accidents since there are always unexpected factors that can cause derailment, such as human error and natural disasters. Thus, it is necessary to develop technology to reduce damage due to derailment. This technology can be categorized as derailment protection [2–4].

To reduce the damage from derailment accidents, protection facilities can be installed in the railway track. In Korea, guard rails to prevent derailment are used at sharp curves, bridges, and switches, as shown in Figure 1a. Containment walls are also installed on bridges for high-speed railways (where the minimum speed of the line is 200 km/h), as shown in Figure 1b [5–7].

In European countries, three different types of derailment containment provisions (DCPs) are used (DCP types I, II, and III), as shown in Figure 2 [8]. The guard rail is one example of a DCP type I facility, where the DCP is installed inside the track gauge. The wheel of the derailed train comes into direct contact with the facility. DCP type II is similar to DCP type I, but is installed outside of the track gauge. DCP type III facilities are installed outside of the track, similar to DCP type II. However, they are different from DCP type II since the axis of the wheel or bogie of the derailed train impacts this type of DCP.

**Figure 1.** Examples of derailment protection facilities: (**a**) guard rail and (**b**) containment wall.

**Figure 2.** Concept of derailment containment provisions (DCPs).

To verify the performance of such derailment protection facilities, tests and simulations must be conducted, including investigation of the post-derailment behavior of the train. Some researchers have conducted derailment and post-derailment simulations by using 3D finite element analysis. Researchers in Sweden [9,10] analyzed the post-derailment behavior of the wheel of a derailed train colliding with a concrete railroad sleeper. Some researchers in China [11–13] investigated the post-derailment dynamic behavior of a railway vehicle under earthquake excitations. Also, there have been studies simulating collision with a protective facility after the derailment of a high-speed train [2–4].

Full-scale derailment testing is the most reliable method to evaluate the post-derailment behavior of a train and the performance of derailment protection facilities. However, it is very difficult, and only a few such studies can be found in the literature. Wu et al. [14] performed derailment testing, but it was lab-scale testing and the speed was limited to 16 km/h. The test methods relating to train derailment, such as the derailing method and data filtering, are not well established and must be investigated.

In this study, full-scale train derailment testing was conducted. For realistic simulation of the derailment situation and reliable data acquisition from the test, a derailment device and data acquisition system were proposed. Then, from the test, the post-derailment behavior of the train was analyzed.

#### **2. Experimental Methodologies**

#### *2.1. Test Site and Track*

Full-scale train derailment testing requires a large test site. The test site consisted of acceleration, test, and braking regions, as shown in Figure 3. To increase the speed of the train at the point of derailment, a sufficient acceleration region is needed. Also, adequate test and braking regions should be provided to investigate the post-derailment behavior and to ensure safety, respectively. The lengths of the acceleration, test, and braking regions were 1200 m, 400 m, and 400 m, respectively. In the 400 m test region, a region of 100 m was used to construct a concrete track. The total length of the test line was about 2000 m. Once the target speed of the bogie is reached, a test bogie is released at the end of the acceleration region. The test bogie is derailed at the start of the test region and the post-derailment behavior is observed in the test region. A braking region is provided at the end of the test region. This is a margin region to provide safety after the unexpected behavior of derailed bogie or train. The derailment test site was constructed using a closed railway line to save on costs and replicate actual railway operating conditions.

**Figure 3.** Overview of the test site.

In this study, the focus was on the post-derailment behavior on a concrete track. A Rheda 2000 concrete track, shown in Figure 4, was constructed in the test region where the post-derailment behavior was observed. The Rheda 2000 concrete (ballastless) track was used for the first time in Germany in 2000 as a pilot project on the new rail line between Erfurt and Halle-Leipzig. It was also installed on the high-speed railway in Korea in 2004. Apart for this 100-m concrete track, ballast tracks were used. It should be noted that the derailment containment wall was installed on the left side of the track, as shown in Figure 4. The distance between the wall and the center of the track was 2700 m, considering the geometric condition of the axle of the test vehicle after the derailment.

**Figure 4.** The concrete track used in the test region.

#### *2.2. Test Bogie and Acceleration Method*

In this study, a bogie-level test was conducted. The bogie used in this study was a casting bogie for a freight train, as shown in Figure 5. The total weight of the test bogie was 40.91 kN and the specifications of the test bogie are listed in Table 1. The frame structure was made of three pieces of cast steel. The fixed wheel base and wheel diameter were 1676 mm and 860 mm, respectively. A suspension system with coil springs was used between the bolster and side frame.

**Figure 5.** Test bogie (casting bogie for a freight train).

**Table 1.** Specifications of the test bogie.

In order to accelerate the test bogie, several methods can be used, such as a reverse towing system, push system, or remote-controlled system. In this study, the push system shown in Figure 6 was used. The power car was linked to the test bogie with a connector. The test car was accelerated by the power car. Then, the test car was released after reaching the target speed. For this, the releasing system of the connector and braking system of the power car were designed to be controlled by air pressure and a remote controller, as shown in Figure 7. The power car used in this study is shown in Figure 8. The power car had a traction power of 147 kN.

**Figure 6.** Concept of the push system to accelerate the test bogie.

**Figure 7.** Release and braking system for the test.

**Figure 8.** Power car to accelerate the test bogie.

#### *2.3. Derailment Device*

A device to induce derailment was needed for the test. Wheel-climbing derailment is one of the most frequent types of derailment. In this study, a derailment device that can induce wheel-climbing derailment was developed and installed at the derailment point, as shown in Figure 9 [15]. The device consisted of a wheel entrance part and a derailment part. In the wheel entrance part, there is an upward slope, and the entered wheel flange reaches the same height as the surface of a rail head. Then, the wheel climbs the surface of a rail head in the derailment part, shown in Figure 9c. Thus, the train or bogie is derailed due to wheel-climbing derailment.

**Figure 9.** Derailment device: (**a**) overview, (**b**) wheel entrance part, and (**c**) derailment part.

#### *2.4. Data Acquisition System*

The position of the derailed train bogie and related data have to be accurately measured during the test. The acceleration and angular velocity of the bogie were measured. For this, an accelerometer and an angular velocity sensor were installed on the geometrical bogie center as shown in Figure 10b, where the capacities of the accelerometer and angular velocity sensor were 2000 g and ±1500 deg/s full-scale range, respectively. Large impact forces are expected after derailment. Thus, the sensors must be securely fixed to the test bogie. To measure the impact force, special accelerometers such as MEMS accelerometers [16] can be used. However, it is hard to determine the exact impact point after derailment. Thus, the impact force to the rail was estimated from the acceleration data of the test bogie. A data logger was installed on the test bogie, as shown in Figure 10c. The data logger used was a shock-resistance data logger (high shock rating of 500 g) for collision testing and it has high sampling capabilities (max. sampling rate of 100 k samples/s/channel). The sampling rate of the data was 1/10,000 s (10,000 Hz) in this test.

**Figure 10.** Data acquisition system for the bogie: (**a**) layout, (**b**) accelerometer and angular velocity sensor, and (**c**) a shock-resistance data logger for collision testing.

The velocity of the test bogie at the derailment point is important information for the analysis of post-derailment behavior. Since the test bogie was released before the derailment point, it was necessary to set up an additional velocity measurement system at the derailment point, as shown in Figure 11. In this study, a photoelectric tube speedometer was used. Two lights and two receiver sensors were installed at a specific distance. As the test bogie passes, the passing time between the two lights can be obtained. Then, the initial derailment velocity of the test bogie can be calculated.

**Figure 11.** Speed measuring system: (**a**) system conceptual diagram; (**b**) photoelectric tube speedometer.

In this study, a high-speed camera was also used to investigate the post-derailment behavior of the test bogie in depth. Image data are very useful to evaluate the trace of the bogie after derailment. The high-speed camera used in this study is shown in Figure 12. According to its specifications, it can store 3.6 s at 1000 fps (frame/s) when the image format is 1024 × 1024. In total, three high-speed cameras were used. One camera was used to take the top view with a crane, as shown in Figure 12. The others were used to take the side views. For each second, 500 frames were obtained (500 fps). Also, some ordinary digital cameras were installed to obtain image data from various viewpoints.

**Figure 12.** High-speed cameras: (**a**) side view; (**b**) top view.

#### *2.5. Simultaneous Trigger System*

Trigger systems are used to synchronize several different pieces of measurement equipment. By synchronizing the equipment, data analysis is convenient since the initial time of data recording is the same for all equipment. In this study, the three high-speed cameras and data logger were synchronized by a trigger system with a tape switch, as shown in

The tape switch method is simple and free from spurious operation. The trigger for the high-speed cameras was set up on the test track at the derailment point, as shown in Figure 13a. The trigger for the data logger, shown in Figure 13b, was installed on the test bogie. When the test bogie passes the derailment point, the trigger on the test bogie contacts the trigger for the high-speed camera. Then, data recoding is started. It should be noted that the trigger on the bogie was made of flexible materials since the post-derailment behavior could be affected by high stiffness and contact force between the triggers. Figure 13.

**Figure 13.** Simultaneous trigger system by a tape switch: (**a**) trigger for the high-speed cameras; (**b**) trigger for the data logger.

#### *2.6. Running Diagram of the Power Car and Test Bogie*

It was necessary to determine the running diagram of the power car and test bogie to achieve the target speed of the test bogie and to ensure the safety of the power car after releasing the test bogie. Figure 14 shows a schematic view of the running diagram of the power car and test bogie. The blue solid line and red dashed line represent the running diagrams of the power car and the test bogie, respectively. Before releasing the test bogie, the power car and the test bogie move together. The velocity at the release time is larger than the target speed. After releasing the test bogie, the power car reduces its speed and stops before the derailment point. The test bogie reaches the target velocity at the derailment point and is derailed.

**Figure 14.** A schematic view of the running diagram of the power car and test bogie.

The distance to the release point depends on the target speed and acceleration capacity of the power car. The running diagram was constructed from the results of several preliminary tests. In this study, the target speeds of the test bogie were 30 km/h and 55 km/h. Figure 15 shows a representative running diagram of the power car through the preliminary tests. An acceleration distance of about 480 m was required for the stationary power car to accelerate to about 56 km/h, and a braking distance of about 150 m was required for subsequent stops. The deceleration to the derailment point after the release of the test bogie was approximately 1–2 km/h. Therefore, the running distance of the power car to ensure safety during the experiment with a target speed of 50 to 55 km/h was determined to be 650 m. In addition, an acceleration distance of about 140 m was required for the power car to accelerate to about 32 km/h, and a braking distance of about 90 m was required for subsequent stops. The deceleration to the derailment point after the release of the test bogie was approximately 2–3 km/h. Thus, the running distance of the power car to ensure safety with a target speed of 25 to 30 km/h was determined to be 250 m.

**Figure 15.** Examples of running diagrams of the power car.

#### *2.7. Post-Processing of Data*

It is necessary to conduct a filtering process on raw data since raw data may contain unexpected noise, especially in the case of impact or collision testing. In the case of car collision testing, the data are analyzed after filtering the raw data obtained at a 1000 Hz or 10,000 Hz frequency. Usually, the moving average method is used for the filtering [17,18]. The moving average method is used not only for the smoothing of the data, but also for the evaluation of equivalent static design force. A previous study reported that the moving average method gave a reasonable estimation of the equivalent static design load [19]. Thus, the 50 ms moving average method was applied with 10,000 Hz data sampling in this study.

#### **3. Experimental Results**

#### *3.1. Test Cases*

After derailment, wheel–rail interaction disappears, and the train runs along the top surface of the track. In this case, various contact and impact conditions, such as gear box–rail, wheel–rail fastener, and wheel–sleeper, arise. These various conditions affect the post-derailment behavior. The post-derailment behavior is an important consideration for the design of derailment protection facilities. In this study, full-scale train bogie derailment tests were conducted. The test bogie and test methodologies are detailed in Section 2. The main test parameter of this study was the derailment velocity. Two different derailment velocities were considered. For Cases #1 and #2, the derailment speeds were 28.08 km/h and 55.05 km/h, respectively.

#### *3.2. Case #1 (Derailment Velocity of 28.08 km*/*h)*

For Case #1, the target speed range was 25–30 km/h. The measured derailment velocity was 28.08 km/h. Figure 16 shows a front view of the post-derailment behavior in Case #1. It can be seen that the bogie ran in the right lateral direction after derailment (0.410–1.810 s), and the left wheel continuously contacted the hump of the concrete track (1.810–3.570 s). Then, the bogie moved to the left direction without contact between the right rail and left wheel (3.570–5.630 s). The direction was changed due to the effect of the hump in this case. In this study, the Rheda 2000 type concrete track was used, and the hump in this track affected the post-derailment behavior.

The acceleration data were analyzed to calculate the velocity and traveling distance of the test bogie. The analysis results were compared with the results obtained in the image data from the high-speed cameras for cross-validation of the data. It is known that the image analysis offers a certain potential for the dynamic investigation [20]. The velocity and traveling distance of the test bogie were also calculated by target mark tracking from the image data, as shown in Figure 17. The results are shown in Figure 18 as a solid line. It can be seen that the velocity obtained from the high-speed camera showed initial fluctuation. However, the overall trend of the velocity was similar to that from the accelerometer. The velocity decreased almost linearly, as shown in Figure 18a. The slope of Figure 18a represents the deceleration, and it was approximately −0.97 m/s2. By integrating the acceleration data twice, the traveling distance could be calculated, as shown in Figure 18b. The traveling distances obtained from the accelerometer and high-speed cameras were similar to each other. The discrepancy may come from accumulated error during the integration or distortion of the image by the use of a wide-angle lens.

**Figure 16.** Post-derailment behavior of the bogie for Case #1, front view.

**Figure 17.** Example of image analysis by target mark tracking.

**Figure 18.** Data analysis results, Case #1: (**a**) velocity, (**b**) traveling distance.

The longitudinal, lateral, and vertical acceleration data of the test bogie were analyzed with the video images together. The first impact was observed at 0.410 s after derailment, as shown in Figure 19. It can be seen that there were sudden changes in acceleration at the time of the first impact, as shown in Figure 19a. The maximum acceleration was found in the vertical direction, and it was 5.07 g (49.78 m/s2). If the total mass of the bogie is included in the calculation of the impact force, it is equivalent to approximately 207 kN of impact force. The first impact occurred between the left wheel and hump, as shown in Figure 19b. Figure 19c shows the damage to the hump and rail fastener due to the first impact.

**Figure 19.** First impact, Case #1: (**a**) acceleration vs. time (0.3–0.6 s); (**b**) test bogie at first impact; (**c**) damage to the hump and rail fastener due to the first impact.

Similar to Figure 19, information about the second and third impacts is presented in Figure 20. The second and third impacts were observed at 0.750 s and 0.810 s, respectively. In the second impact, the right and left rear wheels impacted the fourth hump, as shown in Figure 20b,c. The maximum vertical acceleration was approximately 5.57 g—larger than that from the first impact. This may be attributed to the pitching moment of the bogie after the first impact. In the case of the third impact, the front right wheel contacted the sixth and seventh humps in succession. At the third impact, the peak of the vertical acceleration was markedly decreased.

**Figure 20.** Second and third impacts, Case #1: (**a**) acceleration vs. time (0.65–0.95 s); (**b**) test bogie at the second and third impacts; (**c**) damage to the hump and rail fastener due to the second and third impacts.

After the third impact, the test bogie ran the track without any further significant impact with track components. However, considerable lateral displacement was observed. For example, at 1.170 s, the front left wheel was located at almost the center of the track, as shown in Figure 21a. The lateral displacement continuously increased at 1.170–1.810 s, and the front left wheels contacted the humps in

the left rail at 1.810 s. Then, the test bogie returned to the center of the track at 3.750–5.630 s, as shown in Figure 21b.

(**b**)

**Figure 21.** Position of the test bogie: (**a**) 1.170–1.810 s, (**b**) 2.150–5.630 s.

In summary, for Case #1, where the derailment velocity was 28.08 km/h, the derailed test bogie ran the track with considerable lateral movement. The lateral movement was restrained by contact with the humps in the concrete track. Then, the test bogie returned to the center of the track. The whole trace of the test bogie is shown in Figure 22. The maximum vertical acceleration was approximately 5.57 g at the second impact with the humps. The maximum lateral acceleration was approximately 2.62 g at the first impact with the humps. From the results, it can be seen that the effect of the humps on the post-derailment behavior was significant when the Rheda 2000 track system was used. Thus, the type of track must be considered when evaluating the post-derailment behavior of a derailed train.

**Figure 22.** Trace of the left and right wheels for Case #1.

#### *3.3. Case #2 (Derailment Velocity of 55.05 km*/*h)*

In Case #1, the derailed wheel was not guided by the rail or containment wall. However, it was guided by the humps on the track. To identify the effects of the rail and containment wall, the derailment velocity was increased for Case #2. The target derailment velocity range was 50–55 km/h. The measured derailment speed was 55.05 km/h. Figure 23 presents the variation in the velocity and longitudinal traveling distance after derailment for Case #2. The data were obtained by the accelerometer as well as the high-speed cameras. However, the data from the high-speed cameras cannot be interpreted well after 2 s since the test bogie left the camera range due to the increased speed of the test bogie. The deceleration of the test bogie was approximately −1.48 m/s2. The traveling length was approximately 38 m at 3 s.

**Figure 23.** Data analysis results, Case #2: (**a**) velocity, (**b**) traveling distance.

Figure 24 shows the test bogie at derailment. It was found that the test bogie seemed to jump, bypassing the derailment device as shown in Figure 24a. Then, the front left wheel contacted with the third hump and rail fastener at 0.240 s due to the roll of the test bogie, as shown in Figure 24b. When the test bogie passed through the derailment device, the vertical acceleration in Case #2 was much greater than that in Case #1 before the first wheel contact with the track (0.242 s), as shown in Figure 24c. If higher-speed experiments are carried out, the derailment device should be improved because excessive jumping is expected.

**Figure 24.** Test bogie at derailment for Case #2: (**a**) 0.124 s; (**b**) 0.242 s (first impact); (**c**) vertical acceleration when derailed (0–0.20 s).

The main impact to the track was observed around from 0.4 s to 0.7 s, as shown in Figure 25. At 0.460 s, the front right wheel came into contact with the eighth rail fastener. Then, the front right wheel impacted the ninth hump and the rear right wheel immediately came into contact with the seventh rail fastener. The maximum vertical acceleration was found at this time, and it was approximately 3.53 g. The maximum acceleration in Case #2 was smaller than that in Case #1. This is because increased velocity is related with the longitudinal direction and vertical impact force is mainly a function of the velocity in the vertical direction.

**Figure 25.** Impact of the test bogie after derailment for Case #2: (**a**) acceleration vs. time (0.4–0.7 s), (**b**) 0.460 s; (**c**) 0.534 s.

After 0.7 s, the test bogie ran the track with considerable lateral movement, similar to Case #1. However, for Case #2, the wheel flange was in contact with the rail as shown in Figure 26b. Then, the test bogie impacted the containment wall at 1.910 s, as shown in Figure 26c. The lateral acceleration was increased due to the impact with the rail and containment wall, as shown in Figure 26a. The lateral acceleration values were approximately 1.17 g (11.51 m/s2) and 0.75 g (7.36 m/s2) for the rail and containment wall impact, respectively. Thus, it was expected that the impact force on the containment wall would be reduced approximately 36% by the first guide by the rail for Case #2. Finally, the test bogie returned to the track center direction after containment wall impact, as shown in Figure 27.

**Figure 26.** Impact of the test bogie after derailment for Case #2: (**a**) acceleration vs. time (1.7–2 s), (**b**) 1.780 s; (**c**) 1.910 s.

**Figure 27.** Trace of the left and right wheels for Case #2.

As mentioned in the introduction, the basic concept of DCP type III is that the derailed train is first guided by the rail and the containment wall, preventing excessive lateral movement [8,21,22]. The containment wall installed in this test is DCP type III, and it was seen that the derailed train was effectively guided by DCP type III in test Case #2.

#### **4. Conclusions**

This study presented full-scale train bogie derailment test methodologies and post-derailment behavior based on the test results. The derailment tests were conducted at the bogie level, and the test speeds were 28.08 km/h and 55.05 km/h. The major findings of this study are as follows:

**Author Contributions:** Investigation, Conceptualization, Methodology, Data Curation, Project Administration, Supervision, Writing (original draft, review & editing) from H.-U.B., J.M. and N.-H.L.; Experiment, Data processing from S.-J.L. and J.-C.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by a grant (19RTRP-B122273-04) from the Railway Technology Research Program, funded by the Ministry of Land, Infrastructure, and Transport of the Korean government.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Prediction of Friction Resistance for Slurry Pipe Jacking**

#### **Yichao Ye, Limin Peng, Yang Zhou, Weichao Yang \*, Chenghua Shi and Yuexiang Lin**

School of Civil Engineering, Central South University, Changsha 410075, China<EMAIL_ADDRESS>(Y.Y<EMAIL_ADDRESS>(L.P<EMAIL_ADDRESS>(Y.Z<EMAIL_ADDRESS>(C.S<EMAIL_ADDRESS>(Y.L.)

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +86-137-8723-2438

Received: 18 November 2019; Accepted: 23 December 2019; Published: 26 December 2019

### **Featured Application: The new approach established in this paper can provide accuracy prediction of friction resistance for slurry pipe jacking with various soil conditions, which lays a good foundation for better future design and less construction costs.**

**Abstract:** Friction resistance usually constitutes one of the two main components for the calculation of required jacking force. This paper provides a new approach to predict the friction resistance of slurry pipe jacking. First, the existing prediction equations and their establishment methods and essential hypotheses used were carefully summarized and compared, providing good foundations for the establishment of the new model. It was found that the friction resistance can be uniformly calculated by multiplying an effective friction coefficient and the normal force acting on the external surface of the pipe. This effective friction coefficient is introduced to reflect the effect of contact state of pipe-soil-slurry, highly affected by the effect of lubrication and the interaction of pipe-soil-slurry. The critical quantity of pipe-soil contact angle (or width) involved may be calculated by Persson's contact model. Then, the equation of normal force was rederived and determined, in which the vertical soil stress should be calculated by Terzaghi's silo model with parameters proposed by the UK Pipe Jacking Association. Different from the existing prediction models, this new approach has taken into full consideration the effect of lubrication, soil properties (such as internal friction angle, cohesion, and void ratio), and design parameters (such as buried depth, overcut, and pipe diameter). In addition, four field cases and a numerical simulation case with various soils and design parameters were carefully selected to check out the capability of the new model. There was greater satisfaction with the measured data as compared to the existing models and the numerical simulation approach, indicating that the new approach not only has higher accuracy but is also more flexible and has a wider applicability. Finally, the influence of buried depth, overcut, and pipe diameter on the friction resistance and lubrication efficiency were analyzed, and the results can be helpful for the future design.

**Keywords:** slurry pipe jacking; friction resistance; effective friction coefficient; pipe-soil-slurry interaction; lubrication efficiency

#### **1. Introduction**

In many parts of the world, the numerous constructions of municipal tunnels are creating unforeseen problems, such as blocking of roads, existing pipelines failure, and buildings subsidence. This has motivated attempts at the development of trenchless construction technology, such as pipe jacking, especially in metropolitan cities [1,2]. Pipe jacking is defined as a trenchless excavation technique, which employs hydraulic jacks to thrust specially made pipes through the ground behind a jacking machine, from a drive shaft to a reception shaft, as illustrated in Figure 1. It has many technical merits, such as a short time limit, high security, low environmental effect, and little traffic

disturbance [3–6]. Because of that, pipe jacking has been widely used in the construction of infrastructure for traffic and transportation systems in cities [7,8].

In pipe jacking, the jacking force is a critical factor that determines the thickness of pipe and reaction wall, selection or design of jacking machine and lubricant requirements [9]. The accuracy of prediction of jacking force is directly related to the structural safety and construction cost.

The main component of the jacking force is due to frictional resistance. Application of a lubricant such as bentonite slurry in pipe jacking (so-called 'slurry pipe jacking') is essential to reduce the friction resistance and, therefore, the jacking force [3,10–13]. However, the use of slurry makes it more complex to calculate or predict the friction resistance because of the change in contact conditions between the pipe and soil and lubricant slurry. The new contact state, which is due to the pipe-soil-slurry interaction, is affected by factors such as pipe diameter, soil properties, overcut [3,9], lubrication efficiency [9], pipeline misalignments [10,14–16], and stoppages [9,14,17–19]. The existing prediction models have not fully taken these factors into consideration, leading to an overestimation or underestimation of the friction resistance [9,20–23]. It is therefore obvious that a new prediction approach or model is imperatively needed to be established to solve the problem in slurry pipe jacking [24].

**Figure 1.** Schematic of slurry pipe jacking.

#### **2. Overview of the Existing Prediction Models of Friction Resistance**

Numerous models that calculate the friction resistance of pipe jacking have been proposed by authors from all over the world. The superposition principle is usually used, which holds that the comprehensive outcome of two or more linear factors of a system is equal to the accumulation of the effect of each factor. In pipe jacking, the linear factors to generate friction resistance are due to the weight of pipe (*fW*), soil pressure (*fs*), slurry pressure (*fm*), pipe-soil cohesion (*fsc*), and pipe-slurry cohesion (*fmc*). Their equations can be expressed as [5,6,9,10,14,19,21,22,25,26].

$$f\_W = \mu\_s \mathbf{W} \tag{1}$$

$$f\_s = \mu\_s N\tag{2}$$

$$f\_m = \mu\_m \text{N} \tag{3}$$

$$f\_{\text{sc}} = c\_s B\_s \tag{4}$$

$$f\_{mc} = c\_m B\_m \tag{5}$$

where μ*s* and μ*m* are the kinematic friction coefficient of pipe-soil and pipe-slurry, respectively; *W* is the weight of pipe per unit length, kN/m; *N* is the total normal force acting on the pipe, kN/m; *cs* and *cm* are the pipe-soil cohesion resistance and pipe-slurry cohesion resistance, respectively, kPa; *Bs* and *Bm* are the pipe-soil contact width and pipe-slurry contact width, respectively, m.

Some hypotheses have been made to establish the prediction models, by which one or some of the items listed above should be considered. Typical hypotheses are:

**Hypothesis 1.** *The excavated tunnel is self-stable, the pipeline simply slides along the bottom of the tunnel due to its own weight (see Figure 2a) [9,10,14].*

**Hypothesis 2.** *The angular space due to overcut is completely filled with lubricant slurry, and the excavated tunnel is stable under the slurry pressure (see Figure 2b) [22,25].*

**Hypothesis 3.** *The excavated tunnel is unstable, the surrounding soil collapses and is in full contact with the whole area of the jacking pipes (see Figure 2c) [5,6,9,21,26].*

**Hypothesis 4.** *The excavated tunnel is stable under the pressure of slurry, and part of the pipe comes in contact with the surrounding soil (see in Figure 2d) [3].*

**Figure 2.** The models to calculate friction resistance according to different hypotheses: (**a**) Hypothesis 1; (**b**) Hypothesis 2; (**c**) Hypothesis 3; (**d**) Hypothesis 4.

According to Hypothesis 1, there are two kinds of prediction models, in which the item of *fw* has to be taken into consideration. The first one assumes that the friction resistance is only due to the weight of pipe [9,21]. It is given as

$$f = \mu\_s \mathcal{W} \tag{6}$$

$$
\mu\_s = \tan(q/2) \tag{7}
$$

The second one also takes the pipe-soil cohesion resistance (item of *fsc*) into account [10,14,21], which is given by

$$f = \mu\_s \mathcal{W} + \mathfrak{c}\_s B\_s \tag{8}$$

where the pipe-soil contact width *Bs* is calculated by Hertzian contact model, as [10,14,21]

$$B\_s = 1.6(\mathrm{Pk}\_d \mathbb{C}\_\varepsilon)^{1/2} \tag{9}$$

$$k\_d = \frac{D\_c D\_p}{D\_c - D\_p} \quad , \quad \mathbb{C}\_c = \frac{1 - v\_p^2}{E\_p} + \frac{1 - v\_s^2}{E\_s} \tag{10}$$

where *Dc* and *Dp* are the internal diameter of cavity and external diameter of pipe, respectively, m; ν*p* and ν*s* are the Poisson's ratio for pipe and soil material, respectively; *Ep* and *Es* are the elasticity modulus of the pipe and soil, respectively, kPa; and *P* is the effective external force acting on the center of the pipe, kN/m, usually it is considered as equal to the weight of pipe per unit length *W*.

According to Hypothesis 2, the friction resistance is related only to the properties of lubricant slurry, and the only model found was one that takes both the items of *fm* and *fmc* into account [22,25].

$$f = \mu\_m \mathcal{N} + \mathfrak{c}\_m B\_m \tag{11}$$

$$N = \pi D\_p P\_m \tag{12}$$

$$B\_m = \pi D\_p \tag{13}$$

where *Pm* is the mud slurry pressure.

It is noted that most of the studies completed to date have focused exclusively on the prediction models established by Hypothesis 3. This may be attributed to the assumption of full contact of the pipe and soil that leads to a large prediction value of friction resistance, or in other words, Hypothesis 3 is conservative. Because of that, this kind of model is widely accepted by authorities and standards from all over the world, such as Japan Sewage Association (JSA) [22], UK Pipe Jacking Association (PJA) [27], Chinese Trenchless Technology Association (CTTA) [28], and Germany Standard (AVT A-161) [29]. This kind of model can be summarized and divided into the following four categories.

$$f = \mu\_{\text{s}} \mathcal{N} \tag{14}$$

$$f = \mu\_s \mathcal{N} + \mu\_s \mathcal{W} \tag{15}$$

$$f = \mu\_s \mathcal{N} + \mu\_s \mathcal{W} + \mathcal{c}\_s B\_s \tag{16}$$

$$f = \beta \mu\_s \mathcal{N} + \mu\_s \mathcal{W} + c\_s B\_s \tag{17}$$

The fourth kind of model (Equation (17)) introduces an empirical constant β (smaller than 1) on the basis of the third model (Equation (16)) to reflect the effect of lubrication.

For the models summarized above, the item of *fs* (=μ*sN*) has to be taken into consideration. Thus, the key problem for this kind of model is exactly focused on the calculation of soil pressure *N.* CTTA suggests *N* to be calculated by Rankine's formula, which gives that [28]

$$N = 2(1+K)D\_p \sigma\_v \tag{18}$$

$$
\sigma\_{\mathcal{V}} = \gamma h \tag{19}
$$

$$K = \tan^2(\pi/4 - \phi/2)\tag{20}$$

where σ*v* is the vertical soil stress; *K* is the coefficient of soil pressure above the pipe; γ is the unit weight of soil; *h* is the overburden depth of the pipeline; and ϕ is the internal friction angle of soil.

However, JSA suggests using Terzaghi's silo model, the expression of *N* is then expressed as [22]

$$N = \pi D\_p \sigma\_v \tag{21}$$

$$\sigma\_{\upsilon} = \frac{b\gamma - 2c\_s}{2K \tan(\delta)} \left( 1 - e^{-2K \tan(\delta) \hbar/b} \right) \tag{22}$$

where *cs* is the cohesion of soil; δ is the friction angle between the pipe and soil; *b* is the influencing silo width of soil above the pipe, and the other symbols have the same meanings as before.

Although a lot of prediction equations of friction resistance have been proposed and even some of them have been applied to engineering practice, it is obvious that their hypotheses are quite different, and even for the same hypothesis, models and parameters can be different too. Thus, it is bound to make the prediction friction resistances vary greatly. Furthermore, apart from Equations (11) and (17), the other models completely ignore the effect of lubrication, which is very important to slurry pipe jacking.

In the design philosophy of slurry pipe jacking, the angular space due to overcut is expected to be completely filled with lubricant slurry, to reduce the friction resistance with maximum efficiency, creating a 'filter cake' layer around the cavity and is then pressurized to the support pressure required for the soil (see Figure 2b) [3,10]. In this case, the friction resistance should be only related to the slurry pressure and the friction coefficient between slurry and the pipe. From this point of view, the expression Equation (11) seems a convincing explanation here.

However, the more general case is that the excavated tunnel is stable under pressure of slurry and part of the pipe inevitably comes in contact with the soil (see Figure 2d) [3]. The reasons for the occurrence of pipe-soil contact can be complex, such as insufficient design and control of grouting amount of slurry, the pipeline deviates from the intended line and level, irregular deformation of the surrounding soil, and interpenetration between the soil and slurry. Thereby, the state of contact can change from 'pipe-slurry' into 'pipe-soil-slurry' (Figure 2b,d). In such a case, a simple form of Equation (17) seems more suitable to reflect the effect of lubricant slurry. However, it is logical that the greater the contact width between the pipe and the soil, the smaller the effect of lubrication will be, and, therefore, the greater the friction resistance will be. Thus, the value of β should be highly affected by the pipe-soil contact width. For different soils, grouting amount of slurry and design parameters (such as buried depth and overcut) is bound to lead to completely different contact widths of pipe-soil. Thus, β should be in a large range, and it would be rather difficult to pick out a value to use in application.

In fact, a successful prediction model of friction resistance should not only consider the effect of lubrication but also needs to be able to reflect the effect of pipe-soil-slurry interaction. It is based on this understanding that the following model comes into being.

#### **3. Calculation of Friction Resistance for Slurry Pipe Jacking**

The general contact state of pipe-soil-slurry due to the interaction of the pipe, surrounding soil, and lubricant slurry is shown in Figure 3. In the picture, the position of pipe-soil contact is arbitrary, with a contact width of *Bs* and the corresponding contact angle of 2ε.

**Figure 3.** The general contact state of pipe-soil-slurry.

From Figure 3, it is obvious that to calculate friction resistance *f*, both of the items of *fs* and *fm* have to be taken into account, which can be expressed as

$$f = \mu \mathbf{N} = f\_{\mathbf{s}} + f\_{m} = \mu\_{\mathbf{s}} \mathbf{N}\_{\mathbf{s}} + \mu\_{m} \mathbf{N}\_{m} \tag{23}$$

where μ is an effective friction coefficient introduced to reflect the effect of lubrication and the influence of pipe-soil contact width. It is generally accepted that μ*s* = tan(ϕ/2) for the coefficient of kinematic friction between soil and the pipe [21,22]; μ*m* for the coefficient of kinematic friction between mud slurry and the pipe can be taken as 0.01, according to the test result reported by Guo [30]. *Ns* and *Nm* are the total normal force of pipe-soil and pipe-slurry in contact, respectively.

To calculate *Ns* and *Nm* precisely, the location of pipe-soil contact and the magnitude of contact angle (or contact width) and the contact force have to be determined. For various reasons leading to the occurrence of pipe-soil contact, it seems impossible to calculate these quantities in a target section of the pipeline. However, if taking the whole pipeline into consideration, and assuming that the pipe-soil contact can occur at any position of a section of the pipeline with a same probability and the contact force is approximately equal to the soil pressure in the contact area, this problem can be greatly simplified. In this case, we have the following equations:

$$P = N\_s = \frac{B\_s}{\mathbb{C}} N = \frac{\varepsilon}{\pi} N \tag{24}$$

$$N\_m = \frac{B\_m}{C} N\_c \tag{25}$$

where *C* (=π*Dp*) is the external circumference of pipe and ε is the semi-angle of contact (Figure 3). By substituting Equations (24) and (25) in Equation (23), after some algebra, giving that

$$\begin{cases} \mu = \mu\_{\text{s}} \lambda\_{\text{s}} + \mu\_{\text{m}} \lambda\_{\text{m}}\\ \lambda\_{\text{s}} = \frac{B\_{\text{s}}}{C} = \frac{\varepsilon}{\pi} \end{cases}, \quad \lambda\_{\text{m}} = \frac{B\_{\text{m}}}{C} \tag{26}$$

$$B\_m = C - \frac{B\_s}{1+c} \tag{27}$$

where *e* is the void ratio of soil.

By substituting Equation (27) in Equation (26), the expression of μ can be further rewritten as

$$
\mu = \mu\_s \frac{\varepsilon}{\pi} + \mu\_m \left( 1 - \frac{\varepsilon}{\pi(1+\varepsilon)} \right) \tag{28}
$$

According to Equation (28), the calculation of ε is essential to calculate μ*.* Hertzian model provides a simple way for the calculation of the width of contact (or contact angle) as we have mentioned before (see Equations (8) and (10)); however, the Hertzian contact problem is approached only when the applied force is small, or the large radial clearance is large, and the limited angle of contact is smaller than about 30◦ [31]. Due to the technical limitations, most of the pipe jacking projects encounter clay or sandy soils and with a small overcut, it is therefore important that the applicability of the Hertzian contact model should be extremely limited here. Actually, the Hertzian contact model is just a special case of the Persson's contact model with a small contact width (or angle) [31]. If a large contact angle (larger than 30◦) occurs, the more general contact model proposed by Persson should be taken as the first choice. The following singular integro-differential governing equation of contact angle is derived by Persson, as [31,32]

$$B = 4(1 - \beta) - 2(1 - \alpha) \int\_{-\xi}^{+\xi} q(t) \frac{dt}{1 + t^2} - \frac{\pi}{2}(1 + \alpha) \frac{E\_p \Delta R}{(1 - v\_p^2)P} \tag{29}$$

or

$$\frac{\pi(1+a)E\_p\Delta R}{(1-v\_p^2)P} = 4(1-\beta) - 2(1-a)\int\_{-\xi}^{+\xi} q(t)\frac{dt}{1+t^2} - B \tag{30}$$

The involved auxiliary variables are defined as [31]

$$\begin{array}{ll} \Delta \mathsf{R} = \frac{D\_{\mathsf{r}} - D\_{\mathsf{p}}}{2} \\ \Delta \mathsf{R} = \frac{1 - \eta}{1 + \eta} \quad , \quad \mathsf{\beta} = \frac{\lambda}{2(1 + \eta)} \\ \eta = \frac{E\_{\mathsf{p}}}{E\_{\mathsf{s}}} \frac{1 - v\_{\mathsf{s}}^{2}}{1 - v\_{\mathsf{p}}^{2}} \quad , \quad \lambda = \frac{1 - 2v\_{\mathsf{p}}}{1 - v\_{\mathsf{p}}} - \eta \frac{1 - 2v\_{\mathsf{s}}}{1 - v\_{\mathsf{s}}} \end{array} \tag{31}$$

After some approximate treatments, the key terms of Equation (30) have been solved by Michele [32], as follows:

$$\int\_{-\xi}^{+\xi} q(t) \frac{dt}{1+t^2} = \frac{1}{2\pi} \frac{I\_b}{\xi^2(\xi^2+1)} + \frac{\xi^2}{\xi^2+1} \tag{32}$$

$$I\_b = \pi \log(\xi^2 + 1)\tag{33}$$

$$B = \frac{2\xi^4 + 2\xi^2 - 1}{\xi^2(\xi^2 + 1)}\tag{34}$$

By substituting these equations into Equation (29), Michele obtained an approximate form of general contact angle relation.

$$\frac{\pi(a+1)E\_p\Delta\mathcal{R}}{(1-v\_p^2)P} = \frac{(a-1)[\ln(\xi^2+1)+2\xi^4]+2}{(\xi^2+1)\xi^2} - 4\beta\tag{35}$$

As compared with Equation (9), Equation (35) is a far more complex nonlinear equation. It can be further simplified with respect to that the elastic modulus of soil *Es* is generally much smaller than that of pipe *Ep* (the difference between the two can be three orders of magnitude). Thus, from Equation (30), the magnitude of auxiliary variable η should be very large, and, therefore, the following approximate relations can be obtained: π(α+1)*Ep*

$$\begin{cases} \frac{\pi (a+1)E\_p}{(1-v\_p^2)} \approx \frac{2\pi E\_s}{(1-v\_s^2)'}\\ \alpha \approx -1\_\prime\\ \beta \approx \frac{1-2v\_s}{2(1-v\_b)} \end{cases} \tag{36}$$

Using Equation (36), Equation (35) can be then simplified as

$$\frac{\pi E\_{\text{s}} \Delta R}{(1 - \upsilon\_{\text{s}}^2)P} + \frac{1 - 2\upsilon\_{\text{s}}}{1 - \upsilon\_{\text{s}}} = \frac{1 - [\ln(\xi^2 + 1) + 2\xi^4]}{(\xi^2 + 1)\xi^2} \tag{37}$$

From Equation (37), it is essential to calculate *P*, which requires one to calculate the total normal force *N*. It can be gained by integrating the normal stress σ*n* on an element of the pipe surface and is determined on the basis of vertical and horizontal soil stresses.

$$
\sigma\_n = \sigma\_\upsilon \sin \theta + \sigma\_d \cos \theta \tag{38}
$$

$$N = 4\int\_0^{\pi/2} \sigma\_n \frac{D\_p}{2} d\theta \tag{39}$$

where θ is defined as the angle between the corresponding radius line and the horizontal line at each point of the pipe (Figure 4).

By substituting Equation (38) in Equation (39), it is easy to obtain the equation of *N*, which has the same form as Equation (18).

To calculate Equation (18), the vertical soil stress σ*v* has to be first determined. It is noted that at the present time, by far the most commonly used model for soil pressure calculation is Terzaghi's silo model (Equation (22)) [5,6,9,19,21,26]. According to Equation (22), the calculation of the vertical soil stress requires some physical parameters that may be determined with some accuracy, such as the height of cover *h*, the cohesion *cs*, and the unit weight of soil γ, but also some empirical parameters, such as *b*, δ, and *K*. The definition of these empirical parameters varies from one author to another. Here, typical approaches of Terzaghi, Germany Standard ATV-A 161 E-90 [29], Chinese Standard GB 50332-2002 [33], UK Standard BS EN 1594-09 [34], US Standard ASTM F 1962-11 [35], UK PJA [27], Japan JMTA [36], and Japan JSA [22] would be discussed and compared.

For the calculation of silo width *b*, three kinds of boundary planes of wedge failure assumed by different authors and the corresponding equations are clearly illustrated in Figure 5. The width of the boundary plane is related to the 'vault' effect of soil. Generally, a smaller *b* means a lower 'vault' effect of soil, leading to a larger vertical soil stress.

**Figure 4.** The earth pressure and the normal stress acting on the pipe.

**Figure 5.** Boundary planes of wedge failure assumed by different authors. JSA: Japan Sewage Association; PJA: UK Pipe Jacking Association.

For the determination of δ, most of the guidelines, such as PJA, JSA, JMTA, BS EN 1594-90, and GB 50332-2002 assume shear planes as perfectly rough and take an angle of friction in the shear planes δ equal to the soil internal friction angle ϕ*.* However, ATV-A 161 E-90 and ASTM F 1962-11 make a more cautious assumption and only takes into account half the internal friction angle ϕ/2.

For the lateral pressure coefficient *K* above the tunnel, Terzaghi assumes *K* coefficient is equal to 1, which corresponds to the range of values encountered in clayey soils. PJA, ASTM F 1962-11, and GB 50332-2002 suggest *K* = *Ka* (calculated by Rankine's formula of active soil pressure coefficient), while BS EN 1594 and ATV-A 161 assume *K* = *K*0 (calculated by Rankine's formula of soil pressure coefficient at rest). Moreover, according to ATV-A 161, this *K* coefficient is equal to 0.5, which corresponds to an internal soil friction angle of 30◦, a typical value for sandy soils.

Parameters of *b*, δ, and *K* chosen by the different authors have been summarized in Table 1.

**Table 1.** Definition of the empirical parameters in Terzaghi's silo model by different authors [37].

Note: α = π/4 − ϕ/2; β = π/4 + α/2.

It is noted that none of the approaches use the same parameters. Consequently, the vertical soil stress calculated by these approaches would be quite different. Thus, it is not convincing to pick out an approach to use without checking the field data. This work will be carried out in the next section.

Thus far, all the equations needed to calculate friction resistance have been determined. If the parameters needed for the prediction equations are quantified, by using Equations (18) and (22) the total normal force *N* can be determined, then together with Equations (23), (24), (27), (31), and (37), the contact angle 2ε, the effective friction coefficient μ, and the friction force *f* now can be uniquely identified. The flow chart is shown in Figure 6.

**Figure 6.** Flow chart of friction resistance prediction.

Apparently, the effective friction coefficient here is not just related to the interfriction angle of soil ϕ but also the state of pipe-soil contact and the effect of lubrication.

#### **4. The Verification of the E**ff**ectiveness of the Proposed New Approach**

*4.1. Comparison between the Predicted Friction Resistances and the Field Data*

Before comparison analysis, the existing prediction equations were numbered from M 1 to M 13 (Table 2). For M 7~M 13, the Equations (7), (14), (18), and (22) should be used and the parameters (*b*, δ, and *K*) of Terzaghi's silo model (Equation (22)) in the calculation of earth pressure should be determined by their approaches that are shown in Table 1.

**Table 2.** The numbering of the existing friction resistance prediction equations.

As has been mentioned before, the new approach of this paper introduces an effective friction coefficient μ (Equation (28)) to replace the original pipe-soil friction coefficient μ*s* (Equation (7)). The calculation of μ (Equation (28)) should follow the flow chart shown in Figure 2. Equations (18) and (22) are also used to calculate the normal force *N* and vertical earth pressure σ*v*, respectively. The new prediction equations were numbered from M 7 \* to M 13 \* (Table 3).

**Table 3.** The numbering of the new friction resistance prediction equations.

The superscript "\*" is used to distinguish from the number of models in Table 2.

Some parameters required by the predicted equations (such as μ*s*, *Pm*, and *W*) might not be given in the literature. In principle, during the calculation, the parameters given in the field case should be used, and the missing parameters can be evaluated by the following rules:

6. The average values of soil parameters (such as γ*s*, *cs*, ϕ, and *e*) can be obtained from the Geological Engineering Handbook.

Parameters in each of the four cases were finally determined, which are given in Table 4 [2,9,22]. These cases encountered some representative soils, such as silt, clay, sand, and gravels. Furthermore, they have different overburden depths of 2.72~8.5 m, overcut of 0~20 mm, and pipe diameters of 0.66~4.06 m. All of these characteristics provide good foundations for identifying the capability of the prediction models.

**Table 4.** Parameters that are needed to calculate the prediction equations in each case.

In Table 5, for each of the drives, measured frictional force values are presented and compared to values calculated by the approaches of the existing models. One can see that, for most cases, the prediction results of the models (M 1 and M 2) established based on Hypothesis 1 are generally too small. This is because it is correct only when the overcut is stable and the pipeline slides on its base inside the annular gap remaining open. The same problem is encountered in M3 (established based on Hypothesis 2), which is probably more due to ignoring the occurrence of pipe-soil contact decreasing the magnitude of the effective friction coefficient. Obviously, for slurry pipe jacking, the prediction friction based on these two assumptions may be insufficient and unsafe.

**Table 5.** Comparison of friction resistances calculated by the existing models and the measured data.

Note: *f* mea. is the measured friction resistance; *f* cal. is the calculated friction resistance.

The calculated results of M4 are much larger than the measured data, which presumably result from ignoring the 'vault effect' of soil, so that the soil pressure has been overestimated. Apart from M4, the other models established based on Hypothesis 3 (M5~13) have shown some applicability in case 4, in which the overcut of C4 is equal to zero, which is exactly in line with the assumption of the full contact of pipe-soil of Hypothesis 3. Except for C4, the prediction results of other cases are much larger, and the amplitude may be up to 30 times the measured data. Therefore, for slurry pipe jacking, Hypothesis 3 is generally over-conservative. Although it can ensure the structural safety of the design, it may also cause a much higher construction cost.

From what we have analyzed above, the existing models have good prediction results only when the field case is consistent with the basic hypothesis of each model. However, it is also these hypotheses that determine their limited applicability.

Table 6 presents the friction resistances calculated by the approach of this paper, by using the parameters of Terzaghi silo model chosen from different authors (Table 1). Better agreements with the field data in each of these cases were obtained, which indicate that compared with the existing prediction models, the approach of this paper that considers the effect of lubrication and interaction of pipe-soil-slurry not only has higher accuracy but is also more flexible and has wider applicability.

**Table 6.** Comparison of the friction resistances calculated by the new approach and the measured data.

The superscript "\*" is used to distinguish from the number of models in Table 2.

Moreover, the approaches of PJA (M9), ATV A 161 E-90 (M11), ASTM F 1962-11 (M12), and GB 50332-2002 (M13) provide higher prediction accuracy than the other ones (Table 6). However, it is noted that ATV A 161 E-90, ASTM F 1962-11, and GB 50332-2002 suggest the effect of soil cohesion *c* should be neglected (Table 1), but it is because the cohesion of soil has been considered in the calculation that the good prediction results can be obtained, as shown in Table 6. From this point of view, UK PJA (consider the cohesion of soil) provides the best choice of parameters of Terzaghi' silo model in the calculation of earth pressure. Therefore, the UK PJA approach is supposed to be used in the next analyses.

#### *4.2. Comparison between the Predicted Friction Resistances and the Numerical Simulation Result*

Some excellent work of numerical simulation has been conducted for the estimation of the jacking force (or friction resistance) of pipe jacking. For example, Ji et al. [5,40] and Barla and Camusso [41] presented a novel discrete 2D numerical model to evaluate the normal force acting on the pipe. Then, the friction resistance is determined by multiplying the interface friction coefficient by the normal force. Three different drives in a pipe jacking projects were analyzed, where the calculated pattern of jacking force was compared with the measured data, which demonstrated the effectiveness of the proposed approach. Yen and Shou [25] used a model coupling finite element method and a displacement control method to estimate the required jacking force in pipe jacking. The displacement control option in the numerical analysis software ABAQUS (Abaqus Inc., (Palo Alto, CA, USA) 2012) was used to designate the displacement at the end cross section of the pipe in the launch shaft. Accounting for the contact property and the contact range between the pipes and the soil during the jacking process, the stresses exerted on the pipes were used to back-calculate the jacking forces [25].

The numerical simulation strategy of Yen and Shou [25] is quite consistent with the analytical approach of this paper. Therefore, the numerical simulation approach of Yen and Shou will be described in the next section and its estimated result of friction resistance will be discussed and compared with that calculated by the analytical approach of this paper.

One of the analyses focused on a case of slurry pipe jacking in the Taichung Science Park, the basic parameters of this case were summarized in Table 7. In the numerical simulation, the lateral boundaries were fixed by roller, and hinges were used to constrain the bottom boundary. The three dimensional hexahedron element (C3D8R) with eight nodes was used in the simulation. To verify the jacking force obtained from the displacement control simulation, overcut and lubrication were included in the model by setting the contact range (1/3, 1/2, and 1 pipe-soil contact) and the frictional coefficient (Figure 7a–c). The results suggested the 1/3 contact case can estimate the jacking force with a better accuracy towards the middle and the final stage of the pipe jacking process (Figure 8).

**Table 7.** The details of the slurry pipe jacking in Taichung Science Park.

**Figure 7.** The setting of pipe-soil contact range and frictional coefficient in the numerical simulation. (**a**) 1/3 of the pipe-soil contact; (**b**) 1/2 of the pipe-soil contact; (**c**) 1 of the pipe-soil contact.

In Figure 8, the slope of the linear regression for the scattered points represents frictional resistance, and the intercept represents face resistance. For the calculated result of this paper, the measured face resistance is used, and the frictional resistance (*f* cal. = 57.7 kN/m) is calculated by the new approach (Section 3) using parameters of UK PJA in Terzaghi's silo model (Table 1). For the result of numerical simulation, the decrease of jacking force in the initial stage of drives (before 62.5 m) is explained by authors that the weight of the pipe jacking machine (the material of which is steel, featuring a larger unit weight) directly pressing on the soil, causing a larger face resistance. After the drives of 62.5 m, the influence exerted by the machine weight decreased gradually and the increase of jacking force is caused by the accumulation of friction resistance [25]. In other words, the slope of regression from drives of 62.5 to 100 m represents the friction resistance estimated by numerical simulation (*f* num. = 64.4 kN/m).

It is obvious that both the predicted friction resistances of the analytical equations (*f* cal. = 57.7 kN/m) and the numerical simulation (*f* num. = 64.4 kN/m) are acceptable as compared to the measured data (*f* mea. = 41.5 kN/m). However, better accuracy (especially towards the first 80 m of drives) is obtained by the approach of this paper.

**Figure 8.** The comparison of jacking forces from monitoring, the numerical simulation, and the predicted equations of this paper.

The predicted accuracy of friction resistance by the numerical simulation approach is highly affected by the setting of pipe-soil contact angle. The pipe-soil contact angle of 71 degrees is that calculated by the analytical equations for this field case. Thus, it seems that the pipe-soil contact angle of 120 degrees (1/3 contact) in the simulation is set too large, resulting in the calculated result slightly larger than the measured data. From this point of view, the numerical simulation approach can accurately predict the friction resistance, but the contact angle of pipe-soil in the simulation needs to be set reasonably to obtain good prediction results. Conversely, in the approach of this paper, the pipe-soil contact angle is theoretically calculated with respect to soil property, overcut, pipe diameter, etc., while human factors or empirical factors can be eliminated as far as possible.

#### **5. Influence of Design Factors on Lubrication E**ffi**ciency and Friction Resistance**

For a better design of slurry pipe jacking in the future, it is meaningful to study the influence of design parameters (such as buried depth *h*, pipe diameter *Dp*, and the overcut Δ*R*) on lubrication efficiency and friction resistance. To achieve this objective, a set of reference parameters is used, and then, by changing a target parameter according to the design rules, the effect of that parameter on friction resistance and lubrication efficiency can be obtained. The designed cases were shown in Table 8.

#### *5.1. Influence of Design Factors on Friction Resistance*

The influences of design factors (*h*, *Dp*, and Δ*R*) on the critical quantities of effective friction coefficient μ, normal force acting on the pipe *N*, and friction resistance *f* are respectively shown in Figure 9a–c, Figure 9d–f, Figure 9g–i.

**Figure 9.** The influence of factors of *h*, *Dp*, and Δ*R* on the quantities of μ (**a**–**c**), *N* (**d**–**f**)*,* and *f* (**g**–**i**).

It is evident that increasing buried depth and pipe diameter led to a double action for the increasing of friction resistance. Firstly, this increase then increased the possibility of contact between the pipe and soil, and, therefore, increase the effective friction coefficient μ on the interface. Secondly, the increasing of buried depth increases the vertical soil stress, while the increasing of pipe diameter increases the contact area, both effects of them increase normal force *N* acting on the pipes. The main difference between the two is that buried depth causes both of the effective friction coefficient and normal force to slightly increase and then gradually stabilize, while they approximately linearly increase with pipe diameter. Especially for the normal force induced by pipe diameter, which is strongly increased from 29.83 to 991.34 kN, this leads to a notable increase of friction resistance from 0.34 to 48.16 kN/m. Thus, the additional friction is strongly affected by the pipe diameter but appears not to be greatly affected by the buried depth.

Different from the buried depth and pipe diameter, the overcut with small values has no effect on the normal force (Figure 9h). However, it has a strongly negative effect on the effective friction coefficient on the interface (Figure 9g). In fact, it does determine the volume of injected lubricant slurry, which has a significant influence on the occurrence of the pipe-soil contact, and, therefore, determine the magnitude of the effective friction coefficient. Thus, it has a strongly negative effect on the friction resistance.

#### *5.2. Influence of Design Factors on Lubrication E*ffi*ciency*

Except for friction resistance, engineers are also concerned about the lubrication efficiency [9,42,43]. According to Equation (27), the magnitude of μ is between the pipe-slurry friction coefficient μ*m* and the pipe-soil friction coefficient μ*s*. If there is no contact between the pipe and the soil, the angular space due to overcut is completely filled with lubricant slurry, the effective friction coefficient is equal to μ*m*; and if the soil is in full contact with the external surface of the pipe, the effective friction coefficient mainly depends on the pipe-soil nature and approximately equals to μ*s*. Thus, the lubrication efficiency can be defined as

$$
\chi = \left(1 - \frac{\mu - \mu\_m}{\mu\_s}\right) \times 100\% \tag{40}
$$

By substituting Equation (28) in Equation (31), and considering that μ*m* is far smaller than μ*s*, the χ can approximately be expressed by Equation (41).

$$\chi = \left(1 - \frac{\varepsilon}{\pi}\right) \times 100\% \tag{41}$$

If pipe-soil contact angle ε = 0, χ = 100% for maximum lubrication efficiency, and if ε = π, χ = 0% for minimum lubrication efficiency. It is noted that χ = 0% is not going to happen. According to the Passon model (Equation (37)), the result calculated by the left terms of Equation (37) should be not less than zero, while the right term of Equation (37) is a monotonically decreasing function of pipe-soil contact angle. In other words, when the right term of Equation (37) is equal to zero, the pipe-soil contact angle reaches its maximum value, which is solved by ε = ε*max* = 72◦ = 0.4π, corresponding to the minimum lubrication efficiency of 60% (i.e., χ is theoretically between 60% and 100% on the basis of Passon contact model). Although it is not theoretically correct as compared to that counted by Pellet and Kastner as between 45% and 90% [9] and tested by Zhou as between 47.8% and 78.6% [5], it seems practical to estimate the efficiency of lubrication by the approach of this paper.