How can $mu$ be nonzero if grand canonical density operator commutes with $hat{H}$?

Short answer: you incorrectly deduced that the Hamiltonian does not commute with the particle number operator.

Here is the flawed portion:

Now, for a grand canonical ensemble, the particle number is not conserved due to its contact with a particle reservoir, [?̂,?̂]≠0.

I will explain the mistake in the above quote in two different ways. Please keep in mind that ?̂ in the original question is the Hamiltonian for the system, not for the system + reservoir. This is important.

1. First way: by analogy. The mistake is the same as in the following parallel argument (note that I simply replaced "particle number" with energy/Hamiltonian in the above quote):

Now, for a grand canonical ensemble, the energy is not conserved due to its contact with a particle (and heat) reservoir, so [?̂,?̂]≠0.

Of course the conclusion is absurd, the Hamiltonian commutes with itself! So the logic is flawed, and the flaw is exactly the same here as in the original quote.

2. Second way: correct analysis of the commutator. If you actually let the system evolve according to ?̂ (by disconnecting it from the reservoir), then the particle number will be conserved (for the case we're discussing). As usual, if a quantity is conserved, then the corresponding operator commutes with the Hamiltonian, so [?̂,?̂]=0.

Let me first note that $[H,N]=0$ is a purely mathematical relation, which can be verified for the Hamiltonian that you are working with. Derivations in (quantum) statistical mechanics textbooks often assume that it is the case, while leaving the exceptions from the rule to the chapters about superconductivity and superfluidity. Another purely mathematical fact is that $mu$ is a number, which commutes with either $H$ or $N$.

The quantum mechanical meaning of equality $[H,N]=0$ is that the energy and the particle number can be used together as quantum numbers to label the states. (There are obviously more quantum numbers, but these two are privileged in the context of the statistical mechanics.) Switching from the canonical to the grand canonical ensemble doesn't change anything in the Hamiltonian or its states, but changes the states that we include in the statistical sum: in the canonical ensemble they all have the same $N$, whereas in the grand canonical ensemble different values of $N$ are allowed.

In other words non-conservation of particles has different meaning in QM and Stat. Mech.: in the former it means that the energy states cannot be characterized by a particle number, in the latter - that we consider states with different numbers of particles.

Finally, note that this is not a quantum mechanical feature: although there is no semantic confusion in the classical case, the classical Hamiltonians do conserve the number of particles, as they have a fixed number of degrees of freedom. Thus, when we work in the grand canonical ensemble, we use a different Hamiltonian for every $N$, e.g., $$ H_N=sum_{i=1}^Nfrac{mathbf{p}_i^2}{2m} $$