# Theoretic Framework and Finite Element Implementation on Progressive Collapse Simulation of Masonry Arch Bridge

Authors: Weibing Peng; Ruodan Pan; Fei Dai
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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.


---

## Body

## 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, linear-elastic static; nonlinear static; linear-elastic 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 Wadia-Fascetti [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, Astaneh-Asl 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, computer-graphics 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 Figure1.Figure 1
Geometrical nonlinear configuration.At timet0, the coordinate of P is xi0; the coordinate of Q is xi0+dxi0.At timet, the coordinate of P is xit; the coordinate of Q is xit+dxit, i=1,2,3.Regard the change of object configuration as a mathematic change fromxi0 to xit. For a certain time t, this transformation can be defined by(1)xit=xitx10,x20,x30.According to the deformation continuity requirements, the change must be in one-to-one correspondence, namely, the transformation is single-valued continuous. Meanwhile this transformation should have unique inverse transformation; namely, it has single-valued continuous inverse transformation, which can be shown as follows:(2)xi0=xi0x1t,x2t,x3t.In the view of the definition above, dxi0 and dxit can be written as (3)dxi0=∂xi0∂xjtdxjt=xi,jt0dxjt,dxit=∂xit∂xj0dxj0=xi,j0idxj0.Using the formula above, the distance of P and Q at time t0 and time t is(4)d0s2=dxi0dxi0=xi,mt0xi,nt0dxmtdxnt,dts2=dxitdxit=xi,m0txi,n0tdxm0dxn0.After deformation, the distance of the two points changed into(5)dts2-d0s2=xk,i0txk,j0t-δijdxi0dxj0=2εij0tdxi0dxj0,dts2-d0s2=δij-xk,it0xk,jt0dxitdxjt=2εijttdxitdxjt.Equation (5) can be rewritten as(6)εij0t=12xk,i0txk,j0t-δij,εijtt=12δij-xk,it0xk,jt0,where εij0t is the Green-Lagrange strain tensor, which is defined by coordinate before, and εijtt is the Almansi strain tensor, which is defined by coordinate after deformation.To obtain the relation between strain and displacement, the displacement fielduit is given by(7)uit=xit-xi0,where uit is the displacement from configuration before deformation to configuration after deformation.We then have(8)xij0t=δij+ui,j0t,εij0t=δij-ui,jtt.Using (6) and (8) yields the relation between strain and displacement:(9)εijtt=12ui,jtt+uj,itt-uk,ittuk,jtt,εij0t=12ui,j0t+uj,i0t+uk,i0tuk,j0t.When the displacement is very small, the quadratic term of displacement derivative can be ignored. This equation reduces to(10)εij0t=εijtt=εij.Since dts2-d0s2=0, we further obtain(11)εij0t=εijtt=0.Therefore, 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τijt, as shown in Figure 2.Figure 2
Geometrical nonlinear strains.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 planePtQtRtSt after deformation is dTt/dSt; and assume virtual stress on the plane P0Q0R0S0 before deformation is dT0/dS0.In order to maintain the consistency of mathematics, there are usually two kinds of provisions for expressing the relationship betweendT0 and dTt, the schematic diagram is shown in Figure 3.Figure 3
Geometrical nonlinear stresses.(1) Lagrange Provision. The internal force on area element before deformation is equal to the internal force after deformation and reads(12)d0TiL=dtTi.(2) Kirchhoff Provision. The transformation of internal force on area element before and after deformation is consistent with coordinate transformation and reads(13)d0TiK=xi,jt0dtTj.Because τijt is the stress component of the configuration after deformation, define(14)dtTi=τijtνjtdS,Lagrange  provision:d0TiL=Tji0tνj0dS0=dtTi,Kirchhoff  provision:d0TiK=Sji0tνj0dS0=xi,jt0dtTj,where νjt is the direction cosine on area element before deformation; Tjit0 is the Lagrange strain tensor; Sjit0 is the Kirchhoff strain tensor; the left superscript 0 means that it is measured on configuration before deformation; and the left subscript t 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:(15)σij=Dijklεkl,σij=∂W∂εij,W=12Dijklεijεkl,∂σij∂t=Dijkl∂εkl∂t.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, (16)Sij0t=Dijkl0tεkl0t.Using superelasticity representing such material properties in continuous medium mechanics yields(17)Sij0t=ρ0∂W2t∂0tεij∂0tεij,where Sij0t is Kirchhoff stress tensor, namely, the engineering stress; εij0t is Green strain tensor, namely, the engineering strain; Dijkl0t is constant and elastic constitutive tensor; Wt is the strain energy function.If the materials are plastic, then(18)Sij0t=Dijkl0εkl0t,where Dijkl0 is tangent constitutive tensor on time t, referring to configuration on time t0: (19)Dijkl0=Deijkl-Dpijkl0.


#### 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:xi0, ui0.The coordinates and displacement of each point on the object at timet: xit, uit.The coordinates and displacement of each point on the object at timet+Δt: xit+Δt, uit+Δt.The virtual displacement principle corresponding to the equilibrium condition of the configuration at timet+Δt can be written as(20)∫Vt+Δtτijt+Δtδt+Δteijt+ΔtdV=∫Vt+Δttkt+Δtt+Δtδukt+ΔtdS+∫Vt+Δtρt+Δtfkt+Δtt+Δtδukt+ΔtdV,where δuk is the variation of the real displacement; δt+Δteijt+Δt is the variation of the strain; tkt+Δtt+Δt is the surface load measured in the real configuration; ft+Δtt+Δt is volume load measured in the real configuration; ρt+Δt is the mass density of the real configuration.For total Lagrange Formulation, we have(21)∫0VSij0t+Δtδεij00t+ΔtdV=∫0Stk0t+Δtδuk0dS+∫0Vρ0fkδt+Δtuk0dV.Notes: because OpenSees (the Open System for Earthquake Engineering simulation, is an object-oriented, 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 coordinatesxi0 and displacement uit of each element can be obtained by nodal value interpolation:(22)xi0=∑k=1nNkxik0,xit=∑k=1nNkxikt,xit+Δt=∑k=1nNkxikt+Δt,uit=∑k=1nNkuikt,uit+Δt=∑k=1nNkuikt+Δt.Then analyze an element by using the total Lagrange Formulation, so that we finally obtain(23)KL0t+KNL0tu=Qt+Δt-F0t,KL0t=∑∫0VBLT0tD0BL00tdV,KNL0t=∑∫0VBNLT0tS0tBNL00tdV,F0t=∑∫0VB0tLTS^00tdV,where u is the displacement vector of the node; BL0t is the transformation matrix of the linear strain displacement; BNL0t is the transformation matrix of the nonlinear strain displacement; D0 is the material constitutive matrix; S0t 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 Figures4(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.Figure 4
Collapse simulations by small strain elements in BCSS.
(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 Section2.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.Figure 5
Collapse simulations by large strain elements in BCSS.
(a)
Initial failure
(b)
No snap through element during collapsing


## 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 Figure1.Figure 1
Geometrical nonlinear configuration.At timet0, the coordinate of P is xi0; the coordinate of Q is xi0+dxi0.At timet, the coordinate of P is xit; the coordinate of Q is xit+dxit, i=1,2,3.Regard the change of object configuration as a mathematic change fromxi0 to xit. For a certain time t, this transformation can be defined by(1)xit=xitx10,x20,x30.According to the deformation continuity requirements, the change must be in one-to-one correspondence, namely, the transformation is single-valued continuous. Meanwhile this transformation should have unique inverse transformation; namely, it has single-valued continuous inverse transformation, which can be shown as follows:(2)xi0=xi0x1t,x2t,x3t.In the view of the definition above, dxi0 and dxit can be written as (3)dxi0=∂xi0∂xjtdxjt=xi,jt0dxjt,dxit=∂xit∂xj0dxj0=xi,j0idxj0.Using the formula above, the distance of P and Q at time t0 and time t is(4)d0s2=dxi0dxi0=xi,mt0xi,nt0dxmtdxnt,dts2=dxitdxit=xi,m0txi,n0tdxm0dxn0.After deformation, the distance of the two points changed into(5)dts2-d0s2=xk,i0txk,j0t-δijdxi0dxj0=2εij0tdxi0dxj0,dts2-d0s2=δij-xk,it0xk,jt0dxitdxjt=2εijttdxitdxjt.Equation (5) can be rewritten as(6)εij0t=12xk,i0txk,j0t-δij,εijtt=12δij-xk,it0xk,jt0,where εij0t is the Green-Lagrange strain tensor, which is defined by coordinate before, and εijtt is the Almansi strain tensor, which is defined by coordinate after deformation.To obtain the relation between strain and displacement, the displacement fielduit is given by(7)uit=xit-xi0,where uit is the displacement from configuration before deformation to configuration after deformation.We then have(8)xij0t=δij+ui,j0t,εij0t=δij-ui,jtt.Using (6) and (8) yields the relation between strain and displacement:(9)εijtt=12ui,jtt+uj,itt-uk,ittuk,jtt,εij0t=12ui,j0t+uj,i0t+uk,i0tuk,j0t.When the displacement is very small, the quadratic term of displacement derivative can be ignored. This equation reduces to(10)εij0t=εijtt=εij.Since dts2-d0s2=0, we further obtain(11)εij0t=εijtt=0.Therefore, 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τijt, as shown in Figure 2.Figure 2
Geometrical nonlinear strains.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 planePtQtRtSt after deformation is dTt/dSt; and assume virtual stress on the plane P0Q0R0S0 before deformation is dT0/dS0.In order to maintain the consistency of mathematics, there are usually two kinds of provisions for expressing the relationship betweendT0 and dTt, the schematic diagram is shown in Figure 3.Figure 3
Geometrical nonlinear stresses.(1) Lagrange Provision. The internal force on area element before deformation is equal to the internal force after deformation and reads(12)d0TiL=dtTi.(2) Kirchhoff Provision. The transformation of internal force on area element before and after deformation is consistent with coordinate transformation and reads(13)d0TiK=xi,jt0dtTj.Because τijt is the stress component of the configuration after deformation, define(14)dtTi=τijtνjtdS,Lagrange  provision:d0TiL=Tji0tνj0dS0=dtTi,Kirchhoff  provision:d0TiK=Sji0tνj0dS0=xi,jt0dtTj,where νjt is the direction cosine on area element before deformation; Tjit0 is the Lagrange strain tensor; Sjit0 is the Kirchhoff strain tensor; the left superscript 0 means that it is measured on configuration before deformation; and the left subscript t 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:(15)σij=Dijklεkl,σij=∂W∂εij,W=12Dijklεijεkl,∂σij∂t=Dijkl∂εkl∂t.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, (16)Sij0t=Dijkl0tεkl0t.Using superelasticity representing such material properties in continuous medium mechanics yields(17)Sij0t=ρ0∂W2t∂0tεij∂0tεij,where Sij0t is Kirchhoff stress tensor, namely, the engineering stress; εij0t is Green strain tensor, namely, the engineering strain; Dijkl0t is constant and elastic constitutive tensor; Wt is the strain energy function.If the materials are plastic, then(18)Sij0t=Dijkl0εkl0t,where Dijkl0 is tangent constitutive tensor on time t, referring to configuration on time t0: (19)Dijkl0=Deijkl-Dpijkl0.


### 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:xi0, ui0.The coordinates and displacement of each point on the object at timet: xit, uit.The coordinates and displacement of each point on the object at timet+Δt: xit+Δt, uit+Δt.The virtual displacement principle corresponding to the equilibrium condition of the configuration at timet+Δt can be written as(20)∫Vt+Δtτijt+Δtδt+Δteijt+ΔtdV=∫Vt+Δttkt+Δtt+Δtδukt+ΔtdS+∫Vt+Δtρt+Δtfkt+Δtt+Δtδukt+ΔtdV,where δuk is the variation of the real displacement; δt+Δteijt+Δt is the variation of the strain; tkt+Δtt+Δt is the surface load measured in the real configuration; ft+Δtt+Δt is volume load measured in the real configuration; ρt+Δt is the mass density of the real configuration.For total Lagrange Formulation, we have(21)∫0VSij0t+Δtδεij00t+ΔtdV=∫0Stk0t+Δtδuk0dS+∫0Vρ0fkδt+Δtuk0dV.Notes: because OpenSees (the Open System for Earthquake Engineering simulation, is an object-oriented, 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 coordinatesxi0 and displacement uit of each element can be obtained by nodal value interpolation:(22)xi0=∑k=1nNkxik0,xit=∑k=1nNkxikt,xit+Δt=∑k=1nNkxikt+Δt,uit=∑k=1nNkuikt,uit+Δt=∑k=1nNkuikt+Δt.Then analyze an element by using the total Lagrange Formulation, so that we finally obtain(23)KL0t+KNL0tu=Qt+Δt-F0t,KL0t=∑∫0VBLT0tD0BL00tdV,KNL0t=∑∫0VBNLT0tS0tBNL00tdV,F0t=∑∫0VB0tLTS^00tdV,where u is the displacement vector of the node; BL0t is the transformation matrix of the linear strain displacement; BNL0t is the transformation matrix of the nonlinear strain displacement; D0 is the material constitutive matrix; S0t is the Kirchhoff strain matrix.


## 2.1.1. Strain under Large Deformation

Object in Cartesian Coordinates will continuously change its configuration under an external force, as shown in Figure1.Figure 1
Geometrical nonlinear configuration.At timet0, the coordinate of P is xi0; the coordinate of Q is xi0+dxi0.At timet, the coordinate of P is xit; the coordinate of Q is xit+dxit, i=1,2,3.Regard the change of object configuration as a mathematic change fromxi0 to xit. For a certain time t, this transformation can be defined by(1)xit=xitx10,x20,x30.According to the deformation continuity requirements, the change must be in one-to-one correspondence, namely, the transformation is single-valued continuous. Meanwhile this transformation should have unique inverse transformation; namely, it has single-valued continuous inverse transformation, which can be shown as follows:(2)xi0=xi0x1t,x2t,x3t.In the view of the definition above, dxi0 and dxit can be written as (3)dxi0=∂xi0∂xjtdxjt=xi,jt0dxjt,dxit=∂xit∂xj0dxj0=xi,j0idxj0.Using the formula above, the distance of P and Q at time t0 and time t is(4)d0s2=dxi0dxi0=xi,mt0xi,nt0dxmtdxnt,dts2=dxitdxit=xi,m0txi,n0tdxm0dxn0.After deformation, the distance of the two points changed into(5)dts2-d0s2=xk,i0txk,j0t-δijdxi0dxj0=2εij0tdxi0dxj0,dts2-d0s2=δij-xk,it0xk,jt0dxitdxjt=2εijttdxitdxjt.Equation (5) can be rewritten as(6)εij0t=12xk,i0txk,j0t-δij,εijtt=12δij-xk,it0xk,jt0,where εij0t is the Green-Lagrange strain tensor, which is defined by coordinate before, and εijtt is the Almansi strain tensor, which is defined by coordinate after deformation.To obtain the relation between strain and displacement, the displacement fielduit is given by(7)uit=xit-xi0,where uit is the displacement from configuration before deformation to configuration after deformation.We then have(8)xij0t=δij+ui,j0t,εij0t=δij-ui,jtt.Using (6) and (8) yields the relation between strain and displacement:(9)εijtt=12ui,jtt+uj,itt-uk,ittuk,jtt,εij0t=12ui,j0t+uj,i0t+uk,i0tuk,j0t.When the displacement is very small, the quadratic term of displacement derivative can be ignored. This equation reduces to(10)εij0t=εijtt=εij.Since dts2-d0s2=0, we further obtain(11)εij0t=εijtt=0.Therefore, 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τijt, as shown in Figure 2.Figure 2
Geometrical nonlinear strains.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 planePtQtRtSt after deformation is dTt/dSt; and assume virtual stress on the plane P0Q0R0S0 before deformation is dT0/dS0.In order to maintain the consistency of mathematics, there are usually two kinds of provisions for expressing the relationship betweendT0 and dTt, the schematic diagram is shown in Figure 3.Figure 3
Geometrical nonlinear stresses.(1) Lagrange Provision. The internal force on area element before deformation is equal to the internal force after deformation and reads(12)d0TiL=dtTi.(2) Kirchhoff Provision. The transformation of internal force on area element before and after deformation is consistent with coordinate transformation and reads(13)d0TiK=xi,jt0dtTj.Because τijt is the stress component of the configuration after deformation, define(14)dtTi=τijtνjtdS,Lagrange  provision:d0TiL=Tji0tνj0dS0=dtTi,Kirchhoff  provision:d0TiK=Sji0tνj0dS0=xi,jt0dtTj,where νjt is the direction cosine on area element before deformation; Tjit0 is the Lagrange strain tensor; Sjit0 is the Kirchhoff strain tensor; the left superscript 0 means that it is measured on configuration before deformation; and the left subscript t 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:(15)σij=Dijklεkl,σij=∂W∂εij,W=12Dijklεijεkl,∂σij∂t=Dijkl∂εkl∂t.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, (16)Sij0t=Dijkl0tεkl0t.Using superelasticity representing such material properties in continuous medium mechanics yields(17)Sij0t=ρ0∂W2t∂0tεij∂0tεij,where Sij0t is Kirchhoff stress tensor, namely, the engineering stress; εij0t is Green strain tensor, namely, the engineering strain; Dijkl0t is constant and elastic constitutive tensor; Wt is the strain energy function.If the materials are plastic, then(18)Sij0t=Dijkl0εkl0t,where Dijkl0 is tangent constitutive tensor on time t, referring to configuration on time t0: (19)Dijkl0=Deijkl-Dpijkl0.


## 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:xi0, ui0.The coordinates and displacement of each point on the object at timet: xit, uit.The coordinates and displacement of each point on the object at timet+Δt: xit+Δt, uit+Δt.The virtual displacement principle corresponding to the equilibrium condition of the configuration at timet+Δt can be written as(20)∫Vt+Δtτijt+Δtδt+Δteijt+ΔtdV=∫Vt+Δttkt+Δtt+Δtδukt+ΔtdS+∫Vt+Δtρt+Δtfkt+Δtt+Δtδukt+ΔtdV,where δuk is the variation of the real displacement; δt+Δteijt+Δt is the variation of the strain; tkt+Δtt+Δt is the surface load measured in the real configuration; ft+Δtt+Δt is volume load measured in the real configuration; ρt+Δt is the mass density of the real configuration.For total Lagrange Formulation, we have(21)∫0VSij0t+Δtδεij00t+ΔtdV=∫0Stk0t+Δtδuk0dS+∫0Vρ0fkδt+Δtuk0dV.Notes: because OpenSees (the Open System for Earthquake Engineering simulation, is an object-oriented, 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 coordinatesxi0 and displacement uit of each element can be obtained by nodal value interpolation:(22)xi0=∑k=1nNkxik0,xit=∑k=1nNkxikt,xit+Δt=∑k=1nNkxikt+Δt,uit=∑k=1nNkuikt,uit+Δt=∑k=1nNkuikt+Δt.Then analyze an element by using the total Lagrange Formulation, so that we finally obtain(23)KL0t+KNL0tu=Qt+Δt-F0t,KL0t=∑∫0VBLT0tD0BL00tdV,KNL0t=∑∫0VBNLT0tS0tBNL00tdV,F0t=∑∫0VB0tLTS^00tdV,where u is the displacement vector of the node; BL0t is the transformation matrix of the linear strain displacement; BNL0t is the transformation matrix of the nonlinear strain displacement; D0 is the material constitutive matrix; S0t 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 Figures4(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.Figure 4
Collapse simulations by small strain elements in BCSS.
(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 Section2.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.Figure 5
Collapse simulations by large strain elements in BCSS.
(a)
Initial failure
(b)
No snap through element during collapsing


## 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 Figures4(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.Figure 4
Collapse simulations by small strain elements in BCSS.
(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 Section2.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.Figure 5
Collapse simulations by large strain elements in BCSS.
(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 box-type 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 Figure6. And so the map between feature parameters and geometry coordinates has been created.Figure 6
A single-box dual-room section using feature parameter methods.There are many component section types in actual project. Each section has its unique mapping relationship from feature parameters to three-dimensional 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 Space-Partitioning (BSP) Trees

A section is encircled by several directed edges, according to which its binary space-partitioning (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.Figure7 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.Figure 7
Definition of hollow sections with Boolean operation.LetP and Q be two polygons in 2D space. The difference of polygons P and Q, means Q is subtracted from P, which is denoted as P∖Q. And define ¬Q as taking negation of Q. Thus the calculation formula can be written as follows:(24)P∖Q=P∩¬Q.Take the quadrilateral sectionEFGD in Figure 8 as an example. Now it is needed to add a quadrilateral holeefgd in its internal. The BSP tree of this section after digging a hole is as shown in Figure 8.Figure 8
The BSP tree of this hollow section.
(a)
(b)As is shown in Figure8, the operation of digging a quadrilateral holeefgd in sectionEFGD is, in fact, taking C4 as a new 2D spatial region and then partitioning this region using quadrilateral sectionefgd. Accordingly, the BSP tree is to add more new nodes on the basis of the original one. According to formula P∖Q=P∩¬Q, where P is the sectionEFGD and Q is the holeefgd, thus ¬Q is to do negation of sectionefgd, that is, changing directions of the directed edges. As Figure 8 shows, the direction of the directed edges of sectionefgd is already opposite to the one of sectionEFGD. On the basis of the right-hand rule, we can get the section with hole C4∩C5∩C6∩C7, and from the BSP tree of this section we can see C4, C5, C6, and C7 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 Figure9.Figure 9
The display examples in BCSS.
(a)
“1-box and 3-room” girder section
(b)
“2-box and 1-room” 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 ofP onto the plane S, define the vector of plane S as S·n→+d=0, where S(a,b,c),n→(x,y,z), as shown in Figure 10.Figure 10
The projection of pointP onto plane S.In the picture,P′ is the projection of P onto plane S. The line PP′ connecting P and its projection P′, by definition, are parallel to the normal vector n→ of the plane S.Let the length ofPP′ be r. We obtain(25)m→=P-P′,m→=r·n^n^=n→n→n→=x2+y2+z22.So we have(26)P′=P-r·n→.As 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 Figure11, the standard plane is a torus, which is in theXOY plane. And the projection plane is created by three points: P1=(10,0,0), P2=(0,10,0), and P3=(0,0,10). Figure 12 shows the 3D model of box girder and arch by projection.Figure 11
The 3D model established by projection in BCSS.Figure 12
The 3D models by projection in BCSS.
(a)
Box girder model
(b)
Arch model


### 3.4. Organization Form of Master-Slave 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: master-slave; 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 master-slave 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.Figure 13
Constraints on geometry of all components.(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.Figure 14
The organization of node coupling.For so many arch bridge components, slave components are attached to one or more master components by position or geometry parameters. The master-slave component diagram can be established as in Figure15.Figure 15
The master-slave component diagram.In Figure15, 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.


## 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 box-type 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 Figure6. And so the map between feature parameters and geometry coordinates has been created.Figure 6
A single-box dual-room section using feature parameter methods.There are many component section types in actual project. Each section has its unique mapping relationship from feature parameters to three-dimensional 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 Space-Partitioning (BSP) Trees

A section is encircled by several directed edges, according to which its binary space-partitioning (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.Figure7 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.Figure 7
Definition of hollow sections with Boolean operation.LetP and Q be two polygons in 2D space. The difference of polygons P and Q, means Q is subtracted from P, which is denoted as P∖Q. And define ¬Q as taking negation of Q. Thus the calculation formula can be written as follows:(24)P∖Q=P∩¬Q.Take the quadrilateral sectionEFGD in Figure 8 as an example. Now it is needed to add a quadrilateral holeefgd in its internal. The BSP tree of this section after digging a hole is as shown in Figure 8.Figure 8
The BSP tree of this hollow section.
(a)
(b)As is shown in Figure8, the operation of digging a quadrilateral holeefgd in sectionEFGD is, in fact, taking C4 as a new 2D spatial region and then partitioning this region using quadrilateral sectionefgd. Accordingly, the BSP tree is to add more new nodes on the basis of the original one. According to formula P∖Q=P∩¬Q, where P is the sectionEFGD and Q is the holeefgd, thus ¬Q is to do negation of sectionefgd, that is, changing directions of the directed edges. As Figure 8 shows, the direction of the directed edges of sectionefgd is already opposite to the one of sectionEFGD. On the basis of the right-hand rule, we can get the section with hole C4∩C5∩C6∩C7, and from the BSP tree of this section we can see C4, C5, C6, and C7 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 Figure9.Figure 9
The display examples in BCSS.
(a)
“1-box and 3-room” girder section
(b)
“2-box and 1-room” 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 ofP onto the plane S, define the vector of plane S as S·n→+d=0, where S(a,b,c),n→(x,y,z), as shown in Figure 10.Figure 10
The projection of pointP onto plane S.In the picture,P′ is the projection of P onto plane S. The line PP′ connecting P and its projection P′, by definition, are parallel to the normal vector n→ of the plane S.Let the length ofPP′ be r. We obtain(25)m→=P-P′,m→=r·n^n^=n→n→n→=x2+y2+z22.So we have(26)P′=P-r·n→.As 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 Figure11, the standard plane is a torus, which is in theXOY plane. And the projection plane is created by three points: P1=(10,0,0), P2=(0,10,0), and P3=(0,0,10). Figure 12 shows the 3D model of box girder and arch by projection.Figure 11
The 3D model established by projection in BCSS.Figure 12
The 3D models by projection in BCSS.
(a)
Box girder model
(b)
Arch model


## 3.4. Organization Form of Master-Slave 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: master-slave; 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 master-slave 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.Figure 13
Constraints on geometry of all components.(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.Figure 14
The organization of node coupling.For so many arch bridge components, slave components are attached to one or more master components by position or geometry parameters. The master-slave component diagram can be established as in Figure15.Figure 15
The master-slave component diagram.In Figure15, 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 three-dimensional objects. With OpenGL, programmers can control computer-graphics 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 Figure16.Figure 16
Steps on collapse simulation using BCSS.


## 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 rise-span ratio of main arch is 1/5, the arch-axis 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 rise-span 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 Figure17.Figure 17
Elevation of Dixituojiang Bridge.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 Table1.Table 1
The material parameters of Dixituojiang Bridge.
Components
Bulk modulus(Mpa)
Shear modulus(MPa)
Density(kg/m3)
Main arch
11078
5539
2400
Cross wall
6083
3042
2400
Spandrel arch
6083
3042
2400
Bridge pier
9380
4690
2400
Floor system
6083
3042
2400With the material parameters above, a 3D FEM model of Dixituojiang Bridge is established by feature modeling method, as shown in Figure18.Figure 18
The 3D solid simulation model of Dixituojiang Bridge in BSCC.


### 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 zero-length contact element was set between nodes in foundation and in the bottom of the pier.Due to the geometric nonlinear element analysis, and setting zero-length 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 8-node hexahedral elements, and set 24 zero-length 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εmax⁡>3.5×10-2 or εmin⁡<-4.0×10-3; for the other main arches, elements fail when εmax⁡>5×10-2 or εmin⁡<-5×10-3. Then we can get the progressive collapse process similar with the actual situation, as shown in Figure 19.Figure 19
Collapse simulations for Dixituojiang Bridge in BCSS.
(a)
Initial failure is occurred ont=0.5 s
(b)
The first main arch collapses ont=2.3 s
(c)
The second main arch collapses ont=5.1 s
(d)
The third main arch collapses ont=7.3 s
(e)
The forth main arch collapses ont=7.8 s
(f)
The whole bridge collapses ont=9.1 sThe 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 Figures19(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.Aboutt=0.5 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 t=0.9 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.


## 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 Table1.Table 1
The material parameters of Dixituojiang Bridge.
Components
Bulk modulus(Mpa)
Shear modulus(MPa)
Density(kg/m3)
Main arch
11078
5539
2400
Cross wall
6083
3042
2400
Spandrel arch
6083
3042
2400
Bridge pier
9380
4690
2400
Floor system
6083
3042
2400With the material parameters above, a 3D FEM model of Dixituojiang Bridge is established by feature modeling method, as shown in Figure18.Figure 18
The 3D solid simulation model of Dixituojiang Bridge in BSCC.


## 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 zero-length contact element was set between nodes in foundation and in the bottom of the pier.Due to the geometric nonlinear element analysis, and setting zero-length 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 8-node hexahedral elements, and set 24 zero-length 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εmax⁡>3.5×10-2 or εmin⁡<-4.0×10-3; for the other main arches, elements fail when εmax⁡>5×10-2 or εmin⁡<-5×10-3. Then we can get the progressive collapse process similar with the actual situation, as shown in Figure 19.Figure 19
Collapse simulations for Dixituojiang Bridge in BCSS.
(a)
Initial failure is occurred ont=0.5 s
(b)
The first main arch collapses ont=2.3 s
(c)
The second main arch collapses ont=5.1 s
(d)
The third main arch collapses ont=7.3 s
(e)
The forth main arch collapses ont=7.8 s
(f)
The whole bridge collapses ont=9.1 sThe 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 Figures19(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.Aboutt=0.5 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 t=0.9 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.



---

*Source: 707269-2015-05-18.xml*