In heterogeneous networks of coupled oscillators, phase frustration typically prevents the emergence of complete synchronization in the Sakaguchi-Kuramoto (SK) model. In this study, we propose an analytical framework to overcome this barrier and induce complete synchronization in oscillators governed by phase-frustrated bi-harmonic coupling. We derive a general set of natural frequencies correlated with the network's degree heterogeneity, along with the parameters involved in the bi-harmonic coupling function that lead to complete synchronization (r=1) in the presence of the harmonic coupling terms (K₁, K₂ ≠ 0). On top of that, we found hysteresis in the synchronization transition in the case of scale-free networks, indicating a first-order (discontinuous) phase transition, whereas Erdos–Renyi networks exhibit a second-order (continuous) synchronization transition. Furthermore, we use mean-field approximation to determine the critical coupling strength for the synchronization transition in the absence of first-harmonic coupling (K₁=0). Here, the obtained optimal natural frequencies scale linearly with the node degree, and the critical coupling strength for the onset of synchronization is derived analytically from the self-consistent equations. In this specific regime, we observe distinct dynamical disparities: the second harmonic drives an explosive first order (or second order) transition for the second order parameter (r₂), while the first order parameter (r₁) remains suppressed in the forward direction but emerges during the backward transition. These findings remain robust with higher-order harmonic coupling schemes, as well as across a diverse range of synthetic and empirical networks, including scale–free, Erdos–Renyi, Zachary Karate Club and C.elegans neural network, demonstrating their general applicability.