This work explores charmless $B$ decays using the LHCb detector. LHCb is one of the four main experiments at the Large Hadron Collider (LHC) located at CERN, and is designed to perform CP violation measurements and to study rare decays of hadrons containing $b$ or $c$ quarks. Among charmless $B$ decays, the $B^0_s$ decay modes to final states with two light resonances ($\eta$, $\eta^{\prime}$, $\omega$, $\phi$) are particularly interesting in view of time-dependent CP violation studies. More specifically they can be used to measure the CP-violating phase difference between the $B^0_s$--$\bar{B}^0_s$ mixing amplitude and the $b\to s\bar{s}s$ decay amplitude. Among these, $B^0_s \to \phi\phi$ has been exploited by LHCb through an angular analysis of the vector-vector final state. The other modes have lower measured or expected event yields, but don't require an angular analysis. We present the results of a search for the yet unobserved $B^0_s \to \eta^{\prime} \phi$ decay using $3~\rm{fb}^{-1}$ of data collected by LHCb during the LHC Run 1 (2011--2012). The decay $B^{0}_{s} \to \eta^{\prime} \phi$ has been studied in several theoretical frameworks and the predictions for its branching fraction cover a wide range, typically from $0.05 \times 10^{-6}$ to $20 \times 10^{-6}$.\\ In the analysis presented in this thesis the $B^+ \to \eta^{\prime} K^+$ decay is used as normalisation in the computation of the branching fraction for the searched mode. The $B^{0}_{s} \to \eta^{\prime} \phi$ signal yield is obtained from a simultaneous two-dimensional fit of the reconstructed $B$ and $\eta^{\prime}$ invariant masses of the $B^{0}_{s} \to \eta^{\prime} \phi$ and $B^+ \to \eta^{\prime}K^+$ candidates in Run 1 data. No significant signal is found and, for the first time, an upper limit on the $B^{0}_{s} \to \eta^{\prime} \phi$ branching fraction is set: \begin{equation} {\cal B}(B^{0}_{s} \to \eta^{\prime} \phi)< 0.82\,(1.01)\times 10^{-6} \quad \mbox{at 90\% (95\%) CL}\,. \nonumber \end{equation} Although large theoretical uncertainties make most predictions compatible with the result of this analysis, the upper limit is significantly smaller than the central values of most of the predictions, which tends to favour the lower end of the range of predictions.\\ Furthermore, prospect studies using the Run~2 data collected in 2015 and 2016, are presented for $B^{0}_{s} \to \eta^{\prime} \phi$ and for two other decay modes already studied with Run 1 data, $B^{0}_{s} \to \eta^{\prime} \eta^{\prime}$ and $B^+ \to \phi\pi^+$. The $B^+ \to \eta^{\prime} K^+$} and $B^+ \to \phi K^+$ decay modes are used as normalisation channels. The study shows that at least the full Run 2 dataset, to be collected until the end of 2018, will be needed to aim at an observation of the $B^{0}_{s} \to \eta^{\prime} \phi$ and $B^+ \to \phi\pi^+$ decays, taking into account also the wide range of predictions for these modes, while for the already established $B^{0}_{s} \to \eta^{\prime} \eta^{\prime}$ decay the statistics collected by the end of Run 2 will allow a measurement of the $B^{0}_{s}$ lifetime.