Influence of choice of statistics on gas kinetics In the derivation of distributions over energy states, a common assumption made is that under normal conditions (normal from a fluid dynamics standpoint, so > 300K typically) the energy states are very sparsely populated and so the vast majority of them remain empty. Because of this sparsity, the choice between Bose-Einstein and Fermi-Dirac statistics doesn't make a difference in the limiting case (the Boltzmann limit) when computing the distribution over the microstates. 
However, the validity of this limit comes down to the condition:
$$\frac{V}{N}\times\frac{(2\pi m k T)^{3/2}}{h^3} \gg 1$$
which can be violated if the mass is small and number density is large or when T is very low. When this happens, the choice of statistics matters. 
I have only seen derivations when the Boltzmann limit is satisfied and the choice doesn't matter. 
How do the results change if Fermi-Dirac statistics are used? Do we see changes in the macro-scale results such as the specific heat coefficients or the state relations? Or are the resulting equations the same but the specific values are different? Would macro-scale features such as shocks/expansion fans or chemical non-equilibrium behave differently? Or radiation?
 A: First, the derivation. The condition says that it's very unlikely for the particle number to be greater than one – the number of filled one-particle states is a small fraction of the number of all one-particle states and the doubly/multiply filled ones may be considered a higher-order correction that may be neglected. This approximate inequality is independent on which distribution we use: when the fraction of the occupied states is infinitesimal, the multiply occupied ones may be neglected and all the distributions are the same.
Now, at temperature $T$, the kinetic energy is of order $kT$, so the momentum $p$ is roughly given by $p^2/2m \sim kT$ i.e. $p\sim\sqrt{2mkT}$. We may view this value as the uncertainty of the momentum as well, so the uncertainty of the position – the width of the one-particle "box" in the $x$ or $y$ or $z$ direction – is given by $h/\Delta p$, approximately, which is $h/\sqrt{2mkT}$. The volume of the box in the position space is therefore $h^3/(2mkT)^{3/2}$, up to coefficients of order one (such as $\pi$ which differs in my expressions and yours).
If you divide the actual volume in the position space $V$ by this volume of the one-particle box, you get the number of boxes, and if you divide the number of particles by this ratio, you get the average number of particles per box which is $N/V\cdot h^3/(2mkT)^{3/2}$. This must be much smaller than one for the sparsity approximation to hold. You can see that the exponent in the numerator should be $+3/2$, not $-3/2$ as you wrote.
Of course, when the sparsity approximation doesn't work, it doesn't work. The Boltzmann distribution has to be replaced by the Bose-Einstein or Fermi-Dirac one, depending on the context, and all other physical observables that depend on the distribution are of course affected as well if the distribution has to be changed and the change is substantial.
A: I will try to answer it.
So for a simple kinetic theory we have the following derivation. 
The single particle partition function is:
$$Z_1 = \int \frac{d^3x d^3p}{(2\pi \hbar)^3} \exp{-(p_x^2 + p_y^2 + p_z^2)/(2mk_BT)} = V \left(\frac{(\pi*2mk_BT)^{1/2}}{2\pi\hbar}\right)^3$$
$$\implies Z_1 = V \left(\frac{mk_BT}{2\pi\hbar^2}\right)^2 = V/\lambda^3$$
This is a pretty standard derivation and the integrals are standard as well.
The factor of $2\pi\hbar$ can be proven from the phase space argument, but for now, just let's assume the results.
Now a partition function for N particles can be calculated as:
$$Z_N = \frac{Z_1^N}{N!}$$
Now everything boils down to calculating the chemical potential from the following equation:
$$\mu = \frac{\partial F}{\partial N} = - \beta^{-1} \frac{\partial \ln{Z_N}}{\partial N}$$
This is some algebra, but it is not very important for the discussion.
The most important thing is that the chemical potential becomes positive when:
$$V/\lambda^3 = 1$$
Now this can be interpreted as a way to say, that the particles occupy the whole space of the container and there is now an increase in the energy of the system if more particles are introduced.
This can be also interpreted differently, now the particles become much closer and their wavefunctions start to interact and their nature becomes important. They can be classified to either Fermions or Bosons.
Now to the point: the energy distributions for different systems can be calculated considering the Grand Partition function.
For a simple F system, the grand partition function can be written as follows:
$$ \Xi = \sum_n^1 \left(e^{-\beta(\varepsilon_k - \mu)}\right)^n = 1+e^{-\beta(\varepsilon_k - \mu)}$$
Then the particle number can be calculated as follows:
$$ \left< n \right>_{k,F} = \frac{\partial \Xi}{\partial \mu} = \frac{1}{e^{\beta(\varepsilon - \mu)}-1}$$
For B system we have a similar expression:
$$ \Xi = \sum_n^\infty (1-e^{-\beta(\varepsilon_k - \mu)})^n = \frac{1}{1-e^{-\beta(\varepsilon_k - \mu)}}$$
$$ \left< n \right>_{k,B} = \frac{1}{e^{\beta(\varepsilon - \mu)}+1}$$
I do not know if you'd like me to show how to calculate internal energy from this, but it shouldn't be difficult by making some approximations on the way.
You just need to consider the integral:
$$ U = \int_0^\infty g(\varepsilon) n(\varepsilon)\varepsilon d\varepsilon$$
where $g$ is the dispersion relation, which is different for many different systems, but can be easily obtained once the physical system in question is known.
And $n$ is the density of states for some given energy, which is related to the average particle number as follows:
$$\left<n\right>_k \frac{\partial k}{\partial \varepsilon}\delta\varepsilon = n(\varepsilon)\delta\varepsilon$$
