Does Frobenius norm (not operator 2 norm) preserve the positive semidefinite order of matrices?

Yes. Since $0\preceq A\preceq B$, we have $0\le\lambda_i(A)\le\lambda_i(B)$ when the eigenvalues of $A$ and $B$ are arranged in the same (ascending or descending) order. It follows that $\operatorname{tr}(A^2)=\sum_i\lambda_i(A)^2\le\sum_i\lambda_i(B)^2=\operatorname{tr}(B^2)$.

(However, despite $\operatorname{tr}(B^2)\ge\operatorname{tr}(A^2)$, the difference $B^2-A^2$ is not necessarily PSD.)


Here is an alternative proof that is more direct (I think). Since $B - A$ is symmetric positive semidefinite, we know that eigenvalues of $(B - A)$ are nonnegative. In other words, $\mathrm{Tr}((B - A)^2)\ge 0$. Expanding this using the fact that the trace is linear and $\mathrm{Tr}(CD) = \mathrm{Tr}(DC)$, we get \begin{align*} \mathrm{Tr}(B^2)\ge 2\mathrm{Tr}(AB) - \mathrm{Tr}(A^2) & = \mathrm{Tr}(A^2) + 2\Big[\mathrm{Tr}(AB) - \mathrm{Tr}(A^2)\Big] \\ & = \mathrm{Tr}(A^2) + 2\mathrm{Tr}(A(B - A)). \end{align*} The second term is nonnegative since $A$ and $B - A$ are both symmetric positive semidefinite, see here for the proof of this fact.