【正文】 f solutions ??u as ??? can be shown to converge to the solution of the constrained problem (5). 12 In practice, penalty methods are used successfully even when the equilibrium state is not associated with minimization of the total potential energy. Due to the occurrence of peration for finite values of the penalty parameter ? ,a unique definition of the contact surface C and gap function g is not readily available for use in (11).A crucial step in the proposed formulation is the introduction of a welldefined penalty regularization for (10a) and (10b),such that ???? )()( ?? ? gp (12) With the aid of (12),equation (10a) bees }:{ )(2 1 ???? ???? ? dtwdvbwpdvtdi v w q xtt nxxx ??? ? ? ??? ? ???? 0)2()1( 22)2(11)1( ??????????? ?? CC dnwgdnwg ???? (13) Each of the last two integrals of (13) identical in form to the boundary integral emanating from penalty regularization of a Signorini problem. Although in the present context these two integrals are clearly coupled by the definitions (2a) and (2b) of gap functions )(?g , this deposition has important putational implications, as will be discussed in the next section. The developments presented here can be easily extended to enpass the mulibody contact problem by reducing it to a series of coupled twobody problems. 3. APPLICATION TO TWODIMENSIONAL CONTACT The remainder of this article is devoted to the discretization of twodimensional contact problems based upon the penalty formulation and the identification of gap functions on both surfaces, as suggested in Section two key issues to be addressed here are the geometric construction of contact elements and the choice of admissible fields for the finite element approximation. 14 . Discretization of the contact surfaces A continuous discretization of the contact surfaces is advocated. Surfaces )1(C and )2(C are uniquely determined as those on which 0)1( ?g and 0)2( ?g , respectively. Setting aside implementational details, contact surface )1(C is discretized by a series of normal projections (not necessarily closestpoint) from 1t?? to 2t?? ,as shown in Figure 2. An analogous procedure is followed for the discretization of )2(C .Consequently, each contact element )(?eC relates a single spatial element edge on surface )(?C to the opposite surface. The apparent nonsmoothness of the discrete boundaries results in discontinuity of unit normals and gap functions. No attempt is made here towards circumventing this problem, although its effect might not be negligible, especially in problems of rolling contact. For a special discretization that results in smooth boundary representation in twodimensions, see Reference 13. . Finite element fields A specific contact finite element is suggested in this section, based upon equeations (10a) and (12). The choice of finite dimensional fields is guided by previous works especially on the Signorini admissible displacement and pressure fields lead to unilateral contact elements that are able to stably replicate the imperability condition (. they satisfy the underlying LBB condition) and are accurate. Such convergence analysis, although not currently available for the kinematically nonlinear twobody contact problem addressed here, provides a guideline for the selection of admissible fields. Numerical integration is typically employed for all boundary terms on the contact surface. Discounting discretizations that use straightedge spatial elements(. threenode triangles and fournode quadrilaterals),the gap functions )(?g within a single contact element are nonpolynomial with respect to local boundary coordinate systems. Thus, 15 numerical integration on )(?C is generally inexact and all integration rules introduce errors that directly influence the formulation. The proposed contact element is based on displacement fields emanating from standard Q9 (ninenode isoparametric) elements and employs Simpson’s integration rule.