Advanced simulation of damage of reinforced concrete structures under impact

The efficiency of the discrete element method for studying the fracture of heterogeneous media has been demonstrated, but it is limited by the size of the computational model. A coupling between the discrete element and the finite element methods is proposed to handle the simulation of impacts on large structures. The structure is split into sub-domains in each of which the method of analysis is adapted to optimise the modelling of the structure behaviour under impact. The DEM takes naturally into account the discontinuities and is used to model the media in the impact zone. The remaining structure is modelled by the FEM. Proposed combined DE/FE algorithm is implemented in the EUROPLEXUS fast dynamics software and parallelised with MPI formalism. The efficiency of the EUROPLEXUS multi-domain MPI parallel version is tested. La modélisation des structures par éléments discrets est bien adaptée aux problèmes dynamiques mettant en jeu de la fragmentation mais elle est difficile à mettre en œuvre sur des structures de grande taille. Couplée à une méthode aux éléments finis, cette méthode devient performante pour simuler le comportement des ouvrages en béton armé soumis à des impacts autant à l'échelle locale de l'impact qu'au niveau global de la structure. Dans un formalisme multi-domaine, les éléments discrets sont utilises dans la zone localisée de fort endommagement alors que les éléments finis sont employés pour le reste de la structure. L'algorithme combinant éléments discrets et éléments finis est implanté dans le code de dynamique rapide EUROPLEXUS et parallélisé à l'aide du formalisme MPI.


Introduction
To simulate the mechanical consequences of extreme impact loads on reinforced concrete (RC) structures the numerical model has to be capable of representing advanced damage states of the structure until its complete failure. The continuum approaches, such as the finite element method (FEM) well adapted to non-linear analysis of structures before failure, reach their limits when trying to describe macro cracking and fragmentation mechanisms. It is nevertheless possible to generate discontinuities in FEM model by using the element erosion technique (Belytschko & Lin, 1987) but the numerical precision of the whole algorithm becomes difficult to control during the fragmentation process namely for the mass and momentum conservation aspects.
The discrete element method (DEM) is a powerful alternative to the finite element modelling especially when advanced damage states and failure have to be studied. This method does not rely on any assumption regarding where and how a crack or several cracks occur and propagate since the medium is considered as naturally discontinuous. The behaviour of the discrete element model at the structural level results from multiple elementary interactions. The local cohesive laws linking discrete elements have to be identified in order to reproduce the macroscopic behaviour of the material.
Among numerous discrete element techniques proposed in the literature we use the numerical model based on rigid spherical elements of individual diameter and mass allowing a very easy contact treatment (Rousseau, Frangin, Marin, & Daudeville, 2008;Rousseau, 2009a).
To be precise and able to reproduce the macroscopic behaviour of the material, DEM needs a rather fine discretisation of the computational domain. Consequently, DEM cannot be applied to model large industrial structures because of its prohibitive computational cost. Two ways are investigated to reduce the computational time. The first way consists in building combined numerical models. Indeed in most practical cases the accidental-type dynamic loads lead to localised fractures and fragmentation of concrete occurs only in the vicinity of the impact. Far from the impacted area the structure is generally not damaged and continues to behave in an elastic way. A coupled "finite element/discrete element" approach has been developed using DEM only locally to obtain a detailed description of the fracture and fragmentation mechanisms, whereas a standard FEM is applied elsewhere (Frangin, Marin, & Daudeville, 2006;Rousseau, Frangin, Marin, & Daudeville, 2009b;Rousseau, Marin, Daudeville, Potapov, 2010). Another way to improve the performances of DEM consists in parallelising the resolution of the dynamic problem involving discrete elements. Thus, a multi-domain multimodel parallel framework based on MPI (message passing interface) formalism has been implemented in the fast dynamics software EUROPLEXUS (EUROPLEXUS User's Manual) co-developed and used by CEA and EDF in accidental-type simulations. The use of parallel computers allows significantly reducing the computational time (REPDYN Project., 2011). This paper first presents a short theoretical description of used DEM. A numerical algorithm coupling finite element and discrete element models is presented next. Some details concerning the parallel implementation of the combined FE/DE algorithm within the EUROPLEXUS software are then given. Finally, DE/FE simulations of an impact on two RC structures are presented.

Model description
The fundamentals of the used DEM are described in detail in Hentz, Daudeville, and Donze (2004). The discrete elements considered here are rigid spheres of different sizes and masses constituting a disordered poly-disperse assembly with a particular size distribution. It should be noted that the characteristic sizes of discrete elements are not representative of concrete constituents such as aggregates. Indeed, we use a higher scale model that has to reproduce the macroscopic behaviour of concrete ( Figure 1).
The discrete element assembly is required to be isotropic at the beginning of the calculation in order to reproduce the isotropic property of undamaged concrete. A non-uniform discrete element distribution guarantees that, unlike in the case of particles arranged in a regular pattern, no preferential directions exist during the fracturing process.
Translation and rotational motions of discrete elements are obtained by time integration of the second Newton's law involving elements' accelerations and inter-elements forces. For this purpose the explicit central difference scheme is used.

Modelling of concrete
Two interaction modes are considered: a cohesive mode characterised by an interaction distance defined for each element and a contact mode describing interactions between elements non-contacting initially or contacts appearing after destruction of cohesive links.
Cohesion-type interactions between two discrete elements are defined by means of normal K n and tangential K s stiffnesses characterising the elastic behaviour of concrete. Since the strain energy for a given cohesive link depends on the size of the interacting elements, the local interaction stiffnesses are not identical over the sample. Some "micro-macro" relations are used to calculate these local stiffnesses from macroscopic elastic parameters such as the Young's modulus and the Poisson's ratio. Those relations stem from homogenisation models typically used for regular DE assemblies; they have been adapted to take into account both the relative disorder of the DE assembly and of the interaction surface S int depending on the size of interacting elements. Equations (1) show the micro-macro relations applied to determining K n and K s between two elements a and b. D init stands for the initial distance between elements. It is equal before deformation to the sum of the element radii R a and R b . Parameters β, α and γ accounting for the spatial distribution of the discrete elements depend on the sphere packing algorithm used, and have to be identified. All details about parameters identification can be found in Rousseau (2009a).
Non-linear behaviour of the discrete element model has also to be identified, in particular for the fracture threshold and tension/compression dissymmetry. Locally it is defined by two failure criteria (2) associated with softening. The first Mohr-Coulomb type criterion controls shear and behaviour's dissymmetry, and the second accounting for local tractions.
In (2) U i is the friction angle, C o is cohesion, and T is local tensile strength. Local parameters, such as T, C o and softening factor ζ need to be identified from global parameters such as compressive and tensile strengths σ c and σ t and fracture energy G f through the simulation of uni-axial quasi-static tensile and compressive tests. When the tension criterion is reached (Figure 2), the model adopts softening behaviour driven by parameter ζ > 0, which allows reproducing the macroscopic softening observed in experiments. If ζ tends to zero, the behaviour becomes perfectly brittle. If the system is unloaded during the softening phase, some permanent deformations remain. If the distance D between two elements exceeds D max specified by the user, the cohesive link is definitely destroyed, and only contact-type interaction following the classical Coulomb friction constitutive behaviour characterised by a friction angle U c is possible between two considered elements.
When the shear criterion is reached, sliding-type interaction between elements is activated without deleting cohesion, which corresponds to a non-associative flow law ( Figure 2).

Modelling of steel reinforcement
The steel reinforcement bars are also modelled using the DEM formalism. Each longitudinal and transversal bar is modelled by means of aligned single "steel" discrete elements whose diameters correspond to those of the real reinforcement ( Figure 3). The discrete elements of reinforcement are linked to their immediate neighbours through special spring-type links having a bilinear behaviour in traction-compression, bending and shear (Figure 4). The local stiffnesses of links are calculated from the Young's modulus and the Poisson's ratio of steel using the standard theory of beams formulae. Thus, the DE model of steel reinforcement bar is capable of reproducing the classical behaviour of a beam subjected to different types of loading. Behaving elastically until the prescribed elastic limit σ y , the DE bars can undergo plastic deformations with a classical strain hardening. The link is broken in traction if the ultimate stress σ r or the maximum prescribed elongation D max are reached (D max = D init + ɛ max D init with ɛ max the rupture strain).

Steel-concrete link
The concrete-steel interaction laws adopted here are derived from ideas proposed in the finite element context (Domingez, 2005). Thus, the tangential law ( Figure 5, left) reproduces at first a perfect stick of the two materials with an elastic stiffness K s , then a sliding phase followed by damage and finally by a frictional pseudo-contact with a residual effort. To calibrate the parameters, the push-out test is used. The normal law is much simpler (Figure 5, right). It describes an elastic brittle behaviour in traction and elastic perfectly plastic behaviour in compression.

Discrete element-finite element coupling
To optimise the numerical model, discrete elements used in the vicinity of the impact are coupled with finite elements used in the rest of the structure using the bridging domain method proposed by Xiao and Belytschko (2004) to couple FEM with the     A weighted parameter η is used whose value is 1 in the FE domain and 0 in the DE domain ( Figure 6). To guarantee kinematical continuity, the degrees of freedom of both domains within the interface zone are linked using Lagrange multipliers. The coupling algorithm implemented in EUROPLEXUS allows linking discrete elements with 3D and shell-type finite elements (Rousseau et al., 2010).
Discrete element and finite element models can also interact by contact, which is detected using Pinball method (Casadei, 2003). Kinematical constraints are solved using Lagrange multipliers.

Geometrical support for DEM
The geometrical support of DEM model is different from the classical finite element meshes, thus common mesh generation software cannot be used. To generate polydisperse DEM sphere packings with a particular size distribution we use an algorithm described in Jerier, Richefeu, Imbault and Donze (2010). This algorithm uses a new geometric procedure to pad very efficiently the polydisperse spheres in a given tetrahedral mesh (Figure 7). The padding procedure is stopped when a target solid fraction or number of spheres is reached. Implemented initially into YADE (2011) open-source software, this algorithm has been recently inserted into the SMESH mesh generator of the SALOME (2011) open-source platform which was co-developed and used by EDF and CEA for pre-and post-processing of numerical simulations.  A special function of this computational tool allows generating DEM mesh for reinforced concrete structures (Figure 8).
The quality of the DE meshes obtained with the geometric packing algorithm is satisfactory. Figure 9 shows a typical size distribution obtained for a DE sample.

MPI parallel resolution for DE model
To accelerate the resolution of a dynamic problem involving discrete elements, the multi-domain parallel framework of EUROPLEXUS has been extended to account for the discrete element model and its coupling with the finite element model.
Parallel acceleration is based on domain decomposition: the whole numerical model is partitioned into geometric sub-domains, which are assigned to different processors ( Figure 10). Calculation of internal forces and identification of kinematical constraints is thus distributed on different processors available. One should note the specificity of parallel implementation of the discrete element method when compared with the finite element method, which comes from the fact that the computational domain for DE model starts changing with the beginning of the fragmentation process. Indeed, during the calculation some inter-element cohesive links are broken, creating macro-cracks, whereas some new contact-like links are generated to account for crack closure or contacts between fragments created.
When one structural part is modelled with discrete elements and another one with finite elements, interaction between the two models takes place in the recovering zone only. To optimise parallel implementation and avoid additional inter-processor communications, domain decomposition is modified by forcing any discrete element interacting with a finite element to belong to the sub-domain of concerned finite element (Figure 11).
To manage parallel resolution process and update data at each time step, each processor must access all necessary data locally. Thus each calculation sub-domain uses a list of discrete elements of the sub-domain and a list of their links. A specific treatment Figure 10. Example of subdomains definition for the problem of impact of finite element projectile on a discrete element slab. Figure 11. Subdomains in a combined finite element-discrete element model. Figure 12. Definition of a "crossover" neighbourhood for a near-to-border discrete element (large blue circle in centre). must be applied to the discrete elements located in the vicinity of one sub-domain's boundary, since some of their neighbour elements belong to other sub-domains ( Figure 12). Concerned discrete elements are called remote discrete elements for the current sub-domain and inter-processor communications are required in order to identify new links for elements of neighbouring sub-domain and calculate correctly the corresponding inter-element forces.
If computational domain does not change a lot during simulation, extensibility is very good. Otherwise, if many inter-element links are destroyed, the communications between sub-domains become too costly and dynamic domain repartitioning is needed to maintain the resolution efficiency (REPDYN, 2011).

Applications
Two calculations of an impact on RC structures have been simulated with an MPI version of the software EUROPLEXUS (2011).
The first calculation deals with an impact of a rigid projectile on a RC beam with a rectangular cross section and simply supported at its ends. Impact velocity is 8.29 m/s. Reinforcement consists of four bars situated near cross section angles. DE mesh is composed of 47.878 discrete elements. In the second calculation we simulate test no. 5 of the MEPPEN test campaign dealing with impact of soft missiles on RC slabs (Jonas, Meschkat, Riech, & Rüdiger, 1979). The slab of 6.5 by 6 m and 70 cm thickness is impacted by a tube-like missile (6 m length, 60 cm diameter) at 235 m/s. DE model of the slab contains 256701 elements with 193.000 "concrete" and 63.000 "steel" elements. The projectile is modelled using FEM. Figure 13 shows evolution of damage obtained in the first simulation for different time stations. In the DE context the local level of damage is characterised, for a given discrete element, by the number of broken links with the neighbouring elements. Thus, the dark colours correspond to high damage levels. Qualitatively, the results of the numerical simulation are satisfactory: one can observe formation of a typical conical shape of damage corresponding to a particular perforation mode in form of a shear plug Figure 13. Damage maps during the impact of a rigid projectile on a RC beam. delimited by oblique cracks. The same perforation mode is visible on Figure 14 for the RC slab, as observed on the real slab of the MEPPEN test.
To measure the efficiency (defined as the ratio between speed-up and number of processors) of EUROPLEXUS MPI parallel version, computations were performed with 1, 2, 4, 8, and 16 processors. A satisfactory efficiency is obtained until eight processors are used, then acceleration is saturated due to a small number of elements per subdomain with respect to communication operations needed (Table 1).
For the second example, efficiency of the parallel algorithm decreases progressively because many inter-element links are destroyed, and the communications between subdomains become too costly, limiting extensibility to 16 processors. It is necessary to make a new partitioning of the domain to maintain resolution efficiency.

Conclusion
The proposed combined DE/FE algorithm implemented in EUROPLEXUS software and parallelised with MPI formalism opens the way to simulate large reinforced concrete structures under impact. Computational performances of the multi-domain parallel framework of EUROPLEXUS can still be improved by a better load balancing of processors. To achieve this goal, implementing a dynamic domain partitioning approach is suitable. First, this requires the capability of switching from one domain decomposition to another with memory transfers associated to migrating elements. Second, criteria must be provided to assess the updates of domain decomposition to runtime internal load balancing measures. Finally, multi-time step integration can be implemented to improve efficiency for both sequential and parallel algorithms when coupling between finite elements and discrete elements is used, since the time scale associated to finite elements is classically greater than that associated to discrete elements.