Most numerical analysis of waveguide propagation in photonic crystal fibers are based on ideal structures with full discrete rotational symmetry and periodic boundary conditions to reduce the computational domain. However fabrication defects can yield some kind of disorder appearing as small random displacements in the air-hole distribution in the fibre cladding. The effects introduced by this disorder on the eigenmodes and propagation constants can be studied by the numerical solution of the whole cross-section of the photonic crystal fibre. Here, the finite element method is applied to the solution of the two-dimensional scalar Helmholtz equation. Nonsymmetrical meshes obtained by Delaunay triangulation are used, and a perfect matched layer is introduced outside the air-hole distribution in order to reduce the effects of spurious evanescent modes. For monomode fibres, the weak disorder only changes slightly the effective propagation constant and the field. However, for multimode fibre, the field profile of the higher-order modes deforms significantly even in the presence of weak disorder. The field profile of the fundamental mode adapts to the first row of air-holes with only small changes. But for multimode fibres the degeneracy of the high-order mode profiles, which follows from a group theory analysis of the full discrete symmetry of the fibre, is broken by the disordered air-hole distribution. Surprisingly, the effective propagation constant only suffers small changes. In summary, the results are similar to those obtained in recent experiments on multi-mode propagation in photonic crystal fibres.