This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 License.
Introduction
This paper offers a thorough introduction to mathematical tools to describe wave propagation in solids modeled with a wide collection of viscoelastic laws. Before we even attempt a general description of the models we will be addressing, let us emphasize what our goals are and what has not been tackled in the present paper. We aim for a unified mathematical description of a wide collection of known viscoelastic models, including basic well-posedness results. The models will include all classical viscoelastic wave models, fractional versions thereof, and couplings of different models in different subregions. The techniques that we will employ in the first part of the article (Sections 3 through 7) are those of Laplace transforms, understanding the influence of volume forces, normal stresses, and given displacements, as well as the strain-stress relation as transfer functions describing linear distributional processes in the time domain. These techniques are classical, although we will use them in a language that is borrowed from the recent literature of time domain integral equations. In a second part of the paper (Sections 8 to 10), we will introduce and use tools from the theory of strongly continuous semigroups to analyze the three classical models and some transitional situations where, for instance, classical Zener-style viscoelasticity coexists with pure linear elasticity in different subdomains, with smooth or abrupt transition regions. Our goals for the current piece of work are not in the realm of the modeling: we will analyze but not discuss known models, and we will not deal with physical justifications thereof. A particular issue where we will be very restrictive is the fact that we will only deal with solids moving from equilibrium (no displacement, strain, or stress) at time zero. There are practical reasons for this choice (since stress remembers past strain, it is not entirely justifiable to start the clock with known displacement and stress), but we are also restricted because of our analysis goals. In both parts (transfer function analysis and semigroup analysis) we will give a hint at how to deal with initial conditions.
Barring initial and boundary conditions that are needed to fully describe the model, our goal is the study of a linear elasticity equation
$$\begin{array}{}
\rho \ddot{\mathbf u}=\mathrm{div}\,\boldsymbol\sigma+\mathbf f
\end{array} $$
in a bounded domain of d-dimensional space. Here u is the displacement field, upper dots denote time differentiation, σ is the stress tensor and f represents the volumetric forces. Linear strain ε =
$\begin{array}{}
\frac{1}{2}
\end{array} $ (∇ u + (∇ u)⊤) determines stress through a generic convolutional law (we only display the time variable in the following formulas)
$$\begin{array}{}
\displaystyle \boldsymbol\sigma(t)=\int_0^t \mathcal D(t-\tau){\dot{\boldsymbol\varepsilon}(\tau)}\mathrm d\tau,
\end{array} $$
where 𝓓 is a time-dependent, possibly distributional, tensor-valued kernel. Following the careful description given by Francesco Mainardi in his monograph [19], the four classical models of viscoelasticity are named after Zener, Voigt, Maxwell, and Newton. The strain-to-stress relationship is given by a differential equation as follows:
$$\begin{array}{}
\boldsymbol\sigma+a\,\dot{\boldsymbol\sigma}
=\mathrm C_0 \boldsymbol\varepsilon+\mathrm C_1\dot{\boldsymbol\varepsilon}, \qquad
\mbox{(Zener)} \\
\qquad\,\,\,\boldsymbol\sigma
=\mathrm C_0 \boldsymbol\varepsilon+\mathrm C_1\dot{\boldsymbol\varepsilon}, \qquad
\mbox{(Voigt)}\\
\boldsymbol\sigma+a\,\dot{\boldsymbol\sigma}
=\mathrm C_1\dot{\boldsymbol\varepsilon}, \qquad \qquad\,\,\,\,\,
\mbox{(Maxwell)}\\
\qquad\,\,\,\boldsymbol\sigma
=\mathrm C_1\dot{\boldsymbol\varepsilon}, \qquad \qquad\,\,\,\,\,
\mbox{(Newton)}\\
\qquad\,\,\,\boldsymbol\sigma
=\mathrm C_0\boldsymbol\varepsilon. \qquad \qquad\,\,\,\,\,
\mbox{(Linear elasticity)}
\end{array} $$
Here a is a non-negative function and C0 and C1 are four-index tensors satisfying hypotheses similar to the tensor that is used to describe linear elasticity (we will be covering heterogeneous anisotropic solids), with some additional conditions that will be introduced when we explain the models in detail. Formally speaking, the last three models can be considered as particular examples of Zener’s model. However, they have very different properties. Newton’s model is equivalent to a parabolic equation, as can be seen by substituting the formula for σ = C1ε̇ in the equation of conservation of momentum and integrating once in the time variable. This model therefore does not produce waves and will be ignored in the sequel. Voigt’s model also gives an explicit strain-to-stress relation: if we substitute σ = C0ε+C1ε̇ in the dynamical equation we can see that we end up with a PDE of order three (there are terms with two derivatives in space and one in time). The models of Zener, Maxwell, and Voigt dissipate energy.
Let us now give a quick literature review, including some relevant work from the modeling and mathematical communities. Some of the monographs [1, 18, 19] contain a large collection of references that can be used for a deeper introduction to this fascinating area. For a very early introduction to linear viscoelasticity, see [10]. A generalized model that encompasses all classical models describe above four is examined in [1, Chapter 3]. An overview of the physics of viscoelasticity and its relation to rheology (fluid and soft solid flow) can be found in [29] and [23], where the latter also develops numerical methods for solving problems associated to viscoelasticity, while an overview of the mathematical theories and techniques, including the problem of waves propagating in viscoelastic media can be found in [8]. Viscoelastic models have also been formulated in the language of integral equations (see [8, 12, 32]), and electrostatic models like the Cole-Cole model in [16].
We wish to consider waves propagating in viscoelastic solids using both the classical models and those using fractional derivatives. The relation between fractional derivatives and viscoelastic models (including introductions to the Mittag-Leffler functions) can be found in [19, 22], and [20] offers a short survey of the history of development of this theory. Mainardi, who it can easily be seen is at the center of much of the development and communication of waves in viscoelasticity, shows in [21] that the fractional relaxation process with constant coefficients is equivalent to a similar process governed by variable coefficient ODE. This relationship has also been explored in [16], where fractional derivatives are bypassed by using a method of Yuan and Agrawal to convert the equations to a system of first order ODE. The Zener model has been extensively studied in the context of waves with both classical [2, 5, 14] and fractional [26] time derivatives, as well as in the quasistatic case [28, 30]. Other authors have explored the Maxwell and fractional Maxwell models [7], different variations of Maxwell models [6, 9], and the Voigt model [5, 12, 15].
For some of our stability results, we make rigorous use of Laplace transforms for vector-valued distrivutions. Laplace transforms also appear in [2, 8, 9, 10, 26, 28], used in the context of showing existence and uniqueness of solutions, justifications of models, exploration of the constitutive relations, or to simplify numerical implementation. In the context of stability analysis, a Green’s function representation of the solution to the three dimensional wave problem is used in [8]. While we obtain some estimates in the time domain by making use of an inversion theorem for the Laplace transformation (in the spirit of the Payley-Wiener theorems), in some cases can make use of semigroup theory to obtain estimates directly. Similar analysis has been performed for waves in unbounded domains [5] and in bounded domains [14], where semigroup analysis is used to show the existence and uniqueness of solutions for waves in a Zener model.
While not explored in detail in the present work, we are also interested in the analysis and implementation of numerical schemes for the simulation of waves in viscoelastic materials. As mentioned earlier, [23] contains an overview of numerical methods for problems in viscoelasticity including finite element, boundary element, and finite volume formulations. Numerical implementation with finite elements has also been explored in [12] and [35] for the simulation and comparison to real world data, specifically blood flow in [35]. Also in the context of blood flow, [28] uses discontinuous Galerkin methods to simulate a quasistatic nonlinear 1D fractional Zener model. A DG method for a general linear quasistatic viscoelastic model is proposed in [30] and a priori error estimates are derived. Convergence of finite element methods for viscoelasticity is explored [13] and [14, 15], where the first reference focuses on convergence in time while the latter two are concerned with optimal order of convergence in space. Coupling of elastic and viscoelastic subdomains is examined in [24] where boundary elements for the viscoelastic subdomain are coupled with finite elements for the elastic components, whereas in [34] a scheme involving only finite elements are used for the same problem and the two schemes are compared.
Our paper is structured as follows. After introducing the general model (Section 2), we give a general framework for the viscoelastic material law as a transfer function (Section 3) and then move on to prove that the main classical models (Zener, Maxwell, Voigt), fractional versions of them, and combinations of different models in different subdomains, fit in our general framework (Section 4). Sections 5-7 contain the Laplace domain analysis of the model carried out as follows: first we do a transfer function analysis, then we give the general theory of how to understand transfer functions as Laplace transforms of distributional convolutions in the time variable, and finally we give estimates for the case of smooth (in time) data. In Sections 8-10, we start with a semigroup analysis of the classical models. Because we are striving for generality, we make an effort to include models where the classical viscoelastic models can degenerate into classical elasticity, which motivates a careful discussion of the closure process of a normed space with respect to a certain seminorm and how this affects the action of some operators. Section 9 gives a detailed treatment of Zener’s model, using the tools of the previous section and well-known results of the theory of strongly continuous semigroups in Hilbert spaces. Section 10 sketches the main changes that need to be made to the preceding analysis to study Voigt’s and Maxwell’s models. Finally, and just for the sake of illustration, we show some simulations for one and three dimensional models.
Before we proceed with the work at hand, let us give here some quick notational pointers. While the transient models we will be describing and analysing in this paper take values on spaces of real valued functions, the transfer function analysis will require the introduction of complex variables and complex-valued functions. To be on the safe side, all brackets and forms considered in this paper will be linear or bilinear, never conjugate linear or sesquilinear. We will write
$$\begin{array}{}
\displaystyle \mathrm A:\mathrm B=\sum\limits_{i,j=1}^d a_{ij}b_{ij},
\qquad \mathrm A,\mathrm B\in \mathbb C^{d\times d},
\end{array} $$
with no conjugation involved. The upperscript ⊤ will be used for transposition of matrices, without conjugation. Note that A : B = 0 for all A ∈
$\begin{array}{}
\mathbb C^{d\times d}_\text{sym}
\end{array} $ := {A ∈ ℂd × d : A⊤ = A}, and B ∈
$\begin{array}{}
\mathbb C^{d\times d}_\text{skw}
\end{array} $ := {B ∈ ℂd × d : B⊤ = −B}. We will write ∥M∥2 = M : M.
Given two Banach spaces X and Y, we will consider the space 𝓑(X, Y) of bounded linear maps from X to Y with the operator norm. We will shorten 𝓑(X) := 𝓑(X, X).
An introduction to the model problem
The wave propagation problem will be given in an open bounded domain Ω ⊂ ℝd, whose boundary is denoted by Γ. In order to have a well defined trace operator in classical Sobolev spaces, we will assume that Ω is locally a Lipschitz hypograph, although this hypothesis can be relaxed as long as we have a trace operator. We will assume that Γ is decomposed into Dirichlet and Neumann parts, ΓD and ΓN, satisfying
$$\begin{array}{}
\overline{\Gamma_D}\cup \overline{\Gamma_N}=\Gamma, \qquad \Gamma_D\cap\Gamma_N=\emptyset.
\end{array} $$
The inner products in the Lebesgue spaces
$$\begin{array}{}
L^2(\Omega), \quad \mathbf L^2(\Omega):=L^2(\Omega;\mathbb R^d), \quad
\mathbb L^2(\Omega):=L^2(\Omega;\mathbb R^{d\times d}_\text{sym})
\end{array} $$
(note that the latter is a space of symmetric-matrix-valued functions) will be respectively denoted by
$$\begin{array}{}
\displaystyle (u,v)_\Omega:=\int_\Omega u\,v,\qquad
(\mathbf u,\mathbf v)_\Omega:=\int_\Omega \mathbf u\cdot\mathbf v,\qquad
(\mathrm S,\mathrm T)_\Omega:=\int_\Omega \mathrm S\,:\,\mathrm T,
\end{array} $$
and ∥⋅∥Ω will denote the associated norm in all three cases. We will also consider the Sobolev space
$$\begin{array}{}
\mathbf H^1(\Omega)=H^1(\Omega;\mathbb R^d)
:=\{ \mathbf u :\Omega \to \mathbb R^d\,:\, \mathbf u\in \mathbf L^2(\Omega),\,
\nabla \mathbf u \in L^2(\Omega;\mathbb R^{d\times d})\},
\end{array} $$
endowed with the norm
$$\begin{array}{}
\| \mathbf u\|_{1,\Omega}^2:=\| \mathbf u\|_\Omega^2+\|\nabla \mathbf u\|_\Omega^2,
\end{array} $$
and the symmetric gradient operator
$$\begin{array}{}
\mathbf H^1(\Omega)\ni \mathbf u\longmapsto
\boldsymbol\varepsilon(\mathbf u):=\tfrac12(\nabla \mathbf u+(\nabla\mathbf u)^\top)
\in \mathbb L^2(\Omega).
\end{array} $$
We can define a bounded and surjective trace operator γ : H1(Ω) → H1/2(Γ). For simplicity, we will write γDu := γu|ΓD. We then consider the Sobolev spaces:
$$\begin{array}{}
\,\,\,\quad\mathbf H^1_D(\Omega) := \{ \mathbf w\in \mathbf H^1(\Omega)\,:\, \gamma_D \mathbf w=0 \}=\ker\gamma_D,\\
\,\,\,\mathbf H^{1/2}(\Gamma_D) := \{ \gamma_D\mathbf u\,:\, \mathbf u \in
\mathbf H^1(\Omega)\} =\mathrm{range}\,\gamma_D,\\
\,\,\widetilde{\mathbf H}^{1/2}(\Gamma_N) := \{\gamma\mathbf w|_{\Gamma_N}\,:\, \mathbf w\in \mathbf H^1_D(\Omega)\},\\
\mathbf H^{-1/2}(\Gamma_N) := \widetilde{\mathbf H}^{1/2}(\Gamma_N)'.
\end{array} $$
Our exposition will include the cases where either ΓD or ΓN is empty. To be completely precise, H−1/2(ΓN) is the representation of the dual of H͠1/2(ΓN) making
$$\begin{array}{}
\widetilde{\mathbf H}{^{1/2}}(\Gamma_N)
\subset L^2(\Gamma_N;\mathbb R^d)
\subset \mathbf H^{-1/2}(\Gamma_N)
\end{array} $$
a well-defined Gelfand triple. The reciprocal duality product of the two fractional spaces on ΓN will be denoted with the angled bracket 〈⋅, ⋅〉ΓN.
We will consider the space for symmetric-tensor-valued functions
$$\begin{array}{}
\mathbb H(\mathrm{div},\Omega):=
\{ \mathrm S\in \mathbb L^2(\Omega)\,:\,
\mathrm{div}\,\mathrm S \in \mathbf L^2(\Omega)\}.
\end{array} $$
In the above definition, the divergence operator is applied to the rows of the matrix valued function S. Following well-known results on Sobolev spaces, we can define a bounded linear and surjective operator γN : ℍ(div, Ω) → H−1/2(ΓN) so that the following weak formulation of Betti’s formula
$$\begin{array}{}
\langle \gamma_N \mathrm S,\gamma\mathbf u\rangle_{\Gamma_N}
\!\!\!\!& =(\mathrm S,\nabla\mathbf u)_\Omega+(\mathrm{div}\,\mathrm S,\mathbf u)_\Omega\\
&=(\mathrm S,\boldsymbol\varepsilon(\mathbf u))_\Omega+(\mathrm{div}\,\mathrm S,\mathbf u)_\Omega \qquad
\forall \mathrm S\in \mathbb H(\mathrm{div},\Omega)
\quad\forall \mathbf u\in \mathbf H^1_D(\Omega),
\end{array} $$
holds.
Pending a precise introduction of the material law, which we will give in the Laplace domain in Section 3, we are now ready to give a functional form for the viscoelastic wave propagation problem. We look for u : [0, ∞) → H1(Ω) and σ : [0, ∞) → ℍ(div, Ω) satisfying for all t ≥ 0
$$\begin{array}{}
\rho\, \ddot{\mathbf u}(t) =\mathrm{div}\,\boldsymbol\sigma(t)+\mathbf f(t),
\end{array} $$$$\begin{array}{}
\gamma_D\mathbf u(t) =\boldsymbol\alpha(t),
\end{array} $$$$\begin{array}{}
\gamma_N\boldsymbol\sigma(t) =\boldsymbol\beta(t).
\end{array} $$
Here ρ ∈ L∞(Ω), with ρ ≥ ρ0 almost everywhere for some positive constant ρ0, models the mass density in the solid, which at the initial time t = 0 is at rest on the reference configuration Ω
$$\begin{array}{}
\mathbf u(0)=\mathbf 0, \qquad\dot{\mathbf u}(0)=\mathbf 0.
\end{array} $$
Upper dots are used to denote time derivatives. The data are functions
$$\begin{array}{}
\mathbf f:[0,\infty)\to \mathbf L^2(\Omega),
\qquad
\boldsymbol\alpha:[0,\infty)\to \mathbf H^{1/2}(\Gamma_D),
\qquad
\boldsymbol\beta:[0,\infty)\to \mathbf H^{-1/2}(\Gamma_N).
\end{array} $$
What characterizes the viscoelastic model (for small deformations where the strain at time t can be described by ε(u(t))) is the existence of a strain-stress relation of the form
$$\begin{array}{}
\displaystyle \boldsymbol\sigma(t)=\int_0^t \mathcal D(t-\tau)\,\boldsymbol\varepsilon(\dot{\mathbf u}(\tau))
\mathrm d\tau.
\end{array} $$
The viscoelastic law (1e) is formally written in terms of a convolutional kernel 𝓓. In a first approximation, this kernel can be considered as a fourth order tensor (with some symmetric properties) depending on the time variable. As we will see later on, the most interesting examples arise when the causal convolution (1e) is a distributional one and 𝓓 is described as a causal tensor-valued distribution of the real variable.
Viscoelastic laws in the Laplace domain
In this section we are going to give a precise meaning to the general viscoelastic law (1e). Instead of writing the convolutional law (1e) in the time domain, we are going to introduce a Laplace transformed model which we will analyze in detail. This model will then be used to justify a family of distributional models in the time domain. At this point, we consider the formal Laplace transform of the convolutional process (1e) and introduce
$$\begin{array}{}
\mathrm C(s):=s\, \mathcal L\{ \mathcal D\}(s).
\end{array} $$
(Note that multiplication by s takes care of time differentiation.) Laplace transforms will be defined in the complex half-plane
$$\begin{array}{}
\mathbb C_+:=\{ s\in \mathbb C\,:\, \mathrm{Re}\,s \gt 0\}.
\end{array} $$
The viscoelastic material tensor can be described as a transfer function through a holomorphic map
$$\begin{array}{}
\mathbf C: \mathbb C_+\to \mathcal B(\mathbb C^{d\times d};
L^\infty(\Omega;\mathbb C^{d\times d}))\equiv L^\infty(\Omega;\mathbb C^{d\times d\times d\times d}).
\end{array} $$
Given M ∈ ℂd × d we thus have ℂ+ ∋ s ↦ C(s)M ∈ L∞(Ω;ℂd × d).
Some easy observations: conditions (2b) and (2d) imply (2c); if
$\begin{array}{c}
\displaystyle
\mathbf{M}\in \mathbb C^{d\times d}_\text{skw}
\end{array}$
, then C(s)M = 0 for all s; and if M ∈ ℝd×d, then C(s)M ∈ ℝd×d for all s ∈ (0, ∞).
Hypothesis 2
(Positivity). There exists a non-decreasing functionψ : (0, ∞) → (0, ∞) satisfying
Here and in the sequel, we identify C(s) with its associated ‘multiplication’ operator
$\begin{array}{c}
\displaystyle
\mathbf C(s)\in \mathcal{B}(L^2(\Omega;\mathbb C^{d\times d}_\text{sym}))
\end{array}$
, noticing that
A precise time domain description of the strain-stress material law (6) will require the introduction of some tools of the theory of operator valued distributions. We will do this in Section 6.
Examples
Before detailing the main examples covered with our theory, let us introduce a definition that will make our exposition simpler. Let
where c > 0 is a constant. For simplicity, in the future, we will write C ≥ c to refer to the last inequality and we will say that C is a steady Hookean material model, when conditions (11) are satisfied. The constant c > 0 will be called a lower bound for the model. Note that we can apply the model to complex-valued matrices
Ifa ∈ L∞(Ω) satisfiesa ≥ a0 > 0 almost everywhere, thenaCis a steady Hookean model.
Proof
Note first that
$\begin{array}{c}
\displaystyle
\mathbf C\mathbf{M}\in \mathbb C^{d\times d}_\text{sym}
\end{array}$
almost everywhere for all M ∈ ℂd×d. Therefore, if
$\begin{array}{c}
\displaystyle
\mathrm N\in \mathbb C^{d\times d}_\text{skw}
\end{array}$
, then
which implies CN = 0. Property (b) follows from (a). Properties (c) and (d) are straightforward.
Elastic models
In the case where we take a Hookean model C0 and we consider the constant function C(s) ≡ C0, it is simple to see that the Hypotheses 1-3 are satisfied with ψ(x) := c0x (c0 being the lower bound for the model C0) and (4), r = 0, and ϕ(x) := ∥C0∥. The time domain version of this model is the usual linear strain-stress relation
This hypothesis makes Cdiff a non-strict Hookean model. As we will see in the direct time-domain analysis of Section 9, Cdiff is the diffusive part of the elastic model, while C0 acts as a base or ground elastic model. We will make use of the formula
The above exposition of Zener’s model allows for full anisotropy. The isotropic viscoelastic model can be easily described with two variable coefficients. To do that we let λ, μ : ℂ+ → L∞(Ω) be holomorphic functions with the following properties being satisfied almost everywhere in Ω and for all s ∈ ℂ+:
Examples of functions satisfying (14) can be found using a variant of Zener’s model for viscoelasticity: let a, mμ, bμ, mλ, bλ ∈ L∞(Ω) be strictly positive (bounded below by a positive number) and such that
In this section we explore models of the form C(sν) where C(s) is a Zener model and ν ∈ (0, 1). For fractional powers in the complex plane we will always take the principal determination of the argument, i.e., the one with a branch cut at the negative real axis. A fractional Zener model has the form
is a Caputo fractional derivative of order ν. Note that this fractional derivative coincides with the Riemann-Liouville fractional derivative of the same order, if we are assuming homogeneous initial conditions for all variables. This fractional derivative can also be defined as a distributional fractional derivative in the entire real line.
Proposition 3
(Fractional Zener models). LetC(s) be a viscoelastic Zener model and letν ∈ (0, 1).
Ifc0 > 0 is the lower bound for the tensorC0, thenC(sν) satisfies Hypothesis 2 withψ(x) := c0x.
The modelC(sν) satisfies Hypothesis 3 withr = 0 andϕgiven by(12).
Proof
To prove (a), using the same idea as in Proposition 2, we write
$$\begin{array}{c}
\displaystyle
\overline s \mathbf C (s^\nu) \mathbf{M}:\overline {\mathbf{M}}
=
(\mathbf C_0\mathbf{M}:\overline{\mathbf{M}}) \, \overline s +
(|s|^2 \mathbf C_\text{diff}\mathbf{M}:\overline{\mathbf{M}})\,(1+as^\nu)^{-1} s^{\nu-1}.
\end{array}$$
where C1 is a steady Hookean model (with lower bound c1 > 0) and a ∈ L∞(Ω) satisfies a ≥ a0 > 0 almost everywhere, for some constant a0. In the time domain this gives again an implicit strain-to-stress relation
as a viscoelastic parameter model, where C0 is a steady Hookean material model and C1 is a non-negative Hookean model. In those parts of the domain where C1 = 0, Voigt’s model reduces to classical linear elasticity. In the time domain, this model gives an explicit differential expression for the strain-to-stress relationship
which shows that this model is a third order differential equation, although we admit the possibility that the third order terms vanish in some regions.
Proposition 7
Voigt’s model satisfies Hypotheses 1-3 withr = 1, and
Let Ω1, …, ΩJbe non-overlapping subdomains of Ω such that$\begin{array}{}
\displaystyle
\overline\Omega=\cup_{j=1}^J \overline\Omega_j
\end{array}$
. Assume that Cjis a viscoelastic model in the domain Ωj, satisfying the Hypotheses 1-3. Then
This means, in particular, that we can combine all the above models (Zener, Maxwell, and Voigt) in their differential or fractional versions, with different fractional orders in different subdomains subdomains. Because all our definitions are distributional, whenever there is an interface (smooth or not) we are implicitly imposing continuity of the displacement (this is done by assuming u too take values in H1(Ω)) and of the normal stress (since σ takes values in ℍ (div, Ω)).
Newton’s model is a fourth choice among classical viscoelastic models, given by the law
(here g is an antiderivative of f and takes care of non-vanishing initial conditions if needed). Therefore, Newton’s model becomes a parabolic equation for linear elasticity and does not produce waves. Since the interest of this paper is wave models, we will not investigate this simple model any further.
Transfer function analysis
We consider the energy norm, tagged in a parameter c > 0:
These are d uncoupled scalar problems for each of the components of u͠. Using the Bamberger-HaDuong lifting lemma (originally stated in [3], see [31, Proposition 2.5.1] for a rephrasing in the current language), it follows that
The main theorem of this section studies the operator associated to the Laplace transform of problem (1). In fact, problem (23) below is the Laplace transform of (1) for data that as functions of time are Dirac masses at time equal to zero.
Theorem 12
Letf ∈ L2(Ω;ℂ), α ∈ H1/2(ΓD;ℂ), andβ ∈ H−1/2(ΓN;ℂ). For s ∈ ℂ+, the solution of
$$\begin{array}{}
\displaystyle
\mathbf u \in \mathbf H^1(\Omega;\mathbb C),\qquad \gamma_D \mathbf u =\boldsymbol\alpha,
\end{array}$$
for a certain constantCdepending onρand the geometry.
Proof
The solution of (23) for data (f, β, α) can be decomposed as the sum of the solutions for (f, 0, 0), (0, β, 0) and (0, 0, α). For the first one, we note that u ∈
$\begin{array}{}
\displaystyle
\mathbf H^1_D
\end{array}$
(Ω;ℂ) satisfies
where we have used Korn’s inequality and (1). For the final one, we write u = û + u0, where û = û(α, |s|) is the solution of (21) with c = |s| and u0 ∈
$\begin{array}{}
\displaystyle
\mathbf H^1_D
\end{array}$
(Ω;ℂ). Then
In this section we show how the transfer function studied in Section 5 (specifically in Theorem 12) is the Laplace domain transform of the solution operator for a distributional version of the viscoelastic wave propagation problem (28) (cf. Proposition 16 below). We start with some language about vector-valued distributions, borrowed from [31].
Background on operator valued distributions
Given a real Hilbert space X, its complexificationXℂ := X + ıX is a complex Hilbert space that is isometric to X × X
with the product by complex scalars defined in the natural way. The complexification Xℂ has a naturally defined conjugation, which is a conjugate linear isometric involution in Xℂ. The Lebesgue and Sobolev spaces of complex-valued functions that we have used in Section 5 are complexifications of the corresponding real spaces. If X and Y are real Hilbert spaces, the space of bounded linear operators 𝓑(Xℂ, Yℂ) can be understood as the subspace of 𝓑(X2, Y2) formed by matrices of operators of the form
$$\begin{array}{}
\displaystyle
\begin{pmatrix} A & - B \\ B & A \end{pmatrix} \qquad A,B\in \mathcal B(X,Y).
\end{array}$$
This is easily seen to be isomorphic (with an equivalent but not equal norm) to the complexification of the Banach space 𝓑(X, Y). (For the problem of the many possible equivalent complexifications of real Banach spaces, see [25].) If A ∈ 𝓑(Xℂ, Yℂ) we can define A ∈ 𝓑(Xℂ, Yℂ) by
$$\begin{array}{}
\displaystyle
\overline A x:=\overline{A\overline x},
\end{array}$$
where in the right hand side we use the natural conjugations of Xℂ and Yℂ. With this definition Ax = Ax.
An X-valued tempered distribution is a continuous linear map from the Schwartz class 𝓢(ℝ) to X. We say that the X-valued tempered distribution h is causal, when the action of h on any element of 𝓢(ℝ) supported in (−∞, 0) is zero. Causal tempered distributions have a well defined Laplace transform, which is a holomorphic function H = 𝓛{h} : ℂ+ → Xℂ satisfying
where the conjugation on the left-hand side is the one in Xℂ. We will write h ∈ TD(X) whenever h is an X-valued causal tempered distribution whose Laplace transform satisfies
where μ ∈ ℝ and ψ : (0, ∞) → (0, ∞) is non-increasing and at worst rational at the origin, i.e., sup0<x<1xkψ(x) < ∞ for some k ≥ 0. When, instead of a Hilbert space X and its complexification Xℂ, we are dealing with bounded linear operators 𝓑(Xℂ, Yℂ), the conjugation in (24) is the one for operators between complexified spaces. Some pertinent observations and results:
If h ∈ TD(X) and T ∈ 𝓑(X, Y) is a ‘steady-state’ operator, then Th ∈ TD(Y). We can reverse the roles of T and h and show that if an operator valued distribution T ∈ TD(𝓑(X, Y)) acts on a constant h ∈ X, it defines Th ∈ TD(Y).
A simple argument using the formula for the inverse Laplace transform of s−mH(s), where m is an integer chosen so that μ − m < − 1, can be used to characterize all these distributions (see [31, Proposition 3.1.2]): h ∈ TD(X) if and only if there exists an integer m ≥ 0 and a causal continuous function g : ℝ → X with polynomial growth at infinity such that h = g(m), with differentiation understood in the sense of tempered distributions. Moreover, if H : ℂ+ → Xℂ is a holomorphic function satisfying (24) and (25) (with the conditions given for ψ), then H = 𝓛{h} for some h ∈ TD(X).
If h ∈ TD(ℝ) and a ∈ X, then the tensor product a ⊗ h defines a distribution in TD(X).
If X and Y are Banach spaces and h ∈ TD(𝓑(X, Y)), then the convolution producth ∗ λ is well defined for any X-valued causal distribution λ, independently on whether it is tempered or not [33]. In the simpler case where λ ∈ TD(X), we have the convolution theorem
which can be used as an equivalent (and simple) definition of the convolution of h ∗ λ, and we also have h ∗ λ ∈ TD(Y), as can be easily proved from the definition.
Viscoelastic material law and wave propagator
Theorem 13
If C : ℂ+ → L∞(Ω;ℂd×d×d×d) satisfies hypotheses (2), then there exists
such that 𝓛{𝓓M}(s) = C(s)M for all$\begin{array}{}
\displaystyle
\mathrm M\in \mathbb R^{d\times d}_\text{sym}.
\end{array}$For arbitraryu ∈ TD(H1(Ω)), the convolution
such that 𝓛{𝓓} = C. If we fix
$\begin{array}{}
\displaystyle
\mathrm M\in \mathbb R^{d\times d}_\text{sym}
\end{array}$
, we can easily show that the Laplace transform of 𝓓M ∈ TD(𝕃2(Ω)) is C(s)M. Finally, if u ∈ TD(H1(Ω)), then ε(u) ∈ TD(𝕃2(Ω)) and the convolution 𝓓 ∗ ε(u) is well defined.
The expression σ = 𝓓 ∗ ε(u) can be equivalently written σ = 𝓓 ∗ ε(u̇) where 𝓓 ∈ TD(𝓑(𝕃2(Ω))) is the distribution whose Laplace transform is s−1 C(s). In the simplest example (the purely elastic case), C(s) = C0, we can write 𝓓 = C0 ⊗ δ0 and 𝓓 = C0 ⊗ H, where H is the Heaviside function. This yields the usual elastic law σ = C0ε(u).
Proposition 14
Fors ∈ ℂ+, let us consider the solution map
$$\begin{array}{}
\displaystyle
\mathrm S(s):\mathbf L^2(\Omega;\mathbb C)\times \mathbf H^{1/2}(\Gamma_D;\mathbb C)\times \mathbf H^{-1/2}(\Gamma_N;\mathbb C)
\to \mathbf H^1(\Omega;\mathbb C)
\end{array}$$
defined byu = S(s)(f, α, β) being the solution of (23). The function
is analytic. Moreover, the operator
$\begin{array}{}
\displaystyle
\mathbb A_0(s):=\mathbb A(s)|_{\mathbf H^1_D(\Omega;\mathbb C)}
\end{array}$
is invertible for all s ∈ ℂ+ by the coercivity of the bilinear form b given in (19a). For any two Banach spaces, the inversion operator
is 𝓓∞ and therefore the map s → 𝔸0(s)−1 is analytic from ℂ+ to
$\begin{array}{}
\displaystyle
\mathcal B(\mathbf H^1_D(\Omega;\mathbb C)',\mathbf H^1_D(\Omega;\mathbb C))
\end{array}$
. If we now consider a bounded operator L : H1/2(ΓD) → H1(Ω) that is a right-inverse of γD and its natural extension to complex-valued functions, we can easily write
Since div σ = ρs2 S(s)(f, α, β) − f and C is analytic, it is clear that T is analytic. The conjugation property for T and the formula for its inverse are straightforward.
The above estimates give an upper bound for the norm of ∥T(s)∥, which together with Proposition 15 shows the existence of 𝓣 such that T = 𝓛{𝓣}. To prove that (u, σ) = 𝓣 ∗ (f, α, β) solves (28), we can just take Laplace transforms and use the definition of T(s). We next give an alternative proof that will be used for some arguments later on. Note first that σ = 𝓓 ∗ ε(u), as follows from the definition of T(s). Note also that there exists
for all s ∈ ℂ+, where F = 𝓛{f}, A = 𝓛{α}, and B = 𝓛{β}. Since (29) is uniquely solvable, this proves uniqueness of (28).
This general result about the weak (distributional) version of the viscoelastic wave propagation problem includes the possibility of adding non-homogeneous initial conditions for the displacement and the velocity, since they are just included in the volume forcing function f. These non-zero conditions will make the distributional solution of (28) non-smooth at time t = 0 and therefore the estimates of Section 7 will not be valid.
Estimates in time
In this section we translate the estimates for the transfer function given in Theorem 12 into time-domain estimates. The following theorem rephrases [31, Proposition 3.2.2], which is an inversion theorem for the Laplace transform of causal convolution operators.
Theorem 17
Let F : ℂ+ → 𝓑(X, Y) be a holomorphic function valued in the space of bounded linear operators between two Hilbert spaces and assume that
φ:(0, ∞) → (0, ∞) is non-increasing and sup0<x<1xkφ(x) < ∞ for somek ≥ 0.
Ifλ ∈ 𝓓m+1(ℝ;X) is a causal function such thatλ(m+2)is integrable, then the uniqueY-valued causal functionusuch that 𝓛{u} = F 𝓛{λ} is continuous and satisfies
The integrability of the (m + 2)-th derivative of λ can be relaxed to local integrability, but in that case we need to add a hypothesis concerning the growth of ∥λ(m+2)(t)∥X, since we need to be able to take Laplace transforms in ℂ+. In any case, since all the processes that we are dealing with are causal (the solution depends on the past values of the data, and never on the future ones), the result holds in finite intervals of time assuming only local integrability of the last derivative. Note also that if X and Y are complexifications of real spaces, the process in the time domain is real-valued for real-valued data when we have F(s) = F(s) for all s.
Let F be as in Theorem 17 and$\begin{array}{}
\displaystyle
\lambda\in W^{m+2,1}_+(0,\infty;X).
\end{array}$Ifλ͠ : ℝ → Xis the trivial extension ofλto (− ∞, 0), then there is a uniqueu ∈ 𝓓([0, ∞);Y) such thatu(0) = 0 and 𝓛{u͠} = F 𝓛{λ͠}. Moreover, the estimates(31)hold.
The main theorem of this section follows. It will use the antidifferentiation-in-time operator
for some non-negative m(α), m(β), and m(f). We denote with the same symbols their tilde-extensions, i.e., their extensions by zero to negative times. We finally associate the solution (u, σ) of the viscoelastic wave propagation problem. We want to give hypotheses on the data guaranteeing the existence of solutions of the equation that are continuous functions of time. Key quantities to keep in mind are the parameter r and the functions ϕ and ψ in (3)-(4) (hypotheses on the material model) and the derived functions ϕ* and ψ* in Proposition 10. The key theorem to keep in mind is Theorem 12 and the bound for ∥C(s)∥.
The result is a more or less direct consequence of Corollary 18, using Theorem 12 for the Laplace domain estimates, and going through the language of Propositions 14 and 15. Since the problem is linear, it can be decomposed as the sum of problems with data (f, 0, 0), (0, α, 0), and (0, 0, β). We follow the notation of Proposition 14 and define nine instances of spaces and operators to apply Corollary 18: we list the spaces X and Y (before complexification), as well as the value of m and μ and the function φ in the estimate (30). We separate the operator S(s) in Proposition 14 as a sum of three operators
This proves the continuity properties for u and σ as well as the estimates (32a) and (32c). Note that zero initial values hold whenever the output function is continuous (see Corollary 18). To obtain the bounds for ∥u(t)∥Ω we can use a simple shifting argument, since if (f, α, β) ↦ u̇ (i.e., we use the operators multiplied by s and with values in L2(Ω)), then (f(−1), α(−1), β(−1)) ↦ u.
In Section 9, we give a different analysis, based on the theory of C0-semigroups of operators, of the viscoelastic wave propagation for Zener’s model, including situations where part of the domain is described with a purely elastic material. This requires a certain amount of preparatory work, which we will present in full detail. To avoid keeping track of constants, in this section and in the next we will use the symbol ≲ to absorb constants in inequalities that we do not want to display.
Let Cdiff ∈ L∞(Ω;ℝ(d×d)×(d×d)) be a non-negative Hookean model, which will play the role of the diffusive part of the viscoelastic model. We will identify the tensor with the operator
which is bounded, selfadjoint, and positive semidefinite. We also consider a strictly positive function a ∈ L∞(Ω), and the associated multiplication operator
$$\begin{array}{}
\displaystyle
\mathbb{L}^2(\Omega)\ni \mathrm E \longmapsto m_a \mathrm E:= a \mathrm E\in \mathbb{L}^2(\Omega).
\end{array}$$
Note that maCdiff = Cdiffma. We next consider the following seminorm in 𝕃2(Ω)
The result follows then by (33). The property (c) is a direct consequence of (b). Finally, to prove (d) we write
$$\begin{array}{}
\displaystyle
|b\,\mathrm S|_c^2
=(a\mathrm {C_{diff}}(b\mathrm S),b \mathrm S)_\Omega
=(b^2 a \mathrm {C_{diff}}\mathrm S,\mathrm S)_\Omega \\
\displaystyle\qquad~=\int_\Omega b^2 a \underbrace{(\mathrm {C_{diff}} \mathrm S):\mathrm S}_{\ge 0}
\lesssim \int_\Omega a (\mathrm {C_{diff}} \mathrm S):\mathrm S =|\mathrm S|_c^2.
\end{array}$$
This finishes the proof.
The dynamics of the viscoelastic model will allow for subdomains where energy is conserved and subdomains where the diffusive part of the model is active. The fact that we allow for transition areas between purely elastic and strictly diffusive models forces us to take the completion of the space 𝕃2(Ω) with respect to the seminorm |⋅|c. This standard process requires first to eliminate the kernel of the seminorm and move to its orthogonal complement. We thus consider the subspace
which is clearly closed in 𝕃2(Ω). If T ∈ M and |T|c = 0, then T ∈ (kerCdiff)⊥ ∩ kerCdiff = {0} (Proposition 20(c)) and therefore |⋅|c is a norm in M. Note also that if S ∈ 𝕃2(Ω), we can decompose
and CdiffS = CdiffSM. We define the Hilbert space M̂ as the completion of M with respect to the norm |⋅|c. Consider next the canonical injection I : M → M̂ and the operator Π : 𝕃2(Ω) → M that performs the 𝕃2(Ω) orthogonal projection onto M, where in M we consider the norm |⋅|c. (For the sake of precision, it is important to understand that the target space for Π is M and therefore Π is surjective.) We then define
$$\begin{array}{}
\displaystyle
\mathrm R:=\mathrm I\, \Pi:\mathbb{L}^2(\Omega)\to \widehat M,
\end{array}$$
which is continuous by Proposition 20(a). Note next that Cdiff|M : M → 𝕃2(Ω) is bounded by Proposition 20(b) and, therefore, there exists a unique bounded extension
We now consider b ∈ L∞(Ω) and the associated multiplication operatormb : 𝕃2(Ω) → 𝕃2(Ω) that we recall was given by mb E := b E. Since mb Cdiff = Cdiffmb, it follows that mb maps kerCdiff into kerCdiff. Therefore, if S ∈ M, then
and (c) follows by density. Property (d) is straightforward by a density argument.
We end this section with a very simple example of the above construction, where the diffusive tensor Cdiff is strictly positive in a subdomain and vanishes in the complement. Assume that there are two open sets Ω1 and Ω2 such that
$$\begin{array}{}
\displaystyle
\mathrm {C_{diff}} \mathrm T
=\chi_{\Omega_1}\mathrm {C_{diff}}\mathrm T
=\mathrm {C_{diff}}(\chi_{\Omega_1}\mathrm T).
\end{array}$$
If CdiffT = 0, then χΩ1 T : T = 0 by (35a) and therefore T = 0 almost everywhere in Ω1. Reciprocally, if T = 0 in Ω1, then by (36), CdiffT = 0. We have thus proved that
and therefore M̂ = M. In this case R : 𝕃2(Ω) → M̂ is just the restriction to Ω1 of matrix-valued functions defined on Ω,
$\begin{array}{}
\displaystyle
\widehat{\mathrm {C_{diff}}}
\end{array}$
is the restriction of the action of Cdiff to functions defined on
$\begin{array}{}
\displaystyle
L^2(\Omega_1;\mathbb R^{d\times d}_\text{sym})
\end{array}$
, and the same happens to multiplication operators
$\begin{array}{}
\displaystyle
\widehat{m_b}
\end{array}$
. This simple situation arises when the domain is subdivided into two parts: in one part we will deal with a purely elastic medium (Cdiff = 0), while in the other part we will handle a viscoelastic medium that is ‘strictly’ diffusive, as expressed in (35a).
Semigroup analysis of a general Zener model
In the coming two sections we will use some basic results on C0-semigroups in Hilbert spaces, namely the Lumer-Philips theorem (characterizing the generators of contractive semigroups in Hilbert spaces as maximal dissipative operators) and existence theorems for strong solutions of equations of the form
$$\begin{array}{}
\displaystyle
\dot U(t)=\mathcal A U(t)+F(t),
\end{array}$$
where 𝓐 : D(𝓐) ⊂ 𝓗 → 𝓗 is maximal dissipative and F : [0, ∞) → 𝓗 is a sufficiently smooth function. These results are well known and the reader is referred to any classical book dealing with semigroups (for instance, Pazy’s popular monograph [27]) or to the chapters on semigroups in many textbooks on functional analysis.
Consider now the space
$$\begin{array}{}
\displaystyle
\mathcal H:=\mathbf L^2(\Omega)\times \mathbb{L}^2(\Omega) \times \widehat M,
\end{array}$$
where we have progressively applied the boundary condition
$\begin{array}{}
\displaystyle
\gamma_N (\mathrm C_0\mathrm E+\widehat{\mathrm {C_{diff}}}\mathrm S)=0
\end{array}$
, Proposition 21(b), and Proposition 21(c). This proves the dissipativity of 𝓐.
Take now (f, F, G) ∈ 𝓗, solve the coercive variational problem
so that
$\begin{array}{}
\displaystyle
\mathrm S+\widehat{m_a}\mathrm S=\mathrm R\boldsymbol{\varepsilon}({\mathbf u})+\widehat{m_a}\mathrm G
\end{array}$
and therefore
$$\begin{array}{}
\displaystyle
\mathrm S=\widehat{m_a}^{-1}(\mathrm R\boldsymbol{\varepsilon}({\mathbf u})-\mathrm S)+\mathrm G.
\end{array}$$
The key issue for this proof is regularity. Assuming the data regularity of Theorem 23, we have (47), as continuity of u as a function [0, ∞) → H1(Ω) follows from (39b). Note that (48c) is (39d). By (39a), we have that div(C0E+
$\begin{array}{}
\displaystyle
\widehat{\mathrm C_\text{diff}}
\end{array}$S) ∈ 𝓒([0, ∞);L2(Ω)) and u̇(0) = 0, which was the missing initial condition (recall Corollary 24). We also have
Equations (49) and (50) identify continuous functions of t taking values in L2(Ω), H–1/2(ΓN), and 𝕃2(Ω) respectively. We can then differentiate them in the sense of vector-valued distributions of t to obtain (48a), (48d), and (48b). Note that to be entirely precise, the additional regularity we obtain is
Let$\begin{array}{}
\displaystyle
\boldsymbol\alpha\in W^{3,1}_+(0,\infty;\mathbf H^{1/2}(\Gamma_D)), \boldsymbol\beta\in W^{2,1}_+(0,\infty;\mathbf H^{-1/2}(\Gamma_N)),
\end{array}$andf ∈ W1,1(0, ∞;L2(Ω)). Additionally, let (u, E, S) solve (39) andσbe defined by(46). We have
and equations (48) hold for alltwith all derivatives defined in the strong way in the appropriate spaces.
Proof
If we solve problem (39) with (ḟ, α̇, β̇) as data, and we integrate from 0 to t, we obtain the solution of (39) which is therefore an element of 𝓒2([0, ∞);𝓗). This is enough to prove everything else.
The estimates of Corollary 24 greatly improve those of Section 7 (see Theorem 19) in two aspects: less regularity required for the data, and constants independent of the time variable. There are three particular cases included in the analysis of this section that are worth paying special attention to.
If ker Cdiff = {0}, then M = 𝕃2(Ω) and R = I is just the canonical inclusion of 𝕃2(Ω) into its completion with respect to the norm$\begin{array}{}
\displaystyle
(a{\mathrm C_{\text{diff}}}\,\cdot\,,\,\cdot\,)_\Omega^{1/2}.
\end{array}$
there is no need to use the completion process since M = M̂ = 𝕃2(Ω) and then
$\begin{array}{}
\displaystyle
\widehat{\mathrm C_\text{diff}}
\end{array}$ = Cdiff,
$\begin{array}{}
\displaystyle
\widehat{m_a}=m_a,
\end{array}$, and R is the identity operator. The space
which is equivalent to the usual norm. This makes the analysis of this strictly diffusive viscoelastic problem much simpler.
When Cdiff = 0, we have M = M̂ = {0} and the third equation and unknown do not play any role. In this case the operators ±𝓐 are maximal dissipative and therefore, 𝓐 is the generator of a group of isometries in 𝓗, i.e., this model is conservative. This should not be a surprise, since in this case we recover a first order formulation
The introduction of non-zero initial conditions for the most general version of this model is not trivial. When the model is strictly diffusive (case (b) in the above discussion, i.e., when (51) holds), we are allowed to impose initial conditions
which would be the natural ones for the formulation (48). This is done by modifying the system (39), using initial conditions u(0) = u0, E(0) = 0, S(0) = 0, and adding ρv0 to the right-hand side of (39a) and G0, with
to the right-hand-side of (39c). The general case is much more complicated, given that the initial conditions for σ(0) must match C0ε(u0) in purely elastic subregions.
More semigroup analysis
We are now going to take advantage of the preparatory work of Section 8 to give a quick view of the formulations and estimates that can be obtained for the Maxwell and Voigt models.
Maxwell’s model
Maxwell’s model can be understood as the particular case of Zener’s model when C0 = 0 and Cdiff is strictly positive. However, its analysis is not included in the treatment given in Section 9, due to the fact that C0 was used to define the norm of the space 𝓗. We thus start again, with a new space
Using the operator 𝓐 in an equation of the form (43), we can prove that the hypotheses of Theorem 23 are sufficient to provide a solution (u, S) ∈ 𝓒1([0, ∞);𝓗) of the problem
with vanishing initial conditions. The estimates of Corollary 24 hold for this model as well.
A combination of Zener’s and Maxwell’s models is also available. It requires an even more general framework so that C0 and Cdiff can vanish on separate parts of the domain as long as a certain combination stays strictly positive. (See the variational problem (38) that is solved as a starting step to prove maximal dissipativity. As long as
$\begin{array}{}
\displaystyle
\mathrm C_0+m_{1+a}^{-1}{\mathrm C_{\text{diff}}}
\end{array}$ is a Hookean model, everything else will work.) The analysis would require a completion process with respect to the seminorm
$\begin{array}{}
\displaystyle
(\mathrm C_0\cdot,\cdot)_\Omega^{1/2}
\end{array}$ and the corresponding restriction operator. This is a simple (while a little cumbersome) extension that the reader can do to prove their handle of the techniques developed above.
Voigt’s model
The analysis of Voigt’s viscoelastic model (Zener’s model with a = 0, C0 strictly positive and Cdiff ≥ 0), including areas transitioning to classical linear elasticity, follows from a simple modification of the ideas of Section 9. In Voigt’s model Cdiff = C1 plays the role of a dissipative term. The semigroup analysis of this model is slightly different in involving a second order differential operator. Like in Maxwell’s model, there is no need to involve a completion process to handle transitions to classical linear elasticity. We now consider the following ingredients:
and define E := ε(u) + F, it is easy to prove that (u, E) ∈ D(𝓐) and (u, E) – 𝓐(u, E) = (f, F). Therefore, 𝓐 is maximal dissipative.
Going carefully over the proof of Theorem 23, it is easy to see that with the same hypotheses on the data f, α, and β, we have a unique (u, E) ∈ 𝓒1([0, ∞);𝓗), vanishing at zero, and solving
and the resulting pair (u, σ) is a solution to (48) (with a = 0) satisfying also the estimates of Corollary 24.
Some experiments
We now present some numerical experiments of the various viscoelastic models that we described in Section 4. We use finite elements for space discretization and a trapezoidal rule-based convolution quadrature (TRCQ) for time discretization [4, 11, 17]. The numerical experiments will be divided into three groups. First, we investigate 1D uniaxial wave propagation. Through these experiments, we observe how the viscoelastic behavior is dependent upon the choosing of parameters in the constitutive equations. In the second group of experiments we compare the behaviors of 1D uniaxial wave propagation in elastic, classical viscoelastic, fractional viscoelastic, and heterogeneous models by plotting their 2D space-time contour graphs. In the heterogeneous model, we decompose the region into two subdomains with different viscoelastic models, where the reflection and refraction of waves can be observed at the transition interface. Finally, we present the 3D simulation of a viscoelastic rod. Similar to the heterogeneous domain in the previous experiments, the rod is decomposed into two different subdomains. The snapshots we present show how the simulation accurately captures the memory and relaxation effects of the rod under a sudden change in displacement. The numerical analysis of the discretization schemes employed in this section is the goal of future research. Tests have been performed in sufficiently refined space-and-time meshes to obtain some sort of eye-ball convergence to a solution.
1D experiments
For simplicity, in the one-dimensional examples we will use traditional PDE notation, as opposed to the notation of evolutionary equations used throughout the paper. We first present numerical experiments for different fractional models in one dimension
where the constitutive relation that determines the model is defined through σ and ρ. We use the window function displayed in the left of Figure 1 as Dirichlet boundary data at x = 0, while we take a homogeneous Neumann boundary condition at x = 1. The strain-to-stress relation is given by a general formula
for parameters C0, C1, a and ν to be determined. For the discretization in space we use 𝓟4 finite elements on a mesh with 513 subintervals of equal size. Discretization in time is carried out using a TRCQ with 10, 240 time-steps of equal size in the interval [0, 40]. For the first two sets of experiments, we implement the Dirichlet boundary condition g(t) in Figure 1 at x = 0 and observe the evolution of u(1, t). We use: (a) the values given in Table 1, varying C1, to create the results of Figure 2, and (b) the values in Table 2, varying ν, for the results of Figure 3. In Table 1 the values for C1 are chosen so that the parameter Cdiff, which controls the diffusion, takes same values for the three different models, except for Maxwell, where Cdiff = 0 results in a model that is identically zero. To avoid this while still getting comparable results, we use C1 = 0.05 for the first value in the Maxwell model. In Figure 2, as we increase C1, all three models show a faster energy dissipation. Compared to the Zener and Voigt models, the Maxwell model exhibits less oscillations as a response to the Dirichlet boundary condition.
The parameters used to create the plots given in Figure 2. The first choice of C1 for the Zener and Voigt models reduce the model to linear elasticity
Zener
Maxwell
Voigt
Co
1.5
0
1.5
Cl
0.75,1,2.75
0.05,0.25,2
0,0.25,2
a
0.5
0.5
0
ν
1
1
1
ρ
1
1
1
The parameters used to create the plots given in Figure 3.
Zener
Maxwell
Voigt
Co
1.5
0
1.5
Cl
1
1
1
a
0.5
0.5
0
ν
0.05,0.5,0.95
0.05,0.5,0.95
0.05,0.5,0.95
ρ
1
1
1
In Figure 3, we observe that the decreasing the fractional power ν leads to a slower rate of energy dissipation. We also notice that the change of ν does not have any obvious effect on the frequency of the oscillations.
We now change the boundary conditions (52b) and (52c) to
where h is the function in the right of Figure 1. Once again, we vary ν to see the effects that the fractional order of the derivative has on the Zener, Maxwell, and Voigt models respectively. The parameters we choose for this experiment is given in Table 3 and we again plot u(1, t) in Figure 4.
The parameters used to create the plots given in Figure 4.
Zener
Maxwell
Voigt
Co
1
0
1
Cl
1
1
1
a
0.5
0.5
0
ν
0.25,0.5,0.75,1
0.25,0.5,0.75,1
0.25,0.5,0.75,1
ρ
1
1
1
Here, both Zener and Voigt models show exponential rate response to the suddenly applied traction h(t) and then converge to some equilibrium state. The Maxwell model, as expected, shows a linear rate of creep when the traction is constant in time. As the fractional power ν decreases, we observe that the amplitude of oscillations is larger for the Zener and Voigt models, on the other hand we see the creep rate is slower for the Maxwell model.
Spacetime plots
We now study space-time contour plots of the displacement solution of (52) with Zener, Maxwell and Voigt models. We show three simulations focusing on one of these models where in each experiment we compare it with its fractional version, and observe its behavior in a heterogenous domain coupled with an elastic model. When we work on a heterogenous domain we split the interval [0, 1] into [0, 1/2) and [1/2, 1], where the first half is elastic and the second half is one of the models we are comparing: Zener, Maxwell or Voigt. We show space-time plots of the displacement corresponding to different models side by side where elastic model is included in all cases for the sake of comparison. We implement two different signals as Dirichlet boundary condition: a single pulse and a train of pulses (see Figure 5). For each experiment we display eight plots where the first four are the results of a single pulse while the last four are the results of the same experiment but for a train of pulses.
Our first test focuses on the Zener model where we utilize the parameters given in Table 4, and Figure 6 shows the outcome of this experiment. Here, the contrast between the energy conservative elastic model and the dissipative Zener model can easily be seen. We also notice that fractional Zener model displays slower dissipation than the classical model. In the heterogenous domain we observe a reflection and refraction of waves at the interface x = 1/2. Although the elastic model is conservative, in the case of a heterogenous domain we see the dissipation of the Zener model affecting the coupled system with a loss in wave amplitude. We also observe that in the results whose Dirichlet boundary condition is a train of pulses (the four panels on the right), the solution enters quickly into a time-harmonic regime.
The parameters used in in the fractional Zener model simulations to create the 2D spacetime plots shown in Figure 6.
Elastic
Zener
Fractional Zener
Heterogeneous Domain
Co
1.75
1.5
1.5
1.75 (x < 1/2), 1.5 (x ≥ 1/2)
C1
1.75
1.75
1.75
1.75
a
1
0.5
0.5
1 (x < 1/2), 0.5 (x ≥ 1/2)
ν
1
1
0.3
1
ρ
10
10
10
10
We perform a similar comparison for the Maxwell model using the parameters given in Table 5 and collecting the results in Figure 7. When a single pulse is used, we notice that the waves in the Maxwell model exhibit dissipation. Moreover, the fractional Maxwell reveals less dissipation with a little dispersion. In the heterogenous domain the reflections at x = 1/2 are less obvious for both of the input signals when comparing to the corresponding Zener simulations. Lastly, we demonstrate the space-time plots corresponding to the Voigt model in Figure 8 using the parameters from Table 6. Comparing to Zener and Maxwell, we observe that this model displays more dispersion. In particular, this dispersion is on dramatic display in the heterogenous domain where the reflections are clearly seen in the elastic part, while the dispersion occurs in the Voigt subdomain.
The parameters used in the Maxwell model simulations to create the 2D spacetime plots shown in Figure 7.
Elastic
Maxwell
Fractional Maxwell
Heterogeneous Domain
Co
1.75
0
0
1.75(x < 1/2), 0 (x ≥ 1/2)
C1
1.75
1.75
1.75
1.75
a
1
1
1
1
ν
1
1
0.3
1
ρ
10
10
10
10
The parameters used in the fractional Voigt model simulations to create the 2D spacetime plots shown in Figure 8.
Elastic
Voigt
Fractional Voigt
Heterogeneous Domain
Co
1.75
1.75
1.75
1.75
C1
1.75
1.75
1.75
1.75
a
1
0
0
1 (x < 1/2), 0 (x ≥ 1/2)
ν
1
1
0.3
1
ρ
10
10
10
10
3D numerical simulation
We now present a numerical simulation for viscoelastic waves propagating in the parallelepiped Ω = (0, 1) × (0, 10) × (0, 1) with a Dirichlet boundary on one of the small faces ΓD := (0, 1) × {0} × (0, 1). The PDE we are simulating is
where w(t) represents an enforced displacement at the Dirichlet boundary that takes value 0 in [0, 0.5], 1 in [1, 50] and transitions smoothly from 0 to 1 in the interval [0.5, 1]. Initial conditions are set to zero. The material is isotropic and locally homogeneous, with a strain-stress law given by
Therefore when y < 5, the model is purely elastic and when y ≥ 5 the model is a Zener viscoelastic model. For the discretization in space we use 𝓟2 finite elements on a mesh of 30, 720 tetrahedra obtained by partitioning a uniform quadrilateral mesh of 8 × 80 × 8 elements. Discretization in time is done by TRCQ using 500 time-steps over the interval [0, 50].
Figures 9 and 10 show snapshots of the displacement from the simulation, where the coloring in the first one exhibits the norm of the stress (averaged on each tetrahedron). We remark that in both figures, while we choose the same snapshots, the snapshots are not uniform in time. This is so we can highlight some of the more interesting aspects of the simulation which occur earlier in the time interval.
In the top rows of these figures, we observe the elastic waves, generated by the sudden deformation at y = 0, propagating along the rod in the y direction. We also note that the elastic part of the rod (y < 5) responds quickly to the sudden deformation and adjusts to the new displacement (enforced by the Dirichlet boundary condition) in about 70 time-steps, while the viscoelastic part of the rod (y > 5) shows a much slower response to the change of the displacement taking around 400 time-steps to adjust.
Conclusions
We have presented a framework for the analysis of a large class of models for viscoelastic wave propagation, including the classical models of Maxwell, Zener, and Voigt, fractional derivative versions of all of them, and combinations of different models in different subdomains. We have also shown that the general techniques using Laplace transforms can be improved in the classical models using techniques from evolutionary equations. Among the goals of future research is the analysis of fully discrete approximations of model equations in this general framework and the development of strategies for model fitting.