Why do we use grand canonical ensemble for BEC description? The system we consider has constant $N$, $V$ and $T$ (the number of particles, volume and temperature) This is just the thermodynamic variables for the canonical ensemble, why we use fugacity $z$ or chemical potential $\mu$ there?
 A: If you choose to set $N$, $V$ and $T$, then you are choosing the Canonical Ensemble, so $N$ is fixed. In this approach, the chemical potential $\mu$ is just a parameter that you've to determine by imposing that the total number of particles is $N$. 
Viceversa, if you choose the grandcanonical ensemble, you fix $\mu$ and so the number of particles $N$, fluctuates around an average value $\langle N \rangle$. I hope that I've well understood your problem and answered properly.
A: In the thermodynamic limit fluctuations in particle number should tend to $0$ and so we should get the same result whether we use the cannonical or grand canonical ensamble. If we are ultimately interested in a system with a fixed particle number, we can, if we so choose, set up the problem in the grand cannonical ensamble and then invert the equation for $\langle N \rangle$ to find what chemical potential will give us the desired mean particle number. 
For systems of identical quantum particles we generally choose to use the grand cannonical ensamble because it turns out to be dramatically simpler. In such systems it is convenient to use an occupation number representation for states, rather than listing the state of each individual particle, since it can then be assumed that an appropriate (anti-)symmeturization procedure has been applied and not have to think about it. In the cannonical ensamble, however, we now encounter a problem; since the total number of particles is fixed the occupaion numbers of differnet states are not independant, but must satisfy a constraint. This leads to a rather convoluted combinatorics problem to list the set of allowed ocupation numers.
In the grand cannonical ensamble we do not encounter this problem, since all values of $N$ are allowed. We simply consider all possible sets of occupation numbers, but since states with very large or very small particle number are exponentially weighted against, they do not contirbute. In effect we allow the statistics of the problem to "automatically" find the states with the correct particle number for us. 
