A nonlinear mixed finite element formulation based on the Hellinger-Reissner variational principle is developed for planar contact stress analysis. The formulation is based on the updated Lagrangian approach and accounts for geometric nonlinearity. In the mixed model, both displacements and stresses are approximated independently and this approach has in general been found to be more accurate than the displacement finite element model, especially for contact problems since it avoids the extrapolation of stresses computed at the Gauss points to the boundary nodes. An algorithm based on the penalty technique for equality constraints has been developed to handle the interface boundary conditions arising in a contact problem. The algorithm automatically tracks potential contact nodes, detects overlap during any load step and iteratively restores geometric compatibility at the contact surface. The classical Hertz contact problem is solved to validate the algorithm. The mixed formulation algorithm in cylindrical coordinates is applied in conjunction with the penalty based algorithm to solve the contact problem in layered cylindrical bodies. Static condensation techniques are used to condense out the discontinuous components of stresses at the element level. The contact stress distribution and variation of contact area with load is computed for different loading situations. Furthermore, the effect of the difference in the relative magnitudes of the moduli of the layers on the stability of the contact algorithm is investigated. / Master of Science