Here is a solution with some Matlab help for the analysis, but with a clear manual proof path.
Due to homogeneity, we can demand $a^2+ b^2 + c^2 = 1$. Define $m$ to be the mean of $a,b,c$, i.e. $a + b+c = 3m$. Then note that
$$
9 m^2 = (a+b+c)^2 = a^2+ b^2 + c^2 + 2 (ab + bc + ca) = 1 + 2(ab + bc + ca)
$$
Hence the claim can be written
$$
\sum_{cyc} \sqrt{\frac{a}{b+c}+\frac{b}{c+a}}\ge 2+\sqrt{\frac{2}{9m^2-1}}
$$
Now turn to the LHS. Write $a = m +x$, $b = m + y$, $c = m+z$ with $x+y+z=0$ and $1 = a^2 + b^2 + c^2 = 3 m^2 + x^2 + y^2 + z^2$ which gives two conditions for $(x,y,z)$. W.l.o.g. $(x,y,z)$ can then be expressed as
$$
x = \sqrt\frac23 \sqrt{1 - 3m^2}\cos(\phi-2\pi/3)\\
y = \sqrt\frac23 \sqrt{1 - 3m^2}\cos(\phi-4\pi/3)\\
z = \sqrt\frac23 \sqrt{1 - 3m^2}\cos(\phi)
$$
Hence the claim can be written, with these $(x,y,z)$, as
$$
\sum_{cyc} \sqrt{\frac{m+x}{2m-x}+\frac{m+y}{2m-y}}\ge 2+\sqrt{\frac{2}{9m^2-1}}
$$
The LHS is now a function of $\phi$ whereas the RHS is not. For any $m$, a free (unbounded) minimum w.r.t. $\phi$ of the LHS occurs at $\phi = \pi$ which can be shown by varying $\phi$ about $\pi$. [For bounded minima see below.] So we have to inspect the LHS at that minimum and show that
$$
\lim_{(\phi = \pi)} \sum_{cyc} \sqrt{\frac{m+x}{2m-x}+\frac{m+y}{2m-y}}- 2-\sqrt{\frac{2}{9m^2-1}} \ge 0
$$
Since $(a,b,c)$ should be nonnegative, this requires that $c = m + z = m - \sqrt\frac23 \sqrt{1 - 3m^2} > 0$ or $m > \sqrt2 / 3$, this bound corresponds to $(a,b,c) = (\frac{1}{\sqrt2},\frac{1}{\sqrt2},0)$. On the other hand, the maximum possible $m$ occurs when $a = b = c = m$ or, since $a^2+b^2 + c^2 = 1$, at $m = 1/\sqrt3$.
Let's look at the two extreme values for $m$. Indeed we have (using Matlab) that
$$
\lim_{(m = \sqrt2 / 3)} \lim_{(\phi = \pi)} \sum_{cyc} \sqrt{\frac{m+x}{2m-x}+\frac{m+y}{2m-y}}- 2-\sqrt{\frac{2}{9m^2-1}} = 0 \\
\lim_{(m = 1 / \sqrt3)} \lim_{(\phi = \pi)} \sum_{cyc} \sqrt{\frac{m+x}{2m-x}+\frac{m+y}{2m-y}}- 2-\sqrt{\frac{2}{9m^2-1}} = 0
$$
and for all values of $m$ in between the $> 0 $ holds. Below is a plot which illustrates this.
The minimum of the LHS may as well be bounded by the fact that $(a,b,c)$ should be nonnegative. In that case, the bound arises when the smallest variable, say $c$, is zero, and it must be inspected, while keeping $c=0$, until another variable becomes zero. So that bound is given by $0 = c = m + \sqrt\frac23 \sqrt{1 - 3m^2}\cos(\phi)$ or $m = \sqrt{\frac{\frac23 \cos^2(\phi) }{1 + 2 \cos^2(\phi) }}$ and must be inspected for $\frac43 \pi > \phi > \frac23 \pi$ since at $\frac23 \pi$ (or $\frac43 \pi$ ) we have that also $b =0$ or $a =0$ (then the terms diverge, and this case was excluded by the OP). That means we have to look at (with the $(x,y,z)$ as above)
$$
\lim_{m = \sqrt{\frac{\frac23 \cos^2(\phi) }{1 + 2 \cos^2(\phi) }}} \sum_{cyc} \sqrt{\frac{m+x}{2m-x}+\frac{m+y}{2m-y}}- 2-\sqrt{\frac{2}{9m^2-1}}
$$
which is a function of $\phi$. Variation of $\phi$ about $\pi$ already shows local positivity. Here is a plot (where $\phi$ was denoted $x$) which illustrates the overall behavior:
This proves the claim. $\qquad \Box$
