Three-dimensional numerical simulations of laminar flow around a circular cylinder with sinusoidal variation of cross-section along the spanwise direction, named ‘wavy cylinder’, are performed. A series of wavy cylinders with different combinations of dimensionless wavelength (λ/Dm) and wave amplitude (a/Dm) are studied in detail at a Reynolds number of Re = U∞Dm/ν = 100, where U∞ is the free-stream velocity and Dm is the mean diameter of a wavy cylinder. The results of variation of mean drag coefficient and root mean square (r.m.s.) lift coefficient with dimensionless wavelength show that significant reduction of mean and fluctuating force coefficients occurs at optimal dimensionless wavelengths λ/Dm of around 2.5 and 6 respectively for the different amplitudes studied. Based on the variation of flow structures and force characteristics, the dimensionless wavelength from λ/Dm = 1 to λ/Dm = 10 is classified into three wavelength regimes corresponding to three types of wake structures. The wake structures at the near wake of different wavy cylinders are captured. For all wavy cylinders, the flow separation line varies along the spanwise direction. This leads to the development of a three-dimensional free shear layer with periodic repetition along the spanwise direction. The three-dimensional free shear layer of the wavy cylinder is larger and more stable than that of the circular cylinder, and in some cases the free shear layer even does not roll up into a mature vortex street behind the cylinder. As a result, the mean drag coefficients of some of the typical wavy cylinders are less than that of a corresponding circular cylinder with a maximum drag coefficient reduction up to 18%. The r.m.s. lift coefficients are greatly reduced to practically zero at optimal wavelengths. In the laminar flow regime (60 ≤ Re ≤ 150), the values of optimal wavelength are Reynolds number dependent.