Is there mathematical relationship between the FFT and DCT transforms
You have some confusion in  the indices, maybe the correct expression would be
$$Y[k] = \sum_{n=1}^N \left[ 
\sum_{m=1}^N x[m] \cos\left( \frac{(2m+1)n\pi}{2N} \right) \right] e^{-\frac{j2\pi k n}{N}}$$
The DCT is the real part of the odd-indices coefficients of double length DFT
Expressing DCT in terms of complex harmonics
$$z[n] = \sum_{m=1}^N x[m] \left( \exp\left( j \frac{(2m+1)n\pi}{2N} \right) + \exp\left( -j \frac{(2m+1)n\pi}{2N} \right) \right)$$
If you take $X(n)$ being the fourier transform of $x$, in the reals instead of the integers, you could write $z[n] = X(m+1/2)+X(-m-1/2)$, $X(n)$ has Hermitian symmetry and the result agrees with the fact that $z[n]$ is real.
I think that the most efficient way to express the relation between $Y[m]$ and $x[m]$, is that $Y = DFT(DCT(x))$, unless you want to upsample your signal.
The denominator $2N$ in the exponential, the summation is from $1$ up to $N$, and the numerators may be even, so cannot simplify the exponent, the way to go is to compute a $2N$-samples, DFT.
If you create an intermediate signal with odd samples $v_{2i+1} = x_i + x_{N - i}$ and even samples $v_{2i}=0$, and compute $V = DFT^{-1}(v)$, then the first $N$ elements of $V$ equals to $z$, and the $DFT(V) = DFT(DFT^{-1}(v)) = v$