I guess this is not possible in term of standard mathematical functions. Is it solvable at all?
Edit: I think it's possible (by proper substitution) to do this using the Fermi-Dirac integral. But that would change the lower limit and I'd need to numerically solve the incomplete Fermi-Dirac integral for which I haven't found a working algorithm yet. So I'd prefer to keep the lower limit at 0
Assuming $c>0$ we have
$$ f(c)=\int_{0}^{+\infty}\frac{\sqrt{t+c}}{e^t+1}\,dt=\int_{1}^{+\infty}\frac{\sqrt{c+\log u}}{u(u+1)}\,du=\int_{0}^{1}\frac{\sqrt{c-\log v}}{v+1}\,dv $$ and the last integral is simple to approximate numerically.
By the Cauchy-Schwarz inequality he have $$ f(c)\leq \sqrt{\int_{0}^{1}(c-\log v)\,dv\int_{0}^{1}\frac{dv}{(v+1)^2}}\leq\sqrt{\frac{c+1}{2}}$$ and this bound is both simple and pretty accurate for values of $c$ not much larger than $10$.
Additionally, by the integral representation for the $\eta$ (or $\zeta$) function we have $$ \lim_{c\to 0^+}f(c) =\frac{2-\sqrt{2}}{4}\sqrt{\pi}\,\zeta\left(\tfrac{3}{2}\right) = 0.6780938951531\ldots$$ and by the (inverse) Laplace transform $$ f(c) = \sum_{n\geq 1}\frac{2 \sqrt{c n}+e^{c n} \sqrt{\pi } \operatorname{Erfc}\sqrt{c n}}{2 n^{3/2}}(-1)^{n+1}. $$ For an alternative and more accurate numerical approximation, one may replace the complementary error function by a convergent of its well-known continued fraction.
Addendum: by this tail inequality for the normal distribution proved on MO, for large values of $c$ we get
$$ f(c) \approx \sum_{n\geq 1}\frac{2 \sqrt{c n}+2\frac{1}{\sqrt{cn}+\sqrt{cn+2}}}{2 n^{3/2}}(-1)^{n+1}=\frac{1}{2}\sum_{n\geq 1}\frac{\sqrt{cn+2}-\sqrt{cn}}{n\sqrt{n}}(-1)^{n+1}. $$