High-Intensity Focused Ultrasounds (HIFU) technology has several applications in the medical field, such as tumour ablation and drug delivery. HIFU technology uses ultrasound to induce temperature increase and, eventually, cavitation in each location to destroy a particular tissue or activate a given drug. One of the limitations of HIFU is the side damage caused because of the inaccuracy given by the very complex composition of the human body. Computational fluid dynamics analysis can help evaluate the size of the cavitation zone that forms at the focal point and estimate eventual deviation from the expected location. The numerical simulation can also help estimate the temperatures reached at the focal point or design new devices. The numerical model proposed to simulate HIFU is based on the standard compressible Volume of Fluid (VoF) method. The transport equation of the liquid volume fraction is coupled with the continuity and momentum conservation equations to solve for the flow field. Since the size of incipient bubbles formed by ultrasonically induced cavitation is generally smaller than the single cell, sub-grid models are used to describe nucleation, bubble dynamics and collapse unless the critical size of the bubbles or bubble cluster is enough to guarantee the creation of an interface. This last condition coincides with several adjacent cells having a volume fraction 0 (<1e-6). When a multi-cells bubble is created, the sub-grid models are de-activated, and the newly formed bubble is simulated with the compressible VoF method. The nucleation was simulated by adopting an empirical correlation based on experimental evidence reported in several works. However, it is important to point out that developing a nucleation model to build a fully predictive numerical model to simulate HIFU is essential. The Keller--Miksis (KM) equation described the bubbles' dynamics. This equation is derived by assuming spherical symmetry and applying the Navier--Stokes equations projected onto the radial coordinate. Finally, the collapse was determined by a cavitation threshold that was, once again, derived from an empirical correlation based on experimental evidence as the available thresholds, such as the so-called "Blake threshold" are not suitable for the problem under study. The vibration frequency and forcing pressure were used to determine the fluctuating pressure boundary condition for each transducer composing the transducer array used in the experiments. The numerical model was used to simulate several cases present in the literature. This work presents the first step towards realising a predictive model for simulating ultrasonically induced cavitation in biological fluids. The framework was realised with the scope of evaluating bubbles' cloud location, size and cavitation intensity, thus allowing better treatment control.