$\newcommand\tH{\tilde H}$This problem is one of real algebraic geometry, which can be solved purely algorithmically. In Mathematica, such algorithms are implemented by `Reduce` and similar commands. 

Here is a solution with Mathematica: 

[![enter image description here][1]][1]

---

Here is a "more human" proof: Let $f(x,y,z)$ stand for the left-hand side of your desired inequality. We want to show that $f(x,y,z)\ge3/2$ if $x,y,z>0$ and $xyz=1$. Equivalently, we want to show that 
$$h(x,y):=\big(f(x,y,\tfrac1{xy})-\tfrac32\big)
2 x^3 y^3 (x + y) (1 + x^2 y) (1 + x y^2)\ge3/2 \\ 
=2 x^{10} y^{10}+2 x^9 y^8+2 x^8 y^9+2 x^7 y^7-3 x^7 y^6-3 x^6 y^7-3 x^6
   y^4+2 x^6 y^2-6 x^5 y^5+2 x^5 y^3+2 x^5-3 x^4 y^6-3 x^4 y^3+2 x^4
   y+2 x^3 y^5-3 x^3 y^4+2 x^2 y^6+2 x y^4+2 y^5\overset{\text{(?)}}\ge0$$
for $x,y>0$. 

Consider the partial derivatives $p(x,y):=h_x(x,y)$ and $q(x,y):=h_y(x,y)$. Let $r_1(y)$ and $r_2(x)$ denote the resultants of the polynomials $p(x,y)$ and $q(x,y)$ w.r. to $x$ and $y$, respectively. By symmetry, the real roots of $r_1(y)$ and $r_2(x)$ are the same: $z_1:=0$, $z_2:=1$, and a certain algebraic number $z_3\approx0.8180077783$. Therefore the critical points of $h$ in $(0,\infty)^2$ (if any) are of the form $(z_j,z_k)$ for $j,k\in\{2,3\}$. It is straightforward to check that $h(z_j,z_k)\ge0$ for $j,k\in\{2,3\}$. 

So, it remains to check that the boundary values of $h$, near the boundary of the set $[0,\infty]^2$, are $\ge0$. 

We have $h(0,y)=2 y^5\ge0$ for $y\ge0$ and $h(x,0)=2x^5\ge0$ for $x\ge0$. This does it for the boundary pieces $\{0\}\times[0,\infty)$ and  $[0,\infty)\times\{0\}$. 

By symmetry, it remains to show that $\liminf_{x,y\to\infty}h(x,y)\ge0$ and $\liminf_{x\to0,y\to\infty}h(x,y)\ge0$. 
Let 
$$H_{11}(x,y):=h(\tfrac1x,\tfrac1y)x^{10}y^{10}
=2+2 x^{10} y^5+2 x^9 y^6+2 x^8 y^4-3 x^7 y^6+2 x^7 y^5+2 x^6 y^9-3 x^6
   y^7-3 x^6 y^4+2 x^5 y^{10}+2 x^5 y^7-6 x^5 y^5+2 x^4 y^8-3 x^4 y^6-3
   x^4 y^3-3 x^3 y^4+2 x^3 y^3+2 x^2 y+2 x y^2.$$
Obviously, $H_{11}(0+,0+)=2>0$. So, $\liminf_{x,y\to\infty}h(x,y)\ge0$. 

Let 
$$H_{01}(x,y):=h(\tfrac1x,y)x^{10}
=2 x^{10} y^5+2 x^9 y^4+2 x^8 y^6+2 x^7 y^5-3 x^7 y^4-3 x^6 y^6-3 x^6
   y^3+2 x^6 y-6 x^5 y^5+2 x^5 y^3+2 x^5-3 x^4 y^7-3 x^4 y^4+2 x^4
   y^2+2 x^3 y^7-3 x^3 y^6+2 x^2 y^9+2 x y^8+2 y^{10}.$$
Removing from the latter expression terms dominated by other terms whenever $x,y\to0$, we end up with showing that $\liminf_{x\to0,y\to\infty}\tH(x,y)\ge0$, where 
$$\tH(x,y):=(2+o(1)) x^5 + (2+o(1)) x^4 y^2 - (3+o(1)) x^3 y^6 + 2 x y^8 + 2 x^2 y^9 + 2 y^{10},$$
which follows because 
$$(2+o(1)) x^4 y^2 - (3+o(1)) x^3 y^6 + 2 x y^8 \\
\ge(2+o(1)) x^5 y^4 - (3+o(1)) x^3 y^6 + 2 x y^8 \\ 
=x y^8[(2+o(1)) s^4 - (3+o(1)) s^2 + 2 ]>0$$
for all small enough $x,y>0$, where $s:=x/y$. 
So, $\liminf_{x\to0,y\to\infty}h(x,y)\ge0$. $\quad\Box$




  [1]: https://i.sstatic.net/sC1bd.png