# Finite Particle Method-Based Collapse Simulation of Space Steel Frame Subjected to Earthquake Excitation

Authors: Xiao-Hong Long; Rong Yue; Yong-Tao Ma; Jian Fan
Journal: Shock and Vibration
 (2018)

---

## Abstract

In the process of collapse failure of the space steel frame subjected to earthquake excitation, complex behaviors often are involved, including geometric nonlinearity, material nonlinearity, fracture, contact, and collisions. In view of the unique advantages of the finite particle method to analyze complex structural nonlinear problems, this paper utilized the finite particle method as the basic means of analysis and used MATLAB software for computational analysis. This paper first derived a finite particle method-based space steel frame model, conducted static analysis and dynamic response analysis under earthquake excitation, and compared findings with ANSYS analysis results to validate reliability. This paper established the fracture criterion and failure mode of a steel frame member. Theoretical derivation and numerical simulation indicate that the finite particle method is a feasible and effective way to simulate the collapse of space steel frame structures subjected to earthquake excitation. This method provides a new approach to study the collapse and anticollapse seismic design of space steel frame structures subjected to earthquake excitation.


---

## Body

## 1. Introduction

By investigating the damage conditions of buildings in large earthquakes occurring globally, it is evident that it is not uncommon for steel structures to be damaged or even collapse under the excitation of strong earthquakes. In 1994, an earthquake occurred in Northridge, California, USA. The Northridge earthquake caused various degrees of damage to more than 100 steel frame buildings [1, 2]. In 1995, a strong earthquake occurred at the southern end of Binku County in Japan. This earthquake caused severe damage to the steel structural buildings [3]. Related data indicated that 988 modern steel structural buildings were subjected to damage, of which 90 buildings collapsed, and 322 buildings were severely damaged [4]. To date, theoretical analysis, experimental investigation, and numerical simulation have been the major approaches used to study the structural collapse process proposed by various researchers. The failure process of structural collapse often involves a complex nonlinear situation, and it is difficult to carry out pure theoretical analysis and accurately understand the failure mechanism of the structure. The structural collapse experiments also are limited to testing conditions, difficulty in execution, and higher cost, and the experimental process is difficult to precisely control. With the development of computer and structural analysis technologies, various mechanical model-based numerical simulation methods have been proposed and gradually have improved, thus becoming a practical and effective study approach.Scholars have utilized the discrete element method, finite element method, and particle method to numerically simulate structural collapse and have achieved many favorable outcomes that have provided important help into the study of the continuous collapse mechanism [5–7]. Structure collapse, however, is an extremely complex process that involves dynamic nonlinearity of geometry and materials, element fracture, contact, collisions, and other problems. The conventional methods often need special treatment, which significantly limits determination of the problem. Therefore, how to find more effective numerical simulation is an important topic in the study of structural collapse. Facing this difficulty of applying a conventional numerical analysis method to stimulate complex structural behaviors, Ding [8, 9] proposed a new concept of selecting a point value commonly used in the numerical computation method as a means to describe structural behavior. He introduced generalized vector mechanics as the criterion of movement and deformation, and accordingly, he developed vector structure and solid mechanical theory [10] and provided the theoretical foundation for the development of finite particle method. On the basis of principle theories of vector mechanics, Yu et al. [11, 12] developed a finite particle method, which is a structural analysis method focused more on structural engineering. The basic idea of this analysis method is to discrete the analysis field into a series of particles, connect particles through elements, describe the mechanical relationship between particles through the deformation of elements, follow Newton’s motion law for the movement of particles, and track the structural behaviors through the movement of particles. Because the method maintains the energy equilibrium of each particle during analysis, there is no need for iteration or special correction, nor generalized stiffness matrix. The method is not limited by mesh division during simulation of the fracture of the member, and it can effectively deal with large geometric deformation of spatial structures, the constitutive relationship of nonlinear and discontinuous materials, and the movement of the structural body and the rigid body and their mutual coupling behavior, as well as other complex behaviors. Wang et al. [13] tested the collapse of a single-story four-column wall frame structure and compared the experimental results with the results from the finite particle simulation. Their research indicated that compared with the conventional finite element method, the finite particle method can simulate structural collapse in a simpler and more effective way.In treating large structural deformations, large displacement, and fracture collision, the finite particle method exhibits its unique advantages in stability, calculation convergence, and simplification in program preparation. Currently, under the joint efforts of investigators, the finite particle method has achieved a certain popularity and development; however, study of the collapse of the steel frame under earthquake excitation is rare. This paper utilizes the finite particle method as a basic analysis approach, uses MATLAB software as computation analysis tool, combines theoretical derivation and numerical simulation, and carries out a simulation analysis on the collapse failure process of steel frame under earthquake excitation.


## 2. Description of Finite Particle Method for Three-Dimensional Beam Element

The finite particle method is based on vector mechanics and describes(1) a structure spatially by discrete particles and (2) movement of those particles by path elements in time. Accordingly, this method accurately describes the systematic behaviors of a structure. The movement of particles follows Newton’s second law, the interaction between particles is simulated through elements, and the pure deformation of the elements is obtained by eliminating rigid body motion through reverse movement [11].


### 2.1. Establishment of a Particle Motion Equation

In the finite particle method, the particle displacement of the space beam element can be dissected into linear displacements in three directions and three angle displacements, which correspond to three forces and three moments in the coordinate axis, respectively, as shown in Figure1.Figure 1
Loading onto space beam system particles.The motion equation in path elements for particles in the space beam structure can be expressed as shown in (1) and (2):(1)Mαd2dt2dxdydzα=FxintFyintFzintα+FxextFyextFzextα,and(2)Ιαd2dt2θxθyθzα=IxxIxyIxzIyxIyyIyzIzxIzyIzzd2dt2θxθyθzα=MxintMyintMzintα+MxextMyextMzextα,where Mα is the mass of particle α; dxdydzαT is the displacement vector of particle α; FxintFyintFzintαT is the internal force vector of particle; FxextFyextFzextαT is the external force vector; Ια is the mass moment matrix of inertia of particle; θxθyθzαT is the rotational angle vector of particle; MxintMyintMzintαT is the internal moment vector of particle; and MxextMyextMzextαT is the external moment vector of particle.The internal force is from the beam elements connecting to the particles. The internal force reaction produced after element deformation is applied onto the particles with which the element connects. This chapter introduces the process of calculating internal force with the finite particle method using the steps of local coordinate system transformation, deformation calculation, and internal force calculation.


#### 2.1.1. Transformation of Local Coordinate System

Assume the particles at the two ends of an element at the initial time are A and B, respectively. To determine the pure deformation of the element, set up a group of local coordinate systemsO^x^y^z^ as a deformation coordinate system that moves along with the movement of the element. At the initial time t0, the direction in x^ is given as the length direction of the element. The y^ axis and z^ axis are perpendicular to the x^ axis plane as shown in Figure 2. Select a reference point C, assume y^ axis is perpendicular to the plane consisting of AB and AC, and then the unit vector of each principal can be expressed in (3a), (3b), and (3c):(3a)e^x=eAB,(3b)e^y=eAC×eABeAC×eAB,and(3c)e^z=e^x×e^y.Figure 2
Local coordinate system of space beam element.The local coordinate system of the element keeps continuous change during the movement of the structure, but it only changes at the time of cutting two path elements, and the direction within one pathway element does not change. Assume the time section fromta to tb is the pathway element describing particle movement, and then the vectors for element AB at the two continuous time points in the principal direction of deformation coordinate system are e^xa,e^ya,e^za and e^xb,e^yb,e^zb.Given that the space coordinates for the two particles of an element at theta and ta are xAa,xBa and xAb,xBb, respectively, then the particle rotational angles are βAa,βBa and βAb,βBb, respectively. As shown in Figure 3, the displacement Δxi and rotational angle Δβi for the two particles of an element within that pathway element can be expressed by (4)Δxi=xib-xiaΔβi=βib-βiai=A,B.Figure 3
Particle displacement.The principal axis rotational vector of the element consists of out-of-axis rotation and self-rotation of the principal axise^xa during the reverse movement. The out-of-axis rotation of the principal axis e^xa, as shown in Figure 4, can be obtained through the location vectors xia and xia for the particle at the time of ta and tb, and the principal axis e^xa and e^xb can be expressed by(5a)e^xa=xBa-xAaxBa-xAa,and(5b)e^xb=xBb-xAbxBb-xAb.Figure 4
Rotational angle of space beam element in principal axis direction.Then the out-of-axis rotational vectorθba can be expressed as (6a): (6a)θba=θbaeba,(6b)θba=sin-1e^xa×e^xb,and(6c)eba=e^xa×e^xbe^xa×e^xb,where θba represents the magnitude of rotational angle and eba represents the rotational angle direction vector.The self-rotation of the principal axis can be obtained by projecting the rotational angleΔβA of node A within pathway element onto the e^xa axis as presented by(7)Δβ^xA=ΔβA·e^xa.Thus, the total rotational vector of the principal axis is obtained as(8)γ=θba+Δβ^xA,where the magnitude of the rotational vector γ is γ and axial direction is eλ=lγmγnγ. Then, the rotational matrix of the principal axis can be obtained as give in(9)R∗=I+1-cos⁡γA2+sin⁡γA,where I represents a 3×3 unit matrix and matrix A is expressed as(10)A=0-nγmγnγ0-lγ-mγlγ0.Accordingly the principal axise^yb and e^zb of element at tb can be obtained as(11a)e^yb=R∗e^ya,and(11b)e^zb=e^xb×e^yb.


#### 2.1.2. Pure Deformation Calculation of the Space Beam Element

After determining the direction vectors(e^xa,e^ya,e^za) and (e^xb,e^yb,e^zb) corresponding to times ta and tb, we assume the configuration at time ta to be the virtual reference structural shape. The transformation relationship from the global coordinate system to the local coordinate system for the vector variables under the virtual configuration can be expressed as follows: (12a)v^=Ωv,of which (12b)Ω=e^xaTe^yaTe^zaT.For elementΑ′Β′, we assume that reverse translation is -ΔxA and reverse rotation is -γ to Α′′Β′′, as shown in Figure 5; then, the difference between the virtual location Α′′Β′′ for the element through reverse movement and the reference structural shape AB is referred to as the pure deformation of element.Figure 5
Reverse movement of space beam element.We establish deformation coordinates at the principal vector(e^xa,e^ya,e^za) at time ta. Because the element only includes two particles, A and B, and the virtual location Α′′Β′′ is on the same line with AB, the axial elongation and compressive deformation are the change of the axial length of the element as represented by(13)ΔuBd=Δee^xa=lb-lae^xa,of which la and lb are the lengths corresponding to times ta and tb, respectively. To calculate the flexure and rotational angle between the virtual locations Α′′Β′′ and AB, transform the rotational angle displacement of the particle from the global coordinate system to the local coordinate system of the element according to(14)Δβ^i=ΩΔβi,i=Α,Β.Obtain the rotational angles for nodes A and B as follows:(15a)ΔφxΑ=ΔβxΑ-ΔβxΑ=0,and(15b)ΔφxΒ=ΔβxΒ-ΔβxΑ.The out-of-axis rotational vectorθba is represented by principal coordinates, and the corresponding three principal components are presented by(16a)Δθx=θba·e^xa=0,(16b)Δθy=θba·e^ya,and(16c)Δθz=θba·e^za.After the reverse out-of-axis rotation, as shown by the rotation angle relationship in Figure6, the rotation angles for particle A and B are provided by (17a) through (17d):(17a)ΔφyΑ=Δβ^yΑ-Δθy,(17b)ΔφzΑ=Δβ^zΑ-Δθz,(17c)ΔφyΒ=Δβ^yΒ-Δθy,and(17d)ΔφzΒ=Δβ^zΒ-Δθz.Figure 6
Bending angle of node A.Accordingly, we can obtain six pure deformation diagrams for the element (Figure7): axial deformation Δe, principal rotation ΔφxΒ, two flexure deformations ΔφyΑ and ΔφzΑ for node A, and two flexure deformations ΔφyΒ and ΔφzΒ for node B.Figure 7
Pure deformation of element.


#### 2.1.3. Internal Force Calculation of Space Beam Element

From these six pure deformations of the element, we can determine the element’s internal force increment under the local coordinate system for the pathway element within timesta~tb by (18a) through (18d): (18a)
Δ
f
^
Α
=
Δ
f
^
x
Α
Δ
f
^
y
Α
Δ
f
^
z
Α
=
-
E
A
a
l
a
l
b
-
l
a
6
E
I
z
^
l
a
2
Δ
φ
z
Α
+
Δ
φ
z
Β
-
6
E
I
y
^
l
a
2
Δ
φ
y
Α
+
Δ
φ
y
Β
,
a
n
d
(18b)
Δ
f
^
Β
=
Δ
f
^
x
Β
Δ
f
^
y
Β
Δ
f
^
z
Β
=
E
A
a
l
a
l
b
-
l
a
-
6
E
I
z
^
l
a
2
Δ
φ
z
Α
+
Δ
φ
z
Β
6
E
I
y
^
l
a
2
Δ
φ
y
Α
+
Δ
φ
y
Β
,
(18c)
Δ
m
^
Α
=
Δ
m
^
x
Α
Δ
m
^
y
Α
Δ
m
^
z
Α
=
-
G
J
a
l
a
Δ
φ
x
Β
E
I
y
^
l
a
4
Δ
φ
y
Α
+
2
Δ
φ
y
Β
E
I
z
^
l
a
4
Δ
φ
z
Α
+
2
Δ
φ
z
Β
,
a
n
d
(18d)
Δ
m
^
Β
=
Δ
m
^
x
Β
Δ
m
^
y
Β
Δ
m
^
z
Β
=
G
J
a
l
a
Δ
φ
x
Β
E
I
y
^
l
a
2
Δ
φ
y
Α
+
4
Δ
φ
y
Β
E
I
z
^
l
a
2
Δ
φ
z
Α
+
4
Δ
φ
z
Β
.The above internal force increments are added to the full internal force of the element in the local coordinate system at timeta to obtain the full internal force of the element in the local coordinate system from time ta to time tb as presented in (19a)f^bi=f^ai+Δf^i,i=Α,Β,and(19b)m^bi=m^ai+Δm^i,i=Α,Β.Then through virtual positive rotation and coordinate transformation, calculate the element internal force under global coordinate system at timetb (Figure 8) as presented by (20a)fbi=fxbifybifzbi=R∗ΩTf^xbif^ybif^zbi,i=Α,Β,and(20b)mbi=mxbimybimzbi=R∗ΩTm^xbim^ybim^zbi,i=Α,Β.Figure 8
Translation and rotation in positive direction of space beam element.The internal force and moment of element AB at timetb are applied to particles at both ends by reverse superposition to obtain the total force and total moment of particle α at time tb as represented in(21a)Fα=Fαext-Fαint,and(21b)Mα=Mαext-Mαint.


### 2.2. Solution of Particle Motion Equation

The motion equation of particleα can be expressed by(22)Mαd¨α=Fαint+Fαext,and(23)Iαθ¨α=Mαint+Mαext.where d¨α represents the linear acceleration of particle and θ¨α represents the angle acceleration of particle. To solve (22) and (23), implicit solution and explicit integration methods can be used. Of these, implicit solution needs iteration and there is no issue of convergence. The finite particle method uses a time-display integration method for the solution and utilizes central difference; therefore, iteration is not needed. The velocityd and the acceleration d¨ for the particle at thenth step can be expressed by(24a)d˙n=12hdn+1-dn-1,and(24b)d¨n=1h2dn+1-2dn+dn-1.Accordingly, we can obtain the linear displacementdn+1 and angle displacement θn+1 for the particle at then+1 step as represented by(25a)dn+1=h2MαFαint+Fαext+2dn-dn-1,and(25b)θn+1=h2Iα-1Mαint+Mαext+2θn-θn-1,wherehrepresents the time step.From the equation solution processes, we found that the solution of the finite particle method does not need a generalized stiffness matrix. By applying (25a) and (25b), we can directly obtain the displacement of the particle within each pathway element, obtain the space location and deformation state for the structure at the next time step, and complete the description to structural behaviors.When utilizing the central difference method to determine the motion equation of the particle, we need to determine an appropriate time increment step to control the error accumulation resulting from each step in the difference operation to ensure the stability and accuracy of the results. According to the stability principle [14], the increment time step Δt should meet(26)Δt≤Δtcr=2ωmax,where Δtcr represents the critical time step and ωmax represents the maximum natural frequency.If we consider the impact of damping, the critical time step of the central difference method could be expressed by(27)Δtcr≈2lminνcξ2+1-ξ=2lminξ2+1-ξE/ρ,where lmin represents the minimal length of the element in the system; νc represents the transmission velocity of the axial wave in a straight linear element;E represents the elastic modulus of the material; ρ represents the material density; and ξ represents the damping ratio coefficient corresponding to the highest mode frequency of system.


## 2.1. Establishment of a Particle Motion Equation

In the finite particle method, the particle displacement of the space beam element can be dissected into linear displacements in three directions and three angle displacements, which correspond to three forces and three moments in the coordinate axis, respectively, as shown in Figure1.Figure 1
Loading onto space beam system particles.The motion equation in path elements for particles in the space beam structure can be expressed as shown in (1) and (2):(1)Mαd2dt2dxdydzα=FxintFyintFzintα+FxextFyextFzextα,and(2)Ιαd2dt2θxθyθzα=IxxIxyIxzIyxIyyIyzIzxIzyIzzd2dt2θxθyθzα=MxintMyintMzintα+MxextMyextMzextα,where Mα is the mass of particle α; dxdydzαT is the displacement vector of particle α; FxintFyintFzintαT is the internal force vector of particle; FxextFyextFzextαT is the external force vector; Ια is the mass moment matrix of inertia of particle; θxθyθzαT is the rotational angle vector of particle; MxintMyintMzintαT is the internal moment vector of particle; and MxextMyextMzextαT is the external moment vector of particle.The internal force is from the beam elements connecting to the particles. The internal force reaction produced after element deformation is applied onto the particles with which the element connects. This chapter introduces the process of calculating internal force with the finite particle method using the steps of local coordinate system transformation, deformation calculation, and internal force calculation.


### 2.1.1. Transformation of Local Coordinate System

Assume the particles at the two ends of an element at the initial time are A and B, respectively. To determine the pure deformation of the element, set up a group of local coordinate systemsO^x^y^z^ as a deformation coordinate system that moves along with the movement of the element. At the initial time t0, the direction in x^ is given as the length direction of the element. The y^ axis and z^ axis are perpendicular to the x^ axis plane as shown in Figure 2. Select a reference point C, assume y^ axis is perpendicular to the plane consisting of AB and AC, and then the unit vector of each principal can be expressed in (3a), (3b), and (3c):(3a)e^x=eAB,(3b)e^y=eAC×eABeAC×eAB,and(3c)e^z=e^x×e^y.Figure 2
Local coordinate system of space beam element.The local coordinate system of the element keeps continuous change during the movement of the structure, but it only changes at the time of cutting two path elements, and the direction within one pathway element does not change. Assume the time section fromta to tb is the pathway element describing particle movement, and then the vectors for element AB at the two continuous time points in the principal direction of deformation coordinate system are e^xa,e^ya,e^za and e^xb,e^yb,e^zb.Given that the space coordinates for the two particles of an element at theta and ta are xAa,xBa and xAb,xBb, respectively, then the particle rotational angles are βAa,βBa and βAb,βBb, respectively. As shown in Figure 3, the displacement Δxi and rotational angle Δβi for the two particles of an element within that pathway element can be expressed by (4)Δxi=xib-xiaΔβi=βib-βiai=A,B.Figure 3
Particle displacement.The principal axis rotational vector of the element consists of out-of-axis rotation and self-rotation of the principal axise^xa during the reverse movement. The out-of-axis rotation of the principal axis e^xa, as shown in Figure 4, can be obtained through the location vectors xia and xia for the particle at the time of ta and tb, and the principal axis e^xa and e^xb can be expressed by(5a)e^xa=xBa-xAaxBa-xAa,and(5b)e^xb=xBb-xAbxBb-xAb.Figure 4
Rotational angle of space beam element in principal axis direction.Then the out-of-axis rotational vectorθba can be expressed as (6a): (6a)θba=θbaeba,(6b)θba=sin-1e^xa×e^xb,and(6c)eba=e^xa×e^xbe^xa×e^xb,where θba represents the magnitude of rotational angle and eba represents the rotational angle direction vector.The self-rotation of the principal axis can be obtained by projecting the rotational angleΔβA of node A within pathway element onto the e^xa axis as presented by(7)Δβ^xA=ΔβA·e^xa.Thus, the total rotational vector of the principal axis is obtained as(8)γ=θba+Δβ^xA,where the magnitude of the rotational vector γ is γ and axial direction is eλ=lγmγnγ. Then, the rotational matrix of the principal axis can be obtained as give in(9)R∗=I+1-cos⁡γA2+sin⁡γA,where I represents a 3×3 unit matrix and matrix A is expressed as(10)A=0-nγmγnγ0-lγ-mγlγ0.Accordingly the principal axise^yb and e^zb of element at tb can be obtained as(11a)e^yb=R∗e^ya,and(11b)e^zb=e^xb×e^yb.


### 2.1.2. Pure Deformation Calculation of the Space Beam Element

After determining the direction vectors(e^xa,e^ya,e^za) and (e^xb,e^yb,e^zb) corresponding to times ta and tb, we assume the configuration at time ta to be the virtual reference structural shape. The transformation relationship from the global coordinate system to the local coordinate system for the vector variables under the virtual configuration can be expressed as follows: (12a)v^=Ωv,of which (12b)Ω=e^xaTe^yaTe^zaT.For elementΑ′Β′, we assume that reverse translation is -ΔxA and reverse rotation is -γ to Α′′Β′′, as shown in Figure 5; then, the difference between the virtual location Α′′Β′′ for the element through reverse movement and the reference structural shape AB is referred to as the pure deformation of element.Figure 5
Reverse movement of space beam element.We establish deformation coordinates at the principal vector(e^xa,e^ya,e^za) at time ta. Because the element only includes two particles, A and B, and the virtual location Α′′Β′′ is on the same line with AB, the axial elongation and compressive deformation are the change of the axial length of the element as represented by(13)ΔuBd=Δee^xa=lb-lae^xa,of which la and lb are the lengths corresponding to times ta and tb, respectively. To calculate the flexure and rotational angle between the virtual locations Α′′Β′′ and AB, transform the rotational angle displacement of the particle from the global coordinate system to the local coordinate system of the element according to(14)Δβ^i=ΩΔβi,i=Α,Β.Obtain the rotational angles for nodes A and B as follows:(15a)ΔφxΑ=ΔβxΑ-ΔβxΑ=0,and(15b)ΔφxΒ=ΔβxΒ-ΔβxΑ.The out-of-axis rotational vectorθba is represented by principal coordinates, and the corresponding three principal components are presented by(16a)Δθx=θba·e^xa=0,(16b)Δθy=θba·e^ya,and(16c)Δθz=θba·e^za.After the reverse out-of-axis rotation, as shown by the rotation angle relationship in Figure6, the rotation angles for particle A and B are provided by (17a) through (17d):(17a)ΔφyΑ=Δβ^yΑ-Δθy,(17b)ΔφzΑ=Δβ^zΑ-Δθz,(17c)ΔφyΒ=Δβ^yΒ-Δθy,and(17d)ΔφzΒ=Δβ^zΒ-Δθz.Figure 6
Bending angle of node A.Accordingly, we can obtain six pure deformation diagrams for the element (Figure7): axial deformation Δe, principal rotation ΔφxΒ, two flexure deformations ΔφyΑ and ΔφzΑ for node A, and two flexure deformations ΔφyΒ and ΔφzΒ for node B.Figure 7
Pure deformation of element.


### 2.1.3. Internal Force Calculation of Space Beam Element

From these six pure deformations of the element, we can determine the element’s internal force increment under the local coordinate system for the pathway element within timesta~tb by (18a) through (18d): (18a)
Δ
f
^
Α
=
Δ
f
^
x
Α
Δ
f
^
y
Α
Δ
f
^
z
Α
=
-
E
A
a
l
a
l
b
-
l
a
6
E
I
z
^
l
a
2
Δ
φ
z
Α
+
Δ
φ
z
Β
-
6
E
I
y
^
l
a
2
Δ
φ
y
Α
+
Δ
φ
y
Β
,
a
n
d
(18b)
Δ
f
^
Β
=
Δ
f
^
x
Β
Δ
f
^
y
Β
Δ
f
^
z
Β
=
E
A
a
l
a
l
b
-
l
a
-
6
E
I
z
^
l
a
2
Δ
φ
z
Α
+
Δ
φ
z
Β
6
E
I
y
^
l
a
2
Δ
φ
y
Α
+
Δ
φ
y
Β
,
(18c)
Δ
m
^
Α
=
Δ
m
^
x
Α
Δ
m
^
y
Α
Δ
m
^
z
Α
=
-
G
J
a
l
a
Δ
φ
x
Β
E
I
y
^
l
a
4
Δ
φ
y
Α
+
2
Δ
φ
y
Β
E
I
z
^
l
a
4
Δ
φ
z
Α
+
2
Δ
φ
z
Β
,
a
n
d
(18d)
Δ
m
^
Β
=
Δ
m
^
x
Β
Δ
m
^
y
Β
Δ
m
^
z
Β
=
G
J
a
l
a
Δ
φ
x
Β
E
I
y
^
l
a
2
Δ
φ
y
Α
+
4
Δ
φ
y
Β
E
I
z
^
l
a
2
Δ
φ
z
Α
+
4
Δ
φ
z
Β
.The above internal force increments are added to the full internal force of the element in the local coordinate system at timeta to obtain the full internal force of the element in the local coordinate system from time ta to time tb as presented in (19a)f^bi=f^ai+Δf^i,i=Α,Β,and(19b)m^bi=m^ai+Δm^i,i=Α,Β.Then through virtual positive rotation and coordinate transformation, calculate the element internal force under global coordinate system at timetb (Figure 8) as presented by (20a)fbi=fxbifybifzbi=R∗ΩTf^xbif^ybif^zbi,i=Α,Β,and(20b)mbi=mxbimybimzbi=R∗ΩTm^xbim^ybim^zbi,i=Α,Β.Figure 8
Translation and rotation in positive direction of space beam element.The internal force and moment of element AB at timetb are applied to particles at both ends by reverse superposition to obtain the total force and total moment of particle α at time tb as represented in(21a)Fα=Fαext-Fαint,and(21b)Mα=Mαext-Mαint.


## 2.1.1. Transformation of Local Coordinate System

Assume the particles at the two ends of an element at the initial time are A and B, respectively. To determine the pure deformation of the element, set up a group of local coordinate systemsO^x^y^z^ as a deformation coordinate system that moves along with the movement of the element. At the initial time t0, the direction in x^ is given as the length direction of the element. The y^ axis and z^ axis are perpendicular to the x^ axis plane as shown in Figure 2. Select a reference point C, assume y^ axis is perpendicular to the plane consisting of AB and AC, and then the unit vector of each principal can be expressed in (3a), (3b), and (3c):(3a)e^x=eAB,(3b)e^y=eAC×eABeAC×eAB,and(3c)e^z=e^x×e^y.Figure 2
Local coordinate system of space beam element.The local coordinate system of the element keeps continuous change during the movement of the structure, but it only changes at the time of cutting two path elements, and the direction within one pathway element does not change. Assume the time section fromta to tb is the pathway element describing particle movement, and then the vectors for element AB at the two continuous time points in the principal direction of deformation coordinate system are e^xa,e^ya,e^za and e^xb,e^yb,e^zb.Given that the space coordinates for the two particles of an element at theta and ta are xAa,xBa and xAb,xBb, respectively, then the particle rotational angles are βAa,βBa and βAb,βBb, respectively. As shown in Figure 3, the displacement Δxi and rotational angle Δβi for the two particles of an element within that pathway element can be expressed by (4)Δxi=xib-xiaΔβi=βib-βiai=A,B.Figure 3
Particle displacement.The principal axis rotational vector of the element consists of out-of-axis rotation and self-rotation of the principal axise^xa during the reverse movement. The out-of-axis rotation of the principal axis e^xa, as shown in Figure 4, can be obtained through the location vectors xia and xia for the particle at the time of ta and tb, and the principal axis e^xa and e^xb can be expressed by(5a)e^xa=xBa-xAaxBa-xAa,and(5b)e^xb=xBb-xAbxBb-xAb.Figure 4
Rotational angle of space beam element in principal axis direction.Then the out-of-axis rotational vectorθba can be expressed as (6a): (6a)θba=θbaeba,(6b)θba=sin-1e^xa×e^xb,and(6c)eba=e^xa×e^xbe^xa×e^xb,where θba represents the magnitude of rotational angle and eba represents the rotational angle direction vector.The self-rotation of the principal axis can be obtained by projecting the rotational angleΔβA of node A within pathway element onto the e^xa axis as presented by(7)Δβ^xA=ΔβA·e^xa.Thus, the total rotational vector of the principal axis is obtained as(8)γ=θba+Δβ^xA,where the magnitude of the rotational vector γ is γ and axial direction is eλ=lγmγnγ. Then, the rotational matrix of the principal axis can be obtained as give in(9)R∗=I+1-cos⁡γA2+sin⁡γA,where I represents a 3×3 unit matrix and matrix A is expressed as(10)A=0-nγmγnγ0-lγ-mγlγ0.Accordingly the principal axise^yb and e^zb of element at tb can be obtained as(11a)e^yb=R∗e^ya,and(11b)e^zb=e^xb×e^yb.


## 2.1.2. Pure Deformation Calculation of the Space Beam Element

After determining the direction vectors(e^xa,e^ya,e^za) and (e^xb,e^yb,e^zb) corresponding to times ta and tb, we assume the configuration at time ta to be the virtual reference structural shape. The transformation relationship from the global coordinate system to the local coordinate system for the vector variables under the virtual configuration can be expressed as follows: (12a)v^=Ωv,of which (12b)Ω=e^xaTe^yaTe^zaT.For elementΑ′Β′, we assume that reverse translation is -ΔxA and reverse rotation is -γ to Α′′Β′′, as shown in Figure 5; then, the difference between the virtual location Α′′Β′′ for the element through reverse movement and the reference structural shape AB is referred to as the pure deformation of element.Figure 5
Reverse movement of space beam element.We establish deformation coordinates at the principal vector(e^xa,e^ya,e^za) at time ta. Because the element only includes two particles, A and B, and the virtual location Α′′Β′′ is on the same line with AB, the axial elongation and compressive deformation are the change of the axial length of the element as represented by(13)ΔuBd=Δee^xa=lb-lae^xa,of which la and lb are the lengths corresponding to times ta and tb, respectively. To calculate the flexure and rotational angle between the virtual locations Α′′Β′′ and AB, transform the rotational angle displacement of the particle from the global coordinate system to the local coordinate system of the element according to(14)Δβ^i=ΩΔβi,i=Α,Β.Obtain the rotational angles for nodes A and B as follows:(15a)ΔφxΑ=ΔβxΑ-ΔβxΑ=0,and(15b)ΔφxΒ=ΔβxΒ-ΔβxΑ.The out-of-axis rotational vectorθba is represented by principal coordinates, and the corresponding three principal components are presented by(16a)Δθx=θba·e^xa=0,(16b)Δθy=θba·e^ya,and(16c)Δθz=θba·e^za.After the reverse out-of-axis rotation, as shown by the rotation angle relationship in Figure6, the rotation angles for particle A and B are provided by (17a) through (17d):(17a)ΔφyΑ=Δβ^yΑ-Δθy,(17b)ΔφzΑ=Δβ^zΑ-Δθz,(17c)ΔφyΒ=Δβ^yΒ-Δθy,and(17d)ΔφzΒ=Δβ^zΒ-Δθz.Figure 6
Bending angle of node A.Accordingly, we can obtain six pure deformation diagrams for the element (Figure7): axial deformation Δe, principal rotation ΔφxΒ, two flexure deformations ΔφyΑ and ΔφzΑ for node A, and two flexure deformations ΔφyΒ and ΔφzΒ for node B.Figure 7
Pure deformation of element.


## 2.1.3. Internal Force Calculation of Space Beam Element

From these six pure deformations of the element, we can determine the element’s internal force increment under the local coordinate system for the pathway element within timesta~tb by (18a) through (18d): (18a)
Δ
f
^
Α
=
Δ
f
^
x
Α
Δ
f
^
y
Α
Δ
f
^
z
Α
=
-
E
A
a
l
a
l
b
-
l
a
6
E
I
z
^
l
a
2
Δ
φ
z
Α
+
Δ
φ
z
Β
-
6
E
I
y
^
l
a
2
Δ
φ
y
Α
+
Δ
φ
y
Β
,
a
n
d
(18b)
Δ
f
^
Β
=
Δ
f
^
x
Β
Δ
f
^
y
Β
Δ
f
^
z
Β
=
E
A
a
l
a
l
b
-
l
a
-
6
E
I
z
^
l
a
2
Δ
φ
z
Α
+
Δ
φ
z
Β
6
E
I
y
^
l
a
2
Δ
φ
y
Α
+
Δ
φ
y
Β
,
(18c)
Δ
m
^
Α
=
Δ
m
^
x
Α
Δ
m
^
y
Α
Δ
m
^
z
Α
=
-
G
J
a
l
a
Δ
φ
x
Β
E
I
y
^
l
a
4
Δ
φ
y
Α
+
2
Δ
φ
y
Β
E
I
z
^
l
a
4
Δ
φ
z
Α
+
2
Δ
φ
z
Β
,
a
n
d
(18d)
Δ
m
^
Β
=
Δ
m
^
x
Β
Δ
m
^
y
Β
Δ
m
^
z
Β
=
G
J
a
l
a
Δ
φ
x
Β
E
I
y
^
l
a
2
Δ
φ
y
Α
+
4
Δ
φ
y
Β
E
I
z
^
l
a
2
Δ
φ
z
Α
+
4
Δ
φ
z
Β
.The above internal force increments are added to the full internal force of the element in the local coordinate system at timeta to obtain the full internal force of the element in the local coordinate system from time ta to time tb as presented in (19a)f^bi=f^ai+Δf^i,i=Α,Β,and(19b)m^bi=m^ai+Δm^i,i=Α,Β.Then through virtual positive rotation and coordinate transformation, calculate the element internal force under global coordinate system at timetb (Figure 8) as presented by (20a)fbi=fxbifybifzbi=R∗ΩTf^xbif^ybif^zbi,i=Α,Β,and(20b)mbi=mxbimybimzbi=R∗ΩTm^xbim^ybim^zbi,i=Α,Β.Figure 8
Translation and rotation in positive direction of space beam element.The internal force and moment of element AB at timetb are applied to particles at both ends by reverse superposition to obtain the total force and total moment of particle α at time tb as represented in(21a)Fα=Fαext-Fαint,and(21b)Mα=Mαext-Mαint.


## 2.2. Solution of Particle Motion Equation

The motion equation of particleα can be expressed by(22)Mαd¨α=Fαint+Fαext,and(23)Iαθ¨α=Mαint+Mαext.where d¨α represents the linear acceleration of particle and θ¨α represents the angle acceleration of particle. To solve (22) and (23), implicit solution and explicit integration methods can be used. Of these, implicit solution needs iteration and there is no issue of convergence. The finite particle method uses a time-display integration method for the solution and utilizes central difference; therefore, iteration is not needed. The velocityd and the acceleration d¨ for the particle at thenth step can be expressed by(24a)d˙n=12hdn+1-dn-1,and(24b)d¨n=1h2dn+1-2dn+dn-1.Accordingly, we can obtain the linear displacementdn+1 and angle displacement θn+1 for the particle at then+1 step as represented by(25a)dn+1=h2MαFαint+Fαext+2dn-dn-1,and(25b)θn+1=h2Iα-1Mαint+Mαext+2θn-θn-1,wherehrepresents the time step.From the equation solution processes, we found that the solution of the finite particle method does not need a generalized stiffness matrix. By applying (25a) and (25b), we can directly obtain the displacement of the particle within each pathway element, obtain the space location and deformation state for the structure at the next time step, and complete the description to structural behaviors.When utilizing the central difference method to determine the motion equation of the particle, we need to determine an appropriate time increment step to control the error accumulation resulting from each step in the difference operation to ensure the stability and accuracy of the results. According to the stability principle [14], the increment time step Δt should meet(26)Δt≤Δtcr=2ωmax,where Δtcr represents the critical time step and ωmax represents the maximum natural frequency.If we consider the impact of damping, the critical time step of the central difference method could be expressed by(27)Δtcr≈2lminνcξ2+1-ξ=2lminξ2+1-ξE/ρ,where lmin represents the minimal length of the element in the system; νc represents the transmission velocity of the axial wave in a straight linear element;E represents the elastic modulus of the material; ρ represents the material density; and ξ represents the damping ratio coefficient corresponding to the highest mode frequency of system.


## 3. Key Issues of Collapse Simulation

### 3.1. Geometric Nonlinearity of Space Beam Element

When structures fail in fracture collapse, they often involve large deformations in the structural members. When structural displacement is getting large, the second-order effect resulting from this large displacement could result in an apparent variation in structural stiffness. An equilibrium equation must be established under the unknown displacement state. In the geometric nonlinear issue, Jagannathan et al. [15] proposed that the rigid body motion of a structure could cause a change in the structural internal force direction and virtual strain. This makes it difficult to determine the pure deformation of a structure and its internal force. To address this difficulty, a finite particle method introduces the concept of a “pathway element” to dissect the structural movement process into countless tiny time periods and then to utilize “virtual reverse motion” to eliminate the impact of the rigid body on internal force solution.Example  1: Geometric Large Deformation Analysis for Space Cantilever Curved Beam under the Short Time Impact Load. As shown in Figure 9, we used the space cantilever beam under a short time impact load as the calculation example for a large deformation analysis to validate the correctness of a beam element simulating the geometric nonlinear problems proposed in this paper. For comparison with the literature [16], we adopted the English unit. At the initial time, the curved beam was located in the x-y plane, and the corresponding central angle was 45°. The calculation parameters shown in Figure 9 were perpendicular to the plane of the curved beam. We applied a dead load P = 300 lb along the positive direction of the y-axis at the free end of the beam for a duration of 0.3 s. As for a large deformation analysis, a higher particle division density would be beneficial to obtaining more accurate calculation results. This calculation example utilized 20 elements for simulation of the cantilever beam with a time step of Δt=1.0×10-5s.Figure 9
Space cantilever curved beam.Using the finite particle method for analysis, we calculated the displacement change of the free end of the cantilever beam and compared these calculation results with results from a study by Chan [16]. As shown in Figures 10, 11, and 12, the cantilever curved beam would vibrate under the dead load at the free end. The calculated displacement changes of the free end for the cantilever curved beam along x, y, and z directions matched well with the calculation results from Chan [16] and therefore validated the reliability of applying a beam element in the finite particle method to analyze dynamic geometric nonlinear problems.Figure 10
Z-direction displacement of beam free end.Figure 11
X-direction displacement of beam free end.Figure 12
Y-direction displacement of beam free end.Example  2: Geometrical Nonlinearity Analysis of Spatial Elastic Framed Dome. The structural type and cross-sectional dimensions of the space hexagonal rigid frame model are shown in Figure 13, and the same cross-sectional form is adopted for all bars. Elastic modulus and shear modulus of material are E=2.06×1010N/m2and G=8.83×109N/m2, respectively. Ultimate load is Fcr=1.23×108N. The top center node A is subjected to the concentrated force F in the Z-direction. Through the displacement control of node A, the corresponding external load F can be obtained. The time step is Δt=1×10-3s, displacement increment is Δδ=1×10-3m, and duration of analysis is 3.5s. Each component of the structure is divided into n-segment beam elements, and three calculation cases n=2,n=6, and n=10 are considered, respectively. Using the finite particle method for geometrical nonlinearity analysis of the structure, the load-displacement curve at node A is shown in Figure 14. When n=2, the calculation error is too large; however, when n=6 and n=10, the calculation results tend to be stable. It is indicated that proper selection of a larger particle density can improve the accuracy of calculations while dealing with the large deformation problem of the structure. Park [17] used different nonlinear algorithms to analyze the example which belongs to the problem of geometrical large deformation. By the comparison analysis (Figure 14), the calculation results of the finite particle method are in good agreement with Park’s analysis result. It is shown that the geometrical nonlinearity analysis of the structure using the finite particle method is feasible.Figure 13
Spatial elastic framed dome.Figure 14
Load-displacement curve of node A.


### 3.2. Material Nonlinearity of Space Beam Element

In most situations, when a structure gradually approaches its limit load, the properties of the material gradually transform from linear elasticity to plasticity, and they behave according to the nonlinearity of the material. In the theory of finite particle method, the elastic plasticity of a structure is related only to the determination of an element’s internal force. Compared with the elastic analysis in the finite particle method, the elastic-plastic analysis does not need to alter the basic calculation frame, but it does need to consider a yielding condition and the loading-unloading flow for various internal force calculations while also calculating the internal force of the element. This paper adopted a plastic hinge concept to study the nonlinearity of the material and assumed a rapid development of plasticity after the cross-section entered into a plastic state. Therefore, this paper no longer discusses the detailed development process of plasticity on the cross-section of an element.We represented the yielding condition and the complete plastic condition of the element by a group of generalized stress that was the limit of the curved surface in the stress space could be represented by the yield functionΨe=0 of stress resultant and was the maximum value of Equation [12] out of the various components of the generalized stress relating to elastic behaviors. According to the von Mises criterion, the yielding conditions for the element to withstand axial force, moment, and torque can be expressed by(28)Ψe=NNY±MMY2+TTY2-1.Example  1: Material Nonlinear Analysis of Plane Cantilever Beam. As shown in Figure 15, the cantilever had a length of 5 in and a rectangular 0.1in×0.5in cross-section. We selected a double linear strengthening model for analysis. The free end of the cantilever was applied to point loadP along the y-direction, which gradually increased the load with a time step of Δt=1.0×10-5s and without considering the gravity effect. Compared with the reference literature [18, 19], this calculation example utilized five plane beam elements for simulation and analyzed the elastic-plastic behaviors of the cantilever beam under the application of vertical loading at the free end.Figure 15
Cantilever beam subject to joint load at the free end.As for the calculation example, Yang [18] and Tang [19], respectively, carried out elastic and elastic-plastic analysis using the finite element method. For comparison with the literature, this paper utilized the elastic calculation method and an elastic-plastic calculation method from the finite particle method theory to calculate the displacement of the cantilever beam free end while applying vertical loading. The cantilever structure in the example mainly showed bending, and the elastic modulus at the strengthening stage was used to calculate the stress increment after entering plasticity.Figure16 presents the displacement-loading curve of the cantilever free end that resulted from the elastic method. We found that the elastic calculation results from the finite particle method completely matched with the elastic solution by Tang [19] and validated the correctness of the beam element model. Figure 17 presents the displacement-loading curve of the cantilever free end that resulted from the elastic-plastic method. It was found that, with consideration for plasticity, the elastic-plastic calculation results from the finite particle method were slightly different from the elastic-plastic solutions by Yang [18] and Tang [19]. This was probably caused by the fact that the calculation simplification in the finite particle method did not accurately consider the plastic development process of the cross-section, but overall, the displacement-loading curve was well matched.Figure 16
Load-displacement of beam free end (elastic).Figure 17
Load-displacement of beam free end (elastic-plastic).Example  2: Material Nonlinearity of Spatial Right Angle Frame. The angle steel frame is shown in Figure 18, and the force F is applied at the midpoint C of the beam BD in the y-direction. The elastoplastic behavior of the frame under F action was analyzed using an ideal elastoplastic model. The structure is subject to the bending-torsion effect, and (28) is adopted as the yield condition of the component. The finite particle method is used for modeling analysis. The beams AB and BC are divided into 30 spatial beam elements and the time step is Δt=1×10-4s; the corresponding load increment is ΔP=ΔtN.Figure 18
Spatial right angle frame.Shi [20] and Turkalj [21] analyzed this example. The nondimensional load-displacement curve of the node B in the Y-direction of load F is shown in Figure 19. Both Turkalj and this paper use the plastic hinge method to simulate the yield of the structure, the yield process of the structure shows the characteristics of multiple segments, node D first enters plasticity, and then node A and node C yield one after the other. However, due to the different yield conditions, the results are slightly different, and the structural plasticity obtained by Turkalj has a dimensionless value of 4.977, while the corresponding value obtained by the finite particle method is 4.822, which is a difference of 3.5%. Although the finite particle method uses a plastic hinge model that does not consider the plastic development of the section, the calculation results show that the calculation accuracy can meet the requirements of engineering analysis.Figure 19
Nondimensional load-displacement curves of node B in the Y-direction.


### 3.3. Fracture Criterion of Space Steel Member

Space steel frame is a type of space beam system structure. The elements are under the joint actions of axial force, moment, and torque, and they belong to compression-bending or tension-bending members. During design, the calculation to determine the ultimate limit state of the compression-bending members includes strength and stability calculations. To select failure criterion, this paper utilized the ultimate limit state calculation method from the Chinese Steel Structure Design Code [22] and simultaneously considered strength control and stability control. Strength control was adopted for the ultimate loading capacity of the tension member and stability control was adopted for the ultimate loading capacity of compression member.


#### 3.3.1. Strength Failure Criterion of Member

For the two-way compression-bending or tension-bending members with moment applied in the two principal planes, the Chinese Steel Structure Design Code [22] requires that the formula for the cross-section strength calculation follows (29)NAn±MxγxWnx±MyγyWny≤f,where An represents the net sectional area of cross-section; Wnx and Wny represent the net sectional modulus of the cross-section to the x-axis and y-axis, respectively;frepresents the strength design value; and γx and γy are the sectional plastic development coefficients.


#### 3.3.2. Stability Failure Criterion of Members

The overall instability of two-way compression-bending members often is accompanied with torsion deformation, and this stability loading capacity is related to the proportion ofN, Mx, and My. The Chinese Steel Structure Design Code stipulates that, for compression-bending members with two main planes and biaxial symmetry of a solid web I-shape and a box section applied by moment, stability loading capacity is determined by(30)NφxA+βmxMxγxWx1-0.8N/NEx′+ηβtyMyφbyWy≤f,where φx represents the stability coefficient for the axial compression member corresponding to strong axis x-x and the weak axis y-y; φby represents the overall stability coefficient of the bending member subject to the uniform bending loading; Mx and My represent the maximum moment for the strong axis and weak axis within computation range of the member; NEx′ represents the calculation parameter, NEx′=π2EA/1.1λx2; Wx and Wy represent the sectional modulus for the strong axis and the weak axis; βmx represents the stability equivalent moment coefficient in plane; and βty represents the stability equivalent moment coefficient out of plane.


### 3.4. Failure Mode of Space Steel Members

Fracture behaviors of members involve a series of problems, such as the separation of particles, the redistribution of internal forces and masses, and the redefinition of topological relations of elements [23]. To simulate these behaviors in the analysis of the finite particle method, it is necessary to establish the fracture failure model. When the stress reaches the failure criterion of the material, the failure mode of the members may separate from the two ends of the particle simultaneously, and it is also possible to separate the members from one side of the particle. Particles maintain a dynamic balance under internal and external forces. Although axial forces at the both ends of the bar element are the same, the resultant forces of the particles at both sides of the bar element will not necessarily be the same. Thus, the resultant force of the particles at both sides of the bar element needs to be compared, and the fracture of the members will occur at the end of the greater resultant forces of the particles. Thus, fM and fN are defined to be the resultant forces of particles atM end andN end of the bar element, respectively (Figure 20(a)). If fM>fN, the fracture of the member occurs at particleM, the new particle M′ is generated, and the element is connected to the new particle M′ as shown in Figure 20(b). If fM<fN, the fracture of the member occurs at particle N, the new particle N′ is generated, and the element is connected to the new particle N′ as shown in Figure 20(c). If fM=fN, the fracture of the member occurs at both ends of the member, then two new particles M′and N′ are generated simultaneously, and the element is connected to the new particles M′ and N′, Finally, the element is completely separated from the structure as shown in Figure 20(d).Figure 20
Failure modes of steel members.
(a)
Original truss (b)
Fracture occurs at particleM (c)
Fracture occurs at particleN (d)
Fracture occurs at both endsExample  3. Space Steel Frame Collapse Analysis under Impact Loading. For this calculation example, we selected a five-story frame structural model [24] with a story height of 3 m and an overall height of 15 m as shown in Figure 21. We used this example to analyze the fracture behaviors of the structure and to validate the effectiveness of the failure criterion and failure modes proposed in this paper. The structural members in Figure 21 all utilized thin-wall pipe and shared identical material properties; sectional area A=0.002m2, second moment of area I=2.37×10-6m2, material density ρ=7.85×103kg/m3, the elastic modulus of material E=2.06×1011N/m2, and material strength design value fy=235MPa. At the beam-column joint along the left side of each story, we applied the point load along y-direction and in the form of slop loading. The loading time was 1.5 s, and the loading peak value Fn of each story was 1000N, 2000N, 3000N, 4000N, and 5000N.Figure 21
Five-story space steel frame.The unit time step during the analysis wasΔh=1×10-4s, and the analysis time was 5 s. Considering a gravity effect, we adopted an ideal elastic-plastic model, determined the yield function according to (28), and adopted the failure criterion and failure modes defined previously. The established finite particle model is shown in Figure 22(a).Figure 22
Collapse process of space steel frame.
(a)
t = 0.0 s (b)
t = 1.0 s (c)
t = 1.5 s (d)
t = 2.0 s (e)
t = 3.5 s (f)
t = 4.0 sOur analysis tracked the structural configurations at the various times under the application of a gradually increasing load as shown in Figure18. We listed the number of partially fractured elements, fracture time, and fracture configuration. Before 1.41s, the structure was predominated by plastic deformation, as shown in Figure 22(b); between 1.41s and 1.43s, the first batch of members successively fractured. At first, the bottom column member failed by compression-bending instability; subsequently, strength failure of the beams and columns occurred at the second story, as shown in Figure 22(c). Between 1.7 s and 1.86 s, the second batch members failed. The structure dropped as a result of the gravity effect after losing bottom vertical support, and the beam-column member at the structural middle collided with the ground and failed by strength loss, as shown in Figures 22(c), 22(e), and 22(f). In summary, we applied a gradually increasing load along the structural height, causing the structural bottom to be fractured resulting in the loss of upper structural vertical support and collapse caused by the gravity effect.


## 3.1. Geometric Nonlinearity of Space Beam Element

When structures fail in fracture collapse, they often involve large deformations in the structural members. When structural displacement is getting large, the second-order effect resulting from this large displacement could result in an apparent variation in structural stiffness. An equilibrium equation must be established under the unknown displacement state. In the geometric nonlinear issue, Jagannathan et al. [15] proposed that the rigid body motion of a structure could cause a change in the structural internal force direction and virtual strain. This makes it difficult to determine the pure deformation of a structure and its internal force. To address this difficulty, a finite particle method introduces the concept of a “pathway element” to dissect the structural movement process into countless tiny time periods and then to utilize “virtual reverse motion” to eliminate the impact of the rigid body on internal force solution.Example  1: Geometric Large Deformation Analysis for Space Cantilever Curved Beam under the Short Time Impact Load. As shown in Figure 9, we used the space cantilever beam under a short time impact load as the calculation example for a large deformation analysis to validate the correctness of a beam element simulating the geometric nonlinear problems proposed in this paper. For comparison with the literature [16], we adopted the English unit. At the initial time, the curved beam was located in the x-y plane, and the corresponding central angle was 45°. The calculation parameters shown in Figure 9 were perpendicular to the plane of the curved beam. We applied a dead load P = 300 lb along the positive direction of the y-axis at the free end of the beam for a duration of 0.3 s. As for a large deformation analysis, a higher particle division density would be beneficial to obtaining more accurate calculation results. This calculation example utilized 20 elements for simulation of the cantilever beam with a time step of Δt=1.0×10-5s.Figure 9
Space cantilever curved beam.Using the finite particle method for analysis, we calculated the displacement change of the free end of the cantilever beam and compared these calculation results with results from a study by Chan [16]. As shown in Figures 10, 11, and 12, the cantilever curved beam would vibrate under the dead load at the free end. The calculated displacement changes of the free end for the cantilever curved beam along x, y, and z directions matched well with the calculation results from Chan [16] and therefore validated the reliability of applying a beam element in the finite particle method to analyze dynamic geometric nonlinear problems.Figure 10
Z-direction displacement of beam free end.Figure 11
X-direction displacement of beam free end.Figure 12
Y-direction displacement of beam free end.Example  2: Geometrical Nonlinearity Analysis of Spatial Elastic Framed Dome. The structural type and cross-sectional dimensions of the space hexagonal rigid frame model are shown in Figure 13, and the same cross-sectional form is adopted for all bars. Elastic modulus and shear modulus of material are E=2.06×1010N/m2and G=8.83×109N/m2, respectively. Ultimate load is Fcr=1.23×108N. The top center node A is subjected to the concentrated force F in the Z-direction. Through the displacement control of node A, the corresponding external load F can be obtained. The time step is Δt=1×10-3s, displacement increment is Δδ=1×10-3m, and duration of analysis is 3.5s. Each component of the structure is divided into n-segment beam elements, and three calculation cases n=2,n=6, and n=10 are considered, respectively. Using the finite particle method for geometrical nonlinearity analysis of the structure, the load-displacement curve at node A is shown in Figure 14. When n=2, the calculation error is too large; however, when n=6 and n=10, the calculation results tend to be stable. It is indicated that proper selection of a larger particle density can improve the accuracy of calculations while dealing with the large deformation problem of the structure. Park [17] used different nonlinear algorithms to analyze the example which belongs to the problem of geometrical large deformation. By the comparison analysis (Figure 14), the calculation results of the finite particle method are in good agreement with Park’s analysis result. It is shown that the geometrical nonlinearity analysis of the structure using the finite particle method is feasible.Figure 13
Spatial elastic framed dome.Figure 14
Load-displacement curve of node A.


## 3.2. Material Nonlinearity of Space Beam Element

In most situations, when a structure gradually approaches its limit load, the properties of the material gradually transform from linear elasticity to plasticity, and they behave according to the nonlinearity of the material. In the theory of finite particle method, the elastic plasticity of a structure is related only to the determination of an element’s internal force. Compared with the elastic analysis in the finite particle method, the elastic-plastic analysis does not need to alter the basic calculation frame, but it does need to consider a yielding condition and the loading-unloading flow for various internal force calculations while also calculating the internal force of the element. This paper adopted a plastic hinge concept to study the nonlinearity of the material and assumed a rapid development of plasticity after the cross-section entered into a plastic state. Therefore, this paper no longer discusses the detailed development process of plasticity on the cross-section of an element.We represented the yielding condition and the complete plastic condition of the element by a group of generalized stress that was the limit of the curved surface in the stress space could be represented by the yield functionΨe=0 of stress resultant and was the maximum value of Equation [12] out of the various components of the generalized stress relating to elastic behaviors. According to the von Mises criterion, the yielding conditions for the element to withstand axial force, moment, and torque can be expressed by(28)Ψe=NNY±MMY2+TTY2-1.Example  1: Material Nonlinear Analysis of Plane Cantilever Beam. As shown in Figure 15, the cantilever had a length of 5 in and a rectangular 0.1in×0.5in cross-section. We selected a double linear strengthening model for analysis. The free end of the cantilever was applied to point loadP along the y-direction, which gradually increased the load with a time step of Δt=1.0×10-5s and without considering the gravity effect. Compared with the reference literature [18, 19], this calculation example utilized five plane beam elements for simulation and analyzed the elastic-plastic behaviors of the cantilever beam under the application of vertical loading at the free end.Figure 15
Cantilever beam subject to joint load at the free end.As for the calculation example, Yang [18] and Tang [19], respectively, carried out elastic and elastic-plastic analysis using the finite element method. For comparison with the literature, this paper utilized the elastic calculation method and an elastic-plastic calculation method from the finite particle method theory to calculate the displacement of the cantilever beam free end while applying vertical loading. The cantilever structure in the example mainly showed bending, and the elastic modulus at the strengthening stage was used to calculate the stress increment after entering plasticity.Figure16 presents the displacement-loading curve of the cantilever free end that resulted from the elastic method. We found that the elastic calculation results from the finite particle method completely matched with the elastic solution by Tang [19] and validated the correctness of the beam element model. Figure 17 presents the displacement-loading curve of the cantilever free end that resulted from the elastic-plastic method. It was found that, with consideration for plasticity, the elastic-plastic calculation results from the finite particle method were slightly different from the elastic-plastic solutions by Yang [18] and Tang [19]. This was probably caused by the fact that the calculation simplification in the finite particle method did not accurately consider the plastic development process of the cross-section, but overall, the displacement-loading curve was well matched.Figure 16
Load-displacement of beam free end (elastic).Figure 17
Load-displacement of beam free end (elastic-plastic).Example  2: Material Nonlinearity of Spatial Right Angle Frame. The angle steel frame is shown in Figure 18, and the force F is applied at the midpoint C of the beam BD in the y-direction. The elastoplastic behavior of the frame under F action was analyzed using an ideal elastoplastic model. The structure is subject to the bending-torsion effect, and (28) is adopted as the yield condition of the component. The finite particle method is used for modeling analysis. The beams AB and BC are divided into 30 spatial beam elements and the time step is Δt=1×10-4s; the corresponding load increment is ΔP=ΔtN.Figure 18
Spatial right angle frame.Shi [20] and Turkalj [21] analyzed this example. The nondimensional load-displacement curve of the node B in the Y-direction of load F is shown in Figure 19. Both Turkalj and this paper use the plastic hinge method to simulate the yield of the structure, the yield process of the structure shows the characteristics of multiple segments, node D first enters plasticity, and then node A and node C yield one after the other. However, due to the different yield conditions, the results are slightly different, and the structural plasticity obtained by Turkalj has a dimensionless value of 4.977, while the corresponding value obtained by the finite particle method is 4.822, which is a difference of 3.5%. Although the finite particle method uses a plastic hinge model that does not consider the plastic development of the section, the calculation results show that the calculation accuracy can meet the requirements of engineering analysis.Figure 19
Nondimensional load-displacement curves of node B in the Y-direction.


## 3.3. Fracture Criterion of Space Steel Member

Space steel frame is a type of space beam system structure. The elements are under the joint actions of axial force, moment, and torque, and they belong to compression-bending or tension-bending members. During design, the calculation to determine the ultimate limit state of the compression-bending members includes strength and stability calculations. To select failure criterion, this paper utilized the ultimate limit state calculation method from the Chinese Steel Structure Design Code [22] and simultaneously considered strength control and stability control. Strength control was adopted for the ultimate loading capacity of the tension member and stability control was adopted for the ultimate loading capacity of compression member.


### 3.3.1. Strength Failure Criterion of Member

For the two-way compression-bending or tension-bending members with moment applied in the two principal planes, the Chinese Steel Structure Design Code [22] requires that the formula for the cross-section strength calculation follows (29)NAn±MxγxWnx±MyγyWny≤f,where An represents the net sectional area of cross-section; Wnx and Wny represent the net sectional modulus of the cross-section to the x-axis and y-axis, respectively;frepresents the strength design value; and γx and γy are the sectional plastic development coefficients.


### 3.3.2. Stability Failure Criterion of Members

The overall instability of two-way compression-bending members often is accompanied with torsion deformation, and this stability loading capacity is related to the proportion ofN, Mx, and My. The Chinese Steel Structure Design Code stipulates that, for compression-bending members with two main planes and biaxial symmetry of a solid web I-shape and a box section applied by moment, stability loading capacity is determined by(30)NφxA+βmxMxγxWx1-0.8N/NEx′+ηβtyMyφbyWy≤f,where φx represents the stability coefficient for the axial compression member corresponding to strong axis x-x and the weak axis y-y; φby represents the overall stability coefficient of the bending member subject to the uniform bending loading; Mx and My represent the maximum moment for the strong axis and weak axis within computation range of the member; NEx′ represents the calculation parameter, NEx′=π2EA/1.1λx2; Wx and Wy represent the sectional modulus for the strong axis and the weak axis; βmx represents the stability equivalent moment coefficient in plane; and βty represents the stability equivalent moment coefficient out of plane.


## 3.3.1. Strength Failure Criterion of Member

For the two-way compression-bending or tension-bending members with moment applied in the two principal planes, the Chinese Steel Structure Design Code [22] requires that the formula for the cross-section strength calculation follows (29)NAn±MxγxWnx±MyγyWny≤f,where An represents the net sectional area of cross-section; Wnx and Wny represent the net sectional modulus of the cross-section to the x-axis and y-axis, respectively;frepresents the strength design value; and γx and γy are the sectional plastic development coefficients.


## 3.3.2. Stability Failure Criterion of Members

The overall instability of two-way compression-bending members often is accompanied with torsion deformation, and this stability loading capacity is related to the proportion ofN, Mx, and My. The Chinese Steel Structure Design Code stipulates that, for compression-bending members with two main planes and biaxial symmetry of a solid web I-shape and a box section applied by moment, stability loading capacity is determined by(30)NφxA+βmxMxγxWx1-0.8N/NEx′+ηβtyMyφbyWy≤f,where φx represents the stability coefficient for the axial compression member corresponding to strong axis x-x and the weak axis y-y; φby represents the overall stability coefficient of the bending member subject to the uniform bending loading; Mx and My represent the maximum moment for the strong axis and weak axis within computation range of the member; NEx′ represents the calculation parameter, NEx′=π2EA/1.1λx2; Wx and Wy represent the sectional modulus for the strong axis and the weak axis; βmx represents the stability equivalent moment coefficient in plane; and βty represents the stability equivalent moment coefficient out of plane.


## 3.4. Failure Mode of Space Steel Members

Fracture behaviors of members involve a series of problems, such as the separation of particles, the redistribution of internal forces and masses, and the redefinition of topological relations of elements [23]. To simulate these behaviors in the analysis of the finite particle method, it is necessary to establish the fracture failure model. When the stress reaches the failure criterion of the material, the failure mode of the members may separate from the two ends of the particle simultaneously, and it is also possible to separate the members from one side of the particle. Particles maintain a dynamic balance under internal and external forces. Although axial forces at the both ends of the bar element are the same, the resultant forces of the particles at both sides of the bar element will not necessarily be the same. Thus, the resultant force of the particles at both sides of the bar element needs to be compared, and the fracture of the members will occur at the end of the greater resultant forces of the particles. Thus, fM and fN are defined to be the resultant forces of particles atM end andN end of the bar element, respectively (Figure 20(a)). If fM>fN, the fracture of the member occurs at particleM, the new particle M′ is generated, and the element is connected to the new particle M′ as shown in Figure 20(b). If fM<fN, the fracture of the member occurs at particle N, the new particle N′ is generated, and the element is connected to the new particle N′ as shown in Figure 20(c). If fM=fN, the fracture of the member occurs at both ends of the member, then two new particles M′and N′ are generated simultaneously, and the element is connected to the new particles M′ and N′, Finally, the element is completely separated from the structure as shown in Figure 20(d).Figure 20
Failure modes of steel members.
(a)
Original truss (b)
Fracture occurs at particleM (c)
Fracture occurs at particleN (d)
Fracture occurs at both endsExample  3. Space Steel Frame Collapse Analysis under Impact Loading. For this calculation example, we selected a five-story frame structural model [24] with a story height of 3 m and an overall height of 15 m as shown in Figure 21. We used this example to analyze the fracture behaviors of the structure and to validate the effectiveness of the failure criterion and failure modes proposed in this paper. The structural members in Figure 21 all utilized thin-wall pipe and shared identical material properties; sectional area A=0.002m2, second moment of area I=2.37×10-6m2, material density ρ=7.85×103kg/m3, the elastic modulus of material E=2.06×1011N/m2, and material strength design value fy=235MPa. At the beam-column joint along the left side of each story, we applied the point load along y-direction and in the form of slop loading. The loading time was 1.5 s, and the loading peak value Fn of each story was 1000N, 2000N, 3000N, 4000N, and 5000N.Figure 21
Five-story space steel frame.The unit time step during the analysis wasΔh=1×10-4s, and the analysis time was 5 s. Considering a gravity effect, we adopted an ideal elastic-plastic model, determined the yield function according to (28), and adopted the failure criterion and failure modes defined previously. The established finite particle model is shown in Figure 22(a).Figure 22
Collapse process of space steel frame.
(a)
t = 0.0 s (b)
t = 1.0 s (c)
t = 1.5 s (d)
t = 2.0 s (e)
t = 3.5 s (f)
t = 4.0 sOur analysis tracked the structural configurations at the various times under the application of a gradually increasing load as shown in Figure18. We listed the number of partially fractured elements, fracture time, and fracture configuration. Before 1.41s, the structure was predominated by plastic deformation, as shown in Figure 22(b); between 1.41s and 1.43s, the first batch of members successively fractured. At first, the bottom column member failed by compression-bending instability; subsequently, strength failure of the beams and columns occurred at the second story, as shown in Figure 22(c). Between 1.7 s and 1.86 s, the second batch members failed. The structure dropped as a result of the gravity effect after losing bottom vertical support, and the beam-column member at the structural middle collided with the ground and failed by strength loss, as shown in Figures 22(c), 22(e), and 22(f). In summary, we applied a gradually increasing load along the structural height, causing the structural bottom to be fractured resulting in the loss of upper structural vertical support and collapse caused by the gravity effect.


## 4. Space Frame Collapse Simulation under Earthquake Excitation

### 4.1. Model Establishment and Validation

The five-story steel frame structure in example 3 was the analysis object. We used ANSYS to establish the finite element model. All members utilized the BEAM188 element for the simulation. We utilized MATLAB to establish the finite particle method. All members utilized beam element for simulation while simultaneously considering axial force, moment, and torsion moment. Each particle considered three directions of linear displacement and three directions of angle displacement. We assumed the ground was rigid and that the steel frame and the ground had a rigid connection. The models established by these two methods are shown in Figures23 and 24.Figure 23
Finite element model (ANSYS).Figure 24
Finite particle model (MATLAB).To validate the accuracy of the established model, we used the finite element method and finite particle method for static analysis of the steel frame. We considered two working conditions: working condition 1, in which we applied gravity load with gravity accelerationg=9.8m/s2; and working condition 2, in which we applied point load F=500N along the x-axis positive direction at the four beam-column joints of the top story. In the finite particle method, we used a time step of 1×10-4s and compared the calculation results from the two methods.We extracted the vertical reaction force from the four supporting points under the gravity action and the average vertical displacement along various height joints. As shown in Figure25 and Table 1, we observed that the calculation results from the two methods were very close with very small error. We extracted the average displacement in x-direction at various height joints under the action of joint load in x-direction. As shown in Figure 26 and Table 2, we found that the calculation results from these two methods matched well with a maximum error of −1.44%. Through the static analysis results from the finite particle method and finite element method, we found that the steel frame calculation analysis model based on these finite particle methods was reliable.Table 1
Vertical reaction force of supportings and vertical displacement of nodes under earthquake action.
Analysis method Vertical reaction force of various supportings (N) Vertical displacement of various height nodes (mm) #1 #4 #7 #10 6 m high 9 m high 12m high 15 m high FEM 11309 11309 11309 11309 −0.057 −0.076 −0.087 −0.092  FPM 11308.7 11308.7 11308.7 11308.7 −0.057 −0.075 −0.087 −0.092  Relative error 0.00% 0.00% 0.00% 0.00% 0.00% 0.13% 0.11% 0.00%Table 2
Supporting reaction force under horizontal concentrated load.
Analysis method Reaction force of supporting #1 Reaction force of supporting #4 f x N f y N f z N f x N f y N f z N FEM −500.00 0.00 −4345.6 −500.00 0.00 4345.6  FPM −494.1 0.00 −4310.2 −492.8 0.00 4310.1  Relative error −1.18% 0.00% −0.81% −1.44% 0.00% −0.82%Figure 25
Vertical displacement of various heights under gravity.Figure 26
X-direction displacement of various heights under horizontal point loading.


### 4.2. Earthquake Response Analysis of Space Steel Frame

To further validate the reliability of earthquake response analysis using the finite particle method, we compared the dynamic analysis results from the finite particle method with its counterpart from the finite element analysis. We selected the prior 10-s acceleration time history of the EI centro wave for our earthquake dynamic response analysis; the acceleration time history curve is shown in Figure27. We adjusted the peak acceleration to 1.0 m/s2. We selected the time step in ANSYS as 0.02 s, according to the time interval of the earthquake wave; the time step was 1×10-4s in the MATLAB analysis and was input in x-direction. We selected a Rayleigh damping coefficient for damping, and the structural damping ratio was 0.03. We considered only mass damping and determined it to be α=0.6166 through calculation. The time history of the top displacement is shown in Figure 28. The shear time history curve for the left end of element 22 is shown in Figure 29.Figure 27
EI Centro acceleration time history curve.Figure 28
X-direction time history curve of vertex.Figure 29
Z-direction shear time history curve for the left end of element 22.We extracted the maximum displacement values of various height joints (3 m, 6 m, 9 m, 12 m, and 15 m) in the whole time history for the steel frame, which were calculated from the two methods for comparison, as shown in Figure30. With an increase in height, the maximum displacements calculated from the two methods kept increasing and basically remained consistent with small error. This meant that the structural overall displacement time history calculated from these two methods matched well.Figure 30
Maximum displacement values of various heights.We extracted the time history change of the structural internal force that was calculated from the two methods. Figure29 presents the z-direction shear time history curve at the left end of element 22. We observed that the overall variation tendency of the element end shear calculated from the two methods basically remained consistent. Table 3 lists the vertical supporting reaction extreme value of supporting 4 and the extreme value of the z-direction shear at the left end of element 22. We found that the internal force extreme values resulting from these two methods were quite close with a maximum error of 4.8%. Therefore, we can assume that earthquake responses based on the finite particle method and finite element method are consistent, and based on this observation, we can conduct the subsequent earthquake collapse simulation analysis.Table 3
Partial internal force extreme values.
Analysis method Reaction force of supporting #4 (N) Z-direction shear force for the left end of element 22 (N) FEM 9983.00 −11387.50 3280.16 −3209.82  FPM 10465.56 −11251.20 3268.76 −3336.80  Relative error 4.80% −1.20% −0.35% 3.96%


### 4.3. Space Steel Frame Collapse Simulation under Earthquake Excitation

With consideration for gravity effect, we determined the mass damping to be 0.6166, according to two-order circular frequency. The strength design value of the material wasfy=235MPa. The failure of the members was jointly controlled by strength failure criterion and stability failure criterion, which we determined by (28), (29), and (30). We determined the stability coefficient, moment coefficient, and other related parameters of each member according to the corresponding requirements from the Chinese Steel Structure Design Code. We selected the prior 10-s acceleration of the EI Centro wave for time history analysis. To determine the collapse critical peak ground motion (PGA) of the structure and discuss the collapse behaviors under various earthquake PGA actions, we adjusted the earthquake PGA to 4.5m/s2, 5.0 m/s2, 6.0 m/s2, 8.0 m/s2, and 12.0 m/s2. We summarized the members that failed first, failure time, and failure duration from fracture start to complete collapse. The results are given in Table 4.Table 4
Structural collapse time under various PGAs.
Case Case 1 Case 2 Case 3 Case 4 Case 5 PGA (m·s−2) 4.5 5.0 6.0 8.0 12.0  First failure time (s) No failure 2.098 2.069 2.045 1.555  Member number No failure 22 22 22 22  Duration from fracture beginning to complete collapse of structure (s) No failure 2.310 2.276 2.025 1.983The critical earthquake PGA of the steel frame collapse was between 4.5m/s2 and 6.0 m/s2, and with an increase of PGA, the collapse was more dramatic, and the collapse duration was shorter. To better capture the details of the steel frame collapse process, we selected the collapse process of the steel frame at the earthquake PGA of 6.0 m/s2 for analysis.The collisions between fractured members, member, and structure were not taken into consideration. The steel frame began to fail at 2.069 s, and it completely collapsed within 4.347 s with a duration of 2.276 s.We used the finite particle method to simulate the whole collapse process. We captured the structural configuration of each time and selected structural deformation diagrams at several typical times, as shown in Figure31.Figure 31
Collapse diagram of steel frame at various times.
(a)
t = 0.0s (b)
t = 2.1s (c)
t = 2.6s (d)
t = 3.4s (e)
t = 3.8s (f)
t = 4.0sFigure31(b) shows that the two supportings at the bottom story failed by separation between column and supporting; that is, the transverse beam connected with that column fractured and separated from the upper structure. Around 2.6 s, two additional supportings failed. The upper structure lost its supporting connection and dropped as a result of the gravity effect, as shown in Figure 31(c). Around 3.4 s, the partial members departed from the main structure after it fractured because of a collision with the ground, and subsequently, the upper structure gradually dropped, collapsed, and broke (Figure 31(d)). The force-transferring path of the frame was simple with less redundant constraint. After the bottom support member failed, the internal force redistribution caused a chain failure of adjacent members and finally resulted in the continuous collapse of the upper structure, as shown in Figures 31(e) and 31(f).To thoroughly learn the collapse process of the steel frame, we extracted some typical structural internal variation. Figure32(a) shows the axial force of element 2, moment along two directions of the principal axis at the element end, and variations in the time history of element torsion. Note that the axial force, moment, and torsion all were determined under the local coordinate system. The lower part of element 2 connected with supporting 4, and the upper part of element 2 connected with the main structure. Around 2.1 s, that element first failed because of a strength failure at which point it separated from the supporting; however, the upper part remained connected to the main structure. At first, the axial force of element 2 continued to vary with changes in earthquake acceleration. Element 2 fractured around 2.1 s, the bottom separated from the supporting, and the element axial force decayed. That element fractured into a free end and generated a new particle, and the other end remained connected to the upper structure. Around 3.1 s, this element separated from the main structure, but it remained connected to element 2. Element 21 and element 24 collided with the ground and caused a large fluctuation in element internal force. At this time, connected element 2 was subject to internal force redistribution, causing the axial force to increase dramatically and then decay gradually. The whole variation process is shown in Figure 32(a). Figure 32(b) presents the moment of the element 2 end along the two directions of the principle axis. The overall tendency of change was similar to that of axial force. Before 2.1 s, the moment kept changing with earthquake acceleration, and the element fractured at 2.1 s, after which the moment of element ends gradually decayed. Around 3.1 s, the element was subjected to a collision from adjacent elements, and accordingly resulted in internal force redistribution and subsequent decay. Figure 32(c) presents the change in element torsion. Before fracture, element 2 was mainly subject to bending, and the torsion was relatively small at the beginning of loading. Around 3.1 s, because of the distribution of internal force, the element was loaded, and the torsion had a larger fluctuation; after that, the torsion gradually decreased due to damping (Figure 32(d)).Figure 32
Internal force time history curve of element 2.
(a)
Element axial force (b)
In-plane moment ofx^-z^ plane for element end (c)
In-plane moment ofx^-y^ plane for element end (d)
Element torqueWe extracted the vertical reaction change of the supporting as shown in Figure33. The supporting reaction from 0 s to 2.1 s exhibited regular changes along with earthquake acceleration. Supporting 4 failed around 2.07 s and caused the separation of element 2 from supporting 4, which subsequently was unloaded. Around 2.17 s, supporting 1 failed around 2.17 s and caused a separation of element 1 from supporting 1, which subsequently was unloaded.Figure 33
Vertical reaction time history curves for supporting 1 and 4.
(a)
supporting 1 (b)
supporting 4Therefore, the bottom beam-column members of the steel frame formed weak parts as a result of deficiencies in strength and stability loading capacity. These weak parts under the earthquake excitation first failed and subsequently resulted in the loss of effective vertical support to the upper main structure. This was the main reason for the continuous collapse of the structure. Structural collapse failure is a complex behavior, but the finite particle method is able to perform simulation analysis by selecting appropriate calculation parameters, calculation model, and failure criterion. This method can be used to thoroughly discover the collapse mechanism during structural failure process and provides a new concept and method to investigate the collapse of structures as a result of earthquake excitation.


## 4.1. Model Establishment and Validation

The five-story steel frame structure in example 3 was the analysis object. We used ANSYS to establish the finite element model. All members utilized the BEAM188 element for the simulation. We utilized MATLAB to establish the finite particle method. All members utilized beam element for simulation while simultaneously considering axial force, moment, and torsion moment. Each particle considered three directions of linear displacement and three directions of angle displacement. We assumed the ground was rigid and that the steel frame and the ground had a rigid connection. The models established by these two methods are shown in Figures23 and 24.Figure 23
Finite element model (ANSYS).Figure 24
Finite particle model (MATLAB).To validate the accuracy of the established model, we used the finite element method and finite particle method for static analysis of the steel frame. We considered two working conditions: working condition 1, in which we applied gravity load with gravity accelerationg=9.8m/s2; and working condition 2, in which we applied point load F=500N along the x-axis positive direction at the four beam-column joints of the top story. In the finite particle method, we used a time step of 1×10-4s and compared the calculation results from the two methods.We extracted the vertical reaction force from the four supporting points under the gravity action and the average vertical displacement along various height joints. As shown in Figure25 and Table 1, we observed that the calculation results from the two methods were very close with very small error. We extracted the average displacement in x-direction at various height joints under the action of joint load in x-direction. As shown in Figure 26 and Table 2, we found that the calculation results from these two methods matched well with a maximum error of −1.44%. Through the static analysis results from the finite particle method and finite element method, we found that the steel frame calculation analysis model based on these finite particle methods was reliable.Table 1
Vertical reaction force of supportings and vertical displacement of nodes under earthquake action.
Analysis method Vertical reaction force of various supportings (N) Vertical displacement of various height nodes (mm) #1 #4 #7 #10 6 m high 9 m high 12m high 15 m high FEM 11309 11309 11309 11309 −0.057 −0.076 −0.087 −0.092  FPM 11308.7 11308.7 11308.7 11308.7 −0.057 −0.075 −0.087 −0.092  Relative error 0.00% 0.00% 0.00% 0.00% 0.00% 0.13% 0.11% 0.00%Table 2
Supporting reaction force under horizontal concentrated load.
Analysis method Reaction force of supporting #1 Reaction force of supporting #4 f x N f y N f z N f x N f y N f z N FEM −500.00 0.00 −4345.6 −500.00 0.00 4345.6  FPM −494.1 0.00 −4310.2 −492.8 0.00 4310.1  Relative error −1.18% 0.00% −0.81% −1.44% 0.00% −0.82%Figure 25
Vertical displacement of various heights under gravity.Figure 26
X-direction displacement of various heights under horizontal point loading.


## 4.2. Earthquake Response Analysis of Space Steel Frame

To further validate the reliability of earthquake response analysis using the finite particle method, we compared the dynamic analysis results from the finite particle method with its counterpart from the finite element analysis. We selected the prior 10-s acceleration time history of the EI centro wave for our earthquake dynamic response analysis; the acceleration time history curve is shown in Figure27. We adjusted the peak acceleration to 1.0 m/s2. We selected the time step in ANSYS as 0.02 s, according to the time interval of the earthquake wave; the time step was 1×10-4s in the MATLAB analysis and was input in x-direction. We selected a Rayleigh damping coefficient for damping, and the structural damping ratio was 0.03. We considered only mass damping and determined it to be α=0.6166 through calculation. The time history of the top displacement is shown in Figure 28. The shear time history curve for the left end of element 22 is shown in Figure 29.Figure 27
EI Centro acceleration time history curve.Figure 28
X-direction time history curve of vertex.Figure 29
Z-direction shear time history curve for the left end of element 22.We extracted the maximum displacement values of various height joints (3 m, 6 m, 9 m, 12 m, and 15 m) in the whole time history for the steel frame, which were calculated from the two methods for comparison, as shown in Figure30. With an increase in height, the maximum displacements calculated from the two methods kept increasing and basically remained consistent with small error. This meant that the structural overall displacement time history calculated from these two methods matched well.Figure 30
Maximum displacement values of various heights.We extracted the time history change of the structural internal force that was calculated from the two methods. Figure29 presents the z-direction shear time history curve at the left end of element 22. We observed that the overall variation tendency of the element end shear calculated from the two methods basically remained consistent. Table 3 lists the vertical supporting reaction extreme value of supporting 4 and the extreme value of the z-direction shear at the left end of element 22. We found that the internal force extreme values resulting from these two methods were quite close with a maximum error of 4.8%. Therefore, we can assume that earthquake responses based on the finite particle method and finite element method are consistent, and based on this observation, we can conduct the subsequent earthquake collapse simulation analysis.Table 3
Partial internal force extreme values.
Analysis method Reaction force of supporting #4 (N) Z-direction shear force for the left end of element 22 (N) FEM 9983.00 −11387.50 3280.16 −3209.82  FPM 10465.56 −11251.20 3268.76 −3336.80  Relative error 4.80% −1.20% −0.35% 3.96%


## 4.3. Space Steel Frame Collapse Simulation under Earthquake Excitation

With consideration for gravity effect, we determined the mass damping to be 0.6166, according to two-order circular frequency. The strength design value of the material wasfy=235MPa. The failure of the members was jointly controlled by strength failure criterion and stability failure criterion, which we determined by (28), (29), and (30). We determined the stability coefficient, moment coefficient, and other related parameters of each member according to the corresponding requirements from the Chinese Steel Structure Design Code. We selected the prior 10-s acceleration of the EI Centro wave for time history analysis. To determine the collapse critical peak ground motion (PGA) of the structure and discuss the collapse behaviors under various earthquake PGA actions, we adjusted the earthquake PGA to 4.5m/s2, 5.0 m/s2, 6.0 m/s2, 8.0 m/s2, and 12.0 m/s2. We summarized the members that failed first, failure time, and failure duration from fracture start to complete collapse. The results are given in Table 4.Table 4
Structural collapse time under various PGAs.
Case Case 1 Case 2 Case 3 Case 4 Case 5 PGA (m·s−2) 4.5 5.0 6.0 8.0 12.0  First failure time (s) No failure 2.098 2.069 2.045 1.555  Member number No failure 22 22 22 22  Duration from fracture beginning to complete collapse of structure (s) No failure 2.310 2.276 2.025 1.983The critical earthquake PGA of the steel frame collapse was between 4.5m/s2 and 6.0 m/s2, and with an increase of PGA, the collapse was more dramatic, and the collapse duration was shorter. To better capture the details of the steel frame collapse process, we selected the collapse process of the steel frame at the earthquake PGA of 6.0 m/s2 for analysis.The collisions between fractured members, member, and structure were not taken into consideration. The steel frame began to fail at 2.069 s, and it completely collapsed within 4.347 s with a duration of 2.276 s.We used the finite particle method to simulate the whole collapse process. We captured the structural configuration of each time and selected structural deformation diagrams at several typical times, as shown in Figure31.Figure 31
Collapse diagram of steel frame at various times.
(a)
t = 0.0s (b)
t = 2.1s (c)
t = 2.6s (d)
t = 3.4s (e)
t = 3.8s (f)
t = 4.0sFigure31(b) shows that the two supportings at the bottom story failed by separation between column and supporting; that is, the transverse beam connected with that column fractured and separated from the upper structure. Around 2.6 s, two additional supportings failed. The upper structure lost its supporting connection and dropped as a result of the gravity effect, as shown in Figure 31(c). Around 3.4 s, the partial members departed from the main structure after it fractured because of a collision with the ground, and subsequently, the upper structure gradually dropped, collapsed, and broke (Figure 31(d)). The force-transferring path of the frame was simple with less redundant constraint. After the bottom support member failed, the internal force redistribution caused a chain failure of adjacent members and finally resulted in the continuous collapse of the upper structure, as shown in Figures 31(e) and 31(f).To thoroughly learn the collapse process of the steel frame, we extracted some typical structural internal variation. Figure32(a) shows the axial force of element 2, moment along two directions of the principal axis at the element end, and variations in the time history of element torsion. Note that the axial force, moment, and torsion all were determined under the local coordinate system. The lower part of element 2 connected with supporting 4, and the upper part of element 2 connected with the main structure. Around 2.1 s, that element first failed because of a strength failure at which point it separated from the supporting; however, the upper part remained connected to the main structure. At first, the axial force of element 2 continued to vary with changes in earthquake acceleration. Element 2 fractured around 2.1 s, the bottom separated from the supporting, and the element axial force decayed. That element fractured into a free end and generated a new particle, and the other end remained connected to the upper structure. Around 3.1 s, this element separated from the main structure, but it remained connected to element 2. Element 21 and element 24 collided with the ground and caused a large fluctuation in element internal force. At this time, connected element 2 was subject to internal force redistribution, causing the axial force to increase dramatically and then decay gradually. The whole variation process is shown in Figure 32(a). Figure 32(b) presents the moment of the element 2 end along the two directions of the principle axis. The overall tendency of change was similar to that of axial force. Before 2.1 s, the moment kept changing with earthquake acceleration, and the element fractured at 2.1 s, after which the moment of element ends gradually decayed. Around 3.1 s, the element was subjected to a collision from adjacent elements, and accordingly resulted in internal force redistribution and subsequent decay. Figure 32(c) presents the change in element torsion. Before fracture, element 2 was mainly subject to bending, and the torsion was relatively small at the beginning of loading. Around 3.1 s, because of the distribution of internal force, the element was loaded, and the torsion had a larger fluctuation; after that, the torsion gradually decreased due to damping (Figure 32(d)).Figure 32
Internal force time history curve of element 2.
(a)
Element axial force (b)
In-plane moment ofx^-z^ plane for element end (c)
In-plane moment ofx^-y^ plane for element end (d)
Element torqueWe extracted the vertical reaction change of the supporting as shown in Figure33. The supporting reaction from 0 s to 2.1 s exhibited regular changes along with earthquake acceleration. Supporting 4 failed around 2.07 s and caused the separation of element 2 from supporting 4, which subsequently was unloaded. Around 2.17 s, supporting 1 failed around 2.17 s and caused a separation of element 1 from supporting 1, which subsequently was unloaded.Figure 33
Vertical reaction time history curves for supporting 1 and 4.
(a)
supporting 1 (b)
supporting 4Therefore, the bottom beam-column members of the steel frame formed weak parts as a result of deficiencies in strength and stability loading capacity. These weak parts under the earthquake excitation first failed and subsequently resulted in the loss of effective vertical support to the upper main structure. This was the main reason for the continuous collapse of the structure. Structural collapse failure is a complex behavior, but the finite particle method is able to perform simulation analysis by selecting appropriate calculation parameters, calculation model, and failure criterion. This method can be used to thoroughly discover the collapse mechanism during structural failure process and provides a new concept and method to investigate the collapse of structures as a result of earthquake excitation.


## 5. Conclusions

This paper utilized the finite particle method as fundamental analysis means and MATLAB software as computation analysis tool. We carried out simulation and analysis on the collapse process of space steel structures subjected to earthquake excitation.On the basis of the results from our analysis, the following conclusions can be drawn:(i)
We derived the calculation formula and the calculation steps to solve static and dynamic problems applying the finite particle method used in space beam system structural analysis and achieve using MATLAB programming.(ii)
We selected two calculation examples and utilized the finite particle method to carry out comparison between geometric nonlinearity and material nonlinearity. This comparison validated the correctness and reliability of the finite particle method in our structural nonlinear analysis; established the fracture judge criterion and failure modes for the space beam system structural members; simulated the collapse failure process of the space steel frame structure under a horizontal impact loading; and validated the reasonability and effectiveness of the failure criterion, failure mode, and fracture calculation method proposed.(iii)
We comprehensively considered the geometric failure, material failure, and member fracture behaviors of the structure, and we conducted numerical simulations of the structural collapse process under earthquake excitation. The results showed that the larger the earthquake PGA is, the earlier the structural collapse occurred, the more severe the collapse process was, and the shorter the collapse process lasted. We thoroughly analyzed the collapse process of steel frame by applying a one-way earthquake wave, tracked the internal force change and collapse modes of the structure, and provided reference to the structure resistance design against the continuous collapse.The finite particle method is a newly developing method. The collapse of structures under earthquake excitation involves a series of complex problems. In the future, such investigations will be focused on applying earthquake acceleration or displacement time history to the structural foundation to consider the effect of spatial variability. The finite particle method should continue to be developed in the study of structural resistance design against continuous collapse.



---

*Source: 1952050-2018-06-13.xml*