A numerical method to simulate the periodic fibre suspension in viscoelastic fluid is developed with the completed double layer boundary element method (CDLBEM). The periodic summations that arise in the formulation were well handled by Ewald summation technique to speed up the convergence rate in the computation. The formulation for velocity field in periodic fibre suspension in viscoelastic fluid is derived and is used to simulate multiple fibres suspended in a viscoelastic shear flow. Simulations are carried out for various fibre aspect ratios and volume fractions ranging from dilute to concentrated regimes. Numerical results of macroscopic theological properties of the system are compared to available experiments on viscoelastic fibre suspensions, and are found to agree reasonably well with the experimental data. (c) 2008 Elsevier Ltd. All rights reserved.