Advanced composite material simulation

Proceedings: 11th Portugese Conference on Fracture. Pp 3-22. 13th – 15th February 2008. Eds.: R. Martins, C. Moura Branco, J. Teixeiras, F. Fernandes. ISBN: 978-989-95683-0-3 ADVANCED COMPOSITE MATERIAL SIMULATION S. Oller*, X. Martínez*, A. Barbat*, F. Rastellini** *Departamento de Resistencia de Materiales y Estructuras en la Ingeniería (RMEE) Universidad Politécnica de Cataluña **International Center for Numerical Methods in Engineering (CIMNE) Building C1, Campus Nord UPC Jordi Girona, 1-3. 08034-Barcelona, Spain. e-mail: sergio.oller@upc.edu Abstract A computational methodology is presented for modelling the material non-linear mechanical behaviour of composite structures made of FRP (Fibre-Reinforced Polymers) laminates. The model is based on the appropriate combination of the constitutive models of component materials, considered to behave as isolated continua, together with additional ‘closure equations’ that characterize the micromechanics of the composite from a morphological point of view. To this end, any appropriate constitutive model may be selected for each phase. Each component is modelled separately and the global response is obtained by assembling all contributions taking into account the interactions between components in a general phenomenological way. To model the behaviour of a single unidirectional (UD) composite lamina, a Serial-Parallel continuum approach has been developed assuming that components behave as parallel materials in the fibres alignment direction and as serial materials in orthogonal directions. Taking into account the internal morphology of the composite material, it is devised a strategy for decoupling and coupling component phases. This methodology [Rastellini 2006], named “compounding of behaviours”, allows to take into consideration local degradation phenomena (in the constituents materials), like plasticity, etc. in a coupled manner. It is based on the proper management of homogenous constitutive models, already available for each component. In this way, it is used all developments achieved in constitutive modelling of plain materials, what makes possible the transference of this technology to composites. A lamination theory complemented with the proposed UD model is employed to describe the mechanical behaviour of multidirectional laminates. A specific solution strategy for the general non linear case is proposed. It provides quick local and global convergences, what makes the model suitable for large scale structures. The model brings answers on the non-linear behaviour of composites, where classical micro-mechanics formulas are restricted to their linear elastic part. The methodology is validated through several numerical analyses and contrasted against experimental data and benchmark tests. 1. INTRODUCTION The study of composite materials has been one of the major objectives of computational mechanics in the last decade. The numerical simulation of composite materials has been done, traditionally, using orthotropic materials with average properties of their constituents. With this approximation, no model has been found able to work beyond the constituents elastic limit state. Thus, these procedures limit the numerical computation to elastic cases. Different theories have been proposed to solve this problem, taking into account the internal configuration of the composite to predict its behaviour. The two most commonly used are herein remarked. Homogenization theory: This method deals with the global problem of composite material in a two-scale context. The macroscopic scale uses the composite materials to obtain the global response of the structure; composites are treated as homogeneous materials in this scale. The microscopic scale corresponds to an elemental characteristic volume in which the microscopic fields inside the composite are obtained; this scale deals with the component materials. Homogenization theory assumes a periodical configuration of the composite material to relate these two scales [SanchezPalencia, 1987; Oller et al., 2005]. Mixing theory: The first formulation of the mixing theory corresponds to Trusdell and Toupín [Trusdell and Toupín, 1960] and it is based in two main hypothesis: 1. All composite constituents have the same strains. 2. Each constituent collaborates to the composite behaviour according to its volumetric participation. The main problem of the mixing theory is the iso-strain condition, which forces a parallel distribution of the constituents in the composite. Some improvements to the original formulation can be found in [Car et al., 2000]. This work uses an improvement of the mixing theory. The election of the mixing theory instead of a homogenization theory is based in the final aim of the code developed, compute real structures made of composite materials. A homogenization theory requires a micro-model for each point of the structure that becames non-linear. The resolution of a real structure with this procedure generates such a big amount of degrees of freedom that the calculation is beyond the computation capabilities of nowadays personal computers or small servers. On the other hand, the mixing theory does not increase the degrees of freedom of the problem, as it is only present in the constitutive section of the finite element code. The simulations that will be shown are solved with PLCd [PLCd, 1991]. This is a finite element code, developed at CIMNE (International Center for Numerical Methods in Engineering) and UPC (Politechnical Univeristy of Catalonia). The code works with two and three dimensional solid geometries. It can deal with kinematics and material nonlinearities. It contains various constitutive laws to predict the material behaviour: Von-Mises, Mohr-Coulomb, improved Mohr-Coulomb, DruckerPrager, etc. [Lubliner et al., 1989]. It uses different integration algorithms to simulate the material: Elastic, visco-elastic, damage, damage-plasticity, etc. [Oller et al., 1990]. Dynamic analysis are developed using the Newmark method, [Barbat et al., 1997]. The different formulations included in PLCd to deal with composite materials are exposed hereafter. 2.FORMULATION USED TO SIMULATE THE COMPOSITE MATERIAL BEHAVIOUR 2.1 MIXING THEORY 2.1.1 Classical mixing theory The classical rule of mixtures was originally developed by Trusdell and Toupin [1960]. It considers that the interaction between the components in a material point of the composite is done according to the following hypothesis: (i) each infinitesimal volume of the composite contains a finite number of material components; (ii) each component contribution to the global behaviour of the composite is proportional to its volumetric participation; (iii) all components suffer the same strains (closing equation); (iv) the volume of each component is significantly smaller than the composite volume. The third hypothesis, in the case of small strains, can be written as: ε ij = 1ε ij = 2ε ij = … = nε ij c (1) Where, ε ij is the strain tensor for the composite and ε ij is the strain tensor for component k of the composite. According to second hypothesis, the stresses of the composite can be computed as the proportional addition of each component stresses, according to the volumetric participation, thus: c k σ ij = ∑ k =1, n k k ⋅ kσ ij = ∑ k =1, n k S k ⋅ k Cijkl ⋅ k ε ij = ∑ k =1, n k S k ⋅ k Cijkl ⋅ cε ij (2) where parameter k k is the volumetric participation of component k in the composite, defined as: k k= dVk dV0 (3) A more detailed explanation of this theory, as well as the extension to the large deformations case and its numerical implementation can be found in [Car, 2000] 2.1.2 Serial/Parallel mixing theory The main problem of the classical mixing theory is the iso-strain condition, which forces a parallel distribution of the constituents in the composite (Figure 1, iso-strain case). The serial/parallel rule of mixtures improves the classical mixing theory replacing the iso-strain hypothesis for an iso-strain condition in the fibre direction and an iso-stress condition in the transversal directions, allowing to simulate all components distribution in the composite shown in Figure 1. This theory has been developed by Rastellini and an extensive description of it can be found in [Rastellini, 2006] and in the Deliverable 49 of the LESSLOSS project [CIMNE, 2006]. Iso-strain Iso-stress Mixed case Figure 1. Different distribution of components in a composite material The main hypothesis in which is based the numerical model of the Serial/Parallel mixing theory are: i) ii) iii) iv) Composite is composed by two component materials: fibre and matrix Component materials have the same strain in parallel (fibre) direction Component materials have the same stress in serial direction Composite material response is in direct relation with the volume fractions of compounding materials v) Homogeneous distribution of phases is considered in the composite vi) Perfect bounding between components is also considered The equations that define the stress equilibrium and establish the strain compatibility between components arise form the analysis of the model hypothesis. Thus, Parallel behaviour: c ε P = mε P = f ε P σ P = mk mσ P + f k f σ P c Serial behaviour: c ε S = mk mε S + f k f ε S σ S = mσ S = f σ S c (4) (5) where, ε P and ε S are the parallel and serial components of the stress tensor respectively, σ P and σ S are the parallel and serial components of the strain, the superscripts c , m and f represent the m f composite, matrix and fibre materials and k and k is the volumetric participation of fibre and matrix in the composite. In order to verify the compatibility equations, the Serial/Parallel mixing theory proceeds in the following way. 1. First thing to be done is to split the component strain tensor into its serial and parallel parts. 2. 3. 4. Afterwards, a first prediction of the matrix serial strains is computed. With this initial prediction and the strain relation of equation (5), the fibre serial strains can be also obtained. Using each constituent constitutive law, the stresses for matrix and fibre can be calculated, as well as the actualization of each material internal variable. The stress tensors computed in step 3 are split into their serial and parallel parts and the isostress condition defined in equation (5) is checked. If it is verified, the matrix strain prediction was correct and the composite stress can be obtained using equations (4) and (5). If the isostress condition is not verified, the initial strain prediction must be corrected and the process must continue in step 3. The flow diagram implemented in the numerical code corresponding to the exposed procedure is shown in Figure 2. Figure 2. Flow diagram of the Serial/Parallel mixing theory (SP RoM) The main problem of the Serial/Parallel mixing theory is that the composite material can be composed, only, by two material components. However, this inconvenient is solved with the laminate procedure exposed in the following section. 2.2 LAMINATE COMPOSITES Laminated composites are composed by different layers with different fibre orientations. The orientation of the fibre can be defined by the engineer in order to obtain the better performance of the composite according to its application. Under this scope, the limitation of the Serial/Parallel mixing theory to only two materials is not such, as the composites used are usually defined by multiple layers composed by only two components: fibre and matrix. Thus, the SP RoM formulation can be applied to each layer of the composite and, afterwards, compute the composite behaviour composing the performance of each constituent layer. To obtain the laminate behaviour, the classical mixing theory is applied to each lamina. Applying the classical mixing theory to the different layers of the composite implies the assumption that all laminas of the laminate are under the same strain conditions. This assumption can be considered correct, as the different laminas usually have fibre orientation distributions disposed in such a way that provide the laminate with an in-plane homogeneous stiffness. At this point, it must be said that not all the layers have to behave according to the S/P RoM hypothesis. Sometimes, between fibre oriented layers, a randomly oriented one is disposed. Having random oriented fibres, the layer behaves as an isotropic material and the best theory to be applied to predict its behaviour is the Classical mixing theory. Thus, the composite performance is obtained according to the following procedure: 1. For a given strain of the composite, the strain of each composite layer is obtained applying the Classical mixing theory. c 2. (6) Each layer stress is obtained using the Serial/Parallel mixing theory in case of fibre oriented layers or the Classical Mixing theory in case of randomly oriented layers (or in other cases of iso-strain behaviour such as single materials layers). Lk 3. ε = L1ε = L 2ε = … = Lnε ε → Lk σ; (7) using Classic mixing theory or SP RoM The stresses of the composite are obtained composing the stresses in each composite layer. n σ = ∑ Lj k Ljσ c (8) j =1 2.3 OTHER FORMULATIONS DEVELOPED TO SIMULATE CFRP REINFORCEMENTS Having defined the main formulation frame, in which the code works and deals with composite materials, there are other formulations in it that are used to obtain a better performance of the finite element code and that are required to obtain a better approximation of the mechanical behaviour of composite material structures. In this chapter, all of them are briefly described. 2.3.1 Anisotropy using a mapping space theory This theory is based on the transport of all the constitutive parameters and the stress and strain states of the structure, from a real anisotropic space, to a fictitious isotropic space. Once all variables are in the fictitious isotropic space, an isotropic constitutive model can be used to obtain the new structure configuration. This theory allows considering materials with high anisotropy, such as composite materials, using all the techniques and procedures already developed for isotropic materials. σ All the anisotropy information is contained in two fourth order tensors. One of them, Aijkl , relates the stresses in the fictitious isotropic space ( σ ij ) with the stresses in the real anisotropic space ( σ ij ) and ε the other one, Aijkl , does the same with the strains. The relation of both spaces for the strains and the stresses is exposed in equations (9). σ σ ij = Aijkl : σ ij ; ε ε ij = Aijkl : ε ij (9) A representation of these transformations is displayed in Figure 3. A more detailed description of this methodology, the extension to large strains and its numerical implementation can be obtained in [Car et al., 2001 and Car, 2000] Figure 3. Space transformations. Real and fictitious stress and strain spaces in small strain (image obtained form [Car et al., 2001]) 2.3.2 Fibre-matrix debounding The apparition of matrix cracks in a composite material is usually followed by a relative movement between the fibres and the matrix. This lost of adherence implies a stiffness reduction in the composite material. This phenomenon is introduced in the elastic limit of the material as a modification of its yield surface criterion. The new fibre elastic limit becomes: [ { ( f R ) fib = min ( f N ) fib ; ( f N ) mat ; 2 ⋅ ( f N ) fib− mat / r fib R N ]} (10) N Where ( f ) fib is the new fibre strength, ( f ) fib is the nominal fibre strength, ( f ) mat is the N matrix nominal strength and ( f ) fib − mat is the fibre-matrix interface nominal strength. Equation (10) shows that the debounding happens when one of the composite constituents reaches its nominal strength (considering the fibre-matrix interface as a constituent). The nominal resistance values are obtained from the material properties. The numerical implementation of this phenomenon is described in the reference [Car, 2000 and Oller, 2002] 2.3.3 Construction stages algorithm When performing numerical simulations of structures retrofitted with FRP, it is necessary to add the structural reinforcement once the original structure is already damaged. The “Construction Stages Algorithm” implemented in PLCD permits running the numerical simulation during the desired load cases, with only some structural elements active on the structure. At a certain load case, new elements can be added without stopping the calculation process. These new elements are free of strains and stresses when they are activated. The algorithm requires having all elements defined in the structure. The elastic strains are divided in two components, and active and a non-active, equation (11). If the element is not present in the structure, all strains correspond to the non-active part, equation (12), while, if the element is active, all strains corresponding to the non-active situation will be removed, equation (13), from the total strain. Stresses are computed considering only the active elements, equation (14). e ε e = ε Ae + ε NA e ε NA = ε e; (11) ε Ae = 0 (12) e ε Ae = ε e − ε NA (13) σ e = C e ⋅ ε Ae (14) This procedure is explained with more detail in [Martinez et al., 2006; CIMNE, 2006] 2.3.4 Tangent constitutive stiffness tensor Depending on the constitutive equation used in a composite constituent material, the tangent constitutive tensor cannot be obtained analytically. One solution is to use, in these materials, the initial stiffness matrix, which will lead to the equilibrium state but will require a large amount of structural iterations. Thus, in order to obtain a fast and reliable algorithm, the expression of the tangent constitutive tensor is required. To obtain it, when no analytical expression exists, a numerical derivation using a perturbation method is performed. The definition of the tangent constitutive tensor is: σ = Ct : ε (15) where, σ = [σ 1 σ 2 … σ n ] ; ε = [ε1 ε 2 … ε n ] and ⎡ c11t … c1tn ⎤ ⎢ ⎥ Ct = ⎢ ⎥ t ⎥ ⎢ cnt 1 … cnn ⎣ ⎦ (16) The definition of the tangent constitutive tensor, equation (15), shows that the variation of stresses due t to an increment in the value of the j element of ε depends on the values of the j column of C . t Thus, writing the j column of C as, ctj = ⎡⎣c1t j c2t j … cnjt ⎤⎦ T (17) the stress variation is: σ j = ctj ⋅ ε j (18) t being c j the unknown The perturbation method consists of applying a small perturbation to the strain vector and, using the constitutive equation of the material, determine the variation that will be obtained in the stress tensor due to this perturbation. At this point, the j column of the tangent constitutive tensor can be computed as: ctj = σj εj (19) The smaller the value applied to the perturbation, the better the approximation obtained for the tangent constitutive tensor. Having defined a perturbation value, ε j , the perturbed stress is computed using the constitutive equation of the material applying the following input strain: ε = ⎡⎣ε1 … ε j + ε j … ε n ⎤⎦ T → σˆ (20) And the stress variation due to the perturbation is obtained subtracting the original converged stress from the computed one: σ j = σˆ − σ (21) This procedure must be repeated for all strain components in order to obtain the complete expression of the tangent constitutive tensor. Hence, the numerical cost of using a perturbation method is rather high. However, this procedure allows obtaining an accurate approximation to this tensor for any constitutive equation used, ensuring the convergence of the numerical process in few steps. 3.NUMERICAL EXAMPLES OF THE FORMULATION PROPOSED. 3.1 BENDING REINFORCEMENT OF A RC BEAM In this section, a numerical simulation of a bending reinforced beam is presented and used to validate all the formulations proposed and exposed in previous section. This case shows the efficiency of the serial/parallel rule of mixtures to deal with this sort of structural problems, as it is able to reproduce the complex mechanical behaviour found in the beam with an acceptable computational cost. The numerical results are validated with experimental values. The studied beam is defined in the paper by Spadea [Spadea et al. 1998]. Its geometry and the reinforcements applied to it are shown in next figure: Figure 4. Geometry and reinforcement of the beam studied The red (thick) line displayed in the bottom of the beam corresponds to the FRP reinforcement. This is made of carbon fibres embedded in a polymeric matrix. The content of fibres is 60% and the composite thickness is 1.2 mm. The finite element model developed to simulate the beam reinforcement is shown in Figure 5. This is a 3D finite element made with linear hexahedrons. The S/P rule of mixtures allows considering all the reinforcement details found in the beam using a coarse mesh. In Figure 5 it is also included the composite materials composition. As can be seen, a single finite element contains, in this particular case, up to three different component materials. The steel reinforcements are considered as fibres, whose orientation is defined by the bar direction. The FRP reinforcement has been included adding new finite elements to reproduce with more accuracy its position in the beam. The results obtained with this simulation are compared with the experimental results reported by Spadea [Spadea et al., 1998]. Figure 6 shows the capacity curve of the beam for the numerical and the experimental simulations. This is, vertical displacement of the beam, in the point where the force is applied, against the load value. This figure shows the agreement between numerical and experimental results, which proves the ability of the method to perform this sort of simulations. Figure 6 also includes the results obtained with a numerical simulation of the same beam without FRP. The comparison between the results obtained for the reinforced and for the non-reinforced beam shows the improvement obtained in the beam performance when it is reinforced with FRP. MAT-01: Concrete (100%) MAT-02: Concrete (57%) Longitudinal Steel (42%) Vertical Steel (1%) MAT-03: Concrete (99%) Vertical Steel (1%) MAT-04: Concrete (98%) Vertical Steel (1%) Horizontal Steel (1%) MAT-05: Concrete (99%) Horizontal Steel (1%) MAT-06: Polymeric Matrix (34%) Long.Carbon Fibr. (66%) Figure 5. Finite element model developed to realize the numerical simulation Figure 6. Force-displacement graph comparing the experimental and the numerical results One of the main advantages of the proposed finite element formulation is that it allows obtaining the structural behaviour of each composite component, providing all information required to know the failure cause of the structure. In example, for the simulation developed, the failure modes reported by Spadea are: Tension steel yielding and concrete crushing, in the non reinforced model; and sudden and total loss of load capacity due to an explosive debonding of CFRP plate, in the reinforced model. The failure mode in the non-reinforced beam model is the same as the one obtained in the experimental test. This can be seen in figure Figure 7, which shows the damage to the concrete and the plastic damage to the steel reinforcement for the final calculation step and for the most severely damaged beam section. These figures show that the steel has started to yield and that the compressed concrete has also reached its limit stress (onset of concrete crushing). Under these conditions the transversal section cannot develop more stresses and the code cannot determine a valid solution for the current load step. This situation can be interpreted as beam failure. Figure 7. Damage parameter in concrete and steel in the non-reinforced beam at failure 3.2 CFRP RETROFITTING OF A BEAM Two different numerical models have been developed to study the effect of retrofitting a structure, depending on the existing damage in the beam when the CFRP reinforcement is applied to it: • • Sp3D-Rt2: The CFRP reinforcement is applied when damage starts in concrete material. Sp3D-Rt3: CFRP reinforcement is applied when steel starts yielding. Results obtained with these two models are compared with those obtained when the beam is not reinforced (Sp3D-R0 model) and when the beam is reinforced from the beginning of the loading process (Sp3D-R1 model). The capacity curve obtained for each model is shown in Figure 8: These results show that the structural stiffness does not depends on when the reinforcement is applied to the structure. The structure stiffness obtained when the CFRP reinforcement is applied after the steel yielding (Sp3D-Rt3 model) does not differ significantly form the structure stiffness obtained after steel yielding in the reinforced model (Sp3D-R1). However, retrofitting a structure implies that, when the CFRP reinforcement starts collaborating, the deformation and damage of the structure is larger than if it had been reinforced from the beginning. Damage reduces the load capacity of the beam while deformation can make the serviceability state unacceptable (i.e. when a load of 25 kN is applied to the structure, beam deformations are 45% larger in the retrofit model, Sp3D-Rt3, than in the reinforced one, Sp3D-R1). Figure 8. Comparison between the CFRP reinforcement and retrofitting by using capacity curves 3.3 REINFORCEMENT OF A FRAMED STRUCTURE USING CFRP The main objective of the present simulation is to use the developed formulation to verify the capability of CFRP reinforcements to increase the strength of concrete frames when seismic loads are applied to them. Concrete framed structures are common in building and civil engineering; one of the most stressed zones of these structures, under seismic loads, are the connecting joints between beams and columns. In many occasions, these joints show a lack of strength that can be improved with CFRP. The developed models reinforce the frame joint with two different CFRP configurations, in order to study the ability of the reinforcements to increase the frame strength and to find out which configuration offers better results. 3.3.1 Model description The concrete frame to be studied is defined with the geometry and the steel reinforcement commonly used in buildings. Figure 9 shows the geometry considered and Figure 10 shows the steel reinforcement and the description of the CFRP reinforcement that will be applied to the frame joint. The beam cross section, as well as its steel reinforcement, is dimensioned in such a way to assure the structural failure in the beam, near the joint, in order to increase the effects of CFRP in the frame joint. The structure is loaded with an horizontal force P applied in the middle of the frame joint (Figure 9). Figure 9. Geometric definition of the framed structure considered in the simulation Figure 10. Reinforcements applied to the concrete frame Two dimensional and three dimensional models have been developed for the concrete frame. 2D models have been used to calibrate the mesh since it require less computational effort than 3D models. 2D models results are compared with the 3D ones to validate the accuracy of each simulation. Three different structures are considered to study the effect of CFRP reinforcements on the frame joint: 2D-noR and 3D-noR: Two and three dimension concrete frame without CFRP reinforcements 2D-R and 3D-R: Two and three dimension concrete frame with the upper and lower CFRP reinforcements defined in Figure 10 2D-LR and 3D-LR: Two and three dimension concrete frame with the upper, lower and lateral CFRP reinforcements All composite materials existing in the concrete frame are defined by combination of four different basic materials, which are defined in Table 1. CFRP reinforcement is 1.2mm thick and is composed of 66% of carbon fibres and 34% of polymeric matrix. In the case of the upper and lower reinforcements, the fibres are oriented following the structure longitudinal axis. In lateral reinforcements, two layers of CFRP are applied to the frame, in which the fibres are oriented at +0º and +90º respect the horizontal. Yield criterion Material Young Modulus [ MPa ] Poisson Modulus Yield Stress [ MPa ] Fracture Energy [kPa·m] Compr. Tensile Compr. Tensile Mohr-Coulomb Concrete 2.5·104 0.20 30.0 3.0 50.0 0.5 5 Von-Mises Steel Reinf. 2.1·10 0.00 270.0 270.0 2000.0 2000.0 Mohr-Coulomb Polymeric matrix 1.2·104 0.20 87.5 29.2 36.0 3.0 Von-Mises Carbon Fibres 1.5·10 0.00 2300.0 2300.0 2000.0 2000.0 Table 1. 5 Mechanical characteristics of the constituent materials used to define the composite materials existing in the framed structure. 3.3.2 2D results The structural behaviour of the frame joint for the different reinforcements applied to the structure is studied using the capacity curves obtained for each model (Figure 11). The displacement considered corresponds to the horizontal displacement suffered by the point where the load is applied. This displacement depends on the column, beam and joint stiffness. As the column and the beam are not modified in the different models, if the joint stiffness is increased with the different applied reinforcements, the force-displacement graph will reflect this increment. Figure 11. Capacity curves obtained with the 2D models Figure 11 shows that the upper and lower CFRP reinforcements (2DF-R model) do not improve significantly the frame behaviour. This improvement is only found when lateral reinforcements are applied to the concrete frame. All three curves present a region where the load reduces, to start increasing again afterwards. These points correspond to the development of a plastic hinge in the structure. At this load step the structure adopts a new strength mechanism and can increase its load capacity. Comparing the load applied to the structure until the development of the plastic hinge, the lateral reinforcements (2DF-LR model) increase in a 25% the structural load capacity when compared with the non-reinforced model (2DF-noR). This increment is only a 4% if the structure is reinforced only with upper and lower CFRP. A better comprehension of the effects of each reinforcement can be obtained studying the points where the plastic hinges are formed. Figure 12 shows the longitudinal strains of each model at their last computed step. The cross sections where the plastic hinges are formed are the ones with larger strains. Figure 12 shows that applying only the upper and lower CFRP reinforcements the plastic hinge moves from the beam to the inner part of the joint, where no reinforcements are applied. Thus, the presence of CFRP does not modify substantially the beam behaviour and, once the hinge has been formed, both structures behave similarly (as it can be seen in Figure 11). On the other hand, when the lateral reinforcement is applied to the structure, it restrains damage in the frame joint and moves the plastic hinge to the cross section where no CFRP reinforcement is applied, what allows the structure to increase its load capacity and its stiffness. 3.3.3 3D results The study by means of three dimensional models is also performed, as in the case of 2D models, by means of the capacity curves (Figure 13). The main difference found when comparing these results with the 2D ones is that the 3D models are stiffer and can reach larger loads than the 2D models. This is because the concrete confinement is better reproduced in this case, as steel stirrups are modelled taking into account their 3D distribution and not only in one of their directions. Hence, concrete can reach larger stresses, which increase the stiffness and strength of the structure. (a) (b) (c) Figure 12. Plastic hinges in the concrete frame. 2D models. a) model without CFRP reinforcement, b) model with upper and lower CFRP, c) model with upper, lower and lateral CFRP Figure 13 shows that plastic hinges in the non-reinforced model (3DF-noR) and in the reinforced model (3DF-R) appear for the same load and displacement as in the reinforced beam model (2DF-LR), as a consequence of the increment of concrete strength. However, in the three dimensional simulation, plastic hinges appear before in the reinforced model than in the non-reinforced one. The explanation of this effect is shown in Figure 14 (maximum strains in the unreinforced beam model (3DF-noR) before and after the formation of the plastic hinge) and in Figure 15 (the same results in the case of the 3DF-R model). Figure 13. Capacity curves obtained with the 3D models According to these figures, the cross section in which the plastic hinge develops is nearly the same in both models. But, as this cross section is closer to the initial damage in the reinforced model than in the non-reinforced one, it is easier to simulate the plastic hinge when the beam is reinforced. Thus, even if the CFRP reinforcement increases the joint stiffness, in this case the frame plastic hinge appears for lower loads when this reinforcement is applied than when the joint is not reinforced. More differences are found when comparing the three dimensional model with the two dimensional one, in the case in which lateral CFRP reinforcements are applied to the frame joint (3DF-LR model). The first difference is that the formation of a plastic hinge is not visible in the capacity curve. This is because no section is completely damaged when the algorithm loses its convergence. However, the main difference is found when looking at the most damaged section. The strains in the lateral sections of the frame joint (Figure 16a) have a similar behaviour to that of the 2D case: strains are larger in the cross section where the CFRP reinforcement finishes than in the frame joint. But when the strains in a longitudinal section of the structure are studied (Figure 16b), they show that the plastic hinge is developed in the frame joint, as happens in the 3DF-R model. Two dimensional models consider the CFRP reinforcement applied along the whole cross section while actually it is applied only to the lateral surfaces. Thus, the reinforcement can avoid structural cracks on the surface of the structure but cannot restrain them inside the joint. These last figures show that the structure presents the same structural failure independently of the CFRP reinforcement configuration applied to it. Thus, it can be concluded that the only effect of the lateral CFRP reinforcement over the frame joint is to delay the apparition of cracks in it and the following plastic hinge. However, this delay is enough to allow a load 20% larger in the frame when the horizontal displacement in it is of 3.0 cm, increment more than sufficient to consider this reinforcement typology the best option to improve concrete framed structures seismic strength. The results also show the necessity of using three dimensional elements depending on the simulation to be performed. When the distribution of CFRP reinforcements is not uniform along the whole cross section, 2D simulations can provide incorrect results. However, even if the effectiveness of lateral reinforcements is reduced in the 3D model, the structural strength improvement is significant enough to consider this reinforcement configuration as the best option to reinforce the column-beam joint of RC frame structures. Figure 14. Crack evolution in the 3DF-noR model (model without CFRP reinforcement) Figure 15. Crack evolution in the 3DF-R model (model with upper and lower CFRP) Figure 16. Plastic hinge in the 3DF-LR model. Lateral view 3.4 DELAMINATION AND INTERLAMINAR CRACK PROPAGATION The developed formulation can be used also to simulate effects such as delamination and interlaminar crack propagation in laminated composites. The present simulation shows the capability of the code to deal with these sorts of problems. This concept is also usefullfor retrofitting applications. 3.4.1 Experimental test description The End Notch Flexure test is based in the flexure of a beam with an initial crack in one of its endings. The test has been applied to a composite made of carbon fibres with an epoxi polymeric matrix. Fibres are oriented in the longitudinal direction of the beam and the initial crack is created introducing an insert in the laminate during its fabrication. The dimensions of the beam, as well as the dimensions considered for the numerical simulation, are shown in Figure 17. The experimental test applies a vertical displacement to the beam until the initial crack starts its propagation. The imposed displacement is applied until the crack progression stops and the beam recovers its linear behaviour. At this point, the sample is unloaded. Main results obtained from this test are two: The force-displacement graph, which shows the structural performance of the composite beam, and the final length of the initial crack. These two results are the ones that will be compared with the numerical model developed. Figure 17. Dimensions of the beam used in the ENF test The exact properties of the composite material used in the experimental simulations were unknown when the experimental tests were performed [Corbella et al., 2004]. However, the composite is known to be made of carbon fibres and an epoxy polymeric matrix from Hexcel composites. For the numerical simulation, the mechanical values considered to define the composite are the ones described in Table 2, obtained form Hexcel Product data description. The fibre (AS4) and matrix considered are the ones found in HexPly 8552 UD carbon prepegs. Table 2. Sample geometry used for the ENF test 3.4.2 Numerical Model Description Two diferent numerical models have been used to simulate the “End Notch Flexure Test”. One test is in two dimensional formulation plain stress and a the second one using a three dimensional formulation. The 2D model has been defined with 627 linear quadrilateral elements while the 3D model has 5016 linear brick elements. The mesh defined for the three dimensional model is shown in Figure 18. Figure 18. Three dimensional model developed. Mesh description Two diferent materials have been defined in the numerical simulation. One of them corresponding to the composite material and another one corresponding to the insert material. The composite material is defined by the properties of the epoxy matrix and the carbon fibres shown in Table 2. Fibre material is defined as an elastic material. Matrix material is characterized by a damage law, like the one defined on reference [Martinez et al. 2007; Lubliner et al. 1989; Rastellini et al. 2006; Oliver et al. 1990]. The damage model used requires knowing the relation between the compression strength and the tension strength of the material in order to obtain the correct damage evolution. As these parameters are unknown, what is done in the present document is to consider that both strengths the same and to define the fracture energy of the material as the mode II fracture energy obtained from the experimental results. The value of this energy, for the 3M101 beam, is: GII = 1.02 kJ / m 2 . The definition of the insert material properties has been done taking into account its structural performance. The main efect of this material in the beam is allowing the sliding of the section found above the insert along the section found below it. To do so, a material with a shear modulus nearly zero has been defined (it has not been defined as strictly zero to avoid numerical instabilities during the simulation). 3.4.3 Comparisson between the numerical and the experimental results. The numerical and the experimental results are compared with the force-displacement graph obtained for both cases. The displacement represented corresponds to the vertical deflection of the point where the load is applied. This graph is shown in Figure 19, in which the results for the 2D and 3D numerical simulations and the experimental test are represented. Figure 19. Force-displacement. Graph obtained for the different models. Figure 19 shows that the two dimensional simulation provides exactly same results as the three dimensional one, thus, for these kind of problems, 2D simulations are preferable, as the computer cost is much lower. However, the most important result shown in Figure 19 is the agreement between the numerical and the experimental results. The beam initial elastic stifness obtained in the numerical simulation is nearly the same that is obtained from the experimental test. And this agreement between results is even better when comparing the beam maximum load capacity or failure point. The only result that difers a bit in the numerical simulation is the final beam stifness, when the crack has reached its maximum length. In this case, the numerical beam is a 6 % stifer than the experimental one (the stifness obtained in each case is, respectively, 1146 N/mm and 1076 N/mm). The other result to be compared is the final crack length. The experimental values obtained for this final crack length is of 50.34mm. In the numerical simulation, the crack points correspond to those in which the damage parameter, in matrix material, is equal to one. These points have a matrix stifness equal to zero. This implies that the composite serial stifness is also zero, due to the iso-stress condition imposed by the Serial/Parallel mixing theory. Those points with matrix completely damaged cannot develop any shear strength. Hence, the final crack length can by obtained by finding the point, closer to the beam mid-span, with a value of the damage parameter, in matrix material, equal to one. Figure 20 shows the damage parameter in matrix in the load step in which the beam reaches its maximum deflection. In this figure can be sen that the crack length obtained with the numerical simulation also nearly reaches the mid-span section. Figure 20. Damage in matrix material when the maximum deflection has been reached The exact value of the damage parameter is shown, for the points represented in Figure 21a, in Figure 21b. In this figure can be seen that point 13 (corresponding to mid-span) reaches a damage value of 0.6, while the value of point 12 is approximately 0.98. Considering this ast value close enough to one and thus, the section completely broken, the numerical crack length obtained is of 48mm. The point found at 49mm of the support has a damage value in matrix material of 0.89, which is also close enough to one to consider that the numerical results are exactly the same as the experimental ones. Figure 21. Evolution of the damage parameter in the beam 4.CONCLUSIONS In this work have been exposed the different numerical procedures developed to solve in a reliable and efficient way the problem of reinforcement and/or retrofitting of reinforced concrete (RC) structures using fibre reinforced polymers (FRP) and composite structural behaviour with delamination and interlaminar crack propagation. Due to the complexity of the problem to be solved, efforts have been directed not only to the numerical procedures that allow performing the structural simulation but also to the efficiency of the code and the way the user interacts with it. Of all the formulations developed to solve the delamination problem and FRP reinforcement of RC structures, is of special relevance the general formulation of the mixing theory to deal with laminated composites. This is the division of the composite in its different components until reaching the constitutent material, which will be the one that will provide the structural behaviour of the composite. This formulation can be understood as a manager of constitutive equations and is the one that allows dealing with the reinforcement problem, taking into account all its particularities, without increasing the problem numerical size beyond the computational limits. The good performance of this theory has been tested with the different simulations presented in section 3 of the present document. Besides the mixing theory formulation, other numerical procedures have been developed to deal with the particular case of FRP reinforcements. The most relevant of them are the anisotrophy using a mapped space theory and the fibre matrix debounding. Special attention must be paid in the construction stages algorithm, which allows performing simulations of retrofitted structures (section 3.2). These simulations have shown that, even the FRP performance do not vary significantly if it is applied as a reinforcement or as a retrofit, the structural deformations and the stresses are larger if the FRP is applied when the structure is already damaged. Finally, it has to be said that the compression strength formulation has not been tested yet in a FRP reinforcement simulation, however, first results obtained permit be optimistic about its performance and validity. Also, the paper shows that the Serial/Parallel mixing theory is able to simulate a composite laminate behaviour ant its delamination process by using the appropriate constitutive equations to predict the material behaviour of the composite components. No other help or formulation is needed. This afirmation has been proved with the simulation of the End Notch Flexure (ENF) test. The results obtained with the numerical simulation are practically identical to the results obtained from the experimental tests. The structural behaviour of the materials, obtained from the numerical simulation, shows that the delamination phenomenon occurs as a result of the lost of stifness in matrix material due to the damage produced by the shear stresses in it. This lost of stifness in matrix material implies that no other component material can develop serial stresses, due to the iso-stress condition of the Serial/Parallel mixing theory in the serial direction. This is translated in a lost of stifness of the composite in all serial directions. In the case considered, the serial directions are all directions non coincident with fibre orientation. Thus, the composite is unable to avoid the shear deformations produced by the external loads. The conclusion that all all formulations exposed are valid and usefull to solve composite laminated structures of FRP. However, the objective of these simulations was not only to verify the formulations implemented in the finite element code, but also to study in which way the simulations of structural reinforcements of CFRP must be developed to obtain accurate results. The simulation that has provided more information in this aspect has been the framed structure one. 5.ACKNOWLEDGEMENTS This work has been supported by CEE-FP6 (LESSLOSS project, ref. FP6-50544(GOCE)) and by the Spanish Goverment through the Ministerio de Ciencia y Tecnolog\'ia (RECOMP project, ref. BIA200506952 and the Ministerio de Fomento (project “Retrofitting and reinforcement of reinforced concrete structures with composite materials. Numerical and experimental developments applied to joint of bars and composites anchorage proposal”). The fellowship of X. Martinez has been provided by CIMNE (International Centre for Numerical Methods in Engineering). All this support is gratefully acknowledged. 6.REFERENCES [1] Barbat, A.H., Oller, S., Oñate, E., Hanganu, A. [1997] “Viscous damage model for Timoshenko beam structures” International Journal of Solids and Structures, Vol.34, No 30, pp. 3953-3976 [2] Car, E., Oller, S., Oñate, E. [2001] “A Large Strain Plasticity for Anisotropic Materials – Composite Material Application,” International Journal of Plasticity, Vol. 17, No. 11, pp. 1437-1463. [3] Car, E., Oller, S., Oñate, E. [2000] “An anisotropic elasto-plastic constitutive model for large strain analysis of fiber reinforced composite materials” Computer Methods in Applied Mechanics and Engineering, Vol.185, No 2-4, pp. 245-277 [4] CIMNE [2006] LessLoss Deriverable 49: Integration of knowledge on FRP retrofitted structures, LessLoss project (GOCE-CT-2003-505488), Sub-Project 7 [5] Corbella B, Vicens J, Costa J [2004]. Informe de los resultados de los ensayos de propagación de grieta en modo II del "round robin test" del proyecto GRINCOMP (mat2003-09768-c03). Technical Report 2004-CYT-01-IT01, CIMEP. [6] Lubliner J., Oliver, J., Oller S., Oñate, E. [1989] “A plastic damage model for concrete” International Journal of Solids and Structures, Vol. 25, No 3, pp. 299-326 [7] Martinez, X., Oller, S., Barbat, A., Rastellini, F. [2007] “New procedure to calculate compression strength of FRP using the Serial/Parallel mixing theory”, Accepted at 8th FRPRCS-8, Patras, Greece [8] Martinez, X., Oller, S., Barbat, A. [2006] “Numerical tool to study structural reinforcement of steel reinforced concrete (RC) structures under seismic loads using fibre reinforced polymers (FRP)”, ECEES-1 Proceedings, Geneva, Switzerland [9] Martinez X., Oller S., Barbero E. (2007). Study of delamination in composites by using the serial/parallel mixing theory and a damage formulation. Composites 2007, U. Porto. 12th to 14th Septemer/2007. [10] Oliver J., Cervera M., Oller S., Lubliner J. (1990). A Simple Damage Model For Concrete, Including Long Term Efects. Second International Conference on Computer Aided Analysis And Design of Concrete Structures. Vol. 2 pp. 945-958. Zell Am See, Austria. Viena. 1990. [11] Oller, S., Miquel, J., Zalamea, F. [2005] “Composite material behaviour using a homogenization double scale method” Journal of Engineering Mechanics, Vol.131, No 1, pp. 65-79 [12] Oller, S. [2002], Análisis y cálculo de estructuras de materiales compuestos, CIMNE, Barcelona, Spain [13] Oller, S., Oñate, E., Oliver, J., Lubliner, J. [1990] “Finite element non--linear analysis of concrete structures using a plastic-damage model” Engineering Fracture Mechanics, Vol. 35, No 1-3, pp. 219-231 [14] PLCd Manual [1991-now] PLCd. Non-linear thermomechanic finite element code oriented to PhD student education. Code developed at CIMNE [15] Rastellini, F. [2006] Numerical modelling of the constitutive non-linearity of composite laminates, PhD thesis, Departament de Resistència de Materials i Estructures a l’Enginyeria (RMEE) – UPC, Directors: S. Oller y E. Oñate. [16] Rastellini F., Salomon O., Oller S., Oñate E. (2006). Non-Linear mechanical damage modelling for long fibre-reinforced laminates. The e-Journal and Exhibition of Non-destructive Testing. Electronic Publication. www.ndt.net (CDCM 2006). [17] Sánchez-Palencia, E.S. [1987], Homogenization techniques for composite media, SpringerVerlag, Berlin, Germany. Chapter: “Boundary layers and edge effects in composites”, pp. 121-192 [18] Spadea, G., Benicardino F., Swamy, R. [1998] “Structural behaviour of composite RC beams with externally bonded CFRP” Journal of Composites for Construction, Vol.2, No 3, pp. 132-137 [19] Trusdell C., Toupín, R. [1960] The classical field theories, Handbuch der physic iii/I, Edition, Springer-Verlag, Berlin, Germany [20] Zienkiewicz, O. C., Taylor, L.R., [1991] The finite element method, McGraw-Hill, London, England