Fault systems are restricted by the elasticity equation under the mechanical constraints of both the weakening stiffness under cyclic rupturing and constant fault rupture strength. In addition, according to existing experimental data, the ratio θ of stiffness after rupture to stiffness before rupture under the cyclic rupturing of rocks is constant and is not always one. So, based on these constraints on rock fracture mechanics, we propose a new spring network model which has various θ values to simulate the evolution of fault systems. Our model is supported not only by fracture mechanics but also by natural observations of fault systems (e.g., "time-predictable recurrence" for large earthquakes). The result of our simulation shows clearly that varying θ affects patterns of spatial and temporal evolution of fault systems. An original ANSI-C program is presented for the simulation algorithm of our model. This simulation is a valuable tool for investigating fault system evolution in tectonically active regions of the upper lithosphere.