The following solution is not entirely elementary, but is seems different from the approach cited by user1551. Moreover, the non-elementary part, namely, the part that identifies the points of maximum of a linear functional on the intersection of a cube and a hyperplane, can be made completely elementary, but somewhat detail-crowded.
As you have noticed, every symmetric matrix can be diagonalized by an orthonormal matrix, so the problem is reduced to the case when $A$ is a diagonal matrix. This is because if $A=Q^tDQ$ for some orthonormal matrix, then
$$\langle Au_j,u_j\rangle = \langle DQu_j,Qu_j\rangle$$
and $\{Qu_j\}_{j=1}^k$ is yet another choice of $k$ orthonormal vectors. Moreover, every choice of $k$ orthonormal $\{u_j\}_{j=1}^k$ vectors is obtained by applying $Q$ to the orthonormal vectors $Q^{-1}u_j$, so in effect the problem is reduced to the case when $A$ is a diagonal matrix. That is, we want to prove that whenever $D$ is a diagonal matrix whose diagonal is, say, $d_1,\dots, d_n$, with $d_1\geq d_2\cdots\geq d_n$, (*) then for every choice of orthonormal vectors $\{u_j\}_{j=1}^k$, with $k\geq 2$, we have:
$$\sum_{j=1}^k\langle Du_j,u_j\rangle \leq d_1+\cdots + d_k$$
We clearly have equality when $u_j=e_j$, where $e_j$ are the standard basis vectors.
For each $j$, write
$$u_j=(u_{j1},u_{j2},\dots,u_{jn})$$
Then for each $j=1,2,\dots k$:
$$\langle Du_j,u_j\rangle=\sum_{i=1}^n d_i u_{j i}^2$$
So, summing these up and interchanging the order of summation on the r.h.s:
$$
\sum_{j=1}^k\langle Du_j,u_j\rangle=\sum_{i=1}^nd_i\left(\sum_{j=1}^ku_{ji}^2\right)\enspace\enspace (1)
$$
Observe that since $\{u_j\}_{j=1}^k$ are orthonormal, we have:
$$\sum_{i=1}^n\left(\sum_{j=1}^ku_{ji}^2\right)=\sum_{j=1}^k||u_j||_2^2=k$$
and the summands are clearly non-negative, and by orthonormality are also bounded above by $1$. So if we put $x_i=\sum_{j=1}^ku_{ji}^2$, we can write (1) as $\sum_{i=1}^nd_ix_i$ with $x_i\in [0,1]$ and $\sum_{i=1}^nx_i=k$. Consider the maximum of the linear functional $\Lambda(x)=\sum_{i=1}^nd_ix_i$ subject to the constraints $\sum_{i=1}^kx_i=k$ and $x_i\in[0,1]$. These two constraints define a compact and convex subset of $R^n$, namely, the (obviously non empty) intersection of the cube $[0,1]^n$ and the hyperplane $\sum_{i=1}^nx_i=k$, so our linear functional attains a maximum at some extreme point of that set. Every extreme point of the intersection has all it's coordinates either zero or one, and since the sum is $k$, we have exactly $k$ ones and $n-k$ zeros. As $d_1\geq d_2\geq\cdots\geq d_n$, the maximum will be obtained when all the $1$'s appear at the first $k$ coordinates.
Going back to our problem, we deduce that, for every choice of orthonormal $u_1,\dots,u_k$:
$$\sum_{j=1}^k\langle Du_j,u_j\rangle\leq \sum_{i=1}^kd_i$$
In our particular problem the diagonal of $D$ consists of the eigenvalues of the matrix $A$, so for every choice of orthonormal $u_1,\dots,u_k$:
$$\sum_{j=1}^k\langle Au_j,u_j\rangle\leq\sum_{j=1}^k{\lambda_j}$$
Equality holds when $u_1=De_1,\dots, u_k=De_k$, i.e., when the $u_j$'s are orthonormal eigenvectors of $A$.
(*) We can assume without loss of generality that the elements of the diagonal of $D$ are non increasing, becuase if they are not, we conjugate $D$ by a suitable permutation matrix, which is also orthogonal, and the relevant inner product remains unchanged.
