Let us fix first some notations:
$f=ab+ac+bc-1$
$g=\sqrt{5a+5b+8ab}+\sqrt{5a+5c+8ac}+\sqrt{5b+5b+8bc}$
$\mathcal{D}=\{(a,b,c)\in\mathbb{R}^3,a\geq0,b\geq0,c\geq0\}$
$\mathcal{D}_f=\{(a,b,c)\in\mathcal{D},f=0\}$
$\mathcal{C}=\{(a,b,c)\in\mathbb{R}^3,0\leq a \leq 20,0\leq b \leq 20,0\leq c \leq 20\}$
$\mathcal{C}_f=\{(a,b,c)\in\mathcal{D},f=0\}$
For any triple $(a,b,c)\in\mathcal{D}$ which is not in $\mathcal{C}$, one of the real numbers $a$, $b$ and $c$ is greater than 20 and obviously we have:
$g\geq \sqrt{5a}+\sqrt{5b}+\sqrt{5c}>10>3\sqrt{2}+2\sqrt{5}$.
We are left to proving the inequality for $(a,b,c)\in\mathcal{C}_f$. Since $g$ is continuous and differentiable function of the variables $a$, $b$ and $c$ and since $\mathcal{C}_f$ is a compact subset of $\mathbb{R}^3$, there is a minimum for $g$ on $\mathcal{C}_f$ which is attained for at least one point.
If the point where minimum is achieved lies on the boundary, one of its coordinates is $0$ or $20$ (in which case g is at least 10). Let us suppose we have a coordinate equal to $0$, without loss of generality (by symmetry), we can assume $c=0$. When then need to have $ab=1$, so $a$ is not $0$ and $b=1/a$ which yields:
$g = \sqrt{5a}+\sqrt{5/a}+\sqrt{5a+5/a+8}$
It is easily determined that the minimum of this function for $a>0$ is met when $a=1$ and thus $g=3\sqrt{2}+2\sqrt{5}$.
We are left with points on the interior of $\mathcal{C}_f$. If the minimum is met in the interior, it has to be a local extremum. We will rely on the method of Lagrange Multiplier to search for extremal points. Such a point must verify that the gradients of $f$ and $g$ (as function of $a$, $b$ and $c$) are colinear. These gradients can be computed:
$\mathrm{Grad}(f)=(b+c,a+c,b+c)$
$\mathrm{Grad}(g)=(\frac{8b+5}{2r_1}+\frac{8c+5}{2r_2},\frac{8a+5}{2r_1}+\frac{8c+5}{2r_3},\frac{8a+5}{2r_2}+\frac{8b+5}{2r_3})$, where
$r_1 = \sqrt{5a+5b+8ab}$
$r_2 = \sqrt{5a+5c+8ac}$
$r_3 = \sqrt{5b+5c+8bc}$
The fact that both gradients are colinear translate into their cross product being null, which yields three equations (one per coordinate of the cross product). These equation involve square roots, so we apply the classical trick to multiply each one with its seven conjugates (i.e. applying all possible sign choices for the square roots). Last as noted by Claude Lebovici, one has $c=(1-ab)/(a+b)$ and we can substitute this value in the equations. All in all, we obtain three polynomial equations in $a$ and $b$ by taking the numerators. One can note that $a-b$ appears as a factor in one equation while $a^2+2ab-1$ appears in two others equations and corresponds to case where $b=c$. The other factors imply too many terms to be copied here but are provided by any algebra package on a modern computer.
For each pair of equations, using a resultant, one can eliminate the variable $b$ to get a polynomial in $a$. We can then take a polynomial gcd between the three univariate polynomial to conclude that $p(a)=0$ if $(a,b,c)$ is an extremum of $g$ on $\mathcal{C}_f$. When factoring over $\mathbb{Q}[a]$, the polynomial $p$ has 11 irreducible factors with degrees ranging from 1 to 5. We can then inject the factors (one after each other) in the three equations and a tedious case study (keeping only positive real numbers) identify the following local extrema:
$a=b=c=1/\sqrt{3}$
$a=b$ are roots of the degree 5 five polynomial $512x^5+320x^4-300x^3-100x-125$ and $c=(1-ab)/(a+b)$
The extrema found are higher that the minimum previously identified. This concludes the identification of the potential minima of $g$ on $\mathcal{C}_f$. The inequality is thus proved since it corresponds to the minimum of $g$ on $\mathcal{C}$.
As a further note, this proof is very computational and not elegant. However, I believe the method can be applied to the other questions you submitted on this site.
Problem. Let $a, b, c$ be positive real numbers. Prove that:
$$\frac{a}{\sqrt{a^2+8bc}}+\frac{b}{\sqrt{b^2+8ca}}+\frac{c}{\sqrt{c^2+8ab}}<2.$$
Proof. It suffices to prove that
$$\frac{a}{\sqrt{a^2+2bc}}+\frac{b}{\sqrt{b^2+2ca}}+\frac{c}{\sqrt{c^2+2ab}}<2.$$
Let $x := \frac{2bc}{a^2}, y := \frac{2ca}{b^2}, z := \frac{2ab}{c^2}$. Then $xyz = 8$. We need to prove that
$$\frac{1}{\sqrt{1 + x}} + \frac{1}{\sqrt{1 + y}} + \frac{1}{\sqrt{1 + z}} < 2.$$
WLOG, assume that $x \le y \le z$.
We split into two cases:
Case 1. $x + y \ge 6$
We have $y \ge 3$. Thus,
$$\mathrm{LHS} < 1 + \frac{1}{\sqrt{1 + 3}} + \frac{1}{\sqrt{1 + 3}} = 2.$$
Case 2. $x + y < 6$
Using AM-GM, we have
$$\frac{1}{\sqrt{1 + x}} = \frac{2\sqrt{1 + x}}{2(1 + x)} \le \frac{(1 + x) + 1}{2(1 + x)} = \frac{x + 2}{2 + 2x}$$
and
$$\frac{1}{\sqrt{1 + y}} = \frac{2\sqrt{1 + y}}{1 + y} \le \frac{(1 + y) + 1}{2(1 + y)} = \frac{y + 2}{2 + 2y}.$$
It suffices to prove that
$$\frac{x + 2}{2 + 2x} + \frac{y + 2}{2 + 2y} + \frac{1}{\sqrt{1 + z}} < 2$$
or
$$\frac{x}{2 + 2x} + \frac{y}{2 + 2y} > \frac{1}{\sqrt{1 + z}}.$$
Using AM-GM, it suffices to prove that
$$ 2\sqrt{\frac{x}{2 + 2x} \cdot \frac{y}{2 + 2y}} > \frac{1}{\sqrt{1 + z}}$$
or
$$4\cdot \frac{x}{2 + 2x} \cdot \frac{y}{2 + 2y} > \frac{1}{1 + z}$$
or (using $z = \frac{8}{xy}$)
$$\frac{xy(7 - x - y)}{(1 + x)(1 + y)(xy + 8)} > 0$$
which is true.
We are done.
Best Answer
Denote the expression by $f(a, b, c)$.
Using $2uv \le u^2 + v^2$, we have \begin{align*} &a\sqrt{a^2+2bc}+b\sqrt{b^2+2ca}+c\sqrt{c^2+2ab}\\ \le{}& a\sqrt{a^2+b^2 + c^2}+b\sqrt{b^2+c^2 + a^2}+c\sqrt{c^2+a^2 + b^2}\\ ={}& (a + b + c)\sqrt{a^2 + b^2 + c^2}. \tag{1} \end{align*}
Also, by AM-GM, we have $$(a + b + c)(ab + bc + ca) \ge 3\sqrt[3]{abc}\cdot 3\sqrt[3]{ab \cdot bc \cdot ca} = 9abc. \tag{2}$$
Using (1) and (2), letting $x = ab + bc + ca \in [0, 1/9]$ (using $(a + b + c)^2 \ge 3(ab + bc + ca)$), we have \begin{align*} f(a, b, c) &\le 3(a + b + c)(ab + bc + ca) + (a + b + c)\sqrt{a^2 + b^2 + c^2}\\ &= \sqrt 3\, (ab + bc + ca) + \frac{1}{\sqrt 3} \sqrt{\frac13 - 2(ab + bc + ca)}\\ &= \sqrt 3\, x + \frac{1}{\sqrt 3} \sqrt{\frac13 - 2x}\\ &\le \frac{2\sqrt 3}{9} \end{align*} where we use $a^2 + b^2 + c^2 = (a + b + c)^2 - 2(ab + bc + ca) = 1/3 - 2(ab + bc + ca)$, and $$\left(\frac{2\sqrt 3}{9} - \sqrt 3\, x\right)^2 - \left(\frac{1}{\sqrt 3} \sqrt{\frac13 - 2x}\right)^2 = \frac{1}{27}(9x - 1)^2 \ge 0.$$
Also, when $a = b = c = \frac{1}{3\sqrt 3}$, we have $f(a, b, c) = \frac{2\sqrt 3}{9}$.
Thus, the maximum of $f(a,b,c)$ is $\frac{2\sqrt 3}{9}$.