Fiber-matrix Contact Stress Analysis for Elastic 2 D Composite Solids 1

This paper presents a finite element formulation for the analysis of two dimensional reinforced elastic solids developing both small and large deformations without increasing the number of degrees of freedom. Fibers are spread inside the domain without the necessity of node coincidence. Contact stress analysis is carried out for both straight and curved elements via two different strategies. The first employs consistent differential relations and the second adopts a simple average calculation. The development of all equations is described along the paper. Numerical examples are employed to demonstrate the behavior of the proposed methodology and to compare the contact stress results for both calculations.


Fiber-matrix Contact Stress Analysis for Elastic 2D Composite Solids 1 INTRODUCTION
Composites are made of more than one material in order to take advantage of complementary characteristics.Three widely used composite materials that falls in the fiber reinforced category are the reinforced concrete, the reinforced rubber and the fiber-carbon-epoxy.The reinforced concrete combines the low cost of concrete and its strength to compression with the ductility and strength to traction of the steel.The reinforced rubber combines the large deformability of rubber for dynamic energy absorbing with the steel strength and stiffness.The fiber-carbon composite uses the matrix as a bounding among the carbon fibers.Mechanical analysis of fiber-reinforced composites falls in three main levels: the macro-level, the meso-level and the nano-level.The first is interested in the overall behavior of structural components.The meso-level deals with the interdependent behavior of fiber and matrix, i.e. interfacial stress and slip.Finally the nano-level is interested in the nano-scale constitution of fibers and matrix by themselves and in its influence at meso and macro-levels.
In this paper we intend to collaborate in meso and macro-levels of fiber reinforced modeling via Finite Elements.We propose a way to represent short or long fibers immersed in elastic continuum Latin American Journal of Solids and Structures 12 (2015) 583-611 domains by means of fiber finite elements without increasing the number of degrees of freedom and a way to calculate, with good precision, the contact stresses (shear and normal) between fiber and matrix without using auxiliary bounding layer strategies.
In FEM literature different ways are developed to incorporate fibers inside matrix domain.The reader is invited to consult the works (Radtke et al., 2011(Radtke et al., , 2010a(Radtke et al., , 2010b;;Hettich et al., 2008;Chudoba et al., 2009;Oliver et al., 2008) where field enrichments are imposed inside the 2D domain in order to model the fiber-matrix coupling.These enrichments are based on general known behavior of fiber-matrix connections.They are mostly based on the so called Partition of Unity FEM (Melenk and Babuska, 2006;Duarte andOden, 1996a, 1996b;Oden et al., 1998;Duarte et al., 2000;Babuska and Melenk, 1997).These works are very elegant and well posed; however the pre-known enrichment field is of difficult achievement when, for example, curved fibers are present.Readers are invited to consult the works Schlangen et al. (1992); Bolander and Saito (1997); Liz et al., (2006) in which authors employ lattice strategy to model composites from micro-structures.Moreover, other approaches that adopt slip degrees of freedom to represent the fiber reinforced body can be found in Balakrishnan and Murray (1986); Désir et al. (1999).Works that employ the Boundary Element Method should also be mentioned (Coda, 2001;Leite et al., 2003).
This study presents an alternative geometrically nonlinear formulation to analyze 2D solids reinforced by fibers.The 2D solid finite element applied here to discretize the continuum is isoparametric of any order (Coda, 2009;Coda and Paccola, 2008;Pascon and Coda, 2013).Curved high order fiber elements are developed to be embedded in the continuum.To calculated contact stresses without using slip degrees of freedom two approaches are developed and compared, an average stress calculated from the transfer fiber force resultants and a differential relation among normal fiber internal stress and contact forces.The adopted nodal parameters are positions, not displacements, which is adequate to model curved elements and large deformations due to the natural presence of a numerical chain rule.The formulation is classified as total Lagrangian and the Saint-Venant-Kirchhoff constitutive law is chosen to model the material behavior (Ciarlet, 1993;Ogden, 1984).Therefore, the Green strain and the Second Piola-Kirchhoff stress are adopted.
Fiber elements are introduced in matrix by means of nodal kinematic relations.This strategy directly ensures the adhesion of fibers nodes to the matrix without increasing the number of degrees of freedom and without the need of nodal matching (Sampaio et al., 2013;Vanalli et al., 2008).To solve the resulting geometrical nonlinear problem we adopt the Principle of Stationary Total Potential Energy (Tauchert, 1974).From this principle we find the nonlinear equilibrium equations and the Newton-Raphson iterative procedure (Luenberger, 1989) is used to solve the nonlinear system.External loads are considered conservative and incrementally applied.
The paper is organized as follows.Section 2 describes the general nonlinear solution process, indicating the important variables that will be developed in subsequent sections.Section 3 describes the procedure used to model the two-dimensional continuum.Section 4 presents the developed any order fiber finite element and describes the chain rule applied to generalize the inclusion of fibers into high order 2D solid finite elements without increasing the number of degrees of freedom.Section 5 presents the proposed fiber-matrix contact stresses calculations.Section 6 presents the numerical examples comparing and analyzing the behavior of the proposed formulations.Finally, conclusions are presented in section 7. Structures 12 (2015) 583-611

THE NONLINEAR SOLUTION
In this section, the strategy adopted to solve the reinforced 2D solid geometrically nonlinear equilibrium is described.It is important to clarify that next sections explain how the fiber degrees of freedom are related to the main unknown of the process, i.e., the 2D solid nodal positions without increasing the number of degrees of freedom.
The nonlinear analysis starts writing the total potential energy ( Π ) as follows: where U is the strain energy including matrix and fiber contributions written regarding solid nodal positions and Ω is the potential energy of external conservative applied forces given by: where j F is the vector of external forces and j Y is the current position vector.The Principle of Stationary Total Potential Energy (Tauchert, 1974) is applied writing the equilibrium equations as the derivative of total energy regarding nodal positions (2D solid for instance), as: where int j F is the internal force vector or the strain energy gradient vector calculated regarding solid nodal positions.The nodal current positions are the unknowns of the problem, so, when adopting a trial position in Eq. (3) j g is not null and becomes the unbalanced force vector of the Newton-Raphson (Luenberger, 1989) procedure.Expanding the unbalanced force vector around a trial solution 0 Y , one has: which can be rewritten, neglecting higher order terms as: ( ) ( ) where is the correction of position and is the Hessian matrix or tangent stiffness matrix.
The trial solution is improved by: Latin American Journal of Solids and Structures until k Y ∆ or j g become sufficiently small applied if one wants to describe the equilibrium path of the analyzed structure.

ISOPARAMETRIC 2D SOLID FINITE ELEMENT
As we are interested in composite analysis the procedure described in continuum part of the composite, i.e., the matrix.After achieving the strain energy of the matrix the next section is concerned with the introduction of the fiber strain energy into the mechanical system.

Kinematical approximation and positional mapping
By means of the illustration of a quadratic finite element, Figure 1 shows the 2D solid (matrix) mapping from the initial configuration (not deformed) al., 2000).This mapping is done using a dimensionless auxiliary configuration The initial configuration 0 B whose points have coordinates space 1 B with coordinates i ξ using shape functions of any order, of the nodes l in the initial configuration, Similarly, the current configuration B pression: become sufficiently small (Luenberger, 1989).The load level can be incrementally applied if one wants to describe the equilibrium path of the analyzed structure.

ISOPARAMETRIC 2D SOLID FINITE ELEMENT -ISOTROPIC CONTINUUM
As we are interested in composite analysis the procedure described in this section is applied to the continuum part of the composite, i.e., the matrix.After achieving the strain energy of the matrix the next section is concerned with the introduction of the fiber strain energy into the mechanical pproximation and positional mapping By means of the illustration of a quadratic finite element, Figure 1 shows the 2D solid (matrix) mapping from the initial configuration (not deformed) 0 B to its current configuration B (Bonet et .This mapping is done using a dimensionless auxiliary configuration 1 B .
Initial and current configuration mappings.
whose points have coordinates i x is mapped from the dimensionless using shape functions of any order, ( ) , l φ ξ ξ , and by the coordinates in the initial configuration, l i X , such as: B is mapped from the dimensionless space 1 B by the e ( ) . The load level can be incrementally this section is applied to the continuum part of the composite, i.e., the matrix.After achieving the strain energy of the matrix the next section is concerned with the introduction of the fiber strain energy into the mechanical By means of the illustration of a quadratic finite element, Figure 1 shows the 2D solid (matrix) Bonet et is mapped from the dimensionless , and by the coordinates correspond to coordinate directions.The deformation function f that maps the initial configuration 0 B to the current configuration B can be written as a composition of mappings 0 f and 1 f as: ( ) The deformation gradient A can be derived directly from 0 A and 1 A as (Bonet et al., 2000;Coda and Paccola, 2007): Equation ( 10) can be understood as a numerical chain rule because the initial mapping gradient 0 A is a known numerical quantity.The solid element has any order of approximation and N is the number of nodes given as a function of the approximation order GP as:

Continuum strain energy
Without loss of generality, to simulate the continuum portion of the composite (matrix) we adopted the Saint-Venant-Kirchhoff specific strain energy function (Ciarlet, 1993;Ogden, 1984), as: where ijkl C is the elastic fourth-order tensor and E is the Green-Lagrange second-order strain expressed respectively by: ( ) The variables ⋅ t C = A A and δ δ δ δ are the right Cauchy-Green stretch tensor and the Kroenecker delta, respectively.In Eq. ( 13), G is the shear modulus and ν is the Poisson's ratio.
The strain energy accumulated in the continuum part of the composite (matrix) is calculated by integrating the specific strain energy over the initial volume, i.e., Considering solids with unitary thickness and writing Eq. ( 15) as a function of dimensionless coordinates ( 1 ξ and 2 ξ ) results: where 0 J is the Jacobian of the initial mapping, i.e., with 0 A given by Eq. ( 10).The matrix strain energy can be derived directly regarding the solid finite element positions finding the conjugate internal force vector, as: in which α is the direction and β is the node.The derivative inside the integral term of Eq. ( 18) can be developed as: where E is the second Piola-Kirchhoff stress tensor.From the definition of the right Cauchy stretch and Eq. ( 10) one writes: From Eq. ( 8) results To complete the necessary variables of the solution process (section 2) it is necessary to calculate the second derivative of strain energy regarding nodal positions, resulting into the Hessian matrix as: in which Latin American Journal of Solids and Structures 12 (2015) 583-611 and It should be noted that, for solid elements, the second derivative of 1 A regarding nodal parameters is null, simplifying expression (24) to: Next section presents the necessary expressions to introduce fibers into the composite formulation.

ELASTIC FIBER REINFORCEMENT -KINEMATICS AND ENERGY CONSIDERATIONS
This section is divided into two subsections.The first describes the strain energy of general curved fibers and the second describes the strategy used to introduce fiber energy in the composite solution without increasing the number of degrees of freedom.

Any order curved fiber element
To guaranty total adherent fiber-matrix coupling when using high order solid elements it is also necessary to adopt high order fiber elements (Sampaio et al., 2013).Figure 2 shows the nondeformed initial configuration 0 B , the current configuration B and a non-dimensional auxiliary configuration 1 B for the curved fiber finite element of any order.The initial configuration 0 B whose points have coordinates i x is mapped from the dimensionless space 1 B with coordinates ξ using shape functions of any order, ( ) P φ ξ , and by the coordinates of nodes P , P i X , at the initial configuration, such as: ( ) Latin American Journal of Solids and Structures The current configuration B is mapped from the dimensionless space where i y are the coordinates of points in the current configuration ordinates of fiber nodes.In Eqs. ( 26) and ( 27) index ly, the fiber finite element nodes and coordinate directions.
The tangent vector of the fiber and its modulus are calculated at the initial configuration as 0 ( ) It is important to mention that 0 B T ration one finds: ( ) Using Eqs. ( 28) and ( 29) one write the one or in its expanded form, 1 2 is mapped from the dimensionless space 1 B by: Mapping of the fiber finite element -initial and current configurations.
( ) are the coordinates of points in the current configuration B and P i Y are the current c ordinates of fiber nodes.In Eqs. ( 26) and ( 27) index 1,..., P N = and 1, 2 i = represent, respectiv ly, the fiber finite element nodes and coordinate directions.
The tangent vector of the fiber and its modulus are calculated at the initial configuration as and 0 2 2 2 1 2 ( ) ( ) is the differential Jacobian of 0 i f .For the current config and (28) and ( 29) one write the one-dimensional Green strain as . For the current configu- Using the Saint-Venant-Kirchhoff constitutive law one writes the specific strain energy at a point of the fiber as: where E is the elastic modulus and ( ) E ξ is the Green strain measure defined by Eq. ( 30) or (31).The strain energy of a curved fiber is given by integrating equation (32) over its initial volume 0 V as: In order to proceed with the equilibrium analysis it is necessary to know the first derivative of strain energy regarding positions.Based on the energy conjugate concept the natural internal fiber force vector related to node j and direction k ( ) is calculated regarding fiber parameters as: From Eqs. ( 31) and ( 32) follows in which S is the one-dimensional Second Piola-Kirchhoff stress.Considering S constant over the cross section area A of the fiber one transforms 0 f dV into a simple expression, resulting: Where, as mentioned before, The Hessian matrix components for the fiber element are obtained by the second derivative of the strain energy, i.e.: Latin American Journal of Solids and Structures 12 (2015) 583-611 Developing the necessary calculations one achieves: or Integrals ( 36) and ( 40) are solved using Gauss-Legendre quadrature.
The parameters used to find the internal force and Hessian matrix for fibers are not suitable to be directly applied in the solid solution process; in next section the necessary transformations are presented.

Coupling strategy -kinematical fiber matrix coupling
The procedure adopted here to embed fibers at any position of the domain without increasing the number of degrees of freedom is an extension of the works of Vanalli et al. (2008) where linear elements and linear elasticity were adopted.

General fiber / solid connection
The requisite to start the procedure to embed curved fibers in curved solid elements is to know the solid dimensionless coordinates related to fiber nodes coordinates.This is done solving the pair of dimensionless solid variable 1 2 ( , ) p p ξ ξ associated to the physical fiber node position in the following nonlinear system, 1 2 ( , ) where l φ are the shape functions of the solid element, P i X are the known physical coordinates of fiber nodes (generated independently of solid mesh) and l i X are the know solid nodes coordinates.To solve Eq. ( 41) one expands it in Taylor series until the first order and starts with a trial dimensionless coordinate, 1 2 ( , ) pt pt ξ ξ , i.e.: ) Latin American Journal of Solids and Structures 12 (2015) 583-611 in which Pt i X is a trial position of the fiber node calculated from the solid element geometry and the trial dimensionless coordinates and ij H is a two dimensional matrix.The correction of the trial dimensionless coordinates i ξ ∆ is calculated solving the following linear system of equation: The procedure is a simple and fast Newton-Raphson nonlinear solver that relates all fiber nodes to the connected solid element revealing the pair of dimensionless variables 1 2 ( , ) p p ξ ξ .From this information one also knows the current position of fiber nodes as a function of solid nodes positions, i.e., 1 2 ( , ) in which l i Y are the current positions of solid nodes.Equation ( 44) ensures the connection among nodes of fibers and the matrix 2D elements.
In the next item it will be necessary to differentiate the fiber strain energy regarding the solid nodal coordinates using the chain rule.To make it possible one has to differentiate Eq. ( 44) regarding a generic nodal solid coordinate, as: If the fiber node P belongs to the solid element then 1 l β δ = and if direction α (solid) is equal to direction i (fiber) then ( , ) , otherwise it results zero.

General internal force
The strain energy stored in a reinforced body is the sum of the strain energies stored in the matrix and fibers, where mat U is the strain energy stored in the 2D solid finite elements used to discretize the matrix and f U is the strain energy stored in the fiber finite elements.Therefore, the internal force at a node β following direction α , considering both fiber and matrix contributions is found by the con- jugate energy concept, such as: where Eqs. ( 18), ( 36) and ( 45) have been used and there is no summation over P .

Hessian Matrix
Proceeding as described for the calculation of internal forces, we develop the second derivative of strain energy of the reinforced finite element regarding the solid nodal parameters, as follows ( ) The first integral at the last term of Eq. ( 48) is known, Eq. ( 22).However, it is necessary to observe that the kernel of the last integral is the specific strain energy of a fiber derived twice regarding the solid nodal parameters.As Eq. ( 39) gives its value when derived regarding fiber parameters one has to apply the chain rule twice, described by Eq. ( 45), over Eq. ( 48), that is: where f h is the fiber Hessian matrix kernel, Eq. ( 38).In Eq. ( 50) summation is not implyied.Integrating (50) over fiber volume gives: Using Eq. ( 51) into (48) results the consistent spreading of fibers contribution over the matrix properties, i.e.:

CONTACT STRESS CALCULATION
In the proposed formulation the contact stresses are not used to achieve the equilibrium position.Therefore, the contact stress calculation is done after solving the problem for positions.Two strategies are presented in order to do this calculation.The first is based on differential relations and is presented in subsections 5.1 and 5.2.The second is an average calculation, shown at subsection 5.3, for which the transferred nodal force from fiber to matrix is decomposed following the tangential and normal directions of fiber node and divided by a influence portion of the element contact area.

Shear contact stress -differential formulation
The differential equilibrium depicted in Figure 3 is used to calculate the contact shear stress ( ) q ξ between fiber and matrix as: where ds is a differential of the curvilinear coordinate along the deformed fiber and t is the fiber thickness.ds is calculated as a function of dimensionless coordinate ξ as: Using the chain rule results: Considering constant thickness the Cauchy stress is calculated from Piola Kirchhoff stress as: With ( ) E ξ given by equation ( 31).Instead of differentiate expression (58) regarding ξ and substituting into ( 56) and ( 53) we prefer to calculate the normal force nodal values using ( 31) and ( 58) and make: Substituting ( 60) into ( 56) and ( 56) into (53) results: It is important to mention that the differential formulation cannot be applied for linear fiber element, as it results null values of ( ) q ξ .

Normal contact stress -differential formulation
For an infinitesimal length ds there is a curvature center and a curvature radius, see Figure 4.
From this figure one writes the following geometrical relation: Remembering that t is the body thickness, the equilibrium equation following the fiber orthogonal direction is given by: . . ..( .) 2. . . 2 58), it is necessary to calculate 1 ( ) R ξ , given by: ( ) Substituting ( 65) into (64) results the final expression: Latin American Journal of Solids and Structures 12 (2015) 583-611 ( ) By the same reason described in the previous item, expression (66) cannot be applied for linear fiber elements.In section 6 examples are used to demonstrate the use of presented expressions.

Average contact stress
A generic nodal transfer force (from matrix to fiber) int F is depicted in Figure 5 together with the tangential and normal unit vectors ( n and t ) calculated at the same point.The tangential and normal components of the transfer force are: The normal and tangential forces are divided by the surface influence area ( inf A ) depicted in Figure 5 and the average values results: It is worth noting that for connecting nodes the resulting value is the average between values calculated for each element.

NUMERICAL EXAMPLES
Five numerical examples are chosen to show the behavior of the proposed formulation regarding the overall behavior of reinforced structural members and the estimated contact stress accuracy.A convergence analysis is carried out for both displacement and contact stresses.

Simple supported beam
This example is used to certify that the mechanical coupling between fiber and matrix is working Latin American Journal of Solids and Structures 12 (2015) 583-611 properly.Both displacement and contact stress are compared to a simple analytical solution limited to small displacements and strains in order to certify the coherence of results.
It is a simple supported reinforced beam, depicted in Figure 6, subjected to a uniformly distributed load 10 N cm q = . The simple supports are modeled by vertical distributed loads of 100 N/cm r q = and the adopted geometrical properties are: , see Figure 6.The Young modulus and the Poisson´s ratio of the matrix are  In order to check displacement convergence, see Figure 7, we use 160 linear fiber elements to discretize the horizontal reinforcement and three different meshes of triangular third order finite elements, i.e., mesh (a) 8x80 elements and 12050 dof, (b) 16x160 elements and 47138 dof and (c) 32x320 elements and 186434 dof.As one can see the displacement difference from the second to the third discretization is less than 0.3% characterizing convergence.Adopting the mid position of reaction as the span position one achieves the technical reference value ( ) ℓ obviously smaller than the achieved numerical value.The convergence analysis for shear contact stress is made using linear fibers and the average technique, see Figure 8.We adopt 80 and 160 linear fibers for meshes (a), (b) and (c).As one can see there is no significant difference in results characterizing convergence.However, it is important to note that when one increases the number of fiber finite elements it is also necessary to increase the continuum mesh.Adopting the average calculation formulation and using mesh (b), Figure 9 compares the shear contact stress for equally spaced nodes using 160 linear, 80 quadratic and 53 third order fiber approximations.Figure 10 shows the same results of Figure 9 when the differential procedure is adopted to calculate the shear stress for second and third order approximations (remembering that the differential procedure is not applied for first order approximation).From Figures 9 and 10 one concludes that when using high order fiber elements a perturbation of the shear contact stress appears; moreover the use of differential formulation, expected to be the most precise one, introduces even more oscillations.At the beginning we did not know why this spurious behavior appears, however after thinking over the subject we conclude that, due to the nodal characteristic of internal force transfer from matrix to fibers a not recommended finite element operation is taking place.This not recommended operation is the application of concentrated forces at central points of higher order finite elements which leads to well-known spurious displacements and stiffness distributions.To illustrate the reason of this perturbation the horizontal displacement of node B of the simple example depicted in Figure 11 is solved using two first order bar finite elements and one second order bar finite element.Using two linear elements one achieves exactly the analytical result, i.e., ( ) , however when using one second order element the result is ( ) , a wrong solution that reveals the inappropriate force transfer and stiffness distribution of high order fiber elements immersed in the continuum matrix.
Therefore, one concludes that, although achieving a better adherence between fibers and matrix when using high order fiber elements (as demonstrated by Sampaio et al. (2013)) the contact stress distribution is better described by simple linear fiber elements.Moreover, using a good discretization for linear fiber elements an adequate adherence is achieved, see Sampaio et al. (2013).In order to be complete, Figure 12 shows the stress distribution in the continuum matrix using mesh (c) and 160 linear fiber elements.Figure 13 shows the normal stress along the reinforcement.developed contact stress.Due to symmetry only a half of the problem is discretized, as depicted in Figure 14b.The adopted matrix mesh ( 8x240 ) is chosen after a convergence analysis.For 2kN P = Figure 15 compares the analytical shear stress and the ones achieved using the average technique for 240 and 120 linear fiber elements.Obviously that the analytical solution is only a reference value as the Euler-Bernoulli hypothesis excessively simplifies the problem.In Figure 16 one finds the same result using 240 quadratic and 480 cubic elements, revealing the same spurious behavior detected in the previous example for high order elements.Figure 17 shows the normal contact stress for the same load level using linear fiber elements.As expected the normal contact stress increases near the loaded region.In Figure 18 we multiply by 100 the result of Figure 15 and compare it to the shear contact stress achieved when effectively using 200 kN P = . This procedure indicates the influence of large displacements in results.Figure 19 shows the deformed configuration for 200 kN P = with the Cauchy stress distribution inside the matrix.

Curved reinforcement
In this example the same beam modeled in example 1 is solved considering the curved reinforcement as depicted in Figure 20.The matrix adopted discretization is mesh (b) 16x160 elements and 47138 dof of example 1.In Figure 21 the horizontal component of the contact stresses is depicted using 155 linear fiber elements along the straight part of the reinforcement and 12 linear fiber elements along the curved part.These values are compared to the shear contact stress of example 1.The difference in position is obvious as the curved bar is longer than the horizontal one.From this comparison one may note that when using curved reinforcement a faster transfer of the reinforcement normal force to the continuum occurs.Moreover, a change of sign of this transfer appears at the free extremity of the curve, corresponding to a normal contact stress among fiber and matrix, see Figure 22.In Figure 23 the vertical component of the contact stresses is depicted.As it is expected, integrating the area above and under the zero line of this graphic results zero.

Stretch of a reinforced bar (long fiber)
This example is the stretch of a reinforced bar as depicted in Figure 24.The length of the bar is , the height is 20 cm h = and the thickness is 1 cm b = . The reinforcement has area and is placed at the center of the cross section distant 20 cm from the bar extremities.Due to symmetry a half of the problem is solved, see Figure 24.The matrix property is the same as the adopted in the first example.The reinforcement elastic modulus is varied from  The adopted matrix discretization is 32x320 , mesh (c) of example 1.In Figure 25 the contact shear stress is evaluated for 5 2 210x10 N cm f E = and using 80, 160 and 320 linear fiber elements.Due to the previous results the average strategy is assumed.As expected the elastic problem is mesh dependent because the shear stress presents unbound value at the extremity of the bar.
Figure 26 shows the maximum displacement of the matrix as a function of the number of fiber elements, the relative difference among the last two values is less than 0.05% which characterizes that overall displacement presents low sensitivity to the discretization of fiber.In Figure 27 we present the shear contact stress behavior regarding the elastic modulus of the reinforcement.The behavior of the shear stress does not present a proportionality regarding the stiffness of reinforcement, but becomes more singular, i.e., presents a more pronounced maximum at extremities as the elastic modulus difference between fiber and matrix grows.Figure 28 illustrates the reinforcement influence on the matrix 11 σ stress behavior.Almost all stress is transferred to the bar along less than 25 cm.It is important to mention that the limits adopted in legend of Figure 28 are taken on purpose to show the transition that is hidden if the total stress range is assumed.The matrix horizontal stress distributions for different fiber elastic modulus are presented in Figure 30.The same maximum and minimum values have been adopted to facilitate comparisons.As the reinforcement elastic modulus grows the average matrix stress reduces.It is obvious that at the extremities of reinforcements localized high stress appears, as depicted in Figure 31.

CONCLUSIONS
In this study an alternative methodology to analyze elastic reinforced solids, fiber matrix coupling, including the calculation of contact stresses is proposed.The procedure is applied in 2D domains and does not increase the number of degrees of freedom of the original 2D mesh.Two methodologies to calculate the contact stress between fiber and matrix are proposed and implemented.The first is applicable to high order fiber elements and is based on differential relations among the reinforcement (or fiber) internal normal stress and the contact stresses, tangential and normal.The other methodology is based on the division of the internal transfer force by a contact influence area.The last is also applicable to straight linear fiber elements.
From the examples we conclude that the overall displacement behavior is almost insensible to both order and number of the adopted fiber finite element discretization.However, an unexpected surprise reveals that the contact stresses for high order elements present spurious behavior for both average and differential calculations.Moreover no convergence is achieved if the number of fiber elements is increased.We concluded that this unexpected and undesirable behavior is due, though hidden, to the improper application of concentrated forces on intermediate node of high order finite elements used to discretize reinforcements.The application of random short fibers is also successfully tested.
Latin American Journal of Solids and Structures 12 (2015) 583-611 In future works we intend to apply the successful developments of this paper to consider slip among fiber and matrix without increasing the number of degrees of freedom.

Figure 1 :
Figure 1: Initial and current configuration mappings.

Figure 2 :
Figure 2: Mapping of the fiber fin stress analysis for elastic 2D composite solids Latin American Journal of Solids and Structures 12 (2015) 583-611

Figure 4 :
Figure 4: Infinitesimal part of a fiber finite element.

Figure 5 :
Figure 5: Transfer force and area of influence.
while the Young modulus and the cross-sectional area of the reinsymmetry, only a half of the problem is solved as depicted by Figure6.

Figure 9 :
Figure 9: Contact shear stress for high order elements (average calculation).

Figure 10 :
Figure 10: Contact shear stress for high order elements (differential calculation).

Figure 11 :
Figure 11: Simple example to reveal the spurious contact shear stress behavior of high order reinforcement finite elements.

Figure 12 :
Figure 12: Matrix stress distribution in

Figure 13 :
Figure 13: Normal stress along the reinforcement.

Figure 18 :Figure 19 :
Figure 18: Large displacement influence in contact shear stress distribution.

Figure 21 :
Figure 21: Horizontal component of contact stresses.

Figure 23 :
Figure 23: Vertical component of contact stresses.
N/cm is applied at the extremities.