Transcranial focused ultrasound is a novel technique for the noninvasive treatment of brain diseases. The success of the treatment greatly depends on achieving precise and efficient intraoperative focus. However, compensating for aberrated ultrasound waves caused by the skull through numerical simulation-based phase corrections is a challenging task due to the significant computational burden involved in solving the acoustic wave equation. In this article, we propose a promising strategy using the coupling of the boundary integral equation method (BIEM) and the finite element method (FEM) to overcome the above limitation. Specifically, we adopt the BIEM to obtain the Robin-to-Dirichlet maps on the boundaries of the skull and then couple the maps to the FEM matrices via a dual interpolation technique, resulting in a computational domain including only the skull. Three simulation experiments were conducted to evaluate the effectiveness of the proposed method, including a convergence test and two skull-induced aberration corrections in 2D and 3D ultrasound. The results show that the method's convergence is guaranteed as the element size decreases, leading to a decrease in pressure error. The computation times for simulating a 500 kHz ultrasound field on a regular desktop computer were found to be 0.47 ± 0.01 s in the 2D case and 43.72 ± 1.49 s in the 3D case, provided that lower-upper decomposition (approximately 13 s in 2D and 2.5 h in 3D) was implemented in advance. We also demonstrated that more accurate transcranial focusing can be achieved by phase correction compared to the noncorrected results (with errors of 1.02 mm vs. 6.45 mm in 2D and 0.28 mm vs. 3.07 mm in 3D). The proposed strategy is valuable for enabling online ultrasound simulations during treatment, facilitating real-time adjustments and interventions.Copyright © 2023 Elsevier B.V. All rights reserved.