## 3 Results and Discussion
### Fractional Zener Model  
Now examine the behavior of three-element fractional Zener model. The key model variable, \(\xi\), in this case becomes:  
\[\xi^{2}=\frac{\rho\omega^{2}}{E}\Big{[}\frac{1}{(E_{0}/E)(i\omega\tau)^{\gamma}+ \frac{(i\omega\tau)^{a+\beta}}{(i\omega\tau)^{a}+(i\omega\tau)^{\beta}}}\Big{]}, \tag{28}\]  
or alternatively in nondimensional form:  
\[\tilde{\omega}=\frac{\tilde{\omega}^{2}}{\psi(i\tilde{\omega})^{\gamma}+\frac {(i\tilde{\omega})^{a+\beta}}{(i\tilde{\omega})^{a}+(i\tilde{\omega})^{\beta}}}, \tag{29}\]  
where \(\psi=\frac{E_{0}}{E}\) is an elastic parameter. By definition elastic parameter \(\psi\) can be rewritten as:  
\[\psi=\frac{E_{3}\tau_{3}^{\gamma}}{E_{1}\tau_{1}^{a}}\tau^{\gamma-a} \tag{30}\]For classical Zener model (\(\alpha=1\), \(\beta=0\), and \(\gamma=0\)), \(\psi\) reduces to: \(\psi=E_{3}/E_{2}\), i.e., The ratio of [PERSON]'s moduli, thus justifying the notion introduced above. Similar to fractional Maxwell model, we impose additional conditions for fractional orders \(\alpha\), \(\beta\), and \(\gamma\). From two inequalities: \(\alpha-\beta\geq 0\) and \(\gamma+\alpha-\beta\geq 0\), we end up with the first one. The second one is satisfied automatically, provided the first one is and \(\gamma\) being non-negative. We first examine the dynamics of centerline velocity. Different from two-element models, where we considered the influence of three parameters (\(\alpha\), \(\beta\), and \(\bar{a}\)) on the system behavior, here we need to account for five parameters (\(\bar{a}\), \(\alpha\), \(\beta\), \(\gamma\) and \(\psi\)). Let us fix \(\alpha=0.5\) (as previously), \(\gamma=0.5\) and \(\psi=1\). Since \(\alpha-\beta\geq 0\), we only consider \(0\leq\beta\leq 0.5\). Corresponding centerline velocity profiles are given in Figure 7a,d,g. Centerline velocity changes its behavior from monotonically increasing to resonant with increasing \(\bar{a}\). Moreover, relative velocity amplitudes increase with increasing \(\bar{a}\). In addition to it, an interesting phenomenon is observed in Figure 7g. Relative velocity amplitudes tend to decrease with increasing fractional-order \(\beta\) up to a certain value (\(\bar{\omega}\sim 150\)). At this point, centerline velocities for all fractional orders become close, and then the trend inverses, namely velocity amplitudes start to increase with increasing fractional-order \(\beta\). Next, let us fix \(\alpha=\beta=0.5\), \(\psi=1\) varying \(\gamma\) and \(\bar{a}\). The corresponding plots for centerline velocity are shown in Figure 7b,e,h. Here we can observe switching from monotonically increasing to resonant behavior of the system with \(\bar{a}\) increasing and relative velocity amplitude decrease with increasing fractional order \(\gamma\). Moreover, for larger values of \(\bar{a}\) centerline velocity behavior changes from resonant to almost linearly increasing with increasing fractional order \(\gamma\). Finally, we fixed \(\alpha=\beta=\gamma=0.5\) and varied \(\bar{a}\) and \(\psi\) as shown in Figure 7c,f,i. As can be seen, the system behavior changes from monotonically increasing to resonant with increasing \(\bar{a}\). Moreover, resonant curves become less dispersive, with increasing \(\bar{a}\) and more dispersive with increasing \(\psi\). Relative amplitudes of centerline velocity also decrease with increasing \(\psi\). The dynamics of shear stresses are given in Figure 8. The first column in Figure 8 (Figure 8a,d,g) corresponds to \(\alpha=\gamma=0.5\), \(\psi=1\) and varying \(\tilde{a}\) and \(\beta\). In Figure 8g there exist two values of \(\tilde{\omega}\) where shear stresses become very close for different values of fractional order \(\beta\). Upon reaching the one corresponding to lower value of \(\tilde{\omega}\), the trend reverses with relative stress amplitude decreasing for increasing fractional order, \(\beta\). For large values of \(\tilde{\omega}\), shear stress becomes almost independent of the fractional order, \(\beta\).