This work devotes to solving a class of delay fractional partial differential equations that arises in physical, biological, medical, and climate models. For this, a numerical scheme is implemented that applies operational matrices to convert the main problem into a system of algebraic equations; then, solving the resultant system leads to an approximate solution. The two-variable Chebyshev polynomials of the sixth kind, as basis functions in the proposed method, are constructed by the one-variable ones, and their operational matrices are derived. Error bounds of approximate solutions and their fractional and classical derivatives are computed. With the aid of these bounds, a bound for the residual function is estimated. Three illustrative examples demonstrate the simplicity and efficiency of the proposed method.