This paper applies the theory of contact to masonry joints. Based on the continuity rule of strain and/or displacement field, a plane finite element is derived for the nonlinear analysis of masonry structures subjected to monotonic, cyclic, and dynamic loadings, as well as for limit-state analysis. It consists of six degrees of freedom (DOF) and considers the opening and closing mechanism of masonry joints with or without jointing material. Hereby, an incremental formulation is used. Its implantation in a finite element (FE) program is briefly shown. By investigating the interaction between axial and bending forces, it is shown, that the local loss of joint stiffness is insufficiently represented by the classical theory of plastic hinges. The fall of a temple due to seismic loading is used to illustrate the effectiveness of the applied model for structural analysis. The obtained results are compared with the response of the structure with prestressing of the columns, which increases the horizontal stability of the structure. INTRODUCTION Since the beginning of modern structural design in the 17th century, the analysis of masonry structures has been an exclusive subject traced throughout a rich literature over a period of few hundred years (Heyman 1984). Papers refer to this tradition. by interpreting the historical analytical tools in terms of upper and lower bound of yield design (Salencon 1983, 1990). In the latter, an ultimate thrust line for arches (Heyman 1988), or thrust surface for shells (Oppenheim et al. 1989) in the geometrical frontiers of the structure, establishes a static admissible stress field, which satisfies equilibrium and the material yield criterion (the lower-bound theorem). The masonry is treated as a unilateral material, capable of carrying compressive forces but not tension. With respect to the upper bound, a rigid-body collapse mechanism (Oppenheim et al. 1989) invokes a kinematical admissible velocity field, for which the power of the loads applied to the system is larger than the power that can be dissipated inside the system during its movement (the upper-bound theorem). The lower-bound or inner approach neglects the compatibility of deformation, as well as the boundary's conditions, while the upper bound's estimate does not consider the continuity of the stress field. Vilnay and Cheung (1986) investigate the effect of a distortion in the arch configuration on the balance of energy. Their variational approach leads to an energy function, whose convexity is used as the criterion for stability. However, very little is known about the dynamic behavior of masonry arches. Vilnay (1988) studies the free vibration frequencies by considering two modes of response: the elastic mode and the dynamic response of the arch after being transformed into a mechanism. Nonetheless. the aforementioned approaches can only give a first estimate of the critical or ultimate state of the structure. Continuous formulation of the local loss of stiffness, in which the theory of poroelastoplasticity is applied to masonry structures (Carrere et al. 1990), fails to consider the opening and closing mechanism of masonry joints. In view of these facts, the present paper applies, on a macroscale, the theory of contact to masonry joints. We will start with a brief review of the theory of contact and its uniaxial application to masonry structures, which will lead us to an incremental formulation of the joint mechanism. Here, the rigid-plastic behavior is expressed by the continuity rule of strain and/or displacement field. The differences to the classical theory of plastic hinges will be briefly shown. We present the fall of a temple due to seismic loading as an illustrative example. THEORY OF CONTACT The contact between two solids is locally described by the continuity of strain field and/or displacement field. i.e. |Mathematical Expression Omitted~ |Mathematical Expression Omitted~ |Mathematcal Expression Omitted~ where |Mathematical Expression Omitted~ indicates the derivation normal to the contact surface |gamma~ (Lions and Faurre 1981). Eq. (2) assures that the two solids do not penetrate into each other. By decomposing the contact surface in its adhesive and its dismantled part, i.e. |gamma~ = ||gamma~.sub.ad~ + ||gamma~.sub.dis~ (4) we obtain the general kinematical relation (Lions and Faurre 1981) on |Mathematical Expression Omitted~ on |Mathematical Expression Omitted~ Herein, ||gamma~.sub.ad~ and ||gamma~.sub.dis~, the frontiers of the active and dismantled contact surface, have to be determined. Fig. 1 demonstrates these basic relations for the uniaxial case, for which the contact surface reads as the joint's contact length being a single variable function depending only on y. This approach avoids invoking material laws. It describes the geometrical nonlinear behavior by formulating the kinematic of the displacement field. Herein, all further geometrical nonlinearities due to great deflection and rotation are neglected. APPLICATION TO MASONRY STRUCTURES Capable of carrying compressive forces but incapable of tension (Poleni 1748), the contact between the two masonry solids is assured if the surface is under compression. Therefore, we have on |Mathematical Expression Omitted~ on |Mathematical Expression Omitted~ where M and N = bending moment and axial force, respectively (positive on tension). The actual contact surface ||gamma~.sub.ad~ in the geometrical frontiers of the structure ||gamma~.sub.0~, reads then as follows: |Mathematical Expression Omitted~ |Mathematical Expression Omitted~ |Mathematical Expression Omitted~ By deriving (7) with respect to a kinematical time, we obtain the joint mechanism (for M |is not equal to~ 0 |wedge~ N |is not equal to~ 0) as follows: |Mathematical Expression Omitted~ with |Mathematical Expression Omitted~ = time dependent increments of the bending moment and the axial force. Herein, we suppose a reversible character of the presented model. The opening (closing) is coupled with a loss (gain) of stiffness in the direction of the concerned degrees of freedom (DOF), which represents the discontinuity (continuity) of strain field. For rigid-plastic behavior, this loss (gain) of stiffness is total. It is, therefore, convenient to apply the penalty method to the concerned DOF. Here, we use a 2-node element with 6 DOF, as shown in Fig. 2. The continuity of the displacement field u = u + v invokes the rigid plastic frictional joint behavior. Hence, the stiffness matrix reads as follows: |Mathematical Expression Omitted~ where ||delta~.sub.i1~ and ||delta~.sub.j1~ = Kronecker deltas and consider the continuity of displacement and/or strain field; and |omega~ = penalty factor, which is calculated in the classical manner (Kardestuncer et al. 1987), taking into account the actual global rigidity of the structure to avoid numerical instability. The following flowchart illustrates the algorithmic transformation of the presented formulation. It is part of the calculation of the residual vector in a nonlinear algorithm. Concerning iteration schemes, any explicit or implicit nonlinear algorithm capable of being implanted in a step-by-step integration procedure (Cariou 1988), can be applied. It is used in the (finite element) (FE) program CESAR-LCPC, developed by the Laboratoire des Ponts et Chaussees, Paris (Humbert 1989). FLOWCHART The calculation of residual vector is performed as follows. At time t + dt: 1. Determine M, N at the connecting beams' ends. 2. Then, for each joint element determine joint mechanism |Mathematical Expression Omitted~ using (8). 3. Actualize contact surface |Mathematical Expression Omitted~ |Mathematical Expression Omitted~ 4. Continuity control: (Determine i, j of Kronecker deltas for stiffness matrix of next iteration or time step, respectively) IF |Mathematical Expression Omitted~ ELSE i = 1 IF |Mathematical Expression Omitted~ ELSE j = 0 END 5. Repeat steps 2-5 for each joint element. Assemble vector of unbalanced forces. M/N-INTERACTION As just mentioned, each joint's opening presents a local discontinuity of the displacement field (5). The structure converts into a mechanism such that motion of the system takes place through rotations at the joints. This imposes the question, to which extend the local loss of joints stiffness can be idealized by using the classical theory of plastic hinges (Chen and Han 1988). Consider the velocity field at the edge of opening, as shown in Fig. 3. The corresponding power of applied bending moment M and axial force N is given by |Mathematical Expression Omitted~ with |Mathematical Expression Omitted~ as compatibility condition. The corresponding power that can be dissipated inside the system is |Mathematical Expression Omitted~ |alpha~ depends on the applied constitutive law for |sigma~. It should be carefully, chosen between |Mathematical Expression Omitted~ where the lower bound corresponds to an elastic material behavior and the upper bound to rigid-plastic behavior. By equating (12) to (14) and applying the compatibility condition (13), we obtain the fundamental equation of eccentricity per unit length dx |Mathematical Expression Omitted~ with |Mathematical Expression Omitted~ where |M.sub.0~, and |N.sub.0~ = bending moment and axial force, at which the opening of the joint occurs. With respect to the applied elasto-plastic approach, the corresponding loading surface reads then in terms of the nondimensionalized stress resultants m = M/|M.sub.0~, and n = N/|N.sub.0~ as follows F(m, n) = m/n - (|alpha~ - 1) = 0 (18) The interaction curve of (18) is drawn in Fig. 4 and compared with typical axial and bending interaction curves of different sections (Chen and Han 1988). The differences in shape, curvature, and orientation illustrate clearly that the local loss of joints stiffness is insufficiently represented by the classical theory of plastic hinges. By shape, a Mohr-Coulomb criterion with zero cohesion is necessary to idealize the local joint's behavior. However, it could not take into account the reversible character of the joint mechanism and would, therefore, overestimate the energy dissipation due to joint plastification. This underlines the chosen approach by the theory of contact. EXAMPLE To illustrate the effectiveness of the proposed model, we treat the vertical cross section of a temple subjected to seismic loading. The geometry is taken from the Temple of Zeus in Olympia (Andronicos 1987). Fig. 5 shows the geometrical details (Mueller and Vogel 1990), the schematic element mesh, and the material properties. The core of the temple, which plays a major role concerning horizontal stability, is neglected, as well as all other nonlinear geometrical effects due to vertical displacement. Hence, the treated structure consists of a plane massive panel on six marble columns. The structure is idealized using triangle elements for the frieze and plane beam elements for the columns connected by the proposed joint elements. Herein, the angular inertia of the columns is neglected. The axial loading in the columns due to the own weight of the frieze and the columns is idealized by imposing an axial compressive deflection on the columns. These prestressing displacements are only considered for the calculation of the generalized stresses in the joints. The structure is subjected to a lateral ground motion. The accelerogram of the earthquake, taken from ISPRA (Cooperative 1990), is drawn in Fig. 6. The response of the structure is shown in Fig. 7. The joints open and the structure converts into a mechanism leading to failure, i.e. to the singularity of the system's global stiffness matrix. We compare it with the response of a structure with increased prestressing in which the joints remain closed. Herein, the augmentation of axial compressive deflection does not correspond to an increase of mass. This could be of great importance for the engineering earthquake design of existing masonry structures, with prestressing increasing the horizontal stability. Furthermore, by binding the joints together, it increases the axial stiffness of the system. The axial frequencies augment and will therefore be less activated by vertical accelerations, which were not considered in the present study, but which have a significant influence on the stability of the structure. CONCLUSION The present paper sought to approach the nonlinear behavior of masonry joints by the theory of contact. Hereby, an incremental formulation for the opening and closing mechanism of a plane finite joint element from the continuity of strain and displacement field was derived. This was numerically treated by using the penalty method. By investigating the axial and bending interaction curve of the proposed element, it was shown that the local loss of stiffness due to joint opening is insufficiently represented by the classical theory of plastic hinges. With the example of a collapsing temple subjected to seismic loading, the effectiveness of the applied theory of contact for the earthquake engineering design of masonry structures was illustrated. The effectiveness of the proposed model corresponds, however, to the scale of the engineering problem to which it is applied, a question which remains to be answered by the consulting engineer. APPENDIX I. REFERENCES Andronicos, M. (1987). Olympia: Ausgrabungen und Museum. Ekdotike Athenon S.A., Athens, Greece (in German). Cariou, D. (1988). 'Calculs par elements finis de structures poutres bidimensionnelles et tridimensionnelles en grands deplacements elastoplastiques sous sollicitations dynamiques,' PhD thesis, Ecole Nationale des Ponts et Chaussees, Paris, France (in French). Carrere, A., Coussy, O., and Fauchet B. (1990). 'Stabilite des barrages-poids: Apports de la mecanique des milieux poreux.' Ann. Ponts Chaussees, 3, 10-30. Chen, W. F., and Han, D. J. (1988). Plasticity for structural engineers. Springer Verlag, New York, N.Y. 'Cooperative research on the seismic response of reinforced concrete structure.' (1990). Final rep., ISPRA. Ispra, Italy. Heyman, J. (1984). The masonry arch. Ellis Horwood, Chichester, England. Heyman, J. (1988). 'Poleni's problem.' Proc., Inst. Civ. Engrs., Part 1, 84, 737-759. Humbert, P. (1989). 'CESAR-LCPC: Un code general de calcul par elements finis.' Bull. Liaison. Lab. Ponts Chaussees, 160, 112-116. Kardestuncer, H., Norrie, D., and Brezzi, F. (1987). Finite element handbook. McGraw-Hill, New York, N.Y. Lions, J. L., and Faurre, P. (1981). Cours d'analyse numerique. Deuxieme partie: Notes d'optimisation. Ecole Polytechnique, Palaiseau, France. Mueller. W., and Vogel, G. (1990). dtv-Atlas zur Baukunst. Tafeln und Texte. Band 1, 8. Auflage, Deutscher Taschenbuchverlag, Muenchen, Germany. Oppenheim, I. J., Gunaratnam, D. J., and Allen, R. H. (1989). 'Limit state analysis of masonry domes.' J. Struct. Engrg., 115(4), 868-882. Poleni, G. (1748). Memorie istoriche della gran cupola del Tempio Vaticano. 1748, Padua, Italy (in Italian). Salencon, J. (1983). Calcul a la rupture et analyse limite. Presses de l'Ecole Nationale des Ponts et Chaussees, Paris, France. Salencon, J. (1990). 'An introduction to the yield design theory and its application to soil mechanics.' Eur. J. Mech. A/solids, 9(5). 477-500. Vilnay, O. (1988). 'Dynamical behaviour of three-voussoir arch.' J. Struct. Engrg., 114(5), 1173-1186. Vilnay, O., and Cheung, S.S.(1986). 'Stability of masonry arches.' J. Struct. Engrg., 112(10), 2185-2199. APPENDIX II. NOTATION The following symbols are used in this paper: b = joint width; F(m, n) = loading surface depending on m, n; K = element stiffness matrix; M = bending moment; |M.sub.0~, |N.sub.0~ = bending moment and axial force, at which joint opening occurs; m = M/|M.sub.0~ = nondimensionalized stress resultant of bending moment; N = axial force; n = N/|N.sub.0~ = nondimensionalized stress resultant of axial force; |P.sup.ex~ = power of external applied forces; |P.sup.in~ = power dissipated inside system; u = displacement field; u, v = axial and bending components of displacement field u; |alpha~ = material constant; |gamma~ = contact surface, contact length respectively for uniaxial case; ||gamma~.sub.ad~ = adhesive part of contact surface; ||gamma~.sub.dis~ = dismantled part of contact surface; ||gamma~.sub.0~ = geometrical frontiers of structure; ||delta~.sub.ij~ = Kronecker delta; |epsilon~ = strain field; ||epsilon~.sub.0~ = strain in reference axis; |sigma~ = stress field; ||phi~.sub.1~, ||phi~.sub.2~ = rotation DOF at joint's edges; |chi~ = curvature; |omega~ = penalty factor; |Mathematical Expression Omitted~ = joint mechanism factor; (*) = d/dt = time derivative; and |Mathematical Expression Omitted~ = normal derivative.