Multiple seasonality with ARIMA?

I know that ARIMA can't detect multiple seasonality, but it is possible to use fourier functions to add a second seasonality.

I need to forecast gas consumption composed by a daily, weekly (week days-weekend), yearly seasonality. Does it make sense to apply three times the STL decomposition by LOESS? The reason is that I applied the fourier method and I have bad results but I don't know if it is only because I applied it wrong.

I'm interested in the theoretical explanation, but here you find also the code:

ARIMA + 2 STL:

b - ts(drop(coredata(dat.ts)), deltat=1/12/30/24, start=1)
fit - stl(b, s.window="periodic")
b - seasadj(fit)
dat.ts - xts(b, index(dat.ts))

# The weekdays are extracted
dat.weekdays - dat.ts[.indexwday(dat.ts) %in% 1:5]
dat.weekdaysTS - ts(drop(coredata(dat.weekdays)), frequency=24, start=1)
fit - stl(dat.weekdaysTS, s.window="periodic")
dat.weekdaysTS - seasadj(fit)

arima - Arima(dat.weekdaysTS, order=c(3,0,5))

With fourier:

dat.weekdays - dat.ts[.indexwday(dat.ts) %in% 1:5]
dat.weekdaysTS - ts(drop(coredata(dat.weekdays)), frequency=24, start=1)
z - fourier(ts(dat.weekdaysTS, frequency=365.25), K=5)
arima - Arima(dat.weekdaysTS, order=c(3,0,5),xreg=z)