The issue is that the equation
$$ \mathbf{u} \cdot \mathbf{n}=0$$
or
$$ \mathbf{u} \cdot \mathbf{t}=0$$
is a constraint equation rather than a normal flux, which is what enters naturally in the weak form of the equation.
The weak form derivation goes something like this. Multiply by $\mathbf{v}$:
$$-\int_\Omega \nabla \cdot \nabla \mathbf{u} \cdot \mathbf{v} dV=\int_\Omega \mathbf{f} \cdot \mathbf{v} dV$$
Now, perform partial integration:
$$-\int_{\partial\Omega} (\mathbf{n} \cdot \nabla \mathbf{u}) \cdot \mathbf{v}dS+\int_\Omega \nabla \mathbf{u} \cdot \nabla \mathbf{v} dV=\int_\Omega \mathbf{f} \cdot \mathbf{v} dV$$
The interesting term here is the boundary integral:
$$-\int_{\partial\Omega} (\mathbf{n} \cdot \nabla \mathbf{u}) \cdot \mathbf{v}dS$$
which shows that we control the flux/force $\mathbf{g}$ on the boundary on the following form:
$$-\int_{\partial\Omega} \mathbf{g} \cdot \mathbf{v}dS$$
On strong form this becomes:
$$\mathbf{n} \cdot \nabla \mathbf{u}=\mathbf{g}$$
So the weak form shows that, on the boundary, we do not get direct control over $\mathbf{u}$ and its normal and tangential projections. Instead we only get "access" to the normal projection of the gradient:
$$\mathbf{a}_n=\mathbf{n} \cdot \nabla \mathbf{u}$$
(Note that this vector can be "leaning" and doesn't have to be in the normal direction.)
One formal way to get control over constraints is to add a Lagrange multiplier to the boundary flux. Something like:
$$\mathbf{n} \cdot \nabla \mathbf{u}=\mathbf{P} \vec{\lambda}$$
where $\mathbf{P}$ is a matrix, and find a vector $\vec{\lambda}$ so that the constraint is fulfilled. This vector represents the flux that you need in order to fulfill the constraint (to every constraint there is a corresponding flux).
Now we are getting into theory of Lagrange multipliers.
Let's now partition the boundary up in two (possibly more as needed) parts. One $\partial \Omega_1$ with a Neumann condition and one $\partial \Omega_2$ with a constraint (general Dirichlet) condition.
We can now use the following energy functional (Dirichlet's principle):
$$E=\frac{1}{2}\int_\Omega \nabla \mathbf{u} \cdot \nabla \mathbf{u} dV-\int_{\Omega} \mathbf{f} \cdot \mathbf{u} dV-\int_{\partial \Omega_1} \mathbf{g} \cdot \mathbf{u} dS-\int_{\partial \Omega_2} (P\mathbf{u}-b) \cdot \mathbf{\lambda} dS$$
where $\lambda=(\lambda_1,\lambda_2)$ is a Lagrange multiplier vector.
To simplify matters, let's only consider part of the functional:
$$E=\frac{1}{2}\int_\Omega \nabla \mathbf{u} \cdot \nabla \mathbf{u} dV-\int_{\partial \Omega_2} (P\mathbf{u}-b) \cdot \mathbf{\lambda} dS$$
Now, consider the stationarity of the first variation of this functional:
$$\delta E=\int_\Omega \nabla \mathbf{u} \cdot \nabla \delta\mathbf{u} dV-\int_{\partial \Omega_2} P\delta \mathbf{u} \cdot \mathbf{\lambda} dS-\int_{\partial \Omega_2} (P\mathbf{u}-b) \cdot \delta \mathbf{\lambda} dS=0$$
If we set $\mathbf{v}=\delta \mathbf{u}=(\delta u_1,\delta u_2)=(v_1,v_2)$ and $\mu=\delta \lambda=(\mu_1,\mu_2)$ this becomes:
$$\int_\Omega \nabla \mathbf{u} \cdot \nabla \mathbf{v} dV-\int_{\partial \Omega_2} P\mathbf{v} \cdot \mathbf{\lambda} dS-\int_{\partial \Omega_2} (P\mathbf{u}-b) \cdot \mathbf{\mu} dS=0$$
If we now have constraints in both tangential and normal direction, assuming 2D and $\mathbf{b}=0$, then we have the constraint relationship:
$$P\mathbf{u}=\begin{bmatrix}
n_x & n_{y}\\
t_{x} & t_{y}\\
\end{bmatrix} \begin{bmatrix}
u_1\\
u_2\\
\end{bmatrix}= \begin{bmatrix}
0\\
0\\
\end{bmatrix}$$
Let's say we now only have a normal constraint, then we get:
$$P\mathbf{u}=\begin{bmatrix}
n_x & n_{y}\\
0 & 0\\
\end{bmatrix} \begin{bmatrix}
u_1\\
u_2\\
\end{bmatrix}= \begin{bmatrix}
n_x u_1+n_y u_2\\
0\\
\end{bmatrix}=\begin{bmatrix}
0\\
0\\
\end{bmatrix}
$$
Here, the constraint matrix is:
$$A=\begin{bmatrix}
n_x & n_{y}\\
0 & 0\\
\end{bmatrix}$$
The two Lagrange multiplier terms are now:
$$\int_{\partial \Omega_2} P\mathbf{v} \cdot \mathbf{\lambda} dS=\int_{\partial \Omega_2} (n_x v_1+n_y v_2) \lambda_1 dS$$
and
$$\int_{\partial \Omega_2} P\mathbf{u} \cdot \mathbf{\mu} dS=\int_{\partial \Omega_2} (n_x u_1+n_y u_2) \mu_1 dS$$
Note that whether you put the constraint in the first or second row of the constraint matrix is somewhat arbitrary.
For the equation that is not acted on by the Lagrange multiplier you can apply a Neumann condition.
I didn't find any quite satisfactory references. These are some that might be helpful:
https://users.oden.utexas.edu/~oden/Dr._Oden_Reprints/1982-008.finite_element.pdf
Also the appendix of this book:
https://www.amazon.com/Incompressible-Finite-Element-Isothermal-Laminar/dp/0471492507
These references are far from perfect. I will put in more here if I find any.