Is there an analytical solution for fluid flow in a square duct?

The equation is correct--- the (laminar) flow at small Reynolds number is given by making the flow be along the pipe, and substituting into the Navier Stokes equations, which reduces to your thing. The one issue is the sign--- $\Delta P$ is negative if you mean the flow is going to be in the positive z direction. I will absorb the constants and consider the problem on the box [-1,1]x[-1,1].

There is no analytic solution in elementary functions for this, because the problem is equivalent to solving Laplace's equation with certain Dirichlet boundary conditions on the square. But there is a simple and rapidly convergent series which gives you the answer.

The equation is

$$ \nabla^2 \phi = -A $$

Where A is the (negative) pressure gradient over the viscosity, in length units where the size of the box is 2. To solve this, first note that the quadratic function

$$ \phi_0(x,y) = {A\over 4} (2 - x^2 - y^2)$$

works, but doesn't satisfy the boundary conditions. This flow (plus a constant) gives the parabolic cylinder laminar flow profile, it satisfies no-slip on the circle of radius $\sqrt{2}$, just touching the corners of the square. You replace $\phi$ by $\phi_0 - {A\over 4} \phi$, and the new $\phi$ satisfies Laplace's equation:

$$ \nabla^2 \phi = 0$$

with the boundary conditions $$\phi(x,1) = 1-x^2 $$, so as to subtract out the nonzero velocity on the boundary square. This is the Dirichlet problem.

So you need to solve the Dirichlet problem on the square. In principle, the interior of the unit square can be conformally mapped into the circle, but the tranformation is ugly. So it is best to give a direct approximation.

Write $\phi$ as the real part of an analytic function $f(z)$, where $z=x+iy$. The symmetry of the problem tells you that the real part of $f(iz)$ is the real part of $f(z)$, so that (by analyticity) $f(iz)=f(z)$ and the analytic function f is an expansion in powers of $z^4$.

$$ \phi = \mathrm{Re} f(z) $$ $$ f(z) = a_0 + a_1 z^4 + a_2 z^8 + a_3 z^{12} $$

Then you know that on the boundary, $f(1+iy) = 1 - y^2 $, and this fixes the coefficients. To lowest nontrivial order, you keep only the constant and the $z^4$ term, and you find

$$ \mathrm{Re} f(1+iy) = a_0 + a_1 (1 - 6 y^2 + y^4) = 1-y^2 $$

Which gives (just by setting the lowest order terms equal) $a_1 = 1/6$ and $a_0 = 5/6$. The flow is then, to quartic order:

$$ V_z(x,y) = {A\over 4}( 7/6 - x^2 - y^2 - {1\over 6}(x^4 - 6 x^2 y^2 + y^4 )) $$

This is not a great approximation, but you can go to order 8, order 12, or order 16 and do the same thing to get polynomial approximations to the flow of any order.

I should add that there is a slowly converging solution by expansion in box modes, and it obeys the boundary conditions but it is inferior to the analytic method above--- the Fourier series of a constant function only falls of as 1/n.

Yes, there is an exact analytical solution to this Poisson partial differential equation.The well known problem is the laminar steady-state flow in rectangular ducts. The solution is given in the following references : Shah RK and London (1978) Laminar forced convection in ducts, Academic Press. Timoshenko and Goodier (1970) Theory of elasticity Mc Graw-Hill