Continuous Modeling Technique of Fiber Pullout from a Cement Matrix with Different Interface Mechanical Properties Using Finite Element Program 1

Fiber-matrix interface performance has a great influence on the mechanical properties of fiber reinforced composite. This influence is mainly presented during fiber pullout from the matrix. As fiber pullout process consists of fiber debonding stage and pullout stage which involve complex contact problem, numerical modeling is a best way to investigate the interface influence. Although many numerical research works have been conducted, practical and effective technique suitable for continuous modeling of fiber pullout process is still scarce. The reason is in that numerical divergence frequently happens, leading to the modeling interruption. By interacting the popular finite element program ANSYS with the MATLAB, we proposed continuous modeling technique and realized modeling of fiber pullout from cement matrix with desired interface mechanical performance. For debonding process, we used interface elements with cohesive surface traction and exponential failure behavior. For pullout process, we switched interface elements to spring elements with variable stiffness, which is related to the interface shear stress as a function of the interface slip displacement. For both processes, the results obtained are very good in comparison with other numerical or analytical models and experimental tests. We suggest using the present technique to model toughening achieved by randomly distributed fibers.


Latin American Journal of Solids and
can significantly increase the material fracture strength.We can rate the main advantages of fiber reinforcement composites with the brittle matrix as (1) improvement in fracture toughness and tensile strength of the composite and (2) inhibition of crack-open and propagation.These desirable features have been obtained with a low fiber volume fraction, generally less than 2% (Fantilli and Vallini, 2008).
The fiber-matrix interface property directly affects the mechanical behavior of the composite.The most common way to observe the interface effect is through fiber pullout tests, which consists of debonding process and followed pullout process (Figure 1).For these two processes, the adhesion between fiber and matrix is the key (Kim and Mai, 1998).Engineers and researchers always seek for new techniques of characterizing interface and numerical methods to determine the strengths of the composite reinforced by the fibers with the characterized interface in order to obtain better composites with expected mechanical performance.For this, a number of numerical and computational models has been developed.Without a doubt, the finite element method is one of the best computational tools to solve this problem.With a focus on local mixed-mode at the interface, Becker and Lauke (1997) developed a finite element model using a fracture mechanics criterion for the interface debonding process of a single fiber pullout test.Also using this kind of test, Liu et al. (1999) investigated the effects of fiber pullout rate, thermal residual stress, friction coefficient and fiber volume fraction.In their work, the local shear strain criterion was adopted for interface debonding.Recently, Wei et al. (2012) implemented a cohesive damage model in finite element modeling to simulate fiber-matrix interface debonding in order to investigate the interface shear stress distribution and effects of shear stress transfer across the debonding interface.Pochiraju et al. (2001), Zhang X. et al. (2004) and Tsai et al. (2005) simulated fibers pullout from an epoxy matrix by the finite element method using interface and contact elements, the cohesive zone model and Coulomb friction.Rodrigues et al. (2015) used interface finite elements with high aspect ratio and modeled the interaction between the steel bars and concrete.Structures 13 (2016) 1937-1953 Since many parameters are involved in the pullout process such as the mechanical properties of fiber, matrix, and interface, as well as the geometry of the fiber (Wang and Friedrich, 2013), it is still difficult to continuously model a single fiber pullout test.Others factors such as mesh and contact problems increase the difficulty.Although the modeling of fiber pullout problem has been studied for almost twenty years, we note that the reports on the research advance are decreased substantially recently.The reason would rather be attributed to the great difficulty in the numerical convergence of modeling than that the problem has been resolved.Also, we note that there is scarcely the numerical modeling that can systematically deal with whole fiber pullout process with desired interface mechanical properties.No one of the reported works can continuously model whole fiber pullout with desired interface mechanical properties without interrupting and modifying the modeling program.

Latin American Journal of Solids and
This work presents continuously modeling technique without interruption to simulate the whole process of fiber pullout from cement matrix using popular finite elements program ANSYS® interacting with the MatLab®.At debonding stage, we suggest using interface elements with the cohesive zone model (CZM) presented by Xu and Needleman (1994) and obtain a satisfied relationship between the tractions acting on the fiber surface and the displacements jump across the interface.At fiber pullout stage, we propose employing nonlinear spring elements to model the debonded interface.The nonlinear spring elements are of variable stiffness and able to relate interface shear behavior with fiber slip displacement relative to the matrix.The number of the nonlinear spring elements on the interface decreases with the fiber slip displacement.The switching from debonding process to pullout process is made automatically through a computational subroutine implemented in the software MatLab®.

PROPOSED FINITE ELEMENT MODEL FOR FIBER PULL-OUT PROCESS
The axisymmetric model representing a single fiber inserted in a concentric cylindrical matrix is shown in Figure 2. By use of the software ANSYS, all geometric models of fiber and matrix used in this work were discretized with 2700 elements named as PLANE182 by ANSYS.The element PLANE182 is defined by four nodes and each node has two degrees of freedom: translations in the nodal x and y directions.The element has plasticity, stress stiffening, large deflection, and large strain capabilities.To simulate different processes, two interface modes were used (Figure 2): interface elements for debonding process and spring elements for pullout process.For the interface element, the nodes on both sides (fiber and matrix) were duplicated.300 asymmetric elements were used to represent the cohesive interface, each with 4 nodes and zero thickness.The interface element INTER202 of ANSYS has four nodes and two degrees of freedom at each node.The spring element COMBIN40 of ANSYS with nonlinear stiffness k is made up of a pair of nodes.One of them is on the fiber and another in the matrix.In total 150 spring elements were used.The COMBIN40 element has two nodes and is a combination of a spring-slider and damper in parallel, coupled to a gap in series.A mass may be associated with one or both nodal points.The element has one degree of freedom at each node, either a nodal translation, rotation, pressure or temperature.The debonding process is uniquely dependent on the required properties of the cohesive interface; meanwhile, the pullout process is dependent on the spring stiffness.When considering shear stress as constant, known as constant- model (Marshall et al. 1985, Gopalaratman and Shah 1987, Li 1992), we can write it as

Latin
or where Pmax is the maximal pullout load before the debonding process is over, d the fiber diameter and Le the embedded length.If  is constant, then the right term in Equation 2 is also constant.
We denote the initial stiffness of the spring at the end of debonding process as which is corresponding the stiffness of the null slip displacement 0 s  .When 0 s  we suggest the spring stiffness as: where n(s) is the number of effective springs.Using the spring stiffness determined in previous load step, we compute current slip displacement s and know the position of the fiber end node or how much nodes on the matrix side has been exceeded by the fiber end node, then, the number of effective springs can be accounted as the total nodes number on the matrix side minus the nodes number exceeded by the fiber end node (Figure 3).In this way, we can get the springs stiffness, ki, for the next load step through Equation 4.However, the constant  model is not suitable to describe the fiber pullout, since there is ample evidence which shows that the shear stress varies with the interface slip displacement and causes slip-hardening or slip-softening (Figure 4) as reported by Li and Stang (1997), Lin and Li (1997) and Wang et al. (1988).The variation may be attributed to several factors such as fiber surface abrasion, interface roughness, matrix plasticity, fragmentation of fibers coatings at the interface, etc. Latin American Journal of Solids and Structures 13 (2016) 1937-1953 In these cases, the spring stiffness is also a function of the variation of the shear stress and thus, we propose the spring stiffness as Using this relation, we can simulate the pullout of fibers in different types, for example, steel fiber with slip softening behavior, as seen in the experimental works of Naaman and Shah (1976) and Li et al. (1991) and synthetic fibers with slip hardening behavior commonly observed in this fiber type, where the interfacial stress increases with slip distance (Wang et al. (1988), Li et al. (1995) and Ting et al. (2015)).This hardening behavior is due to wear between the surfaces, occurring at the moment when low hardness fibers try to slide over the stiffer matrix, a resistance to movement is imposed by the accumulation of the fiber debris, and thus increases the interface shear stress.These characteristics bring us a comprehensive model for simulating a great number of cases and helping to evaluate others factors such as fiber bridging stress and toughening.

COMPUTER ANALYSIS PROCEDURE
The MATLAB is a powerful program that offers not only different optimization algorithms already implemented but also the possibility of conducting an optimization algorithm defined by the user.Thus, the interactive use of the MATLAB and the ANSYS allows the user to fully control and optimize the simulation in order to get the best results as soon as possible.
Figure 5 shows the computational procedure made in this work.We performed simulations using the commercial software ANSYS® 11.0.With the aid of MATLAB® R2012a software, we shortened the simulation time and facilitated switching between debonding and pullout processes and taking the loop of the displacement δ applied at the right end of the fiber (Figure 3).
Following the procedure, first, we start the subroutine pullout.m in Matlab and create the text file "debonding.txt" by introducing an initial small displacement δ about 0.01 mm.The txt file contains the commands in ANSYS Parametric Design Language (APDL), which generate the geometric model, mesh and displacement restrictions of the model and finite element solution in batch mode.The displacement restrictions are imposed on the boundary of matrix cylinder.
As a result, the ANSYS yields the file "Results.txt" in two columns: the first for the applied displacement δ and the second for the reaction force obtained as the sum of the reactions at the displacement prescript nodes.If debonding is over, the reaction force will drop to zero.Therefore, after the computation of current load is accomplished, we check whether the reactions sum is zero in order to verify whether debonding has been finished.
If not yet, we increase the displacement by an amount same as the initial displacement applied.The new results are stored as a new line in the "Results.txt"file.
If the condition of zero reactions sum is met, i.e. the debonding has been completed, the subroutine loads the new file "Pullout.txt",and then the interface elements are switched to spring elements and the characterized behavior of shear stress on the interface is introduced.Start again with the displacement δ applied and the initial stiffness ki both obtained immediately before the debonding completion and calculate the displacement of the left end of the fiber, which is sliding relatively to the matrix i.e. the slip displacement s.The sliding quantity determines the number of effective springs n(s) as shown in Equations 4 or 5 and alters the spring stiffness, depended on what shear stress behavior is expected.Following this, a new larger displacement δ is added.The modeling is kept on until no effective spring is available, which means that the pullout process is over.

Debonding Process
The cohesive zone model consists of a constitutive relationship between tractions T acting on the interface and the corresponding separation distances i.e. displacement jump Δ across the interface (Lin et al. 2001).The use of cohesive zone allows a better control over the debonding process.As soon as the final separation is accomplished, the pullout process will start.
In the present work, the fiber-matrix interface was based on the exponential failure model of Xu and Needleman (1994)  interface parameters used, are listed in Table 1.Since there is little information on the tested specimens, we assume that the matrix is 1.43 mm in diameter and 12 mm in length, approximately corresponding to fiber volume fraction 2%.The matrix is characterized by the stress-strain curve in compression, as shown in Figure 6, with the maximum compressive strength equal fc.We suggested a curve and implanted it into the ANSYS in order for the modeling to be more realistic.The obtained results are compared with the experimental test of Li et al. (1991) and the two numerical models of Mobasher and Li (1995) and Li and Mobasher (1998) (Figure 7).These two numerical models dealt the interface in different manners.The first developed a theoretical model of interface failure based on the criterion of maximum strain energy release rate applied to a partially debonding length; meanwhile, the second dealt the interface as a third linear elastic phase of the composite, where the debonding criterion is given by integrating maximal normal and shear stresses into the finite element method.
From Figure 7 we see that the proposed model produces a good approximation to the experimental results, especially when the slip is large.Regardless of there is still a bit of deviation from the experiment curve on the initial debonding stage, in comparison with the other two numerical models, we get a better curve of pullout force -slip displacement, because the curves resulted from the other numerical models can not extend to a sufficient slip displacement and the debonding process is ended too early.Furthermore, our pullout force peak appears at an instant closest to the experimental one.The maximal peak load obtained in the debonding process is a very important parameter for the pullout modeling, as shown in Equation 3, because it determines the initial spring stiffness that would compromise the results.To check whether the slip displacement step size affects the pullout force, we repeated the computation and plot Figure 8 bellow with slip displacement in different step sizes.From Figure 8 we see that the slip displacement step size doesn't affect the pullout force.This coincides with our intuition that the slip displacement step size shouldn't influence the pullout force.It has been verified that implanting the curve of stress-strain for concrete in compression into the ANSYS is indispensable.If the matrix phase is considered only as linear elastic, the results will be overestimated.Figure 9 shows the influence of three different concrete compressive curves on the pullout load.As observed, although the compressive curves are different, there is the unnoticeable difference in the responses obtained.It is not yet clear why there is little difference among the responses of the three compressive curves.Perhaps, the response is mainly depended on the interface performance parameters max  , n  (normal separation characteristic length) and t  (tangential sep- aration characteristic length) as well as the matrix stiffness.However, there is an obvious difference between the response of linear elastic curve and the ones of the real compressive curves, which have the softening characteristic that inhabits continue raise of pullout load after the peak load reaches Pmax.The curve 2 with a lower residual strength after the peak stress presents a slightly lower response than the others do. Figure 10 shows the influence of the interface strength ( max  ) during the debonding process.We note there is a significant increase in the maximum pullout force with increasing max  .
The gain in the pullout force is more than 100% when max  varies from 1.5 MPa to 3.5 MPa.
The proper adhesion between fiber and matrix helps to enhance the composite strength.However if the adhesion is too strong to allow the fiber debonding, the stress in the fiber would exceed the fiber tensile strength and lead to the fiber breakage, decreasing the number of effective fibers, consequently, degrading contribution to improvement in the mechanical properties of the composite especially in the toughness, as demonstrated as some experiments reported by Li et al. (2001).As the debonding predominantly occurred in the tangential direction, we also investigated the influence of tangential separation characteristic length t  on the pullout force.The curves of the pullout force versus slip displacement show an inverse relationship between the maximal pullout force and the interface characteristic length (Figure 11).This relationship may be attributed to the cohesive law of Xu and Needlemann (1994) used by the ANSYS, which imposes the conservation of the works realized by normal and tangential separations, resulting in the increase of shear stress and pullout force when the interface characteristic length decreases.

Slip Softening Modeling
As seen, the debonding process is represented very well by exponential failure model.Since the experimental test of Li et al. (1991) did not show the pullout process, thus, we validated the pullout model (Equation 4) by means of the experimental test of Leung and Shapiro (1999) who investigated the influence of steel fibers with different yield strengths ( y f ).To simulate the slip softening behavior commonly observed in the composite reinforced by steel fibers, an appropriable interface shear stress model (s)  should be introduced into the Equation 4.
In this work, after debonding process is accomplished, we adopted the relationship (s)  suggest- ed by Fantilli and Vanilli (2003) (Figure 12a) as: where, max  is the maximum shear stress (Equation 1); s the slip displacement, which is approxi- mately equal to  the applied displacement; 1 s the applied displacement just immediately before debonding process was completed; fin  is the asymptotic value of shear stress and kc a coefficient.These latter two parameters given in Figure 12b are dependent on fiber type and its manufacture process.There is the difference between the present work and of Fantilli and Vanilli in the determination of parameters 1 s and max  .Fantilli and Vanilli (2003) estimated 1 s and defined max  as function of fiber diameter and the concrete compressive strength.We got 1 s and max  in a more natural manner through simulating debonding process and identifying when the debonding process was over and the respectively applied displacement and maximal pullout force.It is clear that all these parameters rely on the CZM and the results obtained from the modeling of the debonding process, which is fundamental to complete the following modeling of fiber pullout process.Fantilli and Vallini (2003).Leung and Shapiro (1999) tested five specimen groups.For convenience, we simulated the case of the fiber type b.All parameters used are listed in Table 2, where t f and c f are the tensile and compressive strengths of the matrix respectively.The required stress-strain relationship for concrete in compression is the same as in Figure 6, used for the debonding process validation.The composite has 0.05% of fiber volume fraction, which is defined here by the ratio of the volume of a single fiber to the matrix volume mentioned in section 4.1.Leung and Shapiro (1999) where the elastic modulus was 30 GPa; b Estimate; c By Li and Mobasher (1998).
Table 2: Mechanical and geometrical properties of the specimen by Leung and Shapiro (1999).
The comparison of results from the pullout process modeling is demonstrated in Figure 13.It is clear that the constant- model causes large discrepancies and the model of considering the interface shear as a function of the slip displacement, associated to maximum pullout load obtained in the debonding process, is more suitable for modeling pullout process.The model proposed in this work presents a good agreement with experimental result, especially when the slip displacement is larger than 2 mm.Therefore, the use of spring elements for simulating interface is able to control interface performance.
In comparison, our result and the analytical result of Fantilli and Vanilli (2007) are similar.With kc=2 as used by Fantilli and Vanilli, our result prevails against the one of Fantilli and Vanilli only when the slip displacement is larger than 2 mm.However with kc=5, our result is closer to the experiment curve over the whole pullout process.Fantilli and Vallini (2007) and the experiment made by Leung and Shapiro (1999).Solids and Structures 13 (2016) 1937-1953

Slip Hardening Modeling
To exemplify interface slip hardening behavior we referred the experimental test of Li et al. (1995) where polyethylene fibers were added into a cement matrix.The geometric and mechanical properties of fiber, cement matrix, and interface parameters are listed in Table 3.The composite has 2% of fiber volume fraction.For the comparison, the same experimental test was referred to the analytical model of Lin and Li (1997) with a simple interface constitutive relationship adopted from Bao and Song (1993) to quantify the slip-hardening interface behavior, which is: where the df is fiber diameter, 0  the frictional sliding shear stress at the tip of debonded zone where no slip occurs and  a dimensionless hardening parameter.The latter two parameters need to be determined empirically.
In terms of Equation 7, Lin and Li (1997) derived the pullout load: where / , cosh 1 This work adopted the slip hardening model created by Wang et al. (1998).The relationship between interface shear stress  and slip displacement s in Equation 4 was taken as a quadratic function of s: where the constants a0, a1 and a2 are empirically determined, so that the theoretical curves remain reasonably close to the experimental curves.We have chosen three load values to fit the experimental curve, they are P1 (the load when the debonding was completed), P2 (the maximal load of the experiment) and P3=0 when the fiber was pulled out from the matrix.
Figure 14 gives the comparison of the results.Clearly, the proposed model produces an alike but better fiber pullout curve than the curves of Lin and Li (1997), who used an optimal parameter 0.0125

 
, before the load hits its maximum.After the load attains the maximum, all the modeled curves deviate from the experimental curve, since the experimental curve is not quadratic.As seen, the  parameter is very sensible.To adequate value, many trials should be done.This is an uneasy matter. with a value larger than 0.015 or less than 0.005 leads larger deviation from the experi- ment.With 0.0125

 
, you get a good approximation to the experiment before the peak load but large deviation after then; with 0.0085

 
, you have the inverse result.In comparison, our model has no need of trials but requires three load values.The first value is available automatically when the debonding is complete; the second value is fit to the maximal experimental one, and the third one is P=0 when the full fiber is pulled out from the matrix.Among the three values, only the second one is obligated to the experiment, the others are acquired naturally.With these three values, we got a better approximation to the experiment.Therefore, we can say that the presented model integrated into the ANSYS is suitable for the modeling because it is more simple but capable of describing very well the fiber pullout process.and the analytic model of Lin and Li (1997) with different  values.

REMARKS
As satisfactory results were obtained in the modeling, the following techniques are suggested to realize continuous modeling of fiber pullout process: • For the debonding process, the employment of interface element with cohesive surface traction and exponential failure behavior is highly recommended.
Latin American Journal of Solids and Structures 13 (2016) 1937-1953 • The compression test curve of the matrix material with a residual strength after peak load must be integrated into the modeling because linear elasticity behavior brings about super stiffness and leads modeling to fail.• For pullout process, the use of spring elements with variable stiffness, which is associated to interface shear stress as a function of interface slip displacement through Equation 4 or 5 proposed in this work, prevails over contact elements and allows for modeling different interface shear behavior such as slip softening or hardening and constant shear.Though we have realized the modeling by means of the ANSYS, the points given above may be included in any finite element programming for continuous modeling of fiber pullout process.The suggested computational procedure with the interaction of the ANSYS and the MATLAB can reduce substantially the time consumption in computation, which is of great advantage when there are a large number of parameters involved.With the presented technique, we are able to model toughening achieved by randomly distributed fibers.

Figure 2 :
Figure 2: Axisymmetric concentric-cylinder model, mesh, and modes of the interface.
and the modeling was accomplished through commercial finite element software ANSYS 11.0 where for interface separation, only three parameters are needed: (1) maximum normal stress max  at the interface; (2) normal characteristic length n  and (3) tangential charac- teristic length t  .The experimental test of Li et al. (1991) is generally used for the validation of the cohesive zone model of the debonding process.The authors Li et al. paid great attention to the process and their work has become a typical reference for the study on the behavior of steel fibers inserted into cement matrix.The geometric and mechanical properties of steel fiber and cement matrix, as well as Latin American Journal of Solids and Structures 13 (2016) 1937-1953

Figure 6 :
Figure 6: Stress-strain curve in compression suggested for the concrete.

Figure 7 :
Figure 7: Comparison of pullout force vs slip displacement of present work, Mobasher and Li (1995), and Li and Mobasher (1998) with the experimental test of Li et al. (1991).

Figure 8 :
Figure 8: Influence of slip displacement step size Δδ on pullout force.

Figure 9 :
Figure 9: Influence of stress-strain compression curve of concrete on the pullout force.(a) Different concrete compressive curves.(b) Responses of pullout force.

Figure 10 :
Figure 10: Influence of the interface strength max 

Figure 11 :
Figure 11: Influence of the tangential separation characteristic length t  , considering

Figure 13 :
Figure 13: Comparison of the results from the proposed models, the analytical solution of

Figure 14 :
Figure 14: Comparison of results of the proposed model with the experimental test ofLi et al. (1995)

Table 1 :
Mechanical and geometrical properties of the specimen considered.

Table 3 :
Mechanical and geometrical properties of the specimen considered.