Let us recall Quadratic Mean - Arithmetic Mean - Geometric Mean inequalities.
For $n\in\mathbb{N}$ and positive real numbers $a_i$, $i=1,\ldots,n$, it holds
$$\sqrt{\frac{a_1^2+a_2^2+\ldots+a_n^2}{n}}\geq \frac{a_1+a_2+\ldots+a_n}{n}\geq \sqrt[n]{a_1a_2\cdots a_n},$$ with equality of and only if $a_1=a_2=\ldots =a_n$.
The proposed problem is well-posed for $n=1$ and the solution follows directly from the QM-GM inequality: $$1=\sqrt{\frac{x^2+y^2+z^2}{3}}\geq\sqrt[3]{xyz}.$$
Let $n\geq 2$. Condition $$\sum_{i=1}^n(x_i^2+y_i^2+z_i^2)=3n$$can be written as $$X+Y+Z=3n,\quad X=\sum_{i=1}^nx_i^2,\quad Y=\sum_{j=1}^ny_j^2,\quad Z=\sum_{k=1}^nz_k^2.$$
Without the loss of generality let $X\geq Y\geq Z$.
Function
$$f(x_1,\ldots,x_n,y_i,\ldots,y_n,z_1,\ldots,z_n)=\sum_{n|i+j+k}x_iy_jz_k$$
cannot attain its maximal value if some $x_i'y_j'z_k'$ is negative, so it is sufficient to only consider the case $x_i,y_j,z_k\geq0$.
Let us start with the case when all numbers are positive.
For fixed $i\in\{1,\ldots,n\}$ we notice that $x_i$ appears in exactly $n$ terms in $f$ and if $x_iy_jz_k$ is the part of $f$ for some $j$ and $k$, then $x_iy_kz_j$ is also the part of $f$.
By $\displaystyle{s_{x_i}}$ we denote the sum of all terms in $f$ which contains $x_i$.
Then
\begin{equation}\label{eq1}s_{x_i}=x_i\left(\sum_{j+k=n-i} y_jz_k+\sum_{j+k=2n-i} y_jz_k+\sum_{j+k=3n-i} y_jz_k\right),\quad j,k\in\{1,\ldots,n\}.\end{equation}
To each term of each sum in $s_{x_i}$ we apply squared GM-QM inequality and obtain
$$s_{x_i}\leq \frac{x_i}{2}\left(\sum_{j+k=n-i}(y_j^2+z_k^2) +\sum_{j+k=2n-i} (y_j^2+z_k^2)+\sum_{j+k=3n-i} (y_j^2+z_k^2)\right)=\frac{x_i}{2}(Y+Z)=\frac{x_i}{2}(3n-X).$$
Since $$s_{x_1}+s_{x_2}+\ldots+s_{x_n}=\sum_{n|i+j+k}x_iy_jz_k,$$
we have
$$\sum_{n|i+j+k}x_iy_jz_k\leq\frac{3n-X}{2}(x_1+x_2\ldots+x_n).$$
If we have $x_1,\ldots,x_p>0$ and $x_{p+1}=\ldots=x_n=0$, for some $1\leq p<n$, then
$$\sum_{n|i+j+k}x_iy_jz_k\leq\frac{3n-X}{2}(x_1+x_2+\ldots+x_p)=\frac{p}{2}\frac{x_1+\ldots+x_p}{p}(3n-X)\leq \frac{p}{2}\sqrt{\frac{X}{p}}(3n-X),$$
where the last inequality follows from AM-QM inequality.
We get
$$\sum_{n|i+j+k}x_iy_jz_k\leq\frac{1}{2}\sqrt{pX}(3n-X)<\frac{1}{2}\sqrt{nX}(3n-X).$$
Let us define $$g(X):=\frac{1}{2}\sqrt{nX}(3n-X).$$
The first derivative of $g$ is
$$g'(X)=\frac{3\sqrt{n}}{4}\frac{(n-X)}{\sqrt{X}}.$$
From the fact that $X\geq n$ (otherwise, from $X\geq Y\geq Z$ we would have $X+Y+Z<3n$) we conclude that $g'(X)\leq 0$ so $g$ is decreasing function.
Since $g'(X)=0$ if and only if $X=n$, we conclude that $g$ attains its maximum value at $X=n$ which is equal to $n^2$.
Finally, we conclude that $$\sum_{n|i+j+k}x_iy_jz_k\leq n^2.$$
We can easily find real numbers for which equality holds. This will be the case when all $x_i$, $y_j$ and $z_k$ are equal to one, but there are other cases too. For example all $x_i$ and $y_j$ are equal to $-1$, and all $z_i$ are equal to 1...
Remark: the number of terms of the sum in $f$ can be easily calculated with the help of generating functions. We are searching for the total number of solutions of equations $i+j+k=n$, $i+j+k=2n$ and $i+j+k=3n$, where $1\leq i,j,k\leq n$. After we build generating function $F$ of variable $a$, we need to sum up the coefficients of $a^n$, $a^{2n}$ and $a^{3n}$ in $F$. The result will be $n^2$.