Mucociliary clearance is the first defense mechanism of the respiratory tract against inhaled particles. This mechanism is based on the collective beating motion of cilia at the surface of epithelial cells. Impaired clearance, either caused by malfunctioning or absent cilia, or mucus defects is a symptom of many respiratory diseases. Here, by exploiting the Lattice Boltzmann particle dynamics technique, we develop a model to simulate the dynamics of multiciliated cells in a two-layer fluid. Firstly, we tuned our model to reproduce the characteristic length- and time-scales of the cilia beating. We then check for the emergence of the metachronal wave as a consequence of hydrodynamic mediated correlations between beating cilia. Finally we tune the viscosity of the top fluid layer to simulate the mucus flow upon cilia beating, and evaluate the pushing efficiency of a carpet of cilia. With this work, we build a realistic framework that can be used to go explore several important physiological aspects of mucoliary clearance.