As this is very long, I'll put the final numerical answer up front:
$$C(n) = (n-1)\left(\frac{2}{n\cdot (n+1)!}\right)^{\frac1{n-1}}$$
This is, fundamentally, an optimization problem. We're not going to think of it as an inequality; we're going to be looking for a maximum or minimum.
And no, this has nothing to do with Hardy's inequality.
First, generalize. Those coefficients $\frac1k-\frac{2}{n(n+1)}$ and $k^2$ are a bit messy, and we don't want to deal with them yet. So we don't; we replace them with arbitrary positive sequences $u_k$ and $v_k$ respectively.
Next, the constraint $\prod_k a_k$ is inconvenient - so we kill it. Replace $\frac1{a_k}$ with $\frac{\prod_i a_i}{a_k}=\prod_{i\neq k} a_i$ in the sum on the right-hand side. Both sides of the new comparison
$$\sum_{k=1}^n u_k a_k \ge C(n)\left(\sum_{k=1}^n v_k\prod_{i\neq k} a_i\right)^{\frac1{n-1}}$$
are homogeneous of degree $1$, so the constraint no longer matters and we are free to let the $a_k$ range over all nonnegative (not all zero) values, or to set a new constraint somewhere more convenient.
What's convenient? I wobbled around somewhat in my work, but settled on the affine constraint $\sum_{k=1}^n u_ka_k = n$. Raise both sides of the inequality to the $n-1$ power, and we seek the largest $C(n)$ so that
$$\left(\frac{n}{C(n)}\right)^{n-1} \ge \sum_{k=1}^n v_k\prod_{i\neq k} a_i = F(\vec{a})$$
or equivalently the maximum of that right hand side, subject to the constraint $\sum_k u_ka_k = n$.
This is now a constrained optimization problem, to which we can apply familiar techniques. The region is a $n-1$-dimensional simplex in $\mathbb{R}^n$ with vertices on the coordinate axes. All of the facets up to dimension $n-3$ are in subspaces where at least two of the $a_k$ are zero - which then makes every term of $F(\vec{a})$ zero. We won't get a maximum there. The $n-2$-dimensional faces are each in a coordinate hyperplane $x_k=0$. There, all but one term of $F(\vec{a})$ is zero - but one term is enough to call for a closer look. A quick application of AM-GM finds the lone critical point where $u_ia_i=\frac{n}{n-1}$ for each $i\neq k$, and $F(\vec{a})=\left(\frac{n}{n-1}\right)^{n-1}\cdot\frac{v_k}{\prod_{i\neq k}u_i}$. For extremely skewed patterns of the $u_i$ and $v_i$, where some $u_kv_k$ is much larger than all the others (more than $n-1$ times their sum), this can be the global maximum - but it's not going to happen in the case we're worried about, at least for large $n$.
That leaves finding the critical point inside the simplex. For this, we go to Lagrange multipliers. The constraint equation is easy enough to handle, with $j$th partial derivative $u_j$. Then
$$\frac{\partial F}{\partial a_j} = \sum_{k\neq j}v_k\prod_{i\neq j,k}a_i = \frac1{a_j}\sum_{k\neq j}v_k\prod_{i\neq k}a_i =\frac1{a_j}\left(F(\vec{a})-v_j\prod_{i\neq j}a_i\right)$$
Moving the $\frac1{a_j}$ factor, the Lagrange multiplier system for the critical point is (repeating for all $j$)
$$\lambda u_ja_j = F(\vec{a})-v_j\prod_{i\neq j}a_i$$
Sum them all up, and
$$\lambda \sum_j u_ja_j = n\lambda = (n-1)F(\vec{a})$$
Now, go back to that $j$th equation. Multiply and divide by $a_j$ to complete the product, and substitute our new formula for $F(\vec{a})$:
$$u_ja_j -\frac{n\lambda}{n-1} +\frac{v_j}{a_j}\prod_{i=1}^n a_i = 0$$
Let that product of all the $a_i$ be $A$. Treat it as a constant (hey, at least it doesn't depend on $j$), and solve the quadratic for $a_j$:
$$a_j = \frac{\frac{n\lambda}{n-1} \pm\sqrt{\left(\frac{n\lambda}{n-1}\right)^2 - 4Au_jv_j}}{2\lambda u_j} = \frac{n\pm\sqrt{n^2 -4(n-1)^2Bu_jv_j}}{2(n-1)u_j}$$
where $B=\frac{A}{\lambda^2}$ is a new constant.
OK, now it's time to get specific. We have $u_j=\frac1j-\frac{2}{n(n+1)}$ and $v_j=j^2$, so $u_jv_j=j-\frac{2j^2}{n(n+1)}$ and the equation becomes
$$u_ja_j = \frac{n\pm \sqrt{n^2 - 4(n-1)^2jB + \frac{8(n-1)^2j^2B}{n(n+1)}}}{2(n-1)}$$
Which sign should we take on those square roots? Since our constraint is that the sum of the $u_ja_j$ is $n$, most of them should be reasonably close to $1$ - and that means taking plus signs. Assume for now that they're all plus signs; we can revisit later if it doesn't work out.
And now... a miracle happens. This was labeled as contest math. I have faith that there is a solution - so I'm going to guess that it's nice. What would be nice? Having those square roots all come out rational would be nice - and that calls for $B=\frac{2n}{(n-1)^2(n+1)}$ to complete the square. Then the quantity inside the square root is $\left(n-\frac{4j}{n+1}\right)^2$ and
$$u_ja_j = \frac{2n-\frac{4j}{n+1}}{2(n-1)}=1+\frac1{n-1}-\frac{2j}{n^2-1}$$
This gives $\sum_j u_ja_j = n$ as it should be - check.
What about the product? There's no constraint there; $B$ came from the combination of the product $A$ of the $a_i$ and a freely varying parameter $\lambda$, so it'll be consistent no matter what. We have indeed found the critical point with this guess.
Next, we need to find $F(\vec{a})$ at this critical point. Fortunately, the formula for $a_j$ is simpler than that for $u_ja_j$; it comes out to just $a_j = \frac{jn}{n-1}$. Then $v_j\prod_{i\neq j}a_j = j^2\cdot \frac{n!}{j}\cdot\left(\frac{n}{n-1}\right)^{n-1}$ and
$$F(\vec{a}) = \sum_{j=1}^n j^2\cdot \frac{n!}{j}\cdot\left(\frac{n}{n-1}\right)^{n-1} = \frac{n!\cdot n^{n-1}\cdot n(n+1)}{(n-1)^{n-1}\cdot 2}$$
Set this equal to $\left(\frac{n}{C(n)}\right)^{n-1}$ to find the answer we seek. That comes out to
$$C(n) = (n-1)\left(\frac{2}{n\cdot (n+1)!}\right)^{\frac1{n-1}}$$
One last loose end - could there be any other critical points? If we increased $B$, that would decrease all of the $a_j$, and even taking plus signs on all the square roots wouldn't be enough for the $u_ja_j$ to have a large enough sum. Nope. If we decreased $B$, we'd have to flip a square root to the negative for some $j$. In that case, note that our quadratic formula gives $u_ja_j\le \frac{n}{n-1}$ with equality only if $B=0$ and we take the plus sign. Set $B=0$ and we get a solution with some $a_j$ zero and each other $u_ja_j=\frac{n}{n-1}$. We've already dismissed these face critical points as maximum possibilities - the values of $F(\vec{a})$ increase as we move into the interior from there - but they've shown up anyway. That's it for critical points.
On that note, we do have to pay attention to the points on the faces for small $n$. For $n=2$, the condition of having one $u_kv_k$ at least $n-1$ times the rest is automatically reached, so we look closer - and with $u_1=\frac23$, $u_2=\frac16$, $v_1=1$, $v_2=4$, they're actually equal. This leads to the left and right sides of the original inequality both being constant multiples of the same sum, and the inequality becomes an equation. It doesn't matter where we sample $F$, since $F$ is constant on that segment. While the logic of the argument comes out a little odd, we have the right value there.
For $n=3$, $u_1v_1=\frac56\cdot 1$, $u_2v_2=\frac13\cdot 4$, and $u_3v_3=\frac16\cdot 9$. None of them dominates, and the points on the faces result in smaller values of $F$. From there on out, the advantage of the interior point becomes even stronger.