In this study, a new traction‐separation based constitutive model for use in finite element simulation of masonry joints under complex loading conditions is developed for cohesive elements. The proposed model is formulated using damage parameters and plastic deformation with mutual couplings, and can accurately simulate the complex nonlinear behaviors of masonry joints considering hardening or softening of strength and stiffness degradation. To enhance the numerical stability of the model, plasticity and damage are separated algorithmically and implemented in two phases. In the first phase, the plastic deformations are treated using a multi‐surface plasticity model composed of a smooth hyperbolic yield surface for tension‐shear mixed‐mode failure and an elliptical cap primarily for the compressive failure. This is implemented in effective stress space and helps restrict the evolution of yield surfaces with no softening, significantly enhancing the efficiency of stress return mapping by the closed point projection method. In addition, an adaptive sub‐stepping scheme is adopted to further improve the robustness of the numerical implementation. In the second phase, nominal stresses are computed from the effective stresses using damage parameters. The evolution of these damage parameters is defined in terms of plastic work which is defined by a polynomial form, and is recommended in this study for a better calibration capability. Improvements are made in the formulation of compressive cap including incorporation of hardening of strength and stiffness degradations as these are ignored in existing interface models. This approach helped improve simulation of masonry under cyclic loads with tension‐compression transitions. For the structural level applications, the interface model is implemented within a finite element program, which is utilized to simulate failure of a number of masonry specimens under in‐plane/out‐of‐plane monotonic/cyclic loading. The simulated results are rigorously validated with existing experimental data that shows a good potential in modeling masonry structures.