A coupled implicit material point-finite element method for modeling reinforced materials

Ahmad Chihadeh (Institute for Structural Analysis, Technische Universität Dresden, Dresden, Germany)
Michael Kaliske (Institute for Structural Analysis, Technische Universität Dresden, Dresden, Germany)

Engineering Computations

ISSN: 0264-4401

Article publication date: 14 April 2022

Issue publication date: 5 July 2022

383

Abstract

Purpose

This paper aims to introduce a method to couple truss finite elements to the material point method (MPM). It presents modeling reinforced material using MPM and describes how to consider the bond behavior between the reinforcement and the continuum.

Design/methodology/approach

The embedded approach is used for coupling reinforcement bars with continuum elements. This description is achieved by coupling continuum elements in the background mesh to the reinforcement bars, which are described using truss- finite elements. The coupling is implemented between the truss elements and the continuum elements in the background mesh through bond elements that allow for freely distributed truss elements independent of the continuum element discretization. The bond elements allow for modeling the bond behavior between the reinforcement and the continuum.

Findings

The paper introduces a novel method to include the reinforcement bars in the MPM applications. The reinforcement bars can be modeled without any constraints with a bond-slip constitutive model being considered.

Originality/value

As modeling of reinforced materials is required in a wide range of applications, a method to include the reinforcement into the MPM framework is required. The proposed approach allows for modeling reinforced material within MPM applications.

Keywords

Citation

Chihadeh, A. and Kaliske, M. (2022), "A coupled implicit material point-finite element method for modeling reinforced materials", Engineering Computations, Vol. 39 No. 7, pp. 2553-2580. https://doi.org/10.1108/EC-10-2021-0623

Publisher

:

Emerald Publishing Limited

Copyright © 2022, Emerald Publishing Limited


Highlights

  1. Modeling of reinforced materials by the generalized interpolation material point (GIMP) method is introduced.

  2. The coupling of the GIMP method with the truss finite elements is shown.

  3. Bond elements are used to link the truss finite elements to the continuum elements.

1. Introduction

Large deformation problems in solid mechanics are often a challenge for engineers. Such problems, e.g. plastic forming, soil collapse and high- speed load applications like impact scenarios, require numerical methods, which are capable of modeling the structural response with satisfying accuracy and acceptable computational cost.

The finite element method (FEM), which is the most widely used numerical method in engineering applications for solid mechanics, fails when the deformations of the mesh are extremely high and the elements are strongly distorted. For such cases, the MPM is an alternative numerical method.

The MPM has been first introduced by Sulsky et al. (1994, 1995) under the name particle method, which was later named as the MPM in Sulsky and Schreyer (1996). The idea of the MPM is inspired by the fluid implicit particle (FLIP) method (Brackbill and Ruppel, 1986) for fluid mechanics with two main differences; first, the constitutive equations are solved at the material points instead of the cell center; second, the MPM formulation is consistent with the FEM, which yields an efficient combination of both methods (Lian and Zhang, 2018; Lian et al., 2011a; Zhang et al., 2016). Subsequently, the Generalized Interpolation Material Point (GIMP) method was proposed by Bardenhagen and Kober (2004) to overcome the problem of cell crossing noise which occurs when the material points cross elements in the background mesh.

Formerly, explicit time integration has been employed for most of the MPM developments, while the implicit MPM has been recently considered (Charlton et al., 2017; Coombs et al., 2020; Wang et al., 2016). The advantages of the implicit MPM over the explicit MPM are the larger time steps, better stability and higher guarantee of accurate results (see Guilkey and Weiss, 2003).

In this paper, the implicit GIMP method for large deformations is addressed for quasi static applications.

As the MPM is used to model large deformation problems, such structures may include reinforced materials. This combination is applied in different fields as reinforced soil fill, reinforced rubber or reinforced concrete. Reinforced concrete, for instance, is used to construct civilian buildings, nuclear structures, dams, etc, which is important to investigate its responses to blast and impact loading. One example is ongoing research on reinforced concrete under impact, see Figure 1 (Hering et al., 2020), for which a penetration problem is addressed.

Including the reinforcement into the continuum can be achieved by many approaches. Discrete approach: The domain is discretized in which the reinforcement bars, modeled by truss elements, are connected to the continuum element nodes. Comparing to other approaches, this is the most accurate one. Bond slip can be introduced (see Häussler-Combe, 2015; Ngo and Scordelis, 1967; Schäfer, 1975). This approach has the main disadvantage, that the discretization of the continuum is dependent on the reinforcement layout. Smeared approach: This formulation is suitable to model structures with uniformly distributed reinforcement bars, e.g. reinforced plates. In the smeared approach, the reinforcement bars are smeared as a layer of equivalent thickness placed at the position of the bars (Häussler-Combe, 2015). The effect of dowel action can be included (see Zhang et al., 1994). Embedded approach: The main advantage of this approach is that the discretization of the continuum is independent of the reinforcement layout. The reinforcement bars are modeled explicitly by truss elements and embedded in the continuum elements (Elwi and Hrudey, 1989; Hartl, 2002). Truss elements and continuum elements are coupled by means of a special type of elements. This element is a so-called bond element, which is very stiff in the lateral direction and has a bond law in the longitudinal direction (see Häussler-Combe et al., 2019, 2020).

In this work, coupling of the GIMP method with truss- finite elements is introduced using the embedded approach. The coupling is achieved by means of bond elements, which allow for coupling freely distributed truss elements with the structured background mesh. This approach is well suited for the MPM as it gives flexibility for discretizing the truss elements independently of the background mesh. Truss elements can move over the computational mesh, and it allows modeling the bond behavior between the reinforcement and the continuum.

Coupling MPM with FEM in general has been investigated by many researches (see Chen et al., 2015; Lian et al., 2011a; Zhang et al., 2006). It is of interest because of the attempt to take advantage of both methods. MPM is more robust than FEM for applications with large deformations, as the problem of element distortion is avoided in the regions with severe deformations. On the other side in the regions of small deformations, FEM is more accurate, as the accuracy of Gauss quadrature in FEM is higher than material point quadrature and is computationally more efficient, as no mapping of data is required and it needs less computational storage. Coupling the MPM with other numerical methods has been introduced as well. For instance, coupling of MPM with discrete element method (DEM) has been introduced in Liang and Zhao (2020).

Introducing truss elements into the MPM is a special type of coupling MPM–FEM and has been proposed in Lian et al. (2011b). The proposed approach in Lian et al. (2011b) differs from the one in this work at hand. It is implemented in explicit MPM, bond slip between the truss elements and the continuum elements is not included and all data from both material points and truss element nodes are mapped to the grid nodes. However, in the presented work, bond slip is included and truss element data are not mapped to the grid nodes, as the truss elements are included explicitly in the system and have their own nodes in the global stiffness matrix.

The paper at hand is organized as follows: In Section 2, the basics of the MPM and the formulation of the non-linear implicit MPM are presented. The GIMP method is presented in Section 3. The coupling of the GIMP method with truss- finite elements is introduced in Section 4. Numerical examples are shown and commented in Section 5. Finally, in Section 6, conclusive remarks are given.

The authors of this paper restrict themselves to two-dimentional problems for simplicity and clarity, however, the extension to three-dimensional is straight forward.

2. Implicit MPM

2.1 Basics

The main feature of the MPM is moving of the material points, carrying all material properties, over a continuously reset mesh. The computational cycle can be divided into three phases. In the first phase, all state variables are mapped from the material points to the nodes, known as the Lagrangian phase. In the second phase, the system is solved for the unknown displacements, usually in an iterative manner, e.g. by the Newton–Raphson method. In the third phase, the state variables are mapped back to the material points and the mesh is reset to its undeformed configuration, this phase is known as convective phase. Figure 2 illustrates the three phases.

The mapping of the state variables is achieved using shape functions. For instance, the nodal force vector fn of an element in the mesh is given by

(1)fn=p=1npN(ξp)Tfp,
where N is the matrix of the shape functions, ξp is the local coordinate of material point p, fp is the material point force vector. On the other side, the displacement vector up of a material point is given by
(2)up=N(ξp)un,
where un is the nodal displacement vector of the corresponding element in the mesh.

2.2 Formulation of non-linear implicit MPM

In non-linear finite element analysis, the reference configuration is taken either the initial configuration or the current configuration. In the MPM, due to mesh resetting and also changing the elements each time step, considering the initial configuration as the reference configuration is not a good choice. In Coombs et al. (2020), taking the reference configuration at the previously converged state (the state at the beginning of a time step) is proposed as it is computationally more efficient than the current configuration formulation. Hence, in this paper, the reference configuration is taken at the previously converged state. The discretized weak form of the equilibrium is given by

(3)VtBT(Δu)σ̃dVt=VtNTbdVt+AtNTtdAt,
with the matrix of shape function derivatives, B, function of the displacement, Δu, starting from the beginning of the time step, t, σ̃ is the second Piola Kirchhoff stress tensor pushed forward to the previously converged state, body forces b acting over the volume at the previously converged state Vt, and ∂At is the part of the domain surface, where tractions t are applied (Bathe, 2006; De Borst et al., 2012).

Eq. (3) is non-linear in terms of the nodal displacements. Using the Newton–Raphson method, and after linearization, the global vector of internal forces and the global stiffness matrix are given by (Ghorbani et al., 2009)

(4)f=VtBLT(Δu)σ̃dVt,
(5)K=VtBLT(Δu)CTBL(Δu)+BNLTSBNLdVt,
respectively. Those quantities are built by assembly of all material points contributions in a similar way to the Gauss points in FEM. The vector of internal forces and the stiffness matrix of a material point are
(6)fp=BLT(Δu)σ̃Vp,
(7)Kp=BLT(Δu)CTBL(Δu)+BNLTSBNLVp,
where BL and BNL are matrices of the shape function derivatives, S is the second Piola Kirchhoff stress tensor pushed forward to the previously converged state, and CT is the tangential material matrix, i.e. tensor. In the presented work, the simple St. Venant–Kirchhoff material model (Ogden, 1997) is used when elastic material is described. Thus, in this case, CT is a constant material tensor C0 mapped to the previously converged state and Vp is the material point volume at the previously converged state.

The abovementioned components are given as

(8)BL(Δu)=BL0+BL1(Δu),
where
(9)BL0=N1,x0N2,x0N3,x0N4,x00N1,y0N2,y0N3,y0N4,yN1,yN1,xN2,yN2,xN3,yN3,xN4,yN4,x,
(10)BL1(Δu)=lxxN1,xlyxN1,xlxxN2,xlxyN1,ylyyN1,ylxyN2,ylxxN1,y+lxyN1,xlyxN1,y+lyyN1,xlxxN2,y+lxyN2,xlyxN2,xlyxN4,xlyyN2,ylyyN4,ylyxN2,y+lyyN2,xlyxN4,y+lyyN4,x,
with
(11)lij=k=14Nk,jΔuik,
where Δu equals to 0 for the first iteration in each time step. The tangential material matrix mapped to the previously converged state is given by
(12)CT=1det(Ft)TTC0T,
with the transformation matrix (Bathe, 2006; De Borst et al., 2012)
(13)T=Ft(1,1)2Ft(2,1)2Ft(1,1)Ft(2,1)Ft(1,2)2Ft(2,2)2Ft(1,2)Ft(2,2)2Ft(1,1)Ft(1,2)2Ft(2,1)Ft(2,2)Ft(1,1)Ft(2,2)+Ft(2,1)Ft(1,2).
Ft is the deformation gradient at the previously converged state given by
(14)Ft=utx+1,
where ut is the material point displacement vector at the previously converged state and x is the nodal coordinate vector of the element in the background mesh. The second term in Eq. (7) is required in order to include geometrical non-linearity, with
(15)BNL=N1,x0N2,x0N3,x0N4,x0N1,y0N2,y0N3,y0N4,y00N1,x0N2,x0N3,x0N4,x0N1,y0N2,y0N3,y0N4,y,
and
(16)S=S¯xxS¯xy00S¯yxS¯yy0000S¯xxS¯xy00S¯yxS¯yy,
with
(17)S¯=1det(Ft)FtS¯0FtT,
(18)S¯0=s¯0xxs¯0xys¯0xys¯0yy,
(19)s0=C0ϵ.
ϵ is the Green–Lagrange strain tensor. σ̃ in Eq. (6) is then
(20)σ̃T=S¯xxS¯yyS¯xy.
Finally, the volume assigned to a material point is
(21)Vp=det(Ft)Vp0,
with Vp0 as the initial volume.

For the formulation shown above, the local coordinates of the material points in the background mesh elements are required. Those coordinates are changing each time step and evaluated using the Newton–Raphson method. Given material point global coordinates, xp, and element coordinates, x, material point local coordinates are calculated as shown in Figure 3.

3. Implicit GIMP method

As mentioned previously, GIMP method is used to reduce the cell crossing noise. This is achieved by replacing the standard shape function N used in the MPM with a weighting function Svp. Each material point is assigned a domain. The weighting function, as proposed in Bardenhagen and Kober (2004), is given by

(22)Svp=1VpΩpΩχp(ξ)Nv(ξ)dξ,
where Ωp is the influence domain of the material point, Ω is the problem domain, and χp is a characteristic function in the form
(23)χp(ξ)=1,ifξΩp0,otherwise.
Note that the standard shape function can be recovered from Eq. (22) if χp(ξ) is as follows
(24)χp(ξ)=Vp,ifξ=ξp0,otherwise.
The gradient of the weighting function ∇Svp is given by
(25)Svp=1VpΩpΩχp(ξ)Nv(ξ)dξ,

The weighting function gradient, ∇Svp, of node n is equal to zero at node position, thus, the contribution to the internal force vector of node n is 0 from a material point, when it is located exactly at node n, see Equations (4) and (6). In order to calculate the element stiffness matrix in the same manner as in the FEM, Charlton et al. (2017) introduced a modification to the weighting function proposed in Bardenhagen and Kober (2004), Equations (22) and (25). The proposed modification is based on considering only the part of the material point's domain overlapping with the element. The main idea is that the summation of the function values from all elements sharing the same node is equal to the original function value proposed by Bardenhagen and Kober (2004). The weighting functions of the one-dimensional element with Nodes 1 and 2 as given by Charlton et al. (2017) are

(26)Svp1=14lp2ξ2ξ222ξ1+ξ12,Svp2=14lp2ξ2+ξ222ξ1ξ12.

In the same manner, the gradient of the weighting functions are constructed and given by

(27)Svp1=ξ1ξ22lpandSvp2=ξ2ξ12lp,
with
(28)ξ1=1,ifξplp2<1ξplp2,ifξplp2>1,ξ2=1,ifξp+lp2>1ξp+lp2,ifξp+lp2<1,
where lp is the material point influence domain length. For 2D and 3D problems, the weighting function is simply the multiplication of the 1D weighting function in each direction.

The length of the material point influence domain, lp, is updated at each time step. Different methods for updating lp are proposed in Coombs et al. (2020). In this work, for 2D problems, lp is updated by

(29)lpi=lpi0Uii(no implied sum oni),
where U is the stretch tensor of the deformation gradient, lpi0 is the domain length in the initial configuration and i indicates the direction (x, y). The algorithm of the implicit GIMP method is shown in Appendix.

4. Coupling the GIMP method with truss finite elements

In many applications, reinforcements are used to strengthen materials, e.g. reinforced concrete and reinforced soil. In the FEM, truss elements are widely used to model reinforcements. However, beam elements and continuum elements are also employed, depending on the purpose of the simulation and the required accuracy (Fleischhauer et al., 2015; Fuchs et al., 2020; Mahmoud, 2016).

In this work, truss elements are implemented and coupled with the GIMP method by means of bond elements to link or to bond the truss elements with the continuum elements in the computational mesh. In the following sub-sections, the implementation of bond elements as well as truss elements in the GIMP method is illustrated.

4.1 Bond element

The main feature of the bond element is the slip between truss element and continuum element. Each bond element has a point called anchor point located along the corresponding truss element at which slip values are calculated, see Figure 4.

Slip at an anchor point is given by

(30)s=uat(ξat)uac(ξac)=Nt(ξat)utNc(ξac)uc=Nt(ξat)Nc(ξac)utuc=Bb(ξa)utuc,
where subscript a indicates an anchor point and subscripts t and c refer to truss element and continuum element, respectively. Bond stresses induced by the slip vector are given by
(31)τxτy=Cx00Cysxsy,
where Cx is the constitutive characteristic relating bond stress in the longitudinal direction to the slip value sx. In general, Cx is non-linear with ∂τx/∂sx = 0 at maximum bond strength. However, for simplicity in this work, a constant Cx is assumed giving a linear relationship between τx(sx) and sx. On the other side, lateral slip, sy, is avoided by a penalty term (assigning to Cy a large value). The bond element has no thickness, hence, no geometric non-linear term exists in the formulation of its stiffness matrix. The stiffness matrix, Kb, and vector of internal forces, fb, of the bond element are
(32)Kb=Bb(ξa)TCbBb(ξa)Ab,
(33)fb=Bb(ξa)TτbAb,
with Ab as the area of the bond element, which is equal to the corresponding truss cross-sectional circumference multiplied by length of bond element. Kb and fb are assembled to K, Eq. (5), and f, Eq. (4), respectively, based on the nodal indices of the corresponding truss and continuum elements. Due to the fact that embedded truss elements have no external constraints or boundary conditions, a minimum of two bond elements is required to provide stability. The more bond elements are used, the more accurate is the result but also the more computational cost has to be paid.

The local coordinates of the anchor points in the continuum elements are calculated using the Newton–Raphson method in the same manner shown in Section 2.2, Figure 3. By means of this bond element, truss elements can be embedded in the continuum with own nodes independent of the continuum discretization. This is the main advantage, which allows coupling the MPM with truss elements as at each time step truss elements are moving over a fixed mesh.

4.2 Truss element

The truss element formulation including geometrical non-linearity with small strain is used assuming that the strain is small enough to keep the truss cross-sectional area and length constant.

(34)AA0andll0.
Eq. (34) states that the cross-sectional area and the length in the initial configuration (A0 and l0, respectively) are approximately equal the cross-sectional area and the length at any configuration (A and l, respectively). Hence, the virtual work equation of the truss element considering only the axial stress component, σxx, yields
(35)A0l0δϵxxσxxdx=δuTfext.

Eq. (35) is non-linear in terms of the nodal displacements. After linearization, the truss element vector of internal forces, ft, and stiffness matrix, Kt, are given by

(36)ft=A0l0BL(u)σxx,
(37)Kt=EA0l0BL(u)BLT(u)+A0σxxl0M,
where
(38)BLT=1l01+u2u1l0v2v1l01+u2u1l0v2v1l0,
where E is the Young's modulus, u1 and u2 are the horizontal displacements at Nodes 1 and 2, respectively, and v1 and v2 are the vertical displacements at Nodes 1 and 2, respectively. u1, u2, v1, and v2 result from the previous iteration. The matrix M is given by
(39)M=1010010110100101.

For more details on the truss element formulation (see De Borst et al., 2012).

4.3 Coupling truss elements with the GIMP method

In Sections 4.1 and 4.2, a brief explanation of the formulation of bond and truss elements applied in FEM is shown. For application in the MPM, two main features have to be addressed:

  1. Truss elements are moving over a fixed mesh.

  2. Global displacement vector u is set to 0 every time step, and new vector with new size and new nodal indices for elements are created.

Figure 5 shows an example of truss elements moving over a fixed mesh from time t to t + Δt. After solving the system in each iteration, the nodal displacement vectors of the truss elements are updated, which means that the positions of the truss elements and the corresponding anchor points change, as a result, the continuum elements of the fixed background mesh in which the truss elements are embedded may change or at least the local coordinates of the anchor points in the continuum elements will change and need to be evaluated each time step. Slip values at the anchor points from the previous time step should be added to the current slip values since the global displacement vector is set to zero at each time step and the current slip values are only calculated from the beginning of the current time step. As mentioned above that the size of the global displacement vector changes each time step as new continuum elements are activated and others are deactivated, the elements' global indices change as well. For this reason, the nodal displacements of the truss elements from the previous time step become history variables and are added to the current nodal displacements.

The description above is summarized by the following steps:

  1. Background element for each anchor point has to be found every time step.,

  2. Local coordinates of the anchor points in the corresponding element have to be evaluated every time step.

  3. Slip values from previous time step are added to the current values.

  4. For truss elements, displacement values from previous time step are added to the current values.

5. Numerical examples

In order to validate modeling reinforced materials by the GIMP method, four numerical examples are presented. The first one models a cantilever beam, the second example describes a clamped beam, the third simulation investigates a reinforced material tension bar, and the fourth model problem evaluates a beam-column frame structure.

5.1 Cantilever beam

The cantilever beam is subjected to its self-weight. The initial beam length and depth are 16 cm and 4 cm, respectively. Initially, the beam consists of 224 elements each with equally distributed 3 × 3 material points which gives a total of 2016 material points. Plane stress 2D continuum elements are used. The beam is described by a Young's modulus, E = 100 N/cm2 and Poisson's ratio, ν = 0, with mass density, ρ = 0.0165 kg/cm3. 32 truss elements are used with a cross-sectional area of 0.02 cm2 placed at a distance of 1 cm from the top of the beam. The reinforcement bar is characterized by Young's modulus, E = 10,000 N/cm2. The self-weight of truss elements is not included. The constitutive characteristic of the bond elements is Cx = 2500 N/cm3, see Eq. (31). The load is increased incrementally with a load step equal to 0.05 up to load factor of 1.0. On the other side, the same model is simulated by the FEM with 4-node and 9-node 2D plane stress continuum elements. Figure 6 shows the cantilever beam modeled by the GIMP method and the FEM in its final configuration. The crosses indicate the supports and the blue cells in the MPM model are the activated elements.

The absolute average displacements of the beam tip in vertical and horizontal directions and the axial stresses of the truss elements are shown in Figures 7 and 8, respectively. The figures show that the results of the GIMP method coupled with truss- finite element fit very well the results of the FEM. The obtained displacements yield a very small shift compared to the FEM, as the material points, for which the displacements are plotted, are located inside the beam and not exactly at the beam tip such as the nodes in the FEM. The plot shown in Figure 8 exhibits a small shift of the stress values in the left part of the beam, up to truss element Number 10, and, then, it overlies the stresses obtained by FEM. This small shift results from the high variation of the number of material points per elements in the left part of the beam, whereas the right part is almost straight, less curvy than in the right and the number of material points per element is more uniform. Figure 9 presents the distribution of the material points' stresses within the cantilever beam. As can be seen, the stress variation of the material points is smooth and no stress oscillations are observed. In Table 1, the number of iterations for load steps are shown. As can be seen, the convergence rate is almost quadratic, which is an indication of a correct implementation.

In Figure 10, the average vertical displacements of the material points at the cantilever beam tip are compared with the average displacements obtained from FEM at the same positions for each case. The FEM is taken as a reference, with 4-node 262,144 (256 × 1,024) continuum elements, with 1/element size equals to 64. The number of truss elements equals the number of continuum elements in the x-direction for all cases. The average material point displacements have been normalized to the displacements given by the reference FEM. The convergence is investigated using different grid element sizes each with different numbers of material points per element. Figure 10 shows that by reducing the grid element size, the normalized displacement values approach one. The displacements in the finer meshes are very similar for different numbers of material points. Same result has been shown in Charlton et al. (2017), which indicates that reducing element size has more influence on the accuracy of the result than changing the number of material points per element.

5.2 Clamped beam

A reinforced clamped beam is subjected to a uniform load applied at the mid span of the top surface with a magnitude of 500 N/cm. The loaded material points and the mapped forces in the initial configuration are shown in Figure 11. The reinforcement bars at the top and bottom are depicted as well. The crosses indicate the supports and the blue cells are the initially activated elements.

The initial beam length and depth are 60 cm and 12 cm, respectively. Initially, the beam consists of 720 elements each with equally distributed 3 × 3 material points, which gives a total of 6,480 material points. Plane stress 2D continuum elements are used. The beam is described by a Young's modulus, E = 5000 N/cm2 and Poisson's ratio, ν = 0.3. Top and bottom reinforcements are placed at a distance of 1.75 cm from the top and bottom of the beam, respectively, each with 59 truss elements of a cross-sectional area of 0.01 cm2. The reinforcement bars are associated with a Young's modulus, E = 10,000 N/cm2. The self-weight of truss elements is excluded. The bond elements exhibit a longitudinal stiffness, Cx = 2500 N/cm3. The load is increased incrementally with a load step of 0.05 up to load factor of 1.0.

On the other side, the same model is simulated by the FEM with 4-node and 9-node 2D plane stress continuum elements. Figure 12 shows the cantilever beam simulated by the GIMP method and the FEM in its final configuration. The red and green color of truss elements indicate tension and compression, respectively.

The axial stresses of truss elements for the GIMP method and the FEM are shown in Figure 13. The stresses for the upper and lower reinforcement bars are plotted showing the variation between tension and compression. The stresses obtained from the GIMP method overly the stresses from FEM. The small variation of the FEM with 9-node continuum elements is due to the difference in the way the load is applied, as for the GIMP and FEM with 4-node elements the load is applied on the nodes at the corners; however, in the FEM with 9-node elements, the load is also applied on the nodes at the midpoint of the element edge as well. Hence, the stress concentration at the region of applied load is different.

It is shown for this example as well that the results of the GIMP method are in a good agreement with the results of the FEM.

For this example, as for the previous one, the stress variation of the material points is smooth and no stress oscillations are observed, see Figure 14. Also, the convergence rate is quadratic, as can be seen in Table 2, which emphasizes that the inclusion of reinforcement is implemented correctly.

The influence of changing the grid elements size as well as the number of material points per element is investigated in this example. In Figure 15, the displacement values of the material points close to the mid span are compared to the displacement at the same positions for each case obtained from the FEM. The FEM is taken as a reference, with 4-node 184,320 (192 × 960) continuum elements, with 1/element size equals to 16. The length of truss finite elements in the top and bottom reinforcement bars equals the length of continuum elements in the x-direction for all cases. The displacements of the GIMP models have been normalized to the displacements given by the reference FEM. Figure 15 shows the same result as in the previous example, that by reducing the grid elements size, the normalized displacement values approach one. The same conclusion is made that reducing the element size influences the result more than changing the number of the material points.

5.3 Elasto-plastic reinforced material tension bar

In the previous examples, coupling of the truss element with the MPM is validated using the FEM and the convergence rate is observed. In this present example, a reinforced material tension bar is modeled in different situations using an elasto-plastic constitutive model to further show the interaction between the reinforcement and the material points. Three models are used and depicted in Figure 16. Model 1 is an unreinforced bar, which is described by two materials, elastic material in the outer parts and elasto-plastic material in the middle segment. Model 2 is a reinforced bar that consists of two elastic parts connected by the reinforcement rebar. Model 3 is a reinforced bar which consists of two materials, elastic material in the outer parts and elasto-plastic material in the middle segment with a reinforcement rebar. The bar length and depth are 60 cm and 7 cm, respectively, with thickness of 3 cm. The element size used is 1 × 1 cm with equally distributed 3 × 3 material points. Plane strain 2D continuum elements are used. The bar is described by a Young's modulus, E = 50,000 N/cm2 and Poisson's ratio, ν = 0.0. The yield stress of the elasto-plastic part is 15 N/cm2. The reinforcement rebar is placed at the middle of the bar with total length of 40 cm made of 40 truss elements of a cross-sectional area of 1 cm2. The reinforcement rebar is associated with a Young's modulus, E = 50,000 N/cm2. A large value is given to the bond elements' longitudinal stiffness assuming a rigid bond connection. The bar is fixed at the left edge and prescribed displacement is applied at the right edge. The displacement is increased incrementally with a step of 0.05 up to 1.0.

The variation of the reaction force with the applied displacement is depicted in Figure 17. The plot for Model 1 shows the typical elasto-plastic behavior. The reaction force of Model 2 increases linearly with the applied displacement as linear elastic material is used. The plot of Model 3 consists of two parts: the first part with the stress value below the yield stress. The stiffness is the same as for Model 2 as the stiffness contribution of the reinforcement rebar is small comparing to the stiffness of the bar itself. The second part of the plot with the stress value exceeds the yield stress and the bar enters the plastic regime. As the stiffness of the plastic material drops to almost zero, the stiffness in the second part of the plot of Model 3 is almost equal to the stiffness of Model 2 and the two lines are parallel.

The axial stresses of the truss elements in Model 2 at the final load step are presented in Figure 18. The axial stresses in the unbonded truss elements are equal to the final reaction force value in Figure 17 for Model 2, with the cross-sectional area is equal to 1 cm2.

The results in Figures 17 and 18 indicate correct implementation of the truss finite element in the MPM.

5.4 Beam-column frame structure

A frame structure is subjected to uniform load in the x-direction at the top left corner with magnitude of 3 MN/m. The frame is uniform in thickness of 1 m with 6 m total height and 7 m total width. Initially, the frame consists of 612 elements, each with equally distributed 3 × 3 material points which gives a total of 5,508 material points. Plane stress 2D continuum elements are used. The frame is described by a Young's modulus, E = 100 MN/m2 and Poisson's ratio, ν = 0.3. A total of 216 truss elements are used with a cross-sectional area of 0.002 m2 placed at a distance of 0.2 m from the frame surfaces. Figure 19 depicts the reinforced frame in the initial configuration. The reinforcement bars exhibit a Young's modulus, E = 20,000 MN/m2. The constitutive characteristic of the bond elements is Cx = 5000 MN/m3. The load is increased incrementally with a load step of 0.1 up to load factor of 1.0. On the other side, the same model is simulated by the FEM with 4-node 2D plane stress continuum elements. Figure 20 shows the reinforced frame modeled by the GIMP method and the FEM in its final configuration. The red and green color of truss elements indicates tension and compression, respectively.

The axial stresses of truss elements for the GIMP method and the FEM are shown in Figure 21. The stresses for the reinforcement bars at the outer surface of the frame are plotted showing the variation between tension and compression. The stresses obtained from the GIMP method are in good agreement with the stresses obtained from the FEM. A main difference in the stress value of the truss element at the top left corner is seen at the region of applied load. This difference is due to the manner the applied load is distributed over the nodes in the GIMP method and the FEM, see Figure 20. As the applied load is distributed in different ways, the resulting stress in the continuum and the truss elements in the near field will be different as well. The distribution of material points' stresses within the frame structure is depicted in Figure 22. The figure shows smooth stress variation and no oscillations.

In order to show the effect of the reinforcement, the frame structure is modeled in different scenarios: using elastic material and the Von-Mises elasto-plastic constitutive formulation with yield stress of 5 MN/m2, both with and without reinforcement. The frame is subjected to 5.4 MN/m at the left upper corner as shown previously. The horizontal displacements in all cases at the upper right corner are presented in Figure 23. As shown in the figure, for all cases, the presence of the reinforcement increases the stiffness of the structure and reduces the displacement. However, the influence of the reinforcement is larger in case of the elasto-plastic material as the stiffness of the continuum reduces significantly after yielding. The deformed shapes of the frame structure using the elasto-plastic constitutive approach with and without reinforcement are shown in Figure 24. The unreinforced frame undergoes extremely large deformation in comparison to the reinforced frame which illustrates the influence of the reinforcement on the global structural response.

6. Conclusions

Modeling reinforced materials discretized by the GIMP method is presented in this paper. The reinforcements are described using truss- finite elements. Beam elements can be applied as well if the bending stiffness is dominating. Truss elements are coupled to continuum elements in the background mesh using bond elements. By means of the bond elements, truss elements, which are moving over the fixed mesh, are linked to continuum elements at each time step. The displacement values are added from the previous time step as in the GIMP method the global displacement vector is reset at each time step. The same is also applied to the bond elements, where the slip values are added from previous time step.

Looking at the presented numerical examples, the calculated axial stresses of the truss elements coupled with the GIMP method are in good agreement with the results from the FEM. Small differences have been seen in the areas with large deformations, where a severe change of the material points number per elements is present and from the other side where the finite elements in the FEM are distorted.

The numerical examples have shown quadratic convergence during the simulations. It has been shown as well that by reducing the background mesh size, the result approaches toward the correct solution. Both outcomes indicate a correct implementation of the coupling of the truss finite elements in the GIMP method.

The numerical examples used in this paper are chosen in a way that they can be modeled by both FEM and GIMP method to serve in the validation of coupling the truss- finite elements in the GIMP method. This paper is introducing the method of coupling GIMP-truss FE that will be used for applications of severe deformations for which the FEM fails to yield correct results.

The constitutive relations used in this paper are linear elastic and simple elasto-plastic material models. This simplification is chosen because the focus of this publication is on coupling the GIMP method to truss elements regardless the material description used. Using different constitutive formulations will be presented in future works. The discussion is made for 2D problems for simplicity; however, the extension to 3D is straight forward.

Figures

Reinforced concrete plate under impact

Figure 1

Reinforced concrete plate under impact

Phases in a computational cycle in implicit MPM

Figure 2

Phases in a computational cycle in implicit MPM

Material point local coordinates algorithm

Figure 3

Material point local coordinates algorithm

Truss elements embedded in continuum elements

Figure 4

Truss elements embedded in continuum elements

Coupling the GIMP method with truss finite elements

Figure 5

Coupling the GIMP method with truss finite elements

Reinforced cantilever beam model

Figure 6

Reinforced cantilever beam model

Cantilever tip displacements

Figure 7

Cantilever tip displacements

Axial stresses of truss elements

Figure 8

Axial stresses of truss elements

Material points' stresses in x-direction of the cantilever

Figure 9

Material points' stresses in x-direction of the cantilever

Convergence of displacement with changing element size

Figure 10

Convergence of displacement with changing element size

Clamped beam

Figure 11

Clamped beam

Clamped beam in the final configuration

Figure 12

Clamped beam in the final configuration

Axial stresses of truss elements

Figure 13

Axial stresses of truss elements

Material points' stresses in x-direction of the beam

Figure 14

Material points' stresses in x-direction of the beam

Convergence of displacement with changing element size

Figure 15

Convergence of displacement with changing element size

Different models of the tension bar

Figure 16

Different models of the tension bar

Reaction force of the tension bar

Figure 17

Reaction force of the tension bar

Axial stresses of the truss elements in Model 2

Figure 18

Axial stresses of the truss elements in Model 2

Beam-column frame structure

Figure 19

Beam-column frame structure

Frame structure in the final configuration

Figure 20

Frame structure in the final configuration

Axial stresses of truss elements

Figure 21

Axial stresses of truss elements

Material points stresses magnitude of the frame structure

Figure 22

Material points stresses magnitude of the frame structure

Displacements in x-direction of the frame structure

Figure 23

Displacements in x-direction of the frame structure

Deformed shape of elasto-plastic frame structure

Figure 24

Deformed shape of elasto-plastic frame structure

Algorithm of implicit GIMP method

Figure A1

Algorithm of implicit GIMP method

Algorithm to evaluate K and f

Figure A2

Algorithm to evaluate K and f

Convergence rate of reinforced cantilever beam

StepIterRes. [N]StepIterRes. [N]
0.5511.2e−000.811.7e−00
0.5521.7e−010.821.5e−01
0.5531.3e−040.831.1e−04
0.5544.1e−080.844.7e−08
0.612.7e−000.8512.9e−00
0.621.8e−010.8521.6e−01
0.631.4e−040.8531.3e−04
0.645.3e−080.8546.0e−08
0.6512.0e−000.911.7e−00
0.6521.7e−010.921.4e−01
0.6531.3e−040.931.0e−04
0.6545.0e−080.944.2e−08
0.712.4e−000.9512.7e−00
0.721.7e−010.9521.5e−01
0.731.4e−040.9531.0e−04
0.745.7e−080.9544.2e−08
0.7511.5e−001.012.3e−00
0.7521.6e−011.021.7e−01
0.7531.2e−041.031.2e−04
0.7545.0e−081.044.5e−08

Convergence rate of reinforced clamped beam

StepIterRes. [N]StepIterRes. [N]
0.5513.9e+010.813.5e+01
0.5522.2e+000.822.2e+00
0.5536.7e−040.831.4e−03
0.5541.5e−100.841.3e−09
0.613.7e+010.8513.9e+01
0.622.0e+000.8522.1e+00
0.637.2e−040.8531.3e−03
0.643.2e−100.8541.3e−09
0.6517.1e+010.918.2e+01
0.6522.7e+000.921.2e+01
0.6531.3e−030.933.5e−02
0.6541.0e−090.945.9e−07
0.714.8e+010.9515.5e+01
0.724.0e+000.9523.2e+00
0.733.2e−030.9533.4e−03
0.746.7e−090.9547.1e−09
0.7514.1e+011.015.7e+01
0.7522.3e+001.023.6e+00
0.7531.2e−031.035.8e−03
0.7546.9e−101.042.3e−08

Appendix GIMP method algorithm

Figure A1 summarizes the main steps of the computation within a time step. Note that in the “update material point” step, the following computations are performed

  1. Update the displacement, Eq. (2).

  2. Update the volume, Eq. (21).

  3. Update CT of each material point, Eq. (12).

  4. Update material point domain, łp, Eq. (29).

Figure A2.

References

Bardenhagen, S. and Kober, E. (2004), “The generalized interpolation material point method”, Computer Modeling in Engineering and Sciences, Vol. 5, pp. 477-496.

Bathe, K.-J. (2006), Finite Element Procedures, Prentice-Hall, Englewood Cliffs.

Brackbill, J.U. and Ruppel, H. (1986), “Flip: a method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions”, Journal of Computational Physics, Vol. 65, pp. 314-343.

Charlton, T., Coombs, W. and Augarde, C. (2017), “igimp: an implicit generalised interpolation material point method for large deformations”, Computers and Structures, Vol. 190, pp. 108-125.

Chen, Z., Qiu, X., Zhang, X. and Lian, Y. (2015), “Improved coupling of finite element method with material point method based on a particle-to-surface contact algorithm”, Computer Methods in Applied Mechanics and Engineering, Vol. 293, pp. 1-19.

Coombs, W., Augarde, C., Brennan, A., Brown, M., Charlton, T., Knappett, J., Motlagh, Y. and Wang, L. (2020), “On Lagrangian mechanics and the implicit material point method for large deformation elasto-plasticity”, Computer Methods in Applied Mechanics and Engineering, Vol. 358, 112622.

De Borst, R., Crisfield, M., Remmers, J. and Verhoosel, C. (2012), Nonlinear Finite Element Analysis of Solids and Structures, Wiley, Chichester.

Elwi, A.E. and Hrudey, T.M. (1989), “Finite element model for curved embedded reinforcement”, Journal of Engineering Mechanics, Vol. 115, pp. 740-754.

Fleischhauer, R., Qinami, A., Hickmann, R., Diestel, O., Götze, T., Cherif, C., Heinrich, G. and Kaliske, M. (2015), “A thermomechanical interface description and its application to yarn pullout tests”, International Journal of Solids and Structures, Vol. 69, pp. 531-543.

Fuchs, A., Curosu, I. and Kaliske, M. (2020), “Numerical mesoscale analysis of textile reinforced concrete”, Materials, Vol. 13, p. 3944.

Ghorbani, M., Khiavi, M.P. and Moghaddam, F.R. (2009), “Nonlinear analysis of shear wall using finite element model”, International Journal of Civil and Architectural Engineering, Vol. 3, pp. 375-379.

Guilkey, J. and Weiss, J. (2003), “Implicit time integration for the material point method: quantitative and algorithmic comparisons with the finite element method”, International Journal for Numerical Methods in Engineering, Vol. 57, pp. 1323-1338.

Hartl, H. (2002), “Development of a continuum-mechanics-based tool for 3D finite element analysis of reinforced concrete structures and application to problems of soil-structure interaction”, PhD thesis, Graz University of Technology, Institute of Structural Concrete.

Häussler-Combe, U. (2015), Computational Methods for Reinforced Concrete Structures, Ernst & Sohn, Berlin.

Häussler-Combe, U., Chihadeh, A. and Shehni, A. (2019), “Modelling of discrete cracks in reinforced concrete plates with the strong discontinuity approach (sda)”, Advances in Engineering Materials, Structures and Systems: Innovations, Mechanics and Applications, Cape Town.

Häussler-Combe, U., Shehni, A. and Chihadeh, A. (2020), “Finite element modeling of fiber reinforced cement composites using strong discontinuity approach with explicit representation of fibers”, International Journal of Solids and Structures, Vol. 200, pp. 213-230.

Hering, M., Bracklow, F., Scheerer, S. and Curbach, M. (2020), “Reinforced concrete plates under impact load—damage quantification”, Materials, Vol. 13, p. 4554.

Lian, Y. and Zhang, X. (2018), “A coupled finite element material point method for large deformation problems”, Advances in Computational Coupling and Contact Mechanics, Vol. 11, pp. 251-288.

Lian, Y., Zhang, X. and Liu, Y. (2011a), “Coupling of finite element method with material point method by local multi-mesh contact method”, Computer Methods in Applied Mechanics and Engineering, Vol. 200, pp. 3482-3494.

Lian, Y., Zhang, X., Zhou, X. and Ma, Z. (2011b), “A femp method and its application in modeling dynamic response of reinforced concrete subjected to impact loading”, Computer Methods in Applied Mechanics and Engineering, Vol. 200, pp. 1659-1670.

Liang, W. and Zhao, J. (2020), “Coupled MPM/DEM multiscale modelling geotechnical problems involving large deformation”, 16th Asian Regional Conference on Soil Mechanics and Geotechnical Engineering, ARC, 2019.

Mahmoud, A.M. (2016), “Finite element modeling of steel concrete beam considering double composite action”, Ain Shams Engineering Journal, Vol. 7, pp. 73-88.

Ngo, D. and Scordelis, A.C. (1967), “Finite element analysis of reinforced concrete beams”, ACI Journal Proceedings, Vol. 64, pp. 152-163.

Ogden, R.W. (1997), Non-linear Elastic Deformations, Dover Publications, New York.

Schäfer, H. (1975), “A contribution to the solution of contact problems with the aid of bond elements”, Computer Methods in Applied Mechanics and Engineering, Vol. 6, pp. 335-353.

Sulsky, D. and Schreyer, H. (1996), “Axisymmetric form of the material point method with applications to upsetting and Taylor impact problems”, Computer Methods in Applied Mechanics and Engineering, Vol. 139, pp. 409-429.

Sulsky, D., Chen, Z. and Schreyer, H. (1994), “A particle method for history-dependent materials”, Computer Methods in Applied Mechanics and Engineering, Vol. 118, pp. 179-196.

Sulsky, D., Zhou, S.-J. and Schreyer, H. (1995), “Application of a particle-in-cell method to solid mechanics”, Computer Physics Communications, Vol. 87, pp. 236-252.

Wang, B., Vardon, P., Hicks, M. and Chen, Z. (2016), “Development of an implicit material point method for geotechnical applications”, Computers and Geotechnics, Vol. 71, pp. 159-167.

Zhang, Y.-G., Lu, M.-W. and Hwang, K.-C. (1994), “Finite element modeling of reinforced concrete structures”, Finite Elements in Analysis and Design, Vol. 18, pp. 51-58.

Zhang, X., Sze, K. and Ma, S. (2006), “An explicit material point finite element method for hyper-velocity impact”, International Journal for Numerical Methods in Engineering, Vol. 66, pp. 689-706.

Zhang, X., Chen, Z. and Liu, Y. (2016), The Material Point Method: A Continuum-Based Particle Method for Extreme Loading Cases, Academic Press, London.

Acknowledgements

The presented work is funded by the German Research Foundation (DFG) within Research Training Group GRK 2250 Mineral bonded composites for enhanced structural safety, Sub-project B4. This financial support is gratefully acknowledged by the authors.

Corresponding author

Michael Kaliske can be contacted at: Michael.Kaliske@tu-dresden.de

Related articles