This paper presents a novel Matlab computational program designed for simulating the time-dependent behavior of compressed fracture surfaces. Based on the minimum complementary energy principle and the dependence of relaxation modulus on time, the program analyzes the variations over time in normal creep deformation, creep rate, and damage area of fracture surfaces under different normal stress conditions during compression. Comparative analysis is conducted on the changes in the damaged area of fracture surfaces before and after 100 h of creep. It’s shown that the normal closure velocity gradually slows down with the increase of time, which is caused by the rapid decay of the relaxation modulus in the initial stage of creep and the slow decline in the later stage. The damage area accumulates continuously with the increase of time, but the rate is decreasing. The novel approach considers the interactions between multiple asperities, ensuring the uniqueness of the solution to the contact problem; it improves the calculation efficiency and reduces the storage space. It provides a new method for studying the time-dependent law of creep and damage of rock fractures under normal stress.