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.
There is one famous example in which there is no solution for the Stokes' flow case: Stokes flow around a cylinder, which is approriately named Stokes' paradox. In this case it is impossible to match the boundary conditions both at infinity (uniform flow) and at the cylinder surface (no-slip) with Stokes flow dynamics. See e.g. paragraph 6.4 of the Fluid Dynamics Lecture Notes of Jacques Lewalle. The breakdown of a solution means that there will always be some inertial effect, regardless how small the Reynolds number is.
An approximate solution to this issue was first proposed by Oseen who introduced a linearized inertial term to account for inertial contributions in the far field.
Later, Proudman and Pearson calculated a more precise solution through asymptotic expansions and matching of the far field solution and the near-cylinder solution.
So to answer your question: existence for Stokes' flow is not guaranteed even though the criterion $Re << 1$ is satisfied. A pretty good explanation for the precise reasons of non-existence is given in chapter 7 of the Fluid Dynamics I lecture notes by Prof. Childress. In the same document they also show (in paragraph 7.2) that Stokes' flow does exhibit uniqueness for non-trivial cases (i.e. $\textbf{u}\neq0$).
For more details on conditions for solvability: there is quite a lot of mathematical fluid dynamics literature on the solvability of Stokes' flow. Nazarov and Pileckas reference a number of them in their paper with the telling title "On the Solvability of the Stokes and Navier-Stokes Problems in the Domains That Are Layer-Like at Infinity"
Michiel's answer is more from the aspects of physics. Here is the pde style answer.
Short answer: The Stokes flow's variational problem is well-posed (uniqueness and existence) in certain Hilbert spaces pair which relies on inf-sup condition.
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'. $$
Stokes problem is well-posed when the follpwing operator is an isomorphism: $$ \mathscr{S}:(\b{v},q)\in X\times Y \mapsto (A\b{v}+B'q,B\b{v})\in X'\times Y'. $$
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):
$B: \mathrm{ker}(B)^{\perp}\subset X \to Y'$ is an isomorphism and $\|\b{v}\|_X \leq \beta \|B\b{v}\|_{Y'}$.
$B': Y \to (\mathrm{ker}(B)^{\perp})'$ is an isomorphism and $\|q\|_Y \leq \beta \|B'q\|_{X'}$.
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).
While the answer by Shuhao Cao is fully correct, there is a simpler (though slightly less robust) way to look at things.
Consider the Stokes problem: \begin{alignat}{2} - \mu \mathbf{\Delta} \mathbf{u} + \nabla p ={}& \mathbf{f}&\quad&\textrm{in }\Omega, \label{stokes:eq:stokes.momentum}\\ \nabla\cdot \mathbf{u} ={}& 0&\quad&\textrm{in }\Omega, \label{stokes:eq:stokes.mass}\\ \mathbf{u} ={}& \mathbf{0}&\quad&\textrm{on }\partial\Omega, \label{stokes:eq:stokes.no.slip}\\ \int_{\Omega} p ={}& 0, &\quad&\label{stokes:eq:stokes.pressure} \end{alignat} for some source term $\mathbf{f}\in L^2(\Omega)^d$ and bounded Lipschitz domain $\Omega\subset\mathbb{R}^d$.
A typical weak formulation of this problem is as follows: find $\mathbf{u}\in H_0^1(\Omega)^d$, $p\in L_0^2(\Omega)$ such that \begin{alignat}{2} \mu\mathrm{a}(\mathbf{u}, \mathbf{v}) + \mathrm{b}(\mathbf{v}, p) ={}& \int_{\Omega}\mathbf{f}\cdot\mathbf{v}\quad&\forall& \mathbf{v}\in H_0^1(\Omega)^d, \\ -\mathrm{b}(\mathbf{u}, q) ={}& 0 \quad&\forall& q\in L_0^2(\Omega), \end{alignat} where $H_0^1(\Omega)^d$ denotes the space of functions with square-integrable gradient and zero trace on $\partial\Omega$, $L_0^2(\Omega)$ denotes the space of square integrable functions with zero mean value (i.e. $\int_\Omega q = 0$), and the bilinear forms are given by \begin{equation} \mathrm{a}(\mathbf{u}, \mathbf{v}) := \int_{\Omega}\nabla\mathbf{u} : \nabla\mathbf{v} \quad;\quad\mathrm{b}(\mathbf{v},q) := -\int_{\Omega}(\nabla\cdot\mathbf{v}) q. \end{equation} However, if one considers a stricter regularity of $p\in H^1(\Omega)$, the above formulation is fully equivalent to the following (harder to implement but easier to analyse) formulation: find $\mathbf{u}\in \mathbf{V}$, $p\in Q$ such that \begin{alignat}{2} \mu\mathrm{a}(\mathbf{u}, \mathbf{v}) ={}& \int_{\Omega}\mathbf{f}\cdot\mathbf{v}\quad&\forall& \mathbf{v}\in \mathbf{V}, \tag{1}\label{stokes:eq:stokes.weak.form.momentum}\\ \int_{\Omega}\nabla p \cdot \nabla q ={}& \int_{\Omega}\mathbf{f}\cdot\nabla q - \mu\mathrm{a}(\mathbf{u}, \nabla q) \quad&\forall& q\in Q,\tag{2}\label{stokes:eq:stokes.weak.form.mass} \end{alignat} where I define \begin{equation} \mathbf{V} = \{\mathbf{v}\in H^1(\Omega)^d : \textrm{div }\mathbf{v} = 0,\, \mathbf{v}=\mathbf{0}\textrm{ on }\partial\Omega\}. \end{equation} and \begin{equation} Q = \{q\in H^1(\Omega) : \int_{\Omega}q = 0\}. \end{equation} So, in the second formulation, the zero divergence constraint is enforced strongly. However, the pressure does not appear in equation \eqref{stokes:eq:stokes.weak.form.momentum}. Moreover, due to the zero BCs and a Poincar'e inequality, $\mathrm{a}(\cdot,\cdot)$ defines an inner-product on $\mathbf{V}$. Therefore, by the Riesz representation theorem there exists a unique velocity field $\mathbf{u}\in\mathbf{V}$ satisfying equation $\eqref{stokes:eq:stokes.weak.form.momentum}$. Similarly, by the zero mean value condition and a Poincar'e--Wirtinger inequality, the bilinear form $\int_{\Omega}\nabla p \cdot \nabla q$ defines an inner-product on $Q$. Therefore, by the Riesz representation theorem there exists a unique pressure field $p\in Q$ satisfying equation $\eqref{stokes:eq:stokes.weak.form.mass}$.
Of course, one can relax the extra regularity on the pressure by simply considering the velocity defined by equation \eqref{stokes:eq:stokes.weak.form.momentum} and defining the pressure by \begin{equation} \nabla p = \mathbf{f} + \mu \mathbf{\Delta} \mathbf{u} \end{equation} which requires the solution $\mathbf{u}$ to \eqref{stokes:eq:stokes.weak.form.momentum} to have (at least in a weak sense) a Laplacian. I am fairly certain that even for cracked domains one has $\mathbf{\Delta}\mathbf{u}\in L^{\frac32}(\Omega)^d$.
