To begn with, your Boundary-Value Problem (BVP) is under-determined, because it lacks one boundary condition: because the PDE is elliptic and fourth-order, you need two boundary conditions, not only one. Because you insist on working with $H^2\cap H^1_0$, it seems that the hidden BC is
$$u=0\qquad\hbox{on }\partial\Omega.$$
So let us assume that your BVP includes this condition. Then your calculation works whenever $v\in C^2(\overline\Omega)$ (you should avoid $C^2_0$, because its elements satisfy $\partial_\nu v=0$ as well).If $v$ vanishes on the boundary (but not necessarily its normal derivative), a solution $u$ of the BVP does satisfy
$$B(u,v)=0.$$
If $\Omega$ is bounded with a smooth boundary, then
$$|\partial_\nu w|_{L^2(\partial\Omega)}\le C|w|_{H^2}\le C'|\Delta u|_{L^2},$$
where the second inequality is a kind of Poincar\'e inequality. Therefore $B$ is continuous over $H^2(\Omega)$ and $B(u,v)=0$ extends to all elements $v\in H^2\cap H^1_0$ by density.
The existence follows from Lax--Milgram, using the Poincar\'e inequality to show that $B$ is coercive over $H^2\cap H^1_0$.
Conversely, suppose that $u\in H^2\cap H^1_0$ and $B(u,v)=0$ for every $v\in H^2\cap H^1_0$. Then making the calculations backward, you find
$$\langle \Delta^2 u,v\rangle_\Omega+\int_{\partial\Omega}(\Delta u+\rho\partial_\nu u)\partial_\nu vds(x)=0,$$
where the first term is duality between $H^{-2}(\Omega)$ and its dual. Taking all $v\in{\mathcal D}(\Omega)$, you first obtain $\Delta^2 u=0$ in the sense of distributions. There remains therefore
$$\int_{\partial\Omega}(\Delta u+\rho\partial_\nu u)\partial_\nu vds(x)=0,$$
for every $v\in H^2\cap H^1_0$. Because $v\mapsto\partial_\nu v$ is onto over $H^{1/2}(\partial\Omega)$, this gives you $\Delta u+\rho\partial_\nu u=0$. Finally, $u=0$ on the boundary is ensured because $u\in H^1_0$.
You are using the wrong space, that's all. The correct space is
$$
X=\{u \in H^1_{\textrm{loc}}(\mathbb R^2): u(\cdot+n)=u(\cdot) \quad \forall n=(n_1,n_2)\in \mathbb Z^2\}.
$$
Then you apply Lax-Milgram and you are good to go. Look at
Asymptotic Analysis for Periodic Structures, by Bensoussan, Lions, Papanicolaou (1979) chapter 1, for example.
Answering your comment. $f$ is likewise extended periodically, namely $f$ is extended by zero into $f_0$, and then identified with $\sum_{n\in \mathbb Z^2} f_0(\cdot + n)$ which is periodic on $\mathbb R^2$ and $L^2_{\textrm{loc}} (\mathbb R^2)$.
Now, on say $(-2,2)^2$, $-\Delta u + u=f$ as a weak solution, and therefore interior regularity shows $u\in H^2((-\frac32,\frac32)^2)$. All in all, $u\in H^2_{\textrm{loc}} (\mathbb R^2)$. Consequently, $\partial_x u \in X$, and $\partial_y u \in X$.
Incidentally, because all coefficients are constant, you can write write $u$ explicitely as a (double) Fourier series, its coefficients being that of $f$, divided by $n_1^2+n_2^2+1$. It saves you the trouble of using Lax-Milgram.
Best Answer
The first observation is that the $u$ above satisfies $\nabla u=0$ on $\partial \Omega$ if and only if $f$ is orthogonal to all harmonic functions $v$ in $\Omega$, continuous up the the boundary. In fact, $\int_{\Omega} fv=\int_{\Omega} (\Delta u) v=\int_{\Omega} u \Delta v=0$, by the boundary conditions. Conversely, if this holds for $f$, then, since $u=0$ on $\partial \Omega$, $0=\int_{\Omega} (\Delta u) v=\int_{\partial \Omega} v \frac{\partial u}{\partial n}$ for every harmonic $v$. This gives $\nabla u=0$ at the boundary, since $v$ is arbitrary on $\partial \Omega$ and the tangential derivatives of $u$ are zero by the boundary conditions.
The second remark is that such $u$ vanishes in a neighborhood of the boundary where $f$ is zero (hence $u$ is harmonic). In fact, the equation and the boundary conditions imply that all the derivatives up to the second order vanish on the boundary. Continuing $u$ by zero across the boundary we obtain an harmonic function vanishing in an open set and hence in any connected set containing the boundary where $f=0$.
This observation allows to change $\Omega$ to a ball $B_R$ containing it. In fact, if $u$ solves the problem above in $\Omega$, since it is zero in a neighborhood of the boundary, it solves the same probelm in $B_R$ and conversely if it solves in $B_R$ then it is 0 in $B_R \setminus \Omega$, by the above argument and $u, \nabla u=0$ at $\partial \Omega$.
Let us take therefore $\Omega=B_R$. By the above discussion and by density, the problem has a solution iff $f$ is orthogonal to all harmonic polynomials $P$, that is $\int_{\mathbb R^n} P(x)f(x)=0$ or $\left (P(iD)\hat f\right )(0)=0$.