I am looking for a reference or proof of a lemma (if it's true) or a counter-example otherwise. It goes as follows:
Let $B_1$ and $B_2$ are two concentric balls of radius $1$ and $2$ in some $n$-dimensional Euclidean space. Then for any $f$ and $F$, there exists a $C>0$ such that $$ \text{if }\;\Delta f=\text{div}F\; \text{ weakly, then }\; \|f\|_{L^p_1(B_1)}\leq C\big(\|F\|_{L^p(B_2)}+\|f\|_{L^p(B_2)}\big). $$
I would call this estimate "the classical Calderon-Zygmund estimate" but indeed it is hard to track down a statement for the right-side in divergence form. Usually it is stated as an estimate from $L^p$ to $W^{2,p}$, when what you want is $W^{-1,p}$ to $W^{1,p}$. This is unfortunate, because the latter is the more natural and useful/important estimate! These statements are of course related, and you can get from one to the other.
However, for a direct proof of the statement you want, the only place I know to find it is in our book: see Proposition 7.3 (which actually proves a slightly more general statement than the one you want). The arxiv version is here.
Note that the statement is of course only true for $1<p<\infty$, there are counterexamples given in our book as well.
Here is a sketch of the proof. By assumption $\int_{B_2} f \Delta \phi=-\int_{B_2} F\cdot D\phi$ for every $\phi \in C_c^\infty (B_2)$. Fix $\eta \in C_c^\infty(B_2)$, $\eta=1$ in $B_1$ and $\psi \in C_c^\infty (\mathbb R^n)$. Applying the equality above to $\phi=\eta \psi$ we get with $v=\eta f$ $$ \left |\int_{\mathbb R^n} v(\psi-\Delta \psi)\right | \leq C(\|F\|_p+\|f\|_p)\|\psi\|_{1,p'}. $$ Inserting the difference quotient $\psi=D_{-h}\phi:=|h|^{-1}(\psi(\cdot-h)-\psi(\cdot))$ one arrives at $$ \left |\int_{\mathbb R^n} D_h v(\phi-\Delta \phi)\right | \leq C(\|F\|_p+\|f\|_p)\|\phi\|_{2,p'}. $$ By density this inequality extends to $\phi \in W^{2,p'}(\mathbb R^n)$. Choosing $\phi \in W^{2,p'}(\mathbb R^n)$ such that $$\phi-\Delta \phi=D_hv|D_h v|^{p-2}, \quad \|\phi\|_{2,p'} \leq C\|D_hv\|_p^{p-1}$$ one obtains $\|D_hv\|_p \leq C(\|F\|_p+\|f\|_p)$ and the estimate follows letting $h \to 0$, since $f=v$ in $B_1$.
Suppose $p>1$. For each $i\in \{1,\ldots,n\}$, consider $\varphi_i\in W^{2,p}(B_2)\cap W^{1,p}_0(B_2)$ be the solution of $$ -\Delta \varphi_i = F_i \textrm{ in } B_2, $$ and let $\Phi=(\varphi_i)_{1\leq i \leq n}$.
Alternatively, you could set $\varphi_i = \Gamma\star F_i$, where $\Gamma$ is the fundamental solution of the Laplacian in free space ($C|x|^{n-2}$ for $n>2$, $C\ln |x|$ when $n=2$) and $F$ is extended by zero outside $B_2$ .
In both cases, it satisfies $\| \phi \|_{W^{2,p}(B_2)} \leq C\|F\|_{L^p(B_2)}$.
Now consider $f+\textrm{div} \Phi$ on $B_{3/2}$. It is harmonic, and therefore satisfies $$ \| f + \textrm{div} \Phi \|_{C^1(\overline{B_1})} \leq C \| f + \textrm{div} \Phi \|_{L^p(\overline{B_{3/2}})} \leq C \left( \| f \|_{L^p(B_{2})}+ \|F\|_{L^p(B_2)} \right). $$ Finally, $$\| f \|_{W^{1,p}(B_1)} \leq C \left(\| f + \textrm{div} \Phi \|_{C^{1}(\overline{B_1})} + \| \textrm{div} \Phi \|_{W^{1,p}(B_1)}\right)\leq C \left( \| f \|_{L^p(B_{2})}+ \|F\|_{L^p(B_2)} \right).$$
The unproven fact here is the first statement, regarding the Dirichlet problem for the Laplacian. This is however true, and can be done for example using the explicit form the Green function on the ball, I believe. That doesn't hold for $p=1$ though.
