Numerical Simulation of the Daikai Station Subway Structure Collapse due to Sudden Uplift during Earthquake

In this paper, a numerical simulation of the Daikai Station subway structure influenced by a single upward displacement wave of short duration was performed with a focus on the main section of the station, which suffered severe damage during the Hyogoken-Nanbu earthquake. Fig. 4 shows the numerical analysis model, and Fig. 4(a) shows the discretization of the symmetrical cross section of the whole structure, including the surrounding soils. The soil layers were discretized in the region from the base stratum to the ground surface and laterally more than twice the width of the main section of the station on both sides. The surrounding soils consist of seven layers, among which diluvial deposits (Pleistocene clay, sand, and gravel) have been overlain by alluvial deposits (Holocene clay and sand) (Iida et al. 1996). The decomposed granite soil was used as backfill material for the sidewalls.

The artificial finite boundaries that were previously defined must be set up for the stress waves to ensure that they are not reflected similarly to the infinite domain. In LS-DYNA (Livermore Software Technology Corporation 2016), due to the application of the normal and shear stress waves in Eqs. (1) and (2) to the boundaries, the nonreflecting boundaries for the stress waves were accomplishedwhere $\rho$, $C_d$, and $C_s$ = material density, dilatational wave speed, and shear wave speed, respectively, of the transmitting soil layer; the magnitude of these stresses is proportional to the particle velocity of the normal, $V_{\mathrm{normal}}$, and tangential, $V_{\tan }$, respectively.

A contact algorithm with the ability to consider sliding without friction and separating actions was employed for the contacted surfaces of the surrounding soil to all of the outer surfaces of the station structure. Fig. 4(b) shows the 3D discretization of the main section of the station without platforms for the numerical simulation. Because the central columns in the station were arranged at a center-to-center spacing of 3.5 m, a part of the station supported by one column was taken into consideration in the longitudinal direction. Therefore, the length of the station model in the longitudinal direction was 1.75 m considering symmetry. To accurately investigate the dynamic nonlinear behavior of the central column, the concrete for the column and the corresponding upper and lower girders was discretized using cubic elements of approximately 25 mm in length. The rebar pieces were separately discretized by using beam elements, which were perfectly bonded with the concrete elements based on the actual structure of the station before collapse. The rebar and concrete elements were bonded by coupling.

Regarding the boundary conditions, the nodal points on the symmetrical surface and end surface in the longitudinal direction were restrained in the normal directions. The damping factor was ignored because the maximum duration of the forced displacement was only 50 ms and the transient dynamic response characteristics of the station developed during the stress phase due to the propagation of reflection waves. Self-weight was also considered in the numerical analyses.

The applicability of calculation method employed in this study to the transient response analysis of the soil–structure system, which was subjected to a sudden uplift loading, was confirmed by conducting a preliminary elastic analysis and by checking the propagation state of the stress waves and dynamic behavior of the displacement and stress at the free surface and/or nonreflecting boundary of the soil layers. Regarding the failure behaviors of concrete members, by conducting drop-weight impact loading tests of RC beams and numerical simulation of these tests with LS-DYNA, a comparison between the experimental results and numerical results confirmed that the time history of the displacement and crack patterns of the RC beams can be adequately simulated by setting the longitudinal element size of concrete to a length range of 25–35 mm (Kishi et al. 2009, 2011; Kishi and Bhatti 2010). Therefore, the dynamic behavior of the central column may be accurately analyzed in this study.

Fig. 5 shows the stress-strain relationship for each material used in this numerical analysis: concrete, steel rebar, and overburden and backfilled soils. However, the surrounding soil layers, excluding the remolded soils described previously, were assumed to be elastic. The shear deformation in these layers may be very small due to the lack of lateral loading. The material properties were based on geological survey results obtained by boring and borehole observations, as listed in Table 1 (Iida et al. 1996; Sato Kogyo Co. Ltd. 1997). In this table, the layer numbers refer to those in Fig. 4(a).

A variety of constitutive models for concrete are provided in LS-DYNA (Livermore Software Technology Corporation 2016). In this study, the Karagozian and Case concrete model (KCC model or MAT072R3 in LS-DYNA), as shown in Fig. 5(a), was employed (Malvar et al. 1997) to better regenerate the axial compression failure of concrete. This model can effectively capture key concrete behaviors, including postpeak softening and confinement effects (Wu et al. 2012). Additionally, the model can automatically generate parameters by inputting the compressive strength $f_c^{\prime }$, Poisson’s ratio ${\nu }_c$, density of concrete ${\rho }_c$, and the element size (Malvar et al. 2000). These parameters are applied to determine the relationship between the pressure and the volumetric strain. Three kinds of shear strength surfaces, which are yield, maximum, and residual strength surfaces, are utilized to evaluate the damage level, which will be subsequently described. The compressive strength of the concrete was set to $f_c^{\prime }=37\mathrm{MPa}$ based on Schmidt hammer tests conducted after the earthquake (Sato Kogyo Co. Ltd. 1997). Poisson’s ratio and the density were assumed to be ${\nu }_c=1/6$ and ${\rho }_c=2.3\times {10}^3\mathrm{kg}/{\mathrm{m}}^3$, respectively.

By inputting these values and discretizing the concrete into cubic elements with a length of approximately 25 mm, the tensile strength was estimated; the value was approximately $1/10$ of the compressive strength, and the compressive stress and tensile stress were released at levels of approximately $-0.9\%$ strain and 0.2% strain, respectively. The elastic wave speed of concrete was estimated as $C_c=\mathrm{3,458}\mathrm{m}/\mathrm{s}$. The strain-rate effects for the concrete were ignored because the numerical analyses conducted in this study concentrated on investigating the qualitative dynamic response characteristics of the structure for upward loading.

Fig. 5(b) shows the stress-strain relationship for the axial and shear rebar used in this study. A bilinear type of constitutive model was applied considering the plastic hardening effect after yielding. The yield strength of the rebar was set to $f_y=306\mathrm{MPa}$. The density ${\rho }_s$, elastic modulus $E_s$, and Poisson’s ratio ${\nu }_s$ for the rebar were assumed using the nominal values of ${\rho }_s=7.85\times {10}^3\mathrm{kg}/{\mathrm{m}}^3$, $E_s=206\mathrm{GPa}$, and ${\nu }_s=0.3$, respectively. It was assumed that the yielding of the rebar was controlled by the von Mises yield criterion, and the plastic hardening coefficient $H_s^{\prime }$ was assumed to be 1% of the Young’s modulus $E_s$.

Fig. 5(c) shows the stress-strain relationship for the overburden soils at the station established with a cut-and-cover method and for backfill on the outer side of the walls. The model was assumed to be a perfect elastoplastic body, and the elastic modulus at unloading was $E_{\mathrm{Gul}}=10\mathrm{GPa}$; at this value, the stress was perfectly released as soon as the soil unloaded, and no resistance of the tensile force occurred. The movement of the overburden soils should follow the movement of the ceiling slab when the slab drops due to the collapse of the central column.

Assuming a state with a single pulselike uplift motion and a period $T_0$ at the base stratum, the bottom surface of the stratum was forcibly and linearly moved up to the maximum displacement $d_{i,\max }$ during the half time of period $T_0$ and was then kept at this displacement location. Fig. 6 shows the schematic displacement $d_i$-time $t$ and velocity $V_i$-time $t$ relationships. From this figure, the velocity $V_i$ of the movement is obtained as $V_i=2d_{i,\max }/T_0$, and after passing the time of $T_0/2$, the velocity $V_i$ drops to zero. Generally, at the beginning of loading, when the end of an elastic body vertically impacts a rigid body with velocity $V$, a wave propagates in the elastic body with an axial stress that is given by the following equation:where $C$ = elastic wave speed of the body.

Figs. 7 and 8 show comparisons of the time histories of the response displacement $d$ and relative displacement $d_r$ at each point. Figs. 7(a and b) and 8(a and b) are the results for the loading durations of $T_0/2=5\mathrm{ms}$ and $T_0/2=50\mathrm{ms}$, respectively. In these figures, to investigate the propagation characteristics of a wave passing through the station structure from the base stratum, each time history is indicated for points along the axes of the central column and sidewall. Reference time histories for a wave that is not influenced by the substructure are also indicated along the vertical line 10 m inside the free ends of the soil layers.

From the distribution of the reference time histories of the displacement for the soil layers shown in Fig. 7(a), it can be observed that the amplitude of the wave crest gradually decreases up to a similar depth level as that of the lower surface of the base slab and increases near the subsurface layer. The amplitude at the ground surface is almost twice that of the input displacement. This result may be due to the effect of the reflection wave when the axial stress is released at the ground surface.

By comparing the time histories of the waves between the lower point of the central column/sidewall and the upper surface point of the soil layer just below the base slab shown in Fig. 7(a), the waves at the column/sidewall monotonically increase, even though the waves behave similarly to $t=20\mathrm{ms}$ after the start of uplift loading (referred to as time $t$). However, the waves at the soil layer do not increase in amplitude in a similar manner; additionally, the soil layers vibrate alternatively with a period of approximately 60 ms.

Based on the time histories at two column points, when the input velocity is equal to or less than $V_i=1\mathrm{m}/\mathrm{s}$, the wave may propagate continuously. However, by increasing the input velocity $V_i$ by $2\mathrm{m}/\mathrm{s}$, the displacement at the lower point of the column increases up to approximately $d=130\mathrm{mm}$, but at the upper point, this value is less than $d=20\mathrm{mm}$. From Fig. 8(a), these findings are clearly observed from the relative displacement $d_r$, which refers to that at the upper point of the soil layer just below the base slab. Thus, the column collapses under axial compression failure because the average axial compressive strain of the column reaches more than $\epsilon =2.8\%$ at a 3.82-m column height. By increasing the input velocity to more than $V_i=2\mathrm{m}/\mathrm{s}$, because the amplitude of the wave at the lower point of the column increases more than that at $V_i=2\mathrm{m}/\mathrm{s}$, the column will suffer severe damage. Because the time histories of the waves between two points on the sidewall are approximately the same, the average axial strain in the wall may be very small, and the wall may not reach axial compressive failure, thus differing from the column state. This discrepancy may be inferred from the difference in axial stiffness between the column and wall.

From Fig. 8(a), because the relative displacement $d_r$ at the ground surface on the axis of the sidewall, which refers to that of the ceiling slab, is distributed in the negative region and reaches a maximum of approximately 150 mm, the displacement at the ground surface above the sidewall may be restrained significantly compared with that of the sidewall. Alternatively, the ground surface on the axis of the central column moves slightly upward because the maximum relative displacement $d_{r,\max }$ at the ground surface, which refers to that at the ceiling slab, is approximately 20 mm.

Fig. 7(b) shows the results for the case of a loading duration of $T_0/2=50\mathrm{ms}$. The distribution of the reference time histories of the waves for the soil layers show that although the amplitude of the wave at the ground surface is increased due to the effect of the reflection wave, in each soil layer, this amplitude is similar to that of the input displacement.

The distribution of the time histories of displacement along two axes shows that the wave at the upper surface of the soil layer just below the base slab cannot propagate continuously to the lower points of the central column/sidewall, as in the case of a loading duration of $T_0/2=5\mathrm{ms}$. The soil layers separate from the base slab and change to a vibration state with a period of approximately $T=150\mathrm{ms}$, and the period is prolonged from $T=60$ to 150 ms due to the loading duration extending from 5 to 50 ms.

The results of the time histories of the displacements at both the upper point and lower point of the central column indicate that when the input velocity is $V_i=0.5\mathrm{m}/\mathrm{s}$, the wave propagates continuously in the upward direction. However, when the input velocity $V_i$ is increased above $0.5\mathrm{m}/\mathrm{s}$, the amplitude of the wave at the lower point increases abruptly, but at the upper point, the amplitude does not reach half that at the lower point and the relative displacement is distributed almost in the negative region to time $t=200\mathrm{ms}$, as shown in Fig. 8(b). The strain in the column at time $t=50\mathrm{ms}$ is estimated when the input velocity $V_i=1\mathrm{m}/\mathrm{s}$ and the displacement amplitudes at the lower and upper points of the column are approximately $d=70$ and 10 mm, respectively. Notably, the average axial compressive strain is greater than $\epsilon =1.5\%$, and the column may collapse under axial compression failure. Therefore, by increasing the input velocity $V_i$, the column may experience increasingly severe damage.

Regarding the dynamic response of the sidewall, because the wave continuously propagates in the upward direction along the wall, severe damage does not occur, irrespective of the duration of loading.

From Fig. 8(b), it is observed that the relative displacements $d_r$ at the ground surface on the two axes, which refers to each distribution at the ceiling slab, are approximately 50 mm larger than those in the case of a loading duration of $T_0/2=5\mathrm{ms}$ [Fig. 8(a)]. However, the distribution characteristics are almost similar to time $t=100\mathrm{ms}$. When the input velocity is equal to or greater than $V_i=2\mathrm{m}/\mathrm{s}$, because the relative displacement $d_r$ at the bottom of the overburdened soils on the axis of the central column, which refer to that of the ceiling slab, abruptly increases at time approximately $t=250\mathrm{ms}$, the overburden soils separate from the ceiling slab and begin to substantially rise.

Fig. 9 shows comparisons of the time histories of the particle velocity $V$ at each point, which are similar to those for the displacement $d$ shown in Fig. 7. Fig. 9(a) indicates the results for a loading duration of $T_0/2=5\mathrm{ms}$. The distribution of the reference time histories for the soil layers shows that at the beginning of wave propagation, the amplitude of the wave remains near the input value in each soil layer, except at the base of the second subsurface layer and at the ground surface; the values at these locations are largely influenced by the reflection wave due to the axial stress reaching zero at the ground surface. Then, the positive-sign reflection wave propagates in the downward direction.

From the distribution of the time histories of the wave along the two other axes, because the displacement wave cannot propagate continuously at the boundary between the bottom of the base slab and the upper surface of the soil layer, as mentioned previously, wave discontinuity can be observed. Although the reflection wave with a positive sign propagates in the downward direction in the soil layers, the upper surface of the soil layer below the base slab acts as a free boundary. Because the time histories at the bottom surface points of the overburden soils and the upper points of the column/sidewall are different, the overburden soils may separate slightly from the ceiling slab at some time.

Based on a comparison of the time histories of the particle velocity $V$ at the two points on the column, when the input velocity is equal to or greater than $V_i=2\mathrm{m}/\mathrm{s}$, the amplitude of the wave at the lower point maintains an approximately similar value after separation from the soil layer and then gradually decreases. Additionally, at the upper point, the amplitude remains less than $V=1.5\mathrm{m}/\mathrm{s}$ and decreases to a negative value. This result suggests that due to the comparatively low axial stiffness of the column, the column may reach axial compression failure, and the stress wave at the beginning of loading may not directly propagate toward the upper part of the column. The axial stress in the column due to particle movement was estimated. The maximum particle velocities $V_{\max }$ in the case of the input velocities $V_i=2$, 3, and $4\mathrm{m}/\mathrm{s}$ are approximately $V_{\max }=4.5$, 7, and $10\mathrm{m}/\mathrm{s}$, respectively, and the stress that corresponds to each maximum velocity $V_{\max }$ is evaluated, assuming that the concrete column is an elastic body, using Eq. (3) with $\sigma =36$, 56, and 60 MPa. From these results, as discussed in the previous section, it is confirmed that the column may collapse with axial compression failure because the compressive strength of the concrete was $f_c^{\prime }=37\mathrm{MPa}$.

When the input velocity is $V_i=1\mathrm{m}/\mathrm{s}$, the particle velocity $V$ at the upper point of the column is approximately $V=1\mathrm{m}/\mathrm{s}$ at approximately $t=25\mathrm{ms}$ and then decreases approximately linearly to zero at time $t=40\mathrm{ms}$. Moreover, at the lower point, even though the velocity $V$ at time $t=20\mathrm{ms}$ is greater than $V=1.5\mathrm{m}/\mathrm{s}$, the velocity decreases to zero around $t=30\mathrm{ms}$. This change suggests that after time $t=25\mathrm{ms}$, the reflection wave from the upper point may not directly propagate to the lower point, potentially because horizontal circumferential cracks may occur along the column, and the stress wave may disappear, as supported by the axial stress trend at a particle velocity of $V=0.5\mathrm{m}/\mathrm{s}$ and tensile stress of $\sigma =4\mathrm{MPa}$ at time $t=30\mathrm{ms}$.

Fig. 9(b) shows comparisons of the time histories of the particle velocity $V$ for a loading duration of $T_0/2=50\mathrm{ms}$. The reference time histories of the wave with depth in the soil layers show that the distribution characteristics of the wave are similar to those for a loading duration of $T_0/2=5\mathrm{ms}$. However, the amplitude of each wave during uplift loading abruptly increases, which implies that the reflection wave propagating in the downward direction from the ground surface influences the time history for a long duration of loading. After unloading, the soil layers reach the free vibration state.

From the distributions of the time histories of the wave along the two axes, the response characteristics of the waves are similar to those in the case of a loading duration of $T_0/2=5\mathrm{ms}$, except that the amplitude of the wave at each point in the soil layers below the base slab abruptly increases, thus mimicking the reference time histories.

In the LS-DYNA code, the KCC model is applied as the constitutive model for concrete, and the damage level of a concrete element can be evaluated using an index (Malvar et al. 1997; Markovich et al. 2011). In the KCC model, three independent strength surfaces are formulated as follows:where $i=y$, $m$, and $r$, i.e., the yield strength surface, maximum strength surface, and residual strength surface, respectively; $p$ = pressure; and $a_{ji}$ ($j=0$, 1, 2) = parameters calibrated from the test data, which are set as default parameters in LS-DYNA (Livermore Software Technology Corporation 2016). The failure surface is interpolated between the maximum strength surface and either the yield strength surface or the residual strength surface. The failure surface is a function of the internal damage parameter $\lambda$, which is definedwhere the effective plastic strain increment $d{\epsilon }^p$ is given bywhere $f_t$ = tensile strength of concrete; $b_1$ and $b_2$ = damage scaling parameters for the case of uniaxial compression and tension, respectively; $r_f$ = dynamic increase factor that accounts for the strain rate effects; and $s$ = user-defined scale of the damage measure. These parameters excluding $s$ are set as default parameters in LS-DYNA (Livermore Software Technology Corporation 2016). The parameter $s$ can range from 0 to 100. In the case of $s=0$, the strain-rate effects are omitted; for $s=100$, the strain rate effects are fully included. In this study, the strain rate effects were omitted, as previously described. When the failure surface is distributed in the region between the yield strength surface and the maximum strength surface and between the maximum strength surface and the residual strength surface, the damage index was defined to move linearly from 0 to 1 and from 1 to 2, respectively. More details have been provided by Malvar et al. (1997) and Markovich et al. (2011).

In this study, it is assumed that the element reaches the ultimate state when the damage index is equal to or greater than 1.98, and elements that reach this value are colored red.

Fig. 10 shows the changes in the deformation and damage patterns in the cross section of the structure for each input velocity $V_i$ at time $t=100$ and 400 ms and loading duration $T_0/2=5$ and 50 ms; additionally, the enlargement factor for the deformation is set to 1.

Fig. 10(a) shows the results in the case of the loading duration $T_0/2=5\mathrm{ms}$. From this figure, when the input velocity $V_i=1\mathrm{m}/\mathrm{s}$, horizontal circumferential cracks form over the entire height of the central column, and the ceiling and base slab are subjected to positive and negative bending moments, respectively, based on the flexural crack patterns. Horizontal circumferential cracks may occur due to the downward propagation of the reflection wave, as mentioned in the previous section, and the cracks observed in the RC bridge piers (Fig. 11) after the Hyogoken-Nanbu earthquake may be due to the same phenomenon associated with these crack patterns. At an input velocity of $V_i=2\mathrm{m}/\mathrm{s}$, the central column collapses around midheight, and the base slab tends to move up and separate from the soil layer. Upon increasing the input velocity $V_i$, the magnitude of the axial compression failure of the column increases, and the area tends to move in the downward direction. As mentioned previously, the central column partially fails, but the sidewall does not, and only flexural cracks occur.

From Fig. 10(b), for a loading duration of $T_0/2=50\mathrm{ms}$, horizontal circumferential cracks occur at an input velocity of $V_i=0.5\mathrm{m}/\mathrm{s}$, even though an indication of this phenomenon is not clearly stated in the previous sections. Additionally, the central column reaches axial compression failure at $V_i=1\mathrm{m}/\mathrm{s}$, and the damage level of the column may increase due to the long duration of loading. Both the base slab and ceiling slab are in a state of separation from the soil layer when the input velocity is greater than $V_i=1\mathrm{m}/\mathrm{s}$.

Fig. 12 shows comparisons of the time histories of the displacement at the ground surface point and upper point of the central column, the relative displacement at the bottom of the overburden soils and upper surface of the ceiling slab, and relative displacement at the upper and lower points of the column along the column axis at the input velocity of $V_i=4\mathrm{m}/\mathrm{s}$ and loading duration of $T_0/2=5\mathrm{ms}$. Notably, the ground surface rises to a height of 90 mm because the base slab is lifted and then gradually falls to more than 600 mm below the original ground level due to the axial compression failure of the column and the effect of the self-weight of the soils. Additionally, the overburden soils separate from the ceiling slab with a maximum distance of approximately 100 mm and contact the slab again due to the effect of gravity after a time of approximately $t=500\mathrm{ms}$. Because the relative displacement between the upper point and lower point of the column increases monotonically to more than 500 mm due to uplift loading, the axial compression failure of the column proceeds continuously. Even the relative displacement maintains a constant value when the overburden soils separate from the ceiling slab. After the soils contact the slab again due to gravity, the displacement gradually increases, and the compression failure of the column proceeds.

Fig. 13 shows the temporal transitions of the deformation and damage patterns in the cross section of the station structure and overburdened soils 1,000 ms from the beginning of loading based on the case presented in Fig. 12. This figure confirms that due to uplift loading, the base slab is lifted, the central column reaches axial compression failure, and the overburden soils separate from the ceiling slab at a time of approximately $t=400\mathrm{ms}$. Due to the effect of gravity, the overburden soils contact the ceiling slab again, the cross section of the structure gradually falls with the overburden soils, and the compression failure of the column proceeds.

As shown in Fig. 3(b), the damage state of the Daikai Station subway structure reached during Hyogoken-Nanbu earthquake is described as follows: (1) the central columns were shortened due to axial compression failure; (2) the upper areas of the sidewalls and ceiling slab near the central column perfectly broke in flexure and folded out sharply and symmetrically; and (3) a maximum subsidence of the ground surface above the central column reached 2.5 m. The numerical results indicate that even though the sidewalls and ceiling slab were damaged significantly and deformed symmetrically, they have not perfectly broken in flexure and the maximum subsidence of the ground surface did not reach 2.5 m but reached approximately 0.6 m. However, because the failure characteristics of the station structure obtained from this study are approximately similar to those of the actual structure during the Hyogoken-Nanbu earthquake, the magnitude of the actual damage level may be reproduced due to the increase in the input velocity and/or the duration of displacement loading. Thus, the potential for collapse of the central column of the Daikai Station subway structure due to impulse uplift loading may be confirmed.