Nonlinear analysis of engineering structures by combined finite-discrete element method

A numerical model for dynamic analysis of dry-stone masonry and reinforced-concrete structures, based on the combined finite-discrete element method (FEM/DEM), is presented in the paper. This model describes behaviour of such structures when exposed to dynamic load, crack initiation and propagation, energy dissipation mechanisms due to non linear effects, inertial effects due to motion, contact impact and achievement of the state of rest as a consequence of energy dissipation mechanisms in the system.


Introduction
In modern civil engineering, many engineering structures are realized either as reinforced-concrete or masonry structures.Despite differences in their behaviour, they also present some common features, particularly with regard to crack initiation.The nonlinear analysis of reinforced-concrete and masonry structures exposed to static and dynamic load has imposed itself over the past several decades as a technique that enables us to monitor behaviour of such structures all the way until failure, with an accuracy that varies depending on the model or method used for this purpose.Crack initiation is recognised as one of significant causes of nonlinear behaviour of reinforced-concrete and masonry structures.Most models described in literature, by which the behaviour of reinforced-concrete and masonry structures is simulated, are based on the finite element method.In reinforced-concrete structures, the cracking of material is described using the smeared crack model [1][2][3][4][5][6][7] or the discrete crack model [8].In the smeared crack model [1][2][3][4][5][6][7], the cracked concrete is modelled as an elastic orthotropic material with a reduced elastic modulus in the direction perpendicular to the crack.In this approach, local discontinuities in the displacement zone are distributed over the finite element, while an averaged relation between the stress and strain values is described with constitutive laws.In the discrete crack model, cracks are modelled as a geometrical discontinuity by separating finite elements [8] or by introducing basic functions through which the discontinuity in the displacement zone can be described [9,10].In the masonry structure modelling context, the distinction can be made between the micro modelling and macro modelling.Macro modelling is the masonry structure modelling approach that is now dominantly used in practice [11][12][13][14].In this approach, the structure is simulated as an orthotropic continuum in which the mean relation between the stress and strain is obtained experimentally or through homogenization techniques.A drawback of the concrete and masonry structure modelling by means of continuum is that the occurrence of large discontinuities in the displacement range can not be described.One of the ways to overcome this drawback is to introduce contact elements into the finite element mesh.In this approach, the behaviour of material in the finite element is most often linear-elastic, while the material nonlinearity is modelled in contact elements.The constitutive law of contact element behaviour is based on the theory of plasticity [15][16][17][18] or damage mechanics [19,20].A drawback of models based on the finite element method is the impossibility to simulate mechanical interaction between several bodies, which is important during analysis of structures exposed to impact load, and also during analysis of progressive collapse of structures.Models based on the discrete element method have been developed for the analysis of problems involving mechanical interaction between several bodies which might have large rotations and displacements.This method was initially used for simulating sliding and detachment of rock mass along joints defined in advance [21], and has later on found a suitable application in the analysis of masonry structures [22][23][24] and reinforced-concrete structures [25].The main feature of the discrete element method that has enabled its application in the analysis of structures, masonry structures in particular is the presentation of a structure as a set of discrete elements connected with each other by contact elements.This approach enables us to simulate structural collapse due to rotation, sliding, and impact load.A shortcoming of models based on the discrete element method is their impossibility to describe the state of stress and strain within discrete elements, which is highly significant in the analysis of the initiation and development of cracks.Many attempts have been made in recent times to make use of advantages presented by the finite and discrete element method [26][27][28].In the scope of this method, the behaviour of material until initiation of cracks is modelled as in the finite element method, while the discrete crack appears at the moment when the tensile strength of material is exceeded.The appearance of cracks and fragmentation of discrete elements is comprised in contact elements that are modelled between finite elements.An advantage of the combined finite-discrete element method lies in the possibility to describe phenomena such as the behaviour of the structure due to dynamic action in linear-elastic phase, initiation and propagation of cracks, inertia effects due to motion, interaction resulting from dynamic contact and, finally, achievement of the state of rest which occurs as a consequence of energy dissipation in the structure.For that reason, this method has imposed itself as a very appropriate tool for the development of a new numerical model for concrete and reinforcement that enables simulation of response of reinforced-concrete structures exposed to dynamic load [30,31].Studies have also been made to demonstrate advantages of using the combined finite-discrete element method in the non-linear analysis of dry stone masonry structures [32].

Combined finite-discrete element method
The combined finite-discrete element method (FEM/DEM) [33][34][35] is based on the simulation of behaviour of a great number of discrete elements that can be in an interaction.Each discrete element is discretized with its own mesh of finite elements and hence its deformability is enabled.The material nonlinearity, including the occurrence and development of cracks and, finally, fragmentation of discrete elements, is enabled through the model of contact elements that are implemented between finite elements.In order to take into account all these effects, appropriate algorithms have been developed in the scope of this method.These algorithms cover, in every time interval, the detection and interaction of a contact, monitoring of stress and strain in a finite element, Nonlinear analysis of engineering structures by combined finite-discrete element method initiation and propagation of cracks, integration of equations of motion in time which includes large displacements and rotations and visualisation of the mentioned effects.The stress and strain relationship in the finite triangular element is taken into account through the Hooke's law, according to the following expression where T is the Cauchy distribution tensor, E is the elastic modulus, u is the Poisson's ratio, d E  is a part of the Green-St.Venant strain tensor relating to the change of shape, s E  is a part of the strain tensor relating to the change in volume, m is the damping ratio, and D is the strain rate tensor [34].

Contact detection and interaction
The objective of the contact detection algorithm is to find pairs of neighbouring discrete elements that are in contact, and to eliminate pairs that are too far away and can no longer remain in contact.The NBS (no binary search) contact detection algorithm has been implemented in the FEM/DEM model [33,34].The total time needed to detect all contact pairs is proportional to the total number of discrete elements.Once the discrete element pairs have been detected, the contact forces between two discrete elements in contact are defined.One of these elements is the contactor, and the other is the target (Figure 1).In the interaction algorithm, the distributed contact forces are defined by means of the penalty method which is based on the principle of potential contact forces.When in contact, the contactor and target overlap across the surface S which is bounded by the external edge G βm∩βk .Then the total differential contact force at the contactor df k is defined as where P m and P k are the target and contactor overlapping points, while j is the corresponding function of the potential.The total contact force is obtained by integration (2) over the entire overlapping surface S which can also be written as where n G is the unit external normal to the edge G of the overlapping surface S.
The Coulomb dry friction model for shear forces is implemented in the scope of the contact force algorithm as follows where f t is the tangential elastic contact force, k t is the penalty coefficient for friction, and δ t is the tangential vector of displacement between two elements [36].If f t is greater than the maximum friction force defined by the Coulomb law, |f t |>µ|f n | then the elements slide one along the other while the shear force between them is defined through the elastic normal force f n based on where µ is the friction coefficient.

Crack model
The crack model implemented into the combined finitediscrete element method is used for simulating the initiation and propagation of cracks in brittle materials subjected to load modes I and II.The crack initiation in tension is described in mode I, while the crack in shear is described in mode II.The model is based on approximation of the experimental stressstrain curve for concrete in direct tension [37].The area under the curve of stress-strain in tension is divided into two parts, as shown in Figure 2. In this model, the part "A" is implemented into the behaviour of finite elements in a standard way through the constitutive law of material behaviour [31,34].The part "B" presents tensile softening where the stress reduces with an increase in strain [31,37].This is modelled with the discrete crack model, shown in Figure 3, in which it is assumed for simplicity that the crack is aligned with the edge of the finite element.Separation of edges of two Hrvoje Smoljanović, Nikolina Živaljić, Željana Nikolić neighbouring finite elements induces stress which is taken as the separation function δ, as shown in Figure 2b.
The surface under the stress-strain curve from the moment when the stress falls to zero represents the fracture energy G f .This is the work that has to be applied to create a crack of a unit area.Theoretically, the separation of edges of two neighbouring finite elements should be equal to zero all the way until the tensile strength of material is attained, which would mean that δ t =0.In the presented model, the separation of neighbouring edges of two finite elements is ensured by topology of finite elements in such a way that not a single node belongs to two finite elements.The continuity between finite elements until the time the tensile strength is ensured by means of the penalty method [33].Springs of great stiffness are modelled at the edge of a finite element in the direction of the normal (Figure 4a), and so the relation δ t =δ p is valid.The following relation is valid for the separation of δ < δ p (Figure 2b) where: is the separation in the moment when the stress corresponds to the tensile strength of material f t , h is the size of the finite element, while p 0 is the penalty coefficient by which the relative error of numerical solution is controlled.In the limit case when the separation of edges of two neighbouring finite elements is equal to zero, which corresponds to the moment in which the tensile strength of material f t has been achieved.The analysis of influence of penalty coefficient on the relative error in normal displacements shows that the error is less than 1 % for p 0 =100 E, where E is the modulus of elasticity of material.The stress between finite element edges falls with an increase in separation δ>δ p (Figure 2b), and the stress becomes σ c =0 at the moment of δ=δ c .The following relationship between stress and displacement has been adopted for the area ofδ c >δ>δ p .
where z is the function of the experimental curve that describes the behaviour of concrete in tension [37] with coefficients c 1 =3.00 and c 2 =6.93.Cracks subjected to load in mode II are assumed to behave in a way similar to that shown for mode I. Until the moment the shear strength of material is achieved, the edges of two neighbouring finite elements are supported with shear stresses that are calculated by means of the penalty method [33].Edges are supported with shear springs (Figure 4b) as shown in the expression below Nonlinear analysis of engineering structures by combined finite-discrete element method where: is the separation in the moment when the stress corresponds to the shear strength of material f s , h is the size of the finite element, and p 0 is the penalty coefficient.
In the limit case when the sliding of edges of two neighbouring finite elements is equal to zero, which corresponds to the moment in which the shear strength of material f s has been achieved.The stress between finite element edges falls with an increase in sliding t>t p , and the stress becomes t c =0 at the moment of t=t c .The following relationship between the stress and sliding has been assumed for the area of t c >t>t p where D is defined by the relation The complete relation describing the relationship τ c -t in mode II can be presented as follows In case the crack is subjected to load in mode I and mode II, then the same expressions as previously explained are used to calculate normal or shear stresses, but the damage factor D is adopted in this case.This factor is defined as ( )/ ( for and for or for and The fracture criterion is defined by

Numerical model of reinforcement in the FEM/ DEM model
In the scope of the FEM/DEM method, structures are discretized through the mesh of finite and contact elements, as shown in Figure 5. Reinforced-concrete structures are modelled by modelling the concrete using triangular finite elements, while the reinforcement is modelled by two-node bar elements [38,39] embedded in finite elements of concrete.
In the first step, each reinforcing bar is introduced as one bar.
Once the structure is divided into a certain number of finite elements, the search is made to identify reference points or intersects of reinforcing bars and ends of finite elements of concrete.This is how linear finite elements of reinforcement are formed (Figure 5).
Prior to appearance of cracks in concrete, the structure is in a linear-elastic range, while triangular finite elements of concrete and two-node elements of reinforcement act as a single body.The deformation of the triangular element influences deformation of the finite reinforcement element.This causes stresses in reinforcement which in turn results in forces whose influence is taken into account in form of equivalent forces in triangular-element nodes.This is how the bond between the concrete and steel in the linear-elastic range is established.The appearance of cracks, i.e. the separation of edges of neighbouring finite elements, is enabled by describing neighbouring edges of triangular finite elements of concrete, and end points of neighbouring linear finite elements of reinforcement, by means of different nodes.The initiation and propagation of cracks in concrete takes place in the 2D contact element of concrete placed in between edges of triangular finite elements.At the same time, the reinforcement within the contact element of concrete is deformed, and its nonlinear behaviour is modelled with the linear contact element of reinforcement inserted in between neighbouring nodes of finite elements of reinforcement.In this model, it is assumed in the finite linear element of reinforcement that the relation between stress and strain is linearly elastic, and that there is no sliding, which means that the deformation of reinforcement and concrete in the finite element is the same [30].Nonlinear behaviour of reinforcement is modelled within the contact element of reinforcement.The model of reinforcement in the contact element [30] is divided into the part before occurrence of crack and the part after occurrence of crack.Before occurrence of crack, the contact element of reinforcement maintains the continuity between linear elements of reinforcement, which means that theoretically the separation of d neighbouring points of the contact element (Figure 6) should be equal to 0. In actual implementation, this has been achieved by means of penalty method [30].
The model of reinforcement behaviour in contact element is based on experimental curves that describe the state of deformation of reinforcing bar in crack, taking also into account plastic deformations due to cyclic load [40].In the contact element of reinforcement, the influence of curvature along the reinforcing bar in bending zone was modelled, as well as the influence of crack spacing [30].The improved Kato's material model was adopted for reinforcing steel as it enables simulation of cyclic behaviour of steel [41].

Numerical examples
The validation of the developed numerical model based on the combined finite-discrete element method is presented in this section, using reinforced-concrete beam and stone wall exposed to a monotonically increasing load.Literature examples with known results of physical experiments are used.In addition, the analysis also covers examples showing possibilities for using the developed model in the analysis of reinforcedconcrete and dry stone masonry structures subjected to impact load, and in the incremental dynamic analysis which was in this case conducted by incremental increase of the real earthquake amplitude until failure of the structure.

Bresler-Scordelis reinforced-concrete beam
The Bresler-Scordelis reinforced-concrete beam taken from literature [42] was selected in order to validate the implemented model for the case of monotonically increasing load, and for the adopted nonlinear properties of material.Geometrical properties of the analysed reinforced-concrete beam are shown in Figure 8. Properties of materials used in numerical analysis, taken from literature, are presented in Table 1.

Table 1. Properties of materials as input parameters
The beam is exposed to the action of the concentrated monotonically increasing force in mid span.The load that causes structural failure, i.e.F=258 kN, was obtained by experiment.Nonlinear analysis of engineering structures by combined finite-discrete element method Concrete is discretized with 2355 triangular finite elements, and each reinforcing bar with 168 two-node elements.The load is applied all the way until failure.The discretization of structure is shown in Figure 9. Results obtained by experiment, and by numerical analyses using the 3D FEM method, are presented in literature [42].The nonlinear numerical analysis with threedimensional discretization of the structure is conducted in the cited literature, and the displacement of point in mid span until failure is observed.Results of these analyses, and comparisons with results obtained by numerical model based of the FEM/ DEM method, are presented in diagram shown in Figure 10.The above diagram shows displacement of the point in the midspan of the beam as related to a given concentrated force F. The biggest deviation of results obtained by the FEM/DEM method is 3.8 %.The 0.8 % accuracy was obtained for the limit load which amounts to 258.1 kN according to experiment, and 260.2 kN according to the FEM/DEM method.Compared to 3D analysis with finite elements, the value of displacement immediately prior to failure obtained by FEM/DEM method shows a better correspondence with experimental results.The diagram presented in Figure 11 shows dependency between the stress in reinforcement and concentrated force F. The element below the point of the F force action in the bottom bar of the reinforcement was selected.It can be noticed that the liquid limit in reinforcement has not been achieved.It is significant to note that the full structural failure, both in the experiment and in numerical analysis conducted according to the developed model, occurs due to yield of concrete in compression.This demonstrates that this model describes very well the nonlinear behaviour of reinforcedconcrete structures for the case of failure along the concrete.

Dry stone wall exposed to monotonically increasing shear load
The validity of the developed numerical model is in this case analyzed during description of the stone wall behaviour at shear in its own plane.The experiment conducted by Oliveira [43] to study shear behaviour of stone walls was selected for this purpose.Numerical results obtained using the model based on the FEM/DEM method is compared with the results of the above mentioned experiment.The experimental program consisted of a series of quasi-static monotonic tests conducted on a small stone wall sample whose geometry is shown in Figure 13a.The wall discretization was used in numerical analysis presented in Figure 13b.The wall consisted of stone blocks of regular dimensions.Average values of mechanical characteristics of granite used in the experiment are taken from literature and are presented in Table 2.
Table 2. Average values of mechanical properties of stone used in the experiment [43] The shear behaviour of walls was analysed for two vertical longitudinal force values (30 kN and 100 kN), which corresponds to the pre-compression stress of 0.15 MPa and 0.5 MPa.The coefficient of friction between stone blocks was obtained by experiment, and it amounted to µ=0.62 .Average values of elastic modulus for wall, dependent on pre-compression load, and used in the numerical analysis, are presented in Table 3..Table 3.Average elastic modulus values for wall [43] Results obtained by experiments and calculations based on numerical model developed by Lourenco and Rots [17] are presented in the cited literature [43].This model is based on the finite element method containing also contact elements whose constitutive law of behaviour is based on the plastic theory.Results of these analyses, and comparison with results obtained by numerical model based on the FEM/DEM method, are presented in diagram shown in Figure 14.
A good correspondence can be observed between results obtained by model based on the FEM/DEM method, and numerical results obtained by Oliveira [43].The correspondence between numerical and experimental results can be considered satisfactory because the experiment was conducted with the dry wall made of natural stone.
In such walls, the influence of irregularities between blocks has a considerable effect on wall behaviour, and this effect is very hard to model by numerical procedure.The comparison of numerical results and results obtained by physical experiment is presented in Figure 15, for pre-compression stress values of 0.15 MPa and 0.5 MPa.The comparison of numerical results obtained by the FEM/DEM model with physical experiments shows that the failure mode obtained numerically is similar to that obtained by experiment.It can also be seen that the failure of stone blocks does not occur in case of low values of pre-compression stress, while the full structural failure occurs due to overturning.However, the cracking of stone blocks does occurs in case of greater precompression stress values, which imposes the need to take into account possibilities of occurrence and development of cracks in stone blocks in this kind of analysis.This is why this possibility was taken into account in this example (Figure 15d).Nonlinear analysis of engineering structures by combined finite-discrete element method

Reinforced-concrete beam exposed to impact load
The possibility of using numerical model based on the FEM/DEM method in the analysis of reinforced-concrete structures exposed to impact load is presented in this example.The simple supported reinforced-concrete beam exposed to projectile impact was selected as an example.Geometrical properties of the beam and projectile are presented in Figure 16.The projectile weighing 37.5 km was exposed to initial velocities of v 1 =60 m/s and v 2 =80 m/s in numerical analysis.

Figure 16. Geometrical properties of reinforced-concrete beam
Material properties of concrete and reinforcement used in the beam are presented in Table 5.

Table 5. Material properties as input parameters
Previous analyses conducted during use of this method have shown that the value of damping coefficient m adopted in this numerical analysis corresponds to the restitution coefficient of 0.18.The restitution coefficient value was determined on stone samples.It was assumed that the concrete used in this analysis has the same restitution coefficient.The numerical analysis was conducted for three cases: a) unreinforced concrete beam, b) beam reinforced with longitudinal reinforcement, and c) beam reinforced with longitudinal reinforcement and stirrups.
The structure discretization for all three cases is shown in Figure 17.Nonlinear analysis of engineering structures by combined finite-discrete element method Cracks occurring over time are shown in Figure 20 for the simple supported beam reinforced with longitudinal reinforcement and stirrups, for velocities of v=60 m/s and v=80 m/s.The total kinetic energy of the projectile and beam as a function of time is shown in Figure 21.If potential beam energy is neglected, which is possible as there are no oscillations prior to and after impact, then the difference between the initial and kinetic energy represents the energy dissipation during the impact.
It can be seen in Figure 18 that the projectile fully penetrates through the beam in case of an unreinforced beam.If the beam is reinforced with longitudinal reinforcement only, the beam splits above the bottom reinforcing bar (Figure 19), while in case of beam reinforced with longitudinal reinforcement and stirrups great damage to concrete and reinforcement is registered, but the projectile does not pierce the beam (Figure 20).This can also be seen from the kinetic energy diagrams (Figure 21) where the maximum energy dissipation occurs in case of beam reinforced with longitudinal reinforcement and stirrups.

Dry stone masonry bell tower exposed to impact and seismic load
The possibility of using the developed numerical model in the analysis of dry stone structures exposed to impact and seismic incremental load is presented in this example.
The dry stone masonry bell tower structure was selected for the analysis.The geometry and discretization of the bell tower structure is shown in Figure 22.
The dry stone masonry bell tower consists of stone blocks of regular dimensions.Material properties of stone used in numerical analysis are shown in Table 6.The damping coefficient value m adopted in this numerical analysis corresponds to the restitution coefficient of 0.18 as obtained experimentally on stone samples.The coefficient of friction amounting to µ=0.60 is used in the numerical analysis.The fracture energy adopted for tension is typical for granite, and has been taken from literature [44], while the fracture energy for shear is assumed to be similar to that used for tension.The dry stone bell tower exposed to projectile impact at the projectile speed of v=50 m/s is presented in this analysis.The structure exposed to such load in individual time increments is presented in Figure 23.Here we can see the advantage of this model which can simulate impact load and predict the way in which full structure failure will be occurred.This failure depends on the energy dissipation during the impact, initiation and propagation of cracks, and inertia effects of individual parts of the structure.The dry stone masonry bell tower presented in this section is analyzed using the method of incremental dynamic analysis which has been used in recent times for analyzing response of structures exposed to seismic load [45].It is based on incremental increase of load (in this case of the real earthquake amplitude), and it enables one to monitor behaviour of structures exposed to seismic load over time, all the way until failure.This enables analysis of the way in which structure failure has occurred, determination of bearing capacity of the structure, determination of behaviour factor, and monitoring other structure-ductility parameters.The incremental dynamic analysis of dry stone bell tower was spend by exposing the structure to horizontal ground acceleration (Figure 24) which was registered on 15 April 1979 in Dubrovnik in rocky soil during an earthquake whose epicentre was in Petrovac (Montenegro).The accelerogram was first scaled to peak acceleration a g =0.22 g which is valid for Split, and was then gradually increased until the peak acceleration which causes full structural failure.Block displacements for individual peak acceleration values are presented in Figure 25.
Figure 26 shows block displacements over time for peak acceleration of a g =1.50g at which the full structural failure was registered.
Here we can see the advantage of this model which can simulate full structural failure caused by big displacements and rotation of stone blocks subjected to friction force acting in between them.

Conclusions
The possibility of using the combined method of finite-discrete elements in the nonlinear analysis of engineering structures is presented in the paper.Numerical results obtained by the developed model based on the FEM/DEM method are compared in the paper with numerical and experimental results available in literature, for reinforced-concrete and dry stone masonry structures.
In the analysis of the reinforced-concrete beam exposed to monotonically increasing load, numerical results for failure load obtained by this model, as compared to numerical results of the nonlinear analysis software based on the finite element method, present greater accuracy in relation to physical experiment.
In the analysis of dry stone wall exposed to monotonically increasing load, the results obtained by the model based on the FEM/DEM method show a good correspondence with numerical results of the nonlinear model based on the finite element method.In this case, the correspondence of numerical results with experimental results can be considered satisfactory because the experiment was conducted with the dry wall made of natural stone and, in such walls, the influence of irregularities between blocks is significant, and it is very hard to model by numerical procedures.Some possibilities for using the method in the analysis of structures exposed to impact load, and in incremental dynamic analysis, are also presented in the paper.The analysis of unreinforced and reinforced concrete beam has shown that the initiation and propagation of cracks, and the corresponding dissipation of energy, is dependent on the reinforcement method applied.The analysis of the dry stone masonry bell tower exposed to impact load has pointed to the advantages of this model in simulation of impact load, and in predicting of the way in which the full structural failure will be operated.In addition,

Figure 1 .
Figure 1.Contact differential force in the vicinity of points P m and P k

Figure 3 .
Figure 3. Model of a crack for tensile softening, shown in the stress -strain relation complete relation showing the relationship δ-σ c in the mode I can be presented as follows

Figure 5 .
Figure 5. Discretization of RC structure

Figure 6 .
Figure 6.Contact element of reinforcement

Figure 8 .
Figure 8. Geometry of a simple RC beam

Figure 10 .
Figure 10.Diagram of beam mid-span displacement

Figure 11 .
Figure 11.Diagram showing dependence between stress in reinforcement and concentrated force F

Figure 13 .
Figure 13.Schematic view of a stone wall: a) geometry and load; b) discretization of structurecije

Figure 17 .Figure 18 .
Figure 17.Discretization of simple supported beam: a) unreinforced concrete beam; b) beam reinforced with longitudinal reinforcement; c) beam reinforced with longitudinal reinforcement and stirrups Cracks occurring over time are shown in Figure 18 for the simple unreinforced concrete beam, for velocities of v=60 m/s and v=80 m/s.Cracks occurring over time are shown in Figure 19 for the simple supported concrete beam reinforced with longitudinal reinforcement, for velocities of v=60 m/s and v=80 m/s.

Figure 22 .
Figure 22.Schematic view of stone bell tower: a) geometry of structure and projectile; b) discretization of structure