Start with the time-independent Schrödinger equation for the harmonic oscillator:
$$ -\frac{\hbar^2}{2m} \frac{d^2\psi(x)}{dx^2} + \frac{1}{2}m\omega^2x^2\psi(x) = E\psi(x) $$The classical momentum $p(x)$ is defined as:
$$ p(x) = \sqrt{2m(E - V(x))} = \sqrt{2m\left(E - \frac{1}{2}m\omega^2x^2\right)} $$Assume the WKB form of the wavefunction:
$$ \psi(x) = A(x)\exp\left(\pm \frac{i}{\hbar} \int p(x') dx'\right) $$where $A(x)$ is a slowly varying amplitude.
Identify the classical turning points $x_1$ and $x_2$, where $E = V(x)$. For the harmonic oscillator, these are:
$$ x_1 = -\sqrt{\frac{2E}{m\omega^2}}, \quad x_2 = \sqrt{\frac{2E}{m\omega^2}} $$Apply the quantization condition:
$$ \oint p(x) dx = 2\pi n \hbar, \quad n = 0, 1, 2, \dots $$Here, $\oint p(x) dx$ represents the integral over one complete cycle between the turning points.
For the harmonic oscillator:
$$ p(x) = m\omega \sqrt{\frac{2E}{m\omega^2} - x^2} $$Evaluate the integral:
$$ I = \oint p(x) dx = 4 \int_0^{x_2} m\omega \sqrt{\frac{2E}{m\omega^2} - x^2} dx $$Substitute $x = x_2 \sin(\theta)$, and solve:
$$ I = 4 m\omega \int_0^{\pi/2} x_2^2 \cos^2(\theta)d\theta = 4 m\omega x_2^2 \cdot (\pi/4) = \pi m\omega x_2^2 = \pi m\omega \left(\frac{2E}{m\omega^2}\right) = \frac{2\pi E}{\omega} $$Equating $I = \frac{2\pi E}{\omega} = 2\pi n \hbar$, solve for $E_n$:
$$ \frac{2\pi E}{\omega} = 2\pi(n+1/2) \hbar $$ $$ E_n = (n + \frac{1}{2})\hbar \omega $$The WKB method yields:
$$ E_n = \left(n + \frac{1}{2}\right)\hbar \omega $$WKB method gives approximate result, but in this case it matches the exact quantum mechanical result for the harmonic oscillator. You might want to think about why it does so.