Let $\Omega \subseteq \mathbb{R}^3$ be the complement of a closed ball $\overline{B(0,R)}$, with radius $R > 0$, and $u \in C^2(\Omega) \cap C^0(\overline{\Omega})$ a harmonic function, such that the integral (i.e. corresponding "energy" on the domain $\Omega$) $$\int_\Omega \left\lVert\nabla u\right\rVert^2 \, \mathrm{d}V$$ is finite.
Is limit $\lim_{\left\lVert x \right\rVert \to \infty} u(x)$ necessarily well-defined and finite?
$\newcommand{\R}{\mathbb{R}}$ $\newcommand{\intt}{\int\limits}$
The answer to your question is yes, and even with some decent polynomial bound. My proof works for all dimensions $d\ge 3$, for $d = 2$ I didn't think too hard but I think the answer is negative (at least, if we replace complement of a disk by a half-plane $\Re z > 0$ then there is a counterexample $u(z) = \log (\log (z + 10))$.)
First, let me change the notation a bit. Number $R$ is more or less irrelevant for us since we are interested in the limit as $||x||\to \infty$ so I will leave the radius of the initial ball unspecified. Also let me change $u$ to $v$ and denote $u:=\nabla v$, since this is the function we will be mostly working with and I prefer my main harmonic function to be called $u$ (note that $u:\R^3 \to \R^3$!).
Our starting point is the following Lemma
$\textbf{Lemma 1}$ Let $u$ be a harmonic function on $B(0, 1)$ and $\int_{B(0, 1)} ||u(x)||^2 dx \le 1$. Then for all $x\in B(0, 1/2)$ we have $|u(x)| \le C$ where $C$ is some absolute constant.
This lemma can be proved in a hundred different ways, for example by considering a ball of radius $1/2$ around $x$ and using mean value theorem and the Cauchy-Schwartz inequality.
From this lemma by a simple rescaling we can prove the following Corollary
$\textbf{Corollary 1}$ Let $u$ be a harmonic function on $B(p, R)$ and $\int_{B(p, R)} ||u(x)||^2 dx \le A$. Then for all $x\in B(p, R/2)$ we have $|u(x)| \le C\frac{A^{1/2}}{R^{3/2}}$ where $C$ is some absolute constant.
With this corollary at hand we are ready to go. To prove the existence of the limit we will prove that our sequence is Cauchy. That is, we pick two points $x, y\in \R^3$ with $\min (||x||, ||y||)\to \infty$ and we want to show that $|v(x) - v(y)|\to 0$.
Let $\gamma$ be any piecewise smooth curve connecting $x$ to $y$. Then we have $|v(x) - v(y)| \le \int_\gamma ||u(t)||dt$.
Let $||x|| = R_1$ and $||y|| = R_2$. Without loss of generality we assume that $R_1 \le R_2$. Denote by $z$ the intersection of the sphere $S(0, R_1)$ and the segment $[0y]$. Let us connect $x$ and $z$ by a piecewise linear curve $x=p_1, p_2, p_3, \ldots p_{n-1}, p_n=z$ such that $p_k\in S(0, R_1)$ for all $k$ and the distance between any two consecutive points is at most $R_1/10$. It's easy to see that this can be achieved with $n\le 100$.
So, our curve $\gamma$ will be the broken line $p_1 \ldots p_n y$. The problem is now separated into two quite different ones: a) estimate the difference $|v(p_k) - v(p_{k+1})|$ for $k = 1, \ldots, n-1$ and b) estimate $|v(z) - v(y)|$. We begin with the first one
So, let $p, q\in S(0, R_1)$ with $dist(p, q)\le R_1/10$. Consider the ball $B = B(p, R_1/5)$. If $R_1$ is big enough then it lies entirely inside of the support of $u$ and $v$ (that is, it doesn't intersect the removed ball from the OP). We have $$\intt_B ||u(s)||^2ds \le A,$$ where $A$ is the value of the integral which we assumed is finte, therefore for all $t\in [pq]$ we have $|u(t)|\le C\frac{\sqrt{A}}{R_1^{3/2}}$. Integrating this over the interval $[pq]$, whose length is at most $R_1/10$, we get the final bound $|v(p)-v(q)| \le \frac{C_1}{R_1^{1/2}}$. Summing these estimates (as we said, $n\le 100$) we get $|v(x)-v(z)|\le \frac{100C_1}{R_1^{1/2}}$.
Note that with a bit more care this would work even for $d = 2$. It hints that the second part would require a new idea. And indeed, we will use a good old dyadic decomposition.
So, let us pick the points $z = q_1, q_2, \ldots, q_m = y$ on the interval $[zy]$ such that $\frac{||q_{k+1}||}{||q_k||}=1.01$ for all $k = 1, \ldots, m-2$ and $1\le \frac{||q_m||}{||q_{m-1}||}\le 1.01$.
Essentially repeating the above argument we can get that $|v(q_{k+1})-v(q_k)|\le \frac{C_1}{||q_k||^{1/2}}$. Finally, we note that $\frac{C_1}{||q_k||^{1/2}}$ is a decreasing geometric progression which is known to converge to something constant away from the first term, so summing these estimates up we get $$|v(z) - v(y)| \le \frac{C_2}{||z||^{1/2}} = \frac{C_2}{R_1^{1/2}}.$$
Doing one last addition we get $|v(x) - v(y)| \le \frac{100C_1 + C_2}{R_1^{1/2}},$ which tends to zero as $\min(||x||, ||y||) = R_1$ tends to infinity.