Regularizing a solenoidal vector field $u\in L^p(\Omega)^N$. Let $\Omega\subset\mathbb{R}^N$ be a bounded smooth domain and suppose that $u\in L^p(\Omega)^N$, $p\in (1,\infty)$. Assume that in the sense of distributions, $\operatorname{div}u=0$ where $\operatorname{div}u=\sum_{i=1}^N\frac{\partial u_i}{\partial x_i}$. 
Given $\epsilon>0$, can I find $v\in C^1(\overline{\Omega})$ such that $\|v-u\|_p<\epsilon$ and $\operatorname{div}v=0$?
Thank you
 A: The paper Extension and Representation of Divergence-free Vector Fields
on Bounded Domains may be helpful here.
I also add a couple of remarks.
Remark 1. In the Banach space of $L^p$-vector fields, the subspaces of conservative fields and of divergence-free fields are closed. Moreover, in the reflexive range the space of divergence-free $L^p$ fields is precisely the annihilator of conservative fields in $L^q$, the dual space.
Remark 2. For $1<p<\infty$, for smooth domains, the above subspaces are complemented in $L^p$: that is, they admit bounded projections. This comes from the boundedness of singular integrals on $L^p$. The key term is $L^p$ Hodge decomposition.  
Having a bounded projection allows you to mollify a divergence-free field and then project it back onto the subspace, without  going too far from the original field. 
As a starting point, look at this paper and references given there.

From an older answer: 
Remark 3, on reflecting a divergence-free field in the hyperplane $x_0=0$. Given $x=(x_1,\dots,x_n)$, denote $x^*=(x_1,x_2,\dots,-x_n)$. Here is how the reflection of divergence-free vector field $u$ should be defined:


*

*Write $u$ in components: $u=(u_1,\dots,u_n)$. 

*When $j<n$, use odd reflection $u_j(x^*)=-u_j(x)$. 

*For $j=n$ use even reflection: $u_n(x^*)=u_n(x)$. 


I claim that $\int u\cdot \nabla \varphi=0$ for any test function $\varphi$. Indeed, any $\varphi$ can be decomposed into even part 
$\frac12(\varphi(x)+\varphi(x^*))$ and odd part $\frac12(\varphi(x)-\varphi(x^*))$. 


*

*If $\varphi$ is even, then $u\cdot \nabla \varphi$ is odd by construction, and therefore integrates to zero. 

*If $\varphi$ is odd, then $u\cdot \nabla \varphi$ is even, therefore it suffices to integrate it over one half-space. And since $\varphi$ is odd, it vanishes on the boundary. By the divergence-free property  of the original field $u$, the integral of  $u\cdot \nabla \varphi$ is zero.

