Fault damage deteriorates the rock mass quality which is a significant factor influencing rock mass stability. This study uses the discrete element method to simulate the formation of fault core and damage to host rocks, and to reveal the structure of the damage zone.The numerical simulation was performed using the PFC 2D 5.0 software, employing commonly used parameters for stratigraphic scale simulations by previous studies.By analyzing the fracture count-displacement, strain energy-displacement,and principal stress-displacement curves, fault damage can be classified into four distinct stages: the diffusing damage stage, the clustering damage stage, the nucleation damage stage, and the sliding damage stage. During the diffusing damage stage, prior to the arrival of the fault core, damage primarily yielded under the far-field stress field. Subsequently, fault gradually formed in a clustered manner along the Riedel shear which is the orientation of the potential fault. The nucleation stage represents the primary phase of damage development during which Mode I fractures within the damaged rock were connected by Mode II fractures, while stress and strain energy underwent dramatic fluctuations during accumulation and fracture initiation, continuously linking the remaining intact rock bridges and facilitating the formation of fault core through yield shear deformation. As a result, the fracture count at the fault tip increased rapidly. When the tip damage zone was “abandoned” upon the passage of the fault core, it transformed into wall damage zone on both sides of the fault. Following this, active side experienced abrupt slide along the newly formed fault core, leading to elastic rebound and a rapid increase in fault throw. After the fault core had passed, the irregular boundary of the fault core induced erosion of the host rock and created new damage. This process continuously gave rise to small-scale stress concentration, fault core expansion and stress release. On the whole, the fracture frequency decreases with the spacing to fault core. The damage zone can also be divided by the varying frequencies of fracture in damage zone. Our field survey of the Baini-Shawan fault in the Pearl River Delta reveals that its fault damage feature, such as fault core width, is consistent with the simulation results.The field work from previous studies also align well with the simulation presented in this paper. This simulation study elucidates the fault damage process and reconstructs the structural features of the fault damage zone, which holds significant implications for understanding the causes of hazards associated with rock mass quality deterioration in fault damage zone as well as for implementing effective risk prevention measures.