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 (but each state in itself may have a well-defined number 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}^N\frac{\mathbf{p}_i^2}{2m}
$$
Update in response to the bounty message
The flawed part of the argument in the OP is indeed this:
Now, for a grand canonical ensemble, the particle number is not conserved due to its contact with a particle reservoir, $$[\hat{H},\hat{N}]\neq 0,$$ the only option is $\mu=0$.
The eigenstates of the Hamiltonian can be characterized by quantum number $n$, corresponding to the number of particles, and other quantum numbers - describing the state of n-particle system, which I denote by $j$. Thus, we have eigenstates
$$H|j,n\rangle = \epsilon|j,n\rangle.$$
These states are also eigenstates of the particle number operator, with quantum number $n$:
$$N|j,n\rangle=n|j,n\rangle.$$
We see that the Hamiltonian and the number operator have the same eigenstates and therefore commute - one could directly verify this by calculating the matrix elements of the commutator:
$$\langle j,n|[H,N]|j',n'\rangle=0.$$
Now, in canonical ensemble we would fix the number of particles calculate averages and the partition function using the subset of the eigenstates corresponding to this number of particles, i.e., our density matrix would be
$$
\rho^{n}_{j,j'}=\frac{1}{Z_n}\langle j,n|e^{-\beta H}|j',n\rangle, Z_n=\sum_{j}\langle j,n|e^{-\beta H}|j,n\rangle=\text{tr}_n e^{-\beta H}
$$
(where $j,j'$ do not necessarily refer to the Hamiltonian eigenstates, but the only states used have the same particle number $n$.)
In the Grand canonical ensemble we would consider probabilities of finding the system in states with different numbers of particles:
$$
\rho_{j,n;j',n'}=\frac{1}{Z}\langle j,n|e^{-\beta (H-\mu N)}|j',n'\rangle,\\
Z = \sum_j\sum_n \langle j,n|e^{-\beta (H-\mu N)}|j,n\rangle=
\text{tr}e^{-\beta (H-\mu N)}=\sum_ne^{\beta\mu n}\text{tr}_n e^{-\beta H}
$$
If we are in an eigenstate with a certain number of particles, this number would not change with time - it is conserved. However, in the grand canonical ensemble we have finite probability of eigenstates with different numbers of particles. Using term conserved to describe averaging over different states, each of which conserves the number of particles, is misleading.