. This paper presents a computational rational model to predict the ultimate and optimized load capacity of reinforced concrete (RC) beams strengthened by a combination of longitudinal and transverse fiber reinforced polymer (FRP) composite plates/sheets (fl exure and shear strengthening system). Several experimental and analytical studies on the confinement effect and failure mechanisms of fiber reinforced polymer (FRP) wrapped columns have been conducted over recent years. Although typical axial members are large -scale square/ rectangular reinforced concrete (RC) columns in practice, the majority of such studies have concentrated on the behavior of small -scale circular concrete specimens. A high performance concrete, known as polymer concrete, made up of natu ral aggregates and an orthophthalic polyester binder, reinforced with non -metallic bars (glass reinforced polymer) has been studied. The material is described at micro and macro level, presenting the key physical and mechanical properties using different e xperimental techniques. Furthermore, a full description of non -metallic bars is presented to evaluate its structural expectancies, embedded in the polymer concrete matrix. In this paper, the mechanism of mechanical interaction of smooth and lugged FRP rods with concrete is presented. A general modeling and application of various elements are demonstrated. The contact parameters are defined and the procedures of calculation and evaluation of contact parameters a re introduced. The method of calibration of the calculated parameters is presented. Finally, the numerical results are obtained for different bond parameters which show a good agreement with experimental results reported in literature. Keywords : fiber-reinforced polymers; reinforced concrete (FRP); mechanical interaction; bond; contact; friction Khosro Shahpoori Arani et al. electromagnetically corrosion (Sinaei et al. 2011). The acceptance of FRP reinforcing bars in structural engineering has been inhibited partly due to lack of design criteria, particularly with regard to bond (Feroldi and Russo 2016) . Several techniques have been developed to inhibit or eliminate corrosion, including epoxy-coated reinforcing bars, synthetic membranes, latex concrete, Cathodic protection, special paints, and Sealants (Benmokrane and Tighiouart 1996). The FRP material is transversely anisotropic. As a result, their mechanical properties in longitudinal direction differ from the other two orthogonal directions. Many experimental studies have demonstrated that the key properties of the physical and mechanical behavior of FRP rebars are qualitively and quantitively different from the well-known properties of steel rebars (Cosenza et al. 1997) . In recent years many experimental studies have been conducted in order to investigate bond between FRP rebars and concrete. However, fewer researches have been carried out by scientists on modelling and simulation of bond of FRP rebars to concrete. The aim of this paper is to propose a reliable FE modelling of mechanical interaction of smooth rods with the surrounding concrete by evaluation, calibration and validation of contact parameters. Retrofit of reinforced concrete (RC) members using externally bonded fibre reinforced polymer (FRP) components has become a very common practice, widely recognized by modern design codes (Naderpour et al. 2011, Sinaei et al. 2011, Ashour and Kara 2014, Mohebi et al. 2016, Naderpour et al. 2016). In particular, the flexural strength of a beam can be significantly increased by application of carbon (CFRP) or glass (GFRP) FRP plates/sheets adhesively bonded to the beam tension face. This paper focuses on efficient and accurate finite element (FE) modelling of RC beams flexurally retrofitted with externally bonded FRP plates/sheets. The finite element modeling has been followed those of previously conducted researches (Kheyroddin 2008b, Arabnejad Khanouki et al. 2010, Arabnejad Khanouki et al. 2011, Daie et al. 2011, Jalali et al. 2012, Aghakhani et al. 2015, Bazzaz et al. 2015, Khanouki et al. 2016, Mansouri et al. 2016, Safa et al. 2016, Toghroli et al. 2016, Khorami et al. 2017a, Mansouri et al. 2017, Bazzaz 2018, Heydari and Shariati 2018, Hosseinpour et al. 2018, Ismail et al. 2018, Mehrmashhadi et al. 2018a, Mehrmashhadi et al. 2018b, Sedghi et al. 2018, Toghroli et al. 2018a, Toghroli et al. 2018b, Zhao et al. 2018). FRP-retrofit of RC beams for flexural strengthening presents numerous advantages compared to other flexural strengthening techniques, e.g., strengthening using steel plates. Some of these advantages are: small increase in structural size and weight, easy transportation, speed and simplicity of in-situ application, and good resistance to corrosion and other degradation processes due to hostile environmental conditions. This technique has found numerous applications in retrofitting of bridges and buildings during the last two decades. At the same time, extensive experimental, analytical and numerical research has been carried out to understand and model the structural behavior of FRP-strengthened RC beams. For literature reviews on different aspects of FRP-strengthening of RC structures, the interested reader is referred to (Jalal et al. 2013, Nazerian and Sadeghiipanah 2013, Ashour and Kara 2014, Mohebi et al. 2016). Particular attention has been given to recognizing and understanding the failure modes that RC beams retrofitted with FRP can experience. The experimentally identified failure modes can be grouped as follows: (1) flexural failure by concrete crushing or by steel yielding followed by concrete crushing (flexure failure mode, which is similar to the failure mode of conventional/non-retrofitted RC beams), (2) flexural failure due to FRP rupture (FRP rupture failure mode), (3) flexural failure due to plate end interfacial debonding to concrete cover separation or to intermediate crack induced debonding (debonding failure mode), and (4) shear failure. A more detailed description of these failure modes can be found in (Rahimi and Hutchinson 2001, Naderpour et al. 2016). It is noteworthy that all the failure modes typical of FRP-retrofitted beams (i.e., FRP rupture and debonding) are brittle in nature. In addition, the debonding failure modes correspond to a less than optimal use of the strength capabilities of the FRP material. In order to increase the efficiency of the FRP-retrofit reducing the impact of debonding failure modes, mechanical anchorage techniques have been devised and employed, e.g., use of mechanical devices at the FRP plate/sheet ends, FRP sheets wrapped around the RC member at the FRP plate/sheet ends (U- wrap), and U-shaped FRP plates/sheets along the entire length of the RC beam (U-shape ) (Hadi 2007). In the existing literature, numerous studies deal with modelling of FRP-strengthened RC members for bending behaviour and correlation of experimental results with numerically predicted ultimate flexural strength, e.g., (Hadi 2007). Understanding and modelling of debonding failure modes has been and still is a very active field of research, mainly due to the complexity of the problem at hand. Indeed, there is a strong need of reliable and robust formulations of specialized FEs that could help researchers to model accurately FRP strengthened RC members and structures in order to predict their response under different strengthening configurations, e.g., for design-oriented parametric studies. Previous FE studies of FRP-retrofitted beams involve the use of refined FE meshes of (1) frame elements with an increased number of DOFs per element when compared to a classical Euler-Bernoulli frame element (Triantafillou and Plevris 1992), (2) two-dimensional plate/shell elements (Jerome and Ross 1997), or (3) three-dimensional solid elements (Yang et al. 2009). The high computational cost of structural response analyses based on FE models such as the ones referred above has prompted the development of purely numerical methods (i.e., not based on mechanics) for the analysis and design of FRP-strengthened RC structures (Yang et al. 2009). The research work presented in this paper develops a new nonlinear frame FE, based on the classical Euler-Bernoulli assumptions, and able to model the mechanical behaviour of FRP-strengthened RC beams. This FE, referred to as FRP- FB-beam in the sequel, (1) allows reducing the complexity and computational cost of FE analyses based on existing FE models, (2) provides a 62 Computational optimized finite element modelling of mechanical interaction of concrete with fiber reinforced polymer qh Fig. 1 The tangential and normal behavior of contact element sound mechanical description and interpretation of the phenomena leading to failure of FRP-retrofitted RC beams, and (3) simulates the structural response of the considered structural systems with accuracy satisfactory for practical applications. 2. Mechanical interaction of FRP smooth and lugged rods with concrete Bond strength is more affected by the low transverse stiffness which is dependent mostly on the matrix material. The bond mechanism of FRP reinforcement with concrete consists of two different aspects; the major and minor phenomena. The major phenomena are: friction and mechanical interferences, including bearing between lugs and the blocks of concrete. The former is the main cause of mechanical interaction in smooth rods while the latter mainly happens in the lugged rods. The minor phenomena are chemical adhesion, temperature and moisture effects and cog wheel behavior between the rib and core of reinforcement after shear cut-off takes place. The commonly used bond strength calculated by the division of the pull out force by the area of embedment length decreases as the embedment length increases for a specified db which is the diameter of bar (Cosenza et al. 1997) . The tolerance mechanism of pull-out load in the lugged rods is mainly based on bearing between the ribs and blocks of concrete. In fact, it is the reason of occurrence of general failure in system, such as “concrete splitting ” or “rib cut- off”. If concrete splitting happens in all concrete blocks below the ribs the main factor to provide further bond strength is the friction between different parts of concrete on both sides of the gap. However, in the case of lug shear cut-off the main factor, which provides bond strength, is the so-called “cog wheel ” behavior between the cut lugs and the internal core of reinforcement. It must be noted that decreasing the deformation depth results in some decrease in bond strength. The analytical results show that in the case of large bond lengths (larger than about 15 times the diameter), the analytical bond splitting strength is inversely proportional to the bond length and is clearly influenced by the Young modulus (Sakai et al. 1999). Bakis developed a simplified model of the mechanical interaction without explicitly modelling the geometry of the surface structure (Guo and Cox 2000) . 3. Numerical modelling In order to simulate the mechanical interaction of smooth rods with concrete, the contact elements are employed in the contact area between the reinforcement and concrete cylinder. The contact element involves two stiffness components: The tangential (frictional) and normal components. Another parameter, which is basically applied in the definition of contact, is the coefficient of friction µ. A tension cut-off model is used in the normal direction of interface element, while in compression it behaves linea r with stiffness coefficient Kn through the whole linear zone, as shown in Fig. 1. An elasto-plastic frictional behavior is considered for the tangential behavior of interface, in which the maximum tolerable tangential force; Fs, is µ Fn. When the concrete is casted in the prepared framework, it behaves as a fluid for a while without any capability to exert a pressure on the reinforcement inside. It gradually becomes more stiff and solid. Cement hydration produces a temperature gradient through the whole concrete and results in different temperature of the reinforcement and concrete. Considering different coefficient of temperature expansion the two materials possess, some initial pressure is exerted through the boundary of adjacent materials which is defined by some initial interference of contact element and is equal to radial stiffness Kn times the initial interference. In order to model the chemical adhesion between the concrete and rod, which is one of the suppliers of bond especially at initial stages of loading, some nonlinear springs have been used. The proposed equation for force- displacement curve of nonlinear springs by Al-Zahrani is Ae-Bd, where A and B are variable parameters (Bakis et al. 1998). It means that by increasing the amount of relative displacement, the total tolerable force of nonlinear springs, which indicates the adhesion force, decreases exponentially. In the above formula, d is the relative displacement. The model used in this analysis is axisymmetric so only one half of the experimental specimen is modeled. Fig. 2 illustrates a part of the model applied for simulation of smooth FRP rod pull out test in two modes: the initial and deformed configurations. In the deformed shape, both contact and nonlinear spring elements are shown. It can be observed that the length of these elements increases by exerting a pullout force, which states the opening of the interface, as soon as a longitudinal load is exerted on the rod. 4. Evaluation of contact parameters For the analysis of mechanical interaction between rod and concrete, some appropriate values for the contact element parameters Kn, Ks, µ, int (initial interference) should be evaluated. The procedure for evaluation of Kn, Ks can be made by assuming a predefined value for Kn, following a gradual increase in the value of Ks. For each corresponding pair of values of the parameters Kn and Ks, a separate analysis is made and the amount of final pull-out 63 Khosro Shahpoori Arani et al. Fig. 2 Typical shapes of smooth rod model (Deformed and Unreformed) (c) Fig. 3 (a) slip and load versus ks diagrams - kn evaluation , (b) Slip and load versus kn diagrams - ks evaluation , (c) Slip versus ks diagrams - kn evaluation 64 Computational optimized finite element modelling of mechanical interaction of concrete with fiber reinforced polymer (a) (b) Fig. 4 Bond strength versus slip diagrams: (a) µ=0.21 (b) µ=0.18 load along with the slip at the loaded end and final displacement of the reinforcement are derived. As can be seen in Fig. 3, for a specified point, any further increase in the value of Ks doesn ’t change the relative slip considerably. However, a small decremented change in the value of Ks results in merely high variation in the value of slips or end displacement. A high slope of the curve can be observed. This point is similar to “turned-point ” in mathematics and called here as “critical point ” (Fig. 3). Now a new iterative procedure is begun in which Ks is fixed at the amount derived from previous stage (the amount corresponding to critical point on the graph). As demonstrated in Fig. 3, by increasing Kn gradually, we can obtain a good approximation of Kn at the “critical point ”. A third procedure should be conducted similar to the first procedure. By increasing the value of Ks in an incremental manner, it is easily observed that the obtained value of Ks corresponding to the critical point is quite similar to that gained in the first procedure. This indicates that Ks equals to 5.0E7 and Kn equals to 8.0E8 are acceptable evaluations for the contact element parameters. The initial interference, which was derived from Al-Zahrani test (Bakis et al. 1998) , was about (3.5-3.6) µm. This value can be also concluded theoretically. According to Al-Zahrani test (Bakis et al. 1998), a 40°C decrease in temperature of the whole system must result in 0.7 MPa reductions in total average bond strength (the ultimate load divided by the lateral surface of the part of rod within the bond length). It is clear that the bigger value of coefficient of friction applied in simulation results in greater ultimate load and average bond strength. Thus µ will be increased gradually in the model and its variation in the bond strength due to 40°C decrease of temperature can be Fig. 5 Bond strength versus slip for contact parameters derived from calibration stage measured. The values of Ks and Kn, which have been applied in this analysis are exactly those derived from the iterative procedure explained above. The obtained coefficient of friction is 0.25. However, the proposed value from experiments is approximately between 0.18 and 0.2. As mentioned in previous sections, for a fixed value of Ks, increasing the value of Kn relative to its corresponding amount at the critical point doesn ’t affect the response of system very much. Of course it is true for µ=1 so that the effect of coefficient of friction on the system response is neglected and the tangential force always lies on the linear branch of frictional behavior curve. So Kn is increased up to 2.0E9 and the expected reduction in bond strength is gained for a µ parameter equal to 0.21. This value is very near to the value observed in tests (0.18-0.2) and could be proposed as an appropriate evaluation. 5. Calibration of contact parameters In order to calibrate the evaluated parameters, the load- slip curve has been applied as a suitable criterion. Al- Zahrani derived such a curve for an experimental sample of the model used in this research (Bakis et al. 1998). So parameters calculated from the previous stage ( Ks=5.0E7, Kn=2.0E9, µ=0.21) were inserted into the model and a reanalysis was carried out. Fig. 4(a) describes a comparison between the curves obtained from the theory and experiment. Although the mechanical behavior is similar in both cases, a considerable difference between the resulted maximum pull-out loads in both cases can be observed. To gain a numerical response, which coincides with observations in reality in a better manner, the coefficient of friction is reduced to 0.18. It can be seen that the maximum pull-out loads in both theory and experiments are quite identical. (Fig. 4(b) The best algorithm for a further decrease in the maximum tolerable pull-out load to make the uniform branch of graph as a good average approximation of the experimental data measured after the maximum load carried out is more reduction in the value of µ. So the analysis has been carried out one more time, applying the following parameters: Ks=5.0E7, Kn=2.0E9, µ=0.16 and the obtained load-slip curves for both ends of the rod for both experimental and numerical solutions are compared in Fig. 5. 65 Khosro Shahpoori Arani et al. It can be observed that the resulted curves in both cases are appropriately similar. Presence of an elasto-plastic behavior in the numerical analysis is because of using linear-full plastic material properties for transversely isotropic material of composite in various orthogonal directions due to computational and analytical limitations in the software applied for analysis. Also the mechanical behavior of concrete is assumed linear in modelling. From Fig. 5, it is concluded that these parameters are appropriate evaluations for the analysis of the mechanical interaction of FRP reinforcement with concrete. 6. Conclusions The purpose in the present paper was to propose a reliable finite element modelling of FRP smooth rods mechanical interaction with concrete by evaluation, calibration and validation of contact parameters. Friction is the main supplier of bond in smooth rods while mechanical interference (bearing) is the main cause of bond in lugged rods. Adding deformation increases bond strength and decreasing the deformation depth decreases the amount of bond strength. Modelling of smooth rods interaction with concrete should be done by applying contact elements. These elements are defined by these parameters: tangential stiffness ( Ks), normal stiffness ( Kn), coefficient of friction (µ) and initial interference ( δ). Also some nonlinear springs must be applied to simulate chemical adhesion. Process of proposing appropriate values for bond parameters Kn, Ks includes three stages: evaluation, calibration, validation. The suitable parameters resulted in this research for Ks is 5.0E7 and for Kn is 2.0E9. The suggested value for parameter µ from this research is 0.16. However, the proposed value from experiments is 0.18-0.20. Not all the parameters involved in reality can be considered in modelling. Also here the transversely isotropic material behavior can only be considered as linear-full plastic due to computational limitations so the numerical curves show a smoother and uniform behavior, although the experimental curves are very curved & crooked.