A simple method for non-linear analysis of steel fiber reinforced concrete

This paper proposes a physical non-linear formulation to deal with steel fiber reinforced concrete by the finite element method. The proposed formulation allows the consideration of short or long fibers placed arbitrarily inside a continuum domain (matrix). The most important feature of the formulation is that no additional degree of freedom is introduced in the pre-existent finite element numerical system to consider any distribution or quantity of fiber inclusions. In other words, the size of the system of equations used to solve a non-reinforced medium is the same as the one used to solve the reinforced counterpart. Another important characteristic of the formulation is the reduced work required by the user to introduce reinforcements, avoiding "rebar" elements, node by node geometrical definitions or even complex mesh generation. Bounded connection between long fibers and continuum is considered, for short fibers a simplified approach is proposed to consider splitting. Non-associative plasticity is adopted for the continuum and one dimensional plasticity is adopted to model fibers. Examples are presented in order to show the capabilities of the formulation.


Introduction
The importance of a good representation of fiber reinforced media for engineering analysis can be identified when observing the great amount of effort in studying the phenomenological behavior of this kind of material and various alternatives present in commercial softwares or scientific papers in order to solve this kind of problem (BARZEGAR; MADDIPUDI, 1994;GOMES;AWRUCH, 2001).For microscopic situations, sophisticated schemes have been recently developed to model non-homogeneous materials including fibers (FISH;BELSKY, 1995;LE TALLEC, 1994;HUET, 1990).However, for macroscopic engineering problems, as reinforced concrete the adoption of usual computational softwares 1,2 becomes a difficult task for engineers and applied scientists, due to the existent cumbersome schemes for fiber inclusion into the discrete model.
Acta Scientiarum.Technology Maringá, v. 32, n. 4, p. 367-374, 2010 Engineering finite element analysis of fiberreinforced domains, with arbitrary fiber distribution, for long or short fibers, is limited to two main approaches.The first is based on coincident nodal positions for fiber and continuum discretization.The second is based on special continuum elements that include fiber characteristics (rebar elements) following a simple straight pattern.The first procedure results in an exhaustive discretization of the continuum in order to adapt it to the fiber discretization, resulting in either a huge number of degrees of freedom and difficult mesh and re-mesh strategies or a simple mesh where the user should introduce node by node the reinforcement connection.The second alternative results in difficulties for the connection among elements and may generate truncated fibers (Figure 1).
Another possibility, for Finite Element Method (FEM) analysis, is to treat the reinforced medium as homogeneous, with equivalent properties (HYER, 1997;LEKHNITSKII, 1981).This approach leads to simplified stress distributions that may overestimate the strength of the analyzed solid or structure.This characteristic is worse when dealing with low rate fibers/continuum.Following this reasoning this paper proposes a finite element formulation that avoids the inconvenience of exhaustive discretization and truncate fiber discretization of FEM procedures.
The proposed formulation allows the consideration of short or long fibers, arbitrarily (or not) distributed inside the domain without increasing the number of degrees of freedom of the numerical model (Figure 2).It is based on kinematical considerations that write the displacement of fiber nodes as function of the displacement of continuum elements.In order to be complete, in section 2 the basic equations of FEM are presented, preparing the introduction of fibers.In section 3 the stiffness matrix including fibers is deduced and in section 4 plastic residuals are introduced in the final system.The strategy to allow the fiber meshing to be free from the matrix discretization is described by Vanalli (2004).A simple random fiber generation is specially applied in order to reproduce the fiber-reinforced specimen.Examples are shown, testing the formulation against theoretical solutions and experimental data.

Basic equations
In this section the basic equations necessary to build the proposed technique are briefly described, following the Principle of Virtual Work (PVW) together with the concept of non-linear residual stress.This approach is interesting because one can introduce plasticity or damage mechanics into the equations without changing the basic steps.Index notation will be followed.In order to achieve the PVW statement one may start from the equilibrium equation of an infinitesimal portion of the body, (1) where: b i represents body force and ij  stress.
Comma means partial derivative and repeated indices represent summation.Performing the inner product between equation (1) and the virtual displacement field i u  and integrating it over the analyzed body  results: Using the Gauss' theorem one achieves: where:  is the boundary and j  its outward normal.Applying the known properties , where ij  is the virtual strain, equation (3) turns into: This equation is the PVW for static problems.The first and the last integrals are, respectively, the work variation of surface and domain forces and the second represents the variation of internal energy.The internal energy can be conservative, resulting in an elastic analysis.For this situation one may adopt linear elasticity, given by the Hookes law ( ) turning equation ( 4) into: The total strain kl  is the composition of elastic and plastic parts, i.e.: The elastoplastic constitutive model is defined as: where: p ij  is the accumulated residual stress achieved by any well known elastoplastic procedure (SIMO;HUGHES, 2000).It is worth noting that the formulation described here does not use tangent matrix, but it can be implemented if desired.The residual stress can also be achieved following other non-linear behavior, as damage mechanics for example.Substituting equation ( 7) into equation ( 4), the desired non-linear statement is written as: The stiffness matrix The starting point of the proposed formulation is equation ( 9).Fibers will be considered bounded inside the domain and free of external forces, so the first and last integrals of equation ( 9) will not suffer its influence.The simplified model to consider splitting of short fibers is described in section 5.
The second integral represents the variation of elastic strain energy stored inside the wholly body, so it contains both matrix and fiber's influence, i.e., it can be split out into two parts: where: is the Young´s modulus for fibers and is the Hooke´s elastic tensor for matrix.
Fibers are approached by simple two-node truss elements, so the considered strain is only the longitudinal one and is written as a function of the displacement of the two nodes and the angle (  ) that the element forms with the horizontal (Cartesian) axis (x), as follows: where: u represents displacement in x direction and v represents displacements in y direction.Equation ( 11) is written in a more general way, as follows: where: ) f ( u  are the four degrees of freedom of the fiber element and relates nodal displacements to the longitudinal strain of the truss bar.The adopted solid element (2D case) is the well-known quadratic strain triangle (QST).From usual FEM assumptions similar expression, as equation ( 12), for displacement-strain relation can be written, as: displacements to strains.Substituting equations ( 12) and ( 13) into equation ( 10) results: Performing the indicated integrals over solid and fiber elements and changing representations to vector notation, one writes: Recalling that the fiber element is the two-node truss element, the displacement vector be written as where: is related to the second node node j.For any fiber element, using equations ( 15) and ( 16) one writes: At this point one should accept knowing in which element each fiber node is contained and its non-dimensional co-ordinate (Figure 2).The important kinematical consideration is that the displacements of any node of the fiber element can be written as a function of the displacement of the nodes of the 2D element that contains it, as follows: where:  is a matrix that contains the shape functions of QST element (i or j) calculated for the node location of the reinforcement.The null matrix, 2 x 20 0 , is an auxiliary value introduced to allow each node of any reinforcing element to belong to different (or the same) 2D elements.From equations ( 17), ( 18) and ( 19) one writes: This expression shows that the variation of fiber strain energy can be written as a function of continuum elements that contains it.The upper bar means that the fiber characteristics are written following exactly the degrees of freedom of the continuum.
One should observe that an adequate degree of freedom numbering must be respected in order to build the global stiffness matrix.This numbering is implicit for indices i and j present over matrix 40 x 40

K
, which means lines and rows belonging to different (or same) 2D finite elements.From equation ( 21) and the former reasoning, equation ( 15) is written as Equation ( 22) results into a general reinforced medium strain energy variation and, as a consequence, the resulting stiffness matrix can model as many as desired reinforcing fibers using exclusively the original continuum degrees of freedom.

The plastic residual force
In the previous section the stiffness matrix for reinforced media was developed, in this section a similar procedure is followed to determine the residual force, derived from residual stress, applied only in the positions due to the matrix media.From equation ( 9) one writes the variation of the work of dissipative internal forces as: This equation can be split out into two parts: Taking into account that residual stress and virtual strain are constant along the truss element and using the results of equations ( 12) and ( 13), equation ( 24) is rewritten as: where A is the area of the cross section of the fiber.Equation ( 25) can be written following the vector notation of equation ( 15) and using relations ( 18) and ( 19) as: Equation ( 26) can be written considering all finite elements simply as: where upper bar means that the fiber residual forces are written following exactly the degrees of freedom of the continuum.Introducing equations ( 22), ( 27) into equation ( 9) and considering, for simplicity, only concentrated external forces, the PVW turns into: Considering that the virtual displacement of nodes are arbitrary, equation ( 28) turns into the non-linear system of equation to be solved, i.e.:

 
The non-linear aspect of equation ( 28b) falls into the achievement of plastic residuals, as elastic stiffness matrices are adopted.The achievement of the residual plastic stress follows any wellestablished procedure (SIMO;HUGHES, 2000).
Here we follow the non-associative rule for Tsai-Wu yielding surface presented by VANALLI, 2004.It is important to note that to write equation ( 28) and ( 29) the nodes of fibers are considered free from continuum discretization, but its position regarding solid elements were considered known.In order to guarantee the necessary comfort for users, section 6 presents the way the computational code automatically identifies fiber nodes positions.

Reinforced concrete beam
In this example a reinforced concrete beam is analyzed using the developed formulation.The concrete is considered elastoplastic following two well-known criteria, Drucker-Prager and Tsai-Wu (adapted to run isotropic material).The reinforcement follows one-axial bilinear elastoplastic relation as it is modeled by simple truss elements.The results are compared with experimental values given by Takeya (1972) and with a numerical model based on BEM -FEM coupling given by Coda (2001).
The geometry and the reinforcement arrangement are shown in Figure 3 Symmetry is considered in order to model the problem.A half of the beam is discretized into a homogeneous mesh of 20 x 20 QST finite elements (plane stress) constituting a total of 1992 degrees of freedom.Thirty two truss finite elements, located exactly in the position of the reinforcement, are employed to model the reinforcement.As previously mentioned, no sliding between long reinforcements and concrete is allowed.In Figure 4, the displacement at the center of the beam is plotted against the applied force.Displacement control has been adopted.As one can see the behavior of the proposed formulation is very close to the average of the experimental result and compares well with BEM -FEM formulation.The smooth behavior of the BEM/FEM technique is explained by the poor discretization used by that author to model concrete along vertical direction.Additional information is given in what follows using a mesh of 130 x 30 finite elements.In Figure 5 one can see the behavior of shear forces along the interface of concrete and reinforcement for an elastic situation, considering and not considering stirrups.As one can see the transfer of forces between fiber and matrix is slightly faster in the situation where stirrups are considered.Finally, the Figure 6 shows the evolution of concrete degeneration (in terms of stress) for different levels of imposed displacement.It is possible to verify that at an imposed displacement of 10 mm the beam loses its supporting capacity.

Reinforced concrete beam --experimental analysis
This example compares experimental results from Ashour et al. (2000) with numerical results obtained using the developed formulation.The tested beams are reinforced by long and short fibers.Long fibers are called tensile reinforcement and short fibers are called simply as fibers.A comparison between experimental and numerical results was made for 3 of the 27 tested beams.It is interesting to note that the fibers are randomly spread over the domain in order to reproduce the laboratory conditions.As the numerical model is twodimensional, the total amount of fibers to spread over the domain is 70% of the total used to model the three-dimensional specimen.
As in the experimental analysis, concrete compressive strength of 49 MPa and tensile reinforcement ratios of 1.18% were used and the fiber contents were 0.0, 0.5 and 1.0% by volume.The geometry and reinforcement arrangement are shown in Figures 7 and 8 extracted from Ashour et al. (2000).As can be seen in Figure 8 flexural reinforcements, 2  18 mm and 2  6 mm totaling a Displacements (cm) Force P (kN) tensile reinforcement ratio of 1.18% are used.Perfect plasticity behavior is adopted for this material.Hooked-ends mild carbon steel fibers (short) with average length of 60 mm, nominal diameter of 0.8 mm and yield strength of 1100 MPa are employed in the experimental specimen (ASHOUR et al., 2000).A concrete with compressive strength of 49 MPa composed the matrix media.Again, perfect plasticity behavior was adopted for concrete and softening behavior for fibers reinforcements.Table 1 presents general material data for reinforcement.The adopted effective yielding stress for short fibers is 276.75 MPa.The difference between the analyzed beams is the fiber contents (percentage of fibers) (Table 2): Obs.: fibers area: 0.50265 mm 2 ; Beam volume: 170,000,000 mm 3 .
The numerical results, in Figures 9, 10 and 11, extracted are compared with experimental ones for 0.0, 0.5 and 1.0% of fiber contents, respectively.It is possible to verify the stiffening in the numerical results compared to the experimental ones.These differences between the results can be associated to factors not considered in the numerical formulation as, for example, the slipping of long fibers.

Conclusion
A simple finite element formulation based on kinematics considerations is presented in order to model fibers arbitrarily distributed inside homogeneous media.A two-dimensional implementation has been successfully carried out using QST finite elements to model the continuum and constant strain finite element to model the one dimensional reinforcement (truss).Very good results have been achieved for linear and physical non-linear analysis concerning usual civil engineering materials.The simple splitting of short fibers presented a good behavior and can be easily evaluated.The formulation is promising and further improvements, such as sliding between long fibers and matrix as well as a three dimensional implementation are recommended.

Figure 3 .
Figure 3. Geometry and reinforcement distribution, lengths in centimeters and diameters in millimeters.

Figure 6 .
Figure 6.Evolution of concrete degeneration for prescribed displacements.
. The material properties are the same as the one measured in laboratory, i.e., for steel the Young modulus is E s 

Table 1 .
General material data.
f' c : Compressive strength of concrete (at 28 days); f r : Modulus of rupture of concrete; f' sp : Splitting tensile strength of concrete.

Table 2 .
Characteristics and quantities of steel fibers used in each beam.