Purpose – The purpose of this paper is to establish a peridynamic method in predicting viscoelastic creep behaviour with recovery stage and to find the suitable numerical parameters of peridynamic method. Design/methodology/approach – A rheological viscoelastic creep constitutive equation including recovery and an elastic peridynamic equation (with integral basis) are examined and used. The elasticity equation within the peridynamic equation is replaced by the viscoelastic equation. A new peridynamic method with two time parameters, i.e. numerical time and viscoelastic real time is designed. The two parameters of peridynamic method, horizon radius and number of nodes per unit volume are studied to get their optimal values. In validating this peridynamic method, comparisons are made between numerical and analytical result and between numerical and experimental data. Findings – The new peridynamic method for viscoelastic creep behaviour is approved by the good matching in numerical-analytical data comparison with difference of < 0.1 per cent and in numerical-experimental data comparison with difference of 4-6 per cent. It can be used for further creep test which may include non-linear viscoelastic behaviour and creep rupture. From this paper, the variation of constants in Burger’s viscoelastic model is also studied and groups of constants values that can simulate solid, fluid and solid-fluid viscoelastic behaviours were obtained. In addition, the numerical peridynamic parameters were also manipulated and examined to achieve the optimal values of the parameters. Research limitations/implications – The peridynamic model of viscoelastic creep behaviour preferably should have only one time parameter. This can only be done by solving the unstable fluctuation of dynamic results, which is not discussed in this paper. Another limitation is the tertiary region and creep rupture are not included in this paper. Practical implications – The viscoelastic peridynamic model in this paper can serve as an alternative for conventional numerical simulations in viscoelastic area. This model also is the initial step of developing peridynamic model of viscoelastic creep rupture properties (crack initiation, crack propagation, crack branching, etc.), where this future model has high potential in predicting failure behaviours of any components, tools or structures, and hence increase safety and reduce loss. Originality/value – The application of viscoelastic creep constitutive model on peridynamic formulation, effect of peridynamic parameters manipulation on numerical result, and optimization of constants of viscoelastic model in simulating three types of viscoelastic creep behaviours.