## 1 Introduction  
Fatigue of metallic materials is a cumulative and irreversible process [1; 2; 3]. Under cyclic loading, the new dynamic balances are broken and the voids are generated, so that the damage will be initiated and accumulated [4; 5; 6]. The whole fatigue life is essentially a multiscale phenomenon [7; 8; 9; 10; 11], depending strongly on the material, geometry of the loaded body, the external loaded forces conditions, environmental factors, etc. The life prediction schemes of ferrite and pearlite wheel material under multiaxial loading have been through extensive research both in industry and scientific institution [12; 13; 14]. Modeling of fatigue crack growth in metallic materials encounters many problems [15; 16; 17]. Based on the framework of traditional continuum theory, the finite element method (FEM) models or various modified versions aim to find the kinetic relations among parameters near the fatigue crack tip that characterize the fatigue crack evolution [18; 19; 20]. By redefining the body, the crack is treated as a boundary so that extra criteria, such as fatigue crack growth speed and direction which guide the crack path, are necessary [21; 22; 23]. To solve these problems, Linear Elastic Fracture Mechanics (LEFM) use Cohesive Zone Elements (CZE) to deal with Mode-I crack mode and mixed-mode fracture [24; 25; 26; 27]. The number of cohesive elements increases with decreasing mesh size, yet the size of the continuum region remains the same. Frequent redefining the body is difficult and costly, especially with multiple interacting cracks. The extended finite element method (XFEM), using the local enrich functions, permits the cracks to propagate on any element surface without remeshing in every incremental crack growth [28; 29; 30; 31; 32].Both the FEM model and the XFEM method need extra criteria such as the time of forming fatigue cracks, fatigue propagation speed, direction, branch, arrest, etc. The final fatigue fracture mechanisms are associated with grain boundaries, dislocations, microcracks, anisotropy, etc [33; 34; 35]. Each of them plays an important role at a specific length scale. It is difficult to obtain the location of fatigue initiation in advance, let alone the extra criteria from the experimental data [36; 37; 38]. In contrast, a peridynamic theory (PD) uses the spatial integral equations as opposed to derivatives of the displacement field to compute internal force acting on a material particle [39; 40; 41]. Material damage is part of the constitutive model. Without the need for some special crack growth criteria, peridynamics allows the cracks to propagate at multiple sites with natural paths not only along the element boundary in a consistent framework. Furthermore, the PD theory can link the micro to macro length scales [42; 43; 44]. The fatigue crack model with peridynamics was originally proposed by [PERSON], [PERSON], and [PERSON] [45] and substantially improved by and [PERSON] [46]. Two phases of fatigue failure: crack initiation, fatigue propagation, are included in a consistent fatigue model and the parameters of the model are calibrated separately with experimental data. In ref [47]. the dynamic relaxation method is used to obtain the static solution. Nucleation and growth of a helical fatigue crack are demonstrated by using an aluminum alloy rod. [PERSON] [48] presented the conjugate gradient energy minimization method to obtain the static solution and applied this fatigue model to two-phase composite materials. The results show that the peridynamic fatigue crack model can deal with multiple crack without extra criteria to guide the crack path [49; 50; 51]. However, these peridynamic fatigue models are mainly focus on the simple linear elastic property. And the fatigue crack criteria just use the critical bond elongation technic so that it is difficulty to build effective links between the actual physical parameters and the current models. Moreover, few peridynamic fatigue models take into the cases of multiaxial fatigue failure.