Research Article  Open Access
Theoretic Framework and Finite Element Implementation on Progressive Collapse Simulation of Masonry Arch Bridge
Abstract
The capacity that computer can solve more complex design problem is gradually increased. Progressive collapse simulation of masonry arch bridge needs a breakthrough in the current development limitations and then becomes more accurate and integrated. This paper proposes a theoretic framework and finite element implementation on progressive collapse simulation of masonry arch bridge. It is intended to develop a new large deformation element in OpenSees, which can be used for analyzing the collapse process of masonry arch bridge. A mathematical method for large deformation element is put forward by large deformation element. The feature model for bridge structure allows families of bridge components to be specified using constraints on geometry and topology. Geometric constraints are established in bridge components by feature dependence graph in the feature model for bridge. A bridge collapse simulation software system was developed according to such combined technologies. Results from our implementation show that the method can help to simulate the progressive collapse process of masonry arch bridge.
1. Introduction
The masonry arch bridges, a conventional bridge type, are widely used throughout the world because of their low cost and beautiful appearance. However, there were a lot of arch bridges collapsed due to deficiencies in design, detailing, construction, maintenance, use of materials, and inadequate consideration of external events. The first four deficiencies represent integral roles in the building of a bridge [1]. Furthermore, the deficiencies in construction are concerned more because of complicated construction operation needs and relative collapsed cases.
A masonry arch bridge Dixituojiang Bridge under construction over the Dixituojiang River in Hunan, China, suddenly collapsed on the 13st of August 2007 and there were many victims in this disaster [2]. The report about the collapse of the Dixituojiang Bridge clarified improper installation and inadequate temporary structure to support the bridge components were the main reasons. Stresses redistributing in the remaining structures are triggered by the local deficiencies because of improper installation and inadequate temporary structure and then lead to the failure of the adjacent parts. As the chain destruction carries on, it will eventually lead to the loss of bearing capacity of the masonry arch bridge and then result in partial collapse or even a progressive collapse.
American Society of Civil Engineer (ASCE) defines progressive collapse as “the spread of an initial local failure from element to element, eventually resulting in the collapse of an entire structure or a disproportionately large part of it” [3]. Taciroglu et al. developed an arbitrary Lagrangian Eulerian method for regressing solid domains and interface tracking [4]. Marjanishvili [5] presented four successively more sophisticated analysis procedures for evaluating the progressive collapse hazard, linearelastic static; nonlinear static; linearelastic dynamic; and nonlinear dynamic, and discussed the advantages and disadvantages of each method. Szyniszewski [6] from Florida University presented physics based collapse simulations of moment resisting steel frame buildings, with an emphasis on the development of energy flow relationships. Based on the background of Schoharie Creek Bridge collapse, Lebeau and WadiaFascetti [7] offered a detailed review of the events that resulted in this bridge failure through the use of fault tree analysis. An accelerated Newton algorithm based on Krylov subspaces was applied to solving nonlinear equations of structural equilibrium by Scott and Fenves [8].
Lau and Wibowo [9, 10] carried out nonlinear analyses of reinforced concrete bridges by the Applied Element Method (AEM), studying the progressive collapse phenomenon of structures during earthquakes. Based on a sufficient investigation for the reason of the Mississippi River steel truss bridge collapsing on August 1, 2007, AstanehAsl et al. [11, 12] established three analytical models to simulate the collapse process of this bridge. Accurate stress analysis on steel box girder of long span suspension bridges based on multiscale submodeling method in order to decrease computational complexity of complicated bridges was proposed [13].
In China, there are many aging masonry arch bridges as well and they need prompt inspection, reinforcement, and maintenance. Many studies focused on the progressive collapse and clarified that local deficiencies may trigger collapse [2].
However, there are few researches on how to simulate the progressive collapse process of the masonry arch bridges. There are also few researches about how to use large deformation elements to simulate geometric nonlinearity properties during bridge collapse. In this research, it is intended to develop a new large deformation element, which can be used for analyzing the collapse process of masonry arch bridge. A mathematical method for large deformation element was put forward by large deformation element. Meanwhile, computergraphics technology, feature modeling, and finite element method were combined to develop Bridge Collapse Software System (BCSS) for the progressive collapse process simulation.
2. The Mathematical Methods for Dynamic Calculation for the Collapse Process of Masonry Arch Bridge
2.1. The Large Deformation Equation
2.1.1. Strain under Large Deformation
Object in Cartesian Coordinates will continuously change its configuration under an external force, as shown in Figure 1.
At time , the coordinate of is ; the coordinate of is .
At time , the coordinate of is ; the coordinate of is , .
Regard the change of object configuration as a mathematic change from to . For a certain time , this transformation can be defined byAccording to the deformation continuity requirements, the change must be in onetoone correspondence, namely, the transformation is singlevalued continuous. Meanwhile this transformation should have unique inverse transformation; namely, it has singlevalued continuous inverse transformation, which can be shown as follows:In the view of the definition above, and can be written as Using the formula above, the distance of and at time and time isAfter deformation, the distance of the two points changed intoEquation (5) can be rewritten aswhere is the GreenLagrange strain tensor, which is defined by coordinate before, and is the Almansi strain tensor, which is defined by coordinate after deformation.
To obtain the relation between strain and displacement, the displacement field is given bywhere is the displacement from configuration before deformation to configuration after deformation.
We then haveUsing (6) and (8) yields the relation between strain and displacement:When the displacement is very small, the quadratic term of displacement derivative can be ignored. This equation reduces toSince , we further obtainTherefore, we conclude that the necessary and sufficient condition of an object doing rigid motion is that the components of the Green strain tensor and the Almansi strain tensor to be zero.
2.1.2. Stress under Large Deformation
In the large deformation problem, it is by cutting off a microelement from the deformed object that the equilibrium equation and its equivalent virtual work principle to be established. Accordingly we should cut off an element from the deformed object to define the stress tensor, that is, Euler stress tensor , as shown in Figure 2.
If the Green strain tensor is defined by coordinates before deformation, then its corresponding stress tensor before deformation is required to be defined. Stress on the plane after deformation is ; and assume virtual stress on the plane before deformation is .
In order to maintain the consistency of mathematics, there are usually two kinds of provisions for expressing the relationship between and , the schematic diagram is shown in Figure 3.
(1) Lagrange Provision. The internal force on area element before deformation is equal to the internal force after deformation and reads
(2) Kirchhoff Provision. The transformation of internal force on area element before and after deformation is consistent with coordinate transformation and readsBecause is the stress component of the configuration after deformation, definewhere is the direction cosine on area element before deformation; is the Lagrange strain tensor; is the Kirchhoff strain tensor; the left superscript 0 means that it is measured on configuration before deformation; and the left subscript means that it is measured on configuration after deformation.
2.1.3. Material Constitutive Relation under Large Deformation
As for small deformation and linear elasticity problem under isothermal and adiabatic conditions, the relationship between stress and strain can be described by the following three kinds of equivalent methods:For the small strain elastic problem, the stress and strain is a corresponding relation. In the condition of small strain, the relationship is linear; in the condition of large strain, the relation is nonlinear. That is, Using superelasticity representing such material properties in continuous medium mechanics yieldswhere is Kirchhoff stress tensor, namely, the engineering stress; is Green strain tensor, namely, the engineering strain; is constant and elastic constitutive tensor; is the strain energy function.
If the materials are plastic, thenwhere is tangent constitutive tensor on time , referring to configuration on time :
2.1.4. The Formation of the Control Equation of Geometric Nonlinearity and Virtual Displacement Principle
The coordinates and displacement of each point on the object at time 0: , .
The coordinates and displacement of each point on the object at time : , .
The coordinates and displacement of each point on the object at time : , .
The virtual displacement principle corresponding to the equilibrium condition of the configuration at time can be written aswhere is the variation of the real displacement; is the variation of the strain; is the surface load measured in the real configuration; is volume load measured in the real configuration; is the mass density of the real configuration.
For total Lagrange Formulation, we haveNotes: because OpenSees (the Open System for Earthquake Engineering simulation, is an objectoriented, open source software framework) uses the total Lagrange Formulation, the updated Lagrange formulation will not be discussed in detail.
Disperse the solution domain by using isoparametric element; then the coordinates and displacement of each element can be obtained by nodal value interpolation:Then analyze an element by using the total Lagrange Formulation, so that we finally obtainwhere is the displacement vector of the node; is the transformation matrix of the linear strain displacement; is the transformation matrix of the nonlinear strain displacement; is the material constitutive matrix; is the Kirchhoff strain matrix.
2.2. Geometric Nonlinear Finite Element Theory and Its Realization
2.2.1. Tentative Simulation of Small Strain Elements
Considering the computational cost, bbarBrick Element in OpenSees was applied to simulate the progressive collapse of bridge based on reduced integration. bbarBrick does not work because it cannot be applied under large deflection or large rotation, such as “snap through.” The process of bridge progressive collapse is the geometrically nonlinear problems with large rotation. On Figures 4(a) and 4(b), an element in the abutment at the beginning of collapsing and another element in deck at the end of collapsing have large rotation and “snap through,” which leads to the calculation results being distorted. We conclude that using small strain element cannot simulate geometrical nonlinear problems in the process of bridge collapse.
(a) A snap through element in abutment
(b) A snap through element in deck
2.2.2. Collapse Simulation of Large Strain Elements
By using the mathematical methods mentioned in Section 2.1 and doing a secondary development of TLFD8nBrick element in OpenSees, a finite element that can calculate large deformation problems is established, as shown in Figure 5.
(a) Initial failure
(b) No snap through element during collapsing
3. Feature Modeling Method of Bridge Components
The feature modeling provides a rapid method for parameterization of bridge models, through which a parametric model for bridges can be set up; and influencing factors analysis of bridge collapse can be carried out easily. The Object Oriented Program (OOP) was used to represent bridge elements.
3.1. Sections Organized by Feature Parameters
An arch bridge can be regarded as a collection of various types of components, and every component can be stretched by one or more bridge sections. A section can be defined by changing the characteristic parameters which are familiar to engineers, such as the roof thickness of boxtype section, web thickness, and floor thickness and flange length. Engineers can input few relevant parameters to complete fast the modeling of typical section, as shown in Figure 6. And so the map between feature parameters and geometry coordinates has been created.
There are many component section types in actual project. Each section has its unique mapping relationship from feature parameters to threedimensional coordinate. Moreover, each section class has public functions or the methods, such as area calculation, reasonable judgment, section translation, section rotation, and eccentric conversion. And so a base class (parent class) supporting public functions or methods are established. Each section class (child class) of different types may call the methods and functions of the base class conveniently.
3.2. Geometric Graphics Description Based on Binary SpacePartitioning (BSP) Trees
A section is encircled by several directed edges, according to which its binary spacepartitioning (BSP) tree can be constructed. The commonly used method to construct a BSP tree for a section is to create a series of nodes, making each node represent a partition line that contains an edge of the section. Other polygon edges are separated by the partition line in the position of the nodes. The section geometric graphics definition is carried out by the Boolean operation on the section. By using difference operation, one of the Boolean operations, a section with hole can be defined conveniently, which provides a convenient and visual representation approach for the box section, a wildly used section in bridge engineering.
Figure 7 shows the process of outside region subtract inside region(s), the purpose of Boolean subtract operation is to make a description and definition for the general box section and hollow section.
Let and be two polygons in 2D space. The difference of polygons and , means is subtracted from , which is denoted as . And define as taking negation of . Thus the calculation formula can be written as follows:Take the quadrilateral section EFGD in Figure 8 as an example. Now it is needed to add a quadrilateral hole efgd in its internal. The BSP tree of this section after digging a hole is as shown in Figure 8.
(a)
(b)
As is shown in Figure 8, the operation of digging a quadrilateral hole efgd in section EFGD is, in fact, taking as a new 2D spatial region and then partitioning this region using quadrilateral section efgd. Accordingly, the BSP tree is to add more new nodes on the basis of the original one. According to formula , where is the section EFGD and is the hole efgd, thus is to do negation of section efgd, that is, changing directions of the directed edges. As Figure 8 shows, the direction of the directed edges of section efgd is already opposite to the one of section EFGD. On the basis of the righthand rule, we can get the section with hole , and from the BSP tree of this section we can see , , , and are on the negative side, which is consistent with Figure 8(a).
Examples of box section by “1 box and 3 rooms” and box section by “2 box and 1 room” in graphic platform are shown in Figure 9.
(a) “1box and 3room” girder section
(b) “2box and 1room” girder section
3.3. Forming a 3D Model by Projection
According to the modeling process, the cross wall will be built after the formation of the main arch ring. Since the cross wall is attached to the main arch, the bottom surface of each cross wall and the top surface of the main arch ring must be coplanar. And it is based on these surfaces that every 3D model of bridges including the cross wall can be formed. Make projections perpendicular to the horizontal plane for these surfaces to a specified height to the projection plane and then establish the cross wall model through connecting the top surfaces and the projection planes. To make projection of an objective plane onto assigned plane is essential to get the projection of the vertices of the objective plane onto the assigned plane. That is to say, when you get all of the projection of the vertices, you get the projection of the objective plane.
In order to obtain the projection of onto the plane , define the vector of plane as , where , as shown in Figure 10.
In the picture, is the projection of onto plane . The line connecting and its projection , by definition, are parallel to the normal vector of the plane .
Let the length of be . We obtain
So we haveAs can be seen from the formulas above, the algorithms needed to add in program designing mainly include vector unitization, the distance of two points, product of two vectors, and direction vector of plane. The direction vector of plane can be obtained by product of two different direction vectors in plane. That is, selecting two arbitrary edges in plane, by which we have two vectors to do their product, and the result is the direction vector of the plane.
After getting the projection of the objective plane onto assigned plane, the projection plane is formed, and then its 3D Solid Model can be established on the basis of the two planes. As shown in Figure 11, the standard plane is a torus, which is in the XOY plane. And the projection plane is created by three points: , , and . Figure 12 shows the 3D model of box girder and arch by projection.
(a) Box girder model
(b) Arch model
3.4. Organization Form of MasterSlave Component
Arch bridges are a collection of various types of components, and the process of establishing the arch bridge model is also the organization of these components. There are three basic organizations that have been utilized in the design of organization for bridge components: masterslave; separate executive for each component; and symmetric or anonymous treatment of all components. For most bridge components, the first organization is availably operated in the masterslave mode. This type is certainly the easiest to implement and may often be produced by making relatively simple extensions to an organization that includes full multiprogramming capabilities.
Therefore we need to establish the feature dependence graph between master and slave which includes following two aspects.
(1) Constraints on Geometry of All Components. Each component model needs to be organized to connect with each other. Constraints are created between master and slave by position in space, as shown in Figure 13.
(2) Form a Computable 3D Model. For general finite element calculation program, the model entities are required numbering the element and nodes. And the connection between component and component reacts in the node coupling between two components, as shown in Figure 14.
For so many arch bridge components, slave components are attached to one or more master components by position or geometry parameters. The masterslave component diagram can be established as in Figure 15.
In Figure 15, G1, G2 are the highest levels of master component, and G3, G4, G5, G6, and G7 on the second layer exist depending on the G1 and G2, which are slave components of G1 and G2. And G8~G17 on the third layer are slave components of G3~G7. The lines in the graph represent feature dependence graph between masters and slave. A slave component has feature dependence graph with one or more master components.
4. Development of Bridge Collapse Simulation Software System
Bridge collapse simulation software system (BCSS) is developed by C# and OpenGL. C# is a massively popular programming language used by many thousands of software developers all over the world. The OpenGL graphics system is a software interface to graphics hardware. (The GL stands for Graphics Library.) It allows programmers to create interactive programs that produce color images of moving threedimensional objects. With OpenGL, programmers can control computergraphics technology to produce realistic pictures or ones that depart from reality in imaginative ways.
BSCC includes feature modeling module and bridge collapse visualization module, as shown in Figure 16.
5. Application Examples
At 4:00 p.m. on June 13, 2007, the Dixituojiang Bridge under construction collapsed suddenly in Fenghuang County of Hunan Province, China, and 59 people were killed. The full length of the Dixituojiang Bridge is 328.5 m. And it has four holes, each of which is 65 m in span length. The lane width of the bridge is 12 m. The main arch is equal in the cross section, and the type of curve is catenary. The risespan ratio of main arch is 1/5, the archaxis coefficient is 2.514, and the thickness is 1.35 m. There are seven Associate Arches on the top of each main arch. The span of Associate Arch is 5 m, and the risespan ratio is 1/5. The thickness of the Associate Arch is 0.4 m. It is also the equal section component, and the type of curve is arc curve. The thickness of the cross wall is 1.1 m. The elevation of Dixituojiang Bridge is shown in Figure 17.
The main arch is composed of No. 20 pebble and No. 60 rubble. The Associate Arch is composed by No. 12.5 mortar and No. 30 dressed stone. The No. 12.5 mortar and No. 30 stone are applied to bricklaying of the cross wall.
5.1. Feature Modeling for Masonry Arch Bridge
According to the characteristics of masonry arch bridge components, this feature modeling is divided into three modules: the module of main arch, the module of cross wall and spandrel arch, and the module of abutment and foundation. Take Dixituojiang Bridge as an example, the design parameters of each module are as shown in Table 1.

With the material parameters above, a 3D FEM model of Dixituojiang Bridge is established by feature modeling method, as shown in Figure 18.
5.2. The Collapse Process and Analysis of the Results
By the cause analysis of the accident: the bridge piers adopting spread foundation stand out in the slightly weathered rock. When providing vertical stress, the foundation can only partially constrain the displacement in the horizontal direction of bridge pier. When the horizontal thrust exceeds the allowable value, the pier may have lateral displacement or even overturn, which is different from fixed constraint. Therefore, in order to ensure the accuracy of the analysis, the zerolength contact element was set between nodes in foundation and in the bottom of the pier.
Due to the geometric nonlinear element analysis, and setting zerolength contact element to simulate the boundary of the pier and the foundation, and adopting the birth and death element technology in the process of analysis, all of which lead to the analysis having extremely strong nonlinearity.
The meshes of the original model are transversely coarsened, the final model is divided into 30,356 8node hexahedral elements, and set 24 zerolength contact elements on the right side of the contact surface between the pier and ground, totally has 51086 nodes. The model sets up 950 calculation steps, and each step is 0.008 s. Before the first arch collapsed, 0.15 cm displacement is applied to the left abutment for each iteration step.
Considering the low quality of No. 1 main arch in construction, the ultimate tensile strain and ultimate compressive strain of No. 1 main arch should be reduced properly, where the former decreases 30% and the latter decreases about 20%. Meanwhile, set the rule that elements fail when or ; for the other main arches, elements fail when or . Then we can get the progressive collapse process similar with the actual situation, as shown in Figure 19.
(a) Initial failure is occurred on s
(b) The first main arch collapses on s
(c) The second main arch collapses on s
(d) The third main arch collapses on s
(e) The forth main arch collapses on s
(f) The whole bridge collapses on s
The simulation of collapse process suggests that the continuous arch effect caused by the failure of No. 1 main arch leads to the collapse of No. 2 and No. 3 main arches in succession, which is consistent with the actual collapse process. From Figures 19(a)–19(f), all of the bridge piers fall down towards the direction of No. 0 abutment, which is also consistent with the actual situation.
About s, namely, the time when the displacement of the left abutment is nearly 8~8.5 cm, the floor system and spandrel arch above the left skewback firstly occur tensile failure elements, and then the elements in the connecting part of the left skewback and the left abutment are killed, which suggests the structure internal force is redistributed for the first time, and then it is quickly reflected through killing the elements of the right skewback of the No. 1 main arch. Meanwhile there are also elements of other parts such as spandrel arch and floor system being killed. About s, most elements of the right skewback of No. 1 main arch are killed, and then No. 1 main arch begins to fall down. Then the rest main arches collapse in succession and eventually result in the whole bridge collapsed.
6. Conclusions
This paper proposes a theoretic framework and finite element implementation on progressive collapse simulation of masonry arch bridge. A new large deformation element in OpenSees was developed, which can be used for analyzing the collapse process of masonry arch bridge. A mathematical method for large deformation element is put forward based on large deformation element. The feature model which allows families of bridge components to be specified using constraints on geometry and topology facilitate the modeling of complicated bridge. Geometric constraints are established in bridge components by feature dependence graph in the feature model for bridge. Results from our implementation show that the method can help to simulate the progressive collapse process of masonry arch bridge.
Conflict of Interests
The authors declared that they have no conflict of interests to this work.
Acknowledgments
This research is financially supported by the Zhejiang Provincial Natural Science Foundation of China (Grant no. LY13E080014). This support is gratefully acknowledged. Any opinions, findings, and conclusions or recommendations expressed in this material are those of authors.
References
 K. Wardhana and F. C. Hadipriono, “Analysis of recent bridge failures in the United States,” Journal of Performance of Constructed Facilities, vol. 17, no. 3, pp. 144–150, 2003. View at: Publisher Site  Google Scholar
 Z. Xu, X. Lu, H. Guan, X. Lu, and A. Ren, “Progressivecollapse simulation and critical region identification of a stone arch bridge,” Journal of Performance of Constructed Facilities, vol. 27, no. 1, pp. 43–52, 2013. View at: Publisher Site  Google Scholar
 ASCE, “Minimum design loads for buildings and other structures,” ASCE 705, 2005. View at: Google Scholar
 E. Taciroglu, A. Acharya, A. Namazifard, and I. D. Parsons, “Arbitrary LagrangianEulerian methods for analysis of regressing solid domains and interface tracking,” Computers & Structures, vol. 87, no. 56, pp. 355–367, 2009. View at: Publisher Site  Google Scholar
 S. M. Marjanishvili, “Progressive analysis procedure for progressive collapse,” Journal of Performance of Constructed Facilities, vol. 18, no. 2, pp. 79–85, 2004. View at: Publisher Site  Google Scholar
 S. Szyniszewski, “Dynamic energy based method for progressive collapse analysis,” in Proceedings of the Structures Congress: Don't Mess with Structural Engineers: Expanding Our Role, pp. 1–10, Austin, Tex, USA, AprilMay 2009. View at: Publisher Site  Google Scholar
 K. H. Lebeau and S. J. WadiaFascetti, “Fault tree analysis of schoharie creek bridge collapse,” Journal of Performance of Constructed Facilities, vol. 21, no. 4, pp. 320–326, 2007. View at: Publisher Site  Google Scholar
 M. H. Scott and G. L. Fenves, “Krylov subspace accelerated newton algorithm: application to dynamic progressive collapse simulation of frames,” Journal of Structural Engineering, vol. 136, no. 5, pp. 473–480, 2010. View at: Publisher Site  Google Scholar
 H. Wibowo, S. S. Reshotkina, and D. T. Lau, “Modeling progressive collapse of RC bridges during earthquakes,” in Proceedings of the CSCE Annual General Conference, GC1761 GC17610, May 2009. View at: Google Scholar
 A. Colombo and P. Negro, “Seismic progressive collapse analysis of reinforced concrete bridges by applied element method,” in Earth and Space 2010: Engineering, Science, Construction, and Operations in Challenging Environments, vol. 27, pp. 3019–3026, American Society of Civil Engineers, Reston, Va, USA, 2010. View at: Google Scholar
 A. AstanehAsl, “Progressive collapse of steel Truss bridges, the case of I35W collapse,” in Proceedings of the 7th International Conference on Steel Bridges, pp. 1–10, Guimarăes, Portugal, 2008. View at: Google Scholar
 S. Hao, “A note of the I35W bridge collapse,” Journal of Bridge Engineering, vol. 15, no. 5, pp. 608–614, 2010. View at: Google Scholar
 H. Wang, A. Li, R. Hu, and J. Li, “Accurate stress analysis on steel box girder of long span suspension bridges based on multiscale submodeling method,” Advances in Structural Engineering, vol. 13, no. 4, pp. 727–740, 2010. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2015 Weibing Peng et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.