compute expected value of a expression with normal distribution how compute expected value of the below expression
$$
E[S^TS]
$$
in which:
$$
S = \sum_{i=1}^{\nu} \mathbb{u}_i \mathbb{u}_i^{\top} 
$$
and $\mathbb{u}_i \sim \mathcal{N}(0,V); V \in R^{p \times p}$ is covariance matrix
 A: [I have used $p$ differently: as a running index]
First off, $S = S^T$ because of the definition of $S$.
$E[S^TS]=E[SS]=E[(\sum_i u_i u_i')(\sum_j u_ju_j')]$
$=E[\sum_{i,j} u_i u_i' u_j u_j']=\sum_{i,j} E[u_i u_i' u_j u_j']=E[(u_i \cdot u_j) u_i u_j']$
$\Rightarrow E[(u_i u_i' u_j u_j')_{pq}]=\sum_k E[u_{ik}u_{jk} u_{ip} u_{jq}]$
Each scalar here is the $(p,q)$th element of the result you are looking for.
So each scalar we need to calculate is similar to a 4th moment.
Now, let $i \neq j$. Then, by marginal probability:
$p(u_{ik},u_{jk}, u_{ip}, u_{jq})=\int p(u_{i1},\ldots,u_{in}, u_{j1},\ldots,u_{jn}) du_\text{others}$
(applying independence of sample $i$ from $j$)
$=\int p(u_{i1},\ldots,u_{in}) d(\text{not }u_{ik} \text{ or } u_{ip}) \int p(u_{j1},\ldots,u_{jn}) d(\text{not }u_{jk} \text{ or } u_{jq})=P_i(u_{ik},u_{ip})P_j(u_{jk},u_{jq})$
We can use the joint rule for the multivariate normal to evaluate each.
$P_i = \text{normpdf}([u_{ik},u_{ip}]|0,V([k,p],[k,p]))$
$P_j = \text{normpdf}([u_{jk}, u_{jq}]|0,V([k,q],[k,q]))$
I am using MATLAB notation here, so that $V([a,b],[a,b])=\begin{bmatrix}V_{aa} &V_{ab} \\V_{ba}  &V_{bb} \end{bmatrix}$
So, now that we have the joint over the 4 variables, we can calculate each expectation and perform the sums:
$E[u_{ik}u_{ip}u_{jk}u_{jq}]=\int u_{ik}u_{ip}u_{jk}u_{jq} P_i(u_{ik},u_{ip}) P_j(u_{jk},u_{jq})du_{ik}du_{ip}du_{jk}du_{jq}$
This would involve the tracking of several cases such as $k=i$, $k=p$, $i=j$, etc. and would take some effort to workout. I suspect that things should simplify as you work them out. But in principle at least, you should be able to calculate the final expression.
