An implicit a posteriori finite element error estimation method is presented to inexpensively calculate lower and upper bounds for a linear functional output of the numerical solutions to the three-dimensional Navier-Stokes (N-S) equations. The novelty of this research is to utilize an augmented Lagrangian based on a coarse mesh linearization of the N-S equations and the finite element tearing and interconnecting (FETI) procedure. The latter approach extends the a posteriori bound method to the three-dimensional Crouzeix-Raviart space for N-S problems. The computational advantage of the bound procedure is that a single coupled non-symmetric large problem can be decomposed into several uncoupled symmetric small problems. A simple model problem, which is selected here to illustrate the procedure, is to find the lower and upper bounds of the average velocity of a pressure driven, incompressible, steady Newtonian fluid flow moving at low Reynolds numbers through an endless square channel which has an array of rectangular obstacles. Numerical results show that the bounds for this output are rigorous, i.e. always in the asymptotic certainty regime, that they are sharp and that the required computational resources decrease significantly. Parallel implementation on a Beowulf cluster is also reported.