Three-Dimensional Finite-Discrete Element Framework for the Fracturing of Reinforced Concrete Structures

This paper presents a three-dimensional numerical model based on finite-discrete element method for simulating cracking and predicting failure in reinforced concrete (RC) structures. Discrete representation of concrete cracks is coupled with a non-linear reinforcement model where reinforcing bars interact with concrete and slip due to crack opening and steel plastic deformation. Hence, it enables a simulation of complex material behaviour of concrete and steel in cracked zones and its influence on global structural behaviour. Several numerical examples are used to study the sensitivity of the model to different numerical and physical parameters and its capabilities in the analysis of RC structures.


INTRODUCTION
The structural members in RC structures are often exposed to tri-directional stress state, therefore an analysis of these structures demands three-dimensional (3D) numerical models, both in design and in assessing the state of the existing buildings. Due to the complexity of the material behaviour of concrete and steel in non-linear range, especially crack opening and closing mechanism in concrete, accompanied by yielding, pulling out and buckling of the reinforcement, 3D numerical modelling of these structures is challenging, both for static and for dynamic loads.
The most commonly used computer framework, which takes into account the discontinuous nature of the concrete after cracking, is based on the finite element method (FEM) and smeared representation of the cracks [1,2], usually combined with classical plasticity or damage constitutive models. The advantages of this continuum approach are primarily in simplicity, low computational cost and relatively good approximation of the global structural response; hence, it represents a reliable choice in the nonlinear static analysis up to failure where the cracking propagates gradually with the increasing of the load. On the other side, RC structures exposed to extreme loadings, such as earthquakes, impact and explosions, exhibit rapid damaging and progressive failure of structural elements. In this case, the discrete crack approach provides a more precise description of the discontinuities in the concrete, which is essential for the simulation of energy dissipation due to the cracking process.
Several computational frameworks are available in the discrete modelling of the crack propagation. One is derived from the FEM, and the others are related to the discrete element method (DEM). Cohesive elements across element edges [3,4] and adaptive refinement techniques [5] are examples of discrete crack modelling coupled with FEM. Complexity of these approaches motivated researchers to find other possibilities in the description of discontinuities, like ED-FEM [6,7] or X-FEM [8,9] methods. These approaches have been adapted in quasi-static analysis of RC structures [10][11][12] and in the dynamic fracturing problems [13][14][15]. Microplane models [16] where the stresses are monitored in different predefined directions, can be also used in simulation of the damage caused by cracking.
Contrary to the FEM that starts from the idealisation of the structures as a continuum, the DEM is based on the modelling of the structure as discontinuum. DEM was initially developed in granular material problems, but later it was applied in the fracture simulation of rocks [17,18] and concrete [19]. Discretisation of the structures is given in the form of particles, blocks or Delaunay/Voronoi tessellations, while the connections between the discrete elements were achieved with contact (joint, interface) elements [18], lattice elements [20][21][22][23][24][25][26] or rigid body springs [27,28]. Discrete elements are assumed as rigid or deformable bodies. Their deformability can be simulated by adding strain deformation parameters or by discretising them into finite elements. Latter approach results with numerical models that combine the advantages of finite and discrete element approaches [29][30][31]. Recently, combined FDEM [32] was used in 2D modelling of the fracture and fragmentation process of quasi-brittle materials [33][34][35], but few models of FDEM were also developed for solving the 3D fracturing with possible applications in concrete, rocks and stone masonry [36][37][38]. Most recently, it was also recognised that FDEM offers useful numerical tools for solving fracture problems of RC structures. Our team applied the method in 2D analysis of the classical RC structures [39][40][41][42] reinforced by steel reinforcement. FDEM was also used in 2D modelling of concrete with glass fibre polymer bars [43]. Spherical DEM approach combined with FEM and embedded rebar reinforcement elements was applied in the damage analysis of RC structures [44].
The objective of this study is to present a 3D model for RC structures within the framework of combined FDEM [32]. The fundamental ideas established for 2D analysis of RC structures [39][40][41][42] by our team have been extended in order to simulate their 3D behaviour. The discrete representation of concrete cracks, as a main strength of the FDEM model, was coupled with the non-linear reinforcement model, where the reinforcing bars interact with concrete and slip due to crack opening and high plastic deformation of the steel. After the transition from the elastic zone to the discrete fracture, the contact detection and interaction algorithm is used to simulate the interaction between fracture surfaces, represented through the normal compression and sliding friction. This very important mechanism enables simulation of opening and closing of the cracks during the cyclic loading of the structures.
Although the paper presents the main features of FDEM in the modelling of 3D RC structures, special emphasis is given to the investigation of the influence of mesh refinement, penalty parameter and crack spacing on the accuracy of the solution and verification of the model in various cases of 3D RC structures.
The paper is structured as follows. Section 2 briefly introduces a 3D concrete fracture model with embedded reinforcement in the framework of the FDEM. Formulation and algorithms, which are essential for coupling of reinforcement and concrete, are given through several phases, starting with the continuum modelling of uncracked RC structures, crack opening in the concrete and its influence on the behaviour of the reinforcement in discrete cracks. Sensitivity analysis of the model to various physical and numerical parameters together with applications in analysis of RC beam under cyclic load and monotonically increasing load are provided in section 3.

THREE-DIMENSIONAL CONCRETE FRACTURE MODEL WITH EMBEDDED REINFORCEMENT
Three-dimensional fracture model for RC structures, presented in this study, was derived as a coupled problem connecting two materials with quite different properties. In fact, when the RC structure is exposed to increasing load, due to the low tensile strength, the concrete passes through several different phases, starting with the continuum, microcracking and finally crack opening at the macrolevel. At the same time, the stress in the reinforcing bars placed in the cracking zones increases, initially in the linear range and later in the non-linear range, due to the significant deformation of the bars inside the opened cracks and pulling out of the bars from the surrounding concrete. The bars hold the separated parts of concrete together until the deformation in the steel achieves the critical value, leading to the breaking of the reinforcement.
FDEM enables the simulation of aforementioned effects in the concrete by using finite elements and joint elements implemented within the finite elements [32]. Joint elements enable the cracking of the concrete, managing through the penalty method and contact detection and interaction algorithm, which interact with each other. There are two main regimes: (1) uncracked, where the penalty serves for maintaining compatibility between the elements, and (2) cracked, where the penalty is used for computing interaction forces between the separated elements. In both regimes, the contact algorithm is needed to detect either the neighbouring elements (uncracked regime) or crack the separated elements which interact with each other (cracked regime). The penalty method and the contact interaction algorithm allow the simulation of opening and closing of the cracks in the joint elements, which makes the method applicable in modelling, not only of the structure under the static loading, but of the cyclic loading and unloading as well.
The presented model, derived for the 3D RC structures, uses discretisation by 4-node tetrahedron finite elements with embedded 6-node cohesive joint elements for the simulation of concrete behaviour (Fig. 1). Reinforcement is represented by the reinforcing bars consisting of 2-node line finite elements and 2-node line joint elements. The intersection between the axis of the reinforcing bar with the sides of the concrete tetrahedron finite elements defines the nodes of the reinforcing bar finite elements. The part of the reinforcing bar between the two reinforcing bar elements belongs to the reinforcing joint element (Fig. 1). As the joint elements of concrete allow the opening and the closing of the cracks, similarly, the reinforcing bar joint elements allow the plastic deformation and the breaking of the reinforcement within the cracks. Hence, in the case of the modelling of RC structures, FDEM combines the FEM formulation for representing the behaviour of the concrete before the fracture initiation and reinforcement with the DEM formulation aiming to model fracture process in the concrete and the reinforcement. These mechanisms occur in the finite and joint composite elements, both consisting of concrete and reinforcement elements. In terms of mathematical formulations and applied numerical algorithms, the 3D concrete fracture model intended for the analysis of RC structures contains the following features: • Discretisation of the structure is performed by 3D finite and joint elements intended for the modelling of concrete behaviour and 1D finite and joint elements for the reinforcing bars. Concrete failure criteria generate discrete cracks in the concrete, represented with concrete joint elements, which lead to high plastic deformation and the breaking of the reinforcing steel. Material non-linearity of the concrete and reinforcement is considered in concrete and reinforcement joint elements.
• Material models consider both static and cyclic material behaviour. Their implementation in the finitediscrete element framework, based on dynamic equilibrium equations, enables the simulation of the behaviour of RC structures for quasi-static and dynamic loading.
• Discretised equilibrium equations of motion are given in the form: where M refers to mass diagonal matrices, u  to the acceleration vector of the nodes, F n to the vectors of nodal forces which include the contribution of reinforcement and concrete, contact interaction algorithm, damping introduced into the model and external loading. Therefore, the vector of nodal forces F n is given in the form: is the force vector generated by the contact interaction algorithm which includes the normal compression and the sliding friction. damp F is the nodal force vector generated by introducing the damping into the model or to model quasi-static phenomena by dynamic relaxation. Finally, ext F is the external nodal force contributed by structural loading, consisting of body forces and surface tractions. An explicit time integration scheme is used to solve the discretised equations of motion (1).

Continuum Modelling of Uncracked RC Structures
Concrete structure is discretised by 4-node tetrahedron finite elements in the area without reinforcement, whereas composite finite concrete-reinforcement elements are used for the discretisation when the reinforcing bars cross the tetrahedron concrete elements. Composite finite elements consist of a 4-node tetrahedron concrete element and a 2node reinforcing bar element (Fig. 2).
By taking into account the Young's modulus of elasticity of steel E s , the stress in the bar is obtained as: This yields the forces in the nodes of the reinforcing bar finite element (Fig. 3a) which are given by: where A s is the cross-sectional area of the reinforcing bar.
Forces f 0se and f 1se acting in points P 0 and P 1 are distributed as equivalent nodal forces into the nodes of the corresponding concrete tetrahedron finite element, as shown in Fig. 3b. Strain in the tetrahedron concrete element is calculated from the global coordinates of each node in the initial configuration (x i , y i , z i ) and the current configuration (x c , y c , z c ). Displacements of the element result from the translation and the rotation of the solid body and its shape and volume deformations. Deformation of the element is derived from the deformation gradient F: Green-St. Venant strain tensor E  can be calculated as follows: while Cauchy stress tensor T is obtained according to: where E is the modulus of elasticity, ν is the Poisson's ratio, εν is the volume deformation expressed as: µ is viscous damping coefficient and D is the rate of the deformation tensor [32]: where L is the velocity gradient expressed as: The nodal forces from the corresponding surface of tetrahedron finite element are obtained as:  12) for each of the four surfaces of the tetrahedron and equivalent nodal forces resulting from the 2-node reinforcing bar finite element.

Discrete Crack Model of Concrete
Fracture initiation and propagation of the cracks in the concrete is modelled by 6-node joint elements inserted between the adjacent 4-node tetrahedron elements. The normal and shear separation (δ, t) of the surfaces of adjacent tetrahedron concrete finite elements induces normal and shear bonding stresses. The normal and shear springs are used to model this phenomenon (Fig. 4b). The material model is shown in Fig. 4a, where Gf is energy release rate and defines the area under the softening part of stress-displacement curve. Normal and shear spring model The penalty method [32] ensures compatibility before reaching the tensile strength. For the separation δ ≤ δ t , the bonding stress is given by: where δ t = 2hf t /p is the normal separation which corresponds to the tensile strength f t , h is the size of the finite concrete element, and p is the penalty parameter.
After reaching the tensile strength f t , the stress decreases with an increasing of the normal separation δ, whereas at δ = δ c the bonding stress tends to zero. For the separation δ t < δ < δ c , the bonding stress is given by: where z is the scaling function representing the softening behaviour of the concrete. It is used according to Hillerborg [45] in order to approximate the experimental stressdisplacement curves for the concrete: where a = 0,63, b = 1,8 and c = 6,0, while the damage parameter D t is determined according to the expression: The faces of the two adjacent elements are held together by the shear stress. After reaching the shear strength f s , the stress decreases. For sliding t s <|t|<t c , the shear stress is represented by: where z is given by Eq. (15), but with damage parameter D s , instead of D t , expressed by: If the normal and shear separations act simultaneously, the tension and shear damage parameter D t and D s in Eq. (15) are replaced by:

Reinforcement Model in Discrete Cracks 2.3.1 Reinforcement Strain-Slip Relation in the Crack
Numerical model of the reinforcing bar joint element is based on the assumption that the separation of the adjacent bar finite elements induces the stress in the bar joint element (Fig. 5a). Considering that there is no separation of the nodes of the reinforcing bar joint elements before reaching the tensile strength, the model of the reinforcing bar in the joint element is divided into parts before and after reaching the tensile strength of concrete (Fig. 5b). Before reaching the concrete tensile strength, the reinforcing bar joint element ensures the continuity of displacements of the reinforcing bar and the concrete finite elements. After reaching the tensile strength and the crack initiation in the concrete, the influence of reinforcing bar axial force on concrete between the adjacent cracks is implemented through the bonding between the reinforcement and the concrete. Therefore, the local stress along the bar and stress at the interface are different. That is why strains along the reinforcing bar are not uniformly distributed, and, inter alia, depend on the bar pull-out from the interface of crack S (Fig. 6a). To model that, mechanical model for the deformed reinforcing bar at the RC interface developed by Soltani and Maekawa [46] was used as the basis for the numerical model of reinforcing bar joint element (Fig. 6b). The model assumes that the longitudinal cracks propagate along the bar causing bond deterioration between the bar and the concrete near the crack plane [46,47].
Unbonding between the reinforcement and the concrete causes the local slip S of the bar from the crack interface. The non-dimensional slip s of the bar is obtained according to: where D is bar diameter and f c is compressive strength of concrete (MPa).
Expressions used in presented model for slip-strain relation in pre-yield and post-yield steel region [39] are presented in Fig. 6., where ε s presents the strain at the reinforcing bar in the crack, ε y and ε sh are the yield strain and strain at the onset of hardening of the bar respectively, while f u and f y are tensile strength and yield stress of steel (MPa). Under the assumption that the distribution of strain in the yield zone is linear (Fig. 7), the normalised steel slip s pl is expressed as: where ε max and s max are steel strain and non-dimensional slip instantly after the change from the loading to the unloading, β is the factor obtained from experiments and is approximately equal to 1,0. By substituting the equation (22) in (21), the strain in the reinforcing bar at the crack can be obtained from the known non-dimensional slip s.

Influence of Shear Stress in Reinforcement to Yield Stress Reduction
Although the reinforcing bar in the RC structures primarily adopts axial forces, their deformation within the cracks also induces shear forces due to the deflection and curvature of the bar near the face of the crack (Fig. 8).

Figure 8 Deformation of reinforcing bar within discrete crack
This phenomenon was investigated in the work of Qureshi and Maekawa [48]. Considering the curvature distribution within the curvature-influencing zone, the length of the curvature-influencing zone and the local slip of the reinforcing bar [39], the shear force in the reinforcing bar within the RC faces is given by: where L c presents the length of the curvature-influencing zone while E s is the Young modulus and I s is the moment of inertia of the bar respectively. The shear stress in the reinforcing bar influences the reduction of the yield stress in the bar, which is computed according to the Von Mises yield criterion: where σ y ' is the reduced yield stress, τ s = V s /A s is the shear stress of the reinforcing bar between the crack surfaces, and A s is the cross-sectional area of the reinforcing bar.

Effect of Adjacent Cracks to Steel Slip
Strain in reinforcement depends on the slip of the bar from the crack interface, as mentioned earlier. Developed formulations consider bond deterioration between the concrete and the reinforcement in an implicit way [46,47]. In order to achieve realistic results, the unbonded length of the bar could not exceed the total bar length. This demand can be ensured by accounting the influence of the adjacent cracks through the reduction factor α depending on the distance lcr [46] (Fig. 9). The steel slip s cr is expressed as: where s is defined by expression (21), and presents the nondimensional slip. The reduction factor α is: In FDEM model, finite element edge defines the position of the crack. Therefore, in this model l cr is an input parameter and it is equal to h/2, where h is the average size of the concrete finite element. Thus, the unbonded length of the bar is limited by the distance l cr between the cracks, i.e. average size of the concrete finite element.

Material Model of Steel
The adopted steel material model for the monotonic and the cyclic loading is shown in Fig. 10. Steel bar in RC structures subjected to cyclic loading exhibits hysteresis behaviour enforced through the Kato's stress-strain model [49]. Stress-strain relations are expressed for the cycles of unloading, negative loading, reloading-unloading and reloading (curves 1, 2, 3, 4) in the following form: , E s , f y and f u are modulus of elasticity, yield stress and tensile strength of the steel, ε s and ε sh are strain at the reinforcing bar in the crack and strain at the onset of hardening, while ε y , ε u are yielding and ultimate strain of steel and σ pm is minimum value of σ s .

Figure 10
Stress-strain model of steel

SENSITIVITY ANALYSIS
In this section, sensitivity analysis of the model to mesh size, crack spacing and penalty parameter has been carried out on the examples of concrete and RC structures. Considering that FDEM starts with the discontinuum representation of uncracked RC structures, maintaining compatibility between the tetrahedron finite elements is essential for the reliable modelling of the structures as a continuum. Therefore, the influence of the penalty parameter on the solution in the linear elastic range was analysed on several examples for both concrete and RC beams. As the cracks in FDEM model are predestined with the finite element mesh, the sensitivity analysis of the model to crack spacing and mesh refinement was also investigated.

Sensitivity Analysis of the Model to Mesh Refinement and Penalty Parameter in Tension
This example was selected to analyse the sensitivity of the model to mesh refinement and penalty parameter p. Concrete beam (Fig. 11) was exposed to monotonic increasing tension load represented by a constant velocity v = 1,25×10 −4 t m/s at the beam ends. Modulus of elasticity of concrete was equal to E c = 30 GPa. The analyses were first performed without and with joint elements using two different finite elements mesh refinement (mesh A, mesh B), as shown in Fig. 12. Mesh A comprised 240 finite elements, whereas mesh B comprised 1920 finite elements. In the analyses performed with joint elements, the penalty term was p = 20 E c . The comparison of the analytical and numerical results is presented in Fig. 13. In the case when analyses were performed without the joint elements, the numerical results corresponded to the analytical solution. The introduction of joint elements led to the relative error of numerical results in comparison to the analytical solution. In both cases (with and without joint elements), the numerical model is not sensitive to mesh refinement.
The sensitivity analyses of the model in the linear elastic stage to the penalty parameter p were performed for three values of the penalty parameter (p = 20Ec, p = 60E c and p = 100E c ) with discretization presented in Fig. 12a. The comparison of analytical and numerical results obtained by the presented model is shown in Fig. 14. The influence of the penalty parameter on the relative error is shown in Table 1. It is evident that the relative error decreases with the increasing of the penalty. The values of the penalty parameter higher than 100E c provide a relative error below 3%. These results are very close to those obtained in the 2D case [40]. Of course, by increasing the penalty parameter the solutions are more accurate but require much smaller time steps to ensure numerical stability [32]. Hence, there is a need to make a reasonable choice with the aim of achieving the smallest possible error and the longest possible time step.

Sensitivity analysis of the model to mesh refinement and penalty parameter in bending
Similarly to the previous case where the sensitivity analysis of the model was performed on the concrete beam dominantly exposed to tension load, in this example the influence of mesh refinement and penalty parameter to the solution error was analysed for the beam with dominant bending. The analysis was performed on a simply supported concrete beam (Fig. 15) subjected to gravitation load. The modulus of elasticity and density of concrete was equal to E c = 30 GPa and ρ = 2500 kg/m 3 , respectively. Three meshes were used for discretisation of the beam (Fig. 16). Mesh A is characterised with finite elements equalling h = H/4, while meshes B and C are characterised with finite elements equalling h = H/6 and h = H/8, respectively, where H is the height of the cross-section. The beam oscillates due to its self-weight and subsequently, finds an equilibrium due to damping (Fig. 17). Relative errors of the numerical results obtained with the proposed model are compared to the theoretical solution for Timoshenko's beam [50] and shown in Tab. 2. The beam's midspan deflection for the beam with shear factor k = 2/3 [50] is equal to 0,729 mm. It is observed that, as the number of finite elements increases, the numerical solution converges to the theoretical ones. The analyses show that the minimum eight finite elements at the height of the beam are required to obtain an acceptable numerical error for cases with dominant bending influence. This finite element mesh pattern is adopted as a basis for further analyses.
The sensitivity analyses of the presented model to the penalty parameter for structures with dominant bending were performed for the penalty values p = 20Ec, p = 60E c and p = 100E c and finite element mesh C (Fig. 17c).
The influence of the penalty parameter on the relative errors of the model with joint elements with regard to the solution for Timoshenko's beam [50] is presented in Tab. 3. The analysis shows that the displacement error can be controlled by determining penalty parameter p depending on the elasticity modulus E c . The analysis also shows that the relative error is below 3% when penalty parameter is p = 100E c . By further increasing the penalty parameter, the error was decreased.

Sensitivity of the reinforcing bar model to the Penalty Parameter for Tension
The sensitivity of the reinforcing bar model in the linear elastic range to the penalty parameter p was analyzed on the RC beam shown in Fig. 18.  Strain at the onset of hardening, ε sh = 0,03 Ultimate strain, ε u = 0,1 Strain at the fracture, ε br = 0,12 Penalty parameters (p = 20E c , p = 60E c and p = 100E c ), previously selected in the sensitivity analyses in earlier examples, are implemented here as well. The comparison of analytical and numerical results obtained by the presented model is shown in Fig. 19.
Obtained results presented in Tab. 5 reaffirm that the penalty above 100E c produces a relative error below 3%.

Figure 19
Force-displacement relations -sensitivity to penalty parameter

Sensitivity Analysis of the Model to Mesh Refinement and Crack Spacing
The law of interaction between the concrete and the reinforcement has changed, starting from the perfect bond in uncracked concrete structures to the slip of the reinforcing bar after the opening of the crack. The slip is a function of the reduction factor α (Eq. 26) and depends on the distance between the adjacent cracks lcr and the bar diameter D. As the faces of tetrahedron finite elements predetermine the cracks in the model, the RC beam under monotonic increasing tension load, with the same geometry and material characteristics as in example 3.3, is used to analyse the mesh and crack spacing sensitivity. The constant velocity v = 0,1 m/s at the ends of the RC beam (Fig. 18a)   The influence of the adjacent cracks is considered through the parameter α depending on the crack spacing l cr = h/2. The relationship between the average stress σ=F/A and the average strain ε=Δl/l is shown in Fig. 21. Fig. 22 shows failure patterns for the analysed meshes. The comparative analyses reveal that the used concrete-reinforcement bond model is not significantly mesh sensitive. Hence, the average stress -strain relation is independent of the distance between the cracks. This conclusion applies under condition that approximately linear relationship of α and l cr /D exists, as it is valid for l cr /D ≤ 10 (see Eq. (26) and Fig. 9). This requirement limits

VALIDATON OF NUMERICAL MODEL 4.1 RC Beam Subjected to Cyclic Load
The ability of the model to simulate structural response under cyclic loading conditions was verified on the RC beam experimentally tested by Maekawa et al. [47]. Geometric characteristics of the beam and cyclic loading history are presented in Figs. 23 and 24, respectively, whereas the material characteristics are provided in Table  6. The position of the cracks in the experiment was predestined by notches, placed at a distance of 30 cm.  Ultimate stress, f u = 540 MPa Strain at the onset of hardening, ε sh = 0,0165 Ultimate strain, ε u = 0,1 Strain at the fracture, ε br = 0,019 Two finite element meshes are used to simulate the beam's behaviour (Fig. 25). In the coarse mesh, vertical faces of tetrahedrons are placed at the positions of the notches, therefore the only possible vertical cracks correspond to those experimentally obtained. In the fine mesh, the crack opening was ensured by reducing the tensile strength, at the position of the notches, in the concrete joint elements. The crack patterns of the beam for meshes correspond to those obtained by experiments (see Fig. 26). The relationship between the non-dimensional slip s and the strain ε in the reinforcing steel at the crack for both discretisations is shown in Fig. 27.

Simply Supported RC Beam Exposed to Monotonically Increasing Load
A simply supported RC beam (Fig. 29) was selected to validate the model for monotonically increasing loading.   The discretisation of the structure was performed by 6912 tetrahedron concrete finite elements and by 216 reinforcing bar finite elements. The load was applied incrementally until the collapse. Fig. 31 shows the discretisation of the structure. The FDEM results were compared to those obtained by experiments [51] and by the computer 3D nonlinear FEM program PRECON3D [52,53] for RC structures based on smeared crack approach. Fig. 31 presents the comparison of the mid-span loaddisplacement curve until failure, where a very good agreement between the presented numerical model (FDEM 3D) and results of experiment is obtained. The loaddisplacement curve shows that the largest deviation of results obtained by FDEM 3D in comparison with the experiment is 10,4%, whereas the failure load obtained by the presented model is 1,07% lower than the experimental result. Both numerical models, FDEM 3D and PRECON3D, exhibit a very good match of failure load with the experiments, whereas the FDEM3D model gives a significantly lower error for the mid-span beam displacement at the failure compared to the solution obtained by PRECON3D. In fact, the softening of the concrete in joint elements coupled with the deformability of joint elements of the concrete and the reinforcement, together with the appropriate modelling of the slip of reinforcing bars, provides a realistic description of the structural deformability leading to the solution of the displacement which is very close to the experimental ones. Fig. 32 shows the changes of the tensile stress in the reinforcement in the mid-span of the beam. The increase of the load caused the reaching of the tensile strength in the concrete joint element and the crack openings accompanied by the significant increases of stress in the reinforcement. The steel yielding stress was achieved for the load f = 56 kN and it increased up to the failure. The potential of the presented model is not only exhibited in the precise description of the loaddisplacement curve, failure load and displacement, but also in the modelling of the failure patterns. The cracking starts at the bottom of the beam and propagates to the top causing the transfer of the stress in the cracking zones from the concrete to the reinforcement. Fig. 33 shows the development of the cracks in the concrete just before the final collapse. The collapse of the beam occurred due to the yielding of the reinforcement, which enabled notable cracking and separation of the concrete. The failure pattern of this beam is typical for four-point bending of RC beam with inadequate flexural capacity, i.e. insufficient amount of reinforcement in the tensile zone leads to the failure of the beam over the reinforcement. Due to the constant value of the bending moment and zero value of shear forces between two loading points, the whole area between these points has constant normal stress, while the shear stress does not exist. Therefore, vertical flexural cracks in the whole area between two external forces occur, firstly in the bottom, and with the load increasing, they elongated to the top of the beam. Left and right of the applied external forces, at the locations of the beams where both the bending moment and the shear stress are high in magnitude, flexural shear cracks are formed as a result of exceeding the diagonal tensile stress due to the combination of significant normal and shear stress.

CONCLUSIONS
A complete 3D FDEM formulation for the fracturing of RC structures has been proposed, formulated and implemented in the open-source FDEM software. The discrete representation of concrete cracks, as a main strength of the FDEM model, was coupled with the nonlinear reinforcement model, where the reinforcing bars interact with the concrete and slip due to the crack opening and high plastic deformation of the steel. Concrete fracturing, reinforcement nonlinearity up to the breaking of the bar, and the slip of the bar from the crack faces were modelled in joint elements of concrete and reinforcing bars by using the approximation of the experimental curves for the concrete and the reinforcement in the interface of crack.
The solution of the presented model is provided in a dynamic framework by using an explicit integration scheme. Therefore, the model enables the simulation of the behaviour of RC structures, both for the static and the dynamic load.
As the cracks in the FDEM model are predetermined with the finite element mesh, special emphasis in the paper is given to the sensitivity analysis of the model to different mesh sizes, penalty parameters and crack spacing. The analyses reveal that the model is not considerably mesh sensitive if the ratio between the distance of cracks lcr and the diameter of bar D meets the requirement l cr /D ≤ 10. This requirement limits the maximum tetrahedron finite element size h max to 20D in order to avoid the mesh size sensitivity.
The performed analyses demonstrated that the presented model has a potential to capture complex behaviour of RC structures, both in predicting the loaddisplacement path and to estimate the failure pattern. However, additional validations are required to test the model for various load conditions. In fact, this paper was primarily focussed on the validation of the 3D model for monotonic increasing load and cyclic load, however, considering the original intention of FDEM to capture dynamic structural response and the fact that the cyclic loading-unloading material models were already built in the code, the testing of the presented model under dynamic load is one of the following aims.