## 3 Results and Discussion  
### General Solution  
Here we follow the approach proposed earlier by [PERSON] [61] for Maxwell model and then generalize it for all other models. Consider a specific type of pressure gradient given as follows:  
\[\frac{\partial p}{\partial z}=P_{0}\mathrm{e}^{i\omega t} \tag{15}\]  
Introduce Fourier and inverse Fourier transforms of flow velocity and pressure:  
\[U_{z}(r,\omega)=\int_{-\infty}^{+\infty}u_{z}(r,t)\mathrm{e}^{-i\omega t}dt, \quad u_{z}(r,t)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}U_{z}(r,\omega)\mathrm{ e}^{i\omega t}d\omega \tag{16}\]  
\[P(r,\omega)=\int_{-\infty}^{+\infty}p(r,t)\mathrm{e}^{-i\omega t}dt,\quad p(r,t)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}P(r,\omega)\mathrm{e}^{i\omega t}d\omega \tag{17}\]  
Recall Fourier transform rule for fractional derivative:  
\[\mathcal{F}\{_{a}D_{t}^{\alpha}u_{z}\}=(i\omega)^{\alpha}U_{z}(r,\omega), \quad\mathcal{F}\{_{a}D_{t}^{\alpha-\beta}\frac{\partial p}{\partial z}\}=(i \omega)^{\alpha-\beta}\frac{\partial P(r,\omega)}{\partial z}, \tag{18}\]  
where appropriate integration limits are implied. In particular, lower terminal value is set to: \(a=-\infty\). Then the general algorithm to find flow velocity is the following: (1) express stress and its derivative from the constitutive eq-n; (2) plug stress and its derivative in momentum eq-n; (3) eliminate stress and its derivative and get modified momentum eq-n; (4) perform Fourier transform of modified momentum eq-n; (5) solve corresponding ODE for Fourier transform of velocity (\(U_{z}(r,\omega)\)); (6) change variables with \(\xi\) for simplicity; (7) perform inverse Fourier transform of \(U_{z}(r,\omega)\) to get \(u_{z}(r,t)\). Regardless of the model considered, \(z\)-component of the velocity reads as:  
\[u_{z}(r,t)=\frac{i}{\rho\omega}\frac{\partial p(z,t)}{\partial z}\Big{[}1- \frac{J_{0}(\xi r)}{J_{0}(\xi R)}\Big{]}, \tag{19}\]  
where \(J_{n}(\cdot)\) is the Bessel function of the first kind of order \(n\). The only difference that defines behavior of the system is hidden in \(\xi\). At the same time, the form of the solution reproduces a well-known classical result (see, for example Reference [73]). Thus it represents a natural extension of integer-order models using apparatus of FC. As long as velocity profile was defined, we could also get expression for the stress. The direct integration of (3) accounting for (15) results in:  
\[\sigma_{rz}=\frac{\partial p}{\partial z}\frac{J_{1}(\xi r)}{\xi J_{0}(\xi R )}, \tag{20}\]  
where again the individual properties of the model are hidden in parameter \(\xi\). Physical meaning of \(\xi\) is worth consideration. First, of all, \([z]=[1/\mathrm{m}]\), so that it can be considered as an \"effective\" wave number. Moreover, as shown below, \(\xi\) is related to complex compliance. Introducing \(\sigma_{0}=\frac{\partial p}{\partial z}R\), nondimensional stress \(\tilde{\sigma}=\sigma_{rz}/\sigma_{0}\) is given by:  
\[\tilde{\sigma}=\frac{J_{1}(\tilde{\xi}r)}{\xi RJ_{0}(\xi R)} \tag{21}\]For practical purposes and analysis simplification, it is more convenient to introduce the following nondimensional quantities and relationships:  
\[\bar{r}=\frac{r}{\sqrt{\
u\tau}},\quad\bar{a}=\frac{R}{\sqrt{\
u\tau}},\quad\bar {\omega}=\omega\tau,\quad\phi=\frac{\bar{r}}{\bar{a}},\quad u_{0}=\frac{1}{4 \rho\omega}\frac{\partial p}{\partial z},\quad r^{2}\xi^{2}=\bar{\omega}\bar{r} ^{2}, \tag{22}\]  
where \(\
u\) is a kinematic viscosity and all other quantities were defined earlier. Then velocity and stress profiles are given by:  
\[\bar{u}=\Big{|}\frac{u_{z}}{u_{0}}\Big{|}=\left(1-\frac{J_{0}(\sqrt{\bar{\omega }}\bar{a}\phi)}{J_{0}(\sqrt{\bar{\omega}}\bar{a})}\right),\quad\bar{\sigma}= \Big{|}\frac{\sigma_{rz}}{\sigma_{0}}\Big{|}=\frac{J_{1}(\sqrt{\bar{\omega}} \bar{a}\phi)}{\sqrt{\bar{\omega}}\bar{a}J_{0}(\sqrt{\bar{\omega}}\bar{a})} \tag{23}\]  
As velocity profiles and shear stresses have been defined, we can now proceed with considering specific fractional viscoelastic models and outlining their specifics.