A finite difference scheme, based upon the Crank–Nicolson scheme, is applied to the numerical approximation of a two-dimensional time fractional non-Newtonian fluid model. This model not only possesses a multi-term time derivative, but also contains a special time fractional operator on the spatial derivative. And a very important lemma is proposed and also proved, which plays a vital role in the proof of the unconditional stability. The stability and convergence of the finite difference scheme are discussed and theoretically proved by the energy method. Numerical experiments are given to validate the accuracy and efficiency of the scheme, and the results indicate that this Crank–Nicolson difference scheme is very effective for simulating the generalized non-Newtonian fluid diffusion model.