Using
$$ K(ik) =
\frac{1}{\sqrt{1+k^2}}K\left(\sqrt{\frac{k^2}{k^2+1}}\right) $$
and a substitution $t^2 = \frac{k^2}{1+k^2}$, rewrite the integral as
$$ \int_0^\infty K(i k)^3\,dk = \int_0^1 K(t)^3\,dt. $$
There is a paper "Moments of elliptic integrals and critical L-values"
by Rogers, Wan and Zucker (http://arxiv.org/abs/1303.2259; also one of
the authors' earlier papers: http://arxiv.org/abs/1101.1132), and the
authors, by relating this integral to an L-series of a modular form (their theorems 1 and 2),
show that
$$ \int_0^1 K(k)^3\,dk = \frac{3}{5}K(1/\sqrt{2})^4 =
\frac{3\Gamma(\frac14)^8}{1280\pi^2}, $$
using $K(1/\sqrt{2}) = \frac14 \pi^{-1/2}\Gamma(\frac14)^2$.
Let us start with a warm-up exercise. Introduce the functions
$$g_{\pm}(x)=\operatorname{Ai}(x)\pm i\operatorname{Bi}(x).$$
Computing the Wronskian of these two solutions of the Airy equation, one can check that
$$\frac{1}{\operatorname{Ai}^2(x)+\operatorname{Bi}^2(x)}=\frac{\pi}{2i}\left[\frac{g_+'(x)}{g_+(x)}-\frac{g'_-(x)}{g_-(x)}\right]$$
This gives the integral $\mathcal{K}(0)$ as
$$\mathcal{K}(0)=\pi\left[\arg g_+(\infty)-\arg g_+(0)\right]=\pi\left[\pi-\frac{\pi}{3}\right]=\frac{\pi^2}{6}.$$
To compute the integral $\mathcal{K}(3n)$, we will need to develop a more sophisticated approach. First note that (see here)
$$g_{\pm}(x)=-2e^{\mp 2\pi i/3}\operatorname{Ai}\left(e^{\mp2\pi i/3}x\right).$$
Therefore
\begin{align}\mathcal{K}(3n)&=\frac{\pi}{2i}\int_0^{\infty}x^{3n}\left[\frac{g_+'(x)}{g_+(x)}-\frac{g'_-(x)}{g_-(x)}\right]dx=\\
&=\frac{\pi}{2i}\lim_{R\rightarrow\infty}\int_{S_R}z^{3n}\frac{\operatorname{Ai}'(z)}{\operatorname{Ai}(z)}dz,
\end{align}
where the contour $S_R$ in the complex $z$-plane
is composed of two segments: one going from $Re^{2\pi i/3}$ to $ 0$ and another one going from $0$ to $ Re^{-2\pi i/3}$.
It is a well-known fact that the Airy function $\operatorname{Ai}(z)$ has zeros (i.e. our integrand has poles) on the negative real axis only. Therefore by residue theorem our integral is equal to
$$\mathcal{K}(3n)=-\frac{\pi}{2i}\lim_{R\rightarrow \infty}\int_{C_R}z^{3n}\left[\ln\operatorname{Ai}(z)\right]'dz,\tag{1}$$
where $C_R$ is the arc of the circle of radius $R$ centered at the origin going counterclockwise from $Re^{-2\pi i/3}$ to $Re^{2\pi i/3}$.
The limit (1), on the other hand, can be computed using the asymptotics of the Airy function as $z\rightarrow\infty$ for $|\arg z|<\pi$:
\begin{align}
\ln\operatorname{Ai}(z)\sim -\frac23 z^{3/2}-\ln2\sqrt{\pi}-\frac14\ln z+
\ln\sum_{k=0}^{\infty}\frac{(-1)^k\left(\frac16\right)_k\left(\frac56\right)_k}{k!}\left(\frac43 z^{3/2}\right)^{-k}.\tag{2}
\end{align}
Note that if we introduce instead of $z$ the variable $s=\frac43z^{3/2}$, then the integration will be done over the circle of radius $\Lambda=\frac43 R^{3/2}$, i.e. a closed contour in the complex $s$-plane. The corresponding integral can therefore be computed by residues by picking the necessary term in the large $s$ expansion of $\ln \operatorname{Ai}(z)$.
More precisely, we have the following formula:
\begin{align}
\mathcal{K}(3n)=-\frac{\pi}{2i}\lim_{\Lambda\rightarrow \infty}\oint_{|s|=\Lambda}
\left(\frac{3s}{4}\right)^{2n}d\left[-\frac16\ln s+\ln\sum_{k=0}^{\infty}\frac{(-1)^k\left(\frac16\right)_k\left(\frac56\right)_k}{k!}s^{-k}\right]\tag{3}
\end{align}
To compute the residue, it suffices to expand the logarithm-sum up to order $2n$ in $s^{-1}$. Note that the Pochhammer symbol coefficients are in fact some rational numbers.
In the simplest case $n=0$, the residue is determined by the term $-\frac16\ln s$ and we readily reproduce the previous result
$$\mathcal{K}(0)=-\frac{\pi}{2i}\cdot 2\pi i\cdot\left(-\frac16\right)=\frac{\pi^2}{6}.$$
The general formula for arbitrary $n$ would look a bit complicated (but straightforward to obtain) due to the need to expand the logarithm of the sum.
Example. The calculation of the corresponding values $M(n)=\mathcal{K}(3n)$ in Mathematica can be done using the command
\begin{align}\mathtt{\text{ M[n_] := -Pi^2 SeriesCoefficient[
Series[(3 s/4)^(2 n) D[-Log[s]/6 +}} \\ \mathtt{\text{ Log[Sum[(-1)^k
Pochhammer[1/6, k] Pochhammer[5/6, k]/(k! s^k), }} \\
\mathtt{\text{{k, 0, 2 n}]], s], {s, Infinity, 1}], 1]}} \end{align}
This yields, for instance,
$$M(0)=\frac{\pi^2}{6},\quad M(1)=\frac{5\pi^2}{32},\quad M(2)=\frac{565\pi^2}{512},$$ $$\ldots, M(10)=\frac{2\,660\,774\,144\,147\,177\,521\,025\,228\,125\,\pi^2}{2\,199\,023\,255\,552},\ldots$$
and so on.
Added: The large $s$ expansion (2) can also be found directly by using Airy equation. Moreover, by transforming it into an equation for $\ln \operatorname{Ai}(z)$, one can avoid reexpanding the logarithm of the sum. The price to pay will be that the expansion coefficients will be determined by a nonlinear recurrence relation instead of explicit formulas.
Best Answer
Here's one way to go.
First, note that $$\begin{eqnarray*} \int_0^\pi\frac{3\cos x+\sqrt{8+\cos^2 x}}{\sin x}x\ \mathrm dx &=& \int_0^\pi\frac{3x(1+\cos x)}{\sin x} \mathrm dx +\int_0^\pi\frac{3x}{\sin x} \left(-1+\sqrt{1-\frac{\sin^2x}{9}}\right)\ \mathrm dx. \end{eqnarray*}$$ For now I'll simply claim that \begin{equation*} \int_0^\pi\frac{3x(1+\cos x)}{\sin x} \mathrm dx = \pi\log 64.\tag{1} \end{equation*} (I would be surprised if this integral has not been handled somewhere on this site.) But $$\begin{eqnarray*} \int_0^\pi\frac{3x}{\sin x} \left(-1+\sqrt{1-\frac{\sin^2x}{9}}\right)\ \mathrm dx &=& \int_0^\pi\frac{3x}{\sin x} \sum_{k=1}^\infty {1/2\choose k} \frac{(-1)^k}{3^{2k}} \sin^{2k}x \ \mathrm dx \\ &=& \sum_{k=1}^\infty {1/2\choose k} \frac{(-1)^k}{3^{2k-1}} \int_0^\pi x \sin^{2k-1}x \ \mathrm dx \\ &=& \sum_{k=1}^\infty {1/2\choose k} \frac{(-1)^k}{3^{2k-1}} \frac{\pi^{3/2}\Gamma(k)}{2\Gamma(k+1/2)} \\ &=& -\pi \sum_{k=1}^\infty \frac{1}{3^{2k-1}2k(2k-1)} \\ &=& -\pi \log \frac{32}{27}. \end{eqnarray*}$$ (The last sum can be found by standard methods. Schematically, $\sum \frac{a^{2k-1}}{2k(2k-1)} = \sum \int {\mathrm da} \frac{a^{2k-2}}{2k}$.) Thus, the integral is $\pi \log 54$ as claimed.
Proof of (1): We have $$\begin{eqnarray*} \int_0^\pi \frac{3x(1+\cos x)}{\sin x} \ \mathrm dx &=& \int_{0^+}^\pi \frac{3x(1+\cos x)}{\sin x} \ \mathrm dx \\ &=& 3\int_{0^+}^\pi x \cot\frac{x}{2} \ \mathrm dx \hspace{5ex}\textrm{(double angle formulas)} \\ &=& 12 \int_{0^+}^{\pi/2} t\cot t \ \mathrm dt \hspace{5ex} (t = x/2) \\ &=& -12\int_{0^+}^{\pi/2} \log\sin t \ \mathrm dt \hspace{5ex}\textrm{(integrate by parts)} \\ &=& -6\int_{0^+}^{\pi/2} \log\sin^2 t \ \mathrm dt \\ &=& -6\int_{0^+}^{\pi/2} \log(1-\cos^2 t) \ \mathrm dt \\ &=& 6 \int_{0^+}^{\pi/2} \sum_{k=1}^\infty \frac{1}{k}\cos^{2k}t \ \mathrm dt \hspace{5ex}\textrm{(series for log)} \\ &=& 6\sum_{k=1}^\infty \frac{1}{k} \int_{0^+}^{\pi/2} \cos^{2k}t \ \mathrm dt \hspace{5ex} \textrm{(Tonelli's theorem)}\\ &=& 6\sum_{k=1}^\infty \frac{1}{k} \frac{\sqrt{\pi}\Gamma(k+1/2)}{2\Gamma(k+1)} \\ &=& 3\pi \sum_{k=1}^\infty {1/2 \choose k}(-1)^{k+1}\frac{2k-1}{k} \\ &=& \pi \log 64. \end{eqnarray*}$$ Note that $$\begin{eqnarray*} 6\pi \sum_{k=1}^\infty {1/2 \choose k}(-1)^{k+1} &=& -6\pi \left[\sum_{k=0}^\infty {1/2 \choose k}(-1)^{k} - 1\right] \\ &=& -6\pi[(1-1)^{1/2} - 1] \\ &=& 6\pi \end{eqnarray*}$$ and $$\begin{eqnarray*} -3\pi \sum_{k=1}^\infty {1/2 \choose k}(-1)^{k+1} \frac{1}{k} &=& 3\pi \sum_{k=1}^\infty {1/2\choose k}(-1)^k \int_0^1 x^{k-1} \ \mathrm dx \\ &=& 3\pi \int_0^1 \frac{1}{x} \left[ \sum_{k=0}^\infty {1/2\choose k}(-1)^k x^{k} -1 \right] \ \mathrm dx \\ &=& 3\pi \int_0^1 \frac{1}{x} \left( \sqrt{1-x} -1 \right) \ \mathrm dx \\ &=& 3\pi(-2+\log 4) \\ &=& -6\pi + \pi\log 64. \end{eqnarray*}$$