# Three-dimensional coupled discrete element - Finite element model: Parameter identification and coupling with shells

Jessica Rousseau, Emmanuel Frangin, Philippe Maurice Marin, Laurent Daudeville, Serguei Potapov

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

## To cite this version:

Jessica Rousseau, Emmanuel Frangin, Philippe Maurice Marin, Laurent Daudeville, Serguei Potapov. Three-dimensional coupled discrete element - Finite element model: Parameter identification and coupling with shells. CST 20085 - 9th International Conference on Computational Structures Technology 2008, Sep 2008, Athenes, Greece. 10.4203/ccp.88.16. hal-02004492

HAL Id: hal-02004492

https://hal.science/hal-02004492v1

Submitted on 30 Jun 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Three-Dimensional Coupled Discrete Element - Finite Element Model: Parameter Identification and Coupling with Shells

J. Rousseau<fnref n="1" />, E. Frangin<fnref n="1" />, P. Marin<fnref n="1" />, L. Daudeville<fnref n="1" /> and S. Potapov<fnref n="2" />

doi:10.4203/ccp.88.16

## Abstract

This paper focuses on concrete structures submitted to impact loading and is aimed at predicting local damage in the vicinity of an impact zone as well as the global response of the structure. The Discrete Element Method (DEM) seems particularly well suited in this context for modeling fractures. An identification process of DEM material parameters from macroscopic data (Young's modulus, compressive and tensile strength, fracture energy, etc.) will first be presented for the purpose of enhancing reproducibility and reliability of the simulation results with DE samples of various sizes.

The modelling of a large structure by means of DEM may lead to prohibitive computation times. A refined discretization becomes required in the vicinity of the impact, while the structure may be modeled using a coarse FE mesh further from the impact area, where the material behaves elastically. A coupled discrete-finite element approach is thus proposed to reduce computation times.

**Keywords:** Discrete Element, Reinforced concrete, Shell, Transient dynamics.

## 1 Introduction

The design of some civil engineering structures must take into account the risk of severe loadings due to natural or manmade hazards, such as rockfalls, aircraft or missile impacts. On a sensitive concrete structure, such loadings may have disastrous consequences, and an efficient and accurate prediction of damage proves to be of prime importance. Generally speaking, these severe loadings lead to fractures and fragmentation localized in one part of the concrete structure. The Discrete Element Method (DEM) (Cundall and Strack [1]) is appropriate for modeling such discontinuities. The model developed uses a disordered assembly of spherical elements of different sizes and masses in order to reproduce isotropic and

<footnote n="1">¹Laboratory Soils, Solids, Structures - Risks (3S-R), University of Grenoble, France</footnote>

<footnote n="2">²LaMSID UMR EDF/CNRS, Clamart, France</footnote>
homogeneous behavior at the macroscopic scale. This method does not rely on any assumption regarding where and how a crack or several cracks occur and propagate since the medium is naturally discontinuous and very well adapted to dynamic problems. In order to offer a predictive model, the DE model must properly depict concrete behavior.

The identification of the DEM material parameters from macroscopic data (Young's modulus, compressive and tensile strength, fracture energy...) will be presented in the first part.

Nevertheless, computation time increases with the number of DE, and the analysis of large structures with DEM proves difficult. Use of the FE method further from the impacted area represents one way to minimize this constraint since in most cases, severe degradation phenomena are localized right in the vicinity of the impact. In addition, the FE method is widely used, and efficient mesh generation software can dramatically reduce modeling duration, with the potential for faster calculations than when applying a full DE approach thanks to the facility of handling various discretization sizes. These facts naturally lead to proposing a coupled FE/DE approach. Thus a coupled discrete-finite element approach is proposed: the impact zone is modelled by means of DE and elastic FE are used for the rest of the structure. The proposed approach will be described. An application of this coupling method will be shown by means of the response of a concrete slab impacted by a rocky block. Finally, a coupling method between DE and Shells will be described and some static tests will be shown.

# 2 Identification of the Discrete Element Model parameters

## 2.1 Discrete Element Model

Discrete methods can be split into two groups: lattice models and discrete element models. Lattice models use simple deformable elements like beams or bars organized into a structural network; these models are widely used to represent concrete at the microscopic or mesolevel scale (Leite et al. [2]; Prado and Van Mier [3], Cusatis and Cedolin [4]).

The second group of methods, i.e. using discrete elements, is based on modeling the continuum by means of rigid particles. The interaction laws between discrete elements serve to determine the macroscopic constitutive behavior. During early developments for non-cohesive materials like sands (Cundall and Strack [1]), particle interactions were described by friction laws, but the interaction laws for cohesive materials have since become refined (Hentz et al. [5] [6]).

To guarantee reasonable computation times, a model based on Distinct Elements (Cundall and Strack [1]) with rigid spheres has been chosen herein. Two types of
interactions are defined. The initial interactions between two elements are generally link interaction (the two elements are not necessarily in contact). Additional interactions of contact type can be added during the simulation.

The used Discrete Element method is fully described in Hentz et al [5] [6]. Interactions between two spheres are defined with normal $K_n$ and tangential $K_s$ stiffness characterizing the elastic behaviour of concrete. “Micro-macro” relations (1) are supposed in order to compute the stiffness $K_n$ and $K_s$ between elements a and b from Young’s Modulus and Poisson’s ratio. Parameters $\alpha, \beta, \gamma$ of these relations were identified for a given average number of interactions per discrete element. Stiffnesses $K_n$ and $K_s$ are used to compute the normal and tangential forces, $F_n$ and $F_s$.

$$ \left\{ \begin{array}{l} E = \frac{D_{init}^{a,b}}{S_{int}} K_n \frac{\beta + \gamma \frac{K_s}{K_n}}{\alpha + \frac{K_s}{K_n}} \quad \text{or} \quad \left\{ \begin{array}{l} K_n = E \frac{S_{int}}{D_{init}^{a,b}} \frac{1+\alpha}{\beta(1+\nu) + \gamma(1-\alpha\nu)} \\ K_s = K_n \frac{1-\alpha\nu}{1+\nu} \end{array} \right. \\ \nu = \frac{1-\frac{K_s}{K_n}}{\alpha + \frac{K_s}{K_n}} \end{array} \right. \quad (1) $$

The aim of this study is to model very large structures or parts of structures. The interaction laws are therefore stated in order to best represent concrete behavior at the macroscopic scale. A modified Mohr-Coulomb criterion associated with softening is used to model the non-linear behaviour of the material (see Figure 1. and 2.). Local parameters, as local tensile strength $T$, cohesion $C_o$ and softening factor $\zeta$ need to be identified from global parameters as compressive and tensile strength $\sigma_c$ and $\sigma_t$ and fracture energy $G_f$. A classical friction constitutive behaviour is used between elements in contact.

<figure><img src="image_3.png" /><figcaption>Figure 1: Softening law</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Interaction laws for DE</figcaption></figure>

## 2.2 Parameters identification of linear elastic constitutive behaviour

Simulations of a compressive test with a DE sample are realized with different values of Ks/Kn. Numerical values of E and v are computed and values of α, β,γ are then identified with (1) in order to fit numerical points of E and v versus Ks/Kn (Figure 3.).

<figure><img src="image_5.png" /><figcaption>Figure 3: Comparison of Young's modulus and Poisson's ratio computed and derived by approximation</figcaption></figure>

Tests have been carried out with 15 different samples (10 of them have identical sizes but different disorder and the 5 others have different geometries). The main assumption of this approach is that samples are representative of modelled concrete material in further simulations of impact on concrete structures. Samples have to be isotropic. With the identified values of α, β,γ parameters, average uncertainty is about 2% (maximal error is 7%) for Young's modulus and about 3% (maximal error is 6%) for the Poisson's ratio.
## 2.3 Parameters identification of non-linear constitutive behaviour.

In this part, the point is to identify local parameters to model macroscopic values as compressive strength $σ_c$, tensile strength $σ_t$ and fracture energy $G_f$. Only local tensile strength $T$, cohesion $C_o$ and the softening factor $ζ$ will vary in this study, friction parameters will be fixed values.

The aim is to realize traction and compression simulations in order to simulate tests performed in 3S-R laboratory (Gabet *et al* [7]). We identified parameters according to Hentz's procedure. In Figure 4., we can observe the result of the traction simulation and the compressive tests on Figure 5.

<figure><img src="image_6.png" /><figcaption>Figure 4: Uniaxial Traction : Stress-Strain curve</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Uniaxial Compression : Stress-Strain curve</figcaption></figure>

# 3 Coupling FE/DE method

## 3.1 Description of the method

This third part concerns the coupling between continuum and discrete domains. Discrete Element Method is convenient for modeling discontinuities and fracture but it may lead to prohibitive durations of calculation with a large structure. During an impact, though the impact zone needs a refined mesh, the main part of the structure far away from the impact zone can be modeled by means of coarse finite elements. We propose to couple discrete elements in the vicinity of the impact with elastic finite elements in the rest of the structure. It is chosen to use a bridging domain in order to couple DE and FE degrees of freedom.

Many methods have been already developed for multi-scaled problems. Xiao and Belytchko [8] proposed to use a bridging domain where the Hamiltonian is a linear
combination of discrete and finite ones. Equation (2) shows DE and FE Hamiltonian are weighted by a parameter α equals to 1 in the FE domain and equals to 0 in the DE domain.

$$H = \alpha H_{FE} + (1 - \alpha) H_{DE} \qquad (2)$$

To guarantee kinematic continuity, the degrees of freedom of both domains within the interface zone must be linked. Several approaches could be considered. Xiao and Belytschko [8] proposed directly linking the Discrete Element degrees of freedom with Finite Element degrees of freedom using Lagrange multipliers. Eq. (3) presents the cinematic conditions on bridging domain r, with d DE displacements and u FE displacements.DE rotations ωr are linked with the antisymmetric part of the FE displacement gradient. Details of the method are explained in Frangin et al [9].

$$\vec{d}_r = \bar{k} \cdot \vec{u}_r \quad \text{and} \quad \bar{\omega}_r = \bar{h} \cdot \vec{u}_r \qquad (3)$$

The difference in discretization size between FE and DE may induce wave reflection at the interface, an effect that can be mitigated in different ways, e.g. damping, yet the choice of damping coefficient is not straightforward. We have proposed a method to attenuate the reflection by introducing a reduction parameter for Lagrange multiplier influence. This method leads to a temporal relaxation of the cinematic constraints and is equivalent to a penalty method with an automatic process for optimizing the penalty parameter. All of these methods are discussed in Frangin et al [9]. In the next section, the application presented introduces just such a relaxation parameter, which serves to divide the influence of Lagrange multipliers in reducing wave reflection.

## 3.2 Impact on a concrete structure

To quantify the efficiency of this method, we propose to compare damage and displacement response of a concrete slab under a rocky impact (Figure 6.) between combined FE/DE method and the full DE method (Figure 7.).

This application consists of a rock impact on a concrete slab. A cubic rock block with a 30-cm side length impacts a concrete slab 2.5 m long, 2.0 m wide and 0.28 m thick. The velocity of the impacting object is set at 40 m/s. The two opposite sides are locked in the direction perpendicular to the medium slab plane. The slab is first divided into two parts: the center modeled by DE, and the sides modeled by FE with three bridging domain layers (Figure 6.). In a second model, only the DE are used.
<figure><img src="image_8.png" /><figcaption>Figure 6: Combined FE/DE model, and a side view of the DE part</figcaption></figure>

Displacements and damage of the slab are very close at each part of the slab. The combined method is 10 times faster than the fully DE method without a multi-time step algorithm.

<figure><img src="image_9.png" /><figcaption>Figure 7: Comparison of displacement between DE (red) method and coupled method (blue)</figcaption></figure>

## 3.3 Coupling DE with Shells

### 3.3.1 Method

Finite Elements of shell types are now widely used to model very large structures such as sensitive concrete buildings. Coupling shells with DE becomes unavoidable. This part described the approach for coupling shells (Batoz and Dhatt [10]) using the
coupling method defined in section 3.1. The idea is to use the previously implemented algorithm of coupling DE with 3D FE. The shell element used is a 4 nodes element based on Reissner/Mindlin model (Reissner [11]).

So, the first step is to build a fictive 3D element from the real shell mesh. For each shell node, a normal vector $n$ is computed as the average of all the normal vectors of the connected elements. The average method guarantees to cover the entire domain without hole or superposition. The normal vector $n$ and the thickness $h$ are later used to define fictive nodes, Equation (4). The fictive nodes are on both sides of the real node.

$$ \bar{x}_{fictive} = \bar{x}_{real} \pm \frac{h}{2}\bar{n} \qquad (4) $$

Discrete Elements are then coupled with the fictive nodes using the previously described method (section 3.1), identically with 3D Finite Elements.

Finally, the transition between fictive nodes displacements $u_{fictive}$ and real nodes displacements $u_{real}$ and rotations $\theta$ are computed using Equation (5), supposing the thickness segment keeps perpendicular and rigid (Reissner [11]).

$$ \bar{u}_{fictive} = \bar{u}_{real} + z\bar{\theta} \wedge \bar{n} \quad \text{with} \quad z = \pm \frac{h}{2} \qquad (5) $$

### 3.3.2 Preliminary tests

To validate this new coupling, panels of elastic test were computed. First, a tensile test was computed to verify if the kinematics continuity is correctly established, and to verify the elastic behaviour of the model. A concrete bar (0.85m x 0.25m x 0.25m) is made in one half with FE Shells and in the other half of DE with 3 layers of bridging domain. One end is blocked and a displacement is applied on the other end of the sample. Figure 8 shows the displacement of a DE and the nearest FE.

<figure><img src="image_10.png" /><figcaption>Figure 8: Displacement of a DE and an adjacent FE in the bridging domain</figcaption></figure>
The Young Modulus measured is about 31Gpa, that is to say an error of about 3%.
This test seems to show that the elastic behaviour in tension is correctly reproduced.

<figure><img src="image_11.png" /><figcaption>Figure 9: Flexion test</figcaption></figure>

The second test is a flexion test of a concrete beam (4m x 0.2m x 0.2m). A force of 9N is applied at one extremity when the other is blocked. The final shape of the beam (Figure 9.) seems to correspond to the static determination of the beam. Globally, the elastic responses of these tests are quite good and the coupling method between Shells and DE seem to be acceptable. Obviously, these tests have to be completed with dynamic tests before its use in a very large structure modelling.

**3.4 Conclusion**

The identification of parameters for concrete behavior gives very accurate results.
First, the elastic behavior seems to be reliable. And secondly, uniaxial tensile and
compressive behaviors can be easily computed.

This identification process combined with an efficient coupled DE/FE (3D elements
or Shells) method give solution to simulate large RC structure under impact.

**References**

[1] Cundall, P. A. and Strack, O.D.L. "A discrete numerical model for granular assemblies", Géotechnique, 29(1) 47-65, 1979.

[2] Leite, J.P.B., Slowik, V. and Mihasshi, H., "Computer simulation of fracture processes of concrete using mesolevel models of lattice structures", Cem. Concr. Res., 34, 1025-1033, 2004.
[3] Prado, E.P. and Van Mier, J.G.M. " Effect of particle structure on mode I fracture process in concrete", Engng. Fract. Mech., 70, 1793-1807, 2003.

[4] Cusatis, G. and Cedolin, L. "Tow-scale study of concrete fracturing behavior", Engng. Fract. Mech., 74, 3-17, 2007.

[5] Hentz S., Daudeville L., Donzé F-V. "Identification and validation of a discrete element model for concrete", ASCE Journal of Engineering Mechanics. 130(6): 709-719, 2004.

[6] Hentz S., Daudeville L. and Donzé F. "Discrete Element Modelling of Concrete Submitted To Dynamic Loading at High Strain Rate", Comput. Struct., 82, 2509-2524, 2004.

[7] Gabet, T. Malécot, Y. and Daudeville, L. "Triaxial behaviour of concrete under high stresses: Influence of the loading path on compaction and limit states", Cem. Concr. Res., 38(3), 403-412, 2008.

[8] Xiao, S.P. and Belytschko, T. "A bridging domain method for coupling continua with molecular dynamics", Comput. Methods Appl. Mech. Engrg. 193 1645-1669, 2004.

[9] E. Frangin, P. Marin, L. Daudeville "On the use of combined finite/discrete element method for impacted concrete structures", Journal of Physics IV, 2006.

[10] J-L. Batoz and G. Dhatt, "Modélisation des structures par éléments finis", Hermes editions, volume 3: Coques, 1992.

[11] E. Reissner, "Linear and nonlinear theory of shells", in Thin Shell Structures, (Fung and Sechler Eds.), Prentice Hall, p. 29-44, 1974.