Hybrid Framework for Simulating Building Collapse and Ruin Scenarios Using Finite Element Method and Physics Engine

: Reliable and high-ﬁdelity virtual ruin scenarios for collapsed buildings are essential for post-earthquake emergency search and rescue training. However, the existing research on the distribution of ruins caused by building collapse is insu ﬃ cient for supporting post-earthquake rescue training. Therefore, this paper proposes a hybrid framework for simulating building collapse and ruin scenarios, using a ﬁnite element (FE) model and a physics engine. Based on this framework, the following methods are proposed: (1) geometric model conversion from the FE model to the physics engine; (2) determination of the initial moment of collapse; and (3) data mapping of the FE simulation results. In addition, a corresponding program, Finite Element Method to Rigid Body Dynamics (FEM2RBD), is developed for the hybrid framework. The proposed framework simulates the entire process of building collapse and the distribution of ruins. The accuracy of the framework is validated using a shaking table test of a three-story reinforced concrete frame. The collapse process and ruin scenario of a real-world library building is simulated as a case study. The results show that the proposed framework combines the advantages of the FE model during the small-deformation stage with the advantages of physics engines during the large-deformation stage. The proposed framework can be valuable in simulating building collapse and ruin scenarios for post-earthquake rescue training.


Introduction
Earthquakes can cause building collapse, and people may be buried by ruins [1,2]. Relevant statistics show that being buried is the most common cause of earthquake-related injuries [3]. Therefore, it is necessary to conduct training to improve rescue efficiency and mitigate casualties after earthquakes. To achieve this, the ability to accurately simulate ruin scenarios is required. The reliability of simulated ruins is also an essential factor for efficient training [4].
To date, several studies on the collapse simulation of buildings have been conducted using the finite element (FE) method. Such studies typically focus on the collapse resistance of structures, rather than the formation of ruins, such as the analyses of the progressive collapse resistance caused by accidental load [5][6][7][8], the simulations of seismic collapse processes, etc. [9][10][11][12][13]. Few studies that simulate the distribution of ruins induced by structural collapses have been conducted, due to the high computational workload, and the difficulty of simulating large deformations and nonlinear behaviors using the FE method [14,15].
Although the physics engine-based building collapse simulation methods have advantages in simulating large structural deformations, the simulation results at the small-deformation stage are less reliable. For example, due to the difficulty of determining constraint damping in Bullet, and the oversimplification of the rebar model in BCB, the simulation results of the BCB and FE methods are quite different at the small-deformation stage. However, the process of the structural collapse is often sensitive to the initial state of the collapse. Therefore, ensuring the accuracy of collapse mode simulations using physics engine-based methods is difficult.
To address the challenges mentioned above, a hybrid framework for simulating building collapse and ruin scenarios, based on an FE model and a physics engine, is proposed. For dynamic responses during the small-deformation stage, an FE model is used because of its accuracy in such simulations. For the large-deformation stage, a physics engine-based method is used. In this work, the MSC.Marc software [33] is used to perform the FE modeling and simulations, and the BCB is used to simulate the ruin scenario. The proposed framework is also applicable to other FE software and physics engines.
In Section 2, the proposed hybrid framework is introduced. Section 3 will illustrate three key methods adopted in the proposed framework: (1) conversion of the geometric data from the FE method to the physics engine; (2) determination of the initial moment of collapse; and (3) data mapping of the FE simulation results. Section 4 will validate the feasibility and accuracy of the proposed framework using a shaking table test of a three-story reinforced concrete frame. The collapse process and ruin scenario of a real-world library building are simulated as a case study to demonstrate the high-fidelity collapse process during an earthquake. The outcome of this research contributes to the construction of virtual building collapse and ruin scenarios, and supports post-earthquake emergency rescue training. Figure 1 shows the proposed hybrid framework, where the software programs are in bold italics, and the detailed functions of the program are in italics. The framework consists of four modules:

•
Module 1: FE model establishment and conversion The FE model of the target structure is established using the FE software MSC.Marc [33]. The geometric model is then converted into the Appl. Sci. 2020, 10, 4408 3 of 19 rigid body model in Blender, using the geometry mapping function that is integrated into the Blender Python script "Finite Element Method to Rigid Body Dynamics (FEM2RBD)". In the rigid body model, the elements are rigid bodies, and no internal stresses and strains are considered. The deformations of the reinforced concrete structures are considered via the constraints/springs among different rigid bodies [29].

Data Integration and Rendering
Rendering the Result using Blender   • Module 2: Small deformation analysis and initial collapse moment determination A nonlinear time history analysis is conducted using the FE model of the target building, to analyze the structural behavior during the small-deformation stage. Subsequently, the initial moment of collapse is determined, and the displacement and velocity data of each time step at, and before, the moment of initial collapse are extracted from the FE model. These data will be used in the large-deformation simulation based on the physics engine in Module 3 and the collapse process visualization in Module 4.
• Module 3: Data mapping and large deformation analysis Based on the mapping method of the FE results proposed in Section 3, the displacement and velocity data from the moment of initial collapse are mapped into the BCB geometric model, using the displacement and velocity mapping functions of FEM2RBD. Subsequently, this model is modified using several BCB processing functions, in the following order: (1) establish the rigid body ground; (2) remove the overlapping portions of rigid bodies; (3) map the remaining ground motion; (4) define the constraint parameters; and (5) build the constraints among rigid bodies. The analysis model used for the physics engine simulation can then be established [29]. After establishing this model, simulations are performed using the Bullet physics engine to determine the structural dynamic behavior during the large-deformation stage.
• Module 4: Data integration, rendering and visualization After the large-deformation stage simulation, the render engines of Blender are used to render the data of the small and large-deformation stages. Subsequently, the rendered videos are integrated to visualize the entire structural collapse process.

Key Techniques
To implement the proposed framework, the following three key issues need to be addressed.

•
Conversion of the geometric model from the FE method to the physics engine. The model in the physics engine-based simulation component is a rigid body spring model, the elements of which are significantly different from the elements of the FE model. Therefore, it is necessary to propose a method for converting the geometric data from the FE model to the physics engine • Determination of the initial moment of collapse. To achieve the proposed hybrid framework, the switching moment (collapse keyframe) from the FE simulation to physics engine simulation needs to be determined based on the results of the FE simulation.

•
Mapping the FE simulation data. After the FE simulation, the data computed for the small-deformation stage need to be mapped into the physics engine.

Geometric Model and Material Mapping Method
Fiber beam-column and multi-layered shell elements have been widely used to simulate beams, columns and shear walls during the nonlinear dynamic analyses of structures [10,34,35]. Compared with models that use solid elements, models that use fiber beam-column and multi-layered shell elements reduce the number of elements and the computational workload required, which makes the simulations of large-scale structures easier to conduct. In this study, soil-structure interaction (SSI) effects are neglected, assuming that that the soil underneath the foundations is infinitely rigid. It is worthwhile noting that previous studies have illustrated that SSI effects can significantly affect the lateral displacement demand and collapse capacity of building structures [36][37][38]. The classical Rayleigh damping with a damping ratio of 5% based on the initial stiffness of the building is used to model the inherent damping of the structure. Future studies should consider incorporating more accurate approaches to modeling inherent damping, as a few research works have identified drawbacks associated with the use of Rayleigh damping based on initial stiffness [39][40][41][42]. In the FE model, the geometric model of a fiber beam-column element is a line [a one-dimensional (1D) element], and the geometric model of a multi-layered shell element is a plane (a 2D element). However, the geometric model in the physics engine is a solid (a 3D object). Because the geometric models of the physics engine and FE model are different, a method is proposed to convert the FE models of the fiber beam-column and multi-layered shell elements into solid objects for the physics engine.

Solid Model Establishment Method in Blender
The information available from FE models consists of geometric information (including the coordinates of element nodes and sections) and material information about different elements. After data preparation, a physics engine model for large deformation analysis is established in Blender, as shown in Figure 2.
model the inherent damping of the structure. Future studies should consider incorporating more accurate approaches to modeling inherent damping, as a few research works have identified drawbacks associated with the use of Rayleigh damping based on initial stiffness [Error! Reference source not found.-Error! Reference source not found.]. In the FE model, the geometric model of a fiber beam-column element is a line [a one-dimensional (1D) element], and the geometric model of a multi-layered shell element is a plane (a 2D element). However, the geometric model in the physics engine is a solid (a 3D object). Because the geometric models of the physics engine and FE model are different, a method is proposed to convert the FE models of the fiber beam-column and multi-layered shell elements into solid objects for the physics engine.

Solid Model Establishment Method in Blender
The information available from FE models consists of geometric information (including the coordinates of element nodes and sections) and material information about different elements. After data preparation, a physics engine model for large deformation analysis is established in Blender, as shown in Figure 2. The element establishment method in Blender is different from that in FE software. In Blender, "primitive" mesh shapes should be established first, then all the elements should be transformed from the "primitive" mesh shapes. However, as the shape and position of "primitive" mesh shapes do not coincide with the corresponding FE elements, the shape and position transformations are performed to get the target elements in FE models. Blender provides several "primitive" mesh shapes to begin modeling, including cubes, planes and cylinders [Error! Reference source not found.]. Because a cube is a shape that best approximates engineering components and FE model meshes, it is selected as the basic mesh shape in Blender. The geometric model mapping is then performed in two steps.
• Step 1: Shape and position transformation (including translation and rotation).
The shape and position transformations are performed on the cube mesh provided in Blender [Error! Reference source not found.], to create geometric objects with the same shapes and positions as the elements in the FE model. The rotation method shown in Figure 3 is used to model inclined elements in FE models. Specifically, the rotation angle is obtained by first determining the endpoint coordinates of the longest edge of an FE element (line segment l1) and the projection of l1 in the XOY plane in the global coordinate system (line segment l2). The angle α between l2 and the X axis, and the angle β between l2 and l1, are then calculated using these endpoint coordinates. Subsequently, the shape-transformed object is rotated by angle α around the Z axis in the global coordinate system. The object is then rotated by angle β around the Y' axis in the rotated coordinate system. The rotated coordinate system is the local coordinate system of the cube after the first rotation, where the X'O' axis is parallel to l2.  The element establishment method in Blender is different from that in FE software. In Blender, "primitive" mesh shapes should be established first, then all the elements should be transformed from the "primitive" mesh shapes. However, as the shape and position of "primitive" mesh shapes do not coincide with the corresponding FE elements, the shape and position transformations are performed to get the target elements in FE models. Blender provides several "primitive" mesh shapes to begin modeling, including cubes, planes and cylinders [43]. Because a cube is a shape that best approximates engineering components and FE model meshes, it is selected as the basic mesh shape in Blender. The geometric model mapping is then performed in two steps.

•
Step 1: Shape and position transformation (including translation and rotation). The shape and position transformations are performed on the cube mesh provided in Blender [43], to create geometric objects with the same shapes and positions as the elements in the FE model. The rotation method shown in Figure 3 is used to model inclined elements in FE models. Specifically, the rotation angle is obtained by first determining the endpoint coordinates of the longest edge of an FE element (line segment l 1 ) and the projection of l 1 in the XOY plane in the global coordinate system (line segment l 2 ). The angle α between l 2 and the X axis, and the angle β between l 2 and l 1 , are then calculated using these endpoint coordinates. Subsequently, the shape-transformed object is rotated by angle α around the Z axis in the global coordinate system. The object is then rotated by angle β around the Y' axis in the rotated coordinate system. The rotated coordinate system is the local coordinate system of the cube after the first rotation, where the X'O' axis is parallel to l 2 .

•
Step 2: Add rigid bodies to the cubes. Convert the cubes into rigid bodies after the shape and position transformations to obtain the rigid body model for the Bullet physics engine simulation.
Convert the cubes into rigid bodies after the shape and position transformations to obtain the rigid body model for the Bullet physics engine simulation.

Material Mapping Method in the BCB
The BCB will group objects according to their names in Blender, and assign the same material properties to the objects in each group during the constraint building process. When establishing a geometric model in Blender, the object names need to be chosen so that objects of the same material in the FE model are mapped into the same material group by BCB. The objects in Blender are named according to the rule of "element type: material type: element number." For example, an object derived from a beam element with material number 1 and element number 234 is named "Beam: 1: 234." More details are available in the BCB manual [Error! Reference source not found.].

Determination of the Initial Moment of Collapse
According to the proposed hybrid framework, the FE simulation is conducted before the initial moment of collapse Tc, and the physics engine simulation is conducted after Tc. Therefore, the determination of Tc is essential. Several collapse criteria have been proposed [Error! Reference source not found.-Error! Reference source not found.]. There are various structural collapse modes, including lateral collapse, vertical collapse, and so on. No matter which collapse mode is, the vertical displacements of structural members will exceed a specific value when the structures collapse [Error! Reference source not found.,Error! Reference source not found.]. Therefore, the definition of collapse depends on the moment when "the deformation of the structure is insufficient to maintain a safe use space" [Error! Reference source not found.]. Besides, Zhao et al. [Error! Reference source not found.] studied four collapse criteria, including the following: Criterion 1) According to Chinese seismic code [Error! Reference source not found.], the earthquake-induced inter-story drift ratio of reinforced concrete frames should not be greater than 1/50; Criterion 2) Federal Emergency Management Agency (FEMA) 356 [Error! Reference source not found.] recommended that the inter-story drift ratio limit for the collapse-prevention performance level of concrete structures be 4%; Criterion 3) According to Villaverde [Error! Reference source not found.], when the tangent slope of the incremental dynamic analysis (IDA) curve is lower than 20% of the slope of the initial elastic stage, the structure will collapse; Criterion 4) "The vertical collapse displacement of the main structural members exceed [

Material Mapping Method in the BCB
The BCB will group objects according to their names in Blender, and assign the same material properties to the objects in each group during the constraint building process. When establishing a geometric model in Blender, the object names need to be chosen so that objects of the same material in the FE model are mapped into the same material group by BCB. The objects in Blender are named according to the rule of "element type: material type: element number." For example, an object derived from a beam element with material number 1 and element number 234 is named "Beam: 1: 234." More details are available in the BCB manual [29].

Determination of the Initial Moment of Collapse
According to the proposed hybrid framework, the FE simulation is conducted before the initial moment of collapse T c , and the physics engine simulation is conducted after T c . Therefore, the determination of T c is essential. Several collapse criteria have been proposed [44][45][46][47][48][49]. There are various structural collapse modes, including lateral collapse, vertical collapse, and so on. No matter which collapse mode is, the vertical displacements of structural members will exceed a specific value when the structures collapse [47,48]. Therefore, the definition of collapse depends on the moment when "the deformation of the structure is insufficient to maintain a safe use space" [47]. Besides, Zhao et al. [46] studied four collapse criteria, including the following: Criterion 1) According to Chinese seismic code [44], the earthquake-induced inter-story drift ratio of reinforced concrete frames should not be greater than 1/50; Criterion 2) Federal Emergency Management Agency (FEMA) 356 [45] recommended that the inter-story drift ratio limit for the collapse-prevention performance level of concrete structures be 4%; Criterion 3) According to Villaverde [49], when the tangent slope of the incremental dynamic analysis (IDA) curve is lower than 20% of the slope of the initial elastic stage, the structure will collapse; Criterion 4) "The vertical collapse displacement of the main structural members exceed[ing] 1 m", is considered as Criterion 4 [46]. The work of Zhao et al. [46] shows that the conventional collapse criteria (i.e., Criteria 1-3) significantly underestimate the collapse resistance of the frames. Compared with the other three collapse criteria, Criterion 4 can better reflect the extreme nonlinear behavior of the collapse limit state. Therefore, collapse Criterion 4, proposed by Zhao et al. [46], is selected as a rule for determining T c in this work.

FE Simulation Data Mapping Method
The published research on combining FE models and physics engines is limited [26,28]. These studies also lack applicable data mapping methods for converting FE-based numerical models that use fiber beam-column and multi-layered shell elements into physics engine models. Therefore, a method for mapping simulation data from the FE model into the Blender model is proposed.

Displacement Mapping Method
The keyframe method of Blender [43] is used to achieve the displacement mapping. A keyframe is a location on a timeline that marks the beginning or end of a transition. For example, the displacement keyframe defines the beginning and endpoints of the displacement. Meanwhile, a kinematic keyframe defines the information about the kinematic type. There are two kinematic types in Blender: (1) the animation system, in which the motion of the object is controlled by the displacement keyframe, and (2) the physics engine system, in which the motion of the object is controlled by the simulation results obtained using the physics engine. Displacement mapping requires only the displacement keyframe. Before the initial moment of collapse T c , the motions of objects in Blender are controlled by the displacement keyframes derived from the FE results. This means that only Blender is used to map the FE results; the Bullet physics engine is not involved.
Specifically, the displacement mapping method is divided into three steps, as shown in Figure 4. Step 1: Extract FE Model Data Step 2: Insert Displacement Keyframe Step 3: Change Kinematic Type

Velocity Mapping Method
When the FE simulation is converted to the Bullet physics engine simulation at the initial moment of collapse Tc, both displacement and velocity should be consistent between the two simulations. Consequently, velocity mapping is necessary. The velocity mapping method requires both the displacement and the kinematic type keyframes [Error! Reference source not found.]. In Blender, the velocity of an object is defined as follows: At frame i, convert the kinematic type of an object to the animation system and insert the kinematic type as a keyframe. At frame j (where i < j), convert the kinematic type of the same object to the physics engine system, and insert the kinematic type as a keyframe again. Consequently, during frame i and frame j, the movement of the object will

•
Step 1. Extract the FE model data. After the FE simulation, the displacement data of each vertex at each time step are extracted, and then the deformed coordinates of each vertex are calculated from the displacement data and the initial coordinates. The deformed centroid coordinates of each object at each time step can be calculated from the deformed vertex coordinates of each object.

•
Step 2. Insert the displacement keyframe. Initially, the keyframe that corresponds to the first time step in the FE model on the Blender timeline is selected. The deformed centroid coordinates obtained from the FE results for one object are then calculated. Subsequently, the deformed centroid coordinates are inserted as the displacement keyframe. When all the objects in the structure have been processed, the simulation moves to the next time step, and the aforementioned procedure is repeated until the initial moment of collapse T c .

•
Step 3: Change the kinematic type. After Step 2, the kinematic type of each object needs to be converted to the animation system, otherwise the object motions will be controlled by the physics engine simulation, rather than by the FE result.

Velocity Mapping Method
When the FE simulation is converted to the Bullet physics engine simulation at the initial moment of collapse T c , both displacement and velocity should be consistent between the two simulations. Consequently, velocity mapping is necessary. The velocity mapping method requires both the displacement and the kinematic type keyframes [43]. In Blender, the velocity of an object is defined as follows: At frame i, convert the kinematic type of an object to the animation system and insert the kinematic type as a keyframe. At frame j (where i < j), convert the kinematic type of the same object to the physics engine system, and insert the kinematic type as a keyframe again. Consequently, during frame i and frame j, the movement of the object will be controlled by the linear interpolation of the displacement at frame i and frame j, respectively. After frame j, the movement of the object will be controlled by the Bullet physics engine, and the object will gain speed v j at frame j, as follows: where i and j (i < j) denote the keyframe number; dt represents the time interval between two adjacent keyframes, and in Blender, the default value is dt = 1 24 s; x i and x j represent the centroid coordinates of an object at frame i and frame j, respectively. The speed v j is the average velocity represented by the secant line, rather than the actual instantaneous velocity represented by the tangent line, as shown in Figure 5.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 20 next frame is determined by the velocity, position and force of the adjacent previous frame, as shown in Equations (2) and (3).
where Fext denotes external forces defined by the user; Fc denotes the constraint forces determined by the position and the velocity of the object in the adjacent previous frame.
If the average velocity j v is used instead of the instantaneous velocity vj at the initial moment of collapse Tc, an error due to approximating the velocity will be introduced. To improve the accuracy of the velocity mapping, this work proposes a velocity mapping method that uses a virtual displacement vector xvirtual. The virtual displacement is obtained using Equation (4), based on the velocity and displacement at the collapse keyframe j. The proposed mapping method is shown in Figure 5. The blue curve is the displacement curve simulated by the FE method; ti and tj are the times corresponding to the displacement keyframe insertion points; xi and xj are the displacements corresponding to the displacement keyframe insertion points; vj is the instantaneous velocity, and xvirtual is the virtual displacement vector. Consider the virtual displacement vector: where xj and vj represent the displacement vector of the FE model and the velocity vector of the In Blender, the Bullet physics engine adopts the symplectic Euler algorithm to update the motion of a rigid body [50], where the position of the object in the next frame is determined by the velocity, position and force of the adjacent previous frame, as shown in Equations (2) and (3).
where F ext denotes external forces defined by the user; F c denotes the constraint forces determined by the position and the velocity of the object in the adjacent previous frame. If the average velocity v j is used instead of the instantaneous velocity v j at the initial moment of collapse T c , an error due to approximating the velocity will be introduced. To improve the accuracy of the velocity mapping, this work proposes a velocity mapping method that uses a virtual displacement vector x virtual . The virtual displacement is obtained using Equation (4), based on the velocity and displacement at the collapse keyframe j. The proposed mapping method is shown in Figure 5. The blue curve is the displacement curve simulated by the FE method; t i and t j are the times corresponding to the displacement keyframe insertion points; x i and x j are the displacements corresponding to the displacement keyframe insertion points; v j is the instantaneous velocity, and x virtual is the virtual displacement vector.
Consider the virtual displacement vector: where x j and v j represent the displacement vector of the FE model and the velocity vector of the collapse keyframe j, respectively; k is the keyframe number of x virtual , and the value of k is defined by Equation (5). After testing, the velocity cannot be mapped successfully in Blender when n is less than four. The real instantaneous velocity and displacement of the FE model are mapped in the Bullet physics engine by replacing x i with the virtual displacement vector x virtual in Equation (1). In summary, the velocity mapping method can be divided into the following three steps, which are schematically diagrammed in Figure 6.

•
Step 2: Obtain xvirtual using Equations (4) and (5), and set the displacement of the frame k to xvirtual.

•
Step 3: At frame k, convert the kinematic type to the animation system, and then insert the displacement and kinematic type keyframes for each object. Convert the frame from frame k to frame j after the keyframe insertion at frame k. At frame j, convert the kinematic type to the physics engine system, and then insert the displacement and kinematic type keyframes for each object. Keyframe k(k<j) Figure 6. Flowchart of the velocity mapping method.

•
Step 1: Extract the FE results at time step t j , including velocity v j and displacement x j . • Step 2: Obtain x virtual using Equations (4) and (5), and set the displacement of the frame k to x virtual . • Step 3: At frame k, convert the kinematic type to the animation system, and then insert the displacement and kinematic type keyframes for each object. Convert the frame from frame k to frame j after the keyframe insertion at frame k. At frame j, convert the kinematic type to the physics engine system, and then insert the displacement and kinematic type keyframes for each object. To validate the reliability of the proposed hybrid framework, a 3D shaking table test was conducted with a three-story reinforced concrete frame [20] as a case study. The collapse of the structure was simulated using three methods: (1) FE simulation, (2) BCB simulation, and (3) the proposed hybrid framework.

Validation Using a 3D Shaking
The details of the structure are shown in Figure 7. The concrete information is shown in Table 1, and the rebar information is shown in Table 2. The structure was subjected to the amplitude-scaled El-Centro ground motion record during the shaking table test, and the load cases are shown in Table 3. The 5%-damped pseudo-acceleration spectra of Load Case 5 are shown in Figure 8. More details can be found in [20]. The first mode is a translation with a fundamental period (elastic period) T 1 = 0.316 s; the second mode is a translation with a fundamental period T 2 = 0.316 s; and the third mode is the planar torsion, with a fundamental period T 3 = 0.178 s.

Comparison of Simulation Results
Based on the information above, an FE model with fiber beam-column and multi-layered shell elements was established in MSC.Marc. Details of these element models have been reported in the published work [35]. The seismic responses of the structure under Load Cases 1-3 are used to validate the rationality of the FE model. Figure 9 shows the comparison of the maximum horizontal displacement in the X-direction between the test results and the FE simulation. The results indicate that the FE simulation agrees well with the test results.
In Load Case 5, a convergence problem occurred in the FE simulation (t = 4.16 s); therefore, the distribution of the ruins could not be obtained. The collapse process predicted by the FE simulation is shown in Figure 10. The colored contours in Figure 10 represent the longitudinal reinforcement yielding in the elements. The numerical simulation reveals that failures of the frame are initiated at both ends of the columns where the bending moment is large. Subsequent failures take place on the ground story, where the columns carry the largest lateral force. Large displacements occurred on the ground story, which leads to the collapse of the entire structure. In Figure 10c, the deformed shape of the structure at 4.16 s is shown.
Based on the information above, an FE model with fiber beam-column and multi-layered shell elements was established in MSC.Marc. Details of these element models have been reported in the published work [Error! Reference source not found.]. The seismic responses of the structure under Load Cases 1-3 are used to validate the rationality of the FE model. Figure 9 shows the comparison of the maximum horizontal displacement in the X-direction between the test results and the FE simulation. The results indicate that the FE simulation agrees well with the test results. In Load Case 5, a convergence problem occurred in the FE simulation (t = 4.16 s); therefore, the distribution of the ruins could not be obtained. The collapse process predicted by the FE simulation is shown in Figure 10. The colored contours in Figure 10 represent the longitudinal reinforcement yielding in the elements. The numerical simulation reveals that failures of the frame are initiated at both ends of the columns where the bending moment is large. Subsequent failures take place on the ground story, where the columns carry the largest lateral force. Large displacements occurred on the ground story, which leads to the collapse of the entire structure. In Figure 10c, the deformed shape of the structure at 4.16 s is shown. From the FE results, the initial moment of collapse is 3.88 s, using the collapse criterion in Section 3.2. A ¼-scale model is adopted; thus, the vertical collapse displacement of the collapse criterion is 0.25 m.
Subsequently, collapse simulation was performed using the proposed hybrid method. Figure  11b shows the distribution of ruins, which shows excellent agreement with the test results.
When BCB is used for collapse simulation, the structure collapses vertically, which agrees poorly with the test results, due to the low accuracy of the BCB simulation of the small-deformation From the FE results, the initial moment of collapse is 3.88 s, using the collapse criterion in Section 3.2. A 1 4 -scale model is adopted; thus, the vertical collapse displacement of the collapse criterion is 0.25 m.
Subsequently, collapse simulation was performed using the proposed hybrid method. Figure 11b shows the distribution of ruins, which shows excellent agreement with the test results. The following observations are made from the comparison and analysis of the simulation methods:

•
For FE models, a convergence problem occurs during the collapse simulation; therefore, the distribution of ruins cannot be easily obtained.

•
The BCB simulation method is not sufficiently accurate for the small-deformation stage. Therefore, the collapse mode and the distribution of ruins differ substantially from the test results.
• The proposed hybrid method achieves the most satisfactory simulation of the structural collapse mode and the distribution of ruins. When BCB is used for collapse simulation, the structure collapses vertically, which agrees poorly with the test results, due to the low accuracy of the BCB simulation of the small-deformation stage. As shown in Figure 11c, the final positions of the three slabs of the three-story reinforced concrete frame almost entirely overlap.
Compared with the BCB results (Figure 11c), the distribution of the ruins obtained by the proposed hybrid method (Figure 11b) shows better agreement with the test results (Figure 11a). The following observations are made from the comparison and analysis of the simulation methods:

•
For FE models, a convergence problem occurs during the collapse simulation; therefore, the distribution of ruins cannot be easily obtained.

•
The BCB simulation method is not sufficiently accurate for the small-deformation stage. Therefore, the collapse mode and the distribution of ruins differ substantially from the test results.

•
The proposed hybrid method achieves the most satisfactory simulation of the structural collapse mode and the distribution of ruins.

Collapse and Ruins Simulation of a Real-World Library Building
To demonstrate the applicability of the proposed hybrid method to real-world complex structures, a case study of a multi-story reinforced concrete library building is performed. The building is a seven-story frame-shear wall structure, with a total height of 19.8 m and a building area of 15,938 m 2 . The total mass of the library building is 3.17 × 10 7 kg. All reinforcement consists of HRB400 reinforcing bars, whose strength is 400 MPa. The material properties and dimensions of the main structural members are shown in Table 4. Specifically, the library is located on a site class III in GB50011-2010 [44], with an approximate equivalent shear wave velocity of 200 m/s for 30-m soil (VS30). The characteristic period of the site is 0.45 s. The building has an 8-degree seismic design intensity, where the peak ground acceleration is 200 cm/s 2 for a design basis earthquake (DBE) with a return period of 475 years. There is no R-factor, coefficient of displacement or overstrength factor in Chinese seismic design code. However, according to Lu et al. [51] and American Society of Civil Engineers (ASCE) 7-10 [52], the design information of the building is similar to the following seismic design parameters: R factor = 4.5; overstrength factor = 2.5 and deflection amplification factor = 4. Note that the sizes of walls mean X-direction/Y-direction × Z-direction.
Three element types are used in this model: (a) fiber beam elements for beams and columns, (b) multilayer shell elements for shear walls, and (c) membrane elements for floor slabs. Details of these element models have been reported in the published work [35]. Figure 12a illustrates the three-dimensional FE model. The first mode of the model is a translational mode, with a fundamental period (elastic period) T 1 = 0.37 s.
The widely used El-Centro ground motion record, the same as that in Section 4.1, is adopted. To simulate the collapse behavior, the PGAs of the input motions along the X-and Y-directions are scaled to 4000 cm/s 2 . The nonlinear time history analysis results show that the maximum vertical displacement of the structural component reaches 1 m at 1.76 s, and the corresponding maximum lateral displacement is 0.97 m. The corresponding deformed shape of the structure is shown in Figure 12b. The colored contours in Figure 12 represent the longitudinal reinforcement yielding in the elements. The amount of shear walls on the third story is much lower than on the first and second stories, which results in a sudden change of stiffness. Therefore, the third story is the weak story of the building. The geometric and material information of the FE model, and its displacement and velocity results, are mapped into Blender using the mapping methods discussed in Section 3. After pre-processing with BCB, a collapse simulation is performed using the Bullet physics engine. The collapse process of the library building is then determined, as shown in Figure 13.
The damage first occurs at the soft story of the structure, as shown in Figure 12b. The structural deformation then increases substantially. Excessive deformation finally causes the collapse of the entire structure. The distribution of ruins is shown in Figure 14. The geometric and material information of the FE model, and its displacement and velocity results, are mapped into Blender using the mapping methods discussed in Section 3. After pre-processing with BCB, a collapse simulation is performed using the Bullet physics engine. The collapse process of the library building is then determined, as shown in Figure 13. The geometric and material information of the FE model, and its displacement and velocity results, are mapped into Blender using the mapping methods discussed in Section 3. After pre-processing with BCB, a collapse simulation is performed using the Bullet physics engine. The collapse process of the library building is then determined, as shown in Figure 13.
The damage first occurs at the soft story of the structure, as shown in Figure 12b. The structural deformation then increases substantially. Excessive deformation finally causes the collapse of the entire structure. The distribution of ruins is shown in Figure 14.

Conclusions
In this work, a hybrid framework for simulating building collapse and ruin scenarios using an FE method and physics engine is proposed, and the corresponding program FEM2RBD is developed. Shaking table tests of a three-story reinforced concrete frame and a real-world complex library building are performed to simulate the ruin scenarios. The following conclusions can be drawn:

•
In the proposed hybrid framework, the FE method simulates structural behaviors during the small-deformation stage, and the physics engine simulates structural behaviors during the large-deformation stage. The proposed framework efficiently combines the advantages of the FE method and the physics engine.
• Using a shaking table test of a three-story reinforced concrete frame, the proposed hybrid simulation method is demonstrated to be more accurate than an approach based on the physics engine alone. The case study of a real-world complex library building shows the high-fidelity of the collapse simulation.
• The collision of structural components and the distribution of ruins after the collapse are considered in the proposed hybrid method. The proposed framework has great significance for simulating building collapse and ruin scenarios for post-earthquake rescue training.
Note that the proposed framework is also applicable to other FE software and physics engines, besides the MSC.Marc software and the Bullet physics engine that were used in this work.
Author Contributions: Conceptualization, X.L. and Z.Z.; data curation, X.L. and Z.Z.; formal analysis, Z.Z.; funding acquisition, X.L.; investigation, Y.T. and Z.Z.; methodology, Y.T. and Z.Z.; project administration, X.L.; The damage first occurs at the soft story of the structure, as shown in Figure 12b. The structural deformation then increases substantially. Excessive deformation finally causes the collapse of the entire structure. The distribution of ruins is shown in Figure 14.

Conclusions
In this work, a hybrid framework for simulating building collapse and ruin scenarios using an FE method and physics engine is proposed, and the corresponding program FEM2RBD is developed. Shaking table tests of a three-story reinforced concrete frame and a real-world complex library building are performed to simulate the ruin scenarios. The following conclusions can be drawn:

•
In the proposed hybrid framework, the FE method simulates structural behaviors during the small-deformation stage, and the physics engine simulates structural behaviors during the large-deformation stage. The proposed framework efficiently combines the advantages of the FE method and the physics engine.
• Using a shaking table test of a three-story reinforced concrete frame, the proposed hybrid simulation method is demonstrated to be more accurate than an approach based on the physics engine alone. The case study of a real-world complex library building shows the high-fidelity of the collapse simulation.
• The collision of structural components and the distribution of ruins after the collapse are considered in the proposed hybrid method. The proposed framework has great significance for simulating building collapse and ruin scenarios for post-earthquake rescue training.
Note that the proposed framework is also applicable to other FE software and physics engines, besides the MSC.Marc software and the Bullet physics engine that were used in this work.

Conclusions
In this work, a hybrid framework for simulating building collapse and ruin scenarios using an FE method and physics engine is proposed, and the corresponding program FEM2RBD is developed. Shaking table tests of a three-story reinforced concrete frame and a real-world complex library building are performed to simulate the ruin scenarios. The following conclusions can be drawn:

•
In the proposed hybrid framework, the FE method simulates structural behaviors during the small-deformation stage, and the physics engine simulates structural behaviors during the large-deformation stage. The proposed framework efficiently combines the advantages of the FE method and the physics engine. • Using a shaking table test of a three-story reinforced concrete frame, the proposed hybrid simulation method is demonstrated to be more accurate than an approach based on the physics engine alone. The case study of a real-world complex library building shows the high-fidelity of the collapse simulation.

•
The collision of structural components and the distribution of ruins after the collapse are considered in the proposed hybrid method. The proposed framework has great significance for simulating building collapse and ruin scenarios for post-earthquake rescue training.
Note that the proposed framework is also applicable to other FE software and physics engines, besides the MSC.Marc software and the Bullet physics engine that were used in this work.