In spite of the fact that the use of contact elements in modeling crack surface contact, and thus in modeling of crack closure is inherently natural, almost no effort to incorporate it in the finite element analysis of fatigue crack growth has been reported in the literature. The traditional method is based on placing truss elements on the crack surface nodes In this research contact elements are used to model crack surface contact and similar to mentioned method the determination of crack opening stress is accomplished by monitoring the state of contact elements during incrementally applied loading. Instead of determining crack opening stress of every load cycle, an algorithm that makes it possible to find crack opening stress at predetermined load cycle intervals is developed. With the developed algorithm it was possible to analyze crack closure behavior during a larger number of load cycles with less execution time. Due to simulation of relatively high number of load cycles, final stabilization of crack opening stress values after a subsequent decay in the initially "stabilized" opening stress values is observed.