Sliding frame‐solid interaction using BEM/FEM coupling

In this work it is presented a coupling between the Boundary Ele‐ ment Method BEM and the Finite Element Method FEM for two‐ dimensional elastostatic analysis of frame‐solid interaction. The BEM is used to model the matrix while the reinforcement is modeled by the FEM. Regarding the coupling formulation a third degree polyno‐ mial is adopted to describe the displacement and rotations of the reinforcement, while a linear polynomial is used to describe the contact force among the domain matrix and the reinforcement. Perfect bounding contact forces are improved by means of redun‐ dant equations and Least squares method. Slip‐bounding with two and three paths written as function of relative displacement are used to calculate the transmitted contact forces. Examples are used to demonstrate that the proposed slip‐bounding procedure regularizes the contact force behavior.


INTRODUCTION
The combination of the finite element method FEM and boundary element method BEM to solve structural analysis is attractive because it allows for an optimal exploitation of the respective advantage of the methods Zienkiewicz et al., 1977 .The main strength of the BEM for boundary-value problems governed by linear, homogenous, and elliptic differential equations with constant coefficients is the reduction of the dimensionality of the problem by one unit for linear constitutive relations Brebbia, 1978;Brebbia, 1980 .Particularly, BEM is useful to model special situations such as very large or unbounded domains, geometrical singularities e. g. cracks or to obtain very accurate results in regions of complicated shape Aliabadi, 1997;Bonnet, 1999;Frangi et al., 2002 .Thus, coupling the BEM and the FEM allows exploiting their complementary advantages.By the other hand, the FEM is appropriate to solve a lot of problems, including e. g. those with heterogeneous or non-linear constitutive properties, or finite deformations.
The idea of combining these two methods goes back to Zienkiewicz et al. 1977 .One branch of BEM/FEM coupling is the iterative coupling in which the individual sub-domains are treated independently by either method.The procedure starts with an initial guess of the interface unknowns that will be improved by solving each sub-domain and returned to interface.This procedure repeats until an error tolerance is achieved.Although this iterative coupling is very attractive to software design, its convergence commonly depends on relaxation parameters which are rather empirical Estorff and Hagen, 2005 .For this reason, a direct coupling approach is adopted in this work.
Standard BEM formulations to deal with solids stiffened by bars or fibers are derived by combining the BEM and FEM algebraic equations.The domain continuum or matrix is analyzed by BEM, while finite elements are used to represent inclusions bars for instance .The coupling is always done by enforcing displacement compatibility and traction equilibrium at interface nodes.Practical applications for non-slipping or slipping coupling using BEM/FEM coupling are presented, for instance, in Beer and Watson 1996, Coda and Venturini 1995, Coda 2001, Leite et al. 2002, Botta and Venturini 2005, 2003, Leonel 2009, Rocha 2010 .As contribution, the present work is able to capture bending effects in the analysis of frame-solid interaction, which includes bending stiffness in reinforcements immersed in 2D continuum.For this purpose, in this study it is adopted a Rissner-Mindlin type frame finite element with third order of approximation for both displacement and rotations resulting in four nodes and 12 degrees of freedom.However, the contact force is modeled by linear approximation resulting in only 4 independent values which leads to a not square force matrix.These approximations displacements and tractions were used because it is the lowest order that satisfies the differential equation governing the problem.In this approach, the boundary element force lines are build in a compatible way, that is there are 4 source points generating a third order approximation for displacements, but the contact force approximation is linear, also resulting in a not square force matrix.
The least Squares technique is used to eliminate the dependent equations due to the above mentioned difference in approximation order for displacement and tractions.Moreover some authors, Botta and Venturini 2005, 2003and Leonel 2009 , claim that this procedure reduces contact force oscillations.
This paper is organized as follows.It is presented in section 2 the FEM formulation to model the frame structure, which is shown the kinematics adopted.In section 3 is shown the BEM formulation adapted to domain modeling.In section 4 it is presented the proposed coupling formulation between BEM and FEM considering both perfect bonding and debonding cases.This section is divided in subsection 4.1 in which there are presented the basic equations to perfect bonding and an example to verify this formulation.In sub-section 4.3 it is presented the coupling formulation considering the slip between reinforcement and domain.Debonding models, basic equations and the non-linear formulation to solve slipping are presented in sub-sections 4.3.1, 4.3.2 and 4.3.3,respectively.The sub-section 4.4 presents two examples.The first simulates a pullout test and the second solve a soil-structure interaction case.Finally, in section 5 the final remarks and conclusions are given.

THE FRAME ELMENT MODELING -FEM FORMULATION
As mentioned before, the FEM is used to model frame elements and structures.Here, elements which have three degrees of freedom per node and cubic approximation for displacement and rotation are employed.This way, the elements have four nodes and these nodes present two translations vertical and horizontal and one rotation.Moreover, the distributed applied forces will follow linear approximation.

Kinematics
For any point on the structure, the horizontal and vertical components of displacements are given by, respectively: where x and y are the reference system in the center of the layer, as shown in figure 1 and p U and p V are the displacements of point P .From equation 1 , one can apply the differential operator to obtain the linear strain components: Applying the constitutive law for the isotropic materials, the stress components at the point "p" are obtained: where E and G are the longitudinal and shear elastic moduli, respectively.
To write the equilibrium equation it is used the Principle of Minimum Total Potential Energy.Using equations 2 and 3 one writes the Total Potential Energy equation as, where e U and p U are the internal and potential energy of external forces, x t and y t are the components of the distributed loading contact tractions applied to the structure, L and A are the length and cross sectional area of the frame element, respectively.For approximate unknowns 0 U , 0 V e 0 q cubic independent approaches were used, as shown: q being the nodal values unknown .Since ( ) i j x and ( ) j y x are shape functions: Minimizing the energy functional, equation 4 , one finds the algebraic equilibrium system given by: where: E K is stiffness matrix, E G is equivalence force matrix, E U unknown vector with displace- ments and rotations, E f is the vector containing the nodal values of the distributed load, F is the concentrated force.Labels NF and extr NF are node numbers for the displacements and rotations , four per element, and node numbers for forces two per element , respectively

THE DOMAIN MODELING -BEM FORMULATION
Let us consider the domain W and its boundary G .For an elastic body defined by the domain W , the equilibrium equation, written in terms of displacements, is given by: ( ) where u represents the displacement vector, G is the shear modulus and n is the Poisson's ratio.
For a domain W with boundary G , the integral representation of displacements is derived by ap- plying reciprocity theorem or Green's second identity .with r being the distance between the source and field points, ( ) r are components of the vector r , n is the boundary normal unit vector and I is the second order identity tensor or Kronecker delta .
For numerical solution, u and p are approximate by polynomial functions over boundary ele- ments, in this work linear polynomials are used for boundary elements internal force lines are further discussed , and the integral equation 9 is converted into an equivalent algebraic system as follows: where matrix H is obtained from the left terms in equation 9 and matrix G from the terms on the right side.U is a vector which contains the nodal values of displacements for all boundary nodes and P is the nodal traction vector.
After substitution of the prescribed boundary conditions, the algebraic equations may be written as The vector X contains all the unknown boundary displacements and tractions, A is a coefficient ma- trix which is usually non-symmetric and densely populated, and B is a matrix which contains the coef- ficients corresponding to the prescribed boundary conditionsY .
One can differentiate equation 9 to derive the integral representation of strains and then apply the Hooke's law to obtain the stress integral equation, written for internal points, as follows, Where * S and * D are well known third-order tensor for the stress equation obtained by applying the Hooke's law on the fundamental solution at the source point Brebbia and Domingues, 1992 .

BEM/FEM COUPLING FORMULATION
In this section, additional terms inserted in the classical formulations of BEM and FEM as well as the coupling between BEM and FEM are shown.

Basic equations -perfect bonding
In this subsection the perfect bounding situation is described.It implies a direct compatibility between displacements and contact force equilibrium or continuity .Thus: Where R U and D U are vectors containing nodal displacements for frame element and domain, re- spectively; R f and D f are nodal distributed force vector applied on the frame finite elements and on the force line boundary element in the 2D domain , respectively.Once vector R U contains three com- ponents two translations and one rotation and vector D U contains two components two transla- tions , the correlation between these two vectors is done by Q matrix.
The components of the vectors , , For 3D solids the unknown forces, D f , would be distributed over internal surfaces.For 2D solids this forces appears distributed along internal lines, which, roughly, work as internal "boundaries" dedicated only for tractions.Thus, the integral equation 9 is modified to include the additional term: where D f is the internal force acting along the interface, G R , between the two materials and represents the fiber effect applied in the domain.Similarly, the integral equation 13 including this additional effect is written as: Selecting a proper number of collocation points source points at the boundary and at the frame element, two sets of algebraic equations are, respectively, written from equation 16 , as follows: where index b denotes boundary and index r denotes frame, respectively.Additionally, the first index denotes source points and the second index denotes field points, which, are on boundary or load line.
The superimposed frame element reinforcement and load line representations can be seen in  Therefore, the displacement compatibility is rewritten as: Where T is the matrix that relates the nodal positions and direction of the D U vector with the R U vector.Furthermore, the  T matrix has dimension 3NF columns by int 2N rows, for which int N is the amount of BEM interface points.
According to figure 3, the amount of finite element nodes is equal to the amount of boundary interface nodes, i. e.

NF N
= , therefore the determination of the coupling parameters and the boundary values can be summarized by the following system of equations: The 1 N values of b U and 2 N values of b P are known on 1 G and 2 G ( ) hence there are only 2 N unknowns in the system of equations 21 .As usual, to introduce these boundary conditions into 21 one has to rearrange the system by moving columns of bb H , rb H with bb G , rb G from one side to the other, respectively.Once all unknowns are passed to the left-hand side and applying conditions of equations 14 and 20 one can write the new system of equations as: Or in matrix form, As described above, at the fiber element level, the number of displacement values is larger than the number of bonding or contact force nodal values.This occurs because in equation 19 the number of algebraic relations is much larger than the number of force values in D f .To reduce the number of equations to the same as the number of unknowns one can apply the Least Square Method LSM .In this work the LSM is applied over equation 21 or 19 , as follows:

25
Where the matrix G G .Therefore, the system of equations 23 turns into a square system as:

Numerical example -perfect bonding
In this section, one numerical example is analyzed to check the performance and accuracy of the proposed BEM/FEM coupling for two-dimensional reinforced solids.
The reinforced simple supported beam subjected to homogeneous transversal surface load    The above results make evident that even the poorest BEM/FEM coupling discretization has good accuracy when compared with the reference result.
Regarding tractions at the interface, the BEM/FEM results can be seen in Figures 8 and 9.According to Figures 8 and 9, for perfect bonded elastic reinforcement, it's observed that for the three BEM/FEM discretizations there is a perturbation in the traction values at the ends of the fiber.The behavior of axial forces, as shown in Figure 8, had been observed by Botta and Venturini 2005 , however, the behavior shown in Figure 9 has not been reported before.Moreover, the use of LSM partially smoothes the behavior of contact force when compared with results that do not apply this strategy for the coupling, Leite et al. 2003 .At this point an important result of this study should be advanced; when debondig is allowed always occurs in a small vicinity of infinity contact stress the above perturbation disappears.
When using nonlinear constitutive relations for matrix, there are evidences of smooth solutions for this kind of problems which were reported by Coda 2001 .Moreover, one may note, in Figure 9, that the extension of the contact force perturbation reduces as discretization increases.

Coupling formulation -debonding
In this section, the previous coupling is extended to accomplish debonding.The adopted models to represent debonding as well as the nonlinear formulation to simulate the slip between domain and frame elements are presented.The feasibility of this formulation is shown through numerical examples.

Debonding models
Fibers embedded in the domain play an important role to improve solid stiffness and loading capacity if enough internal forces along the interface can be sustained.Sliding along the interface may be allowed when a certain amount of strength is preserved.The ideal situation for which perfect bonding is assumed, as shown in previous section, is impossible in practice; at least in the vicinity of fiber ends, as the interface forces approach the infinity Radtke et al, 2011 .Therefore, a certain amount of sliding occurs according to the bonding carrying capacity.
To model the slip that may occur in the fiber-domain interface, a debonding criterion should be considered.In this work, two models were implemented together with the proposed BEM-FEM coupling.
The curves shown in Figure 10 and 11 represent the debonding criterion that relates the bonding force f with the relative displacement at interface slip s .The following parameters define the two models: model 1 depends only on the maximum bonding force max f and model 2 depends on the maximum force max f , residual bonding force res f and slip characteristic values 1 s and 2 s .
From the Figure 10 and 11, the following relationships are written for the adopted models, respectively: Where vector S contains the nodal relative displacement values, T relates the nodal positions D U with S .Furthermore, T matrix has dimension 2 extr NF columns by int 2N rows.The other terms in equations 31 and 32 have already been explained, but now tractions on interface depend on relative displacements s .In this work, D U and R U have been approximated by cubic polynomial and S by linear polynomial.
From the introduction of relative displacements, the equilibrium equation 19 of the boundary element method may be rewritten as: Therefore, the BEM coupling equation can be rewritten as,: Applying boundary conditions, the equation 34 results: By the other hand, if S is known, one writes: ( ) As can be seen, in equations 36 and 37 , there are more equations than unknowns, since . Thus, to reduce the number of equations to be equal to the number of unknowns the Least Square Method LSM is applied over internal point equations 35 .This way: Where  rr G is the transpose matrix of rr G .Therefore, equations 36 and 37 become square, as follows: 39 Or: ( ) U in first and second equation, respectively, of 41 , one has: which can be replaced in equation 41 , resulting: Where I is the identity matrix.Linearizing equation 44 and using the first term of the Taylor's expansion, results: The derivative that appears in equation 47 is directly obtained from equation 44 using the debonding model relationships given by equations 27 -30 .Then, one has: The matrix é ù ë û CTO W , in equation 48 , is the consistent tangent operator of the proposed algorithm.
The derivatives on the right hand side of equation 48 depend on the updated slip value, computed appropriately according to the adopted model defined in equations 27 -.These derivatives are locally defined by: To model 1 Reaching the convergence in equation 44 for the time increment i t D after n iterations, one has to compute the slip variable s to start the next increment, as follows: Figure 14 shows the displacement curves, in meter, at the bar-domain interface.It is possible to verify a decrease on the slope of displacements for the points located near the load application, as the loading progresses.The evolution of decoupling is illustrated in figure 16 wherein the domain displacements are unequal to the bar displacements.At the first displacement increment it is verified that almost all nodes are perfectly coupled, except the ends nodes.As the free end displacement is increased, other nodes begin to decouple until the ninth increment situation in which the bar displacements continue increasing and domain displacement decreasing.One important aspect shown by this example, more evident in figure 13, is that the limitation of adherence force by the debonding process regularizes the shear contact force and other variables.

Example 2
The structure analyzed in this example is shown in Figure 17, it is important to note that this example cannot be solved without considering bending stiffness and transverse contact forces as did in this work.This structure is a 2D deep foundation, i.e., a pile embedded in an infinite soil.In order to simulate the presence of nearby structures and a rigid supporting rock mass, some soil displacement restrictions are imposed.The pile is considered inclined at an angle of 10  ( ) a .A soil region, with 5 me- ters length ( ) L and 20 meters depth ( ) H , is considered.The pile has 4 meters in length ( ) 0 L and a prescribed load of 0.11  Figure 17 shows the analyzed structure and the reference nodes for which results are presented.Figure 18 shows curves of the pile/soil interface forces at reference nodes.This figure also shows the development of those interface forces as the loading steps are increased for both slip's models.As can be seen in figures 19-21, the pile behavior for model 1 presents severe changes.These changes occur when the maximum load capacity is reached increment 35 ; then slip occurs without further gain of resistance.It is important to mention that over the load capacity the solution is unstable and the system loss objectivity.
For model 2 as the interface forces do not reach the residual part of the model the maximum load capacity is not achieved and no abrupt change occurs in the pile behavior.However, if the total load capacity is reached all contact points reach the residual part of the model the collapse occurs.

CONCLUSIONS
In this work a BEM/FEM coupling among frame bars and 2D continuum is successfully developed and implemented.The domain is modeled by BEM and reinforcement is modeled by FEM.The combination of the two methods is made by writing displacement compatibility and interface equilibrium.The most important feature of the formulation is the consideration of sliding at frame/continuum interface.This procedure is able to model, for example, the progressive failure of pile-soil interaction until reaching the collapse load, or the progressive failure of fiber reinforced bodies considering the influence of shear and normal contact forces.
Regarding the solution behavior an important conclusion should be stated.As boundary elements are able to model high stress concentrations, the values of contact forces at the beginning or ending of any perfect bounding region present a strong perturbation.The use of redundant algebraic equations and the least squares method are tested here and result in a small improvement of this phenomenon.The increasing of discretization reduces the extension of perturbation but increases the near singular stresses.
The complete solution for this problem results when using a more realistic model that allows the natural stress relaxation at singularities.The developed non linear behavior of contact forces, considering the sliding or decoupling between reinforcement and continuum, completely regularizes the contact force behavior, leading to reliable solutions for low or high load situations.Further developments are the consideration of non-linear behavior for both continuum and frame media.
the symbols "*" is used to indicate fundamental solution equation 10 and p represent bound- ary traction values.

Figure 2 :
Figure 2: Force and displacement approximation at frame interface.
figure 3 in which n elements finite and line are employed.The displacement nodes for each finite element are represented by squares and crosses represent the collocation points of equation 19 where displacements are calculated.It is observed that for fiber ends BEM equations are not written in the same position of the finite element nodes to avoid singularities; however displacements are extrapolated to the nodal positions in order to make them compatible.Moreover, this modeling allows the reinforcement frame ends to reach the body boundary without interfering in basic BEM equations.

Figure 4 :
Figure 4: Geometry of the structure under analysis.The beam is discretized by 120 linear boundary elements with the same length.Three discretizations are employed to model the reinforcement, with 25, 50 and 100 finite elements.The same number of force line BEM elements is employed to model the interface.The results are compared with the commercial software ANSYS employing 100 BEAM3 elements to model the reinforcement and 2800 PLANE42 2D solid elements to model the domain.Figures 5, 6 and 7compare results axial, transverse displacements and rotations achieved using the BEM/FEM coupling discretizations and ANSYS.
Figure 10: Deboding Model 1 see, equations 39 and 40 are non linear regarding slip s .To solve them, one has to take into account the non-linear relationship described by the debonding model presented in item 4.3.1, in which the relation between the debonding force f and the slip s is established.The equilibri- um equation 39 is then rewritten in terms of the variable increments, as follows:

Figure 14 :
Figure 14: Evolution of the domain displacements along the interface, BEM/FEM.

Figure 15
Figure 15 shows the relative displacements between bar and domain.According to the definition given in equation 32 , the bar displacements are obtained by subtracting the results of figure 14 from figure 15.The evolution of decoupling is illustrated in figure16wherein the domain displacements are unequal to the bar displacements.At the first displacement increment it is verified that almost all nodes are perfectly coupled, except the ends nodes.As the free end displacement is increased, other nodes begin to decouple until the ninth increment situation in which the bar displacements continue increasing and domain displacement decreasing.

Figure 15 :
Figure 15: Relative displacement s . Loads increase from zero to the reference values in 50 equal increments.A mesh of 300 linear elements is used to model the soil and 100 finite elements line elements are used to model the pile.The soil elastic properties are: .For model 1 it is adopted: max 2
Figures 19-21 show the displacement in the axial and transverse directions, as well as the rotation along the pile for each reference node.

Figure 19 :
Figure 19: Axial displacement at the interface pile/soil.

Figure 20 :
Figure 20: Transverse displacement at the interface pile/soil.

Figure 22
Figure 22 a and b show the displacements in the y direction, in meters, for the soil internal points considering the models 1 and 2, respectively.Figures 23 a and b show the stress values, y s ,

Figure 22 :
Figure 22: Soil displacement in the y direction

Figure 23 :
Figure 23: Soil stress, y s , for different adherence models.
F.C. Rocha, W. S. Venturini and H.B. Coda / Sliding frame-solid interaction using BEM/FEM coupling 1391 Latin American Journal of Solids and Structures 11 2014 1376-1399Equation 44 represents a non-linear system of equations given in terms of the slip increment It can be solved by applying the iterative Newton-Raphson scheme.Then, from the iteration n