Context. Eclipsing, spectroscopic double-lined (SB2) binaries remain the prime source of precise and accurate fundamental properties of stars. Furthermore, high-cadence spectroscopic observations of the eclipse phases allow us to resolve the Rossiter-McLaughlin (RM) effect, whose modelling offers the means to probe spin-orbit misalignment in binaries. Aims. We aim to develop a method that provides precise and accurate measurements of radial velocities (RVs) of both binary components, including the in-eclipse orbital phases where line profiles are subject to large distortions due to the RM effect. We also intend to separate spectral contributions of the primary and secondary components in the velocity space in order that a time series of the separated spectroscopic signals can be obtained throughout the binary orbit, preserving any line-profile variability (LPV) that might be present in either or both of those spectroscopic contributions. Methods. In this study, we provide a generalisation of the least-squares deconvolution (LSD) method to SB2 systems. Our LSDBinary algorithm is capable of working with both in-eclipse and out-of-eclipse spectra as input, and delivers the LSD profiles, LSD-based model spectra, and precise RVs of both binary components as output. We offer an option to account for the RM effect in the calculation of the initial guess LSD profiles and components' flux ratio, such that the effect can be modelled within the algorithm itself. In that case, the algorithm delivers both the LSD profiles and RVs, which are no longer distorted by the RM effect. Otherwise, when geometry of the RM effect is ignored in the calculation of the initial guess, the LSDBinary algorithm delivers an RV curve that contains contributions from both the orbital motion of the star and spectral line distortions due to the RM effect. Results. In this study, we provide an extensive test of the LSDBinary software package on simulated spectra of artificial binaries resembling Algol-type systems and detached binaries with similar components. We study the effects of signal-to-noise ratios (S/N) of input spectra, the resolving power of the instrument, uncertain atmospheric parameters of stars, and orbital properties of the binary system on the resulting LSD profiles and RVs measured from them. We find that atmospheric parameters have a negligible effect on the shape of the computed LSD profiles while affecting mostly their global scaling. High-resolution (R 60 000) spectroscopic observations are required in order to investigate the RM effect in detail, although a medium resolving power of R ≈25 000-30 000 might suffice when the amplitude of the effect is large. Our results are barely sensitive to the S/N of the input spectra provided they contain a sufficient number of spectral lines, such as in A-type and later stars. Finally, the orbital inclination angle and the components' radii ratio are found to have the largest effect on the shapes of the LSD profiles and RV curves extracted from them. Conclusions. The LSDBinary algorithm is ...