Invariant differential operators on real Grassmannians

In this particular case, the invariant theory is pretty simple, so there is an explicit description.

Recall that, in general, for a homogeneous space $M=G/H$, a $G$-invariant, $m$-th order linear differential operator $L:C^\infty(M)\to C^\infty(M)$ is specified by its value on the $m$-jets at the identity coset $eH$ of functions in $C^\infty(M)$. In particular, the vector space of such operators is dual to the vector space of $\mathrm{Ad}(H)$-invariant elements in $J^m_{eH}(M,\mathbb{R})$. Taking the inverse limit as $m$ goes to infinity, one sees that the vector space $\mathcal{L}$ of $G$-invariant, linear differential differential operators from functions to functions is dual to the vector space of $\mathrm{Ad}(H)$-invariant elements of the symmetric algebra generated by ${\frak{m}}^*$, where ${\frak{m}}=T_{eH}M$. In other words, as a vector space, $\mathcal{L}$ is isomorphic to the $\mathrm{Ad}(H)$-invariant polynomials on ${\frak{m}}^*$.

In the case of the real Grassmannian $\mathrm{O}(n)/\bigl(\mathrm{O}(k)\times\mathrm{O}(n{-}k)\bigr)$, you are asking for the $\mathrm{O}(k)\times\mathrm{O}(n{-}k)$-invariant polynomials on the $k$-by-$(n{-}k)$ matrices, and this follows from the singular value decomposition for such matrices. Supposing, as one can, that $k\le n{-}k$, if one lets $A$ be such a matrix, then one can define polynomials $p_j(A)$ (of degree $2j$) by the rule $$ \det(A A^T - t\ I_k) = p_k(A) - p_{k-1}(A) t + p_{k-2}(A) t^2 - \cdots + (-1)^kp_0(A)\ t^k. $$ Then these $k{+}1$ polynomials generate the ring of $\mathrm{Ad}(H)$-invariant polynomials. (Of course, $p_0(A) \equiv 1$.) Thus, $p_j$ corresponds to an invariant linear differential operator $L_{2j}$ of order $2j$, and these generate the ring of invariant operators on the Grassmannian.