I had a look at the paper you linked to, here's the best I can do.
Eaves defines the inner product in the relevant space as
$$
\langle\mathbf{a}, \mathbf{b}\rangle = \frac{1}{V}\int_\Omega \mathbf{a}\cdot\mathbf{b} \enspace\text{d}V.
$$
In the problem you are looking at, the Lagrangian $\mathcal{L}$ is a functional of $\mathbf{u}$ involving the inner product, so it boils down to being able to take the functional derivative of the inner product with respect to one of its arguments.
The calculation is as follows (see wikipedia for a derivation of the functional derivative)
\begin{align}
\int_\Omega \left(\frac{\delta \langle\mathbf{a}, \mathbf{b}\rangle }{\delta \mathbf{b}}\right)_\mathbf{a} \enspace \text{d}V &= \frac{1}{V}\int_\Omega \left(\frac{\delta}{\delta \mathbf{b}}\right)_\mathbf{a}\mathbf{a}\cdot\mathbf{b} \enspace\text{d}V\\
&=\frac{1}{V}\int_\Omega \left(\frac{\partial}{\partial \mathbf{b}} - \nabla\cdot\left(\frac{\partial}{\partial (\nabla\mathbf{b})} \right) \right)_\mathbf{a}\mathbf{a}\cdot\mathbf{b} \enspace\text{d}V\\
& =\frac{1}{V} \int_\Omega \mathbf{a} \enspace \text{d}V
\end{align}
where I have used
$$
\frac{\partial}{\partial (\nabla\mathbf{b})}(\mathbf{a}\cdot\mathbf{b}) = \mathbf{0}
$$
and
$$
\left[\frac{\partial}{\partial \mathbf{b}}(\mathbf{a}\cdot\mathbf{b})\right]_i = \frac{\partial}{\partial b_i} a_j b_j = a_j\underbrace{\frac{\partial}{\partial b_i}}_{\delta_{ij}} = a_i=[\mathbf{a}]_i.
$$
Thus, omitting the normalising factor $1/V$,
$$
\left(\frac{\delta \langle\mathbf{a}, \mathbf{b}\rangle }{\delta \mathbf{b}}\right)_\mathbf{a} = \mathbf{a}
$$
Eaves also defines the adjoint operator as
$$
\langle\mathbf{a}, \mathcal{A}\mathbf{b}\rangle=\langle\mathcal{A}^\dagger\mathbf{a}, \mathbf{b}\rangle.
$$
This means that
$$
\left(\frac{\delta \langle\mathbf{a}, \mathcal{A}\mathbf{b}\rangle }{\delta \mathbf{b}}\right)_\mathbf{a} = \left(\frac{\delta \langle\mathcal{A}^\dagger\mathbf{a}, \mathbf{b}\rangle }{\delta \mathbf{b}}\right)_\mathbf{a}=\mathcal{A}^\dagger\mathbf{a}.
$$
The important part of the Lagrangian you are considering is
$$
\mathcal{L} = \Bigg[\mathbf{v}, \underbrace{\left(\frac{\partial}{\partial t}+\mathcal{A}-\frac{1}{\text{Re}}\nabla^2\right)}_\mathcal{B}\mathbf{u} -\nabla p\Bigg] = \int_0^T \langle\mathbf{v},\mathcal{B}\mathbf{u}-\nabla p\rangle\text{d}t.
$$
where I have defined a new operator $\mathcal{B}$. Upon variation with respect to $\mathbf{u}$, the term containing $p$ vanishes as it does not depend on $\mathbf{u}$, and the integral with respect to time is unimportant, so from above we are left with
$$
\frac{\delta \mathcal{L}}{\delta\mathbf{u}} = \mathcal{B}^\dagger \mathbf{v}.
$$
Finally,
$$
\mathcal{B}_{ij} = \left(\frac{\partial}{\partial t}+\mathcal{A}-\frac{1}{\text{Re}}\nabla^2\right)_{ij} = \delta_{ij}\left(\frac{\partial}{\partial t}-\frac{1}{\text{Re}}\nabla^2\right) + \mathcal{A}_{ij}
$$
so when the adjoint of $\mathcal{B}$ is taken, $\delta_{ij}$ is symmetric so remains the same, but $\mathcal{A}_{ij}$ is not, so its adjoint must be computed.
