Let $C$ be the integration path of a semicircle of radius $R$ denoted by $S(R)$ with a straight line in the real axis from $-R$ to $R$. Then we have that:
$$\oint_{C} \frac{1}{z^4+z^2+1}dz = \left(\int_{0}^{R} +\int_{S(R)} + \int_{-R}^{0}\right) \frac{1}{z^4+z^2+1}dz.$$
Note that:
$$\int_{0}^{R} \frac{1}{z^4+z^2+1}dz = \int_{-R}^{0}\frac{1}{z^4+z^2+1}dz$$
with the substitution $u = -z$. Now,by residue theorem we have that:
$$\oint_{C}\frac{1}{z^4+z^2+1}dz = 2i\pi \left(\frac{1}{2(1+i\sqrt{3})}+\frac{1}{2(-1+i\sqrt{3})}\right) = \frac{\sqrt{3}\pi}{2}$$
Where the terms in the LHS are $2i\pi \sum_{j=1}^2 Res(f(z_j))$ where $z_1,z_2$ are the two roots of the polynomial that are enclosed by $C$. I allready prooved that:
$$\lim_{R\to\infty} \int_{S(R)} \frac{1}{z^4+z^2+1}dz = 0$$
With all these results I could conclude that:
$$\int_{0}^{\infty} \frac{1}{x^4+x^2+1} = \sqrt{3}\frac{\pi}{4}$$
However, by a numerical simulation, the value seems to be $\pi/2\sqrt{3}$, but I have no idea where I do a wrong move.
It looks like you've made a mistake in the process of computing the residues.
We know $f$ has poles at $\pm \frac{1}{2} \pm \frac{\sqrt{3}}{2} i$, and assuming we're taking the semicircle in the upper half plane this means the relevant poles are $\pm \frac{1}{2} + \frac{\sqrt{3}}{2} i$.
Then we compute the residues to be $\frac{1}{\pm 3 + \sqrt{3}i}$ (which seems to disagree with what you computed in your question), so that the residue theorem gives us
$$ \oint f\ dz = 2 \pi i \left ( \frac{1}{3 + \sqrt{3}i} + \frac{1}{-3 + \sqrt{3}i} \right ) = \frac{\sqrt{3}\pi}{3} $$
Since, as you argue, your integral should be half of this, we get an answer that agrees with your numerical simulation.
I hope this helps ^_^