Fractional Burgers wave equation

TL;DR

This paper develops and solves fractional Burgers wave equations for viscoelastic media, revealing how different model parameters influence wave speed and propagation characteristics.

Contribution

It introduces two classes of thermodynamically consistent fractional Burgers models and derives their analytical solutions using integral transforms.

Findings

First class exhibits infinite wave propagation speed.

Second class exhibits finite wave propagation speed.

Spatial profiles show unexpected features in wave behavior.

Abstract

Thermodynamically consistent fractional Burgers constitutive models for viscoelastic media, divided into two classes according to model behavior in stress relaxation and creep tests near the initial time instant, are coupled with the equation of motion and strain forming the fractional Burgers wave equations. Cauchy problem is solved for both classes of Burgers models using integral transform method and analytical solution is obtained as a convolution of the solution kernels and initial data. The form of solution kernel is found to be dependent on model parameters, while its support properties implied infinite wave propagation speed for the first class and finite for the second class. Spatial profiles corresponding to the initial Dirac delta displacement with zero initial velocity display features which are not expected in wave propagation behavior.

Figures

Equations

$$\frac{\partial }{\partial x}\sigma (x,t)=\rho \frac{{\partial }^2}{\partial t^2}u(x,t),x\in \mathbb{R},t>0,$$

$$\epsilon (x,t)=\frac{\partial }{\partial x}u(x,t),x\in \mathbb{R},t>0;$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{\gamma }\right)\sigma \left(x,t\right)=\left(b_1{}_0{\mathrm{D}}_t^{\mu }+b_2{}_0{\mathrm{D}}_t^{\nu }\right)\epsilon \left(x,t\right),x\in \mathbb{R},t>0,$$

$${}_0{\mathrm{D}}_t^{\xi }y\left(t\right)=\frac{{\mathrm{d}}^{n+1}}{\mathrm{d}{t}^{n+1}}\left(\frac{t^{-\left(\xi -n\right)}}{\Gamma \left(1-\left(\xi -n\right)\right)}\ast y\left(t\right)\right),t>0,$$

$$u(x,0)=u_0(x),\frac{\partial }{\partial t}u(x,0)=v_0(x),\sigma (x,0)=0,\epsilon (x,0)=0,$$

$$\underset{x\to \pm \infty }{\lim }u(x,t)=0,\underset{x\to \pm \infty }{\lim }\sigma (x,t)=0,$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{\gamma }\right)\sigma \left(t\right)=\left(b_1{}_0{\mathrm{D}}_t^{\mu }+b_2{}_0{\mathrm{D}}_t^{\mu +\eta }\right)\epsilon \left(t\right)$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{\beta +\eta }\right)\sigma \left(t\right)=\left(b_1{}_0{\mathrm{D}}_t^{\beta }+b_2{}_0{\mathrm{D}}_t^{\beta +\eta }\right)\epsilon \left(t\right)$$

$$c=\sqrt{{\sigma }_{sr}^{\left(g\right)}}=\frac{1}{\sqrt{{\epsilon }_{cr}^{\left(g\right)}}}=\sqrt{\frac{b_2}{a_3}},$$

$$\sum _{i=1}^n{a}_i{}_0{\mathrm{D}}_t^{{\alpha }_i}\sigma (x,t)=\sum _{j=1}^m{b}_j{}_0{\mathrm{D}}_t^{{\beta }_j}\epsilon (x,t),a_i,b_j>0,{\alpha }_i,{\beta }_j\in \left(0,1\right),$$

$$\left(1+a{}_0{\mathrm{D}}_t^{\alpha }\right)\sigma (x,t)=E\left(1+b{}_0{\mathrm{D}}_t^{\alpha }\right)\epsilon (x,t),0<a\le b,\alpha \in \left(0,1\right),$$

$$\sum _{i=1}^n{a}_i{}_0{\mathrm{D}}_t^{{\alpha }_i}\sigma (x,t)=\sum _{i=1}^n{b}_i{}_0{\mathrm{D}}_t^{{\alpha }_i}\epsilon (x,t),0\le {\alpha }_1\le \dots \le {\alpha }_n<1,\frac{a_1}{b_1}\ge \dots \ge \frac{a_n}{b_n}\ge 0,$$

$$\frac{\partial }{\partial x}\sigma (x,t)=\frac{{\partial }^2}{\partial t^2}u(x,t),\epsilon (x,t)=\frac{\partial }{\partial x}u(x,t),$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{\gamma }\right)\sigma \left(x,t\right)=\left({}_0{\mathrm{D}}_t^{\mu }+b{}_0{\mathrm{D}}_t^{\mu +\eta }\right)\epsilon \left(x,t\right)$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{\beta +\eta }\right)\sigma \left(x,t\right)=\left({}_0{\mathrm{D}}_t^{\beta }+b{}_0{\mathrm{D}}_t^{\beta +\eta }\right)\epsilon \left(x,t\right),$$

$$u(x,0)=u_0(x),\frac{\partial }{\partial t}u(x,0)=v_0(x),\sigma (x,0)=0,\epsilon (x,0)=0,$$

$$\underset{x\to \pm \infty }{\lim }u(x,t)=0,\underset{x\to \pm \infty }{\lim }\sigma (x,t)=0,$$

$$\bar{x}=\frac{x}{\mathcal{U}},\bar{t}=\frac{t}{T^{\ast }},\bar{u}=\frac{u}{\mathcal{U}},{\bar{u}}_0=\frac{u_0}{\mathcal{U}},{\bar{v}}_0=\frac{T^{\ast }}{\mathcal{U}}{v}_0,\bar{\sigma }=\frac{\sigma }{{\sigma }^{\ast }},{\bar{a}}_1=\frac{a_1}{{\left(T^{\ast }\right)}^{\alpha }},{\bar{a}}_2=\frac{a_2}{{\left(T^{\ast }\right)}^{\beta }},$$

$$T^{\ast }={\left(\frac{\rho {\mathcal{U}}^2}{b_1}\right)}^{\frac{1}{2-\xi }},{\sigma }^{\ast }={\left(\frac{b_1^2}{{\left(\rho {\mathcal{U}}^2\right)}^{\xi }}\right)}^{\frac{1}{2-\xi }},{\bar{a}}_3=\frac{a_3}{{\left(T^{\ast }\right)}^{\zeta }},\bar{b}=\frac{b_2}{b_1{\left(T^{\ast }\right)}^{\eta }},$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{\gamma }\right)\sigma \left(t\right)=\left({}_0{\mathrm{D}}_t^{\mu }+b{}_0{\mathrm{D}}_t^{\mu +\eta }\right)\epsilon \left(t\right),$$

$$0\le \alpha \le \beta \le \gamma \le \mu \le 1,1\le \mu +\eta \le 1+\alpha ,b\le a_i\frac{\cos \frac{\left(\mu -\eta \right)\pi }{2}}{\left|\cos \frac{\left(\mu +\eta \right)\pi }{2}\right|},$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{2\alpha }\right)\sigma \left(t\right)=\left({}_0{\mathrm{D}}_t^{\mu }+b{}_0{\mathrm{D}}_t^{\mu +\alpha }\right)\epsilon \left(t\right),$$

$$\frac{1}{2}\le \alpha \le \beta \le \mu \le 1,\frac{a_3}{a_1}\frac{\left|\sin \frac{\left(\mu -2\alpha \right)\pi }{2}\right|}{\sin \frac{\mu \pi }{2}}\le b\le a_1\frac{\cos \frac{\left(\mu -\alpha \right)\pi }{2}}{\left|\cos \frac{\left(\mu +\alpha \right)\pi }{2}\right|};$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{\alpha +\beta }\right)\sigma \left(t\right)=\left({}_0{\mathrm{D}}_t^{\mu }+b{}_0{\mathrm{D}}_t^{\mu +\alpha }\right)\epsilon \left(t\right),$$

$$0\le \alpha \le \beta \le \mu \le 1,\alpha +\beta \ge 1,\frac{a_3}{a_2}\frac{\left|\sin \frac{\left(\mu -\beta -\alpha \right)\pi }{2}\right|}{\sin \frac{\left(\mu -\beta +\alpha \right)\pi }{2}}\le b\le a_1\frac{\cos \frac{\left(\mu -\alpha \right)\pi }{2}}{\left|\cos \frac{\left(\mu +\alpha \right)\pi }{2}\right|};$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{\alpha +\beta }\right)\sigma \left(t\right)=\left({}_0{\mathrm{D}}_t^{\mu }+b{}_0{\mathrm{D}}_t^{\mu +\beta }\right)\epsilon \left(t\right),$$

$$0\le \alpha \le \beta \le \mu \le 1,1-\alpha \le \beta \le 1-\left(\mu -\alpha \right),\frac{a_3}{a_1}\frac{\left|\sin \frac{\left(\mu -\alpha -\beta \right)\pi }{2}\right|}{\sin \frac{\left(\mu -\alpha +\beta \right)\pi }{2}}\le b\le a_2\frac{\cos \frac{\left(\mu -\beta \right)\pi }{2}}{\left|\cos \frac{\left(\mu +\beta \right)\pi }{2}\right|};$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{2\beta }\right)\sigma \left(t\right)=\left({}_0{\mathrm{D}}_t^{\mu }+b{}_0{\mathrm{D}}_t^{\mu +\beta }\right)\epsilon \left(t\right),$$

$$0\le \alpha \le \beta \le \mu \le 1,\frac{1}{2}\le \beta \le 1-\left(\mu -\alpha \right),\frac{a_3}{a_2}\frac{\left|\sin \frac{\left(\mu -2\beta \right)\pi }{2}\right|}{\sin \frac{\mu \pi }{2}}\le b\le a_2\frac{\cos \frac{\left(\mu -\beta \right)\pi }{2}}{\left|\cos \frac{\left(\mu +\beta \right)\pi }{2}\right|}.$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{\alpha +\beta }\right)\sigma \left(t\right)=\left({}_0{\mathrm{D}}_t^{\beta }+b{}_0{\mathrm{D}}_t^{\alpha +\beta }\right)\epsilon \left(t\right),$$

$$0\le \alpha \le \beta \le 1,\alpha +\beta \ge 1,\frac{a_3}{a_2}\le b\le a_1\frac{\cos \frac{\left(\beta -\alpha \right)\pi }{2}}{\left|\cos \frac{\left(\beta +\alpha \right)\pi }{2}\right|};$$

$$\left(1+a_1{}_0{\mathrm{D}}_t^{\alpha }+a_2{}_0{\mathrm{D}}_t^{\beta }+a_3{}_0{\mathrm{D}}_t^{2\beta }\right)\sigma \left(t\right)=\left({}_0{\mathrm{D}}_t^{\beta }+b{}_0{\mathrm{D}}_t^{2\beta }\right)\epsilon \left(t\right),$$

$$0\le \alpha \le \beta \le 1,\frac{1}{2}\le \beta \le \frac{1+\alpha }{2},\frac{a_3}{a_2}\le b\le a_2\frac{1}{\left|\cos \left(\beta \pi \right)\right|};$$

$$\left(1+{\bar{a}}_1{}_0{\mathrm{D}}_t^{\alpha }+{\bar{a}}_2{}_0{\mathrm{D}}_t^{2\alpha }\right)\sigma \left(t\right)=\left({}_0{\mathrm{D}}_t^{\alpha }+b{}_0{\mathrm{D}}_t^{2\alpha }\right)\epsilon \left(t\right),$$

$$\frac{1}{2}\le \alpha \le 1,\frac{{\bar{a}}_2}{{\bar{a}}_1}\le b\le {\bar{a}}_1\frac{1}{\left|\cos \left(\alpha \pi \right)\right|}.$$

$$\hat{f}(\xi )=\mathcal{F}\left[f\left(x\right)\right]\left(\xi \right)={\int }_{-\infty }^{\infty }f(x){\mathrm{e}}^{-\mathrm{i}\xi x}\mathrm{d}x,\xi \in \mathbb{R},$$

Full text

Key words: thermodynamically consistent fractional Burgers models, fractional Burgers wave equation, wave propagation speed

1 Introduction

Fractional Burgers wave equation is written as the system of equations consisting of: equation of motion corresponding to one-dimensional deformable body

[TABLE]

where $u$ and $\sigma$ are displacement and stress, while $\rho$ is constant material density; strain for small local deformations

[TABLE]

and constitutive equation represented by the fractional Burgers model

[TABLE]

having model parameters assumed as: $a_{1},a_{2},a_{3},b_{1},b_{2}>0,$ $\alpha,\beta,\mu\in\left[0,1\right],$ with $\alpha\leq\beta,$ and $\gamma,\nu\in\left[1,2\right],$ while the operator of Riemann-Liouville fractional derivative ${}_{0}\mathrm{D}_{t}^{\xi}$ of order $\xi\in\left[n,n+1\right],$ $n\in\mathbb{N}_{0},$ is defined by

[TABLE]

see [19], where $\ast$ denotes the convolution in time: $f\left(t\right)\ast_{t}g\left(t\right)=\int_{0}^{t}f\left(t^{\prime}\right)g\left(t-t^{\prime}\right)\mathrm{d}t^{\prime},$ $t>0.$

In order to solve the Cauchy problem on the real line $x\in\mathbb{R}$ and $t>0,$ the system of governing equations (1), (2), and (3) is subject to initial and boundary conditions:

[TABLE]

where $u_{0}$ is the initial displacement and $v_{0}$ is the initial velocity.

Considering the rheological scheme of the classical Burgers model, with the dash-pot element replaced by the Scott-Blair (fractional) element, the fractional Burgers model (3) is derived in [27]. Moreover, using the requirement of storage and loss modulus non-negativity, the analysis of thermodynamical consistency for fractional Burgers model (3), conducted in [27], yielded that the orders of fractional derivatives $\gamma,\nu\in\left[1,2\right]$ cannot be independent of the orders of fractional derivatives $\alpha,\beta,\mu\in\left[0,1\right],$ and this led to formulation of eight thermodynamically consistent fractional Burgers models, divided into two classes.

The first class contains five models, written as

[TABLE]

in an unified manner, such that the highest fractional differentiation order of strain is $\mu+\eta\in\left[1,2\right],$ with $\eta\in\left\{\alpha,\beta\right\},$ while the highest fractional differentiation order of stress is either $\gamma\in\left[0,1\right]$ in the case of Model I, with $0\leq\alpha\leq\beta\leq\gamma\leq\mu\leq 1$ and $\eta\in\left\{\alpha,\beta,\gamma\right\},$ or $\gamma\in\left[1,2\right]$ in the case of Models II - V, with $0\leq\alpha\leq\beta\leq\mu\leq 1$ and $\left(\eta,\gamma\right)\in\left\{\left(\alpha,2\alpha\right),\left(\alpha,\alpha+\beta\right),\left(\beta,\alpha+\beta\right),\left(\beta,2\beta\right)\right\}$. The fractional differentiation order of stress is less than the differentiation order of strain regardless on the interval $\left[0,1\right]$ or $\left[1,2\right].$

The second class contains three models, written as

[TABLE]

in an unified manner, such that $0\leq\alpha\leq\beta\leq 1$ and $\beta+\eta\in\left[1,2\right],$ with $\eta=\alpha,$ in the case of Model VI; $\eta=\beta$ in the case of Model VII; and $\alpha=\eta=\beta,$ $\bar{a}_{1}=a_{1}+a_{2},$ and $\bar{a}_{2}=a_{3}$ in the case of Model VIII. Considering the interval $\left[0,1\right],$ the highest fractional differentiation orders of stress and strain are equal, which also holds true for the orders from interval $\left[1,2\right].$

The responses in creep and stress relaxation tests for Models I - VIII are examined in [28]. Recall, creep compliance $\varepsilon_{cr}$ (relaxation modulus $\sigma_{sr}$) is the strain (stress) history function obtained as a response to the stress (strain) assumed as the Heaviside step function. It is found that models’ behavior near the initial time-instant is different for the first and the second model class: Models I - V have zero glass compliance, i.e., $\varepsilon_{cr}^{\left(g\right)}=\varepsilon_{cr}\left(0\right)=0$ and thus infinite glass modulus, i.e., $\sigma_{sr}^{\left(g\right)}=\sigma_{sr}\left(0\right)=\infty,$ while Models VI - VIII have non-zero glass compliance $\varepsilon_{cr}^{\left(g\right)}=\frac{a_{3}}{b_{2}}$ implying the non-zero glass modulus $\sigma_{sr}^{\left(g\right)}=\frac{b_{2}}{a_{3}}$ as well. On the other hand, the equilibrium compliance is infinite, i.e., $\varepsilon_{cr}^{\left(e\right)}=\lim_{t\rightarrow\infty}\varepsilon_{cr}\left(t\right)=\infty$, so that the equilibrium modulus is zero, i.e., $\sigma_{sr}^{\left(e\right)}=\lim_{t\rightarrow\infty}\sigma_{sr}\left(t\right)=0$ for both model classes and therefore all fractional Burgers models describe fluid-like materials. Note, if the equilibrium compliance is finite, then model would represent the solid-like material.

The implication, proved in the present work, is that fluid-like Burgers models belonging to the first class have infinite, while the ones belonging to the second class have finite wave propagation speed

[TABLE]

as in the case of thermodynamically consistent fractional models arising from the general fractional linear model

[TABLE]

obtained and analyzed in [2] for thermodynamical consistency and used in [22] as constitutive equations in wave propagation modeling. Namely, the results of [20, 21], where the wave propagation speed is found via the conic solution support, i.e., $\left|x\right|<ct,$ in the case of the fractional Zener model and its generalization, respectively given by

[TABLE]

are extended in [22], using the same argumentation as in the previous work, to all four classes of thermodynamically consistent linear fractional models and moreover to the power-type distributed-order model assuming that the orders of fractional differentiation do not exceed one. In particular, it is found that both solid-like and fluid-like materials can have either infinite or finite wave speed. Singularity propagation properties of the memory and non-local type fractional wave equations are investigated in [17, 18] using the tools of microlocal analysis, supporting the results obtained in [20].

Wave propagation phenomena in viscoelastic bodies, modeled by integer and fractional order models, including the question of wave speed and energy dissipation properties are analyzed in [8, 9]. The wavefront expansion of solution, due to Buchen and Mainardi, is introduced in [7] to be later used in [11, 12] when considering the wave equation in viscoelastic materials described by the Bessel as well as by the integer and fractional order Maxwell and Kelvin-Voigt models. The Bessel model for viscoelastic body is introduced in [13] and analyzed in [10]. Features of the wave propagation in viscoelastic media, like the asymptotic behavior of fundamental solution near the wavefront, dispersion, and attenuation is examined in [14, 15, 16]. Wave propagation speed, reinterpreted as the fundamental solution’s peak propagation speed is analyzed in [23, 24, 25]. Modeling viscoelastic materials using the fractional order models, as well as dispersion and attenuation effects described by the corresponding wave equations are reviewed in [26].

Fractional wave equations on bounded and semi-bounded domain are considered in [29, 30, 31] for different fractional models including the Zener, modified Zener, and modified Maxwell models, as well as in [4, 5, 6] in the case of power-type distributed-order model. Generalizations of the classical wave equations and corresponding problems are reviewed in [3, 32].

2 Fractional Burgers model in wave propagation

Fractional Burgers wave equation, as the dimensionless system of equations:

[TABLE]

and either

[TABLE]

for the first class of Burgers models, or

[TABLE]

for the second class of Burgers models, subject to initial and boundary conditions

[TABLE]

is obtained by introducing dimensionless quantities

[TABLE]

with $\xi=\mu$ and $\zeta=\gamma$ for the first class of Burgers models, $\xi=\beta$ and $\zeta=\beta+\eta$ for the second class, and $\mathcal{U}=\sup_{x\in\mathbb{R}}\left|u_{0}\left(x\right)\right|,$ into system of governing equations (1), (2) and either (6) or (7), subject to (4), (5), and by subsequent omittance of bars.

Models in dimensionless form, along with the corresponding thermodynamical restrictions, are listed below.

Model I:

[TABLE]

with $\left(\eta,i\right)\in\left\{\left(\alpha,1\right),\left(\beta,2\right),\left(\gamma,3\right)\right\};$

Model II:

[TABLE]

Model III:

[TABLE]

Model IV:

[TABLE]

Model V:

[TABLE]

Model VI:

[TABLE]

Model VII:

[TABLE]

Model VIII:

[TABLE]

Application of the Fourier transform with respect to the spatial coordinate

[TABLE]

and Laplace transform with respect to the time

[TABLE]

with initial (13) and boundary conditions (14) taken into account, transforms the system of governing equations (10) and either (11), or (12) into ($\xi\in\mathbb{R},$ $\mathop{\mathrm{R}e}s>0$)

[TABLE]

with either

[TABLE]

in the case of the first class of Burgers equation (11), or

[TABLE]

in the case of the second class of Burgers equation (11).

It is obtained that

[TABLE]

with

[TABLE]

once the system of equations (31), (32) is solved with respect to displacement $\widehat{\tilde{u}}$, implying the solution to the fractional Burgers equation (10) and either (11), or (12), subject to (13) and (14), in the form

[TABLE]

where $\ast_{x}$ denotes the convolution with respect to the spatial variable: $f\left(x\right)\ast_{x}g\left(x\right)=\int_{-\infty}^{\infty}f\left(x^{\prime}\right)g\left(x-x^{\prime}\right)\mathrm{d}x^{\prime},$ $x\in\mathbb{R},$ after inverting Fourier and Laplace transforms in (35).

In order to calculate the solution kernel $K,$ the inversion of the Fourier transform is performed in (36) using a well-known inversion formula

[TABLE]

implying

[TABLE]

provided that

[TABLE]

which holds for all Models I - VIII, as proved in Appendix A. Further, inverting the Laplace transformation in (39) by the definition

[TABLE]

where $\Gamma_{0}$ is the Bromwich path, the two forms of solution kernel $K$ are obtained in Appendix B depending on the number and position of branching points of function $\tilde{K},$ given by (39), originating from the zeros of function $\Phi_{\sigma},$ since $\Phi_{\varepsilon},$ except for $s=0,$ has no other zeros in the principal Riemann plane, with $\Phi_{\sigma}$ and $\Phi_{\varepsilon}$ given by either (33) or (34). There are three possible cases, since, as shown in [28], function $\Phi_{\sigma}$ can have no zeros, one negative real zero, or a pair of complex conjugated zeros having negative real part. However, the solution kernel has the same form in the first two cases, thus merged into Case 1 below, while the form of the solution kernel differs in the third case, thus being labeled as Case 2.

Case 1. If function $\tilde{K},$ except for $s=0,$ either has no branching points, or has a negative real branching point, then function $K$ is found as

[TABLE]

either having support in $\mathbb{R}\times\left[0,\infty\right)$ for the first class of fractional Burgers models, or having support in the conic domain $\left|x\right|<\sqrt{\frac{b}{a_{3}}}t,$ for the second class.

Case 2. If function $\tilde{K},$ except for $s=0,$ has a pair of complex conjugated branching points with negative real part: $s_{0}=\rho_{0}\mathrm{e}^{\mathrm{i}\varphi_{0}}$ and $\bar{s}_{0}=\rho_{0}\mathrm{e}^{-\mathrm{i}\varphi_{0}}$, then function $K$ is found as

[TABLE]

either having support in $\mathbb{R}\times\left[0,\infty\right)$ for the first class of fractional Burgers models, or having support in the conic domain $\left|x\right|<\sqrt{\frac{b}{a_{3}}}t,$ for the second class.

The solution support properties, in both cases of solution kernel, define the wave propagation speed: infinite if the support is $\mathbb{R}\times\left[0,\infty\right),$ obtained for the first class of Burgers models, and finite if the support is conic domain $\left|x\right|<\sqrt{\frac{a_{3}}{b}}t,$ obtained as

[TABLE]

for the second class of Burgers models. Since $\sigma_{sr}^{\left(g\right)}=\frac{a_{3}}{b},$ see [28, Eq. (57)], the wave propagation speed (44) is exactly the wave propagation speed (8) that is obtained in [22] for the constitutive models having fractional differentiation orders not exceeding one.

3 Numerical examples

Spatial profiles of the solution to the fractional Burgers wave equations, written as the system of equations (10) and either (11), or (12), subject to initial and boundary conditions (13) and (14), with the initial displacement being the Dirac delta distribution and initial velocity being zero, i.e., $u_{0}=\delta,$ and $v_{0}=0,$ implying that the solution is equal to the solution kernel $K,$ are depicted in Figures 1, 2, and 3 for Model V, representing the first class of fractional Burgers models and in Figures 4, 5, and 6 for Model VII, representing the second class. Recall, in the case of constitutive models belonging to the first class the wave propagation speed is infinite, while in the case of the second class the speed is finite and given by (8). Spatial profiles produced by using the analytical formula for solution kernel $K,$ given by either (42), or (43), are compared with the solution kernel numerically calculated by the fixed Talbot numerical Laplace inversion Mathematica function, developed by J. Abate and P. P. Valkó according to [1] and available at: http://library.wolfram.com/infocenter/MathSource/4738/. In each of the numerical examples good agreement between profiles obtained by these two methods is found.

Figures 1, 2, and 3 present spatial profiles for Model V in cases when function $\tilde{K},$ given by (39), except for $s=0$ does not have other branching points, has one negative real, and has a pair of complex conjugated branching points, respectively. Different number and position of the branching points is a consequence of the change of a single parameter $\beta.$ Apart from the main peak originating from the propagation of the initial Dirac delta displacement, there is a noticeable additional peak that is more prominent for small times and ceasing as time passes. As the parameter $\beta$ increases, the change of the nature (number and position) of the branching points from no branching points to a pair of complex conjugated ones, implies the growth of prominence of the additional peak. During the propagation, due to the energy dissipation, height of the main peak decreases, while the width of profile is increasing, while propagation itself is rather slow.

Wave propagation speed is finite for the second class of fractional Burgers models, and in Figures 4, 5, and 6, presenting spatial profiles for Model VII, it is underlined by denoting the ending points of solution support by circles. It is also noticeable that during the propagation, due to the energy dissipation, height of the peak decreases, while its width increases.

Figure 4 presents spatial profiles depending on the nature of the branching points, different than $s=0,$ of function $\tilde{K}$ given by (39) in three cases obtained as a consequence of changing parameter $\beta$: Figure 4a represents case when there are no other branching points, Figure 4b when there is one negative real branching point, and Figure 4c when there is a pair of complex conjugated branching points. For small times, the profile shapes are considerably different, while as time passes the profile shapes become alike. In all cases there are jumps at the ending points of solution support: in Figures 4a and 4b displacement jumps from a positive value to zero, while in Figure 4c displacement jumps from a negative value to zero.

When compared to the profiles from Figure 4a, where the displacement jumps to zero at the ending point of solution support, the displacements plotted in Figure 5, representing also the case when there are no other branching points than $s=0,$ tend smoothly to zero at the ending points of solution support. Profiles from Figure 5 are similar to the profiles obtained in [20, 21, 22] for fractional constitutive models used wave propagation modeling in viscoelastic dissipative media.

Figure 6 presents spatial profiles in another case of model parameters yielding existence of a pair of complex conjugated branching points (apart of $s=0$) which differ from the ones presented in Figure 4c, since it seems that peaks are situated at zero, while displacement seems to converge to infinity at the ending point of solution support.

4 Conclusion

Fractional Burgers wave equations, considered as a dimensionless system of: equation of motion and strain (10), coupled with the constitutive Burgers models either of the first class (11), or of the second class (12), are solved for the Cauchy initial value problem and their solutions as a response to the initial Dirac delta displacement with zero initial velocity are qualitatively analyzed through numerical examples. The method of Fourier, with respect to space, and Laplace transform with respect to time are used in order to obtain analytical solution as a convolution of the solution kernels and initial data. The form of the solution kernel proved to be dependant on model parameters, so that if parameters yield, except for $s=0,$ either no branching points, or one negative real branching point of the Laplace transform of solution kernel, then solution kernel takes the form (42), while if, except for $s=0,$ the Laplace transform of solution kernel has a pair of complex conjugated branching points, then solution kernel takes the form (43).

Arising from the solution support properties, in both cases of solution kernel, the infinite wave propagation speed is obtained for the first class of Burgers models and finite for the second class. Moreover, the obtained wave propagation speed is consistent with the one obtained for the wave equations involving fractional linear models with differentiation orders below one.

Qualitative analysis has shown the dissipative behavior for both classes of Burgers wave equations, as expected from thermodynamically consistent constitutive laws for viscoelastic body. However, spatial profile shapes differs for the different nature of the branching points. The features of spatial profiles include the jumps from finite value of displacement to zero at the ending points of solution support, as well as profiles that are not expected in wave propagation behavior, like occurrence of the additional peaks and peaks situated at zero.

Appendix A Justification for using the Fourier inversion formula

The solution kernel is obtained by the Fourier and Laplace transforms as (36), and in order to apply the Fourier transform inversion formula (38), the condition (40), i.e.,

[TABLE]

must be fulfilled.

Functions $\Phi_{\sigma}$ and $\Phi_{\varepsilon},$ given by (33) in the case of the first, or by (34) in the case of the second model class, are never zero for $\mathop{\mathrm{R}e}s>0.$ Namely, it is well-known that function $\Phi_{\varepsilon},$ except for $s=0,$ does not have other zeros in the principal Riemann branch $\arg s\in\left(-\pi,\pi\right),$ while for function $\Phi_{\sigma}$ it is proved in [28] that if it has zeros, then they lie in the left complex half-plane.

Therefore, it is left to prove that

[TABLE]

It is clear that if $s=\rho>0,$ then

[TABLE]

Further, by substituting $s=\rho\mathrm{e}^{\mathrm{i}\varphi},$ $\varphi\in\left(-\frac{\pi}{2},\frac{\pi}{2}\right),$ into (45) one obtains

[TABLE]

with

[TABLE]

that will for each fractional Burgers model prove to be strictly positive if $\varphi\in\left(0,\frac{\pi}{2}\right)$ implying that $\psi,$ given by (45) cannot be zero for $\mathop{\mathrm{R}e}s>0.$ Since $\mathop{\mathrm{I}m}\psi\left(\rho,-\varphi\right)=-\mathop{\mathrm{I}m}\psi\left(\rho,\varphi\right),$ note that $\mathop{\mathrm{I}m}\psi\left(\rho,\varphi\right)<0$ if $\varphi\in\left(-\frac{\pi}{2},0\right).$

Model I

is obtained for $\eta\in\left\{\alpha,\beta,\gamma\right\},$ so that function $f_{\rho},$ given by (46), reads

[TABLE]

The thermodynamical restrictions (16) imply the positivity of all terms in (50), yielding $f_{\rho}\left(\varphi\right)>0$ if $\varphi\in\left(0,\frac{\pi}{2}\right).$

Model II

is obtained for $\gamma=2\alpha$ and $\eta=\alpha,$ so that function $f_{\rho},$ given by (46), reads

[TABLE]

Consider function $g$ and its first derivative $g^{\prime}$:

[TABLE]

on the interval $\varphi\in\left(0,\frac{\pi}{2}\right).$ Let $0<\zeta<\xi<1.$ Since function $\frac{\tan x}{x}$ is monotonically increasing for $x\in\left(0,\frac{\pi}{2}\right),$ one has $g^{\prime}\left(\varphi\right)>0,$ $\varphi\in\left(0,\frac{\pi}{2}\right),$ implying that function $g$ is an increasing function on the same interval and therefore

[TABLE]

The thermodynamical restriction (18) yields $0<2\alpha-\mu<\mu<1,$ so that by setting $\zeta=2\alpha-\mu$ and $\xi=\mu$ in function $g$ given by (52), using (53) one has

[TABLE]

Therefore, again by (18), one has that $b-\frac{a_{3}}{a_{1}}\frac{\left|\sin\left(\left(\mu-2\alpha\right)\varphi\right)\right|}{\sin\left(\mu\varphi\right)}>0,$ which, along with the positivity of all other terms in (51), implies that $f_{\rho}\left(\varphi\right)>0$ if $\varphi\in\left(0,\frac{\pi}{2}\right).$

Model III

is obtained for $\gamma=\alpha+\beta$ and $\eta=\alpha,$ so that function $f_{\rho},$ given by (46), reads

[TABLE]

The thermodynamical restriction (20) yields $0<\alpha-\left(\mu-\beta\right)<\alpha+\left(\mu-\beta\right)<1,$ so that by setting $\zeta=\alpha-\left(\mu-\beta\right)$ and $\xi=\alpha+\left(\mu-\beta\right)$ in function $g$ given by (52), using (53) one has

[TABLE]

Therefore, again by (20), one has that $b-\frac{a_{3}}{a_{1}}\frac{\left|\sin\left(\left(\mu-\beta-\alpha\right)\varphi\right)\right|}{\sin\left(\left(\mu-\beta+\alpha\right)\varphi\right)}>0,$ which, along with the positivity of all other terms in (54), implies that $f_{\rho}\left(\varphi\right)>0$ if $\varphi\in\left(0,\frac{\pi}{2}\right).$

Model IV

is obtained for $\gamma=\alpha+\beta$ and $\eta=\beta,$ so that function $f_{\rho},$ given by (46), reads

[TABLE]

The thermodynamical restriction (22) yields $0<\beta-\left(\mu-\alpha\right)<\beta+\left(\mu-\alpha\right)<1,$ so that by setting $\zeta=\beta-\left(\mu-\alpha\right)$ and $\xi=\beta+\left(\mu-\alpha\right)$ in function $g$ given by (52), using (53) one has

[TABLE]

Therefore, again by (22), one has that $b-\frac{a_{3}}{a_{1}}\frac{\left|\sin\left(\left(\mu-\alpha-\beta\right)\varphi\right)\right|}{\sin\left(\left(\mu-\alpha+\beta\right)\varphi\right)}>0,$ which, along with the positivity of all other terms in (55), implies that $f_{\rho}\left(\varphi\right)>0$ if $\varphi\in\left(0,\frac{\pi}{2}\right).$

Model V

is obtained for $\gamma=2\beta$ and $\eta=\beta,$ so that function $f_{\rho},$ given by (46), reads

[TABLE]

The thermodynamical restriction (24) yields $0<2\beta-\mu<\mu<1,$ so that by setting $\zeta=2\beta-\mu$ and $\xi=\mu$ in function $g$ given by (52), using (53) one has

[TABLE]

Therefore, again by (24), one has that $b-\frac{a_{3}}{a_{2}}\frac{\left|\sin\left(\left(\mu-2\beta\right)\varphi\right)\right|}{\sin\left(\mu\varphi\right)}>0,$ which, along with the positivity of all other terms in (56), implies that $f_{\rho}\left(\varphi\right)>0$ if $\varphi\in\left(0,\frac{\pi}{2}\right).$

Model VI

is obtained for $\gamma=\alpha+\beta,$ $\mu=\beta,$ and $\eta=\alpha,$ so that function $f_{\rho},$ given by (46), reads

[TABLE]

The thermodynamical restriction (26) yields $b-\frac{a_{3}}{a_{2}}>0,$ which, along with the positivity of all other terms in (57), implies that $f_{\rho}\left(\varphi\right)>0$ if $\varphi\in\left(0,\frac{\pi}{2}\right).$

Model VII

is obtained for $\gamma=2\beta$ and $\mu=\eta=\beta,$ so that function $f_{\rho},$ given by (46), reads

[TABLE]

The thermodynamical restriction (28) yields $b-\frac{a_{3}}{a_{2}}>0,$ which, along with the positivity of all other terms in (58), implies that $f_{\rho}\left(\varphi\right)>0$ if $\varphi\in\left(0,\frac{\pi}{2}\right).$

Model VIII

is obtained for $\gamma=2\alpha,$ $\beta=\mu=\eta=\alpha,$ $a_{1}+a_{2}=\bar{a}_{1},$ and $a_{3}=\bar{a}_{2},$ so that function $f_{\rho},$ given by (46), reads

[TABLE]

The thermodynamical restriction (30) yields $b-\frac{\bar{a}_{2}}{\bar{a}_{1}}>0,$ which, along with the positivity of all other terms in (59), implies that $f_{\rho}\left(\varphi\right)>0$ if $\varphi\in\left(0,\frac{\pi}{2}\right).$

Appendix B Calculation of the solution kernel

In order to obtain the solution kernels, given by (42) and (43), the inverse Laplace transform (41) will be calculated using the Cauchy integral formula

[TABLE]

where $\Gamma$ is a closed curve containing the Bromwich path $\Gamma_{0}$ from the Laplace inversion formula (41) and chosen differently depending on the number and position of the branching points of function $\tilde{K},$ given by (39).

Branching points of function $\tilde{K}$ are points in which the function under the square root is zero, i.e., in (39) either $\Phi_{\sigma}(s)=0$ or $\Phi_{\varepsilon}(s)=0,$ $s\in\mathbb{C},$ with $\Phi_{\sigma}$ and $\Phi_{\varepsilon}$ given by (33) in the case of the first or by (34) in the case of the second model class. Function $\Phi_{\varepsilon},$ except for $s=0,$ does not have other zeros in the principal Riemann plane $\arg s\in\left(-\pi,\pi\right),$ since

[TABLE]

as proved in [22]. Zeros of function

[TABLE]

with $a_{1},a_{2},a_{3}>0,$ $\alpha,\beta\in\left(0,1\right),$ $\gamma\in\left(0,2\right),$ and $\alpha<\beta<\gamma,$ are analyzed in [28], where it is found that if $\gamma\in\left(0,1\right),$ then function $\Phi_{\sigma}$ has no zeros in the complex plane, which is valid for Model I, while if $\gamma\in\left(1,2\right),$ then the number and position of zeros of function $\Phi_{\sigma}$ is as follows:

[TABLE]

where

[TABLE]

with $\rho^{\ast}$ determined from $\mathop{\mathrm{I}m}\Phi_{\sigma}\left(\rho^{\ast}\right)=0,$ i.e.,

[TABLE]

which is valid for Models II - VII. In the case of Model VIII, zeros of function

[TABLE]

are as follows:

[TABLE]

with $\rho^{\ast}$ determined by

[TABLE]

Note that the branching point $s=0$ is due to the differentiation of fractional order and that function $\tilde{K}$ does not have any singularities other than branching points, justifying the use of the Cauchy integral formula.

B.1 Case 1.

Function $\tilde{K},$ except for $s=0,$ has no other

branching points

If function $\tilde{K}$ (39), except for $s=0,$ has no other branching points, then the contour $\Gamma$ appearing in the Cauchy integral formula (60) is chosen as in Figure 7 and parametrized as in Table 1.

[TABLE]

The integrals along contours $\Gamma_{3},$ $\Gamma_{5},$ and $\Gamma_{0},$ calculated as

[TABLE]

yield the solution kernel $K$ in the form (42) when used in the Cauchy integral formula (60), since the integrals along all other contours will prove to be zero.

The following estimates will be used. According to (33), respectively (34), after the substitution $s=\rho\mathrm{e}^{\mathrm{i}\varphi}$ is made, it is obtained that

[TABLE]

and therefore

[TABLE]

The integral along contour $\Gamma_{1}$ reads

[TABLE]

and since $p+\mathrm{i}R\sim R\,\mathrm{e}^{\mathrm{i}\frac{\pi}{2}},$ as $R\rightarrow\infty,$ one has

[TABLE]

The use of (68) and (71) in (72), due to $0<\frac{\mu+\eta-\gamma}{2}<1,$ yields

[TABLE]

for the first model class and choosing $p_{0}=0$

[TABLE]

for the second model class. Similar argumentation is valid for the integral along $\Gamma_{7}$.

The integral along contour $\Gamma_{2}$ takes the form

[TABLE]

so that

[TABLE]

Using (68) and (71) in (73), due to $0<\frac{\mu+\eta-\gamma}{2}<1$ and $\cos\varphi<0$ for $\varphi\in\left[\frac{\pi}{2},\pi\right],$ yields

[TABLE]

for $\left(x,t\right)\in\mathbb{R}\times\left[0,\infty\right),$ in the case of the first model class and

[TABLE]

for the second model class. Similar argumentation is valid for the integral along $\Gamma_{6}$.

The integral along contour $\Gamma_{4}$:

[TABLE]

tends to zero when $r\rightarrow 0$, since

[TABLE]

due to $\beta,\mu<1$ and

[TABLE]

Function $\tilde{K},$ except for $s=0,$ has a negative real

branching point

If function $\tilde{K}$ (39), except for $s=0,$ has a negative real branching point $-\rho^{\ast},$ determined by (61) or (62), then the contour $\Gamma$ appearing in the Cauchy integral formula (60) is chosen as in Figure 8 and parametrized as in Table 2.

[TABLE]

The integrals along contours $\Gamma_{3a}\cup\Gamma_{3b},$ $\Gamma_{5a}\cup\Gamma_{5b},$ and $\Gamma_{0},$ when $r\rightarrow 0$ and $R\rightarrow\infty,$ are the same integrals as (63), (64), and (65), thus yielding the solution kernel $K$ in the form (42) when used in the Cauchy integral formula (60), since the integrals along contours $\Gamma_{1},$ $\Gamma_{2},$ $\Gamma_{4},$ $\Gamma_{6},$ and $\Gamma_{7}$ already proved to be zero, while the integrals along $\Gamma_{8}$ and $\Gamma_{9}$ will prove to be zero.

Namely, the integral along $\Gamma_{8}$ reads

[TABLE]

so that

[TABLE]

since

[TABLE]

because of $-\rho^{\ast}$ being zero of function $\Phi_{\sigma}.$ Similar argumentation is valid for the integral along $\Gamma_{9}$.

B.2 Case 2.

Function $\tilde{K},$ except for $s=0,$ has a pair of

complex conjugated branching points

If function $\tilde{K},$ except for $s=0,$ has a pair of complex conjugated branching points with negative real part: $s_{0}=\rho_{0}\mathrm{e}^{\mathrm{i}\varphi_{0}}$ and $\bar{s}_{0}=\rho_{0}\mathrm{e}^{-\mathrm{i}\varphi_{0}}$, then the contour $\Gamma$ appearing in the Cauchy integral formula (60) is chosen as in Figure 9 and parametrized as in Table 3.

[TABLE]

The solution kernel $K$ in the form (43) is obtained when the integrals along contours $\Gamma_{3a}\cup\Gamma_{3b},$ $\Gamma_{5a}\cup\Gamma_{5b},$ and $\Gamma_{0},$ calculated as

[TABLE]

are used in the Cauchy integral formula (60), since the integrals along contours $\Gamma_{1},$ $\Gamma_{2},$ $\Gamma_{4},$ $\Gamma_{6},$ and $\Gamma_{7}$ already proved to be zero, while the integrals along $\Gamma_{8}$ and $\Gamma_{9}$ will prove to be zero.

The integral along $\Gamma_{8}$ reads

[TABLE]

so that

[TABLE]

since

[TABLE]

because of $s_{0}$ being zero of function $\Phi_{\sigma}.$ Similar argumentation is valid for the integral along $\Gamma_{9}$.

Acknowledgment

This work is supported by the Serbian Ministry of Education, Science and Technological Development under grants $174005$ and $174024$, by the Provincial Secretariat for Higher Education and Scientific Research under grant $142-451-2384/2018$, as well as by FWO Odysseus project of Michael Ruzhansky.