A detailed study of the low frequency behaviour of the phonon spectrum for 3-dimensional tiling models of icosahedral quasicrystals is presented, in commensurate approximations with up to 10 336 atoms per unit cell. The scaling behaviour of the lowest phonon branches shows that the widths of the gaps relative to the bandwidths vanish in the low frequency limit. The density of states at low frequencies is calculated by Brillouin zone integration, using either local linear or local quadratic interpolation of the branch surface. For perfect approximants it appears that there is a deviation from the normal ω²-behaviour already at relatively low frequencies, in the form of pseudogaps. Also randomized approximants are considered, and it turns out that the pseudogaps in the density of states are flattened by randomization. When approaching the quasiperiodic limit, the dispersion of the acoustic branches becomes more and more isotropic, and the two transversal sound velocities tend to the same value. The dynamical structure factor is determined for several approximants, and it is shown that the linearity and the isotropy of the dispersion are extended far beyond the range of the acoustic branches inside the Brillouin zone. A sharply peaked response is observed at low frequencies, and broadening at higher frequencies. To obtain these results, an efficient algorithm based on Lanczos tridiagonalisation is used.