Stability detection of building concrete structure based on discrete element method

Abstract

In order to improve the safety of building structure, the discrete element method is introduced to test the stability of building concrete structure. Taking reinforced concrete as the research object and using Mohr-Coulomb criterion as the concrete failure criterion, determine the stiffness coefficient, damping coefficient and time step of reinforced concrete. Simulation of failure modes of reinforced concrete beams under different shear-span ratios. Determining the collapse process of reinforced concrete structure under the action of strong earthquake El-Centro north-south wave, and testing the stability of reinforced concrete structure. The results show that at the moment of 3.25 s, the bottom spring of the bottom column of the frame structure is damaged, and the upper end of the right column of the middle two layers is deformed and damaged. At this moment, the frame structure collapses. As the shear-span ratio increases, the failure mode of reinforced concrete beams evolves in the order of diagonal pressure, shear pressure and cable tension. Its shear strength is gradually weakened, and the proposed method can effectively simulate the stability of concrete structures under the action of earthquakes and collapse.

Keywords

Discrete element method buildings concrete structure failure criterion collapse stability detection

1 Introduction

As the main form of building structure, concrete structure also plays an important role as the most essential building material in the current construction industry [4]. As we all know, reinforced concrete, as a form of concrete, is also widely used in buildings. Generally, reinforced concrete members work with cracks. In the process of concrete pouring and hardening, small cracks appear on the surface of concrete due to temperature and load. At the same time, during the use of reinforced concrete buildings, factors such as load effect, Cl-1 ion erosion, carbonation, steel corrosion and so on [25], will cause cracks on the surface of concrete structures, and reduce the durability of structures. In the use of concrete materials, the stability of materials will have a great impact on the reliability of buildings. Therefore, in the concrete construction process, it is necessary to do a good job in the concrete stability test. In the specific test, the damage to the original structure should be minimized to ensure the accuracy of the test results [14]. The finite element method is generally used to simulate the mechanical properties and failure of reinforced concrete structures. But the research object of the finite element method is mainly continuum, and the failure process of brittle materials such as concrete and rock under impact, earthquake and other dynamic loads is a process from continuous to discontinuous to failure [9], so it is not suitable to simulate it with the finite element method.

Discrete element method (DEM) is a numerical method developed in the 1970 s to analyze jointed rock. Its basic idea was originally proposed by Cundall. It can be used to simulate the fracture and failure process of non continuous rock. At first, the discrete element method is only suitable for the rock mass, but it can also be used to simulate the failure process of concrete. It has been proved that the discrete element method is very feasible to simulate the continuous brittle materials such as concrete. The model proposed by Saw amoto et al. Simulates the deformation process of brittle concrete materials and obtains effective experimental results, which has a certain effect. However, this model does not explicitly consider the role of reinforcement, it can be said that it is just a “plain concrete” model. In her research, Wang Zhuolin indicated that the establishment of numerical models should consider three key issues: failure criteria, constitutive relations, and methods for determining material parameters. The deformation law of concrete structure is analyzed in detail, but the influence of concrete structure collapse reaction is not considered. In this paper, the effects of concrete and reinforcement are considered at the same time [7]. The mechanical model, calculation model, element stiffness and other coefficients of the discrete element method are emphatically analyzed. The collapse response of the structure of the discrete element method, the failure mode of the reinforced concrete beam under different shear span ratio, and the collapse process of the occurrence of the reinforced concrete structure under the action of the strong earthquake El-Centro north-south wave by using the SCFS program are simulated and analyzed. By analyzing the state of concrete structure damaged and collapsed under earthquake, the stability detection of building concrete structure based on discrete element method is realized [16].

2 Basic definitions

2.1 Basic idea of discrete element method

In the discrete element method, the research body is divided into a set of several elements, which can be separated in the process of motion, that is, one element can contact or be separated from its adjacent elements. The force of interaction between elements can be calculated according to the relationship between force and displacement [2], while the motion of a single element is completely determined according to the magnitude of the unbalanced force and moment of the element according to Newton’s law of motion. The discrete element method, like the finite element method, divides the region into elements. However, the element is controlled by discontinuities such as joints. In the later equations of motion, the element nodes can be separated, that is, one element can be contacted with its adjacent elements or separated [15]. The force of interaction between elements can be calculated according to the relationship between force and displacement, while the motion of individual elements is completely determined according to the magnitude of the unbalanced force and moment of the element according to Newton’s law of motion. Discrete element method is an explicit numerical method. This method is similar to other explicit calculations in the time domain, such as the explicit difference scheme for solving parabolic partial differential equations. “Explicit” refers to the properties of algebraic equations used in numerical calculation of a physical system. In the calculation of explicit method, the quantities on one side of the equation are known, while the quantities on the other side can be obtained by simple substitution method. This is different from implicit method, which must solve simultaneous equations [1 , 18]. When the explicit method is used, it is assumed that in each iteration time step, each block element only produces the force influence on its adjacent block element, so that the time step needs to be small enough to make the explicit method stable. Since the matrix does not need to be formed when using the explicit method, large displacement and nonlinearity can be considered without additional calculation time.

Now the reinforcement is simulated, but the stiffness of the element is different from that of the concrete element. Because the diameter of the element used in the calculation is generally larger than the diameter of the reinforcement, the actual reinforcement element is the combination of the reinforcement and its surface concrete [22, 23]. The slip between the steel bar and its surface concrete is ignored here. When the interaction between the steel bar unit and the concrete unit occurs, it is actually the interaction between the concrete [19]. Therefore, the rigidity of the concrete is taken as the rigidity in the calculation.

2.2 Interactions during the collapse

The interaction states between elements are divided into two types:

State 1: spring is normal, it can bear not only pressure, but also tension, so there are pressure, tension and shear force between units;

State 2: when cracks occur, the cement mortar loses the ability to resist the tension. Only when the cells are squeezed each other, the spring will work. Therefore, there is only pressure and shear force between the cells [5].

Mohr-Coulomb criterion is adopted as the failure criterion of concrete, as shown in Fig. 1, where, F_t and F_c are the tensile strength and compressive strength of the element respectively.

Fig.1

Failure criterion of concrete.

It can be seen from Fig. 1 that when the normal force between the elements exceeds the tensile strength of the material, the state of interaction between them becomes state II; or when the shear force exceeds (C + _ f_n) (C is the viscosity constant, f_n is the normal force between the elements), the state of interaction also becomes state II, and the shear force becomes _ f_n (_ = tanh, h is the friction angle).

2.3 Basic equation

Taking the plane problem as an example, the basic equations are introduced. In the coordinate system as shown in Fig. 1, it is specified that the displacement increment Δu_n in the normal direction is positive in the far away direction, the displacement increment in the tangent direction is positive in the far away direction, and the displacement increment Δu_s in the tangent direction is positive in the counterclockwise direction, with the values as follows: $\begin{matrix} Δ u_{n} \\ = (Δ u_{j} - Δ u_{i}) cos θ \\ + (Δ v_{j} - Δ v_{i}) sin θ, \\ Δ u_{s} \\ = - (Δ u_{j} - Δ u_{i}) sin θ \\ + (Δ v_{j} - Δ v_{i}) cos θ - (r_{i} Δ H_{i} - r_{i} Δ H_{j}) \end{matrix}}$ (1)

In Equation (1), Δu_i, Δv_i, Δu_j, Δv_j, ΔH_i and ΔH_j are the displacement increment of i and j in Δt along the X, Y axis and the rotation direction, respectively, and θ is the angle between the center line of two circles and the coordinate axis at time t. According to Hooke’s law, the normal force f_n (the pressure is positive) and the shear force f_s in tangent direction (the anticlockwise is positive) between the two elements of i and j at time t are as follows: $\begin{matrix} f_{n} : = f_{n} - k_{n} Δ u_{n}, \\ f_{s} : = f_{s} - k_{s} Δ u_{s} \end{matrix}}$ (2)

In Equation (2), k_n and k_s are element normal stiffness coefficient and element tangent stiffness coefficient respectively. “:=” is the assignment symbol, indicating that the variable on the right end is assigned to the variable on the left. In order to consider the energy loss of the element during operation [3 , 21], the damping force should be added to the calculation earlier. It is assumed that the damping force at the contact point is proportional to the increment of displacement: $\begin{matrix} D_{n} = - {Uk}_{n} Δ u_{n}, \\ D_{s} = - {Uk}_{s} Δ u_{s} \end{matrix}}$ (3)

In Equation (3), U is the stiffness damping coefficient. Then, the above forces are transformed into components under the overall coordinates: $\begin{matrix} F_{x} = - (f_{s} + D_{s}) sin θ - (f_{n} + D_{n}) cos θ, \\ F_{y} = (f_{s} + D_{s}) cos θ - (f_{n} + D_{n}) sin θ \end{matrix}}$ (4)

The resultant force of element i can be calculated by summing up all the contacts of element i and adding the external load of element i. $\begin{matrix} F_{xsum} = \sum_{j} F_{x} + F_{xload}, \\ F_{ysum} = \sum_{j} F_{y} + F_{yload}, \\ M_{sum} = r_{i} \sum_{j} (F_{s} + D_{s}) \end{matrix}}$ (5)

In Equation (5), F_xsum, F_ysum and M_sum are the resultant forces and moments of element i in X and Y directions respectively. According to Newton’s second law of motion, the equation of motion of element i can be obtained as follows: $\begin{matrix} {mu}^{″} + {Tmu}^{'} = F_{xsum} \\ {mv}^{″} + {Tmv}^{'} = F_{ysum} \\ {IH}^{″} + {TmH}^{'} = M_{sum} \end{matrix}}$ (6)

In Equation (6), m is the mass of the element, T is the mass damping coefficient, and I is the moment of inertia of the element. After the first-order central difference of Equation (6), it can get:

$\begin{matrix} u^{'} (t + Δ t / 2) = [u^{'} (t - Δ t / 2) (1 - T Δ t / 2) + \frac{F_{x}}{m} Δ t] / (1 + T Δ t / 2), \\ v^{'} (t + Δ t / 2) = [v^{'} (t - Δ t / 2) (1 - T Δ t / 2) + \frac{F_{y}}{m} Δ t] / (1 + T Δ t / 2), \\ H^{'} (t + Δ t / 2) = [H^{'} (t - Δ t / 2) (1 - T Δ t / 2) + \frac{M}{I} Δ t] / (1 + T Δ t / 2) \end{matrix}}$ (7)

Therefore, the displacement and rotation increment of the element from t to t + Δt and the position of the center of the circle at t + Δt are respectively as: $\begin{matrix} Δ u = u^{'} (t + Δ t / 2) Δ t \\ Δ v = v^{'} (t + Δ t / 2) Δ t \\ Δ H = H^{'} (t + Δ t / 2) Δ t \end{matrix}}$ (8) $\begin{matrix} x : = x + Δ u, \\ y : = y + Δ v \end{matrix}}$ (9)

Thus, the velocity, acceleration, displacement and resultant force of all elements at any time can be obtained.

2.4 Selection of unit parameters

2.4.1 Determination of stiffness coefficient

In the analysis of plane problems, circular elements are generally used, and the combination of elements is hexagon [24]. For the plane stress state, the stiffness coefficient is defined as follows: $\begin{matrix} k_{n} = \bar{3} Et / (3 (1 - v)), \\ k_{s} = \bar{3} (1 - 3 v) Et / (3 (1 - v^{2})) \end{matrix}}$ (10)

For plane strain and axisymmetric stress state, the stiffness coefficient is defined as follows: $\begin{matrix} k_{n} = \bar{3} Ek / (3 (1 + v) (1 - 2 v)), \\ k_{s} = \bar{3} (1 - 4 v) Ek / (3 (1 + v) (1 - 2 v)) \end{matrix}}$ (11)

In Equation (11), E is the elastic modulus of reinforcement, t is the thickness of concrete, and the plane strain state k is 1. For axisymmetric stress state, k is 2πb (b is the distance from the symmetry axis to the center of the element). The rigidity of the reinforcement unit is: $K_{B} = \frac{EA}{L}$ (12)

In Equation (12), A is the cross-sectional area of the reinforcement and L is the distance between the centers of two elements.

2.4.2 Determination of damping coefficient

In the discrete element method, the influence of damping must be considered [13]. In the combined spring model, the damping coefficients c_in, c_im and c_is of normal direction, shear and bending are considered. However, in the process of solving, the damping force will produce a coupling term, which increases the complexity of solving. In order to simplify the calculation, the damping coefficient of the whole coordinate is adopted in this paper. $\begin{matrix} c_{ix} = α m_{i} + β_{x} \sqrt{m_{i} k_{ix}} \\ c_{iy} = α m_{i} + β_{y} \sqrt{m_{i} k_{iy}} \\ c_{i θ} = α I_{i} + β_{θ} \sqrt{J_{i} k_{i θ}} \end{matrix}}$ (13)

In Equation (13), m_i, I_i, k_ix, k_iy and k_iθ respectively represent unit mass, section moment of inertia, linear stiffness and rotational stiffness of connecting spring of unit i (which can be calculated by resultant force on unit i, deformation between unit i and adjacent units); α, β_x, β_y and β_θ are all undetermined coefficients, and α = 0.05, β_x = β_y = 1.5 ∼ 2.0, β_θ = 2.5.

2.4.3 Determination of time steps

The dynamic relaxation method is used to calculate the structure, i.e. it is assumed that the elements are calculated one by one in one time step. In fact, this method assumes that the unbalanced force generated by each element when it is relaxed in one time step is only transmitted to its adjacent elements [6], which has no effect on other elements. Therefore, in order to make the calculation result close to the reality, the time step must be small enough. Therefore, the equation of calculation time step determined in this paper can be shown as follows: $Δ t < \frac{l_{min}}{v_{p}}$ (14)

In Equation (14), it is the minimum unit length; v_p is the propagation speed of P wave (longitudinal wave); according to the wave theory, it can be calculated that for reinforced concrete structure v_p = 3400m/s ∼ 3500m/s, the above time step is a theoretical upper limit. In order to ensure the convergence of the results [3], Δt ≤ 10^-5s is adopted in this paper.

2.5 Analysis of structural collapse response by discrete element method – synchronous dynamic relaxation method

The dynamic relaxation method has path dependence in solving the response of the whole structure, that is to say, the solution result is related to the relaxation order of the element, and different relaxation order may have different calculation results.

In order to reduce the errors caused by the distance correlation and to take into account the simplified numerical calculation process [11], the synchronous dynamic relaxation method is adopted in this paper. The basic idea is to fix all elements of the structure in a time step calculation process from the known initial state, and determine the resultant force and resultant moment of each element according to the connection spring force- displacement relationship between elements. Then all the elements are relaxed at the same time, and the motion displacement of a single element is calculated by the central difference method. All the cells are fixed in the new position and enter the next iteration [8]. According to the principle of the synchronous dynamic relaxation method and the characteristics of the discrete element model of the bar section, this paper compiles the calculation flow chart suitable for solving the collapse process of reinforced concrete structure under the earthquake action, as shown in Fig. 2.

Fig.2

Calculation flow chart of synchronous dynamic relaxation method.

3 Results

3.1 Damage analysis of reinforced concrete beam

Based on the VC6.0 programming platform, a discrete element program for plane problem analysis is developed. The program can be used to calculate and simulate the deformation and failure process of different materials under different loads and boundary conditions. The load and boundary conditions of the whole reinforced concrete structure are simulated by adding initial load and displacement constraints to some elements. The above program is used to simulate the failure of oblique section of beams without web reinforcement under different shear span ratios. The height, width and length of the simulated beam are 0.25, 0.12 and 1.00 m respectively; the beam is equipped with two Φ 20 longitudinal bars, and its physical properties of concrete material are the same as the above example; the constraint condition of the beam is: the lower end of the beam is supported by two supports; the loading condition is: adding concentrated load at beam A (corresponding shear span ratio λ= 0.8), B (λ= 1.3), C (λ= 2.0), as shown in Fig. 3, the radius of the discrete element r = 0.015 m, time step Δt = 10-6 s.

Fig.3

Reinforced concrete beam model.

When loading, the load of each level is 10 kN, which increases gradually until the beam is damaged. The maximum failure load in condition (a) is 80 KN; the failure load in condition (b) is 50 kN; and the maximum failure load in condition (c) is 40 kN. Figure 4 shows the final failure mode of the beam.

Fig.4

Failure mode of reinforced concrete beam.

It can be concluded that with the increase of shear span ratio λ, the failure mode of reinforced concrete beam evolves in the order of diagonal pressure, shear pressure and cable-stayed, and its shear strength gradually decreases, which is quite consistent with the test results. The reason is that the shear span ratio λ is essentially related to the ratio of normal stress e to shear stress f on the section. e and f determine the magnitude and direction of the main stress. Therefore, it is bound to have a great impact on the shear performance of the inclined section. This procedure is only applicable to the plane situation, so the reinforcement at the lower end of the beam is similar to the steel plate, which blocks the extension of the crack and makes the crack development at the bottom of the beam not obvious.

3.2 Example of four-story reinforced concrete structure

A reinforced concrete structure with four stories and one span is chosen as an example. In order to make the structure collapse under strong earthquake, the structure is designed as strong beam and weak column. In the design, the site type is class II, the first group with fortification intensity of 8 degrees, the seismic grade of the frame is class 3, the importance coefficient is 1.0, the design section of the frame column is 300 mm×300 mm, the beam section is 250 mm×500 mm, and the plate thickness is 120 mm. C30 concrete is selected, the main reinforcement is grade II, the stirrup is grade I, and the elevation of the structure is shown in Fig. 5.

Fig.5

Elevation of reinforced concrete frame structure with four floors and one sp.

In order to simulate the collapse of the structure set up in this paper, the seismic wave needs to be input. Different seismic waves will cause different seismic responses of structures, so the selection of seismic waves is very important. At present, El-Centro seismic record is widely used at home and abroad. It has a larger acceleration peak, and its waveform can produce a larger seismic response when the acceleration peak is the same. Therefore, north-south wave of EI-Centro is selected in the simulation analysis of this paper.

El-Centroel-Centro wave record is a record of IMPERIALvALLEY, 9 km from the epicenter of the Richter scale 6.4 earthquake in 1940 in California, USA. The peak value of the north-south horizontal acceleration is 341.7 Gal. This record is the basis for the design of response spectrum in various countries, and is widely used in various theoretical analysis and experimental research. The acceleration time history of El-Centro’s north-south wave is shown in Fig. 6. The ordinate is acceleration, and the unit is centimeter per quadratic second (Gal). The abscissa is time, and the unit is second (s). Considering that the purpose of the study is to ensure the collapse of the structure under the action of earthquake, the peak acceleration of El-Centro’s north-south wave is increased to 0.59, and then the wave is multiplied by the corresponding coefficient.

Fig.6

North South wave for EI-Centro.

According to the plane characteristics of the structure, the beams and columns are discretized. The rules of discretization adopt the failure criteria in this paper. Then, the input file example.txt of SCFS program is compiled according to the relevant input format. The initial discrete element model data is input through the open button of the program, as shown in Fig. 7.

Fig.7

Model data input of SCFS program.

After inputting the initial model data, the acceleration file of El-Centro’s north-south wave is input through the calculation menu, and then is calculated. The time required for calculation depends on the number of model units of the structure and the response time of seismic calculation. When the program calculation is completed, the corresponding animation demonstration file will be generated in the current folder. Finally, the seismic collapse response of the four-story reinforced concrete structure under the El-Centro’s north-south wave with a peak value of 0.59 can be viewed through the animation demonstration under the animation menu. The detail is shown in Fig. 8.

Fig.8

Collapse response of reinforced concrete beam structure under El Centro North South wave.

In the diagram of reinforced concrete structure model (Fig. 8a), there are gaps between the elements, representing the virtual spring. Because the rigid body elements cannot penetrate and the rigid body elements cannot penetrate, the existence of gaps allows the relative deformation between the rigid body elements, which is in line with the actual situation. After the local seismic wave is input, the frame structure will vibrate left and right under the action of earthquake. Fig. 9. is the horizontal displacement time history of each floor of frame structure before 3.25 s of earthquake action. At 3.25 s, the bottom spring of the bottom column of the frame structure is damaged, and the upper end of the right column of the middle two floors is also damaged due to large deformation. At this moment, the frame structure collapses, as shown in Fig. 8b. Later, due to the destruction of the spring between the column units, there is a large relative movement between the column units, resulting in the collision phenomenon, which leads to the continued failure of the structural components. The whole reinforced concrete structure collapses and all structural members scatters on the ground, as shown in Fig. 8c–8e.

Fig.9

Horizontal displacement time history of top and bottom floor of frame structure before 3.25 s of earthquake action.

From the acceleration time curve of El-Centro’s north-south wave, we can get that the seismic acceleration reaches the peak acceleration in about 3.25 s, so the corresponding frame structure will have a larger seismic response at this time. If the peak acceleration is large, the structure is prone to collapse. Therefore, in the SCFS program demonstration, the collapse of the frame structure in 3.25 s is more in line with the actual situation.

4 Conclusions

In this paper, the discrete element method is used to simulate the reinforced concrete structure, and the SCFS program is used to simulate the collapse process of reinforced concrete structure under the action of El-Centro’s north-south wave of strong earthquake.

(1) When using the discrete element method to simulate the reinforced concrete, it is found that the shear span ratio of the reinforced concrete beam has a great influence on the shear performance of the inclined section and affects the stability of the reinforced concrete on the basis of the reasonable selection of the calculation parameters; when using the SCFS program to simulate the collapse process of the reinforced concrete structure under the action of El-Centro’s north-south wave of strong earthquake, it is found that at 3.25 s, the seismic acceleration reaches the maximum, and the structure exits the work due to the damage of some structural members, resulting in the collapse of the structure at this moment. At the next moment, due to the collision between the members, secondary damage is resulted, and finally the structure is scattered on the ground.

(2) Therefore, the discrete element method can effectively simulate the process that the concrete structure is destroyed under the earthquake, leading to the instability and collapse of the structure. It can be seen that the method in this paper can be used to simulate the reinforced concrete structure by the discrete element method, and the SCFS program is used to simulate the collapse process of the reinforced concrete structure under the action of El-Centro’s north-south wave of the strong earthquake, so as to realize the stability detection of the building concrete structure.

In the future, the accuracy of the detection of the stability of concrete structures will be improved to contribute to the construction industry.

References

Wei

, Ji

, Ming

, Simulating the workability of fresh self-compacting concrete with random polyhedron aggregate based on DEM, Materials & Structures50(1) (2017), 92–95.

Wei

, Yan

W.S.

, Song

H.F.

, Blocking analysis of fresh self-compacting concrete based on the DEM, Construction & Building Materials168 (2018), 412–421.

Villanueva

, Correa

, JosÉ

, Pacas

, Control of a single-phase cascaded H-Bridge multilevel inverter for grid-connected photovoltaic systems, IEEE Transactions on Industrial Electronics56(11) (2009), 4399–4406.

Smoljanović

, Živaljić

, Željana

and Munjiza

, Numerical model for confined masonry structures based on finite discrete element method, International Journal of Civil Environmental Structural Construction & Architectural Engineering30 (2017), 19–35.

, Qi

, Progressive collapse simulation based on DEM for single-layer reticulated domes, Journal of Constructional Steel Research128 (2017), 721–731.

J.J.

, Research on the function of electronic medical record and related problems in hospital informatization management, Automation & Instrumentation6 (2017), 226–227, 232.

Shi

, Zhao

, Sun

, Parameter calibration of soil particles contact model of farmland soil in northwest arid region based on discrete element method, Transactions of the Chinese Society of Agricultural Engineering33(21) (2017), 181–187.

Niu

L.L.

, Su

, Niu

L.L.

, Shu

, Research on construction technology and risk management of steel structure based on BIM, Automation & Instrumentation Automation & Instrumentation8 (2017), 165–167.

Zhou

L.Q.

, Yang

, Simulation of a hot-rolled h-section steel beam subject to static loading based on discrete element method, Procedia Engineering210 (2017), 312–319.

10.

Horibe

, Sonoda

, Watano

, Scale-up of lubricant mixing process by using v-type blender based on discrete element method, Chemical & Pharmaceutical Bulletin66(5) (2018), 548–553.

11.

Tang

, Yang

, Li

X.F.

, Two-grid finite element discretization methods for a class of poisson-nernst-planck equations, Journal of Jilin University (Science Edition)57(3) (2019), 71–77.

12.

Celigueta

M.A.

, Latorre

, Arrufat

, Accurate modelling of the elastic behavior of a continuum with the discrete element method, Computational Mechanics60(1) (2017), 1–14.

13.

N.J.

, Zhou

, Zheng

J.L.

, Chen

, Preparation and Electrochemical Performance of Porous V_2O_5 Micro-nanospheres, Chinese Journal of Power Sources42(10) (2018), 1444–1446.

14.

Yang

, Gao

, Li

, Simulation of fracture behaviour of hydrogel by discrete element method, Iet Micro & Nano Letters13(6) (2018), 743–746.

15.

Kaewunruen

, Mirza

, hybrid discrete element - finite element simulation for railway bridge-track interaction, IOP Conference Series: Materials Science and Engineering251(1) (2017), 12–16.

16.

Sinaie

, Application of the discrete element method for the simulation of size effects in concrete samples, International Journal of Solids & Structures108(1) (2017), 244–253.

17.

Smirnov

, Volchenkov

, Five years of phase space dynamics of the standard poor’s 500, Applied Mathematics and Nonlinear Sciences4(1) (2019), 209–222.

18.

Kkl

, Three-dimensional discrete element method for the prediction of protoplasmic seepage through membrane in a biological cell, Journal of Biomechanics65(8) (2017), 115–129.

19.

X.F.

, Zhang

Q.B.

, Li

H.B.

, Grain-based discrete element method (gb-dem) modelling of multi-scale fracturing in rocks under dynamic loading, Rock Mechanics and Rock Engineering51(12) (2018), 3785–3817.

20.

X.L.

, Xu

Y.J.

, Zhou

Z.X.

, Sparse event detection based on parallel discrete social spider optimization algorithm and compressed sensing in wireless sensor networks, Journal of China Academy of Electronics and Information Technology12(2) (2017), 202–208.

21.

Sagna

, Multidimensional BSDE with Poisson jumps of Osgood type, Applied Mathematics and Nonlinear Sciences4(2) (2019), 387–394.

22.

Liu

, Economic analysis of methanol production from coal/biomass upgrading, Energy Sources Part B-Economics Planning and Policy13(1) (2018), 66–71.

23.

Liu

, Experimental analysis of an integrated biomass gasification and PEM fuel cell system, Energy Sources, Part A: Recovery, Utilization, and Environmental Effects41(3) (2019), 360–367.

24.

Liu

, What is the future of solar energy? Economic and policy barriers, Energy Sources Part B-Economics Planning and Policy13(3) (2018), 169–172.

25.

Wang

Z.H.

, Yan

X.C.

, The analysis of long-and-short concrete column structure under horizontally seismic process based on finite element method, Key Engineering Materials737 (2017), 448–453.