Perform a fast Fourier transform on the time series data. This decomposes your time series data into mean and frequency components and allows you to use variables for clustering that do not show heavy autocorrelation like many raw time series.
If time series is real-valued, discard the second half of the fast Fourier transform elements because they are redundant.
Separate the real and imaginary parts of each fast Fourier transform element.
Perform model-based clustering on the real and imaginary parts of each frequency element.
Plot the percentiles of the time series by cluster to examine their shape.
Alternately, you could omit the DC components of the fast Fourier transform to avoid your clusters being based on the mean and instead on the series defined by the Fourier transform, which represents the shape of the time series.
You will also want to calculate the amplitudes and phase angles from the fast Fourier transform so that you can explore the distribution of time series spectra within clusters. See this StackOverflow answer on how to do that for real-valued data.
You could also plot the percentiles of time series shape by cluster by computing the Fourier series from the amplitudes and phase angles (the resulting time series estimate will not perfectly match the original time series). You could also plot the percentiles of the raw time series data by cluster. Here is an example of such a plot, which came about from a harmonic analysis of NDVI data I just did today:
Finally, if your time series is not stationary (i.e., mean and variance shift over time), it may be more appropriate to use a wavelet transform rather than a Fourier transform. You would do so at the cost of information about frequencies while gaining information about location.