Axial forces on a solenoid windings I understand that the windings in a solenoid experience a Lorentz force $\mathbf{f} = \mathbf{J} \times \mathbf{B}$, which tend to cause an outward pressure where $\mathbf{B}$ is directed along the solenoid axis (near the middle), and axial compression, where $\mathbf{B}$ has a radial component (near the ends).
In the case of a solenoid that is much longer than it's radius, the outward pressure is easily derived from energy considerations, but the same logic also seems to imply a force that tends to drive the solenoid apart axially.
Radial force
These forces can be derived from energy considerations. Consider an very long solenoid (or, at least, one where we neglect end effects) with radius $r$ and current $I$. Then 
$$\mathbf{B} = \frac{\mu_0 N I}{\ell}\hat{\mathbf{z}}$$
so the energy density inside is
$$u_B = \frac{1}{2} \frac{B^2}{\mu_0} = \frac{1}{2} \frac{\mu_0 N^2 I^2}{\ell^2}$$
and energy is
$$U_B = \frac{1}{2} \frac{\mu_0N^2 I^2}{\ell^2} \pi r^2 \ell $$
from which we derive a radial (generalized) force
$$f_r = -\frac{\partial U_B}{\partial r} = \frac{\mu_0 N^2 I^2}{2\ell^2} 2\pi r \ell $$
Note: The above expression is missing a minus sign. See Art Brown's answer below.

or, equivalently, a pressure 
$$P = \frac{\mu_0N^2 I^2}{2\ell^2} = \frac{B^2}{2\mu_0} $$
Axial force
It seems that, far away from the solenoid ends, there should also be an expansive force, which we can calculate. Consider a small segment of the solenoid with number of turns $\Delta N$ and length $\Delta \ell$:
$$f_z = -\frac{\partial U_B}{\partial (\Delta\ell)} = \frac{\mu_0\Delta N^2 I^2}{2\Delta\ell^2} \pi r^2  $$
It seems counter-intuitive, both from the point of view of the axial compression due to the ends, and thinking of a solenoid as a linear arrangement of attractive current loops. Where am I going wrong?
 A: A couple problems with your development:
1) (minor, now corrected) You lose a couple factors of $\mu_0$ when you substitute for $B$ in the expression for $u_B$.
2) (major)  You lose a minus sign in calculating the radial force; by your expression, the force should incorrectly tend to compress the coil, not expand it. 
What's going on?
When mechanical work is done, the change in total energy $ \Delta U_B$ includes not just the change in potential energy in the field $\Delta U_{pot}$ but also the electrical energy $\Delta U_{elec}$ required to maintain the solenoid current $I$ constant.
$$\Delta U_{B} = \Delta U_{elec} + \Delta U_{pot}$$
For the potential energy, think of the solenoid as a stack of $N$ dipoles, each of moment $m=I \pi r^2$.  The potential of each in the external field $B=\mu_0 NI/l$ produced by the rest of the dipoles is $-mB$ (since each dipole is aligned with $\boldsymbol{B}$).
If the circuit is energized by "turning on" each dipole sequentially, the first results in no change in the potential energy since the external applied field is zero at the time, while the last sees the "full-up" field B.  Integrating over all the moments, the total dipole potential energy is:
$$ U_{pot}=-\frac{1}{2} NmB = -\frac{1}{2} N I \pi r^2 \mu_0 \frac{NI}{l} = - \frac{B^2}{2 \mu_0} \pi r^2 l$$
which matches your expression for the total energy $U_B$ (once you fix the $\mu_0$'s), except for the minus sign.  
Now one gets the correct expansion force when calculating the radial force $f_r= - \frac{dU_{pot}}{dr}$; work is done on the wire.
Why the sign change? We can calculate these various incremental energy changes (mechanical, electrical, and total) for a radial displacement $\Delta r$, holding current $I$ (and therefore field $B$) constant:
$$ \Delta U_{mech} = - \Delta U_{pot} = \frac{B^2}{2 \mu_0} 2 \pi r l \Delta r $$
$$ \Delta U_{elec} = \int VI dt = I \int N \frac{d \Phi}{dt} dt = NI \Delta \Phi = NIB \Delta A= NI B 2 \pi r \Delta r = 2 \frac{B^2}{2 \mu_0} 2 \pi r l \Delta r $$
$$\Delta U_B = \frac{B^2}{2 \mu_o} 2 \pi r l \Delta r $$
showing that:
$$ \Delta U_B = \Delta U_{elec} + \Delta U_{pot} \text{ , and } \Delta U_B = - \Delta U_{pot}$$
The additional electrical energy is twice as large and opposite in sign to the change in potential energy, so the change in total energy is just the opposite of the change in potential energy, or the mechanical work done.
So, when calculating the force from the total energy, one must use the formula $f= + \frac{dU_B}{dr}$.
Finally, the force on an "end" dipole of the stack is just the standard:
$$ \boldsymbol{f = - \nabla (-m \cdot B)} = \boldsymbol{\nabla (m \cdot B)} $$
which pulls the dipole towards higher values of $B$, i.e. a tension on the coil in the axial direction. 

Update:  More generally, it is useful to rewrite the Lorentz force law, for magnetostatic situations, in terms of B only as:
$$ \boldsymbol{j \times B} = \frac{1}{\mu_0} \boldsymbol{ (\nabla \times B) \times B } 
= - \boldsymbol{\nabla} \left( \frac{B^2}{2 \mu_0} \right) + \frac{1}{\mu_0} \left( \boldsymbol{B \cdot \nabla}\right) \boldsymbol{B}$$
The first term corresponds to a magnetic hydrostatic pressure $B^2/(2 \mu_0)$, with the second acting as an additional tension along the field lines.  This expression can be integrated over a volume $V$ to give the total force $\boldsymbol{F}$ on that volume, or converted to a surface integral of the magnetostatic stress tensor:  
$$ \boldsymbol{F} = \int_{\partial V} \boldsymbol{ T \, dS } $$
where
$$ \boldsymbol{T} = \frac{1}{\mu_0} \left( 
\begin{array} {ccc} B_x^2 - B^2/2 & B_x B_y & B_x B_z \\
B_x B_y & B_y^2 - B^2/2 & B_y B_z \\ 
B_x B_z & B_y B_z & B_z^2 - B^2 /2 \end{array} \right)  $$
or 
$$ T_{ij} =  \frac{1}{\mu_0} \left( B_i B_j - \frac{B^2}{2} \delta_{ij} \right) $$
