$\def\eqdef{\stackrel{\text{def}}{=}}$
The $\ S\ $ you're interested in is that of the last $\ m\ $ order statistics of the values shown on the dice's faces. The order statistics are those values $\ V_1\le$$\,V_2\le\dots\le V_n\ $ arranged in increasing order, and so
$$
S=\sum_{t=1}^mV_{n+1-t}\ ,
$$
or, putting $\ B_t\eqdef V_{n+1-t}\ $,
$$
S=\sum_{t=1}^mB_t\ .
$$
The expectation $\ E(n,k,m)\eqdef\Bbb{E}(S)\ $ can be evaluated using the linearity of expectation:
\begin{align}
\Bbb{E}(S)&=\sum_{t=1}^m\Bbb{E}\big(B_t\big)\\ &=\sum_{t=1}^m\sum_{r=1}^kr\,\Bbb{P}\big(B_t=r\big)\tag{1}\label{exp1}
\end{align}
and the probability mass functions $\ \Bbb{P}\big(B_t=r\big)\ $ of the random variables $\ B_t\ $.
In a separate section below, I derive an expression for the probability mass functions of $\ B_t\ $ and provide a Magma script for evaluating $\ E(n,k,m)\ $ using the expression \eqref{exp1}.
Since the order statistics aren't independent, it's not possible to obtain the probability mass function $\ P(n,k,m,x)\eqdef$$\,\Bbb{P}(S=x)\ $ of $\ S\ $ by convolving those of $\ B_m,$$\,B_{m-1},$$\,\dots,B_1\ $. It is possible, however, to express $\ P(n,k,m,x)\ $ in terms of the joint probability mass function of $\ B_m,B_{m-1},\dots,B_1\ $. While there's no simple closed form expression for it, it can nevertheless be expressed as a sum of products of powers, binomial coefficients and restricted composition numbers. The work necessary to evaluate this sum is a polynomial function of $\ n\ $ and so, when $\ n\ $ is large, it requires much less work than that needed to exhaust over all possible combinations, which grows exponentially with $\ n\ $.
The expression for $\ P(n,k,m,x)\ $ is
$$
\sum_{a=1}^{m-1}{n\choose a}\sum_{r=\max\left(1,\left\lceil\frac{x-ka}{m-a}\right\rceil\right)}^{\min\left(k,\left\lfloor\frac{x-a}{m-a}\right\rfloor\right)}p_{x-(m-a)r,a,r+1,k}\ \ \sum_{b=m-a}^{n-a}{n-a\choose b}\left(\frac{1}{k}\right)^{a+b}\left(\frac{r-1}{k}\right)^{n-a-b}
$$
if $\ x\ $ is not evenly divisible by $\ m\ $, where $\ p_{z,a,\ell,u}\ $ is the number of restricted compositions of length $\ a\ $ of the number $\ z\ $ into parts of sizes between $\ \ell\ $ and $\ u\ $ (inclusive), or
\begin{align}
\sum_{b=m}^n&{n\choose b}\left(\frac{1}{k}\right)^b\left(\frac{x-m}{mk}\right)^{n-b}+\\
&\sum_{a=1}^{m-1}{n\choose a}\sum_{r=\max\left(1,\left\lceil\frac{x-ka}{m-a}\right\rceil\right)}^{\min\left(k,\left\lfloor\frac{x-a}{m-a}\right\rfloor\right)}p_{x-(m-a)r,a,r+1,k}\ \ \sum_{b=m-a}^{n-a}{n-a\choose b}\left(\frac{1}{k}\right)^{a+b}\left(\frac{r-1}{k}\right)^{n-a-b}
\end{align}
if $\ x\ $ is evenly divisible by $\ m\ $.
The following Magma script evaluates $\ P(n,k,m,x)\ $ using the above expressions. You can run it by copying and pasting it into the online Magma calculator. For $\ n=4,k=6,m=3\ $ it agrees with the results obtained by Thomas Andrews in his answer. It also contains code to evaluate $\ E(n,k,m)\ $ using the definition
$$
E(n,k,m)=\sum_{x=m}^{km}xP(n,k,m,x)
$$
and obtains the same values for $\ E(n,k,m)$ Thomas B. does with the method he describes in his answer, except for the case $\ n=33, k=22\ $. When $\ n\ $ and $\ k\ $ are both that large, the computation of the restricted composition numbers $\ p\ $ exceeds the time limit imposed by the online Magma calculator.
restcomp:=function(x,n,l,u)
rcold:=ZeroMatrix(Integers(),1,Min(x,u*(n-1)));
rcnew:=ZeroMatrix(Integers(),1,Min(x,u*(n-1)));
if (x lt n*l) or (x gt n*u) then
return(0);
end if;
if n eq 1 then
return(1);
else
for xi in [l..Min(x,u)] do
rcold[1,xi]:=1;
end for;
for ni in [2..n-1] do
for xi in [ni*l..Min(x,ni*u)] do
for v in [Max(l,xi-(ni-1)*u)..Min(xi-1,u)] do
rcnew[1,xi]:=rcnew[1,xi]+rcold[1,xi-v];
end for;
end for;
rcold:=rcnew;
rcnew:=ZeroMatrix(Integers(),1,Min(x,u*(n-1)));
end for;
tot:=0;
for v in [Max(l,x-(n-1)*u)..Min(x-1,u)] do
tot:=tot+rcold[1,x-v];
end for;
return(tot);
end if;
end function;
P:=function(n,k,m,x)
pr:=0;
for a in [1..m-1] do
g1:=Max(1,Ceiling(Rationals()!(x-k*a)/Rationals()!(m-a)));
g2:=Min(k,Floor(Rationals()!(x-a)/Rationals()!(m-a)));
sc:=0;
for r in [g1..g2] do
s:=0;
for b in [m-a..n-a] do
s:=s+Binomial(n-a,b)*Binomial(n,a)*(r-1)^(n-a-b)/k^n;
end for;
sc:=sc+s*restcomp(x-(m-a)*r,a,r+1,k);
end for;
pr:=pr+sc;
end for;
if IsDivisibleBy(x, m) then
for b in [m..n] do
pr := pr+Binomial(n,b)*(x-m)^(n-b)/(m^(n-b)*k^n);
end for;
end if;
return( pr );
end function;
E:= function(n,k,m)
expv:=0;
for x in [m..k*m] do
expv:=expv+x*P(n,k,m,x);
end for;
return(expv);
end function;
for x in [3..18] do
print x, P(4,6,3,x), RealField(7)!P(4,6,3,x);
end for;
print E(2,6,1), E(2,20,1), E(4,6,3), E(4,6,1);
Probability mass function of $\ B_t $
The event $\ \big\{B_t=r\big\}\ $ occurs if and only if the uppermost faces of some number $\ a\ $ between $\ 0\ $ and $\ t-1\ $ of dice bear numbers strictly greater than $\ r\ $, the uppermost faces of some number $\ b\ $ between $\ t-a\ $ and $\ n-a\ $ of the remaining dice show exactly $\ r\ $, and the uppermost faces of the remaining $\ n-a-b\ $ dice bear numbers strictly less than $\ r\ $. The probability of this event is
$$
\sum_{a=0}^{t-1}{n\choose a}\sum_{b=t-a}^{n-a}\left(\frac{k-r}{k}\right)^a{n-a\choose b}\left(\frac{1}{k}\right)^b\left(\frac{r-1}{k}\right)^{n-a-b}\ ,
$$
which is therefore the value of $\ \Bbb{P}\big(B_t=r\big)\ $. In the above expression, $\ {n-a\choose b}\ $ is the number ways the $\ a\ $ dice to show numbers strictly larger than $\ r\ $ can be chosen, and $\ \left(\frac{k-r}{k}\right)^a\ $ is the probability that they will all in fact show numbers strictly larger than $\ r\ $, $\ {n-a\choose b}\ $ is the number of ways $\ b\ $ of the remaining $\ n-a\ $ dice to show exactly the number $\ r\ $ can be chosen, and $\ \left(\frac{1}{k}\right)^b\ $ is the probability that those $\ b\ $ dice will in fact all show the number $\ r\ $. Finally, $\ \left(\frac{r-1}{k}\right)^{n-a-b}\ $ is the probability that the remaining $\ n-a-b\ $, designated as showing numbers strictly less than $\ r\ $ will in fact all show numbers strictly less than $\ r\ $.
The following Magma script to evaluate $\ E(n,k,m)\ $ for arbitrary values of $\ n,k\ $ and $\ m\ $ using the expression \eqref{exp1}. It obtains the same values for $\ E(n,k,m)$ as the procedure given above using the probability mass function of $\ S\ $. Thanks to Thomas B. for picking up the fact that my $\ m\ $ is his $\ n-m\ $, and thus resolving an apparent discrepancy I had been concerned about earlier.
expv:=function(n,k,m)
ev:=0;
for t in [1..m] do
evbt:=0;
for r in [1..k] do
pbter:=0;
for a in [0..t-1] do
s:=0;
for b in [t-a..n-a] do
s:=s+Binomial(n-a,b)*Binomial(n,a)*(r-1)^(n-a-b);
end for;
pbter:=pbter+s*(k-r)^a/k^n;
end for;
evbt:=evbt+r*pbter;
end for;
ev:=ev+evbt;
end for;
return( ev );
end function;
print expv(2,6,1), expv(2,20,1), expv(4,6,3), expv(4,6,1);
print expv(33,22,22), RealField(10)!expv(33,22,22);
print expv(2,6,1), expv(2,20,1), expv(4,6,3), expv(4,6,1);
print expv(33,22,22), RealField(10)!expv(33,22,22);
Addendum
I agree with the expression you derive for $\ E(n,k)\eqdef$$\,E(n,k,n-1)\ $ in your answer. I believe it's possible to derive your expression by substituting the one I give above for $\ \Bbb{P}\big(B_t=r\big)\ $ into the sum $\ \sum_\limits{t=1}^{n-1}\sum_\limits{r=1}^kr\,\Bbb{P}\big(B_t=r\big)\ $. The only way I can see of proving the identity algebraically, however, is to transfer the sum over $\ t\ $ across the other three sums in the expression so that it becomes the inner sum instead of the outer sum. The algebraic manipulations necessary to do this are rather cumbersome and error-prone, and since I'm satisfied your expression is correct, I eventually gave up trying to prove it by this unnecessarily convoluted procedure. I did however write a Magma function to evaluate your expression for $\ E(n,k)\ $ and compare its output with that of $\ E(n,k,n-1)\ $ as produced by the above Magma script. For every pair $\ n,k\ $ I tried, the results were identical. Here they are:
\begin{align}
E(2,6)&=\frac{161}{36}=E(2,6,1)\\
E(2,20)&=\frac{553}{40}=E(2,20,1)\\
E(4,6)&=\frac{15869}{1296}=E(4,6,3)\\
E(33,22)&=\frac{
56692906323757439755684138687075800338253661}{149889224159230635544290231950914661384192}\\
&=E(33,22,32)\ .
\end{align}
