Simulating Collapse Behaviors of Buildings and Motion Behaviors of Indoor Components During Earthquakes Simulating Collapse Behaviors of Buildings and Motion Behaviors of Indoor Components During Earthquakes

A finite element code, which can handle large-scale collapse and motion behaviors of structural and non-structural components of buildings, was developed. The code was developed with a use of an adaptively shifted integration (ASI)-Gauss technique. It provides higher computational efficiency than the conventional code in those problems with strong nonlinearities including phenomena such as member fracture and elemental contact. Several numerical results obtained by using the numerical code are shown in this chapter: first, a seismic pounding analysis of the Nuevo Leon buildings, in which two out of the three buildings collapsed completely in the 1985 Mexican earthquake, then, a continuous analysis of a steel frame building, subjected to a seismic excitation followed by application of tsunami force, and finally collided with a debris. A motion behavior analysis of a gymnasium is followed as a numerical example, showing the behaviors of indoor components such as ceilings, which dropped occasionally due to detachment of clips and screws. Furthermore, numerical results on motion behavior analysis of furniture were validated with some experimental results.


Introduction
The Great East-Japan Earthquake, which occurred on March 11, 2011, had caused catastrophic disasters along the coast side of East-Japan mainly by a huge tsunami driven by the diastrophism. The tsunami carried different kinds of debris such as ships and cars up the stream, which caused additional damages to the buildings in the area. Moreover, the earthquake not only damaged structural components of buildings and facilities but also indoor non-structural components of buildings, which were located several hundred kilometers away from the epicenter. The indoor components, such as ceilings, doors, and furniture, could become fatal obstacles that obstruct people from evacuating and could even injure them. The motion behaviors of indoor components, along with the collapse behaviors of buildings, have become one of the main targets of recent investigations to reduce the number of victims in such events.
A finite element code that can handle both collapse behaviors of buildings and motion behaviors of indoor components is briefly described in this chapter. The code was developed with a use of an adaptively shifted integration (ASI)-Gauss technique [1]. It obtains highly accurate solutions with a small computational cost and can be easily integrated into existing finite-element codes using linear Timoshenko beam elements. It provides a higher computational efficiency than the conventional code and can handle dynamic phenomena such as member fracturing by shifting the numerical integration points and expressing the precise locations of the fractured sections. This technique provides results with high accuracy at a notably low computational cost compared to the conventional code, and it is efficient in applications in the dynamic analysis of large-scale framed structures. Moreover, frictional contact between objects was considered using a sophisticated penalty method [2] to simulate various contact phenomena of indoor components during seismic excitation.
Several numerical results obtained by using the numerical code are shown in this chapter; first, a seismic pounding analysis of the Nuevo Leon buildings, in which two out of the three buildings collapsed completely in the 1985 Mexican earthquake, then, a continuous analysis of a steel frame building, subjected to a seismic excitation followed by application of tsunami force, and finally collided with a debris. A motion behavior analysis of a gymnasium is followed as a numerical example, showing the behaviors of indoor components such as ceilings, which dropped occasionally due to detachment of clips and screws. Furthermore, numerical results on motion behavior analysis of furniture were validated with some experimental results.

ASI-Gauss code
The general concepts of the ASI-Gauss technique, which is a revised version of the ASI technique [3], are explained in this section. Figure 1 shows a linear Timoshenko beam element and the location of its numerical integration point, in comparison with a rigid bodies-spring model (RBSM) and the location of its rotational and shear springs. The relation between the locations of the numerical integration point and the rotational and shear springs (or the stress evaluation point) where a plastic hinge is formed is expressed as [4] r s = - where s is the location of the numerical integration point and r the location of the stress evaluation point, respectively. They are non-dimensional and take values between −1 and 1. The plastic hinge model represented by the RBSM can be expressed simply, by adopting the linear Timoshenko beam element whose numerical integration point is shifted appropriately by using this relation.
In both the ASI and ASI-Gauss techniques, the numerical integration point is shifted adaptively to form a plastic hinge at exactly the section when a fully plastic section is determined to be formed within an element. Moreover, the corresponding numerical integration point is shifted back to its initial location when the plastic hinge is determined to be unloaded. By using this means, the elastoplastic behavior of the element is simulated appropriately, and the converged solution is obtained with only a small number of elements per member. However, in the ASI technique, the initial location of the numerical integration point is at the midpoint of the linear Timoshenko beam element, which is considered to be optimal for one-point integration in the elastically deformed element. It causes the solutions in the elastic range become inaccurate when the number of elements per member is small. Therefore, to improve the accuracy in the elastic range, the two consecutive elements forming a member are considered as a subset, and the numerical integration points of an elastically deformed member are placed such that the stress evaluation points coincide with the Gaussian integration points of the member [1], in the ASI-Gauss technique. This means that the stresses and strains are evaluated at the locations optimal for two-point integration. In this way, one can maintain the accuracy of bending deformation by two-point integration while using one-point integration per element in the actual calculations. Figure 2 shows the locations of the numerical integration and stress evaluation points in the elastic range for the ASI and ASI-Gauss techniques. In these techniques, the elemental stiffness matrix, the generalized strain, and the sectional force increment vectors in the elastic range are given by the following equations, respectively.   (1). For example, the numerical integration point is shifted immediately to the lower end of the element (s=−1) if a fully plastic section is determined to occur at the upper end of an element (r=1), as shown in Figure 3(b). At the same step, the numerical integration point of the adjacent element forming the same member is shifted back to its midpoint (s=0), where it is an appropriate location for a one-point integration. The condition used to determine the yielding is expressed as follows:

[ ] [ ] [ ][ ]
where f y , M x , M y, and N are the yield function, the bending moments around the x-axis, the yaxis, and the axial force, respectively. The subscript 0 indicates the values that result in a fully plastic section in an element if they act independently on the cross section. Once a section of an element is determined to be yielded based on Eq. (4), the elemental stiffness matrix, the generalized strain, and the sectional force increment vectors of the element are calculated in both the ASI and ASI-Gauss techniques using the equations given below: The highly accurate elastoplastic solutions can be obtained by adaptively shifting the numerical integration point of an element according to Eqs. (2), (3) and (5).
Code verification by elastoplastic response analyses of a simple space frame subjected to stepwise loading, also by the ones subjected to seismic excitation, can be found in Refs. [1,5].
The ASI-Gauss technique converged faster than the conventional scheme and the ASI technique; the exact solution was obtained using only two elements per member. The reason why the converged solution was obtained with a minimum subdivision was that the stress evaluation points in the elements were adaptively shifted in both the elastic and plastic ranges. It was concluded that the two-element subdivision per member was sufficient to determine the overall behaviors of framed structures. The verification and validation results of the numerical code can also be found in Ref. [6].
A member fracture is determined only at the fully plastic section. When the fracture is judged to occur, the numerical integration point of the element is shifted back to its midpoint, and the stiffness is reduced to zero to deactivate the element. Sectional forces exerting at the fractured section are released instantly at the next step after fracture has occurred, which is determined by the following fracture condition examining each strain value in a member: where κ x , κ y , γ xz , γ yz , ε z , κ x0 , κ y0 , γ xz0 , γ yz0 , ε z0 are the bending strains around the x and y axes, the shear strains for the x and y axes, the axial tensile strain and the critical values for these strains, respectively. The critical strain values can be actually obtained from some experiments concerning the fracture of joint bolts [7].
Elemental contact is considered by connecting the pairs of colliding elements with gap elements, using the geometrical relations of the four nodes of the colliding elements. The gap elements are automatically eliminated after some time of contact when the mean value of deformations of gap elements is decreased to a certain ratio. More details regarding this numerical code can be found in Refs. [8][9][10][11][12][13].

Seismic pounding analysis of adjacent buildings
Seismic pounding phenomena, particularly the collision of adjacent buildings under longperiod ground motion, are becoming a more significant issue in Japan. In the 1985 Mexican earthquake, many apartment buildings such as the Nuevo Leon buildings in the Tlatelolco district of Mexico City collapsed due to long-period ground motion (Figure 4). The site was approximately 400 km away from the epicenter. The Nuevo Leon buildings consisted of three similar buildings built consecutively with narrow expansion joints between the buildings, and it was estimated that the collisions between adjacent buildings had caused the destruction.  As shown in Table 1, the natural period of the north building model was set to be 25% longer than the other periods by lowering the structural strengths of the columns; the collisions between building models did not occur in the preliminary numerical test when the natural periods were the same. The difference of natural periods was actually observed in same-type buildings built near the site, as shown in Table 2. It was assumed to be caused by the damage from previous earthquakes [14]. The EW, NS, and UD components of the SCT seismic wave, shown in Figure 6, were used as an input ground motion. The intensity period of the seismic wave was approximately 2 s and it happened to be resonant with buildings of this height.  Table 2. Natural period of same-type building built near the site. According to the obtained results shown in Figure 7, a collision first occurs between the north and center buildings due to the difference of the natural periods. Then, plastic regions spread through the beams and columns due to continuous collisions. The columns near the impact point occasionally lose their structural strengths resulting from the continuous pounding sequence, and eventually, the collapse of the center building initiates at the ceiling of the 9th floor at approximately 70 s. The north building collapses a few seconds after the center; however, the south building withstands the collisions and does not collapse.

Collapse analysis of steel frame building subjected to seismic excitation, tsunami, and debris collision
The huge tsunami occurred after the Great East-Japan Earthquake in 2011 were reported to have caused additional damages to the buildings by washing up large debris such as ships on land. Therefore, a steel frame building subjected to seismic excitation, tsunami, and debris collision was analyzed in a single continuous simulation, to see the effects of each force on the collapse behavior of the building.
As shown in Figure 8 (Figure 9), which was observed in the Great East-Japan Earthquake, was applied only to the building as an input. It was applied for 150 s from the start of the seismic activity.  Then, buoyant forces were applied statically, and drag forces were applied dynamically for the next second, at the locations indicated in Figure 8. The buoyant force applied to the debris balances with the weight and it is expressed as b s F gV r = Simulating where ρ S is the density of sea water including debris which is 1200 kg/m 3 , g the gravitational acceleration, and V the volume of water removed by an object. In the case of the building model with walls under the water line, we used the big box volume under water, as the volume of the water removed. However, in the case of the model without walls under the water line, we used only the volume of columns and beams under the water line. The drag force relates to the relative velocity between tsunami and the object and it is expressed as where A is the projected area under water line, C d the drag force coefficient, and U the relative velocity between tsunami and the object. We used all area below the water line as the projected area for the model with walls under the water line, and side face area of columns below the water line as the projected area for the model without walls under the water line. The overall drag force acts on the model with walls becomes about 10 times larger than that of the model without walls.
At the third phase of the analysis, an initial velocity, which is equal to the velocity of the tsunami, was given to the debris model and the collision to the building was simulated. The drag force applied to the debris was initially zero because the relative velocity between the debris and tsunami was zero. On impact, the value reached to about 3.3 MN when the debris perfectly stopped its motion. For other numerical conditions, we set the incremental time of the dynamic analysis to 1 ms, and we used Newmark's β method with numerical damping as the time integration scheme.   Figure 11 shows the behavior of the building subjected to tsunami and debris collision for the model with walls under the water line. Because of the walls, the drag force applied is very large and it nearly gives a critical damage to the building. The debris collision, however, gives the last trigger to the building to be washed away. After the collision, the building and the debris move with the same velocity and reach to the same value with the tsunami itself. The relative velocity becomes small, and eventually, the drag forces of both decrease to zero as shown in Figure 12.   Figure 13 shows the behavior of the building subjected to the tsunami and debris collision for the model without walls under the water line. Since there are no walls under water, the building withstands the fluid force in the first phase. As shown in Figure 14, the drag force applied to the building is smaller than that of the debris in this case. Although a part of the debris separates in the last phase, both models integrate and stop their motion after the impact. Furthermore, some head-on collision cases as shown in Figure 15 were carried out; these cases can make the drag force acting on the debris smaller while the impact force on the building may become more intensive. Figure 16 shows the relationship between maximum story drift angle of the building and the velocity of the tsunami wave in different inundation heights (d) for the model with walls under the water line. The broken line in the figure indicates the estimated story drift angle of severe damage, 1/30, which is determined in the structural design guideline in Japan. If the inundation height is low (Figure 16(a)), the maximum story drift angle of the building becomes larger in every case for the head-on collision. A head-on collision to lower levels of the building gives severe damages to the columns in lower levels and this makes the maximum story drift angle larger. However, there are some cases in higher inundation heights (Figure 16(b)), where the damage becomes larger if the debris hits the building sideways and in higher velocity. This is due to the drag force which acts on a broadly projected area of the debris in a very high speed flow of tsunami. As can be seen in Figure 17, the case of the model without walls under the water line shows similar tendency to that of the former case with walls, although the maximum story drift angle becomes smaller compared to those in same velocities of the former case. In the case of the tsunami with lower inundation heights, the damage of the building becomes larger if a headon collision occurs. In the case of a big tsunami with high inundation height and high velocity, the damage of the building becomes larger if a sideway collision occurs.  By regaining this information, an appropriate design for tsunami refuge buildings should be estimated and a means to avoid the approach of large debris to such buildings should be regarded at the same time. Figure 18(a) shows a ceiling collapse damage occurred inside a gymnasium during the Great East-Japan Earthquake in 2011. One can observe the plasterboards which weigh about 10 kg each fell and scattered all over the floor. The prevention of these phenomena is an important issue not only to save people's lives but also to keep these facilities to be safely used as shelters after earthquakes.

Ceiling collapse analysis of a gymnasium
The components of the suspended ceiling are quite complicated as shown in Figure 18(b). The clips, indicated in red circles in the figure, connect ceiling joists and ceiling joist receivers. These are small and delicate components and they may be detached during repeated excitation. Once there is a local detachment of clips, a change in the load distribution may cause a chain reaction of detachments, which ends in a drop of plasterboards. Furthermore, the detachments of hanging bolts, which are connected to the structural members composing the roof, and failure of screws on plasterboards are assumed to be other main causes of the ceiling collapse. There were some preliminary tests done to see the actual strengths of these components as shown in Figure 19. For example, the opening of clips began at a tensile force of 0.4 kN [15], and the opening of hangers at 2.8 kN [16]. For the screws connecting plasterboards to ceiling joists, they were known to be shear damaged at a value of 0.3 kN [17], and their heads became loose when 0.2 kN of tensile force was applied [18]. These criteria were implemented in the analysis. However, the strengths of clips were overestimated, and larger criteria of 0.8 kN were adopted here because they tend to slip on ceiling joists during excitation.
Figure 20(a) shows a ceiling model in detail. As we used the ASI-Gauss code, all components (except braces and hanging bolts) were modeled with two linear Timoshenko beam elements per a member. Hangers and hanging bolts were modeled in one piece. The plasterboards were assumed as rigid in outer plane direction and only the mass of rock wool boards was considered. Their strengths were neglected. Clips and screws were modeled with minute, small elements. Each plasterboard was modeled separately to consider local contact between plasterboards, which was simulated by modeling the screws slightly apart. Figure 20(b) indicates a gymnasium model with a ceiling. Elastoplastic buckling of braces and hanging bolts were taken into account by modeling them with eight beam elements each and two hinge elements on both ends. The numerical result was validated with the experimental result, which was performed by the E-Defense under an input wave of K-NET Sendai 50% wave. We used the acceleration data obtained on the shake-table floor of the E-Defense, as shown in Figure 21, as an input wave for the analysis.    The acceleration responses and the spectrum obtained at the location (X4-Y3) on the roof indicated in Figure 20(b), for example, matched well with the experimental result as shown in Figures 22 and 23. A 10 Hz high cut filter was used in the post-processing procedure of the data. The displacement responses obtained at the location (X4-Y4), for example, also matched well with the experimental result as shown in Figure 24. The phase and the amplitude agreed well as shown in Figure 25, in the extracted responses for 15-25 s.     Figures 26 and 27 show the collapse behavior of the ceiling from the front view and the global view of the gymnasium. The plasterboards near walls pattered down occasionally at the first peak of the input wave. These are due to detachment of clips and screws caused by collisions with the walls. Then, the clips near rooftop began to loose due to buckling of hanging bolts caused by vertical excitation, which ends, at the second peak of the input wave, in a drop of plasterboards in a wide range.

Motion behavior analysis of furniture during earthquakes
Improperly secured furniture, especially on the upper floors of high-rise buildings under longperiod ground motion, could become dangerous objects for human life, even if the building itself had no damages at all. Many tumbled furniture such as chairs and tables in schools could become fatal obstacles that obstruct children from evacuating. The motion behaviors of furniture, as well as the collapse behaviors of building itself, should be examined to save people's lives during big earthquakes. Here, some motion behavior analyses of furniture are conducted using the ASI-Gauss code with a contact algorithm, and the results are compared with the experimental results performed on a shake-table.
To simulate various contact phenomena between furniture and walls or floors during seismic excitation, frictional contact between objects was considered using a sophisticated penalty method [2]. Figure 28 shows the subjected forces and the geometrical relations between two elements when they approach each other with a relative velocity v. Once the current distance l between the central lines of two elements becomes shorter than the mean value L of the member widths of the elements, the penalty force vector F P , expressed as follows, is assumed to act in the normal direction of the contact plane. (9) where α is the penalty coefficient, q the penalty index, and n the normal vector at the contact surface. Then, the frictional force vector F D is assumed to act not only in the tangential direction of the contact plane but also in the normal direction; the frictional force vector is decomposed into the following two components.
where μ is the dynamic friction coefficient, D c the coefficient related to damping for normal direction, v T the tangent component of v, and v N the normal component of v. The component in the normal direction F N acts as a damping force along the contact depth and contributes to maintaining numerical stability. Selection of contact parameters in the developed code is easier compared to the distinct element method (DEM). The penalty coefficient and dynamic friction coefficients used in Eqs. (9) and (10) were fixed based upon quite simple rules. The penalty coefficient α were set as the weight of each furniture, the dynamic friction coefficients were set as 80% values of each static friction coefficient, and the coefficient related to damping D c were set as 120% values of the penalty coefficient α. We applied these rules throughout all the analysis.
A set of furniture, shown in Figure 29, were used for the numerical analyses and the experiments. There were five furnitures altogether from a heavy cabinet of about 44 kg to a light office chair with casters of about 6 kg. The sizes and locations of the center of gravity of each furniture are shown in Table 3.   Table 3. Sizes and locations of center of gravity of each furniture. Figure 30 shows the arrangement of furniture on a shake-table and their finite element models constructed using linear Timoshenko beam elements. The motions of each furniture were measured by using a motion capture system with six infra-red cameras arranged around the shake-table. The contact parameters for each furniture were set as shown in Table 4. These parameters were determined according to the rules described earlier.

Conclusion
The ASI-Gauss code described in this chapter requires numerical models to be constructed with linear Timoshenko beam elements. Therefore, it has many limitations. However, the computational cost is very low, and one can practically compute the collapse behaviors of buildings and motion behaviors of indoor non-structural components by using any personal computers with smaller memories.
Some of the applications shown in this chapter had succeeded to obtain much new information that no conventional methods had ever acquired. The seismic pounding analysis had shown possibilities of collisions and fatal destruction of adjacent buildings with different natural periods. The numerical result explained a collapse mechanism of the Nuevo Leon buildings which actually collapsed during the 1985 Mexican earthquake. The analysis of a steel frame building in tsunami had shown the effects of seismic force, tsunami force, and impact force of debris to the damage of the building. The numerical result explained the importance of considering tsunami force along with the impact force of debris, in the structural design of a tsunami refuge building. The ceiling collapse analysis had shown complicated behaviors of ceiling in a gymnasium structure by considering local detachments of clips and screws. It could be used to determine the future design of an anti-seismic ceiling. The motion behavior analysis of furniture had shown the possibilities of applying finite element methods in such analysis and of obtaining practical results in various seismic waves by fixing contact parameters based upon very simple rules. It could help people to design new anti-seismic measures for furniture and could also contribute to disaster education.