What are the solution existence and uniqueness conditions for Stokes' flow?
$$\begin{gathered} \nabla p = \mu \Delta \vec{u} + \vec{f} \\ \nabla \cdot \vec{u} = 0 \end{gathered}$$
Maybe you could also provide some articles or books about the topic? Most physics books seem not to care about these details.
Michiel's answer is more from the aspects of physics. Here is the pde style answer.
Functional equations:$\newcommand{\b}{\boldsymbol}$ Suppose $\mu=1$, the Stokesian flow can be written in the following way (I believe this is called the Pressure-Velocity formulation): $$\left\{ \begin{aligned}&-\Delta \b{u} + \nabla p =\b{f}, \\ &\nabla\cdot \b{u} =0.\end{aligned}\right.\tag{1} $$ In operator form this can be written as the following abstract problem: $$ \begin{pmatrix}A & B' \\ B& 0\end{pmatrix}\begin{pmatrix}\b{u}\\p \end{pmatrix} = \begin{pmatrix}\b{f}\\0 \end{pmatrix}, $$ where $A = -\Delta$ is the vector Laplacian, and $B = -\nabla\cdot$ with $B' = \nabla$, $\b{u}\in X$, and $p\in Y$, the operators: $$ A: X\to X',\quad B:X\to Y', \quad \text{and}\quad B': Y\to X'. $$
Normally, the isomorphism is either proved using Lax-Milgram through coercivity to pin down a fixed point, or using Fredholm alternative.
The sufficient condition for this is a weak version of coercivity (you can view it as invertibility of an operator):
Then by closed range theorem, Babuska proved an equivalence of these conditions with the inf-sup condition(in that pdf link 1.1). Whenever that condition holds for certain Hilbert spaces pair $X\times Y$, (1)'s variational problem has a unique solution.
Weak formulation:
The weak formulation for the abstract version of (1) is then:
$$\left\{ \begin{aligned}\langle A \b{u},\b{v}\rangle + \langle p,B\b{v}\rangle =\langle\b{f},\b{v}\rangle,\quad \forall \b{v} \in X&, \\ \langle q,B\b{u}\rangle =0,\;\;\;\quad\quad \forall q \in Y.&\end{aligned}\right.\tag{2} $$ Possible pairs of Hilbert spaces $X\times Y$ mentioned above are: $$\begin{gathered} H^1_0(\Omega)\times \{q\in L^2(\Omega):\int_{\Omega}q=0\}, \\ H(\mathrm{div})= \{\b{v}\in L^2(\Omega):\nabla\cdot \b{v}\in L^2(\Omega)\}\times L^2(\Omega). \end{gathered}$$ Using integration by parts for (2) leads to: $$\left\{ \begin{aligned}\int_{\Omega} \mathrm{tr}\big((\nabla \b{u})^T \nabla \b{v}\big) +\int_{\Omega}p(\nabla \cdot \b{v}) =\int_{\Omega}\b{f}\cdot\b{v},\quad \forall \b{v} \in X&, \\ \int_{\Omega}q(\nabla \cdot \b{u}) =0,\quad\quad\quad \forall q \in Y.&\end{aligned}\right.\tag{3} $$
Problem (3) can be viewed as a constraint minimization problem for the following conjugate functionals also (viewing the pressure $p$ as a Lagrange multiplier): denote $E(\b{v}) = (\nabla \b{v}^T +\nabla \b{v})/2$ (symmetric part of the Jacobian), the stationary strain tensor, then $\mathrm{tr}\big((\nabla \b{v})^T \nabla \b{v}\big)= |E(\b{v})|^2 $ (a notation usually used in elasticity PDEs). Let $$ \mathcal{L}(\b{v},q) = \int_{\Omega}|E(\b{v})|^2 - \int_{\Omega} \b{f}\cdot \b{v} - \int_{\Omega} q\nabla\cdot \b{v}, $$ and $$\mathcal{J}(\b{v}) = \sup_{q\in Y}\mathcal{L}(\b{v},q) ,$$ then our goal is to minimize $\mathcal{J}$ in $X$ (like looking for a saddle point).