Let $D \subset \mathbb{R}^2$ be smooth bounded domain. Let $$G=\{\nabla f : f\in C^2(\overline{D})\}.$$ Is this space dense in $L^\infty(D)^2$ w.r.t the topology of $L^2$? I have the intuition that this is not true. Maybe, it could be true if we add some condition, but I'm looking for a large class so that the result holds.
Any hint or reference would be appreciated.
No, this is not true. Gradient fields have zero curl, but general vector fields do not have this property. You cannot even approximate all $C^1$-vector fields.
Take $g\in C^1(\bar D)^2$. Take $f_n \in C^2(\bar D)$ such that $\nabla f_n \to g$ in $L^2(D)$. Let $v\in C_c^\infty(D)$. Then $$ 0=\int_D \nabla \times \nabla f_n(x) \cdot v(x)dx = -\int_D \nabla f_n(x) \cdot (\nabla \times v(x)) dx \\ \to -\int_D \int_D g(x) \cdot (\nabla \times v(x)) dx = \int_D \nabla \times g \cdot v(x) dx. $$ Hence $\nabla \times g=0$ necessarily. If $g$ would be just in $L^2(D)^2$, the same reasoning shows that the curl of $g$ (in the weak sense) is zero.
In general, one has the following result for Lipschitz domain $D$: $$ L^2(D)^n = \{ \nabla p: \ p\in H^1_0(D)\} \otimes^\perp \{u\in L^2(D)^n : \ div(u)=0, \ u\cdot n=0 \text{ on } \partial D\}. $$ The conditions on $u$ in the second set have to be formulated in a weak sense. For a reference, see Temam: Navier Stokes equations, or any other book on Navier Stokes equations. The first set can be further decomposed. Iirc, elements in these subspaces can be approximated by smooth functions from these spaces.
Hence, gradients of smooth vector fields are always orthogonal to divergence-free vector fields.