Inverse DFT in terms of original signals for product of two DFTs
Again in this problem you can make use of the duality property to derive the output result without much involvement in the DFT sum formula. Before beginning, however, it's necesary to remember the time-frequency relations of two basic operations, namely, the expansion and compression blocks.

The EXPANSION block is the following:

$$ x_N[n] \longrightarrow \boxed{ \uparrow 2 } \longrightarrow y_{2N}[n] $$

shown operating on a time-domain signal $x_N[n]$ of length $N$ and producing the output $y_{2N}[n]$of length $2N$ given by the relation :

$$y_{2N}[n] = \begin{cases}{ x_N[n/2] ~~,~~~n=2m ~,~~~~~~ m=0,1,..,N-1 \\ ~~~~~0~~~~~~~~,~~~ n = 2m+1 ~, m=0,1,..,N-1} \end{cases} \tag{1}$$

The frequency domain relation between $N$-point DFT $X_n[k]$ of $x_N[n]$ and $2N$ point DFT $Y_{2N}[k]$ of $y_{2N}[n]$ is:

$$ Y_{2N}[k] = X_N[ (k)_N] ~~~,~~~ k = 0,1,2,...,2N-1 \tag{2}$$

Verbally stated, $Y_{2N}[k]$ is a cascade (mere repetition) of $X_N[k]$ and can be shown with the following matrix symbolism, each X,Y are (say) row vectors :

$$ Y = [ ~~~-~ X ~ -~~ | ~~~-~  X ~~~-~ ] \tag{3}$$  

The COMPRESSION block is the following:

$$ x_N[n] \longrightarrow \boxed{ \downarrow 2 } \longrightarrow y_{N/2}[n] $$

shown operating on a time-domain signal $x_N[n]$ of length $N$ and producing the output $y_{N/2}[n]$of length $N/2$ given by the relation :

$$y_{2N}[n] = x_N[2n] ~~~, ~~~ n = 0,1,..,N/2-1 \tag{4}$$

The frequency domain relation between $N$-point DFT $X_N[k]$ of $x_N[n]$ and $N/2$ point DFT $Y_{N/2}[k]$ of $y_{N/2}[n]$ is:

$$ Y_{N/2}[k] = 0.5( X_N[ k ] + X_N[k + N/2]) ~~~,~~~ k = 0,1,2,...,N/2-1 \tag{5}$$

Verbally stated, $Y_{N/2}[k]$ is obtained from $X_N[k]$ by folding: adding second half of $X_N[k]$ onto its first half. 

$$ Y = [ ~-~ {}_{N/2}X ~-~] + [~-~ X_{N/2} ~-~ ] \tag{6}$$  

Then the DUALITY principle states that, if frequency domain sequences are processed by those expansion and compression blocks, then their effects in the time-domain seuqences will be the same as above that happened on the frequency-domain sequences of time-domain expanion and compression. Note that, when different sample rates (or number of DFT points) are used in the input and output of the processing blocks, which is the case for the expander and compressor examples, then the duality relation is scaled by the ratio of sampling rates which would cancel the 0.5 at the Eq(5) in our case.

By using these concepts, we can describe the time-frquency domain relation of the problem you have provided. 

Now observe that your output signal $$Z[k]=Y(e^{j\omega})|_{\omega = \frac{2\pi}{N} 2k} = [Y(0) , Y(2), ... , Y(30) | Y(1) , Y(3) , ..., Y(29) ]$$ can be given by the following block diagram (with $N=31$):

$$ Y_N[k] \rightarrow \boxed{ Y_N -|- Y_N } \rightarrow W_{2N}[k] \rightarrow \boxed{ \downarrow 2} \rightarrow Z_N[k] $$  

Notice that the first blcok is given in (3) and related to the expander operation, where as the second block is clearly a compressor. Lets denote the time domain signals by $y_N[n]$, $w_{2N}[n]$ and $z_N[n]$, and proceed from left to right as follows:

First, given the relation beween the DFT sequences $Y_N[k]$ and $W_{2N}[k]$ then their time-domain originators is the expander operation and hence we have from Eq(1):

$$w_{2N}[n] = \begin{cases}{ y_N[n/2] ~~,~~~n=2m ~,~~~~~~ m=0,1,..,N-1 \\ ~~~~~0~~~~~~~~,~~~ n = 2m+1 ~, m=0,1,..,N-1} \end{cases} \tag{7}$$

And secondly, given the compressor block applied on the frequency signal $W_{2N}[k]$, we conclude its effect, by using duality, on the time domain signal $w_{2N}[n]$ as: 

$$ z_{N}[n] = w_{2N}[n] + w_{2N}[n+N/2] ~~~,~~~n=0,1,...,N-1 \tag{8}$$ namely the folding operation on $w_2[n]$. 

Combining Eq(8) and Eq(9) we conclude the relation between N-point sequences $z[n]$ and $y[n]$ as:

$$\boxed{z[n] = \begin{cases}{ ~~~y_N[n/2] ~~~,~~~n=2m ~,~~~ m=0,1,..,N-1 \\ y_N[n/2 + N/2],~~~ n = 2m+1 } \end{cases} }$$

Eventhough the description seems really long and somewhat abstract, the operation is really neat and quick as the following Matlab/OCtave code demonstrates  it:

N = 31;
y = 1:N;   % just some signal

Y = fft(y,N);    % obtain YN[k]
W = [Y, Y];      % W = [ Y | Y] , w[n] = y[n/2]: is expanded verison of y[n]
Z = W(1:2:end);  % Z[k] = W[2k]; compress W[k] to get Z[k]

z = ifft( Z, N); % obtain z[n] from Z[k]

% Part-II: Use the time-domain relation as described in the answer to obtain z[n]
we(1:2:2*N) = y ;           % we[n] = y[n/2]: expanded version of y[n]
we = [we , 0];              % Append (2-1) zeros to the end of we if N is odd.
z2 = we(1:N) + we(N+1:2*N); % folding we[n] to obtain we[n]

% Display:
[z' , z2' ]