In this paper, we compute the constrained QCD effective potential up to two-loop order with finite quark mass and chemical potential. We present the explicit calculations by using the double line notation and analytical expressions for massless quarks are obtained in terms of the Bernoulli polynomials or Polyakov loops. Our results explicitly show that the constrained QCD effective potential is independent on the gauge fixing parameter. In addition, as compared to the massless case, the constrained QCD effective potential with massive quarks develops a completely new term which is only absent when the background field vanishes. Furthermore, we discuss the relation between the one- and two-loop constrained effective potential. The surprisingly simple proportionality that exists in the pure gauge theories, however, is in general no longer true when fermions are taken into account. On the other hand, for high baryon density $mu_B$ and low temperature $T$, in the massless limit, we do also find a similar proportionality between the one- and two-loop fermionic contributions in the constrained effective potential up to ${cal O}(T/mu_B)$.