# Advanced simulation of damage of reinforced concrete structures under impact

S. Potapov, Vincent Faucher, Laurent Daudeville

► To cite this version:

S. Potapov, Vincent Faucher, Laurent Daudeville. Advanced simulation of damage of reinforced concrete structures under impact. European Journal of Environmental and Civil Engineering, 2012, 16 (9), pp.1090-1101. 10.1080/19648189.2012.699747. hal-02002441

HAL Id: hal-02002441

https://hal.univ-grenoble-alpes.fr/hal-02002441

Submitted on 14 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Advanced simulation of damage of reinforced concrete structures under impact

S. Potapov<sup>a*</sup>, V. Faucher<sup>b</sup> and L. Daudeville<sup>c</sup>

<sup>a</sup>EDF R&D, 1 av. du Général de Gaulle 92140 Clamart, France; <sup>b</sup>CEA, DEN, DANS, DM2S, SEMT, DYN, Gif-sur-Yvette F-91191; <sup>c</sup>UJF/INPG/CNRS, 3SR Lab, DU BP53, 38041 Grenoble Cedex 9, France

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.

**Keywords:** discrete elements; finite elements; coupling; fast dynamics; parallelism

**Mots-clés:** éléments discrets; éléments finis; couplage; dynamique rapide; parallélisme

## 1. 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

<footnote>*Corresponding author. Email: serguei.potapov@edf.fr</footnote>
numerical precision of the whole algorithm becomes difficult to control during the frag-
mentation 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 multi-model 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.

## 2. 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.
<figure><img src="image_2.png" /><figcaption>Figure 1. Considered modelling scale.</figcaption></figure>

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.

### 3. 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 $\beta$, $\alpha$ and $\gamma$ 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).

$$ \left\{ \begin{aligned} K_n &= E \frac{S_{int}}{D_{init}^{\frac{1+z}{z}} \beta(1+v) + \gamma(1-zv)} \\ K_s &= K_n \frac{1-zv}{1+v} \end{aligned} \right. \qquad (1) $$

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.

$$ \left\{ \begin{aligned} f_1(F_n, F_s) &= F_s - \tan(\Phi_i)F_n - S_{int}C_o \\ f_2(F_n, F_s) &= S_{int}T - F_n \end{aligned} \right. \qquad (2) $$

In (2) $\Phi_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 $\zeta$ need to be identified from global parameters such as compressive and tensile strengths $\sigma_c$ and $\sigma_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 $\zeta > 0$, which allows reproducing the macroscopic softening observed in experiments. If $\zeta$ 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 $\Phi_c$ is possible between two considered elements.

$$ F_n = \frac{K_n}{\zeta} (D - D_{max}), \quad \text{with} \quad D_{max} = D_{init} + (1 + \zeta) \frac{S_{int}T}{K_n} $$

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).

<figure><img src="image_3.png" /><figcaption>Figure 2. Interaction laws for concrete.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3. View of a DE reinforcement bar.</figcaption></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4. Interaction laws for the steel reinforcement in (a) traction, and (b) bending.</figcaption></figure>

## 4. 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 $\sigma_y$, the DE bars can undergo plastic deformations with a classical strain hardening. The link is broken in traction if the ultimate stress $\sigma_u$ or the maximum prescribed elongation $D_{max}$ are reached ($D_{max} = D_{init} + \varepsilon_{max}D_{init}$ with $\varepsilon_{max}$ the rupture strain).

## 5. Steel-concrete link

The concrete-steel interaction laws adopted here are derived from ideas proposed in the finite element context (Dominguez, 2005). Thus, the tangential law (Figure 5, left) reproduces at first a perfect stick of the two materials with an elastic stiffness $K_y$, 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.

## 6. 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
<figure><img src="image_7.png" /><figcaption>Figure 5. Interaction laws for the steel-concrete link in tangential (left) and normal (right) directions.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Bridging domain.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7. Original tetrahedral mesh and the discrete element mesh produced in SALOME (2011).</figcaption></figure>

molecular dynamics model and adapted to the DE context in Frangin et al. (2006) and Rousseau, Marin, Daudeville, and Potapov (2010). The idea is to write the global Hamiltonian as a linear combination of the discrete and the finite element ones:

$$H = \eta H_{FE} + (1 - \eta) H_{DE}.$$
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.

## 7. 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.

<figure><img src="image_10.png" /><figcaption>Figure 8. Discrete element mesh for a RC structure produced in SALOME (2011).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9. Typical size distribution for a discrete element mesh with DE radii on abscissa axis.</figcaption></figure>
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.

## 8. 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><img src="image_12.png" /><figcaption>Figure 10. Example of subdomains definition for the problem of impact of finite element projectile on a discrete element slab.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11. Subdomains in a combined finite element-discrete element model.</figcaption></figure>

<figure><img src="image_14.png" /></figure>

White and blue discrete elements belong to sub-domains SD1 and SD2 respectively.

The influence domain of the discrete element, represented by large blue circle and belonging to sub-domain SD2, includes red remote discrete elements belonging to sub-domain SD1.

Green line represents the fictitious boundary between sub-domains.

<caption>Figure 12. Definition of a “crossover” neighbourhood for a near-to-border discrete element (large blue circle in centre).</caption>
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).

## 9. 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><img src="image_15.png" /><figcaption>Figure 13. Damage maps during the impact of a rigid projectile on a RC beam.</figcaption></figure>
EUROPLEXUS

Meppen test simulation

Reinforced Concrete Slab: Discrete Elements

Missile: Finite Elements - Johnson-Cook material law

Impact velocity: 247.7 m/s

MPI calculation - 6 CPU

<figure><img src="image_16.png" /></figure>

Time: 5.0 ms

<figure><img src="image_17.png" /><figcaption>Figure 14. Impact of a tube-like deformable projectile on a RC slab.</figcaption></figure>

<table><caption>Table 1. Extensibility measures for the first example.</caption><thead><tr><th></th><th>Time elapsed<br>for 10,000 cycles</th><th>Speed-up</th><th>Efficiency</th></tr></thead><tbody><tr><td>1 processor</td><td>4294 s</td><td></td><td></td></tr><tr><td>2 processors</td><td>2537 s</td><td>1.69</td><td>0.85</td></tr><tr><td>4 processors</td><td>1140 s</td><td>3.77</td><td>0.94</td></tr><tr><td>8 processors</td><td>728 s</td><td>5.90</td><td>0.74</td></tr><tr><td>16 processors</td><td>604 s</td><td>7.13</td><td>0.45</td></tr></tbody></table>

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 sub-domain 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 sub-domains become too costly, limiting extensibility to 16 processors. It is necessary to make a new partitioning of the domain to maintain resolution efficiency.

## 10. 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 pro-
cessors. 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 dis-
crete elements is used, since the time scale associated to finite elements is classically
greater than that associated to discrete elements.

References

Belytschko, T., & Lin, J.I. (1987). A three-dimensional impact-penetration algorithm with erosion. *International Journal of Impact Engineering*, 5(1-4), 111-127.

Casadei, F. (2003). *A general impact-contact algorithm based on hierarchic pinballs for the Euro-plexus software system*. Ispra, Italy. European Commission, Institute for the Protection and Security of the Citizen, Joint Research Centre.

Dominguez, N. (2005). *Etude de la liaison Acier-Béton: De la modélisation du phénomène à la formulation d'un élément enrichi de béton armé* (Thèse de doctorat). ENS Cachan.

EUROPLEXUS: A computer program for finite element simulation of fluid-structure systems under transient dynamic loading – user’s manual. (2011). Retrieved from http://europlexus.jrc.ec.europa.eu/.

Frangin, E., Marin, P., & Daudeville, L. (2006). On the use of combined finite/discrete element method for impacted concrete structures. *Journal of Physics IV, 134*, 461–466.

Hentz, S., Daudeville, L., & Donze, F. (2004). Identification and validation of a discrete element model for concrete. *Journal of Engineering Mechanics, 130*(6), 709–719.

Jerier, J.-F., Richefeu, V., Imbault, D., & Donzé, F. (2010). Packing spherical discrete elements for large scale simulations. *Computer Methods in Applied Mechanics and Engineering, 199* (25–28), 1668–1676.

Jonas, W., Meschkat, R., Riech, H., & Rüdiger, E. (1979). Experimental investigations to determine the kinetic ultimate bearing capacity of reinforced concrete slabs subject to deformable missiles. In: *Transactions of the 5th SMiRT Conference*, Berlin, paper no. 8/5.

REPDYN Project. (2011). Retrieved from http://www.repdyn.fr.

Rousseau, J., Frangin, E., Marin, P., & Daudeville, L. (2008). Damage prediction in the vicinity of an impact on a concrete structure: a combined FEM/DEM approach. *Computers and Concrete*, 5(4), 343–358.

Rousseau, J. (2009a). *Modélisation numérique du comportement dynamique des structures sous impact sévère avec un couplage éléments discrets/éléments finis* (Thèse de doctorat). Université Joseph Fournier, Grenoble.

Rousseau, J., Frangin, E., Marin, P., & Daudeville, L. (2009). Multidomain finite and discrete elements method for impact analysis of a concrete structure. *Engineering Structures*, 31(11), 2735–2743.

Rousseau, J., Marin, P., Daudeville, L., & Potapov, S. (2010). A discrete element/shell finite element coupling for simulating impacts on reinforced concrete structures. *European Journal of Computational Mechanics*, 19(1–3), 153–164.

SALOME: the open source integration platform for numerical simulation. (2011). Retrieved from http://www.salome-platform.org/.

Xiao, S., & Belytschko, T. (2004). A bridging domain method for coupling continua with molecular dynamics. *Computer Methods in Applied Mechanics and Engineering*, 193, 1645–1669.

YADE. (2011). Retrieved from https://yade-dem.org.