Multi-scale representation of plastic deformation in fiber-reinforced materials: application to reinforced concrete

Here we present a multi-scale model to carry out the computation of brittle composite materials reinforced with fibers, and we show its application to standard reinforced concrete. The computation is carried out within an operator-split framework on the macro-scale, which allows for different failure mechanisms to develop in separate phases, as both the concrete and the bond-slip exhibit non-linear behavior. The computations on the micro-scale are performed for each constituent separately. The reinforcement is taken to be linear elastic, and the bond-slip is handled as a plastic deformation. The standard elastoplastic procedure is used to compute the bond stresses, combined with the X-FEM methodology to give the global representation of slip. The crack development in concrete, on the other hand, is described with a damage model with exponential softening, where ED-FEM is used to represent localized failure. A numerical example is shown to test the developed methodology.


INTRODUCTION
To model fiber-reinforced composite materials, one has to take into account the behavior of each of the constituents, and their interaction.In the case of fiber-reinforced concrete (FRC), we are dealing with concrete, steel fibers, and the interface between them.We can consider standard reinforced concrete as a special case of such a material, with its specific features.First of all, there are not more than a few reinforcement bars embedded in concrete (compared to a large number of randomly oriented short fibres in FRC), that is usually straight, larger in diameter, and anchored at its ends.This facilitates the numerical implementation, as the steel reinforcement can lie along the finite element's edge, and the anchorage prevents the pull-out of the bar.
One of the most important choices to make when dealing with the interaction of two different materials (steel and concrete in this case), is the treatment of the interface between them.Many authors have opted for the insertion of a special interface element, i.e.Ožbolt et al. (2002), Kohnehpooshi and Jaafar (2017), while others have opted for a layered description of reinforced concrete, where each layer can represent a different material, i.e.Jukić et al. (2014) or Šćulac at al. (2014).An example of the utilisation of a special interface element can be found in the work of Dominguez et al. (2005), where a degenerated Q4 element is used to handle the relative displacements between the steel and concrete.Despite the fact that such a description of slip can provide many benefits, it can be gruesome to handle and implement.A simplification regarding the constitutive law implementation in this kind of Q4 element can be found in Ibrahimbegovic et al. (2010), where a Drucker-Prager non-associative plasticity model is used, but still the issue with introducing a special element persists.
That is why we have decided to simplify the handling of bond-slip, but without sacrificing the model capabilities.Building on top of the work of Ibrahimbegovic et al. (2010), we are proposing a novel and efficient way to take into account the slip between the reinforcing bar and the surrounding concrete.Instead of having a special finite element for every constituent (i.e. a CST element for concrete, a truss bar for steel, and a Q4 element for bond-slip), we have encapsulated all three ingredients in a single finite element, without the need to worry about the nodal connections between the concrete, steel, and interface part.
In the next chapter, the formulation of the model is presented, with its specific micro−macro treatment of the interaction of three different constituents, and after that we show the results of the numerical simulation of a reinforced concrete specimen.

Displacement field approximation
We start by defining the bond-slip as the relative displacement between steel and concrete where   represents the slip on the interface along the reinforcement bar,   is the concrete displacement, and   is the steel displacement.This kind of description allows us to write the displacement field of a single finite element as sum of the standard and the enriched part, according to the X-FEM methodology (Fries and Belytschko, 2010) Equation ( 2) exploits the partition of unity property of standard shape functions   , that allows for any other function (in this case the enrichment function ) to be reproduced exactly on the element domain   , (Melenk and Babuška, 1996).The enrichment function  is equal to the Heaviside function that takes the value one in the elements that contain the reinforcement bar (that is located at coordinate ), and zero in all the other elements Since in our case the fiber coincides with the lower edge of the enriched element, as shown on Figure 1, we can simplify equation ( 2) by exploiting the properties of the enrichment function Figure 1: Finite element with enriched degrees of freedom (the reinforcement bar is shown in bold).
Multi-scale representation of plastic deformation in fiber-reinforced materials: application to reinforced concrete Tea Rukavina et al.
Latin American Journal of Solids and Structures, 2019, 16(7 Thematic Section), e141 3/11 The multi-scale characteristics of our model can be observed in the treatment of the failure mechanisms on two levels, macro and micro.The macro-scale computations are divided in two phases, according to the macro-scale operator-split procedure: the macro-global phase (where the concrete and the steel are taken into account), and the macro-local phase (where the redistribution of slip is computed, with the steel bar acting as a coupling term).But, that is not all, since every constituent is taken care of within a micro-scale computation, which can have its own operator-split solution procedure, as in the case of concrete and bond-slip, giving rise to micro-local and micro-global computations.All of this is important for ensuring the best possible description of the behavior of each constituent, with its own peculiar phases related to specific failure mechanisms.Also, it is important to carry out all the computations in a defined order, based on experimental observations on composite specimens.

Micro-scale computations
To be able to couple the three constituents, each one of them has to provide the stress value that is computed from its own constitutive law and kinematics (Figure 2).
Steel computation.Since the steel reinforcement is taken to be linear elastic, its constitutive equation is In equation ( 5),   is the stress,   is the Young's modulus for steel, and   is the strain that is defined as where the steel displacements are computed according to (1).Concrete computation.The chosen material model for concrete is the damage model with isotropic hardening, that can take into account the formation of micro-cracks in the bulk, combined with a localized failure part that describes the opening of the macro-crack at the discontinuity.We will give here just a short description of the model, since all the details regarding the formulation and implementation can be found in the work of Brancherie and Ibrahimbegovic (2009), or Do et al. (2017) for applications in dynamics.
The damage functions  are used to check the admissibility of the stress state, and for the hardening part, this function is defined as where   is the undamaged elastic compliance tensor for the bulk material.It is computed as the inverse of the elastic constitutive matrix:   = (  ) −1 .In the above equation, ||  ||   represents the norm of   in the stress space,   is Young's modulus for concrete,   is the stress at the first cracking, and  �  is the stress-like hardening variable that controls the damage threshold evolution.
Multi-scale representation of plastic deformation in fiber-reinforced materials: application to reinforced concrete Tea Rukavina et al.
Latin American Journal of Solids and Structures, 2019, 16(7 Thematic Section), e141 4/11 Since the softening part of the response is controlled by an anisotropic multi-surface damage model, we have to consider each direction at the discontinuity surface separately.This kind of model can account for the crack opening in mode I (in the normal direction) and mode II (in the tangential direction).In the following equations the subscripts 1 and 2 denote the direction, and  and  stand for the normal and tangential vector on the discontinuity, respectively.The damage functions for the discontinuity are then where  is the traction at the discontinuity,  �  is the ultimate stress in tension, and  �  is the ultimate stress in shear.
The stress-like softening variable  �  is the same for both directions, and it depends on the chosen softening law, i.e. linear or exponential softening.
The constitutive equation relating the traction at the discontinuity  and the crack opening   is based on a traction-separation cohesive law where  � is the damage compliance tensor for the discontinuity.The handling of the displacement jump inside the element is done through the introduction of incompatible mode functions .The total displacement field of a single element is then a sum of the standard and the incompatible part from which it follows that the strain can be written as where  is the matrix of the shape functions' derivatives, and  ˜ consists of the regular part of the incompatible shape functions' derivatives.
The final value of stress in concrete can then be computed from and the value of the elasto-damage modulus   has to be chosen according to the current phase of the computation, which is determined from the damage functions described in (7-8).
Bond-slip computation.The main novelty in this work is the description of bond-slip through a one-dimensional elasto-plastic computation.Since the bond-slip law is given as a relationship between the bond stress   and the slip   , instead of the plastic strain   that is usually employed in this kind of computations, we will have the plastic slip  , .The computation is carried out within an operator-split solution procedure (see Ibrahimbegovic, 2009) where the evolution equations of internal variables are solved in the local phase, on the level of the Gauss numerical integration point of each element, and the equilibrium equations are solved globally, on the level of the whole structure.We start by an elastic trial step, in which the plastic multiplier   is equal to zero, and the values of all internal variables are frozen (they take the value from the previous time step) Multi-scale representation of plastic deformation in fiber-reinforced materials: application to reinforced concrete Tea Rukavina et al.
Latin American Journal of Solids and Structures, 2019, 16(7 Thematic Section), e141 5/11 In the above equations,   represents the hardening variable,   is the stress-like hardening variable for bond-slip, and  ,ℎ is the hardening modulus.The trial value of bond stress is then computed as where   is the bond-slip tangent modulus.The term in the brackets represents the elastic slip  , , since the total slip can be represented as the sum of the elastic and the plastic part:   =  , +  , .The total slip is computed from the nodal values of slip,    , which are obtained from the local iteration of the macro-level computation The bond-slip shape functions    are actually the product of the standard shape functions for concrete, and the X-FEM enrichment function , as given in equation ( 2).For our case it follows since the enrichment function is a Heaviside function.Here,   () are the linear shape functions for a truss bar.
To check if the trial value of stress is admissible, we introduce the yield function where   is the limit value of the bond stress.We have two cases to consider: if the trial value of the yield function is negative or zero, the step is indeed elastic and the trial values are accepted as final.If not, we have to proceed to the plastic step to correct the value of bond stress due to the plastic slip activation.
In the plastic step, we have to compute a new value for the plastic multiplier which is used to update the internal variables according to The final value of stress is then The elastoplastic tangent modulus in the elastic phase is equal to   , and in the plastic phase, it has the value

Macro-scale computations
By inserting the approximation of the virtual displacement field into the weak form of the equilibrium equations, and after gathering the terms related to the concrete displacements   , and the ones related to the bond-slip displacements   , we are left with two systems of equations Here,  denotes the number of elements, and  the number of enriched elements (the ones containing the steel bar).The global equation ( 22a) is solved with a fixed slip, and in the local one (22b) the redistribution of slip is computed.The external force vector and the internal force vectors for concrete, steel, and bond-slip, are defined in the following manner  ,, = ∫         ,, = ∫     ,       ,, = ∫     , In equation ( 23),   is the matrix of derivatives of the shape functions for the steel bar,   is the cross-sectional area of the bar, and   is the bond-slip area that is computed as the circumference of the bar (it is the part of the steel bar in contact with concrete).In the external force vector,  represents the volumetric forces, and  is the imposed traction acting at the boundary   .
By linearizing ( 22), we obtain the following system of equations where the tangent stiffness matrices are defined as The system ( 24) is handled within the macro-scale operator-split solution procedure, that gives us two equations to solve sequentially where  ̂ is the condensed stiffness matrix Equation ( 26a) is the macro-global equation that gives us the value of the concrete displacement for the fixed value of slip.With the value of   in hand, we proceed to the macro-local equation ( 26b), where we update the value of the slip   .
Multi-scale representation of plastic deformation in fiber-reinforced materials: application to reinforced concrete Tea Rukavina et al.

NUMERICAL RESULTS
The proposed methodology was implemented in FEAP (Finite Element Analysis Program), developed at Berkeley by professor R. L. Taylor (2014).All the described constituents are implemented into a single finite element that allows for an efficient computation of their interaction.
To test the performance of the developed formulation, we have performed a tension test on a concrete specimen of dimensions 400x100x100 mm 3 , reinforced by a steel bar of diameter  =16 mm.Since our model is 2D (as shown of Figure 3.), the area of the bar and the bond-slip area are divided by the thickness of the specimen.We consider the steel bar to be fixed at both ends, as it is anchored in concrete.The material properties used in the model are listed here: • for concrete:   = 45700MPa,   = 0.2,  �  = 3.5MPa,  ℎ = 1000MPa,  �  = 4MPa,  � , = 3.8MPa,   = 20; • for the steel bar:   = 210000MPa,   = 2.01mm 2 ; • for the bond:   = 30N/mm 3 ,   = 6MPa,  ,ℎ = 0.03N/mm 3 ,   = 0.5mm.
Among the properties listed above,   is the Poisson's ratio for concrete, and   is the parameter that controls the exponential softening. � , is the ultimate stress for a strip of weakened elements in the middle of the specimen where the crack will appear.The specimen is fixed at the left-hand side, and there is an imposed displacement  � = 0.5 mm acting on the right-hand side.The finite element mesh shown on Figure 4 consists of 800 CST elements, where the triangle sides of each element have the dimension 10 mm.There are 40 enriched elements along the reinforcement bar.
The force-displacement diagram shown on Figure 5. plots the reaction in the -direction at the left-hand side against the imposed displacement on the right-hand side.We can observe several particular phases of the composite behavior: the linear elastic phase is followed by a hardening phase when micro-cracks start to appear in concrete.After the ultimate load of about 500 N is reached, the crack develops in the weak zone in the middle of the specimen, so we enter the softening phase.After a while, the reaction starts ascending, due to the redistribution of stresses and the reinforcement activation.This resembles the typical diagrams that are found in the literature, which are obtained for a specimen with a single crack.concrete Tea Rukavina et al.
Latin American Journal of Solids and Structures, 2019, 16(7 Thematic Section), e141 8/11 On Figure 6.we can see the distribution of slip along the reinforcement bar, where the largest slips are taking place near the crack, and at the ends the slip is zero.On the left of the crack,   has a positive value, and on the left side it is negative, because the reinforcement is moving (relatively to the concrete) from the ends to the centre (in other words, towards the crack).The bond stress has a very similar distribution (Figure 7a), according to the chosen bond-slip law for the interface.The evolution of bond stress in time gives rise to an interesting plot (Figure 7b), where the dashed lines represent the elements left of the crack, and the dotted lines represent the ones on the right.We can see that the elements on the left have a positive value of bond-stress, and for the elements on the right the bond stress is negative.The bold solid line represents the value of bond stress in the cracked element, that is very near to zero, since it is the inflection point of the curve shown on  Related to that is the plastic slip development shown on Figure 8b, where it can be observed that plasticity occurs only at the end of the analysis, in five elements that are nearest to the crack (two on the left side, and three on the right).The same can be seen on Figure 8a, where the plastic slip is plotted in Gauss points along the reinforcement bar.The plastic slip is equal to zero in most of the enriched element, since they have not yet entered the plastic phase and the interface is behaving elastically.Figure 9. gives the representation of the crack opening in concrete, that reaches the value of   =0.47 mm at the end of the analysis.As it has been already explained, the macro-crack starts to open when the material reaches the ultimate stress, which happens around time 0.25.There are also two small cracks in the elements right next to the main crack, but their size is negligible compared to the main crack.It is interesting to point out that the crack opening in concrete coincides with the slip activation, which can be compared on Figure 7b and Figure 9b.This is in accordance with our model assumption that there is no slip when the concrete and the steel have the same displacement.So, the crack opening in concrete,   , is giving rise to the difference between the concrete and steel strain field, that is, in turn, activating the bond-slip   .Moreover, we can test this notion by plotting the crack opening evolution and the slip evolution on the same diagram, as has been done on Figure 10.We are plotting the values of   in the Gauss point of the cracked element (represented by the solid line), and the values of   in left and right node of the same element (dashed and dotted line, respectively).While we have already noted that the two nodes in question have the displacements of opposite sign, it is even more compelling that, quantitatively, the absolute values of the two slips add up to give nearly the same value as the crack opening at a given time.The sum of two slips is represented by a dotted line on Figure 10.

CONCLUSIONS
The methodology presented in this paper allows us to model the failure mechanisms in reinforced concrete happening on the different scales of the material, by taking into account the behavior of each of the constituents.The bond-slip is handled within a 1D elastoplastic framework, and we have shown that this simplification does not affect the model predictivity; hence it represents an improvement compared to previous works based on interface elements.
We have given the multi-scale set-up for this kinds of problems, where multi-level operator-split solution procedures can be performed during the computation.In the numerical example we have shown the relationship between the crack opening and the bond-slip, which corresponds to the realistic processes happening during the failure of fiber-reinforced composites.
The next step would be to extend the methodology to take into account the complete pull-out of the fiber, in order to capture the failure mechanisms taking place when the fiber's ends are not anchored in the surrounding material.

Figure 2 :
Figure 2: Constitutive law for each material: (a) damage model for concrete; (b) linear elasticity for steel; (c) elastoplasticity for bond-slip.

Figure 3 :
Figure 3: Geometry of the reinforced concrete specimen in 2D.

Figure 4 :
Figure 4: Finite element mesh with enriched elements shown in grey.

Figure 5 :
Figure 5: Force-displacement diagram for the tension test.

Figure 6 :
Figure 6: Distribution of slip   along the reinforcement bar at the end of the analysis (nodal values).

11 Figure 7 :
Figure 7: Bond stress: (a) distribution of   along the reinforcement bar (values at Gauss points); (b) evolution of   in time for all enriched elements.

Figure 8 :
Figure 8: Plastic slip: (a) distribution of  , along the reinforcement bar at the end of the analysis (values at Gauss points); (b) evolution in time for the enriched elements that have entered the plastic phase.

Figure 9 :
Figure 9: Crack opening in concrete: (a) distribution of   along the reinforcement bar at the end of the analysis (values at Gauss points); (b) evolution of   in time for the cracked element in the middle of the specimen.

Figure 10 :
Figure 10: Comparison of the evolution in time of   and   : the crack opening in concrete in the middle of the element is nearly equal to the sum of the absolute values of slip in the left and right node of the same element.