Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 55, 3, pp. 897-910, Warsaw 2017 DOI: 10.15632/jtam-pl.55.3.897 MODELING AND SIMULATIONS FOR QUASISTATIC FRICTIONAL CONTACT OF A LINEAR 2D BAR Mikael Barboteu, Nacera Djehaf Laboratoire de Mathématiques et Physique, University of Perpignan Via Domitia, Perpignan, France e-mail: barboteu@univ-perp.fr; nacera.djehaf@univ-perp.fr Meir Shillor Department of Mathematics and Statistics, Rochester, MI, USA; e-mail: shillor@oakland.edu Mircea Sofonea Laboratoire de Mathématiques et Physique, University of Perpignan Via Domitia, Perpignan, France e-mail: sofonea@univ-perp.fr This work considers amathematical model that describes quasistatic evolution of an elastic 2D bar that may come in frictional contact with a deformable foundation. We present the model and some of its underlying assumptions. In particular, the novelty in the model is that both vertical and horizontal motions are taken into account, whichmakes it especially usefulwhen frictional contact is concerned.Contact is describedwith the normal compliance conditionand frictionwith theCoulomblawofdry friction.We introduceahybridvariational formulation of the problem and a numerical discretization based on a uniform time step and the finite element method in space. The numerical algorithm has been implemented, and we present computer simulations that illustrate themechanical behavior of the systemwith emphasis on frictional aspects of the problem. Keywords: new 2Dbar, contact, friction, computational results, finite element discretization 1. Introduction Aconsiderable progress has been achieved in the last two decades in themodeling,mathematical analysis and numerical simulations of various processes involved in the mechanical contact of solids. As a result, Mathematical Theory of Contact Mechanics (MTCM) has been reaching a state of maturity. The theory is concerned with mathematical structures that underly the modeling of general contact processes with different constitutive laws, i.e., different materials, different possible geometries and different contact conditions, see for instance (Eck et al., 2005; Han andSofonea, 2002;Migórski et al., 2013; Panagiotopoulos, 1993; Shillor et al., 2004; Sofonea and Matei, 2012) and the many references therein. MTCM aims to provide a sound, clear and rigorous framework to construct models for processes involved in contact, and necessary tools and ideas to prove the existence, possible uniqueness and regularity results for the solutions of thesemodels.Moreover, the variational formulation of thesemodels leads directly and naturally to sophisticated numerical methods for computer approximations of the solutions. The interest in contact problems involving thin structures such as bars, rods, beams and plates lies in the fact that their mathematical analysis is simpler as it avoids some of the com- plications arising in 3D settings and often provide considerable insight into the possible types of behavior of the solutions, i.e., behavior of such structures. Often, there is decoupling of some of the equations, which simplifies the analysis and computer simulations. Moreover, one may use such models as tests and benchmarks for computer schemes meant for simulations of more complicated multidimensional contact problems. Models, analysis and computer simulations of 898 M. Barboteu et al. various contact problems for rods and beams can be found in (Ahn et al., 2012; Andrews et al., 2012; Kuttler et al., 2001; Shillor et al., 2001) and the references therein. Amathematical model that describes unilateral contact of a beam between two deformable obstacles was studied in (Barboteu et al., 2012b). In this paper, we introduce a new mathematical model that describes quasistatic frictional contact between a 2Dbar and a foundation.The fullmodel derivation and analysis can be found in (Sofonea and Shillor, 2017). We assume that the foundation is deformable and wemodel the contactwith thenormal compliance contact condition, and friction is describedby the associated Coulomb law. Themodel is two-dimensional and its main unknowns are vertical and horizontal displacement fields, both defined on an arbitrary section of the bar. We state the variational formulation of the model, which includes a set-inclusion related to the friction condition, then present a numerical approach. Finally, we provide numerical simulations that illustrate the mechanical behavior of the solution of the quasistatic frictional contact model. In particular, we study the behavior of the numerical solution with respect to the stiffness coefficient of the normal compliance law. This study clearly shows that the problemwith a unilateral constraint, in which the obstacle is assumed to be completely rigid, may be approached as closely as one wishes by the solution to the contact problemwith normal compliance, with a sufficiently large stiffness coefficient. The rest of paper is structured as follows. In Section 2, we describe our model. Section 3 introduces the variational formulation of the problem, and a fully discrete variational approxi- mation by considering a hybrid formulation. Section 4 describes a special 2D finite rectangular element used in the discretization of the 2D bar. In Section 5, we present numerical results on the contact of the bar with a planar or circular foundation. Finally, in Section 6, we provide a very short summary andmention some further topic for study. 2. The model We consider an elastic 3D rectangular domain B identified with a region in R3 that is the undeformed reference configuration of the body.We denote by (x,y,z) the coordinates and we assume that B is sufficiently long in the direction Oz so that the end effects in this direction are negligible. Thus,B= (0,L)× (−h,h)× (−∞,+∞). Since B is a 3D body, which is infinite in the z direction, we refer to B as a plate, and L and 2h represent its length and its thickness, respectively. We denote in what follows by Ω = (0,L)× (−h,h) the cross section of the plate and, therefore, B = Ω × (−∞,+∞). Moreover, when h ≪ L we refer to Ω as a 2D bar, which is the topic of this paper. Let ΓD = 0× (−h,h), ΓN = (0,L)×{h}, ΓF = L × (−h,h) and ΓC = (0,L)×{−h}. The plate is clamped on ΓD × (−∞,+∞) so the displacement field vanishes there. It is free on ΓF × (−∞,+∞). On the top, ΓN × (−∞,+∞), it is being acted upon by distributed surface tractions of densityp. On the bottom, ΓC×(−∞,+∞), the platemay come in frictional contact with a foundation or obstacle describedby the functiony = Ψ(x),which for the sake of simplicity is assumedtobetime independent.Thephysical setting (thecross section of theplate) isdepicted in Fig. 1. Contact is described with the normal compliance condition (in the vertical direction) and friction with the Coulomb law of dry friction (in the horizontal direction). It is assumed that the forces and tractions vary sufficiently slowly so that the quasistatic approximation is valid. In addition, for the sake of simplicity, body forces are neglected. We denote by ν the unit normal vector on the surface of B and we use the subscripts ν and τ to represent normal and tangential components, respectively, of vectors and tensors. The time interval of interest is [0,T ] (T > 0) and a dot above a variable represents its partial time derivative. We denote by S3 the linear space of the second order symmetric tensors in R3 or, Modeling and simulations for quasistatic frictional contact... 899 Fig. 1. The setting of the problem; ΓC is the potential contact surface and Ψ describes the obstacle or foundation equivalently, the space of symmetric matrices of the order three, while “ ·” and ‖ · ‖ represent the inner products and the Euclidean norms on R3 and S3. Themathematical model for the quasistatic evolution of the state of the elastic plate, on the assumptions described above, is the following (Sofonea and Shillor, 2017). Problem P3D. Find a displacement field u : Ω × (−∞,+∞)× (0,T)→ R 3 and a stress field σ : Ω × (−∞,+∞)× (0,T)→S3 such that: σ= λ(trε(u))I3+2δε(u) in Ω × (−∞,+∞)× (0,T) (2.1) Divσ=0 in Ω × (−∞,+∞)× (0,T) (2.2) u=0 on ΓD × (−∞,+∞)× (0,T) (2.3) σν =0 on ΓF × (−∞,+∞)× (0,T) (2.4) σν =p on ΓN × (−∞,+∞)× (0,T) (2.5) −σν = λnc(uν −g)+ on ΓC × (−∞,+∞)× (0,T) (2.6) ‖στ‖¬ µ|σν| −στ = µ|σν| u̇τ ‖u̇τ‖ if u̇τ 6=0      on ΓC × (−∞,+∞)× (0,T) (2.7) u(0)=u0 in Ω × (−∞,+∞) (2.8) We turn to a short description of equations and conditions (2.1)-(2.8). First, equation (2.1) represents the linear elastic constitutive law inwhich λ and δ are thematerial Lamé coefficients, ε(u) is the linearized strain tensor associated with the displacement field u= (u,w,v), trε(u) denotes its trace and I3 represents the identity tensor ormatrix inS 3. Equation (2.2) represents the balance of the forces. We use it in this simplified version since we assume that the process is quasistatic and we neglect any body forces. Here, Divσ represents divergence of the stress fieldσ. Condition (2.3) describes the clamping on ΓD, and conditions (2.4) and (2.5) represent the traction conditions, described above. Next, we describe the contact conditions. Equality (2.6) represents the so called normal compliance condition in which g =−h−Ψ denotes the gap between the lower surface at−h and the obstacle, measured in the direction of the outward normal, λnc is the normal compliance stiffness coefficient of the foundation and r+ =max{0,r}. The normal compliance conditionwas introduced in (Martins andOden, 1987) andwas studied extensively, see, e.g., Han and Sofonea (2002), Klarbring et al. (1988, 1989), Shillor et al. (2004) and the many references therein, wheremore general normal compliance conditions can also be found. Condition (2.7) represents Coulomb’s law of dry friction in which µ is the coefficient of friction, assumed to be a positive constant. We refer to references (Eck et al., 2005; Han and Sofonea, 2002; Shillor et al., 2004; Sofonea andMatei, 2012), among themany others for explanation of this condition and related generalizations. 900 M. Barboteu et al. Finally, the problem is quasistatic since Coulomb’s condition contains the tangential spe- ed ‖u̇τ‖, therefore, we need to provide initial condition (2.8), in which u0 denotes the given initial displacement. Next, we follow (Gao, 1998; Gao and Russell, 1994) and introduce additional assumptions that allow us to derive a simplified two-dimensional model associated with Problem P3D. We assume that p= [q,f,0] with f = f(x,t) and q = q(x,t) (2.9) In other words, we assume that on the top y = h the plate is subjected to a distributed vertical load f and tangential traction q, which do not depend on z. Also, we assume that the initial displacement is such that u0 = [u0,w0,0] with u0 = u0(x,y) and w0 = w0(x) (2.10) Load (2.9) and initial data (2.10), because of the symmetry, cause deformations of the elastic plate with a displacement field u that is independent of z in the form u= [u,w,0] with u = u(x,y,t) and w = w(x,t) (2.11) Here, u is the horizontal displacement, w is the vertical one, and the displacement in the z-direction vanishes.We note that w does not depend on y, which is also an interesting part of the model. It follows that the components of the stress field donot dependon z. Therefore,we are in the situation when both the data and the unknown of this problem do not depend on z. Thus, the simplified problem with symmetry is two-dimensional and can be formulated in the domain Ω as follows Problem P2D. Under assumptions (2.9) and (2.10), find a displacement field u = (u,w) : Ω × (0,T)→ R2 and a stress field σ : Ω × (0,T)→S3 such that (2.1) and (2.2) hold in Ω×(0,T), (2.3)-(2.5) hold in ΓD×(0,T), ΓF ×(0,T) and ΓN ×(0,T), respectively, (2.6) and (2.7) hold in ΓC × (0,T) and, finally, (2.8) holds in Ω. A detailed description of Problem P2D, together with its variational analysis, including the existence and uniqueness results, can be found in our recent paper (Sofonea and Shillor, 2017). We also refer to Sofonea and Bartosz (2017) where the analysis of a dynamic contact problem for viscoelastic plates was provided. There, the model considered was two-dimensional and was derived from the full three-dimensional problem by using very similar arguments. 3. Variational form and its approximation In this Section, we present a hybrid variational formulation of Problem P2D and then its ap- proximation that leads to our algorithm for its numerical solutions. The method is based on a combination of the penalty method for normal compliance condition (2.6) and the augmented Lagrangian method for the numerical treatment of friction conditions (2.7). The Lagrangemul- tiplier is associated with the tangential frictional traction. Then, we present the approximation of the problem by using a uniform discretization of the time interval and the finite element method in the plane. We use arguments similar to those used in Barboteu et al. (2012a, 2014, 2016) and, for this reason, we skipmany of the details. Modeling and simulations for quasistatic frictional contact... 901 3.1. A hybrid variational formulation To introduce the hybrid variational problem,we need the following function spaces.We note that u = u(x,y,t) is defined on Ω × [0,T ] while w = w(x,t) is defined on [0,L]× [0,T ], which is the peculiarity of our problem, hence we seek the displacement field in the spaces V = {u ∈ H1(Ω) : u(0, ·) = 0} W = {w ∈ H1(0,L) : w(0)= 0} X = V ×W These are real Hilbert spaces with inner products defined by (u,ψ)V = ∫∫ Ω (uw+uxψx+uyψy) dxdy ψ ∈ V (w,ϕ)W = L ∫ 0 (wψ +wxψx) dx w,ϕ ∈ W (u,v)X =(u,ψ)V +(w,ϕ)W u= [u,w], v= [ψ,ϕ]∈ X We seek for the stress field in the space Q = { σ=(σij) : σij = σji ∈ L 2(Ω) } endowed with its canonical inner product. Also, we consider the function f : [0,T ]→ X given by (f(t),v)X = L ∫ 0 q(x,t)ψ(x,h) dx+ L ∫ 0 f(x,t)ϕ(x) dx for all u= [u,w], v= [ψ,ϕ]∈ X and t ∈ [0,T ]. Note that the definition of f is based on Riesz’s representation theorem, under appropriate regularity assumptions on q and f. To deal with the Lagrange multiplier, we introduce the space Xτ = {ψ(x,−h) : ψ ∈ V}, equipped with its usual norm and denote by X′τ its dual. Then, we introduce a function Iτ: Xτ → (−∞,+∞] by Iτ(ψ)= L ∫ 0 |ψ(x,−h)| dx ψ ∈ Xτ Next,wenote that for all t ∈ (0,T), friction condition (2.7) is equivalent to the subdifferential inclusion −ξτ(t)∈−µλnc(−w(t)−g)+∂Iτ(u̇(t)) in X ′ ν (3.1) where∂Iτ denotes the subdifferential of Iτ in the sense of the convex analysis. Inclusion (3.1) le- adsus to consider theLagrangemultiplier that is related to the friction traction and is considered as an additional unknown.Then, the hybrid variational formulation of the contact problemP2D, obtained by multiplying the equations with the relevant test functions and performing integra- tion by parts, is as follows. Problem PV2D. Given u0 ∈ X, find a displacement field u= [u,w] : [0,T ]→ X, a stress field σ : [0,T ]→ Q and a Lagrange multiplier ξτ : [0,T ]→ X ′ τ such that for all t ∈ [0,T ] 902 M. Barboteu et al. σ(t)= λ(trε(u(t)))I2+2δε(u(t)) in Ω (3.2) ∫∫ Ω σ(t) ·ε(v) dx− L ∫ 0 λnc(−w(t)−g)+ϕ dx (3.3) − L ∫ 0 ξτ(t)ψ(x,−h) dx =(f(t),v)X v= [ϕ,ψ]∈ X −ξτ(t)∈−µλnc(−w(t)−g)∂Iτ(u̇(t)) in X ′ ν (3.4) u(0)=u0 in Ω (3.5) 3.2. Numerical approximation We turn to the numerical approximation of ProblemPV2D.We use the finite elementmethod for spatial discretization of the domain Ω by introducing a specific finite rectangle element and a uniform discretization of the time interval [0,T ]. Since Ω is a rectangular domain, we denote by {T h} a regular family of rectangular finite element partitions of Ω that are compatible with the boundary decomposition Γ = ΓD ∪ΓN ∪ΓF ∪ΓC. Here and below, h represents the spatial discretization parameter. In the numerical examples presented in the next Section, we approximate the spaceX by thefinite dimensional space of continuous piecewise affine functions, denoted Xh. The space Q is approximated by the finite element space of piecewise constants, denoted Qh. For the discretization of the Lagrange multiplier ξτ, we consider a discrete space Y hτ ⊂ X ′ τ∩L 2(ΓC).For the timediscretization,weusea collection of discrete times{tn} N n=0which define a uniform partition of the time interval [0,T ] = ⋃N n=1[tn−1, tn] with t0 =0, tn = tn−1+k and tN = T , where N > 0 is an integer and k = T/N is the time step size.We use the notation gj = f(tj), 0¬ j ¬N, for a continuous function g(t) with values in a function space. Additional details about the discretization can be found in (Khenous, 2006a,b). Let u h 0 ∈ X h be a finite element approximation of u0. Then, using the previous notations and an implicit Euler scheme δuhkn =(u hk n −u hk n−1)/k for the approximation of the time direvative u̇(x,−h), the fully discrete approximation of ProblemPV h 2D at the time tn is the following. Problem PV h 2D. Find a discrete displacement field u hk = {uhkn } N n=0 ⊂ X h, a discrete stress field σ hk = {σhkn } N n=0 ⊂ Q h and a discrete tangential traction ξhkτ = {ξτ hk n } N n=0 ⊂ Y h τ such that, for all n =1, . . . ,N σ hk n = λ(trε(u hk n ))I2+2δε(u hk n ) (3.6) ∫∫ Ω σ hk n ·ε(v h) dx− L ∫ 0 λnc(−w hk n −g)+ϕ h dx (3.7) − L ∫ 0 ξτ hk n ψ h(x,−h) dx =(fhkn ,v h)X ∀v h =(ϕh,ψh)∈ Xh −ξτ hk n ∈−µλnc(−w hk n −g)+∂Iτ(δu hk n ) (3.8) u hk 0 =u h 0 (3.9) Concerning the numerical solution of hybrid variational ProblemPV h 2D, we have the following comments: The numerical treatment of condition (3.8) is based on a combination of the penalty method for the normal compliance condition, with an augmented Lagrangian method for the Modeling and simulations for quasistatic frictional contact... 903 friction condition.Then, the numerical approximation ofProblemPV h 2D leads at each time step to the solution of a system of nonlinear equations. Next, the unknownpair (u,ξτ) of this nonlinear system is computed by using a generalized Newton method which leads, at each iteration, to the solution of a linear non-symmetric system.Details on thismethod can be found in (Laursen, 2002). Finally, the successive linear non-symmetric systems are solved by using a Conjugate Gradient Squared Method (CGS) with Incomplete LU factorization preconditioners. We note that the contact and friction terms lead topoor conditioningof thenon-symmetric systemmatrix that is overcome by using an element-by-element factorization, seeAlart et al. (1997) for details. 4. A specific finite rectangular element In this Section, we focus on the presentation of a 2D finite rectangular element used in the discretization of the 2D linear elastic bar considered in this work. Since the domain Ω is rectan- gular, we consider finite rectangular elements Ωe for the spatial discretization of Ω. Due to the fact that the horizontal displacement u depends on x and y, and the vertical displacement w depends only on x, the use of the usual isoparametric formulation is not possible. Therefore, in each element we introduce the following displacement interpolations for u(x,y) and w(x), respectively (see Fig. 2). Fig. 2. 2D finite rectangular element We write u(ξ,η) =N(ξ,η)ue w(ξ) =M(ξ)we (4.1) where the shape functionsN andM are defined in the local coordinate system (ξ,η) by N(ξ,η) = 1 4      (1− ξ)(1−η) (1+ ξ)(1−η) (1+ ξ)(1+η) (1− ξ)(1+η)      T M(ξ)= 1 2 [ 1− ξ 1+ ξ ]T (4.2) In (4.1), we useue andwe to denote the vectors of displacements at the local nodes of the finite rectangular element, that is u e =      u1 u2 u3 u4      w e = [ w1 w2 ] = [ w3 w4 ] (4.3) We note that in contrast with the displacement interpolations used in the literature, the requiredmapping from the local coordinate system (ξ,η) to the global coordinate system (x,y) 904 M. Barboteu et al. is based on the same shape functionN, both for x and y. Then, the finite element mapping can be defined by x =N(ξ,η)xe y =N(ξ,η)ye (4.4) Here, the local vectors xe and ye contain the value of the local coordinates at the nodes of the rectangle, that is x e = [ x1 x2 x3 x4 ]T = [ x1 x1+a x1+a x1 ]T y e = [ y1 y2 y3 y4 ]T = [ y1 y1 y1+b y1+ b ]T where a and b are real positive numbers that represent the width and height of the rectangle, respectively. Based on (4.1) and (4.4), we construct inwhat follows the element stiffnessmatrixKe arising from the elementary discretization of the first term in (3.3). To this end, we consider a special finite rectangular element oriented with its sides parallel to the x,y axes. Then, the Jacobian mappingmatrix J is defined by J(ξ,η) =     ∂x ∂ξ ∂y ∂ξ ∂x ∂η ∂y ∂η     =     ∂N1 ∂ξ ∂N2 ∂ξ ∂N3 ∂ξ ∂N4 ∂ξ ∂N1 ∂η ∂N2 ∂η ∂N3 ∂η ∂N4 ∂η          x1 y1 x2 y2 x3 y3 x4 y4      In our special case, by using the shape function N given in (4.2), we obtain J(ξ,η) = 1 4 [ −1+η 1−η 1+η −1−η −1+ ξ −1− ξ 1+ ξ 1− ξ ]      x1 y1 x1+a y1 x1+a y1+ b x1 y1+ b      As is customary in the finite elementmethod,we use a vectorial notation for the components of both the strain tensor ε and the stress tensorσ, i.e. ε=    εxx εyy 2εxy    =    ux wy uy +wx    =         ∂ ∂x 0 0 ∂ ∂y ∂ ∂y ∂ ∂x         [ u(x,y) w(x) ] (4.5) and σ=    σxx σyy 2σxy    =Eε (4.6) Here, thematrixE takes into account the linear elastic constitutive behavior given in (3.2) and is defined by E=    λ+2δ λ 0 λ λ+2δ 0 0 0 δ    =    E1 E2 0 E2 E3 0 0 0 E6    Modeling and simulations for quasistatic frictional contact... 905 Using (4.1) and (4.2), equalities (4.5) and (4.6) yield ε=B(ξ,η) ·de σ=EB(ξ,η) ·de (4.7) where de = [ue,we]T is the vector of nodal variables and B(ξ,η) is the deformation strain- displacement matrix given by B(ξ,η) =        ∂N1 ∂x ∂N2 ∂x ∂N3 ∂x ∂N4 ∂x 0 0 0 0 0 0 0 0 ∂N1 ∂y ∂N2 ∂y ∂N3 ∂y ∂N4 ∂y ∂M1 ∂x ∂M2 ∂x        (4.8) It follows from forms (4.2) of the shape functionsN andM that B(ξ,η) =        − 1 2a (1−η) 1 2a (1−η) 1 2a (1+η) − 1 2a (1+η) 0 0 0 0 0 0 0 0 − 1 2b (1− ξ) − 1 2b (1+ ξ) 1 2b (1+ ξ) 1 2b (1−ξ) − 1 a 1 a        Using now (4.7), the element stiffness matrixKe of our 2D linear elastic bar can be written as K e = ∫ Ωe B T EB dΩ = 4 ∑ i=1 1 4 [B(ξi,ηi)] T E[B(ξi,ηi)]|J|ωi. Then, after some algebra, we derive the following element stiffness matrix K e = ab 16                     2 (E1 a2 + E6 b2 ) −2 E1 a2 0 −2 E6 b2 2 E6 ab −2 E6 ab −2 E1 a2 2 (E1 a2 + E6 b2 ) −2 E6 b2 0 2 E6 ab −2 E6 ab 0 −2 E6 b2 2 (E1 a2 + E6 b2 ) −2 E1 a2 −2 E6 ab 2 E6 ab −2 E6 b2 0 −2 E1 a2 2 (E1 a2 + E6 b2 ) −2 E6 ab 2 E6 ab 2 E6 ab 2 E6 ab −2 E6 ab −2 E6 ab −4 E6 a2 4 E6 a2 −2 E6 ab −2 E6 ab 2 E6 ab 2 E6 ab −4 E6 a2 4 E6 a2                     (4.9) This form, (4.9), of the element stiffness matrix is the starting point in the construction of the global stiffness matrix of the system to be solved. 5. Numerical simulations Themethod described in the previous Section has been implemented and a number of compu- ter solutions for the contact problem obtained. Here, we present numerical simulations of two settings in which the foundation is either planar or curved. In the first case, the foundation is defined by the function Ψ(x) =−0.1. In the second case, it is a circular arc lying on the circle defined by the function y = Ψ(x) that is given by (x−0.5)2+(y +1)2 =0.92. 906 M. Barboteu et al. Fig. 3. The bar in potential contact with a planar obstacle Fig. 4. The bar in potential contact with a circular obstacle The physical settings of these two configurations are depicted in Figs. 3 and 4. We pay particular attention to the mechanical aspects of the frictional contact conditions (2.6) and (2.7) and we provide a study of the dependence of the numerical solution with respect to the stiffness coefficient of the normal compliance law. In the computations, we use a rectangular mesh composed of a uniform partition of finite rectangular elements introduced in the previous Section. The spatial discretization parameters in both the x and y directions are hx =1/128 and hy =1/240, respectively. The rest of the data are the following: L =1m, h =0.05m, E =1000N/m2, G =0.3, T =1.1s, N =11, k =0.1s, f = [0,−20]N/m2 on ΓN, µ =0.4 (friction) and µ =0 (frictionless), u0 = [0,0]m in Ω. Our numerical simulations are presented in Figs. 5-10 below in which the deformed con- figuration of the bar, at the end of the time interval, is depicted. The arrows in the figures that originate on the contact surface represent the frictional contact interface tractions exer- ted by the bar on the foundation. Moreover, for the simulations presented in Figs. 5-8, we use λnc =100N/m 2. A detailed description of our numerical results is the following. We observe in Fig. 5 that in the case of a planar obstacle, the contact takes place on a large fraction of the contact nodes on the contacting surface y = −h and so the interface tractions are nonzero there. Here, we chose the friction coefficient to be µ =0.4, which is rather large, to make the effects more noticeable. Moreover, it is seen that the nodes that are on the right side of the contact region are in a stick state, a state in which u̇ = 0. The other nodes that are in active contact are in a slip state. This is a state in which u̇ 6=0 and the friction bound has been reached by the friction traction. In the frictionless case, depicted in Fig. 6, when µ =0, we note Modeling and simulations for quasistatic frictional contact... 907 that the contact forces are vertical and all the nodes that are in active contact are in the slip state (since the friction bound vanishes). Fig. 5. Deformedmesh and interface tractions in the frictional case, µ =0.4 Fig. 6. Deformedmesh and interface tractions in the frictionless case, µ =0 We turn to describe the second case of contactwith a circular obstacle, see Fig. 7. In contrast to the first case, here active contact arises in fewer nodes on the contact boundary and the rest of the boundary is in state of separationwhere the interface tractions vanish.We note thatmost of the nodes in contact are in the slip state, and only a few nodes are in the stick state. In the frictionless case, Fig. 8, all the contact nodes are in the slip state and the contact tractions are oriented in the normal direction of the foundation. Fig. 7. Deformedmesh and interface tractions in the case with friction Fig. 8. Deformedmesh and interface tractions in the frictionless case 908 M. Barboteu et al. We next describe our numerical experiments concerning the normal compliance stiffness co- efficient λnc, since we expect the solutions to the problemwith the normal compliance condition to converge, as λnc → ∞, to the solutions of the problem with the Signorini nonpenetration contact condition. The latter describes a perfectly rigid foundation. We present results obta- ined in the case of the planar obstacle in Fig. 9, for four different values of λnc. We plot the deformed configurations as well as the frictional contact interface tractions for λnc = 1N/m 2, 10N/m2, 100N/m2 and λnc =1000N/m 2.We note that for λnc =1N/m 2, a large proportion of the contact nodes are in relatively large penetration into the foundation since, in this case, the foundation is soft and so easily deformable. As the stiffness coefficient λnc becomes larger, the penetration of the bar into the foundation decreases. For λnc =1000N/m 2, we observe that the penetration is negligible. This behavior of the numerical solution shows that for a large stiffness coefficient λnc, the foundation behaves like a rigid one, and shows that we may use the normal compliance as an approximation for very stiff foundations. Fig. 9. Deformedmesh and interface tractions for various values of λnc Fig. 10. Deformedmesh and interface tractions for various values of λnc Finally, inFig. 10,wepresent similar results for the case of the circular obstacle.These results provide the same conclusion: the contact problemwith unilateral constraint may be approched by a contact problemwith normal compliance, with a sufficiently large stiffness coefficient. Our numerical results are sumarized in Table 1. Modeling and simulations for quasistatic frictional contact... 909 Table 1 λnc =1 λnc =10 λnc =100 λnc =1000 λnc =10000 max. penetration 0.242322 0.0285327 0.003274 0.000418 0.000438 (planar obstacle) max. penetration 0.188659 0.032852 0.005161 0.000984 0.000101 (circular obstacle) 6. Conclusions and comments This work presents a two-dimensional model for a long thin bar. It is simpler than a 2D model of an elastic long object in which the vertical displacement depends only on x. This makes the model easier to analyze and computationally simulate. The numerical simulations show that the model is especially useful in dealing with frictional contact. The extension of themodel to includedynamic effects is straightforward. Itmaybeof interest to extend themodel and set it as an optimal control problem by introducing the traction on the top surface as the control. Finally, more numerical simulations with different friction coefficients may be of interest (see, e.g., Barboteu et al., 2012b; Gao and Russell, 1994). References 1. Alart P., Barboteu B., Lebon F., 1997, Solutions of frictional contact problems using an EBE preconditioner,Computational Mechanics, 30, 370-379 2. Ahn J., Kuttler K.L., Shillor M., 2012, Dynamic contact of two Gao beams, Electronic Journal of Differential Equation, 194, 1-42 3. Andrews K.T., Dumont Y., M’Bengue M.F., Purcell J., Shillor M., 2012, Analysis and simulations of a nonlinear dynamic beam, Journal of Applied Mathematics and Physics, 63, 1005-1019 4. BarboteuM.,BartoszK.,KalitaP.,RamadanA., 2014,Analysis of a contact problemwith normal compliance, finite penetration and nonmonotone slip dependent friction, Communications in Contemporary Mathematics, 16, 1, DOI 10.1142/S0219199713500168 5. Barboteu M., Cheng X.-L., Sofonea M., 2016, Analysis of a contact problemwith unilateral constraint and slip-dependent friction,Mathematics and Mechanics of Solids, 21, 791-811 6. Barboteu M., Matei A., Sofonea, M., 2012a, Analysis of quasistatic viscoplastic contact problems with normal compliance,Quarterly Journal of Mechanics and Applied Mathematics, 65, 555-579 7. Barboteu M., Sofonea M., Tiba D., 2012b, The control variational method for beams in contact with deformable obstacles, Zeitschrift für Angewandte Mathematik und Mechanik, 92, 25-40 8. Eck C., Jarušek J., Krbeč M., 2005, Unilateral Contact Problems: Variational Methods and Existence Theorems, Pure and AppliedMathematics, 270, Chapman/CRCPress, NewYork 9. Gao D.Y., 1998, Bi-complementarity and duality: A framework in nonlinear equilibria with ap- plications to the contact problems of elastoplastic beam theory, Journal of Mathematical Analysis and Applications, 221, 672-697 10. Gao D.Y., Russell D.L., 1994, A finite element approach to optimal control of a ‘smart’ beam, [In:] International Conference of Computational Methods in Structural and Geotechnical Engine- ering, P.K.K. Lee, L.G. Tham andY.K. Cheung (Edit.), Hong Kong, 135-140 910 M. Barboteu et al. 11. HanW., SofoneaM., 2002,Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity, Studies in Advanced Mathematics, 30, American Mathematical Society, Providence, RI, Interna- tional Press, Somerville, MA 12. Khenous H.B., Laborde P., Renard Y., 2006a, On the discretization of contact problems in elastodynamics, Lecture Notes in Applied and Computational Mechanics, 27, 31-38 13. Khenous H.B., Pommier J., Renard Y., 2006, Hybrid discretization of the Signorini problem with Coulomb friction. Theoretical aspects and comparison of some numerical solvers, Applied Numerical Mathematics, 56, 163-192 14. Kuttler K.L., Park A., Shillor M., Zhang W., 2001, Unilateral dynamic contact of two beams,Mathematical and Computer Modelling, 34, 365-384 15. KuttlerK.L.,PurcellJ., ShillorM., 2012,Analysis and simulationsof a contactproblemfor a nonlinear dynamic beamwith a crack,Quarterly Journal ofMechanics andAppliedMathematics, 65, 1-25 16. Klarbring A., Mikelic A., Shillor M., 1988, Frictional contact problems with normal com- pliance, International Journal of Engineering Science, 26, 811-832 17. Klarbring A., Mikelic A., Shillor M., 1989, On friction problems with normal compliance, Nonlinear Analysis, 13, 935-955 18. Laursen T., 2002,Computational Contact and Impact Mechanics, Springer, Berlin 19. Martins J.A.C., Oden J.T., 1987, Existence and uniqueness results for dynamic contact pro- blems with nonlinear normal and friction interface laws,Nonlinear Analysis TMA, 11, 407-428 20. Migórski S., Ochal A., Sofonea M., 2013,Nonlinear Inclusions and Hemivariational Inequ- alities. Models and Analysis of Contact Problems, Advances in Mechanics and Mathematics, 26, Springer, NewYork 21. Panagiotopoulos P.D., 1993,Hemivariational Inequalities, Applications in Mechanics and En- gineering, Springer-Verlag, Berlin 22. Shillor M., Sofonea M., Telega J.J., 2004, Models and Analysis of Quasistatic Contact, Lecture Notes in Physics, 655, Springer, Berlin 23. ShillorM., SofoneaM., Touzani R., 2001,Quasistatic frictional contact andwear of a beam, Dynamics of Continuous, Discrete and Impulsive Systems, 8, 201-218 24. Sofonea M., Bartosz K., 2017, A Dynamic contact model for viscoelastic plates, Quarterly Journal of Mechanics and Applied Mathematics, 70, 1, 1-19, DOI:10.1093/qjmam/hbw013 25. SofoneaM.,MateiA., 2012,MathematicalModels inContactMechanics, LondonMathematical Society Lecture Note Series, 398, Cambridge University Press, Cambridge 26. SofoneaM., ShillorM., 2017,Modelling and analysis of a frictional contact problem for elastic bars, submitted to Electronic Journal of Differential Equations Manuscript received January 16, 2017; accepted for print March 31, 2017