Single hadron azimuthal asymmetries of positive and negative hadrons produced in muon semi-inclusive deep inelastic scattering off longitudinally polarised deuterons are determined using the 2006 COMPASS data and also combined all deuteron COMPASS data. For each hadron charge, the dependence of the azimuthal asymmetry on the hadron azimuthal angle ϕ is obtained by means of a five-parameter fitting function that besides a ϕ-independent term includes four modulations predicted by theory: sin ϕ, sin 2 ϕ, sin 3 ϕ and cos ϕ. The amplitudes of the five terms have been extracted, first, for the hadrons in the whole available kinematic region. In further fits, performed for hadrons from a restricted kinematic region, the ϕ-dependence is determined as a function of one of three variables (Bjorken-x, fractional energy of virtual photon taken by the outgoing hadron and hadron transverse momentum), while disregarding the others. Except the ϕ-independent term, all the modulation amplitudes are very small, and no clear kinematic dependence could be observed within experimental uncertainties.