Energy of classical ideal gas in the grand canonical ensemble The canonical partition function for an ideal gas is
$$
Z(N,V,\beta) = \frac{1}{N!} \left(\frac{V}{\lambda^3}\right)^N
$$
where $\lambda = \sqrt{\frac{\beta h^2}{2 \pi m}}$ is the thermal De-Broglie wavelength. It is straightforward to obtain
$$
\langle E \rangle = -\frac{\partial \log Z}{\partial \beta} = \frac{3}{2} N k_B T .
$$
From $Z$ the grand-canonical partition function is
$$
Q(\mu,V,\beta) = \sum_{N=0}^\infty \frac{1}{N!} \left(\frac{e^{\beta \mu} V}{\lambda^3}\right)^N = e^{\frac{e^{\beta \mu} V}{\lambda^3}} .
$$
The average particle number is
$$
\langle N \rangle = \frac{\partial \log Q}{\partial (\beta \mu)} = \frac{e^{\beta \mu} V}{\lambda^3} .
$$
To get the average energy we should do, substituting $\langle N \rangle$,
$$
\langle E \rangle = - \frac{\partial \log Q}{\partial \beta} = \frac{3}{2} \langle N \rangle k_B T
$$
but this is true only if we magically ignore the $e^{\beta \mu}$ factor when taking the derivative, otherwise there is an extra (nonsensical) term. I've checked a few sources and this is the accepted solution (after all, it must be this one to be consistent with the canonical ensemble result), although they mysteriously gloss over the issue, so I'm missing something. Thanks.
 A: Hint: You have an error in your computations.  In particular in the grand canonical ensemble,
\begin{align}
  \langle E \rangle \neq -\frac{\partial \log Q}{\partial \beta}.
\end{align}
Moreover, I just did the whole computation having corrected this error in the appropriate way, and it worked out the way it should.
Addendum, 2019-02-02. Details Beyond the Hint
Step 1. Recall the following definitions of the grand canonical partition function $Q$, the ensemble average energy $\langle E\rangle$ in the grand canonical ensemble, and the ensemble average particle number $\langle N\rangle$ in the grand canonical ensemble. All sums are over states $i$ of the system:
\begin{align}
  Q \equiv \sum_ie^{-\beta(E_i - \mu N_i)}, \qquad \langle E\rangle \equiv \sum_i \frac{e^{-\beta(E_i - \mu N_i)}}{Q}E_i, \qquad \langle N\rangle \equiv \sum_i \frac{e^{-\beta(E_i - \mu N_i)}}{Q}N_i
\end{align}
Step 2. Show that the following identity follows from the definitions in Step 1:
\begin{align}
  \langle E\rangle = -\frac{\partial \ln Q}{\partial \beta} + \mu\langle N\rangle.
\end{align}
Step 3. Show that if we take
\begin{align}
  Q = V\frac{e^{\beta\mu}}{\lambda^3},
\end{align}
then
\begin{align}
  -\frac{\partial \ln Q}{\partial \beta} = \frac{3}{2}\frac{\langle N\rangle}{\beta} - \mu\langle N\rangle.
\end{align}
Step 4. Combine steps 2 and 3 to obtain the desired result.
A: Yes, there is a subtle issue. Recall the general problem, from which those formulae are derived. 
An example first. Thermodynamics say that
$$U=\Phi+TS+\mu N$$
($\Phi$ is the grand-canonical potential)
from which you could obtain 
$$S=\left( -\frac{\partial \Phi}{\partial T} \right)_{V,\mu}$$
And check that, in thermodynamics, we always indicate a dubindex in the derivatives. That subindex $(V,\mu)$, is reminding us that those variables must be kept constsant.
If you think about it, they shoudln't be neccesary, because a partial derivative already assumes implicitly that "all other variables" are kept constant. But thermodynamics can become quite messy, so it is actually really important to keep track on the variables that are kept constant.
Now, if you get this, let's get to the point. When dealing with ensembles, we use the entropic representation. That is working with $S$ as fundamental relation, instead of $U$. Consequently, you must rearrange the expression:
$$S=-\frac{\Phi}{T} + \frac{U}{T} + \frac{\mu \bar{N}}{T}$$
And now, variables are obtained as always:
$$U=\frac{\partial \Phi}{\partial (\lambda_U)}$$
where $\lambda_U$ is the intensive variable associated with $U$, which is $(1/T)$. So basically you have to calculate:
$$U=\frac{\partial \Phi}{\partial (1/T)}$$
But this partial derivative must be done maintaining the other variables constant. That includes maintaining constant $(\mu/T)$
Yes, it is weird: $1/T$ varies, but $\mu/T$ must not vary. That's the thing.
... and since $\Phi=-k_B T \ln {(Q)}$
$$U=\frac{\partial [-k_B\ln{(Q)}]}{\partial (1/T)}$$
And dividing everything by $k_B$, you get
$$U=\frac{\partial [-\ln{(Q)}]}{\partial [1/(k_BT)])}=-\frac{\partial \ln{(Q)}}{\partial \beta}$$
But this derivative is done maintaining $(\mu/T$) constant. If that is constant, dividing by $k_B$ will keep being constant, and hence we write
$$U=-\left(\frac{\partial \ln{(Q)}}{\partial \beta}\right)_{(\beta\mu)}$$
That is, the derivative must be done maintaining $\beta\mu$ constant. That's the explanation of your calculus.
There are interesting articles on the lack of rigour of this notation haha.
A: The expression you chose for $\left<E\right>$ is not consistent with a temperature-independent chemical potential! 
To find its dependence, recall that the constants like $\beta$ arise from solving the stationary points of the entropy functional:
$$
S = -\int \mathrm{d} \Gamma \rho(Q,P) \ln{ \left( \rho \left(Q,P \right) \right) }
$$
subject to various conditions, which give rise to Lagrange multipliers, that are later identified with various thermodynamic properties:


*

*$\lambda = -\ln{(Z_G)}-1$ from normalization constraint;

*$\beta := T^{-1} $ from constant $\left<E\right> ;$

*$\gamma := -\mu\beta$ from constant $\left<N\right> ;$
where $$Z_G = \sum_{N=0}^\infty \int \mathrm{d} \Gamma_N e^{-H\beta -\gamma N}$$
It is now clear that the fugacity $e^{\mu\beta} = e^{-\gamma}$  is truly independent of $\beta$.
Now the expectation value of energy is given by 
$$-\frac{\partial(\ln Z)}{\partial\beta}$$
with the understanding that $\gamma$ is constant. This is the same as the formula for the canonical ensemble.
However, the usual treatment is to treat $\gamma$ as the dependent variable and $\mu$ as independent, since $\mu$ has a more obvious physical interpretation.
Therefore the average energy is 
$$
\begin{align}
\left<E\right> & ~=~ \frac{1}{Z_G}\sum_{N=0}^\infty \int \mathrm{d} \Gamma_N H e^{-H\beta + \beta\mu N} \\[5px]
 & ~=~ \frac{1}{Z_G}\sum_{N=0}^\infty \int \mathrm{d} \Gamma_N \left(-\frac{\partial}{\partial\beta} + \mu N \right) e^{-H\beta + \beta\mu} \\[5px]
 & ~=~ \mu \left<N\right> - \frac{\partial}{\partial\beta}\left(\ln{Z_G}\right)
\end{align}
$$
where the partial derivative holds $\mu$ constant. 
So depending on whether you choose to hold $\gamma$ or chemical potential as your independent constant, you get a different formula for $\left<E\right> .$ The end result of your calculation should be independent of this choice, as long as you are consistent with all your derivations.
This explains the "magic" in ignoring the $e^{\beta\mu}$ factor. Holding its argument constant is consistent with your expression for $\left<E\right>. $ The "correct" (conventional) Grand canonical ensemble expression contains a term that exactly cancels the "extra (nonsensical) term".
A: Commenting on FGSUZ answer, the somewhat odd result that
$$
  \langle E \rangle = -\left(\frac{∂\log Q}{∂β}\right)_{βμ}
$$
is not an accident, in fact it has a quite deep interpretation in terms of probability theory.
Consider the grand canonical distribution
$$
  \begin{aligned}
  &w_{k} = \frac{1}{Q}e^{β(μN_k - E_k)}\\
  \\
  &\text{where }Q = \sum_k e^{β(μN_k - E_k)}
  \end{aligned}
$$
giving the probability of the state $k$ with energy $E_k$ and $N_k$ particles at the given $β,μ$.
In probability theory, given a random variable $E$, one defines its moment generating function as the expectation value of $e^{sE}$, $M(s) = \mathbb{E}[e^{sE}]$.
By expanding the exponential we see why:
$$
  M(s) = \mathbb{E}[e^{sE}]
       = \mathbb{E}\left[\sum_{n=0}^{+∞} \frac{(sE)^n}{n!}\right]
       = \sum_{n=0}^{+∞} \mathbb E[E^n] \frac{s^n}{n!}
$$
So, the n-th derivative of $M$ evaluated at $s=0$ gives the n-th moment $m_n = \mathbb E[E^n]$ of $E$.
In the case of the grand canonical distribution, $M$ can be expressed in terms of partition function $Q(β, μ)$ alone¹:
$$
\begin{aligned}
  M(s) &= \sum_k w_k e^{sE_k} \\
       &= \frac{1}{Q} \sum_k e^{βμN_k - (β-s)E_k} \\
       &= \frac{1}{Q} \sum_k e^{(β-s)/(β-s) ⋅ βμN_k - (β-s)E_k} \\
       &= Q(β-s, βμ/(β-s)) \: \big/ \: Q(β,μ)
\end{aligned}
$$
Taking derivatives of this is quite annoying because $β$ and $μ$ appear mixed in the second argument.
It's beneficial to do a change of variables to decouple them: define $λ=e^{βμ}$, the so-called fugacity².
The partition function then becomes
$$
  Q(β, λ) = \sum_k e^{βμN_k} e^{-βE_k} = \sum_k λ^{N_k} e^{-βE_k}
$$
and $M(s) = Q(β-s, λ)/Q(β,λ)$, which is much easier to deal with.
We can now easily compute the first moment, that is, the average energy:
$$
  \begin{aligned}
    \langle E\rangle
        = M'(0) &= \frac{∂}{∂s} \left[\frac{Q(β-s, λ)}{Q(β,λ)}\right]_{s=0} \\
                &= -\frac{1}{Q} \frac{∂Q}{∂β}(β,λ) \\
                &= -\left(\frac{∂\log Q}{∂β}\right)_λ
  \end{aligned}
$$
For comparison, after a tedious computation $Q(β,μ)$ would have given the formula:
$$
  \langle E\rangle = -\left(\frac{∂\log Q}{∂β}\right)_μ
                     +\frac{μ}{β}\left(\frac{∂\log Q}{∂μ}\right)_β
$$
which is completely equivalent to the above, but uglier and less useful.

Bonus: higher moments can't be so easily expressed as derivatives wrt $β$, but we can compute the
cumulants instead.
These $k_n$ are combinations of the usual moments $m_n$ and are equally good in describing a probability distribution.
They are generated by the function $\log M(s)$, as before:
$$
\begin{aligned}
  K(s) &= \log M(s) = \sum_{n=1}^{+∞} k_n \frac{s^n}{n!} \\
  k_n &= K^{(n)}(s)
\end{aligned}
$$
This choice is very convenient because we have (up to the inconsequential constant $-\log Q(β, λ)$)
$$
  K(s) = \log Q(β-s, λ)
$$
From this all the cumulants (the first three being the average $\langle E\rangle$, fluctuations $\langle ΔE^2\rangle$ and $\langle ΔE^3\rangle$) follow easily:
$$
  k_n =  \left(\frac{∂^n\log Q}{∂β^n}\right)_λ
$$

¹ This is just a consequence of $w_k$ being exponential.
² The exponential is not important, but it's conventionally defined like this and a little bit nicer.
A: You may take derivative by assuming constant fugacity i.e $e^{\beta\mu}$. Because the Main formula for calculation of partition function in grand canonical ensemble is $Q=\sum\,e^{-\beta(E-\mu N)}$. So when you want calculate mean Energy or $\langle E\rangle =\sum p_i E_i =\sum\frac{e^{-\beta(E_i-\mu N_i)}}{Q}E_i $. 
We know that $p_i=\frac{e^{-\beta(E_i-\mu N_i)}}{Q}$. so you should take derivative of $\hbox{ln}(Q)$ in a way that for each term just $E_i$ come beside $p_i$ it mean you should take derivative in constant fugacity or constant $e^{\beta\mu}$
