1. Introduction
The application of fractional derivatives in time in rheological constitutive equations has a long history [
1]. One of the advantages of fractional time-derivatives in linear viscoelasticity is that with their help the damping behavior of viscoelastic media is described by using less parameters, compared with integer-order models. For a review of the main aspects of wave propagation in linear homogeneous viscoelastic media and the simplest and most used fractional constitutive models, we refer to [
2].
In linear viscoelasticity, the rheological properties of a viscoelastic medium are described through a linear constitutive relation between stress
and strain
. Following [
2], we restrict our considerations to the uniaxial case, in which
and
, and consider systems quiescent for all times prior to some starting time
.
The fractional Zener constitutive law is extensively used as a model of solid-like viscoelastic behavior [
2,
3]. The fractional Zener stress–strain relation reads
where
denotes the Riemann–Liouville fractional derivative [
4] of order
and
. Constitutive Equation (
1) is introduced for the first time in 1971 by Caputo and Mainardi [
5,
6]. Another early extensive study of model (
1) can be found in [
7]. A necessary and sufficient condition for physical acceptability of model (
1) is the restriction
. Under this condition, the corresponding relaxation modulus is a completely monotone function, see, e.g., [
2].
Recall that a functions
is said to be completely monotone if it is real-valued and infinitely differentiable on
and
Completely monotonic relaxation moduli are commonly encountered in polymer physics and other branches of rheology and mechanics [
2,
8,
9]. For a discussion on the importance of the complete monotonicity property of the relaxation moduli in linear viscoelasticity, we refer to [
9].
Wave equations in viscoelastic media described by the fractional Zener model are studied in many works, see, e.g., [
3,
10] and the references cited there. Various types of distributed-order generalizations of model (
1) are proposed and studied along with the corresponding wave equations: see, e.g., [
11,
12,
13,
14,
15,
16]. The complete monotonicity of the relaxation modulus for Zener-type constitutive equations with multiple time-derivatives and distributed-order time-derivatives with power-type weight function is established in [
17] under appropriate constraints on the parameters. The proofs are based on a representation of the relaxation modulus as a Laplace transform of a non-negative spectral function. For other applications of this method, we refer to [
2,
18].
Recently, equations with general fractional differential operators have gained increasing interest, see the pioneering work [
19], the review [
20], as well as [
21,
22,
23,
24,
25], to mention only few of many publications from the last years.
In the article [
26], published very recently, a generalization of model (
1) in the context of general fractional calculus is proposed and studied. In this work, the fractional derivative
in the constitutive Equation (
1) is replaced by a general convolutional derivative. For two specific cases of the kernel of this derivative, it is proven that the condition
is sufficient for the dissipation inequality to hold and thus, for the thermodynamic compatibility of the corresponding viscoelastic models.
Motivated by the work [
26], in the present paper, we study the following generalized fractional Zener model
Here,
and
are given parameters, and
is a generalized convolutional derivative in the Riemann–Liouville sense, defined by
where
k is a locally integrable memory kernel,
. For the kernel
, we assume in addition that its Laplace transform
exists for all
and
where
denotes the class of Stieltjes functions. The precise definition of this class is given in the next section, see (
6). The basic particular examples of operators
, defined in this way, are the first-order derivative (in this case,
), the Riemann–Liouville fractional derivative of order
(corresponding to
), as well as linear combinations with positive coefficients of such derivatives.
Let us note that assumptions (
5) are slightly weaker than those required in the definition of the general fractional derivative, introduced in [
19] and studied in detail in [
20,
22], where additional assumptions on the limiting behavior of
are imposed. We also point out that the assumption
, which is essential also in the definition of the general fractional derivative in [
19], allows the use of the convenient Bernstein functions technique developed in [
27].
In this work, based on the Bernstein functions technique, we establish the complete monotonicity of the relaxation modulus for model (
3), provided the constraint
is satisfied. This result is stronger than the thermodynamic compatibility results proven in [
26], since the dissipation inequality is satisfied for systems with completely monotone relaxation moduli. The complete monotonicity property implies also the subordination principle for the generalized fractional Zener wave equation. The propagation function is obtained as the solution of an equation governing wave propagation in a viscoelastic half-space with the generalized fractional Zener model (
3) and its properties are studied. The relation between the propagation function and the subordination kernel is established. Some examples of application of subordination principle are given. Numerical results for specific cases of generalized derivatives are presented and discussed.
The rest of the paper is organized as follows.
Section 2 contains preliminaries and properties of Bernstein functions and related classes of functions, which are necessary for the study. In
Section 3, we establish the complete monotonicity property of the relaxation modulus and discuss the relation of this property with thermodynamic compatibility conditions. Propagation of waves in generalized fractional Zener medium is discussed in
Section 4. A subordination principle is formulated in
Section 5 and applied to two basic problems. Numerical examples for some specific kernels are discussed in
Section 6. The last
Section 7 contains concluding remarks.
2. Preliminaries
This section contains preliminaries on the following functions classes: completely monotone functions (
), Stieltjes functions (
), Bernstein functions (
) and complete Bernstein function (
). Their definitions and properties, necessary for our proofs, are given in the terminology of [
27].
We use the following notations for the Laplace transform of a function
The Laplace transform of a function with respect to the time variable t, where x is considered as a parameter, is denoted by or .
We proceed with the definitions of the four special classes of functions. The definition of the class
of completely monotone functions is already given in the
Section 1, see (
2). The characterization of the class
is given by the Bernstein’s theorem: a function is completely monotone if and only if it can be represented as the Laplace transform of a non-negative (generalized) function.
The class of Bernstein functions consists of all non-negative functions defined and infinitely differentiable on , such that
The class of Stieltjes functions (
) consists of all functions
defined on
which admit the representation (see [
19])
where
,
, and the Laplace transform of
f exists for any
. The constant
can be determined from the identity (see, e.g., [
28], Theorem 2.6)
.
The class of complete Bernstein functions consists of all function defined on , such that
The simplest examples of functions from the classes
and
are given as:
The following properties are satisfied for :
- (P1)
The set is a convex cone: for all and . Moreover, for all . The set is closed under pointwise limits;
- (P2)
The sets , and are convex cones. They are closed under pointwise limits;
- (P3)
Any Stieltjes function is completely monotone. Any complete Bernstein function is a Bernstein function;
- (P4)
Let . Then, the function is completely monotone;
- (P5)
, where ∘ denotes composition of functions from the corresponding sets;
- (P6)
Every function , such that , is the Laplace transform of a function , which is locally integrable on and completely monotone;
- (P7)
if and only if ;
- (P8)
Let . Then, if and only if ;
- (P9)
Let
. Assume
are such that
. Then,
- (P10)
Any function
from the sets
or
admits an analytic extension to the complex plane cut along the negative real axis
, which satisfies
Moreover,
, where * denotes the complex conjugate, and
For details on the proofs, see [
27,
28], Theorem 2.6.
The Mittag–Leffler function is defined as follows:
where
and
. If
and
, then the function of Mittag–Leffler type
is completely monotone for
. The Laplace transform of this function is given by the identity
3. Completely Monotone Relaxation Modulus
Consider the generalized Zener constitutive Equation (
3). We assume that
and
for all
and
. Moreover, we suppose that
.
Let us note that the assumption
is equivalent to the following representation of the kernel
where
,
denotes the Dirac delta function and
is a completely monotone function.
The relaxation modulus
in a linear viscoelastic model is defined as the stress response to a unit step of strain. Specifically, for a quiescent system at
, the following equation is satisfied [
2]:
Applying formally the Laplace transform to (
3) and taking into account the Laplace transform identity
, we deduce
This, together with the identity
(which follows from (
11)), yields
We next discuss the properties of the relaxation modulus
, based on the Laplace transform pair (
12).
In a thermodynamically compatible model, the function
is non-negative and non-increasing. We prove next that these properties can be satisfied only if
. Indeed, since the rate of relaxation is non-increasing,
, it follows
for
. Therefore, taking into account identity (
12), we obtain for
Here, we have used that the initial value theorem for the Laplace transform implies
due to the assumption
as
. The inequality in (
13) implies
, since the denominator is a positive function. In this way, we have proved that condition
is necessary for viscoelastic model (
3) to be physically acceptable.
We prove next that the restriction
together with the assumptions (
5) on the kernel
implies the complete monotonicity of the relaxation modulus for viscoelastic model (
3).
Theorem 1. Assume that and that conditions (5) on the kernel are satisfied. Then, there exists a locally integrable and completely monotone function , , which Laplace transform is given by the identity (12). Proof. For the proof, we use the Bernstein functions technique. For
, it holds
Here, we use the assumption
, which is equivalent to
(see properties (P7) and (P8)). Since
, this implies
. By (P8), this is equivalent to
, which, by noting that
, yields
. Applying again (P7), it follows
. Moreover, (
12) and assumption
as
imply
. Therefore, according to (P6), there exists a locally integrable function
, which Laplace transform is given by (
12). □
By the uniqueness of Laplace transform, it follows that the function
from the above theorem coincides with the relaxation modulus of constitutive model (
3). In this way, we proved that the generalized Zener constitutive Equation (
3) defines a viscoelastic system with a completely monotonic relaxation modulus. It is known (see, e.g., [
8]) that systems with completely monotonic relaxation moduli obey the dissipation inequality, that is, for any
and
, the following holds:
Corollary 1. Under the conditions of Theorem 1, the dissipation inequality (16) is satisfied for model (3). Proof. For completeness, we give a direct proof of (
16). Taking into account (
11), the dissipation inequality (
16) is equivalent to
This is satisfied when
is a function of positive type, that is,
We will use the fact that every completely positive function is a function of positive type, see, e.g., [
29], Theorem 3.1. A function
is called completely positive if it satisfies the following conditions:
and there exist
and a non-negative and non-increasing function
, such that
It is known that any completely monotone function is completely positive, see, e.g., [
29]. For completeness, we will check directly that the relaxation modulus
for the generalized Zener model (
3) is a completely positive function. To this end, let us consider the function
Since
is a complete Bernstein function and
, property (P8) implies
. Moreover,
. Therefore, according to (P6), there exists a locally integrable function
, such that
. Then, identity (
19) is satisfied with
defined in this way and
. This can be easily checked by applying Laplace transform since
Therefore, the dissipation inequality (
16) is satisfied for model (
3). □
Let us note that conditions for thermodynamic compatibility are considered in different form in [
7]. In our notation, these conditions read
For completely monotone functions
, conditions (
20) are always fulfilled. This follows from the fact that in this case,
, which is equivalent to
(see (P7)), and then properties (P10) imply (
20).
Let us point out that since
, the limit
exists and is real (see [
27], Theorem 6.2), which implies
This limit confirms that the class of models (
3) govern solid-like behavior, see, e.g., [
2].
4. Generalized Fractional Zener Wave Equation
The one-dimensional equation governing wave propagation in a viscoelastic medium with the generalized Zener constitutive law is derived by combining constitutive Equation (
3), the one-dimensional equation of motion
and the kinematic equation
where
denotes the particle displacement,
is the body force,
, and
. Combining the three equations, we derive formally the following generalized wave equation of Zener type
It is worth noting that in the limiting case
, Equation (
24) is the classical wave equation, which is well studied. Therefore, in what follows, we consider only coefficients satisfying
Assuming the initial conditions
we rewrite the wave Equation (
24) as a Volterra integral equation as follows. Applying Laplace transform to (
24) and taking into account the initial conditions (
25), we obtain
Denote by
the function with Laplace transform
According to (
12), the function
is well defined and
Taking the inverse Laplace transform in (
26), we derive the following Volterra integral equation
where
Let us note that
implies
In what follows, we study Equation (
28) as a model for wave propagation in a viscoelastic medium with the generalized Zener constitutive law (
3). We consider different possibilities for the domain of the spatial operator.
4.1. Propagation Function
Let us consider first the basic problem of wave propagation of impact waves in a viscoelastic half-space with constitutive Equation (
3):
where the kernel
is defined in (
27) and
is the Heaviside unit step function. The solution
is often referred to as propagation function, since it represents the propagation in time of a disturbance at
.
It is worth noting that in the limiting case , corresponding to the classical wave equation, the propagation function is .
Let us also note that the initial conditions
follow from Equation (
29), taking into account representation (
27) for the kernel
.
By applying Laplace transform with respect to the temporal variable in (
29) and (
30), we obtain the following ODE in
x (with
s considered as a parameter):
where
The solution of problem (
32) is given by
The following property is important in the study of the behavior of function .
Proposition 1. If and satisfies conditions (5), then for the characteristic function , defined in (33), it holds Proof. Applying (
33) and (
12), we obtain the relation
Since
and
(here, we use
, and property (P8)), the function
is a product of two functions from the set
. Then, (
35) follows by applying property (P9) with
. □
According to (
35)
Then, property (P5) implies
with respect to
(and
considered as a parameter). Therefore, the function
is a product of two functions, completely monotone for
, which yields
. This implies by Bernstein’s theorem
. Moreover, the propagation function
is non-decreasing in
t and non-increasing in
x. This can be established in the same way as Theorem 2.2 in [
30] by proving that the Laplace transforms of the functions
and
are completely monotone functions. Indeed, since
, and
, by property (P4) also
. This, together with (
37), yields
as products of two completely monotone functions. Further, since
, it follows from (
34) and (
37)
4.2. Propagation Speed
Applying the general formula for the propagation speed
c of a disturbance in a viscoelastic medium [
2,
31], we obtain for the considered model (
3) with conditions (
5)
This can be seen by inverting the Laplace transform in representation (
34). Therefore, the propagation speed is finite for model (
3) and
for
.
From general theory (see, e.g., [
31], Section 5.4) in the case of finite propagation speed
c, a jump discontinuity of the propagation function
at the planar surface
exists only if
Moreover, this discontinuity is damped exponentially with exponent
, i.e.,
as
and
. In our specific case, it holds that
We find the leading terms in the last expression as
by taking into account the expansion
as
,
, for
and
. In this way, we deduce
Therefore, a jump at the wave front exists only when
, i.e., only when
in the kernel representation (
10). This is, for instance, the case with the composite derivative
, where
denotes the Riemann–Liouville fractional derivative of order
and
. For
, the classical Zener model is recovered. On the other hand, if
(
in the kernel representation (
10) and thus,
), the wave front is smooth. In this case, the coexistence of finite wave speed and absence of jump discontinuity at the wave front is observed.
In the next theorem, we summarize the most important findings in this section concerning the propagation function.
Theorem 2. Assume that and conditions (5) on the kernel are satisfied. Then, the solution of problem (29) and (30) is a non-negative function, non-increasing in and non-decreasing in . Moreover, for If , a jump discontinuity is observed at , otherwise, the function is smooth. Plots of the propagation functions
for specific kernels at different values of the parameters are given in
Figure 1 and
Figure 2. The observed behavior in the figures is in agreement with the established properties of the propagation function.
6. Numerical Results
In this section, we present numerical examples. The function and the subordination kernel are computed for three types of kernels . The numerical computations are based on integral representations of these functions, obtained as follows.
Integral representation for the propagation function
can be derived by inverting its Laplace transform (
34). Following the method of our previous work [
30], Theorem 2.5, which uses solely property (
35) and
as
(and not the specific form of
), we have the following integral representation for the propagation function:
where
In the case of characteristic function
, given in (
33), we obtain for the functions
and
where
and
, for
.
Further, the relation (
51) and the integral representation (
59) imply the following integral representation for
:
Here,
are the functions defined in (
60) and (
61).
Next, we list some examples of kernels
in the generalized derivative
, which satisfy conditions (
5). These conditions guarantee that the results obtained in this work for constitutive model (
3) with
are satisfied. The specific kernels are used in numerical experiments, the results of which are discussed in this section.
Since the numerical experiments are based on the integral representations (
59) for the propagation function
and (
62) for the subordination function
, for each kernel, we provide also the corresponding Laplace transforms and the expressions for the functions
and
, which appear in (
61).
For the sake of brevity, we use the notation
If , we set , , where is the Dirac delta function.
Example 1. Two-term composite derivative: , where and are Riemann–Liouville fractional derivatives with orders satisfying and . In this case,and Example 2. Generalized derivative with a singular Mittag–Leffler-type kernel: Example 3. Distributed-order fractional derivative with uniform distribution over the interval . In this case,and Another example of kernel, satisfying conditions (
5), is the following
where
and
. In fact, in [
26], the generalized fractional Zener model is studied for two types of kernels: (
65) and (
67). Generalizations of the Mittag–Leffler kernel (
65) can also be considered, such as the completely monotone multinomial Mittag–Leffler-type functions, see, e.g., [
35,
36].
In the present work, we treat numerically only the cases in Examples 1–3. All mentioned kernels admit decomposition (
10) and their Laplace transforms satisfy conditions (
5).
In
Figure 1 and
Figure 2, plots of propagation functions are given, while
Figure 3 presents plots of subordination functions for different kernels
. The numerical results provide a numerical verification of the analytical findings in this paper.
In all figures, except
Figure 2b, the parameters
a and
b are taken
, which corresponds to a constant propagation speed,
. Therefore, in these figures,
vanishes for
and
vanishes for
.
In
Figure 1a, the propagation function
is plotted for the kernel (
64) with
, which corresponds to the fractional Zener model (
1) when
. The function
is plotted as a function of
x for
and for different values of the parameter
. The limiting case
corresponds to the classical integer-order Zener model, where a jump discontinuity is observed at the wave front. The step function given with thinner line corresponds to
, recovering the classical elastic wave equation. For
, we observe, as expected, the coexistence of a finite wave speed and smooth wave front.
In
Figure 1b, the propagation function
is plotted for the composite kernel (
64) with
,
and
The function
is plotted as a function of
x for different times
t. Since
, it holds
, see (
40). Therefore, a jump discontinuity at the wave front is expected and observed. The thinner line is the graph of the function
and it marks the height of the jump.
Figure 2a corresponds to the Mittag–Leffler-type kernel (
65) with
and different values of
,
.
Figure 2b corresponds to the distributed-order kernel (
66),
a takes different values in the interval
,
,
. In
Figure 2, the propagation function
is plotted as a function of
x. Since
for the considered kernels, the propagation function is smooth at the wave front.
In
Figure 3, the subordination function
is plotted as a function of
for two types of kernels.
Figure 3a corresponds to kernel (
64) with
and
,
. Different values of
t are considered.
Figure 3b corresponds to the Mittag–Leffler-type kernel (
65) with
and different
,
;
. The plotted functions in
Figure 3 are probability densities as functions of
.