Let $\Omega\subset\mathbb{R}^d$, $d>1$ be an open and bounded Lipschitz domain. It is well-known, that for a given function $f\in L^q(\Omega)$ with $1<q<\infty$ and and such that $\int_\Omega f=0$, there exists a solution $u\in W^{1,q}_0(\Omega;\mathbb{R}^d)$ of the underdetermined PDE: $$\mathrm{div} u=f$$ and the estimate $$\lVert \nabla u \rVert_q\leq C\lVert f\rVert_q$$ holds for a scaling and translation invariant constant $C=C(\Omega,q)>0$. The solution to this PDE is $u:=\mathcal{B}f$, where $\mathcal{B}:L^q(\Omega)\to W^{1,q}(\Omega;\mathbb{R}^d)$ is so-called Bogovskii's operator.
I'm interested if there is a version of this theorem, preferably with the similar estimate, where the divergence operator is replaced by an arbitrary first order PDE operator $\mathcal{A}:=a^i\partial_i$ with some restrictions on coefficients $a^i$. Unfortunately, it seems that the original proof of Bogovskii does not carry over to the general case.