Your answer of $Na^2$ is correct. To get there, one can notice that this is equivalent to the following:
What is the expected value of $|A|^2$ if $A$ is the sum of $n$ points uniformly distributed on the unit circle?
This follows from "complexifying" the problem. In particular, the harmonic addition formula may be derived by noting that
$$\sin(x)=\text{Im}(e^{ix})$$
where $\text{Im}$ is the imaginary part of a number. Note that, since $e^{ix}$ traces out a circle in the complex plane and $\text{Im}(x)$ is a projection onto the imaginary axis, the geometric interpretation of this is that if you move an angle of $x$ away from one axis on a unit circle, then project onto the other axis, you get $\sin(x)$.
Let us ignore $k$ and $a$ for now, as they just act to scale the problem. If we take a sum of $\sin(x+c_k)$ for various $c_k$ uniformly distributed in $[0,2\pi]$, we get
$$\sum_{k=1}^{N}\sin(x+c_k)=\sum_{k=1}^N\text{Im}(e^{i(x+c_k)})=\sum_{k=1}^n\text{Im}(e^{ic_k}e^{ix})=\text{Im}\left(e^{ix}\cdot \sum_{k=1}^N e^{ic_k}\right)$$
Then, if we write things in terms of sines again, we get that if $z= \sum_{k=1}^N e^{ic_k}$, then the sum evaluates to $|z|\sin(x+\text{arg}(z))$, where $\text{arg}(z)$ is the argument of $z$.
To calculate the expected value of $|z|$, we first rewrite it as $z\overline{z}$. Then, noting that $\overline{e^{ic_k}}=e^{-ic_k}$, we expand $\mathbb E(|z|^2)$ to the following:
$$\mathbb E(z\overline{z})=\mathbb E\left(\left(\sum_{k=1}^Ne^{ic_k}\right)\left(\sum_{j=1}^Ne^{-ic_j})\right)\right)=\sum_{k=1}^N\sum_{j=1}^N\mathbb E(e^{i(c_k-c_j)})$$
Then, we know that $\mathbb E(e^{i(c_k-c_j)})=0$ for $k\neq j$, since $c_k-c_j$ will be uniformly distributed mod $2\pi$ so $e^{i(c_k-c_j)}$ is uniformly distributed on the unit circle, so its expectation is the center of the unit circle, which is $0$. Then $\mathbb E(e^{i(c_k-c_j)})=1$ when $k=j$, since then $e^{i(c_k-c_k)}=1$. Thus, the above have $N$ non-zero terms, all of which are $1$, so
$$\mathbb E(|z|^2)=N.$$
All of this algebra just expresses the more general fact that, if you take several independently and randomly distributed vectors, each distribution having the mean at the origin, then the expected value of $|z|^2$ for their sum equals the sum of the expected values of $|v|^2$ in each distribution.
Then, to get to your problem, one just multiplies the sum of sines by $a$, giving a new amplitude of $Na^2$. The value of $k$ does not affect the amplitude, so is irrelevant in the answer.
If one wants to do it directly via the harmonic addition theorem given as (for this case)
$$A^2=\sum_{i=1}^N\sum_{j=1}^Na^2\cos(c_i-c_j)$$
we can just use linearity of expectation:
$$\mathbb E(A^2)=\sum_{i=1}^N\sum_{j=1}^N\mathbb E(a^2\cos(c_i-c_j)).$$
Then, again, this sum vanishes for $i\neq j$, since $c_i-c_j$ is distributed uniformly mod $2\pi$ and $\cos(\theta)$ has a mean of $0$. For $i=j$, we of course have $\cos(c_i-c_j)=1$, so the above sum has $n$ non-zero summands, each of which is $a^2$. So, the sum is
$$\mathbb E(A^2)=Na^2.$$
One should note that this calculation of expectation is precisely what happens when we take the real part of the equation we derived for $|z|^2$ before - indeed, the harmonic sum theorem follows from the work in this post.
Best Answer
$$\mu\left[X(t)\right]=\mu\left[A\sin(\omega_0t+B)\right]=\dfrac{1}{2\pi}\int_{-\pi}^{\pi}A\sin(\omega_0t+B)dB=0\\R\left[X(t+\tau)X^{*}(t)\right]=\dfrac{1}{2\pi}\int_{-\pi}^{\pi}A^2\sin(\omega_0t+B)\sin\left(\omega_0(t+\tau)+B\right)dB\\=\dfrac{1}{4\pi}\int_{-\pi}^{\pi}A^2\cos(\omega_0\tau)dB-\dfrac{1}{4\pi}\int_{-\pi}^{\pi}A^2\cos(2\omega_0t+2B+\omega_0\tau)dB\\=\dfrac{A^2}{2}\cos(\omega_0\tau)$$