Solution 1:

Let us try to model the grid:

hex grid 1

We have six base directions $$ u_k = (x_k, y_k) = d (\cos k \pi/3, \sin k \pi/3) \quad (k\in \{0,\ldots,5\}) $$ where $d$ is the incircle diameter of a hexagon cell. Starting from the origin $(0,0)$ each cell can be reached by walking along one of the six directions, e.g. $$ u = u_2 + 2u_5 - u_1 + 7u_2 + 2 u3 + \cdots $$ This shrinks down to the integer multiples of two directions, because $u_3 =-u_0$ etc and $u_2 = u_1-u_0$: $$ u = c_0 u_0 + c_1 u_1 $$ so we can use $(c_0,c_1)$ as coordinates for each cell.

The vertices of the center cell are $$ v_k = D (\cos (2k+1) \pi/6, \sin (2k+1) \pi/6) \quad (k\in \{0,\ldots,5\}) $$ where $D$ is the diameter of the circumscribed circle. It relates to $d$ via: $$ D = \frac{2}{\sqrt{3}} d $$

We can describe every vertex via $(c_0,c_1,k)$, which is not injective.

Each vertex within the circle fulfills $$ \lVert u + v_k \rVert \le R $$

Here is an image where the red circles (except for the one around the center) have radius $r_c = (2k + 1) r$.

hex grid 2

Here the hexagons were coloured the same, when they have the same circumscribed circle (these are the easy to calculate cases):

hex grid 3

The number of black circles between red circles increases with increasing radius. This is because there are more discrete possibilities to position the hexagons.

hex grid 4

So it is possibly a good idea to let a machine do the counting, it seems to get only more complicated, with the limit of continous radius possibility for large $r$.

hex grid 5

Please take the above graphics with a grain of salt: The solid graph is $f(r) = \frac{A(r)}{A_H} = \frac{\pi r^2}{2 \sqrt{3}}$ which is the area of the circumscribed circle expressed in hexagon cells. Due to the gaps at the rim, the number of hexagons (dots) is less than this ideal number. However the dots with the hexagon count for that radius seem too low, I need to hunt the bug and provide a version with corrected numbers.

Solution 2:

Posting the trivial estimates so that we can measure progress.

Assume that the side of the hexagon has length $1$, and that the radius of the circle is $r$. The area of a single hexagon is then $6\cdot\sqrt3/4=3\sqrt3/2$. Let $N$ be the number of hexagons inside the circle. Their total area is less than that of the circle, so we get the inequality $$ \frac{3\sqrt3 N}2<\pi r^2. $$ On the other hand the hexagons cover the concentric circle of radius $r-2$ completely. Therefore we get another inequality $$ \frac{3\sqrt3 N}2>\pi (r-2)^2. $$ The midway point gives then the ballpark estimate of $$ N\approx\frac{2\sqrt3\pi(r-1)^2}9 $$ together with the bound on the absolute value of the error $$ |\Delta N|<\frac{4\sqrt3\pi r}9. $$ The link by Chris Culter suggests that it may be possible to get an error term of order $r^{2/3}$. The problem considered there is about the number of lattice points, so it is not exactly equivalent.