The flux is given by
$$\iint_{\text{square}} dA \, \mathbf{F} \cdot \hat{\mathbf{n}} $$
where $\hat{\mathbf{n}}$ is the unit normal to the square, which in this case is simply the $z$ direction. Thus, $\mathbf{F} \cdot \hat{\mathbf{n}}$ is the $z$ component of $\mathbf{F}$, or $x y z$. Thus, the flux is
$$a \int_{-a/2}^{a/2} dx \, x \, \int_{-a/2}^{a/2} dy \, y = 0$$
The situation seems basically similar to a "dry clutch" in engineering. The simplest model of that, I think, would be almost identical to a block sliding on a flat surface, subject to frictional forces. And the common assumption there is that the kinetic friction force is proportional to the normal force applied to the block (which might be its weight), and to a coefficient of friction, but it is independent of the speed. This is the angular equivalent of that. The frictional torque between the two discs will act on the relative motion, and will conserve total angular momentum.
So we can start by writing
$$
I_1 \dot{\omega}_1(t) = -\tau , \qquad I_2 \dot{\omega}_2(t) = +\tau
$$
where
the dot represents the time derivative and
$\tau$, the frictional torque, is a constant.
We need to integrate these equations in time, from $t=0$, up to the point where $\omega_1(t)=\omega_2(t)=\omega_f$. Beyond that point, the discs will rotate together.
This is simple enough to do analytically. The angular velocities will change linearly in time
$$
\omega_1(t) = \omega_1(0) -\frac{\tau t}{I_1 } , \qquad
\omega_2(t) = \omega_2(0) +\frac{\tau t}{I_2 }
$$
An easy way to get a solution is to write an equation for the relative angular velocity $\omega=\omega_1-\omega_2$,
set $\omega(t)=0$, and solve it for $\tau t$:
\begin{align*}
\omega(t) &= \omega(0)
- \tau t \left(\frac{I_1+I_2}{I_1I_2}\right) = 0
\\
\quad\Rightarrow\quad
\tau t_f &= \left(\frac{I_1I_2}{I_1+I_2}\right) \omega(0)
= \left(\frac{I_1I_2}{I_1+I_2}\right)
[ \omega_1(0)-\omega_2(0)]
\end{align*}
where I've called this final time $t_f$.
We can check that this gives the correct answers by substituting back in
\begin{align*}
\omega_1(t_f) &= \omega_1(0) -\frac{1}{I_1 }\left(\frac{I_1I_2}{I_1+I_2}\right)
[ \omega_1(0)-\omega_2(0)]
=
\frac{I_1\omega_1(0) + I_2\omega_2(0)}{I_1+I_2}
\\
\omega_2(t_f) &= \omega_2(0) +\frac{1}{I_2 } \left(\frac{I_1I_2}{I_1+I_2}\right)
[ \omega_1(0)-\omega_2(0)]
=
\frac{I_1\omega_1(0) + I_2\omega_2(0)}{I_1+I_2}
\end{align*}
So, both equal to $\omega_f$ at that time $t=t_f$.
Of course,
if you want to make the frictional torque depend on
relative angular velocity in some complicated way,
the solution may require a computer.
But the underlying equations will be similar to the above.
[Edit following OP comments]
Note that I wrote the equation for $t_f$, time needed to reach equal angular velocities, as an expression for the product $\tau\, t_f$:
this product is equal to a function of the initial angular velocities and the moments of inertia.
You need to know $\tau$ before you can calculate $t_f$.
To evaluate $t_f$ for a particular physical case, you need to
multiply my equation for $\tau t_f$ by $1/\tau$ on both sides, i.e.
take $\tau$ over to the right hand side.
This illustrates that $t_f$ is inversely proportional to $\tau$,
if the other parameters are kept constant.
If the friction between the two discs is zero, it will take
an infinite time to reach the same angular velocity,
because the two discs will have no effect on each other.
If the friction is very large,
the time taken will be very short.
In any case, the integration of those equations
should stop at $t=t_f$,
since after that time the two discs are not rotating relative to each other.
So, to solve the problem you want to solve, you must provide the physics of the interaction between the discs. Let me emphasise that my solution is only based on the simplest assumption about this. Your situation might be more complicated. However, some aspects of my solution will still apply, e.g. the torque between the discs must still conserve total angular momentum. For more discussion of friction in general, see for example https://physics.stackexchange.com/questions/2408/does-the-force-of-kinetic-friction-increase-with-the-relative-speed-of-the-objec and https://physics.stackexchange.com/questions/154443/why-is-the-equation-for-friction-so-simple .
Best Answer
Let $S_e$ denote the cross section at the exit of the tube. The exact momentum flux is
$$Q_{\text{momentum}} = \int_{S_e}\rho \mathbf{V} \mathbf{V} \cdot \mathbf{n} \, dA$$
In your solution for momentum flux as the product of the mass flow rate and the bulk-averaged x-velocity at the exit, you are implying that
$$\int_{S_e}\rho \mathbf{V} \mathbf{V} \cdot \mathbf{n} \, dA = \frac{ \int_{S_e}\rho \mathbf{V} \cdot \mathbf{n} \, dA \int_{S_e}\mathbf{V} \, dA}{ \int_{S_e}\, dA}$$
This need not be true as an equality but may serve as a reasonable approximation, for example, in turbulent flow where the velocity profile is relatively flat.
To examine this in the case of laminar flow, assume that cross section is a disk of radius $R$ normal to the x-direction and the velocity profile is that of steady axisymmetric Poiseuille flow,
$$\mathbf{V} = C\left[1- \left(\frac{r}{R}\right)^2\right]\mathbf{e_x}$$
In this case we have $\mathbf{n} = \mathbf{e_x}$ and
$$\int_{S_e}\rho \mathbf{V} \mathbf{V} \cdot \mathbf{n} \, dA = 2\pi \rho C^2\int_0^R\left[1- \left(\frac{r}{R}\right)^2\right]^2 r\, dr \,\,\mathbf{e_x} = \frac{\pi \rho C^2R^2}{3}\, \mathbf{e_x} \\ \int_{S_e}\, dA = \pi R^2 \\ \int_{S_e}\mathbf{V} \, dA = 2\pi C \int _0^R\left[1- \left(\frac{r}{R}\right)^2\right] r\, dr\, \mathbf{e_x} = \frac{\pi C R^2}{2}\, \mathbf{e_x} \\ \int_{S_e}\rho \mathbf{V} \cdot \mathbf{n} \, dA = 2\pi \rho C \int_0^R\left[1- \left(\frac{r}{R}\right)^2\right] r\, dr = \frac{\pi \rho CR^2}{2}$$
Hence,
$$ \frac{ \int_{S_e}\rho \mathbf{V} \cdot \mathbf{n} \, dA \int_{S_e}\mathbf{V} \, dA}{ \int_{S_e}\, dA} = \frac{\pi\rho C^2 R^2}{4} \mathbf{e_x} \neq \frac{\pi\rho C^2 R^2}{3 } \mathbf{e_x} = \int_{S_e}\rho \mathbf{V} \mathbf{V} \cdot \mathbf{n} \, dA $$
Here your approximation is not particularly accurate -- deviating from the exact momentum flux by a factor of $3/4$.