Analysis of models for viscoelastic wave propagation

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{matrix} ρ \ddot{u} = d i v σ + f \end{matrix}$ $$\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{matrix} \frac{1}{2} \end{matrix}$ $\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{matrix} σ (t) = \int_{0}^{t} D (t - τ) \dot{ε} (τ) d τ, \end{matrix}$ $$\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{matrix} σ + a \dot{σ} = C_{0} ε + C_{1} \dot{ε}, (Zener) \\ σ = C_{0} ε + C_{1} \dot{ε}, (Voigt) \\ σ + a \dot{σ} = C_{1} \dot{ε}, (Maxwell) \\ σ = C_{1} \dot{ε}, (Newton) \\ σ = C_{0} ε . (Linear elasticity) \end{matrix}$ $$\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 C₀ and C₁ 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 σ = C₁ε̇ 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 σ = C₀ε+C₁ε̇ 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{matrix} A : B = \sum_{i, j = 1}^{d} a_{i j} b_{i j}, A, B \in C^{d \times d}, \end{matrix}$ $$\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{matrix} C_{sym}^{d \times d} \end{matrix}$ $\begin{array}{} \mathbb C^{d\times d}_\text{sym} \end{array} $ := {A ∈ ℂ^{d × d} : A^⊤ = A}, and B ∈ $\begin{matrix} C_{skw}^{d \times d} \end{matrix}$ $\begin{array}{} \mathbb C^{d\times d}_\text{skw} \end{array} $ := {B ∈ ℂ^{d × d} : B^⊤ = −B}. We will write ∥M∥² = 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{matrix} \bar{Γ_{D}} \cup \bar{Γ_{N}} = Γ, Γ_{D} \cap Γ_{N} = \emptyset . \end{matrix}$ $$\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{matrix} L^{2} (Ω), L^{2} (Ω) := L^{2} (Ω; R^{d}), L^{2} (Ω) := L^{2} (Ω; R_{sym}^{d \times d}) \end{matrix}$ $$\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{matrix} (u, v)_{Ω} := \int_{Ω} u v, (u, v)_{Ω} := \int_{Ω} u \cdot v, (S, T)_{Ω} := \int_{Ω} S : T, \end{matrix}$ $$\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{matrix} H^{1} (Ω) = H^{1} (Ω; R^{d}) := {u : Ω \to R^{d} : u \in L^{2} (Ω), \nabla u \in L^{2} (Ω; R^{d \times d})}, \end{matrix}$ $$\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{matrix} ∥ u ∥_{1, Ω}^{2} := ∥ u ∥_{Ω}^{2} + ∥ \nabla u ∥_{Ω}^{2}, \end{matrix}$ $$\begin{array}{} \| \mathbf u\|_{1,\Omega}^2:=\| \mathbf u\|_\Omega^2+\|\nabla \mathbf u\|_\Omega^2, \end{array} $$

and the symmetric gradient operator $\begin{matrix} H^{1} (Ω) ∋ u ⟼ ε (u) := \frac{1}{2} (\nabla u + (\nabla u)^{⊤}) \in L^{2} (Ω) . \end{matrix}$ $$\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 γ : H¹(Ω) → H^1/2(Γ). For simplicity, we will write γ_Du := γu|_{Γ_D}. We then consider the Sobolev spaces: $\begin{matrix} H_{D}^{1} (Ω) := {w \in H^{1} (Ω) : γ_{D} w = 0} = \ker γ_{D}, \\ H^{1 / 2} (Γ_{D}) := {γ_{D} u : u \in H^{1} (Ω)} = r a n g e γ_{D}, \\ {\tilde{H}}^{1 / 2} (Γ_{N}) := {γ w |_{Γ_{N}} : w \in H_{D}^{1} (Ω)}, \\ H^{- 1 / 2} (Γ_{N}) := {\tilde{H}}^{1 / 2} (Γ_{N})^{'} . \end{matrix}$ $$\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{matrix} \tilde{H}^{1 / 2} (Γ_{N}) \subset L^{2} (Γ_{N}; R^{d}) \subset H^{- 1 / 2} (Γ_{N}) \end{matrix}$ $$\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{matrix} H (d i v, Ω) := {S \in L^{2} (Ω) : d i v S \in L^{2} (Ω)} . \end{matrix}$ $$\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{matrix} 〈 γ_{N} S, γ u 〉_{Γ_{N}} & = (S, \nabla u)_{Ω} + (d i v S, u)_{Ω} \\ = (S, ε (u))_{Ω} + (d i v S, u)_{Ω} \forall S \in H (d i v, Ω) \forall u \in H_{D}^{1} (Ω), \end{matrix}$ $$\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, ∞) → H¹(Ω) and σ : [0, ∞) → ℍ(div, Ω) satisfying for all t ≥ 0 $\begin{matrix} ρ \ddot{u} (t) = d i v σ (t) + f (t), \end{matrix}$ $$\begin{array}{} \rho\, \ddot{\mathbf u}(t) =\mathrm{div}\,\boldsymbol\sigma(t)+\mathbf f(t), \end{array} $$(1a) $\begin{matrix} γ_{D} u (t) = α (t), \end{matrix}$ $$\begin{array}{} \gamma_D\mathbf u(t) =\boldsymbol\alpha(t), \end{array} $$(1b) $\begin{matrix} γ_{N} σ (t) = β (t) . \end{matrix}$ $$\begin{array}{} \gamma_N\boldsymbol\sigma(t) =\boldsymbol\beta(t). \end{array} $$(1c)

Here ρ ∈ L^∞(Ω), with ρ ≥ ρ₀ almost everywhere for some positive constant ρ₀, models the mass density in the solid, which at the initial time t = 0 is at rest on the reference configuration Ω $\begin{matrix} u (0) = 0, \dot{u} (0) = 0 . \end{matrix}$ $$\begin{array}{} \mathbf u(0)=\mathbf 0, \qquad\dot{\mathbf u}(0)=\mathbf 0. \end{array} $$(1d)

Upper dots are used to denote time derivatives. The data are functions $\begin{matrix} f : [0, \infty) \to L^{2} (Ω), α : [0, \infty) \to H^{1 / 2} (Γ_{D}), β : [0, \infty) \to H^{- 1 / 2} (Γ_{N}) . \end{matrix}$ $$\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{matrix} σ (t) = \int_{0}^{t} D (t - τ) ε (\dot{u} (τ)) d τ . \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol\sigma(t)=\int_0^t \mathcal D(t-\tau)\,\boldsymbol\varepsilon(\dot{\mathbf u}(\tau)) \mathrm d\tau. \end{array} $$(1e)

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{matrix} C (s) := s L {D} (s) . \end{matrix}$ $$\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{matrix} C_{+} := {s \in C : R e s > 0} . \end{matrix}$ $$\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{matrix} C : C_{+} \to B (C^{d \times d}; L^{\infty} (Ω; C^{d \times d})) \equiv L^{\infty} (Ω; C^{d \times d \times d \times d}) . \end{matrix}$ $$\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}).

Hypothesis 1

(Symmetry). Almost everywhere in Ω:

$\begin{matrix} C (\bar{s}) \bar{M} = \bar{C (s) M} & \forall s \in C_{+} & \forall M \in C^{d \times d}, \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C(\overline s)\overline{\mathbf{M}}=\overline{\mathbf C(s)\mathbf{M}} & \qquad & \forall s\in \mathbb C_+ & \quad & \forall \mathbf{M}\in \mathbb C^{d\times d}, \end{array}$$(2a)

$\begin{matrix} C (s) M \in C_{sym}^{d \times d} & \forall s \in C_{+} & \forall M \in C^{d \times d}, \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C(s)\mathbf{M} \in \mathbb C^{d\times d}_\text{sym} & & \forall s\in \mathbb C_+ & & \forall \mathbf{M}\in \mathbb C^{d\times d}, \end{array} $$(2b)

$\begin{matrix} C (s) M = C (s) (\frac{1}{2} (M + M^{⊤})) & \forall s \in C_{+} & \forall M \in C^{d \times d}, \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C(s)\mathbf{M} = \mathbf C(s) (\tfrac12(\mathbf{M}+\mathbf{M}^\top)) & & \forall s\in \mathbb C_+& & \forall \mathbf{M} \in \mathbb C^{d\times d}, \end{array}$$(2c)

$\begin{matrix} C (s) M : N = C (s) N : M & \forall s \in C_{+} & \forall M, N \in C^{d \times d} . \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C(s)\mathbf{M}:\mathrm N=\mathbf C(s)\mathrm N: \mathrm M & & \forall s\in \mathbb C_+ & & \forall \mathbf{M},\mathrm N\in \mathbb C^{d\times d}. \end{array}$$(2d)

Some easy observations: conditions (2b) and (2d) imply (2c); if $\begin{matrix} M \in C_{skw}^{d \times d} \end{matrix}$ $\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

$\begin{matrix} inf_{0 < x < 1} x^{- ℓ} ψ (x) > 0, ℓ > 0, \end{matrix}$ $$\begin{array}{c} \displaystyle \inf_{0 \lt x \lt 1} x^{-\ell} \psi(x) \gt 0, \qquad \ell \gt 0, \end{array}$$(3a)

and such that almost everywhere in Ω

$\begin{matrix} R e (\bar{s} C (s) M : \bar{M}) \geq ψ (R e s) ∥ M ∥^{2} \forall M \in C_{sym}^{d \times d}, \forall s \in C_{+} . \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf{R}\text{e}(\overline{s}\mathbf {C}(s)\mathbf{M}:\overline{\mathbf{M}})\ge \psi(\mathbf{R}\text{e}~ s)\, \|\mathbf{M}\|^2 \qquad \forall \mathbf{M} \in \mathbb C^{d\times d}_\text{sym}, \quad \forall s\in \mathbb C_+. \end{array}$$(3b)

For each s ∈ ℂ₊, we take ∥C(s)∥ to be the smallest real number so that, almost everywhere in Ω,

$\begin{matrix} ∥ C (s) M ∥ \leq ∥ C (s) ∥ ∥ M ∥ \forall M \in C_{sym}^{d \times d} . \end{matrix}$ $$\begin{array}{c} \displaystyle \| \mathbf C(s)\mathbf{M}\|\le \|\mathbf C(s)\|\,\|\mathbf{M}\| \qquad \forall \mathbf{M}\in \mathbb C^{d\times d}_\text{sym}. \end{array}$$

Hypothesis 3

(Boundedness). There exists an integerr ≥ 0 and non-increasing functionϕ : (0, ∞) → (0, ∞) such that

$\begin{matrix} sup_{0 < x < 1} x^{k} ϕ (x) < \infty, k \geq 0, \end{matrix}$ $$\begin{array}{c} \displaystyle \sup_{0 \lt x \lt 1} x^k\phi(x) \lt \infty, \qquad k\ge 0, \end{array}$$(4a)

and almost everywhere in Ω

$\begin{matrix} ∥ C (s) ∥ \leq | s |^{r} ϕ (R e s) \forall s \in C_{+} . \end{matrix}$ $$\begin{array}{c} \displaystyle \| \mathbf C(s)\|\le |s|^r \phi(\mathbf{R}\text{e}~ s) \qquad \forall s\in \mathbb C_+. \end{array}$$(4b)

We can equivalently introduce this material tensor with a collection of holomorphic functions (the material coefficients)

$\begin{matrix} C_{i j k l} : C_{+} \to L^{\infty} (Ω; C), i, j, k, l = 1, \dots, d, \end{matrix}$ $$\begin{array}{c} \displaystyle C_{ijkl}:\mathbb C_+\to L^\infty(\Omega;\mathbb C), \qquad i,j,k,l=1,\ldots,d, \end{array}$$

satisfying

$\begin{matrix} C_{i j k l} (\bar{s}) = \bar{C_{i j k l} (s)} & \forall s \in C_{+} & i, j, k, l = 1, \dots, d, \end{matrix}$ $$\begin{array}{c} \displaystyle C_{ijkl}(\overline s)=\overline{C_{ijkl}(s)} & \qquad & \forall s\in \mathbb C_+ & \quad & i,j,k,l=1,\ldots,d, \end{array}$$(5a)

$\begin{matrix} C_{i j k l} (s) = C_{j i k l} (s) = C_{i j l k} (s) = C_{k l i j} (s) & \forall s \in C_{+} & i, j, k, l = 1, \dots, d . \end{matrix}$ $$\begin{array}{c} \displaystyle C_{ijkl}(s)=C_{jikl}(s)=C_{ijlk}(s)=C_{klij}(s) & \qquad & \forall s\in \mathbb C_+ & \quad & i,j,k,l=1,\ldots,d. \end{array}$$(5b)

In this case, we just define

$\begin{matrix} (C (s) M)_{i j} = \sum_{k, l = 1}^{d} C_{i j k l} (s) m_{k l} i, j = 1, \dots, d, \end{matrix}$ $$\begin{array}{c} \displaystyle (\mathbf C(s)\mathbf{M})_{ij}=\sum\limits_{k,l=1}^d C_{ijkl}(s)m_{kl} \qquad i,j=1,\ldots, d, \end{array}$$

and notice that (2) is equivalent to (5). When we write the material tensor in terms of coefficients, we can take

$\begin{matrix} ∥ C (s) ∥_{max} := d^{2} max_{i, j, k, l} ∥ C_{i j k l} (s) ∥_{L^{\infty} (Ω)} \end{matrix}$ $$\begin{array}{c} \displaystyle \| \mathbf C(s)\|_{\max}:=d^2 \,\max_{i,j,k,l} \|C_{ijkl}(s)\|_{L^\infty(\Omega)} \end{array}$$

as an upper bound of ∥C(s)∥ almost everywhere. With this point of view C(s) can be considered as an element of L^∞(Ω;ℂ^d×d×d×d).

The general viscoelastic material law in the Laplace domain is

$\begin{matrix} σ (u) = C (s) ε (u) . \end{matrix}$ $$\begin{array}{c} \displaystyle \boldsymbol{\sigma}(\mathbf u)=\mathbf C(s)\boldsymbol{\varepsilon}(\mathbf u). \end{array}$$(6)

Here and in the sequel, we identify C(s) with its associated ‘multiplication’ operator $\begin{matrix} C (s) \in B (L^{2} (Ω; C_{sym}^{d \times d})) \end{matrix}$ $\begin{array}{c} \displaystyle \mathbf C(s)\in \mathcal{B}(L^2(\Omega;\mathbb C^{d\times d}_\text{sym})) \end{array}$ , noticing that

$\begin{matrix} ∥ C (s) ∥_{L^{2} \to L^{2}} \leq e s s sup ∥ C (s) ∥ \leq ∥ C (s) ∥_{max} \forall s \in C_{+} . \end{matrix}$ $$\begin{array}{c} \displaystyle \|\mathbf C(s)\|_{L^2\to L^2} \le \mathrm{ess}\sup \|\mathbf C(s)\|\le \|\mathbf C(s)\|_{\max} \qquad \forall s\in \mathbb C_+. \end{array}$$(7)

The associated bilinear form is

$\begin{matrix} a (u, w; s) := (C (s) ε (u), ε (w))_{Ω} . \end{matrix}$ $$\begin{array}{c} \displaystyle a(\mathbf u,\mathbf w;s):=(\mathbf C(s)\boldsymbol{\varepsilon}(\mathbf {u}),\boldsymbol{\varepsilon}(\mathbf w))_\Omega. \end{array}$$

It is clear that a : H¹(Ω;ℂ) × H¹(Ω;ℂ) → ℂ is bilinear, bounded, symmetric, and satisfies

$\begin{matrix} | a (u, w; s) | \leq | s |^{r} ϕ (R e s) ∥ ε (u) ∥_{Ω} ∥ ε (w) ∥_{Ω} & \forall u, w \in H^{1} (Ω; C), s \in C_{+}, \end{matrix}$ $$\begin{array}{c} \displaystyle |a(\mathbf u,\mathbf w;s)| \le |s|^r \phi(\mathbf{R}\text{e}~ s) \,\|\boldsymbol{\varepsilon} (\mathbf{u})\|_{\Omega}\|\boldsymbol{\varepsilon} (\mathbf {w})\|_{\Omega} & \qquad &\forall \mathbf u,\mathbf w\in \mathbf H^1(\Omega;\mathbb C), \quad s\in \mathbb C_+, \end{array}$$(8)

$\begin{matrix} a (u, w; s) = (C (s) ε (u), \nabla w)_{Ω} & \forall u, w \in H^{1} (Ω; C), s \in C_{+}, \end{matrix}$ $$\begin{array}{c} \displaystyle a(\mathbf u,\mathbf w;s) =(\mathbf C(s)\boldsymbol{\varepsilon} (\mathbf{u}),\nabla \mathbf w)_\Omega & \qquad & \forall \mathbf u,\mathbf w\in \mathbf H^1(\Omega;\mathbb C), \quad s\in \mathbb C_+, \end{array}$$(9)

$\begin{matrix} R e a (u, \bar{s u}; s) \geq ψ (R e s) ∥ ε (u) ∥_{Ω}^{2} & \forall u \in H^{1} (Ω; C) s \in C_{+} . \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf{R}\text{e}~ a(\mathbf u,\overline{s\,\mathbf u};s) \ge \psi( \mathbf{R}\text{e}~ s)\|\boldsymbol{\varepsilon} (\mathbf{u})\|_\Omega^2 & \qquad & \forall \mathbf u\in \mathbf H^1(\Omega;\mathbb C)\quad s\in \mathbb C_+. \end{array}$$(10)

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

$\begin{matrix} C \in B (R^{d \times d}; L^{\infty} (Ω; R^{d \times d})) \equiv L^{\infty} (Ω; R^{d \times d \times d \times d}) \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C \in \mathcal B(\mathbb R^{d\times d};L^\infty(\Omega;\mathbb R^{d\times d}))\equiv L^\infty(\Omega;\mathbb R^{d\times d\times d\times d}) \end{array}$$(11a)

satisfy almost everywhere in Ω

$\begin{matrix} C M \in R_{sym}^{d \times d} \forall M \in R^{d \times d}, \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C \mathbf{M} \in \mathbb R^{d\times d}_\text{sym} \qquad \forall \mathbf{M}\in \mathbb R^{d\times d}, \end{array}$$(11b)

$\begin{matrix} C M : N = C N : M \forall M, N \in R^{d \times d}, \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C\mathbf{M}:\mathrm N=\mathbf C \mathrm N:\mathbf{M} \qquad\forall \mathbf{M},\mathrm N\in \mathbb R^{d\times d}, \end{array}$$(11c)

$\begin{matrix} C M : M \geq c ∥ M ∥^{2} \forall M \in R_{sym}^{d \times d}, \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C\mathbf{M}:\mathbf{M}\ge c\|\mathbf{M}\|^2 \qquad \forall \mathbf{M} \in \mathbb R^{d\times d}_\text{sym}, \end{array}$$(11d)

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

$\begin{matrix} C (M_{r e} + ı M_{i m}) := C M_{r e} + ı C M_{i m} . \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C(\mathbf{M}_{\mathrm{re}}+\imath \mathbf{M}_{\mathrm{im}}) :=\mathbf C\mathbf{M}_{\mathrm{re}}+\imath \mathbf C\mathbf{M}_{\mathrm{im}}. \end{array}$$

When hypotheses (11) are satisfied with c = 0, we will call C a non-negative Hookean model. For a steady Hookean model C, we will write

$\begin{matrix} ∥ C ∥_{max} := d^{2} max_{i, j, k, l} ∥ C_{i j k l} ∥_{L^{\infty} (Ω)} . \end{matrix}$ $$\begin{array}{c} \displaystyle \|\mathbf C\|_{\max}:=d^2\max_{i,j,k,l}\|C_{ijkl}\|_{L^\infty(\Omega)}. \end{array}$$

Lemma 1

LetCbe a steady Hookean model. Then:

If $\begin{matrix} N \in C_{skw}^{d \times d} \end{matrix}$ $\begin{array}{c} \displaystyle \mathrm N\in \mathbb C^{d\times d}_\text{skw} \end{array}$ , thenCN = 0 almost everywhere.

Almost everywhere in Ω

$\begin{matrix} C M = C (\frac{1}{2} (M + M^{⊤})) \forall M \in C^{d \times d} . \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C\mathbf{M} =\mathbf C (\tfrac12(\mathbf{M}+\mathbf{M}^\top)) \qquad \forall \mathbf{M}\in \mathbb C^{d\times d}. \end{array}$$

Almost everywhere in Ω

$\begin{matrix} C M : \bar{M} \geq c M : \bar{M} = c (∥ M_{r e} ∥^{2} + ∥ M_{i m} ∥^{2}) \forall M \in C_{sym}^{d \times d} . \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C\mathbf{M}:\overline{\mathbf{M}} \ge c \,\mathbf{M}:\overline{\mathbf{M}} = c(\|\mathbf{M}_{\mathrm{re}}\|^2+\|\mathbf{M}_{\mathrm{im}}\|^2) \qquad\forall \mathbf{M}\in \mathbb C^{d\times d}_\text{sym}. \end{array}$$

Ifa ∈ L^∞(Ω) satisfiesa ≥ a₀ > 0 almost everywhere, thenaCis a steady Hookean model.

Proof

Note first that $\begin{matrix} C M \in C_{sym}^{d \times d} \end{matrix}$ $\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{matrix} N \in C_{skw}^{d \times d} \end{matrix}$ $\begin{array}{c} \displaystyle \mathrm N\in \mathbb C^{d\times d}_\text{skw} \end{array}$ , then

$\begin{matrix} 0 = C M : N = M : C N \forall M \in C^{d \times d}, \end{matrix}$ $$\begin{array}{c} \displaystyle 0=\mathbf C\mathbf{M} :\mathrm N=\mathbf{M}: \mathbf C\mathrm N \qquad \forall \mathbf{M}\in \mathbb C^{d\times d}, \end{array}$$

which implies CN = 0. Property (b) follows from (a). Properties (c) and (d) are straightforward.

4.1

Elastic models

In the case where we take a Hookean model C₀ and we consider the constant function C(s) ≡ C₀, it is simple to see that the Hypotheses 1-3 are satisfied with ψ(x) := c₀x (c₀ being the lower bound for the model C₀) and (4), r = 0, and ϕ(x) := ∥C₀∥. The time domain version of this model is the usual linear strain-stress relation

$\begin{matrix} σ (t) = C_{0} ε (u (t)) . \end{matrix}$ $$\begin{array}{c} \displaystyle \boldsymbol\sigma(t)=\mathbf C_0 \boldsymbol\varepsilon(\mathbf u(t)). \end{array}$$

4.2

Zener’s classical viscoelastic model

We now consider the material law

$\begin{matrix} C (s) = (1 + a s)^{- 1} (C_{0} + s C_{1}), \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C(s)=(1+a\,s)^{-1} (\mathbf C_0+ s\,\mathbf C_1), \end{array}$$

where a ∈ L^∞(Ω) is strictly positive, C₀ and C₁ are steady Hookean material models, with

$\begin{matrix} C_{diff} := C_{1} - a C_{0} \geq 0, \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C_\text{diff}:=\mathbf C_1-a\,\mathbf C_0 \ge 0, \end{array}$$

that is, almost everywhere

$\begin{matrix} C_{1} M : M \geq a C_{0} M : M \forall M \in R_{sym}^{d \times d} . \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C_1\mathbf{M}:\mathbf{M}\ge a\,\mathbf C_0\mathbf{M}:\mathbf{M} \qquad \forall \mathbf{M}\in \mathbb R^{d\times d}_\text{sym}. \end{array}$$

This hypothesis makes C_diff a non-strict Hookean model. As we will see in the direct time-domain analysis of Section 9, C_diff is the diffusive part of the elastic model, while C₀ acts as a base or ground elastic model. We will make use of the formula

$\begin{matrix} (1 + a s)^{- 1} (C_{0} + s C_{1}) = C_{0} + s (1 + a s)^{- 1} C_{diff} . \end{matrix}$ $$\begin{array}{c} \displaystyle (1+as)^{-1}(\mathbf C_0+s\mathbf C_1)=\mathbf C_0+s(1+as)^{-1}\mathbf C_\text{diff}. \end{array}$$

The Laplace domain stress-strain relation can be written in implicit form

$\begin{matrix} σ + a s σ = C_{0} ε (u) + s C_{1} ε (u), \end{matrix}$ $$\begin{array}{c} \displaystyle \boldsymbol\sigma+a\,s\,\boldsymbol\sigma=\mathbf C_0\boldsymbol{\varepsilon} (\mathbf{u})+ s \mathbf C_1\boldsymbol{\varepsilon} (\mathbf{u}), \end{array}$$

corresponding to the differential relation in the time domain

$\begin{matrix} σ (t) + a \dot{σ} (t) = C_{0} ε u (t) + C_{1} ε \dot{u (t)}, σ (0) = 0. \end{matrix}$ $$\begin{array}{c} \displaystyle \boldsymbol\sigma(t)+a\,\dot{\boldsymbol\sigma}(t)= \mathbf C_0\varepsilon{\mathbf u(t)}+\mathbf C_1\varepsilon\dot{\mathbf u(t)}, \qquad \boldsymbol\sigma(0)=0. \end{array}$$

Proposition 2

(Laplace domain properties of Zener’s model). LetC(s) be a viscoelastic Zener model.

Ifc₀ > 0 is the lower bound for the tensorC₀, then Hypothesis 2 is satisfied withψ(x) := c₀x.

Hypothesis 3 is satisfied withr = 0 and

$\begin{matrix} ϕ (x) := \frac{(1 + ∥ a ∥_{L^{\infty} (Ω)})}{a_{0}^{2}} (∥ C_{0} ∥_{max} + ∥ C_{1} ∥_{max}) \frac{1}{min {1, x}^{2}}, \end{matrix}$ $$\begin{array}{c} \displaystyle \phi(x):= \frac{(1+\|a\|_{L^\infty(\Omega)})}{a_0^2} (\|\mathbf C_0\|_{\max}+\|\mathbf C_1\|_{\max}) \frac1{\min\{1,x\}^2}, \end{array}$$(12)

wherea₀ = 1/∥a^–1∥_L^∞(Ω)is a lower bound fora.

Proof

To prove (a), note first that

$\begin{matrix} \bar{s} (1 + a s)^{- 1} (C_{0} + s C_{1}) = \bar{s} C_{0} + (1 + a s)^{- 1} | s |^{2} C_{diff}, \end{matrix}$ $$\begin{array}{c} \displaystyle \overline s (1+as)^{-1} (\mathbf C_0+s\mathbf C_1) =\overline s \mathbf C_0+(1+as)^{-1} |s|^2 \mathbf C_\text{diff}, \end{array}$$

and therefore

$\begin{matrix} \bar{s} (1 + a s)^{- 1} (C_{0} + s C_{1}) M : \bar{M} = (C_{0} M : \bar{M}) \bar{s} + (| s |^{2} C_{diff} M : \bar{M}) (1 + a s)^{- 1}, \end{matrix}$ $$\begin{array}{c} \displaystyle \overline s (1+as)^{-1} (\mathbf C_0+s\mathbf C_1)\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)^{-1}, \end{array}$$

where all the bracketed quantities in the right hand side are real. Taking real parts and noticing that

$\begin{matrix} R e (1 + a s)^{- 1} = \frac{1}{| 1 + a s |^{2}} (1 + a R e s) \geq 0, \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf {R}\text{e} (1+as)^{-1}=\frac{1}{|1+as|^2} (1+a\mathbf{R}\text{e}~ s)\ge 0, \end{array}$$

the result follows.

To prove (b), we start with the explicit form of the coefficients

$\begin{matrix} C_{i j k l} (s) = (1 + a s)^{- 1} (C_{i j k l}^{0} + s C_{i j k l}^{1}) . \end{matrix}$ $$\begin{array}{c} \displaystyle C_{ijkl}(s)=(1+as)^{-1} (C_{ijkl}^0+ sC_{ijkl}^1). \end{array}$$

Let then $\begin{matrix} g_{0} = C_{i j k l}^{0} and g_{1} = C_{i j k l}^{1} \end{matrix}$ $\begin{array}{c} \displaystyle g_0=C_{ijkl}^0~~ \text{and}~~ g_1=C_{ijkl}^1 \end{array}$ . An easy computation shows that

$\begin{matrix} ∥ (1 + a s)^{- 1} ∥_{L^{\infty} (Ω)} \leq \frac{1 + | s | ∥ a ∥_{L^{\infty} (Ω)}}{a_{0}^{2} | s |^{2}}, \end{matrix}$ $$\begin{array}{c} \displaystyle \|(1+as)^{-1}\|_{L^\infty(\Omega)} \le \frac{1+|s|\, \|a\|_{L^\infty(\Omega)}}{a_0^2|s|^2}, \end{array}$$

where a₀ = 1/∥a^–1∥_L^∞(Ω) so that a ≥ a₀ almost everywhere. Using

$\begin{matrix} min {1, R e s} max {1, | s |} \leq | s | \forall s \in C_{+}, \end{matrix}$ $$\begin{array}{c} \displaystyle \min\{1,\mathbf{R}\text{e}~ s\}\max\{1,|s|\}\le |s| \qquad \forall s\in \mathbb C_+, \end{array}$$(13)

we easily estimate

$\begin{matrix} ∥ (1 + a s)^{- 1} (g_{0} + s g_{1}) ∥_{L^{\infty} (Ω)} \leq \frac{1 + | s | ∥ a ∥_{L^{\infty} (Ω)}}{a_{0}^{2} | s |^{2}} (∥ g_{0} ∥_{L^{\infty} (Ω)} + | s | ∥ g_{1} ∥_{L^{\infty} (Ω)}) \\ \leq \frac{1 + ∥ a ∥_{L^{\infty} (Ω)}}{a_{0}^{2}} (∥ g_{0} ∥_{L^{\infty} (Ω)} + ∥ g_{1} ∥_{L^{\infty} (Ω)}) \frac{max {1, | s |}^{2}}{| s |^{2}} \\ \leq \frac{1 + ∥ a ∥_{L^{\infty} (Ω)}}{a_{0}^{2}} (∥ g_{0} ∥_{L^{\infty} (Ω)} + ∥ g_{1} ∥_{L^{\infty} (Ω)}) \frac{1}{min {1, R e s}^{2}}, \end{matrix}$ $$\begin{array}{} \displaystyle \| (1+as)^{-1}(g_0+s g_1)\|_{L^\infty(\Omega)} \le \frac{1+|s|\, \|a\|_{L^\infty(\Omega)}}{a_0^2|s|^2} (\| g_0\|_{L^\infty(\Omega)}+|s|\, \|g_1\|_{L^\infty(\Omega)}) \\ \displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\, \le \frac{1+ \|a\|_{L^\infty(\Omega)}}{a_0^2} (\| g_0\|_{L^\infty(\Omega)}+ \|g_1\|_{L^\infty(\Omega)}) \frac{\max\{1,|s|\}^2}{|s|^2} \\ \displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\, \le \frac{1+ \|a\|_{L^\infty(\Omega)}}{a_0^2} (\| g_0\|_{L^\infty(\Omega)}+ \|g_1\|_{L^\infty(\Omega)}) \frac1{\min\{1,\mathbf{R}\text{e}~ s\}^2}, \end{array}$$

which proves the result.

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 ∈ ℂ₊:

$\begin{matrix} λ (\bar{s}) = \bar{λ (s)}, μ (\bar{s}) = \bar{μ (s)}, \end{matrix}$ $$\begin{array}{c} \displaystyle \lambda(\overline s)=\overline{\lambda(s)}, \quad\qquad \mu(\overline s)=\overline{\mu(s)}, \end{array}$$(14a)

$\begin{matrix} R e (\bar{s} λ (s)) \geq 0, R e (\bar{s} μ (s)) \geq μ_{0} R e s (μ_{0} > 0) . \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf {R}\text {e}(\overline s\, \lambda(s))\ge 0, \qquad \mathbf {R}\text {e}(\overline s\,\mu(s))\ge \mu_0 \mathbf{R}\text{e}~ s \qquad (\mu_0 \gt 0). \end{array}$$(14b)

We then define

$\begin{matrix} C (s) M := 2 μ (s) (\frac{1}{2} (M + M^{⊤})) + λ (s) (t r M) I, \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C(s)\mathbf{M}:=2\mu(s) (\tfrac12(\mathbf{M}+\mathbf{M}^\top))+\lambda(s) \,(\mathrm{tr}\,\mathrm{M})\, \text{I}, \end{array}$$

where I is the d × d identity matrix, so that the material law is

$\begin{matrix} σ = 2 μ (s) ε (u) + λ (s) (\nabla \cdot u) I . \end{matrix}$ $$\begin{array}{c} \displaystyle \boldsymbol\sigma=2\mu(s) \boldsymbol{\varepsilon} (\mathbf{u})+\lambda(s)(\nabla\cdot\mathbf u)\,\mathrm{I}. \end{array}$$

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

$\begin{matrix} a m_{μ} \leq b_{μ}, a m_{λ} \leq b_{λ} . \end{matrix}$ $$\begin{array}{c} a\,m_\mu\le b_\mu, \qquad a\, m_\lambda\le b_\lambda. \end{array}$$

Then

$\begin{matrix} λ (s) := \frac{m_{λ} + b_{λ} s}{1 + a s} and μ (s) := \frac{m_{μ} + b_{μ} s}{1 + a s} \end{matrix}$ $$\begin{array}{c} \displaystyle \lambda(s):=\frac{m_\lambda+b_\lambda s}{1+a\,s} \qquad \text{and} \qquad \mu(s):=\frac{m_\mu+b_\mu s}{1+a\,s} \end{array}$$

satisfy (14). To prove the lower bound (14b), note that

$\begin{matrix} R e (\bar{s} μ (s)) = R e (\frac{\bar{s} m_{μ} (1 + a s) + \bar{s} s (b_{μ} - a m_{μ})}{1 + a s}) \\ = m_{μ} R e s + | s |^{2} (b_{μ} - a m_{μ}) \frac{1 + a R e s}{| 1 + a s |^{2}} \geq m_{μ} R e s . \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf {R}\text{e}(\overline s\mu(s)) = \mathbf {R}\text{e} \left( \frac{\overline sm_\mu(1 + as) + \overline s s(b_\mu - am_\mu)}{1+as} \right) \\ \displaystyle\qquad\qquad~ = m_\mu \mathbf{R}\text{e}~ s+ |s|^2 (b_\mu-a m_\mu) \frac{1 + a \mathbf{R}\text{e}~ s}{| 1+ a s|^2} \ge m_\mu \mathbf{R}\text{e}~ s. \end{array}$$

4.3

Fractional Zener models

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

$\begin{matrix} (1 + a s^{ν})^{- 1} (C_{0} + s^{ν} C_{1}), \end{matrix}$ $$\begin{array}{c} \displaystyle (1+a\,s^\nu)^{-1}(\mathbf C_0+s^\nu\mathbf C_1), \end{array}$$

where a, C₀, and C₁ satisfy the same hypotheses as in Section 4.2. In the time domain, this corresponds to

$\begin{matrix} σ (t) + a \partial^{ν} σ (t) = C_{0} ε (u (t)) + C_{1} ε \partial^{ν} (u (t)), \end{matrix}$ $$\begin{array}{c} \displaystyle \boldsymbol\sigma(t)+a\,\partial^\nu\boldsymbol\sigma(t) =\mathbf C_0\boldsymbol{\varepsilon}({\mathbf u(t))}+\mathbf C_1\boldsymbol{\varepsilon}{\partial^\nu(\mathbf u(t))}, \end{array}$$

where

$\begin{matrix} (\partial^{ν} f) (t) = \frac{1}{Γ (1 - ν)} \int_{0}^{t} \frac{\dot{f} (τ)}{(t - τ)^{ν}} d τ \end{matrix}$ $$\begin{array}{c} \displaystyle (\partial^\nu f)(t)= \frac1{\Gamma(1-\nu)}\int_0^t \frac{\dot f(\tau)}{(t-\tau)^{\nu}} \mathrm d\tau \end{array}$$

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).

Ifc₀ > 0 is the lower bound for the tensorC₀, thenC(s^ν) satisfies Hypothesis 2 withψ(x) := c₀x.

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{matrix} \bar{s} C (s^{ν}) M : \bar{M} = (C_{0} M : \bar{M}) \bar{s} + (| s |^{2} C_{diff} M : \bar{M}) (1 + a s^{ν})^{- 1} s^{ν - 1} . \end{matrix}$ $$\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}$$

We thus only need to show

$\begin{matrix} R e (1 + a s^{ν})^{- 1} s^{ν - 1} \geq 0 \forall s \in C_{+} . \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf{R}\text{e} (1+as^\nu)^{-1} s^{\nu-1} \ge 0 \qquad \forall s \in \mathbb C_+. \end{array}$$

To see that, first observe

$\begin{matrix} (1 + a s^{ν})^{- 1} s^{ν - 1} = \frac{1}{s^{1 - ν} + a s}, \end{matrix}$ $$\begin{array}{c} \displaystyle (1+as^\nu)^{-1} s^{\nu-1} = \frac{1}{s^{1-\nu} + as}, \end{array}$$

where 1 > 1 – ν > 0. Since s ∈ ℂ₊, we have

$\begin{matrix} R e (s^{1 - ν} + a s) = R e s^{1 - ν} + R e a s \geq a_{0} R e s > 0, \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf{R}\text{e} (s^{1-\nu} +as)=\mathbf{R}\text{e}~ s^{1-\nu}+ \mathbf{R}\text{e}~ as \geq a_0 \mathbf{R}\text{e}~ s \gt 0, \end{array}$$

which proves the result. The proof of (b) is a direct consequence of Proposition 2(b) and Lemma 4.

Lemma 4

The following inequality holds:

$\begin{matrix} min {1, R e s} \leq R e s^{ν}, \forall ν \in (0, 1), s \in C_{+} . \end{matrix}$ $$\begin{array}{c} \displaystyle \min\{1,\mathbf{R}\text{e}~ s\}\le \mathbf{R}\text{e}~ s^\nu, \qquad \forall \nu\in (0,1), \quad s\in \mathbb C_+. \end{array}$$

Proof

Writing s = r, e^θ with r > 0 and θ ∈ (−π/2, π/2), it is clear that an equivalent form of the result is the inequality

$\begin{matrix} min {1, r \cos θ} \leq r^{ν} \cos (ν θ) r > 0, θ \in (- π / 2, π / 2) . \end{matrix}$ $$\begin{array}{c} \displaystyle \min\{1,r\,\cos\theta\}\le r^\nu \cos(\nu\theta) \qquad r \gt 0, \qquad \theta\in(-\pi/2,\pi/2). \end{array}$$(15)

It is also clear that we only need to prove (15) for θ ∈ [0, π/2). Fix then θ ∈ [0, π/2) and consider the function

$\begin{matrix} f_{θ} (ν) := \cos (ν θ) - (\cos θ)^{ν} . \end{matrix}$ $$\begin{array}{c} \displaystyle f_\theta(\nu):=\cos(\nu\theta)-(\cos\theta)^\nu. \end{array}$$

We have

$\begin{matrix} f_{θ} (0) = 0, f_{θ} (1) = 0, f_{θ}^{″} (ν) = - θ^{2} \cos (ν θ) - (\cos θ)^{ν} \log^{2} (\cos θ) \leq 0, \end{matrix}$ $$\begin{array}{c} \displaystyle f_\theta(0)=0, \quad f_\theta(1)=0, \quad f_\theta''(\nu)= -\theta^2 \cos(\nu\theta) -(\cos\theta)^\nu\log^2(\cos\theta)\le 0, \end{array}$$

and therefore (by concavity) f_θ(ν) ≥ 0 for ν ∈ (0, 1) or equivalently

$\begin{matrix} \cos (ν θ) \geq (\cos θ)^{ν} ν \in (0, 1), θ \in [0, π / 2) . \end{matrix}$ $$\begin{array}{c} \displaystyle \cos(\nu\theta)\ge (\cos\theta)^\nu \qquad \nu\in (0,1), \quad \theta\in [0,\pi/2). \end{array}$$

Finally, this implies

$\begin{matrix} r^{ν} \cos (ν θ) \geq (r \cos θ)^{ν} \geq (\begin{matrix} 1, & if r \cos θ \geq 1, \\ r \cos θ, & if r \cos θ < 1, \end{matrix}) \end{matrix}$ $$\begin{array}{c} \displaystyle r^\nu\cos(\nu\theta)\ge (r\cos\theta)^\nu \ge \begin{cases} 1, & \mbox{if}~~ r\cos\theta\ge 1,\\ r\cos\theta, & \mbox{if} ~~r\cos\theta \lt 1, \end{cases} \end{array}$$

which proves (15) and hence the result.

4.4

Maxwell’s model

Maxwell’s model is given by

$\begin{matrix} C (s) = (1 + a s)^{- 1} s C_{1}, \end{matrix}$ $$\begin{array}{c} \displaystyle \mathbf C(s)=(1+as)^{-1}s\mathbf C_1, \end{array}$$

where C₁ is a steady Hookean model (with lower bound c₁ > 0) and a ∈ L^∞(Ω) satisfies a ≥ a₀ > 0 almost everywhere, for some constant a₀. In the time domain this gives again an implicit strain-to-stress relation

$\begin{matrix} σ (t) + a \dot{σ} (t) = C_{1} ε (\dot{u} (t)) . \end{matrix}$ $$\begin{array}{c} \displaystyle \boldsymbol\sigma(t)+a\dot{\boldsymbol\sigma}(t)=\mathbf C_1\boldsymbol{\varepsilon}({\dot{\mathbf u}(t))}. \end{array}$$

Proposition 5

Maxwell’s model satisfies Hypotheses 1-3 withr = 0 and

$\begin{matrix} ψ (x) := \frac{c_{1} min {1, a_{0}^{3}}}{2 ∥ a ∥_{L^{\infty} (Ω)}^{2}} min {1, x^{3}}, ϕ (x) := \frac{(1 + ∥ a ∥_{L^{\infty} (Ω)})}{a_{0}^{2}} ∥ C_{1} ∥_{max} \frac{1}{min {1, x^{2}}} . \end{matrix}$ $$\begin{array}{c} \displaystyle \psi(x):=\frac{c_1\min\{1,a_0^3\}}{2\| a\|_{L^\infty(\Omega)}^2} \min\{1,x^3\}, \quad \phi(x):=\frac{(1+\|a\|_{L^\infty(\Omega)})}{a_0^2} \|\mathbf C_1\|_{\max} \frac1{\min\{1,x^2\}}. \end{array}$$

Proof

Hypothesis 1 is easy to verify. Hypothesis 3 can be verified using the proof of Proposition 2 taking C₀ = 0. To prove Hypothesis 2 note that almost everywhere

$\begin{matrix} R e (\bar{s} C (s) M : \bar{M}) = (C_{1} M : \bar{M}) | s |^{2} R e (1 + a s)^{- 1} \\ \geq c_{1} ∥ M ∥^{2} | s |^{2} R e (1 + a s)^{- 1} \forall M \in C_{sym}^{d \times d}, s \in C_{+} . \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf{R}\text{e} (\overline s\mathbf C(s)\mathbf{M}:\overline{\mathbf{M}}) =(\mathbf C_1\mathbf{M}:\overline{\mathbf{M}}) |s|^2 \mathbf{R}\text{e} (1+as)^{-1}\\ \displaystyle\qquad\qquad\qquad\qquad \ge c_1 \|\mathbf{M}\|^2|s|^2 \mathbf{R}\text{e} (1+as)^{-1} \qquad\forall \mathbf{M}\in \mathbb C^{d\times d}_\text{sym}, \quad s\in \mathbb C_+. \end{array}$$

Note now that since

$\begin{matrix} \frac{x^{2}}{1 + x^{2}} \geq \frac{1}{2} min {1, x^{2}} \forall x > 0, \end{matrix}$ $$\begin{array}{c} \displaystyle \frac{x^2}{1+x^2}\ge \frac{1}{2}\min\{1,x^2\} \quad \forall x \gt 0, \end{array}$$(16)

then

$\begin{matrix} \frac{| a s |^{2}}{| 1 + a s |^{2}} \geq \frac{1}{2} \frac{| a s |^{2}}{1 + | a s |^{2}} \geq \frac{1}{4} min {1, | a s |^{2}} \\ \geq \frac{1}{4} min {1, (a_{0} R e s)^{2}} \geq \frac{1}{4} min {1, a_{0}^{2}} min {1, (R e s)^{2}}, \end{matrix}$ $$\begin{array}{} \displaystyle \frac{|as|^2}{|1+as|^2} \ge \frac{1}{2} \frac{|as|^2}{1+|as|^2} \ge\frac14\min\{1,|as|^2\}\\ \qquad\qquad \ge \frac14 \min\{ 1, (a_0\mathbf{R}\text{e}~ s)^2\} \ge \frac14\min\{1,a_0^2\} \min\{1, (\mathbf{R}\text{e}~ s)^2\}, \end{array}$$

and therefore, almost everywhere and for all s ∈ ℂ₊

$\begin{matrix} | s |^{2} R e (1 + a s)^{- 1} = \frac{1 + a R e s}{a^{2}} \frac{| a s |^{2}}{| 1 + a s |^{2}} \\ \geq \frac{min {1, a_{0}} (1 + R e s)}{∥ a ∥_{L^{\infty} (Ω)}^{2}} \frac{1}{4} min {1, a_{0}^{2}} min {1, (R e s)^{2}}, \end{matrix}$ $$\begin{array}{} \displaystyle |s|^2 \mathbf{R}\text{e}(1+as)^{-1} = \frac{1+a\mathbf{R}\text{e}~ s}{a^2} \,\frac{|as|^2}{|1+as|^2}\\ \qquad\qquad\qquad\quad~ \ge \frac{\min\{1,a_0\}(1+\mathbf{R}\text{e}~ s)}{\| a\|_{L^\infty(\Omega)}^2} \frac14\min\{1,a_0^2\} \min\{1, (\mathbf{R}\text{e}~ s)^2\}, \end{array}$$

which proves the result.

Proposition 6

(Fractional Maxwell’s model). Ifν ∈ (0, 1) andC(s) = (1 + as)⁻¹sC₁is a Maxwell model, thenC(s^ν) satisfies Hypotheses 1-3 withr = 0 and the functionsψandϕ of Proposition 5.

Proof

Hypothesis 1 is straightforward and Hypothesis 3 follows from the fact that

$\begin{matrix} \frac{1}{min {1, R e s^{ν}}^{2}} \leq \frac{1}{min {1, R e s}^{2}} s \in C_{+}, ν \in (0, 1), \end{matrix}$ $$\begin{array}{} \displaystyle \frac1{\min\{1,\mathbf{R}\text{e}~ s^\nu\}^2}\le \frac1{\min\{1,\mathbf{R}\text{e}~ s\}^2} \qquad s\in \mathbb C_+, \quad \nu\in (0,1), \end{array}$$

as follows from Lemma 4. To verify Hypothesis 2 we first estimate

$\begin{matrix} R e (\bar{s} C (s^{ν}) M : \bar{M}) \geq c_{1} ∥ M ∥^{2} | s |^{2} R e \frac{s^{ν - 1}}{1 + a s^{ν}} \\ \geq \frac{c_{1}}{∥ a ∥_{L^{\infty} (Ω)}^{2}} ∥ M ∥^{2} \frac{| a s |^{2}}{| s^{1 - ν} + a s |^{2}} R e (s^{1 - ν} + a s), \forall s \in C_{+}, \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf{R}\text{e}(\overline s\,\mathbf C(s^\nu)\mathbf{M}:\overline{\mathbf{M}}) \ge c_1 \|\mathbf{M}\|^2 |s|^2 \mathbf{R}\text{e} \frac{s^{\nu-1}}{1+as^\nu}\\ \displaystyle\qquad\qquad\qquad\qquad~\ge \frac{c_1}{\| a\|_{L^\infty(\Omega)}^2} \|\mathbf{M}\|^2 \frac{|as|^2}{|s^{1-\nu}+as|^2}\mathbf{R}\text{e}(s^{1-\nu}+as), \qquad \forall s\in \mathbb C_+, \end{array}$$

almost everywhere. Using (16) and Lemma 4, we can easily bound

$\begin{matrix} \frac{| a s |^{2}}{| s^{1 - ν} + a s |^{2}} \geq \frac{1}{2} \frac{| a s |^{2}}{| s |^{2 - 2 ν} + | a s |^{2}} \geq \frac{1}{4} min {1, \frac{| a s |^{2}}{| s |^{2 - 2 ν}}} \\ \geq \frac{1}{4} min {1, a_{0}^{2}} min {1, | s^{ν} |}^{2} \geq \frac{1}{4} min {1, a_{0}^{2}} min {1, (R e s)^{2}} . \end{matrix}$ $$\begin{array}{} \displaystyle \frac{|as|^2}{|s^{1-\nu}+as|^2} \ge \frac{1}{2}\,\frac{|as|^2}{|s|^{2-2\nu}+|as|^2} \ge \frac14 \min\{ 1,\frac{|as|^2}{|s|^{2-2\nu}}\}\\ \displaystyle\qquad\qquad\quad~ \ge \frac14\min\{1,a_0^2\} \min\{ 1,|s^\nu|\}^2 \ge \frac14\min\{1,a_0^2\} \min\{1,(\mathbf{R}\text{e}~ s)^2\}. \end{array}$$

At the same time

$\begin{matrix} R e (s^{1 - ν} + a s) \geq min {1, R e s} + a_{0} R e s \geq 2 min {1, a_{0}} min {1, R e s}, \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf{R}\text{e}(s^{1-\nu}+as) \ge \min\{1,\mathbf{R}\text{e}~ s\} + a_0 \mathbf{R}\text{e}~ s \ge 2\min\{1,a_0\}\min\{1,\mathbf{R}\text{e}~ s\}, \end{array}$$

(we have used Lemma 4 again) and the proof is finished.

4.5

Voigt’s model

Voigt’s model uses

$\begin{matrix} C (s) := C_{0} + s C_{1} \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf C(s):=\mathbf C_0+s\mathbf C_1 \end{array}$$

as a viscoelastic parameter model, where C₀ is a steady Hookean material model and C₁ is a non-negative Hookean model. In those parts of the domain where C₁ = 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

$\begin{matrix} σ (t) = C_{0} ε (u (t)) + C_{1} ε (\dot{u} (t)) . \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol\sigma(t)=\mathbf C_0\boldsymbol{\varepsilon}({\mathbf u(t))}+\mathbf C_1\boldsymbol{\varepsilon}({\dot{\mathbf u}(t))}. \end{array}$$

Note that this can be plugged into the momentum equation yielding

$\begin{matrix} ρ \ddot{u} (t) = d i v (C_{0} ε (u (t)) + C_{1} ε (\dot{u} (t))) + f (t), \end{matrix}$ $$\begin{array}{} \displaystyle \rho\,\ddot{\mathbf u}(t)= \mathrm{div}\,(\mathbf C_0\boldsymbol{\varepsilon}({\mathbf u(t))}+\mathbf C_1\boldsymbol{\varepsilon}({\dot{\mathbf u}(t))})+\mathbf f(t), \end{array}$$

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

$\begin{matrix} ψ (x) := c_{0} x, ϕ (x) := \frac{∥ C_{0} ∥ + ∥ C_{1} ∥}{min {1, x}}, \end{matrix}$ $$\begin{array}{} \displaystyle \psi(x):=c_0 x, \qquad \phi(x):=\frac{\| \mathbf C_0\|+\|\mathbf C_1\|}{\min\{1,x\}}, \end{array}$$

wherec₀is the lower bound ofC₀.

Proof

It is straightforward proof using the type of inequalities of the proof of Proposition 2.

Proposition 8

(Fractional Voigt’s model). Ifν ∈ (0, 1) andC(s) = C₀ + sC₁is a Voigt model, thenC(s^ν) satisfies Hypotheses 1-3 withr = 1,

$\begin{matrix} ψ (x) := c_{0} x, ϕ (x) := \frac{∥ C_{0} ∥ + ∥ C_{1} ∥}{min {1, x^{2}}} . \end{matrix}$ $$\begin{array}{} \displaystyle \psi(x):=c_0 x, \qquad \phi(x):=\frac{\| \mathbf C_0\|+\|\mathbf C_1\|}{\min\{1,x^2\}}. \end{array}$$

Proof

By Proposition 7

$\begin{matrix} ∥ C (s^{ν}) ∥ \leq | s |^{ν} \frac{∥ C_{0} ∥ + ∥ C_{1} ∥}{min {1, R e s^{ν}}} \leq \frac{| s |}{R e s^{1 - ν}} \frac{∥ C_{0} ∥ + ∥ C_{1} ∥}{min {1, R e s}}, \end{matrix}$ $$\begin{array}{} \displaystyle \| \mathbf C(s^\nu)\| \le |s|^\nu \frac{\| \mathbf C_0\|+\|\mathbf C_1\|}{\min\{1,\mathbf{R}\text{e}~ s^\nu\}} \le \frac{|s|}{\mathbf{R}\text{e}~ s^{1-\nu}}\frac{\| \mathbf C_0\|+\|\mathbf C_1\|}{\min\{1,\mathbf{R}\text{e}~ s\}}, \end{array}$$

where we have used Lemma 4. Using Lemma 4 again we obtain the upper bound for ∥C(s)∥ almost everywhere. For positivity (Hypothesis 2) note that

$\begin{matrix} R e (\bar{s} C (s^{ν}) M : \bar{M}) = (R e s) (C_{0} M : \bar{M}) + | s |^{2} (C_{1} M : \bar{M}) R e s^{ν - 1} \geq c_{0} R e s ∥ M ∥^{2}, \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf{R}\text{e} (\overline s\mathbf C(s^\nu)\mathbf{M}:\overline{\mathbf{M}}) =(\mathbf{R}\text{e}~ s) (\mathbf C_0\mathbf{M}:\overline{\mathbf{M}})+ |s|^2 (\mathbf C_1\mathbf{M}:\overline{\mathbf{M}}) \mathbf{R}\text{e}~ s^{\nu-1} \ge c_0 \mathbf{R}\text{e}~ s \|\mathbf{M}\|^2, \end{array}$$

almost everywhere.

4.6

Coupled models

Proposition 9

Let Ω₁, …, Ω_Jbe non-overlapping subdomains of Ω such that $\begin{matrix} \bar{Ω} = \cup_{j = 1}^{J} {\bar{Ω}}_{j} \end{matrix}$ $\begin{array}{} \displaystyle \overline\Omega=\cup_{j=1}^J \overline\Omega_j \end{array}$ . Assume that C_jis a viscoelastic model in the domain Ω_j, satisfying the Hypotheses 1-3. Then

$\begin{matrix} C (s) := \sum_{j = 1}^{J} χ_{Ω_{j}} C_{j} (s) \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm C(s):=\sum\limits_{j=1}^J \chi_{\Omega_j}\mathrm C_j(s) \end{array}$$(17)

defines a viscoelastic model in the full domain Ω.

Proof

Hypothesis 1 follows readily. Assume now that there exist non-decreasing ψ_j : (0, ∞) →(0, ∞) satisfying

$\begin{matrix} ψ_{j} (x) \geq c_{j} x^{ℓ_{j}} \forall x \in (0, 1], ℓ_{j} \geq 0, c_{j} > 0, \end{matrix}$ $$\begin{array}{} \displaystyle \psi_j(x)\ge c_j x^{\ell_j} \quad \forall x\in (0,1], \qquad \ell_j\ge 0, c_j \gt 0, \end{array}$$

and

$\begin{matrix} R e (\bar{s} C_{j} (s) M : \bar{M}) \geq ψ_{j} (R e s) ∥ M ∥^{2} a.e. in Ω_{j} \forall M \in C_{sym}^{d \times d}, s \in C_{+} . \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf{R}\text{e}(\overline s \mathbf C_j(s) \mathbf{M}:\overline{\mathbf{M}})\ge \psi_j(\mathbf{R}\text{e}~ s) \|\mathbf{M}\|^2 \qquad \mbox{a.e. in}~~ \Omega_j \quad \forall \mathbf{M} \in \mathbb C^{d\times d}_\text{sym}, \quad s\in \mathbb C_+. \end{array}$$

Let now ψ(x) := min{ψ₁(x), …, ψ_J(x)}. If we take c := min{c₁, …, c_J} and ℓ := max{ℓ₁, …, ℓ_J} it follows that ψ is non-decreasing,

$\begin{matrix} ψ (x) \geq c x^{ℓ} \forall x \in (0, 1], \end{matrix}$ $$\begin{array}{} \displaystyle \psi(x)\ge c\,x^\ell \quad\forall x\in (0,1], \end{array}$$

and

$\begin{matrix} R e (\bar{s} C (s) M : \bar{M}) \geq ψ (R e s) ∥ M ∥^{2} a.e. in Ω \forall M \in C_{sym}^{d \times d} . \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf{R}\text{e}(\overline s \mathbf C(s) \mathbf{M}:\overline{\mathbf{M}})\ge \psi(\mathbf{R}\text{e}~ s) \|\mathbf{M}\|^2 \qquad \mbox{a.e. in} ~~\Omega \quad \forall \mathbf{M} \in \mathbb C^{d\times d}_\text{sym}. \end{array}$$

For the upper bounds, consider integers r_j ≥ 0 and non-increasing functions ϕ_j : (0, ∞) → (0, ∞) such that

$\begin{matrix} ϕ_{j} (x) \leq d_{j} x^{- k_{j}} \forall x \in (0, 1], k_{j} \geq 0, d_{j} > 0, \end{matrix}$ $$\begin{array}{} \displaystyle \phi_j(x)\le d_j x^{-k_j} \quad \forall x\in (0,1],\qquad k_j\ge 0, d_j \gt 0, \end{array}$$

and

$\begin{matrix} ∥ C_{j} (s) M ∥ \leq | s |^{r_{j}} ϕ_{j} (R e s) ∥ M ∥ a.e. in Ω_{j} \forall M \in C_{sym}^{d \times d}, s \in C_{+} . \end{matrix}$ $$\begin{array}{} \displaystyle \| \mathbf C_j(s)\mathbf{M}\| \le |s|^{r_j}\phi_j(\mathbf{R}\text{e}~ s)\|\mathbf{M}\| \qquad \mbox{a.e. in}~~ \Omega_j \quad \forall \mathbf{M} \in \mathbb C^{d\times d}_\text{sym}, \quad s\in \mathbb C_+. \end{array}$$

Let then r := max{r₁, …, r_J} and the non-increasing function

$\begin{matrix} ϕ (x) := max {x^{r_{1} - r} ϕ_{1} (x), \dots, x^{r_{J} - r} ϕ_{J} (x)} . \end{matrix}$ $$\begin{array}{} \displaystyle \phi(x):=\max\{ x^{r_1-r}\phi_1(x),\ldots,x^{r_J-r}\phi_J(x)\}. \end{array}$$

If we take d := max{d₁, …, d_J} and k := max{k₁, …, k_J}, we have that

$\begin{matrix} ϕ (x) \leq d x^{- k} \forall x \in (0, 1], \end{matrix}$ $$\begin{array}{} \displaystyle \phi(x)\le d\,x^{-k} \qquad \forall x\in (0,1], \end{array}$$

and

$\begin{matrix} ∥ C (s) M ∥ \leq | s |^{r} ϕ (R e s) ∥ M ∥ a.e. in Ω \forall M \in C_{sym}^{d \times d}, s \in C_{+}, \end{matrix}$ $$\begin{array}{} \displaystyle \|\mathbf C(s)\mathbf{M}\| \le |s|^r \phi(\mathbf{R}\text{e}~ s)\|\mathbf{M}\| \qquad \mbox{a.e. in}~~ \Omega \quad \forall \mathbf{M} \in \mathbb C^{d\times d}_\text{sym}, \quad s\in \mathbb C_+, \end{array}$$

which finishes the proof.

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 H¹(Ω)) 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

$\begin{matrix} C (s) = s C_{1}, \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf C(s)=s\mathbf C_1, \end{array}$$

where C₁ is a steady Hookean model. Since the viscoelastic law in the time domain becomes

$\begin{matrix} σ (t) = C_{1} ε (\dot{u} (t)), \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol\sigma(t)=\mathbf C_1\boldsymbol{\varepsilon}({\dot{\mathbf u}(t))}, \end{array}$$

the model can be simplified to

$\begin{matrix} ρ \dot{u} (t) = d i v C_{1} ε (u (t)) + g (t) \end{matrix}$ $$\begin{array}{} \displaystyle \rho \dot{\mathbf u}(t)=\mathrm{div}\,\mathbf C_1\boldsymbol{\varepsilon}({\mathbf u(t))} +\mathbf g(t) \end{array}$$

(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:

$\begin{matrix} | | | u | | |_{c}^{2} := c^{2} ∥ ρ^{1 / 2} u ∥_{Ω}^{2} + ∥ ε (u) ∥_{Ω}^{2} . \end{matrix}$ $$\begin{array}{} \displaystyle |\!|\!|{\mathbf u}|\!|\!|_c^2:=c^2 \| \rho^{1/2}\mathbf u\|_\Omega^2+\|\boldsymbol{\varepsilon} (\mathbf{u})\|_\Omega^2.\end{array}$$

Recall that the mass density function is ρ ∈ L^∞(Ω) is strictly positive. By (13), it follows that

$\begin{matrix} min {1, R e s} | | | u | | |_{1} \leq | | | u | | |_{| s |} \leq \frac{| s |}{min {1, R e s}} | | | u | | |_{1} \forall u \in H^{1} (Ω; C) s \in C_{+} . \end{matrix}$ $$\begin{array}{} \displaystyle \min\{1,\mathbf{R}\text{e}~ s\} |\!|\!|\mathbf u|\!|\!|_1 \le |\!|\!|\mathbf u|\!|\!|_{|s|} \le \frac{|s|}{\min\{1,\mathbf{R}\text{e}~ s\}} |\!|\!|\mathbf u|\!|\!|_1 \qquad \forall \mathbf u \in \mathbf H^1(\Omega;\mathbb C) \quad s\in \mathbb C_+. \end{array}$$(18)

Finally, consider the bilinear form

$\begin{matrix} b (u, w; s) := a (u, w; s) + s^{2} (ρ u, w)_{Ω} \\ = (C (s) ε (u), ε (w))_{Ω} + s^{2} (ρ u, w)_{Ω} . \end{matrix}$ $$\begin{array}{} \displaystyle b(\mathbf u,\mathbf w;s) := a(\mathbf u,\mathbf w;s)+s^2 (\rho\mathbf u,\mathbf w)_\Omega \\ \displaystyle\qquad\qquad~ = (\mathbf C(s)\boldsymbol{\varepsilon} (\mathbf{u}),\boldsymbol{\varepsilon} (\mathbf {w}))_\Omega+s^2 (\rho\mathbf u,\mathbf w)_\Omega. \end{array}$$

Proposition 10

Ifψ, ϕandrare the functions and integer appearing in Hypotheses 2 and 3 and we define

$\begin{matrix} ψ_{⋆} (x) := min {x, ψ (x)}, ϕ_{⋆} (x) := max {x^{- r}, ϕ (x)}, \end{matrix}$ $$\begin{array}{} \displaystyle \psi_\star(x):=\min\{ x,\psi(x)\}, \qquad \phi_\star(x):=\max\{x^{-r},\phi(x)\}, \end{array}$$

then for alls ∈ ℂ₊

$\begin{matrix} R e b (u, \bar{s u}; s) \geq ψ_{⋆} (R e s) | | | u | | |_{| s |}^{2} \forall u \in H^{1} (Ω; C), \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf{R}\text{e}~ b(\mathbf u,\overline{s\mathbf u};s) \ge \psi_\star(\mathbf{R}\text{e}~ s) |\!|\!|\mathbf u|\!|\!|_{|s|}^2 \qquad \forall \mathbf u\in \mathbf H^1(\Omega;\mathbb C), \end{array}$$(19a)

$\begin{matrix} | b (u, w; s) | \leq | s |^{r} ϕ_{⋆} (R e s) | | | u | | |_{| s |} | | | w | | |_{| s |} \forall u, w \in H^{1} (Ω; C), \end{matrix}$ $$\begin{array}{} \displaystyle |b(\mathbf u,\mathbf w;s)| \le |s|^r \phi_\star(\mathbf{R}\text{e}~ s)|\!|\!|\mathbf u|\!|\!|_{|s|}|\!|\!|\mathbf w|\!|\!|_{|s|} \qquad \forall \mathbf u,\mathbf w\in \mathbf H^1(\Omega;\mathbb C), \end{array}$$(19b)

Ifℓ ≥ 0 andk ≥ 0 are the quantities in (3) and (4), then

$\begin{matrix} inf_{0 < x < 1} x^{- max {1, ℓ}} ψ_{⋆} (x) > 0, sup_{0 < x < 1} x^{max {r, k}} ϕ_{⋆} (x) < \infty . \end{matrix}$ $$\begin{array}{} \displaystyle \inf_{0 \lt x \lt 1} x^{-\max\{1,\ell\}}\psi_\star(x) \gt 0, \qquad \sup_{0 \lt x \lt 1} x^{\max\{r,k\}}\phi_\star(x) \lt \infty. \end{array}$$(20)

Proof

We have

$\begin{matrix} R e b (u, \bar{s u}; s) = (R e s) ∥ ρ^{1 / 2} s u ∥_{Ω}^{2} + R e a (u, \bar{s u}; s) \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf{R}\text{e}~ b(\mathbf u,\overline{s\mathbf u};s)= (\mathbf{R}\text{e}~ s) \|\rho^{1/2} s \mathbf u\|_\Omega^2 +\mathbf{R}\text{e} a(\mathbf u,\overline{s\mathbf u};s) \end{array}$$

and (19a) follows from (10). Similarly, by (8)

$\begin{matrix} | b (u, w; s) | \leq | s |^{r} ϕ (R e s) ∥ ε (u) ∥_{Ω} ∥ ε (w) ∥_{Ω} + ∥ ρ^{1 / 2} s u ∥_{Ω} ∥ ρ^{1 / 2} s w ∥_{Ω} \\ \leq max {| s |^{r} ϕ (R e s), 1} | | | u | | |_{| s |} | | | w | | |_{| s |} \end{matrix}$ $$\begin{array}{} \displaystyle |b(\mathbf u,\mathbf w;s)| \le |s|^r\phi(\mathbf{R}\text{e}~ s) \|\boldsymbol{\varepsilon}({\mathbf u})\|_\Omega\|\boldsymbol{\varepsilon}({\mathbf w})\|_\Omega +\|\rho^{1/2} s \mathbf u\|_\Omega\|\rho^{1/2} s \mathbf w\|_\Omega\\ \displaystyle\qquad\qquad~~\le \max\{ |s|^r \phi(\mathbf{R}\text{e}~ s),1\}|\!|\!|\mathbf u|\!|\!|_{|s|}|\!|\!|\mathbf w|\!|\!|_{|s|} \end{array}$$

and (19b) follows from the fact that 1 ≤ |s|/Res for all s ∈ ℂ₊. The asymptotic bounds (20) can be proved easily.

Lemma 11

There existsC_Ω,ρ > 0 such that for allc > 0 andα ∈ H^1/2(Γ_D), the solution of the variational problem

$\begin{matrix} \hat{u} \in H^{1} (Ω), γ_{D} \hat{u} = α, \end{matrix}$ $$\begin{array}{} \displaystyle \widehat{\mathbf u} \in \mathbf H^1(\Omega), \qquad \gamma_D\widehat{\mathbf u}=\boldsymbol\alpha, \end{array}$$(21a)

$\begin{matrix} (ε \hat{u}, ε (w))_{Ω} + c^{2} (ρ \hat{u}, w)_{Ω} = 0 \forall w \in H_{D}^{1} (Ω), \end{matrix}$ $$\begin{array}{} \displaystyle (\boldsymbol{\varepsilon}{\widehat{\mathbf u}},\boldsymbol{\varepsilon} (\mathbf {w}))_\Omega+ c^2 (\rho\widehat{\mathbf u},\mathbf w)_\Omega=0 \quad \forall \mathbf w\in \mathbf H^1_D(\Omega), \end{array}$$(21b)

satisfies

$\begin{matrix} | | | \hat{u} | | | c \leq C_{Ω, ρ} max {1, c}^{1 / 2} ∥ α ∥_{1 / 2, Γ_{D}} . \end{matrix}$ $$\begin{array}{} \displaystyle |\!|\!|{\widehat{\mathbf u}}|\!|\!|c\le C_{\Omega,\rho} \max\{1,c\}^{1/2} \| \boldsymbol\alpha\|_{1/2,\Gamma_D}. \end{array}$$

Proof

Let E : H^1/2(Γ_D) → H^1/2(Γ) be a bounded extension operator and consider the solution of the elliptic boundary value problem

$\begin{matrix} \tilde{u} \in H^{1} (Ω), γ \tilde{u} = E α, \end{matrix}$ $$\begin{array}{} \displaystyle \widetilde{\mathbf u} \in \mathbf H^1(\Omega),\qquad \gamma\widetilde{\mathbf u} =E\boldsymbol\alpha, \end{array}$$(22a)

$\begin{matrix} - Δ \tilde{u} + c^{2} \tilde{u} = 0 in Ω . \end{matrix}$ $$\begin{array}{} \displaystyle -\Delta \widetilde{\mathbf u}+c^2\widetilde{\mathbf u} =0 \quad\mbox{in}~~\Omega. \end{array}$$(22b)

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

$\begin{matrix} c^{2} ∥ \tilde{u} ∥_{Ω}^{2} + ∥ \nabla \tilde{u} ∥_{Ω}^{2} \leq C_{Ω} max {1, c} ∥ E α ∥_{1 / 2, Γ}^{2} \leq C_{Ω}^{'} max {1, c} ∥ α ∥_{1 / 2, Γ_{D}}^{2} . \end{matrix}$ $$\begin{array}{} \displaystyle c^2 \| \widetilde{\mathbf u}\|_\Omega^2 + \|\nabla\widetilde{\mathbf u}\|_\Omega^2 \le C_\Omega \max\{1,c\} \| E\boldsymbol\alpha\|_{1/2,\Gamma}^2 \le C_\Omega' \max\{1,c\} \| \boldsymbol\alpha\|_{1/2,\Gamma_D}^2. \end{array}$$

However, the solution of (21) minimizes ∥|û∥|_c among all û satisfying γ_Dû = α. Therefore

$\begin{matrix} | | | \hat{u} | | |_{c}^{2} \leq | | | \tilde{u} | | |_{c}^{2} \leq c^{2} ∥ ρ^{1 / 2} \tilde{u} ∥_{Ω}^{2} + ∥ \nabla \tilde{u} ∥_{Ω}^{2} \leq max {∥ ρ ∥_{L^{\infty} (Ω)}, 1} (c^{2} ∥ \tilde{u} ∥_{Ω}^{2} + ∥ \nabla \tilde{u} ∥_{Ω}^{2}) . \end{matrix}$ $$\begin{array}{} \displaystyle |\!|\!|{\widehat{\mathbf u}}|\!|\!|_c^2\le |\!|\!|{\widetilde{\mathbf u}}|\!|\!|_c^2 \le c^2 \|\rho^{1/2}\widetilde{\mathbf u}\|_\Omega^2 + \|\nabla\widetilde{\mathbf u}\|_\Omega^2 \le \max\{\|\rho\|_{L^\infty(\Omega)},1\} \Big(c^2\|\widetilde{\mathbf u}\|_\Omega^2 + \|\nabla\widetilde{\mathbf u}\|_\Omega^2\Big). \end{array}$$

This finishes the proof.

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 ∈ L²(Ω;ℂ), α ∈ H^1/2(Γ_D;ℂ), andβ ∈ H^−1/2(Γ_N;ℂ). For s ∈ ℂ₊, the solution of

$\begin{matrix} u \in H^{1} (Ω; C), γ_{D} u = α, \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u \in \mathbf H^1(\Omega;\mathbb C),\qquad \gamma_D \mathbf u =\boldsymbol\alpha, \end{array}$$(23a)

$\begin{matrix} b (u, w; s) = (f, w)_{Ω} + 〈 β, γ w 〉_{Γ_{N}} \forall w \in H_{D}^{1} (Ω; C), \end{matrix}$ $$\begin{array}{} \displaystyle b(\mathbf u,\mathbf w;s)=(\mathbf f,\mathbf w)_\Omega+\langle \boldsymbol\beta,\gamma\mathbf w\rangle_{\Gamma_N} \quad \forall \mathbf w\in \mathbf H^1_D(\Omega;\mathbb C), \end{array}$$(23b)

satisfies

$\begin{matrix} | | | u | | |_{| s |} \leq \frac{C}{ψ_{⋆} (R e s)} (∥ f ∥_{Ω} + \frac{| s |^{3 / 2 + r} ϕ_{⋆} (R e s)}{min {1, R e s}^{1 / 2}} ∥ α ∥_{1 / 2, Γ_{D}} + \frac{| s |}{min {1, R e s}} ∥ β ∥_{- 1 / 2, Γ_{N}}), \end{matrix}$ $$\begin{array}{} \displaystyle |\!|\!|\mathbf u|\!|\!|_{|s|} \le \frac{C}{\psi_\star(\mathbf{R}\text{e}~ s)}\Big( \|\mathbf f\|_\Omega + \frac{|s|^{3/2+r}\phi_\star(\mathbf{R}\text{e}~ s)}{\min\{1,\mathbf{R}\text{e}~ s\}^{1/2}} \| \boldsymbol\alpha\|_{1/2,\Gamma_D} + \frac{|s|}{\min\{1,\mathbf{R}\text{e}~ s\}} \|\boldsymbol\beta\|_{-1/2,\Gamma_N} \Big), \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{matrix} H_{D}^{1} \end{matrix}$ $\begin{array}{} \displaystyle \mathbf H^1_D \end{array}$ (Ω;ℂ) satisfies

$\begin{matrix} b (u, w; s) = (f, w)_{Ω} \forall w \in H_{D}^{1} (Ω; C), \end{matrix}$ $$\begin{array}{} \displaystyle b(\mathbf u,\mathbf w;s)=(\mathbf f,\mathbf w)_\Omega \quad\forall \mathbf w\in \mathbf H^1_D(\Omega;\mathbb C), \end{array}$$

and we can use w = su as test function. Applying (19a), it follows that

$\begin{matrix} ψ_{⋆} (R e s) | | | u | | |_{| s |}^{2} \leq R e b (u, \bar{s u}; s) = R e (f, \bar{s u})_{Ω} \\ \leq | (ρ^{- 1 / 2} f, ρ^{1 / 2} \bar{s u})_{Ω} | \leq ∥ ρ^{- 1 / 2} f ∥_{Ω} | | | u | | |_{| s |} . \end{matrix}$ $$\begin{array}{} \displaystyle \psi_\star(\mathbf{R}\text{e}~ s) |\!|\!|\mathbf u|\!|\!|_{|s|}^2 \le \mathbf{R}\text{e}~ b(\mathbf u,\overline{s\mathbf u};s)=\mathbf{R}\text{e} (\mathbf f,\overline{s\mathbf u})_\Omega \\ \displaystyle\qquad\qquad\qquad~ \le | (\rho^{-1/2}\mathbf f,\rho^{1/2}\overline{s\mathbf u})_\Omega| \le \|\rho^{-1/2}\mathbf f\|_\Omega |\!|\!|\mathbf u|\!|\!|_{|s|}. \end{array}$$

For the second one, we use the same argument to bound

$\begin{matrix} ψ_{⋆} (R e s) | | | u | | |_{| s |}^{2} \leq | 〈 β, \bar{s} γ \bar{u} 〉_{Γ_{N}} | \leq | s | ∥ β ∥_{- 1 / 2, Γ_{N}} ∥ γ u ∥_{1 / 2, Γ_{N}} \\ \leq C_{1} | s | ∥ β ∥_{- 1 / 2, Γ_{N}} ∥ u ∥_{1, Ω} \leq C_{2} | s | ∥ β ∥_{- 1 / 2, Γ_{N}} | | | u | | |_{1} \\ \leq C_{2} \frac{| s |}{min {1, R e s}} ∥ β ∥_{- 1 / 2, Γ_{N}} | | | u | | |_{| s |}, \end{matrix}$ $$\begin{array}{} \displaystyle \psi_\star(\mathbf{R}\text{e}~ s)|\!|\!|\mathbf u|\!|\!|_{|s|}^2 \le |\langle\boldsymbol\beta,\overline {s}\gamma\overline{\mathbf u}\rangle_{\Gamma_N}| \le |s| \|\boldsymbol\beta\|_{-1/2,\Gamma_N} \| \gamma\mathbf u\|_{1/2,\Gamma_N} \\\displaystyle\qquad\qquad\qquad~ \le C_1 |s| \|\boldsymbol\beta\|_{-1/2,\Gamma_N}\|\mathbf u\|_{1,\Omega} \le C_2 |s| \|\boldsymbol\beta\|_{-1/2,\Gamma_N}|\!|\!|\mathbf u|\!|\!|_1 \\\displaystyle\qquad\qquad\qquad~ \le C_2 \frac{|s|}{\min\{1,\mathbf{R}\text{e}~ s\}}\|\boldsymbol\beta\|_{-1/2,\Gamma_N} |\!|\!|\mathbf u|\!|\!|_{|s|}, \end{array}$$

where we have used Korn’s inequality and (1). For the final one, we write u = û + u₀, where û = û(α, |s|) is the solution of (21) with c = |s| and u₀ ∈ $\begin{matrix} H_{D}^{1} \end{matrix}$ $\begin{array}{} \displaystyle \mathbf H^1_D \end{array}$ (Ω;ℂ). Then

$\begin{matrix} b (u_{0}, w; s) = - b (\hat{u}, w; s) \forall w \in H_{D}^{1} (Ω; C) . \end{matrix}$ $$\begin{array}{} \displaystyle b(\mathbf u_0,\mathbf w;s)=-b(\widehat{\mathbf u},\mathbf w;s) \qquad\forall \mathbf w\in \mathbf H^1_D(\Omega;\mathbb C). \end{array}$$

Taking w = su₀ above, and using (19b), we can bound

$\begin{matrix} ψ_{⋆} (R e s) | | | u_{0} | | |_{| s |}^{2} \leq | b (\hat{u}, u_{0}; s) | \leq | s |^{r} ϕ_{⋆} (R e s) | | | \hat{u} | | |_{| s |} | s | | | | u_{0} | | |_{| s |} . \end{matrix}$ $$\begin{array}{} \displaystyle \psi_\star(\mathbf{R}\text{e}~ s) |\!|\!|{\mathbf u_0}|\!|\!|_{|s|}^2 \le | b(\widehat{\mathbf u},\mathbf u_0;s)| \le |s|^r \phi_\star(\mathbf{R}\text{e}~ s) |\!|\!|{\widehat{\mathbf u}}|\!|\!|_{|s|}\,|s|\, |\!|\!|{\mathbf u_0}|\!|\!|_{|s|}. \end{array}$$

Therefore,

$\begin{matrix} | | | u | | |_{| s |} \leq | | | \hat{u} | | |_{| s |} + | | | u_{0} | | |_{| s |} \leq (1 + \frac{| s |^{r + 1} ϕ_{⋆} (R e s)}{ψ_{⋆} (R e s)}) | | | \hat{u} | | |_{| s |} \\ \leq C_{Ω, ρ} (1 + \frac{| s |^{r + 1} ϕ_{⋆} (R e s)}{ψ_{⋆} (R e s)}) max {1, | s |}^{1 / 2} ∥ α ∥_{1 / 2, Γ_{D}} . \end{matrix}$ $$\begin{array}{} \displaystyle |\!|\!|\mathbf u|\!|\!|_{|s|}\le |\!|\!|{\widehat{\mathbf u}}|\!|\!|_{|s|}+|\!|\!|{\mathbf u_0}|\!|\!|_{|s|} \le \left(1+ \frac{|s|^{r+1}\phi_\star(\mathbf{R}\text{e}~ s)}{\psi_\star(\mathbf{R}\text{e}~ s)}\right) |\!|\!|{\widehat{\mathbf u}}|\!|\!|_{|s|} \\ \displaystyle\qquad\qquad\qquad\qquad\qquad\quad \le C_{\Omega,\rho} \left(1+ \frac{|s|^{r+1}\phi_\star(\mathbf{R}\text{e}~ s)}{\psi_\star(\mathbf{R}\text{e}~ s)}\right) \max\{1,|s|\}^{1/2} \|\boldsymbol\alpha\|_{1/2,\Gamma_D}. \end{array}$$

Using (13) and

$\begin{matrix} (1 + \frac{| s |^{r + 1} ϕ_{⋆} (R e s)}{ψ_{⋆} (R e s)}) \leq \frac{| s |}{ψ_{⋆} (R e s)} (1 + | s |^{r} ϕ_{⋆} (R e s)) \\ \leq \frac{2 | s |^{r + 1}}{ψ_{⋆} (R e s)} max {(R e s)^{- r}, ϕ_{⋆} (R e s)} = \frac{2 | s |^{r + 1}}{ψ_{⋆} (R e s)} ϕ_{⋆} (R e s), \end{matrix}$ $$\begin{array}{} \displaystyle \left(1+ \frac{|s|^{r+1}\phi_\star(\mathbf{R}\text{e}~ s)}{\psi_\star(\mathbf{R}\text{e}~ s)}\right) \le \frac{|s|}{\psi_\star(\mathbf{R}\text{e}~ s)} (1+|s|^r\phi_\star(\mathbf{R}\text{e}~ s)) \\\displaystyle \qquad\qquad\qquad\qquad\qquad\, \le \frac{2|s|^{r+1}}{\psi_\star(\mathbf{R}\text{e}~ s)}\max\{ (\mathbf{R}\text{e}~ s)^{-r},\phi_\star(\mathbf{R}\text{e}~ s)\} =\frac{2|s|^{r+1}}{\psi_\star(\mathbf{R}\text{e}~ s)}\phi_\star(\mathbf{R}\text{e}~ s), \end{array}$$

the result follows.

Distributional propagation of viscoelastic waves

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].

6.1

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

$\begin{matrix} ∥ x_{1} + ı x_{2} ∥_{X_{C}}^{2} := ∥ x_{1} ∥_{X}^{2} + ∥ x_{2} ∥_{X}^{2}, \forall x_{1}, x_{2} \in X, \end{matrix}$ $$\begin{array}{} \displaystyle \|x_1+\imath x_2\|_{X_\mathbb C}^2:=\|x_1\|^2_X+\|x_2\|^2_X, \qquad \forall x_1,x_2\in X, \end{array}$$

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 𝓑(X², Y²) formed by matrices of operators of the form

$\begin{matrix} (\begin{matrix} A & - B \\ B & A \end{matrix}) A, B \in B (X, Y) . \end{matrix}$ $$\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{matrix} \bar{A} x := \bar{A \bar{x}}, \end{matrix}$ $$\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

$\begin{matrix} \bar{H (s)} = H (\bar{s}) \forall s \in C_{+}, \end{matrix}$ $$\begin{array}{} \displaystyle \overline{H(s)}=H(\overline s) \qquad \forall s\in \mathbb C_+, \end{array}$$(24)

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

$\begin{matrix} ∥ H (s) ∥_{X_{C}} \leq | s |^{μ} ψ (R e s) \forall s \in C_{+}, \end{matrix}$ $$\begin{array}{} \displaystyle \| H(s)\|_{X_\mathbb C} \le |s|^\mu \psi(\mathbf{R}\text{e}~ s) \qquad \forall s\in \mathbb C_+, \end{array}$$(25)

where μ ∈ ℝ and ψ : (0, ∞) → (0, ∞) is non-increasing and at worst rational at the origin, i.e., sup_0<x<1x^kψ(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

$\begin{matrix} L {h * λ} (s) = L {h} (s) L {λ} (s) \forall s \in C_{+}, \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal L\{ h*\lambda\}(s)=\mathcal L\{h\}(s) \mathcal L\{\lambda\}(s) \qquad \forall s\in \mathbb C_+, \end{array}$$

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.

6.2

Viscoelastic material law and wave propagator

Theorem 13

If C : ℂ₊ → L^∞(Ω;ℂ^d×d×d×d) satisfies hypotheses (2), then there exists

$\begin{matrix} C \in T D (B (L^{2} (Ω))) \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal C\in \mathrm{TD}(\mathcal B(\mathbb L^2(\Omega)) ) \end{array}$$

such that 𝓛{𝓓M}(s) = C(s)M for all $\begin{matrix} M \in R_{sym}^{d \times d} . \end{matrix}$ $\begin{array}{} \displaystyle \mathrm M\in \mathbb R^{d\times d}_\text{sym}. \end{array}$For arbitraryu ∈ TD(H¹(Ω)), the convolution

$\begin{matrix} C * ε (u) \in T D (L^{2} (Ω)) \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal C* \boldsymbol{\varepsilon} (\mathbf{u}) \in \mathrm{TD}(\mathbb L^2(\Omega)) \end{array}$$

is well defined.

Proof

Let us first recall that for s ∈ ℂ₊, we have defined the operator

$\begin{matrix} L^{2} (Ω; C_{sym}^{d \times d}) ∋ U ⟼ C (s) U \in L^{2} (Ω; C_{sym}^{d \times d}) \end{matrix}$ $$\begin{array}{} \displaystyle L^2(\Omega;\mathbb C^{d\times d}_\text{sym})\ni \mathrm U \longmapsto \mathrm C(s)\mathrm U \in L^2(\Omega;\mathbb C^{d\times d}_\text{sym}) \end{array}$$

and that we have (see (7))

$\begin{matrix} ∥ C (s) U ∥_{Ω} \leq ∥ C (s) ∥_{max} ∥ U ∥_{Ω} \forall U \in L^{2} (Ω; C_{sym}^{d \times d}) . \end{matrix}$ $$\begin{array}{} \displaystyle \| \mathbf C(s)\mathrm U\|_\Omega \le \|\mathbf C(s)\|_{\max} \|\mathrm U\|_\Omega \qquad \forall \mathrm U\in L^2(\Omega;\mathbb C^{d\times d}_\text{sym}). \end{array}$$

Also (by (2a))

$\begin{matrix} \bar{C (s)} U = \bar{C (s) \bar{U}} = C (\bar{s}) U . \end{matrix}$ $$\begin{array}{} \displaystyle \overline{\mathbf C(s)}\mathrm U =\overline{\mathbf C(s)\overline{\mathrm U}} =\mathbf C(\overline s)\mathrm U. \end{array}$$

This means that

$\begin{matrix} C : C_{+} \to B (L^{2} (Ω; C_{sym}^{d \times d})) \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf C:\mathbb C_+\to \mathcal B(L^2(\Omega;\mathbb C^{d\times d}_\text{sym})) \end{array}$$

satisfies the conditions (24) and (25) and therefore there exists

$\begin{matrix} C \in T D (B (L^{2} (Ω))) \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal C\in\mathrm{TD}(\mathcal B(\mathbb L^2(\Omega))) \end{array}$$

such that 𝓛{𝓓} = C. If we fix $\begin{matrix} M \in R_{sym}^{d \times d} \end{matrix}$ $\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(𝕃²(Ω)) is C(s)M. Finally, if u ∈ TD(H¹(Ω)), then ε(u) ∈ TD(𝕃²(Ω)) and the convolution 𝓓 ∗ ε(u) is well defined.

The expression σ = 𝓓 ∗ ε(u) can be equivalently written σ = 𝓓 ∗ ε(u̇) where 𝓓 ∈ TD(𝓑(𝕃²(Ω))) is the distribution whose Laplace transform is s⁻¹ C(s). In the simplest example (the purely elastic case), C(s) = C₀, we can write 𝓓 = C₀ ⊗ δ₀ and 𝓓 = C₀ ⊗ H, where H is the Heaviside function. This yields the usual elastic law σ = C₀ε(u).

Proposition 14

Fors ∈ ℂ₊, let us consider the solution map

$\begin{matrix} S (s) : L^{2} (Ω; C) \times H^{1 / 2} (Γ_{D}; C) \times H^{- 1 / 2} (Γ_{N}; C) \to H^{1} (Ω; C) \end{matrix}$ $$\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

$\begin{matrix} S : C_{+} \to B (L^{2} (Ω; C) \times H^{1 / 2} (Γ_{D}; C) \times H^{- 1 / 2} (Γ_{N}; C); H^{1} (Ω; C)) \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm S:\mathbb C_+\to \mathcal B(\mathbf L^2(\Omega;\mathbb C)\times \mathbf H^{1/2}(\Gamma_D;\mathbb C)\times \mathbf H^{-1/2}(\Gamma_N;\mathbb C); \mathbf H^1(\Omega;\mathbb C)) \end{array}$$

is analytic and S(s) = S(s)for alls ∈ ℂ₊.

Proof

For s ∈ ℂ₊ and u ∈ H¹(Ω;ℂ), we define

$\begin{matrix} A (s) u := b (u, \cdot; s) = s^{2} (ρ u, \cdot)_{Ω} + (C (s) ε (u), ε \cdot)_{Ω} \in H_{D}^{1} (Ω; C)^{'} . \end{matrix}$ $$\begin{array}{} \displaystyle \mathbb A(s)\mathbf u:=b(\mathbf u,\,\cdot\,;s) =s^2(\rho\,\mathbf u,\,\cdot\,)_\Omega+(\mathbf C(s)\boldsymbol{\varepsilon} (\mathbf{u}),\boldsymbol{\varepsilon}{\,\cdot\,})_\Omega\in \mathbf H^1_D(\Omega;\mathbb C)'. \end{array}$$

Using the analyticity of C it is easy to see that

$\begin{matrix} A : C_{+} \to B (H^{1} (Ω; C), H_{D}^{1} (Ω; C)^{'}) \end{matrix}$ $$\begin{array}{} \displaystyle \mathbb A:\mathbb C_+\to \mathcal B(\mathbf H^1(\Omega;\mathbb C),\mathbf H^1_D(\Omega;\mathbb C)') \end{array}$$

is analytic. Moreover, the operator $\begin{matrix} A_{0} (s) := A (s) |_{H_{D}^{1} (Ω; C)} \end{matrix}$ $\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

$\begin{matrix} {T \in B (X, Y) : T is bijective} \to B (Y, X) \end{matrix}$ $$\begin{array}{} \displaystyle \{ T\in \mathcal B(X,Y)\,:\, T ~~\text{is bijective} \} \to \mathcal B(Y,X) \end{array}$$

is 𝓓^∞ and therefore the map s → 𝔸₀(s)⁻¹ is analytic from ℂ₊ to $\begin{matrix} B (H_{D}^{1} (Ω; C)^{'}, H_{D}^{1} (Ω; C)) \end{matrix}$ $\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 : H^1/2(Γ_D) → H¹(Ω) that is a right-inverse of γ_D and its natural extension to complex-valued functions, we can easily write

$\begin{matrix} S (s) (f, α, β) = L α + A_{0} (s)^{- 1} ((f, \cdot)_{Ω} + 〈 β, γ \cdot 〉_{Γ_{N}} - A (s) L α), \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm S(s)(\mathbf f,\boldsymbol\alpha,\boldsymbol\beta) =L\boldsymbol\alpha+ \mathbb A_0(s)^{-1}\left( (\mathbf f,\,\cdot\,)_\Omega+\langle\boldsymbol\beta,\gamma\,\cdot\,\rangle_{\Gamma_N} -\mathbb A(s)L\boldsymbol\alpha\right), \end{array}$$

which shows that S is analytic. Finally (2a) and a simple computation show that 𝔸(s) = 𝔸(s) and therefore S(s) = S(s).

Proposition 15

Fors ∈ ℂ₊, let us consider the solution map

$\begin{matrix} T (s) : L^{2} (Ω; C) \times H^{1 / 2} (Γ_{D}; C) \times H^{- 1 / 2} (Γ_{N}; C) \to H^{1} (Ω; C) \times H (d i v, Ω; C) \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm T(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)\times \mathbb H(\mathrm{div},\Omega;\mathbb C) \end{array}$$

defined by T(s)(f, α, β) := (u, C(s)ε(u)), whereu = S(s)(f, α, β). The function

$\begin{matrix} T : C_{+} \to B (L^{2} (Ω; C) \times H^{1 / 2} (Γ_{D}; C) \times H^{- 1 / 2} (Γ_{N}; C); H^{1} (Ω; C) \times H (d i v, Ω; C)) \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm T:\mathbb C_+\to \mathcal B(\mathbf L^2(\Omega;\mathbb C)\times \mathbf H^{1/2}(\Gamma_D;\mathbb C)\times \mathbf H^{-1/2}(\Gamma_N;\mathbb C); \, \mathbf H^1(\Omega;\mathbb C)\times \mathbb H(\mathrm{div},\Omega;\mathbb C)) \end{array}$$

is analytic and T(s) = T(s)for all s ∈ ℂ₊. Finally

$\begin{matrix} T (s)^{- 1} (u, σ) = (ρ s^{2} u - d i v σ, γ_{D} u, γ_{N} σ) \forall s \in C_{+} . \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm T(s)^{-1}(\mathbf u,\boldsymbol\sigma)=(\rho\,s^2\,\mathbf u-\mathrm{div}\boldsymbol\sigma,\gamma_D\mathbf u,\gamma_N\boldsymbol\sigma) \qquad \forall s\in \mathbb C_+. \end{array}$$

Proof

If u = S(s)(f, α, β) and σ = C(s)ε(u), then

$\begin{matrix} u \in H^{1} (Ω; C), σ \in H (d i v, Ω; C), \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u\in \mathbf H^1(\Omega;\mathbb C), \quad \boldsymbol\sigma\in \mathbb H(\mathrm{div},\Omega;\mathbb C), \end{array}$$(26a)

$\begin{matrix} ρ s^{2} u = d i v σ + f, \end{matrix}$ $$\begin{array}{} \displaystyle \rho\,s^2\mathbf u=\mathrm{div}\,\boldsymbol\sigma +\mathbf f, \end{array}$$(26b)

$\begin{matrix} σ = C (s) ε (u), \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol \sigma =\mathrm C(s)\boldsymbol{\varepsilon}({\mathbf u}), \end{array}$$(26c)

$\begin{matrix} γ_{D} u = α, γ_{N} σ = β . \end{matrix}$ $$\begin{array}{} \displaystyle \gamma_D \mathbf u=\boldsymbol\alpha, \quad \gamma_N\boldsymbol\sigma=\boldsymbol\beta. \end{array}$$(26d)

Since div σ = ρs² 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.

Proposition 16

There exists a distribution

$\begin{matrix} T \in T D (B (L^{2} (Ω) \times H^{1 / 2} (Γ_{D}) \times H^{- 1 / 2} (Γ_{N}), H^{1} (Ω) \times H (d i v, Ω))) \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal T\in \mathrm{TD}(\mathcal B(\mathbf L^2(\Omega)\times\mathbf H^{1/2}(\Gamma_D)\times \mathbf H^{-1/2}(\Gamma_N), \mathbf H^1(\Omega)\times \mathbb H(\mathrm{div},\Omega))) \end{array}$$

such that for all

$\begin{matrix} f \in T D (L^{2} (Ω)), α \in T D (H^{1 / 2} (Γ_{D})), β \in T D (H^{- 1 / 2} (Γ_{N})) \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf f\in \mathrm{TD}(\mathbf L^2(\Omega)), \qquad \boldsymbol\alpha\in \mathrm{TD}(\mathbf H^{1/2}(\Gamma_D)), \qquad \boldsymbol\beta\in \mathrm{TD}(\mathbf H^{-1/2}(\Gamma_N)) \end{array}$$(27)

the pair (u, σ) = 𝓣 ∗ (f, α, β) is the unique solution to

$\begin{matrix} u \in T D (H^{1} (Ω)), σ \in T D (H (d i v, Ω)), \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u \in \mathrm{TD}(\mathbf H^1(\Omega)), \quad \boldsymbol\sigma\in \mathrm{TD}(\mathbb H(\mathrm{div},\Omega)), \end{array}$$(28a)

$\begin{matrix} ρ \ddot{u} = d i v σ + f, \end{matrix}$ $$\begin{array}{} \displaystyle \rho\,\ddot{\mathbf u} =\mathrm{div}\,\boldsymbol\sigma+\mathbf f, \end{array}$$(28b)

$\begin{matrix} σ = C * ε (u), \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol\sigma=\mathcal C*\boldsymbol{\varepsilon} (\mathbf{u}), \end{array}$$(28c)

$\begin{matrix} γ_{D} u = α, γ_{N} σ = β . \end{matrix}$ $$\begin{array}{} \displaystyle \gamma_D\mathbf u =\boldsymbol\alpha, \quad \gamma_N\boldsymbol\sigma=\boldsymbol\beta. \end{array}$$(28d)

Proof

Let (U, Σ) = T(s)(F, A, B). We first estimate

$\begin{matrix} ∥ Σ ∥_{Ω} + ∥ d i v Σ ∥_{Ω} \leq | s |^{r} ϕ (R e s) ∥ ε U ∥_{Ω} + | s |^{2} ∥ ρ U ∥_{Ω} + ∥ F ∥_{Ω} . \end{matrix}$ $$\begin{array}{} \displaystyle \|\boldsymbol\Sigma \|_\Omega+\|\mathrm{div}\,\boldsymbol\Sigma \|_\Omega \le |s|^r\phi(\mathbf{R}\text{e}~ s) \|\boldsymbol{\varepsilon}{\boldsymbol U }\|_\Omega+|s|^2 \|\rho\,\boldsymbol U \|_\Omega+\|\boldsymbol F\|_\Omega. \end{array}$$

By Korn’s inequality, (18) and Theorem 12, we can also bound

$\begin{matrix} ∥ U ∥_{1, Ω} \leq C | | | U | | |_{1} \leq \frac{C}{min {1, R e s}} | | | U | | |_{| s |} \\ \leq ϕ_{1} (R e s) ∥ F ∥_{Ω} + | s |^{3 / 2 + r} ϕ_{2} (R e s) ∥ A ∥_{1 / 2, Γ_{D}} + | s | ϕ_{3} (R e s) ∥ B ∥_{- 1 / 2, Γ_{N}}, \end{matrix}$ $$\begin{array}{} \displaystyle \| \boldsymbol U \|_{1,\Omega} \le C |\!|\!|{\boldsymbol U }|\!|\!|_1 \le \frac{C}{\min\{1,\mathbf{R}\text{e}~ s\}}|\!|\!|{\boldsymbol U }|\!|\!|_{|s|} \\\displaystyle\qquad\qquad\qquad\quad~~ \le \phi_1(\mathbf{R}\text{e}~ s) \| \boldsymbol F \|_\Omega + |s|^{3/2+r}\phi_2(\mathbf{R}\text{e}~ s) \|\boldsymbol A \|_{1/2,\Gamma_D} + |s|\phi_3(\mathbf{R}\text{e}~ s) \|\boldsymbol B \|_{-1/2,\Gamma_N}, \end{array}$$

where

$\begin{matrix} ϕ_{1} (x) := \frac{C_{1}}{ψ_{⋆} (x) min {1, x}}, ϕ_{2} (x) := \frac{C_{2} ϕ_{⋆} (x)}{ψ_{⋆} (x) min {1, x^{3 / 2}}}, ϕ_{3} (x) := \frac{C_{3}}{ψ_{⋆} (x) min {1, x^{2}}} . \end{matrix}$ $$\begin{array}{} \displaystyle \phi_1(x):= \frac{C_1}{\psi_\star(x)\min\{1,x\}}, \quad \phi_2(x):=\frac{C_2\phi_\star(x)}{\psi_\star(x)\min\{1,x^{3/2}\}}, \quad \phi_3(x):=\frac{C_3}{\psi_\star(x) \min\{1,x^2\}}. \end{array}$$

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

$\begin{matrix} E \in T D (B (H^{1} (Ω) \times H (d i v, Ω), L^{2} (Ω) \times H^{1 / 2} (Γ_{D}) \times H^{- 1 / 2} (Γ_{N})) \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal E\in \mathrm{TD}( \mathcal B(\mathbf H^1(\Omega)\times \mathbb H(\mathrm{div},\Omega), \mathbf L^2(\Omega)\times\mathbf H^{1/2}(\Gamma_D)\times \mathbf H^{-1/2}(\Gamma_N)) \end{array}$$

such that 𝓛{𝓔}(s) = T(s)⁻¹ for all s ∈ ℂ₊. In fact,

$\begin{matrix} E * (u, σ) = (ρ \ddot{u} - d i v σ, γ_{D} u, γ_{N} σ) \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal E*(\mathbf u,\boldsymbol\sigma)=(\rho\ddot{\mathbf u}-\mathrm{div}\,\boldsymbol\sigma,\gamma_D\mathbf u,\gamma_N\boldsymbol\sigma) \end{array}$$

for all (u, σ) ∈ TD(H¹(Ω)× ℍ(div, Ω)). Since T(s)⁻¹ T(s) is the identity operator for all s ∈ ℂ₊, it follows that

$\begin{matrix} E * T * (f, α, β) = (f, α, β) \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal E*\mathcal T*(\mathbf f,\boldsymbol\alpha,\boldsymbol\beta)=(\mathbf f,\boldsymbol\alpha,\boldsymbol\beta) \end{array}$$

and therefore (u, σ) solves (28). Finally, if (u, σ) solves (28), then U = 𝓛{u} and Σ = 𝓛{ σ} satisfy

$\begin{matrix} U (s) \in H^{1} (Ω; C), Σ (s) \in H (d i v, Ω; C), \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol U(s)\in \mathbf H^1(\Omega;\mathbb C), \quad \boldsymbol\Sigma(s)\in \mathbb H(\mathrm{div},\Omega;\mathbb C), \end{array}$$(29a)

$\begin{matrix} ρ s^{2} U (s) = d i v Σ (s) + F (s), \end{matrix}$ $$\begin{array}{} \displaystyle \rho\,s^2\boldsymbol U(s)=\mathrm{div}\,\boldsymbol\Sigma(s)+\boldsymbol F(s), \end{array}$$(29b)

$\begin{matrix} Σ (s) = C (s) ε U (s), \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol \Sigma(s)=\mathrm C(s)\boldsymbol{\varepsilon}{\boldsymbol U(s)}, \end{array}$$(29c)

$\begin{matrix} γ_{D} U (s) = A (s), γ_{N} Σ (s) = B (s), \end{matrix}$ $$\begin{array}{} \displaystyle \gamma_D \boldsymbol U(s)=\boldsymbol A(s), \quad \gamma_N\boldsymbol\Sigma(s)=\boldsymbol B(s), \end{array}$$(29d)

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

$\begin{matrix} ∥ F (s) ∥_{X \to Y} \leq | s |^{m + μ} φ (R e s) \forall s \in C_{+}, \end{matrix}$ $$\begin{array}{} \displaystyle \| \mathrm F(s)\|_{X\to Y}\le |s|^{m+\mu}\varphi(\mathbf{R}\text{e}~ s) \qquad \forall s\in \mathbb C_+, \end{array}$$(30)

where:

m ≥ 0 is an integer, andμ ∈ [0, 1),

φ:(0, ∞) → (0, ∞) is non-increasing and sup_0<x<1x^kφ(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

$\begin{matrix} ∥ u (t) ∥_{Y} \leq C_{μ} {(\frac{t}{1 + t})}^{μ} φ (t^{- 1}) \sum_{ℓ = m}^{m + 2} \int_{0}^{t} ∥ λ^{(ℓ)} (τ) ∥_{X} d τ \forall t \geq 0. \end{matrix}$ $$\begin{array}{} \displaystyle \| u(t)\|_Y \le C_\mu \left(\frac{t}{1+t}\right)^\mu \varphi(t^{-1}) \sum\limits_{\ell=m}^{m+2}\int_0^t \|\lambda^{(\ell)}(\tau)\|_X \mathrm d\tau \quad\forall t\ge 0. \end{array}$$(31)

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.

We now consider Sobolev spaces

$\begin{matrix} W_{+}^{m, 1} (0, \infty; X) := {f \in C^{m - 1} ([0, \infty); X) : f^{(m)} \in L^{1} (0, \infty; X), f^{(ℓ)} (0) = 0 ℓ \leq m - 1}, \end{matrix}$ $$\begin{array}{} \displaystyle W^{m,1}_+(0,\infty;X) :=\{ f\in \mathcal C^{m-1}([0,\infty);X)\,:\, f^{(m)}\in L^1(0,\infty;X), \, f^{(\ell)}(0)=0 \quad \ell\le m-1\}, \end{array}$$

where X is any Hilbert space.

Corollary 18

Let F be as in Theorem 17 and $\begin{matrix} λ \in W_{+}^{m + 2, 1} (0, \infty; X) . \end{matrix}$ $\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

$\begin{matrix} g^{(- 1)} (t) := \int_{0}^{t} g (τ) d τ . \end{matrix}$ $$\begin{array}{} \displaystyle g^{(-1)}(t):=\int_0^t g(\tau)\mathrm d\tau. \end{array}$$

Consider now data satisfying

$\begin{matrix} f \in W_{+}^{m (f), 1} (0, \infty; L^{2} (Ω)), \\ α \in W_{+}^{m (α), 1} (0, \infty; H^{1 / 2} (Γ_{D})), \\ β \in W_{+}^{m (β), 1} (0, \infty; H^{- 1 / 2} (Γ_{N})), \end{matrix}$ $$\begin{array}{} \displaystyle ~\mathbf f \in W^{m(f),1}_+(0,\infty;\mathbf L^2(\Omega)),\\\displaystyle \boldsymbol\alpha \in W^{m(\alpha),1}_+(0,\infty;\mathbf H^{1/2}(\Gamma_D)),\\\displaystyle \boldsymbol\beta \in W^{m(\beta),1}_+(0,\infty;\mathbf H^{-1/2}(\Gamma_N)), \end{array}$$

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)∥.

Theorem 19

$\begin{matrix} m (f) = 2, m (α) = 3 + r, m (β) = 3, \end{matrix}$ $$\begin{array}{} \displaystyle m(f)=2,\quad m(\alpha)=3+r, \quad m(\beta)=3, \end{array}$$

thenu ∈ 𝓓¹([0, ∞);L²(Ω)) ∩ 𝓓([0, ∞);H¹(Ω)), and we have the estimates

$\begin{matrix} ∥ \dot{u} (t) ∥_{Ω} + ∥ ε (u (t)) ∥_{Ω} \leq C_{f} (t) \sum_{k = 0}^{2} \int_{0}^{t} ∥ f^{(k)} (τ) ∥_{Ω} d τ + C_{α} (t) \sum_{k = 1 + r}^{3 + r} \int_{0}^{t} ∥ α^{(k)} (τ) ∥_{1 / 2, Γ_{D}} d τ \\ + C_{β} (t) \sum_{k = 1}^{3} \int_{0}^{t} ∥ β^{(k)} (τ) ∥_{- 1 / 2, Γ_{N}} d τ, \end{matrix}$ $$\begin{array}{} \displaystyle \| \dot{\mathbf u}(t)\|_\Omega +\|\boldsymbol{\varepsilon}({\mathbf u(t))}\|_\Omega \le C_f(t) \sum\limits_{k=0}^2 \int_0^t \|\mathbf f^{(k)}(\tau)\|_\Omega\mathrm d\tau + C_\alpha(t) \sum\limits_{k=1+r}^{3+r} \int_0^t \|\boldsymbol\alpha^{(k)}(\tau)\|_{1/2,\Gamma_D}\mathrm d\tau \\ \displaystyle\qquad\qquad\qquad\qquad\qquad~ +C_\beta(t) \sum\limits_{k=1}^3 \int_0^t \|\boldsymbol\beta^{(k)}(\tau)\|_{-1/2,\Gamma_N}\mathrm d\tau, \end{array}$$(32a)

$\begin{matrix} ∥ u (t) ∥_{Ω} \leq C_{f} (t) \sum_{k = - 1}^{1} \int_{0}^{t} ∥ f^{(k)} (τ) ∥_{Ω} d τ + C_{α} (t) \sum_{k = r}^{2 + r} \int_{0}^{t} ∥ α^{(k)} (τ) ∥_{1 / 2, Γ_{D}} d τ \\ + C_{β} (t) \sum_{k = 0}^{2} \int_{0}^{t} ∥ β^{(k)} (τ) ∥_{- 1 / 2, Γ_{N}} d τ . \end{matrix}$ $$\begin{array}{} \displaystyle \|\mathbf u(t)\|_\Omega \le C_f(t) \sum\limits_{k=-1}^1 \int_0^t \|\mathbf f^{(k)}(\tau)\|_\Omega\mathrm d\tau + C_\alpha(t) \sum\limits_{k=r}^{2+r} \int_0^t \|\boldsymbol\alpha^{(k)}(\tau)\|_{1/2,\Gamma_D}\mathrm d\tau \\ \displaystyle\qquad\qquad~~ +C_\beta(t) \sum\limits_{k=0}^2 \int_0^t \|\boldsymbol\beta^{(k)}(\tau)\|_{-1/2,\Gamma_N}\mathrm d\tau. \end{array}$$(32b)

with

$\begin{matrix} C_{f} (t) := \frac{C_{1}}{ψ_{⋆} (t^{- 1})}, C_{α} (t) := \frac{C_{3} max {1, t^{1 / 2}} ϕ_{⋆} (t^{- 1})}{ψ_{⋆} (t^{- 1})} {(\frac{t}{1 + t})}^{1 / 2}, \\ C_{β} (t) := \frac{C_{2} max {1, t}}{ψ_{⋆} (t^{- 1})}, \end{matrix}$ $$\begin{array}{} \displaystyle C_f(t) :=\frac{C_1}{\psi_\star(t^{-1})}, \qquad C_\alpha(t) := \frac{C_3 \max\{1,t^{1/2}\} \phi_\star(t^{-1})}{\psi_\star(t^{-1})} \left(\frac{t}{1+t}\right)^{1/2},\\\displaystyle C_\beta(t) := \frac{C_2\max\{1,t\}}{\psi_\star(t^{-1})}, \end{array}$$

for some constantsC₁, C₂, C₃. If additionally

$\begin{matrix} m (f) = 2 + r, m (α) = 3 + 2 r, m (β) = 3 + r, \end{matrix}$ $$\begin{array}{} \displaystyle m(f)=2+r, \qquad m(\alpha)=3+2r, \qquad m(\beta)=3+r, \end{array}$$

thenσ ∈ 𝓓([0, ∞); 𝕃²(Ω)) and we have the bound

$\begin{matrix} ∥ σ (t) ∥_{Ω} \leq C_{f} (t) \sum_{k = r}^{2 + r} \int_{0}^{t} ∥ f^{(k)} (τ) ∥_{Ω} d τ + C_{α} (t) \sum_{k = 1 + 2 r}^{3 + 2 r} \int_{0}^{t} ∥ α^{(k)} (τ) ∥_{1 / 2, Γ_{D}} d τ \\ + C_{β} (t) \sum_{k = 1 + r}^{3 + r} \int_{0}^{t} ∥ β^{(k)} (τ) ∥_{- 1 / 2, Γ_{N}} d τ . \end{matrix}$ $$\begin{array}{} \displaystyle \|\boldsymbol\sigma(t)\|_\Omega \le C_f(t) \sum\limits_{k=r}^{2+r} \int_0^t \|\mathbf f^{(k)}(\tau)\|_\Omega\mathrm d\tau + C_\alpha(t) \sum\limits_{k=1+2r}^{3+2r} \int_0^t \|\boldsymbol\alpha^{(k)}(\tau)\|_{1/2,\Gamma_D}\mathrm d\tau \\ \displaystyle\qquad\qquad~~ +C_\beta(t) \sum\limits_{k=1+r}^{3+r} \int_0^t \|\boldsymbol\beta^{(k)}(\tau)\|_{-1/2,\Gamma_N}\mathrm d\tau. \end{array}$$(32c)

Proof

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

$\begin{matrix} S (s) = S_{f} (s) + S_{α} (s) + S_{β} (s) . \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm S(s)=\mathrm S_f(s)+\mathrm S_\alpha(s)+\mathrm S_\beta(s). \end{array}$$

We will also use the embedding of H¹(Ω) into L²(Ω). The following table lists all nine operators:

F(s)	X	Y	φ(x)	m	μ	λ	u
sI_H¹→L² S_f(s)	L²(Ω)	L²(Ω)	1/ψ_*(x)	0	0	f	u̇
ε ∘ S_f(s)	L²(Ω)	𝕃²(Ω)	1/ψ_*(x)	0	0	f	ε(u)
C(s)(ε ∘ S_f(s))	L²(Ω)	𝕃²(Ω)	ϕ(x)/ψ_*(x)	r	0	f	σ
s I_H¹→L² S_α(s)	H^1/2(Γ_D)	L²(Ω)	$\begin{matrix} \frac{ϕ_{⋆} (x)}{ψ_{⋆} (x) min {1, \sqrt{x}}} \end{matrix}$ $\begin{array}{} \displaystyle \frac{\phi_\star(x)}{\psi_\star(x)\min\{1,\sqrt x\}} \end{array}$	1 + r	$\begin{matrix} \frac{1}{2} \end{matrix}$ $\begin{array}{} \displaystyle \frac{1}{2} \end{array}$	α	u̇
ε ∘ S_α(s)	H^1/2(Γ_D)	𝕃²(Ω)	$\begin{matrix} \frac{ϕ_{⋆} (x)}{ψ_{⋆} (x) min {1, \sqrt{x}}} \end{matrix}$ $\begin{array}{} \displaystyle \frac{\phi_\star(x)}{\psi_\star(x)\min\{1,\sqrt x\}} \end{array}$	1 + r	$\begin{matrix} \frac{1}{2} \end{matrix}$ $\begin{array}{} \displaystyle \frac{1}{2} \end{array}$	α	ε(u)
C(s)(ε ∘ S_α(s))	H^1/2(Γ_D)	𝕃²(Ω)	$\begin{matrix} \frac{ϕ (x) ϕ_{⋆} (x)}{ψ_{⋆} (x) min {1, \sqrt{x}}} \end{matrix}$ $\begin{array}{} \displaystyle \frac{\phi(x)\phi_\star(x)}{\psi_\star(x)\min\{1,\sqrt x\}} \end{array}$	1 + 2r	$\begin{matrix} \frac{1}{2} \end{matrix}$ $\begin{array}{} \displaystyle \frac{1}{2} \end{array}$	α	σ
s I_H¹→L² S_β(s)	H^−1/2(Γ_N)	L²(Ω)	$\begin{matrix} \frac{1}{min {1, x} ψ_{⋆} (x)} \end{matrix}$ $\begin{array}{} \displaystyle \frac1{\min\{1,x\}\psi_\star(x)} \end{array}$	1	0	β	u̇
ε ∘ S_β(s)	H^−1/2(Γ_N)	𝕃²(Ω)	$\begin{matrix} \frac{1}{min {1, x} ψ_{⋆} (x)} \end{matrix}$ $\begin{array}{} \displaystyle \frac1{\min\{1,x\}\psi_\star(x)} \end{array}$	1	0	β	ε(u)
C(s)(ε ∘ S_β(s))	H^−1/2(Γ_N)	𝕃²(Ω)	$\begin{matrix} \frac{ϕ (x)}{min {1, x} ψ_{⋆} (x)} \end{matrix}$ $\begin{array}{} \displaystyle \frac{\phi(x)}{\min\{1,x\}\psi_\star(x)} \end{array}$	1 + r	0	β	σ

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 L²(Ω)), then (f⁽⁻¹⁾, α⁽⁻¹⁾, β⁽⁻¹⁾) ↦ u.

In general (see Proposition 10)

$\begin{matrix} \frac{1}{ψ_{⋆} (t^{- 1})} ≲ t^{max {1, ℓ}}, ϕ (t^{- 1}) ≲ t^{k}, ϕ_{⋆} (t^{- 1}) ≲ t^{max {r, k}} \end{matrix}$ $$\begin{array}{} \displaystyle \frac1{\psi_\star(t^{-1})}\lesssim t^{\max\{1,\ell\}}, \qquad \phi(t^{-1}) \lesssim t^k, \qquad \phi_\star(t^{-1})\lesssim t^{\max\{r,k\}} \end{array}$$

so all bounds in Theorem 19 polynomial. In Zener’s model and its fractional version (see Proposition 2 and 3) we have r = 0 and

$\begin{matrix} ψ (x) = c_{0} x, & ψ_{⋆} (x) = min {1, c_{0}} x, \\ ϕ (x) = C max {1, x^{- 2}}, & ϕ_{⋆} (x) \leq max {1, C} max {1, x^{- 2}} . \end{matrix}$ $$\begin{array}{} \displaystyle \psi(x) =c_0\,x, \qquad& \psi_\star(x)=\min\{1,c_0\} x,\\ \phi(x)=C\max\{1,x^{-2}\},\qquad & \phi_\star(x) \le \max\{1,C\} \max\{1,x^{-2}\}. \end{array}$$

Technical work towards a time domain analysis

In Section 9, we give a different analysis, based on the theory of C₀-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 C_diff ∈ 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

$\begin{matrix} C_{diff} : L^{2} (Ω) ⟶ L^{2} (Ω), \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm C_\text{diff}:\mathbb L^2(\Omega)\longrightarrow \mathbb L^2(\Omega), \end{array}$$

which is bounded, selfadjoint, and positive semidefinite. We also consider a strictly positive function a ∈ L^∞(Ω), and the associated multiplication operator

$\begin{matrix} L^{2} (Ω) ∋ E ⟼ m_{a} E := a E \in L^{2} (Ω) . \end{matrix}$ $$\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 m_aC_diff = C_diffm_a. We next consider the following seminorm in 𝕃²(Ω)

$\begin{matrix} | S |_{c}^{2} := (a C_{d i f f} S, S)_{Ω} . \end{matrix}$ $$\begin{array}{} \displaystyle |\mathrm S|_c^2:=(a\,\mathrm {C_{diff}} \mathrm S,\mathrm S)_\Omega. \end{array}$$

Proposition 20

The following properties hold:

|S|_c ≲ ∥ S∥_Ωfor all S ∈ 𝕃²(Ω).

∥C_diffS∥_Ω ≲ |S|_cfor all S ∈ 𝕃²(Ω).

|S|_c = 0 if and only if S ∈ kerC_diff.

Ifb ∈ L^∞(Ω), then |bS|_c ≲ |S|_cfor all S ∈ 𝕃²(Ω).

Proof

Since m_a and C_diff are bounded linear operators in 𝕃²(Ω), the proof of (a) is straightforward

$\begin{matrix} | S |_{c}^{2} = (a C_{d i f f} S, S)_{Ω} \leq ∥ m_{a} C_{d i f f} S ∥_{Ω} ∥ S ∥_{Ω} ≲ ∥ S ∥_{Ω}^{2} . \end{matrix}$ $$\begin{array}{} \displaystyle |\mathrm S|_c^2 =(a\mathrm {C_{diff}}\mathrm S,\mathrm S)_\Omega \le \| m_a\mathrm {C_{diff}}\mathrm S\|_\Omega \|\mathrm S\|_\Omega \lesssim \|\mathrm S\|_\Omega^2. \end{array}$$

To prove (b) note first that since m_1/a is bounded, then

$\begin{matrix} ∥ C_{d i f f} S ∥_{Ω} ≲ ∥ a C_{d i f f} S ∥_{Ω} \end{matrix}$ $$\begin{array}{} \displaystyle \|\mathrm {C_{diff}}\mathrm S\|_\Omega \lesssim \| a\mathrm {C_{diff}}\mathrm S\|_\Omega \end{array}$$(33)

and therefore, using (a)

$\begin{matrix} ∥ a C_{d i f f} S ∥_{Ω}^{2} = (a C_{d i f f} S, a C_{d i f f} S)_{Ω} = (S, a C_{d i f f} S)_{c} \\ \leq | S |_{c} | a C_{d i f f} S |_{c} ≲ | S |_{c} ∥ a C_{d i f f} S ∥_{Ω} . \end{matrix}$ $$\begin{array}{} \displaystyle \,\| a\mathrm {C_{diff}}\mathrm S\|_\Omega^2 =(a\mathrm {C_{diff}}\mathrm S,a\mathrm {C_{diff}}\mathrm S)_\Omega =(\mathrm S,a\mathrm {C_{diff}}\mathrm S)_c \\ \displaystyle\qquad\qquad\quad\le |\mathrm S|_c |a\mathrm {C_{diff}}\mathrm S|_c \lesssim |\mathrm S|_c \|a\mathrm {C_{diff}}\mathrm S\|_\Omega. \end{array}$$

The result follows then by (33). The property (c) is a direct consequence of (b). Finally, to prove (d) we write

$\begin{matrix} | b S |_{c}^{2} = (a C_{d i f f} (b S), b S)_{Ω} = (b^{2} a C_{d i f f} S, S)_{Ω} \\ = \int_{Ω} b^{2} a \underset{\geq 0}{\underset{⏟}{(C_{d i f f} S) : S}} ≲ \int_{Ω} a (C_{d i f f} S) : S = | S |_{c}^{2} . \end{matrix}$ $$\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 𝕃²(Ω) 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

$\begin{matrix} M & := {S \in L^{2} (Ω) : (S, T)_{Ω} = 0 \forall T, | T |_{c} = 0} \\ = {T \in L^{2} (Ω) : | T |_{c} = 0}^{⊥} = (\ker C_{d i f f})^{⊥}, \end{matrix}$ $$\begin{array}{} \displaystyle M\!\!\!\! &:=\{\mathrm S\in \mathbb{L}^2(\Omega)\,:\, (\mathrm S,\mathrm T)_\Omega=0 \qquad \forall \mathrm T\,,\, |\mathrm T|_c=0\} \\ &=\{\mathrm T\in \mathbb{L}^2(\Omega)\,:\, |\mathrm T|_c=0\}^\perp =(\ker\mathrm {C_{diff}})^\perp, \end{array}$$

which is clearly closed in 𝕃²(Ω). If T ∈ M and |T|_c = 0, then T ∈ (kerC_diff)^⊥ ∩ kerC_diff = {0} (Proposition 20(c)) and therefore |⋅|_c is a norm in M. Note also that if S ∈ 𝕃²(Ω), we can decompose

$\begin{matrix} S = S_{M} + S_{C}, S_{M} \in M, S_{C} \in \ker C_{d i f f}, \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm S=\mathrm S_M+\mathrm S_C, \qquad \mathrm S_M\in M, \qquad \mathrm S_C\in \ker\mathrm {C_{diff}}, \end{array}$$

and C_diffS = C_diffS_M. 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 Π : 𝕃²(Ω) → M that performs the 𝕃²(Ω) 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{matrix} R := I Π : L^{2} (Ω) \to \hat{M}, \end{matrix}$ $$\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 C_diff|_M : M → 𝕃²(Ω) is bounded by Proposition 20(b) and, therefore, there exists a unique bounded extension

$\begin{matrix} \hat{C_{d i f f}} : \hat{M} \to L^{2} (Ω) . \end{matrix}$ $$\begin{array}{} \displaystyle \widehat{\mathrm {C_{diff}}}:\widehat{M} \to \mathbb{L}^2(\Omega). \end{array}$$

An easy argument about extensions and the fact that M = (kerC_diff)^⊥ prove that

$\begin{matrix} \hat{C_{d i f f}} R = \hat{C_{d i f f}} I Π = C_{d i f f} |_{M} Π = C_{d i f f} . \end{matrix}$ $$\begin{array}{} \displaystyle \widehat{\mathrm {C_{diff}}} \mathrm R =\widehat{\mathrm {C_{diff}}}\,\mathrm I\, \Pi =\mathrm {C_{diff}}|_M \,\Pi =\mathrm {C_{diff}}. \end{array}$$(34)

We now consider b ∈ L^∞(Ω) and the associated multiplication operatorm_b : 𝕃²(Ω) → 𝕃²(Ω) that we recall was given by m_b E := b E. Since m_b C_diff = C_diffm_b, it follows that m_b maps kerC_diff into kerC_diff. Therefore, if S ∈ M, then

$\begin{matrix} (b S, T)_{Ω} = (S, b T)_{Ω} = 0 \forall T \in \ker C_{d i f f}, \end{matrix}$ $$\begin{array}{} \displaystyle (b\mathrm S,\mathrm T)_\Omega =(\mathrm S,b\mathrm T)_\Omega=0 \qquad \forall \mathrm T\in \ker\mathrm {C_{diff}}, \end{array}$$

which shows that bS ∈ M. The bounded linear map m_b|_M : M → M (see Proposition 20(d) can then be extended to a bounded linear selfadjoint map

$\begin{matrix} \hat{m_{b}} : \hat{M} \to \hat{M} . \end{matrix}$ $$\begin{array}{} \displaystyle \widehat{m_b}:\widehat M \to \widehat M. \end{array}$$

Note that since a ∈ L^∞(Ω) is strictly positive, $\begin{matrix} \hat{m_{1 / a}} = {\hat{m_{a}}}^{- 1} . \end{matrix}$ $\begin{array}{} \displaystyle \widehat{m_{1/a}}=\widehat{m_a}^{-1}. \end{array}$

Proposition 21

The following properties hold:

$\begin{matrix} \hat{m_{b}} R = R m_{b} . \end{matrix}$ $\begin{array}{} \displaystyle \widehat{m_b}\mathrm R=\mathrm R m_b. \end{array}$

$\begin{matrix} (\hat{m_{1 / a}} S, S)_{c} \geq 0 \end{matrix}$ $\begin{array}{} \displaystyle (\widehat{m_{1/a}}\mathrm S,\mathrm S)_c\ge 0 \end{array}$for all S ∈ M̂.

( S, RE)_c = $\begin{matrix} (\hat{C_{d i f f}} \hat{m_{a}} S, E)_{Ω} \end{matrix}$ $\begin{array}{} \displaystyle (\widehat{\mathrm {C_{diff}}}\widehat{m_a}\mathrm S,\mathrm E)_\Omega \end{array}$for all S ∈ M̂and E ∈ 𝕃²(Ω).

$\begin{matrix} m_{a} \hat{C_{d i f f}} = \hat{C_{d i f f}} \hat{m_{a}} . \end{matrix}$ $\begin{array}{} \displaystyle m_a\widehat{\mathrm {C_{diff}}}=\widehat{\mathrm {C_{diff}}}\widehat{m_a}. \end{array}$

Proof

It is simple to see that Π m_b = m_b|_M Π and therefore

$\begin{matrix} \hat{m_{b}} R = \hat{m_{b}} I Π = I m_{b} |_{M} Π = R m_{b}, \end{matrix}$ $$\begin{array}{} \displaystyle \widehat{m_b}\mathrm R= \widehat{m_b}\mathrm I\,\Pi=\mathrm I\, m_b|_M\Pi=\mathrm R m_b, \end{array}$$

which proves (a). If S ∈ M and note that

$\begin{matrix} (\hat{m_{1 / a}} S, S)_{c} = (a^{- 1} S, S)_{c} = (C_{d i f f} S, S)_{Ω} \geq 0, \end{matrix}$ $$\begin{array}{} \displaystyle (\widehat{m_{1/a}}\mathrm S,\mathrm S)_c =(a^{-1}\mathrm S,\mathrm S)_c=(\mathrm {C_{diff}}\mathrm S,\mathrm S)_\Omega \ge 0, \end{array}$$

then (b) follows by density. Finally, if S ∈ M, we have

$\begin{matrix} (S, R E)_{c} = (S, Π E)_{c} = (C_{d i f f} m_{a} S, Π E)_{Ω} = (C_{d i f f} m_{a} S, E)_{Ω} \forall E \in L^{2} (Ω), \end{matrix}$ $$\begin{array}{} \displaystyle (\mathrm S,\mathrm R\mathrm E)_c =(\mathrm S,\Pi\mathrm E)_c =(\mathrm {C_{diff}} m_a\mathrm S,\Pi\mathrm E)_\Omega =(\mathrm {C_{diff}} m_a\mathrm S,\mathrm E)_\Omega \quad \forall \mathrm E\in \mathbb{L}^2(\Omega), \end{array}$$

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 C_diff is strictly positive in a subdomain and vanishes in the complement. Assume that there are two open sets Ω₁ and Ω₂ such that

$\begin{matrix} Ω_{1} \cap Ω_{2} = \emptyset, \bar{Ω} = \bar{Ω_{1}} \cup \bar{Ω_{2}} \end{matrix}$ $$\begin{array}{} \displaystyle \Omega_1\cap\Omega_2=\emptyset, \qquad \overline\Omega=\overline{\Omega_1}\cup\overline{\Omega_2} \end{array}$$

and we have

$\begin{matrix} C_{d i f f} M : M \geq c_{0} χ_{Ω_{1}} M : M a.e. \forall M \in R_{sym}^{d \times d}, \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm {C_{diff}}\mathrm M \,:\, \mathrm M \ge c_0 \chi_{\Omega_1} \mathrm M:\mathrm M \qquad \mbox{a.e.} \quad \forall \mathrm M\in \mathbb R^{d\times d}_\text{sym}, \end{array}$$(35a)

$\begin{matrix} C_{d i f f} M = 0 & a.e. in Ω_{2} . \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm {C_{diff}}\mathrm M =0 && \mbox{a.e. in} \Omega_2. \end{array}$$(35b)

We can always write

$\begin{matrix} T = χ_{Ω_{1}} T + χ_{Ω_{2}} T \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm T=\chi_{\Omega_1} \mathrm T+\chi_{\Omega_2} \mathrm T \end{array}$$

and note that by (35b)

$\begin{matrix} C_{d i f f} T = χ_{Ω_{1}} C_{d i f f} T = C_{d i f f} (χ_{Ω_{1}} T) . \end{matrix}$ $$\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}$$(36)

If C_diffT = 0, then χ_Ω₁ T : T = 0 by (35a) and therefore T = 0 almost everywhere in Ω₁. Reciprocally, if T = 0 in Ω₁, then by (36), C_diffT = 0. We have thus proved that

$\begin{matrix} \ker C_{d i f f} = {T \in L^{2} (Ω) : T = 0 in Ω_{1}} \end{matrix}$ $$\begin{array}{} \displaystyle \ker\mathrm {C_{diff}}=\{ \mathrm T\in \mathbb{L}^2(\Omega)\,:\, \mathrm T=0 \mbox{ in} ~\Omega_1\} \end{array}$$(37)

and therefore

$\begin{matrix} M = {T \in L^{2} (Ω) : T = 0 in Ω_{2}} \equiv L^{2} (Ω_{1}) . \end{matrix}$ $$\begin{array}{} \displaystyle M=\{\mathrm T\in \mathbb{L}^2(\Omega)\,:\, \mathrm T=0\mbox{ in}~ \Omega_2\} \equiv \mathbb L^2(\Omega_1). \end{array}$$

However, now, due to (35b), we have

$\begin{matrix} | T |_{c} \approx ∥ T ∥_{Ω_{1}} \forall T \in M, \end{matrix}$ $$\begin{array}{} \displaystyle |\mathrm T|_c \approx \| \mathrm T\|_{\Omega_1} \qquad \forall \mathrm T\in M, \end{array}$$

and therefore M̂ = M. In this case R : 𝕃²(Ω) → M̂ is just the restriction to Ω₁ of matrix-valued functions defined on Ω, $\begin{matrix} \hat{C_{d i f f}} \end{matrix}$ $\begin{array}{} \displaystyle \widehat{\mathrm {C_{diff}}} \end{array}$ is the restriction of the action of C_diff to functions defined on $\begin{matrix} L^{2} (Ω_{1}; R_{sym}^{d \times d}) \end{matrix}$ $\begin{array}{} \displaystyle L^2(\Omega_1;\mathbb R^{d\times d}_\text{sym}) \end{array}$ , and the same happens to multiplication operators $\begin{matrix} \hat{m_{b}} \end{matrix}$ $\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 (C_diff = 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 C₀-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{matrix} \dot{U} (t) = A U (t) + F (t), \end{matrix}$ $$\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{matrix} H := L^{2} (Ω) \times L^{2} (Ω) \times \hat{M}, \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal H:=\mathbf L^2(\Omega)\times \mathbb{L}^2(\Omega) \times \widehat M, \end{array}$$

endowed with the norm

$\begin{matrix} ∥ (u, E, S) ∥_{H}^{2} := (ρ u, u)_{Ω} + (C_{0} E, E)_{Ω} + | S |_{c}^{2} . \end{matrix}$ $$\begin{array}{} \displaystyle \|(\mathbf u,\mathrm E,\mathrm S)\|_{\mathcal H}^2 :=(\rho\,\mathbf u,\mathbf u)_\Omega +(\mathrm C_0\mathrm E,\mathrm E)_\Omega +|\mathrm S|_c^2. \end{array}$$

We then define the operator

$\begin{matrix} A (u, E, S) := (ρ^{- 1} d i v (C_{0} E + \hat{C_{d i f f}} S), ε (u), \hat{m_{1 / a}} (R ε (u) - S)), \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal A(\mathbf u,\mathrm E,\mathrm S) :=(\rho^{-1}\mathrm{div}\,(\mathrm C_0\mathrm E+\widehat{\mathrm {C_{diff}}}\mathrm S), \boldsymbol{\varepsilon}({\mathbf u}), \widehat{m_{1/a}} (\mathrm R\boldsymbol{\varepsilon}({\mathbf u})-\mathrm S)), \end{array}$$

with domain

$\begin{matrix} D (A) := H_{D}^{1} (Ω) \times ((E, S) \in L^{2} (Ω) \times \hat{M} : \begin{matrix} d i v (C_{0} E + \hat{C_{d i f f}} S) \in L^{2} (Ω) \\ γ_{N} (C_{0} E + \hat{C_{d i f f}} S) = 0 \end{matrix}) . \end{matrix}$ $$\begin{array}{} \displaystyle D(\mathcal A) := \mathbf H^1_D(\Omega) \times \left\{ (\mathrm E,\mathrm S)\in \mathbb{L}^2(\Omega)\times \widehat M\,:\, \begin{array}{l} \mathrm{div}(\mathrm C_0\mathrm E+\widehat{\mathrm {C_{diff}}}\mathrm S) \in \mathbf L^2(\Omega)\\ \gamma_N (\mathrm C_0\mathrm E+\widehat{\mathrm {C_{diff}}}\mathrm S)=0 \end{array}\right\}. \end{array}$$

Proposition 22

The following properties hold:

(𝓐U, U)_𝓗 ≤ 0 for allU ∈ D(𝓐).

The operatorD(𝓐) ∋ U ↦ U − 𝓐U ∈ 𝓗 is surjective.

Therefore 𝓐 is maximal dissipative and generates a strongly continuous contractive semigroup in 𝓗.

Proof

If U = (u, E, S) ∈ D(𝓐), then

$\begin{matrix} (A U, U)_{H} = (d i v (C_{0} E + \hat{C_{d i f f}} S), u)_{Ω} + (C_{0} ε (u), E)_{Ω} + ({\hat{m_{a}}}^{- 1} (R ε (u) - S), S)_{c} \\ = - (\hat{C_{d i f f}} S, ε (u))_{Ω} + ({\hat{m_{a}}}^{- 1} R ε (u), S)_{c} - ({\hat{m_{a}}}^{- 1} S, S)_{c} = - (\hat{m_{1 / a}} S, S)_{c} \leq 0, \end{matrix}$ $$\begin{array}{} \displaystyle (\mathcal A\,U,U)_{\mathcal H} = (\mathrm{div}\,(\mathrm C_0\mathrm E+\widehat{\mathrm {C_{diff}}}\mathrm S),\mathbf u)_\Omega +(\mathrm C_0\boldsymbol{\varepsilon}({\mathbf u}),\mathrm E)_\Omega +(\widehat{m_a}^{-1}(\mathrm R\boldsymbol{\varepsilon}({\mathbf u})-\mathrm S),\mathrm S)_c\\\displaystyle\qquad\qquad~~ =-(\widehat{\mathrm {C_{diff}}}\mathrm S,\boldsymbol{\varepsilon}({\mathbf u}))_\Omega +(\widehat{m_a}^{-1}\mathrm R\boldsymbol{\varepsilon}({\mathbf u}),\mathrm S)_c -(\widehat{m_a}^{-1}\mathrm S,\mathrm S)_c =-(\widehat{m_{1/a}}\mathrm S,\mathrm S)_c\le 0, \end{array}$$

where we have progressively applied the boundary condition $\begin{matrix} γ_{N} (C_{0} E + \hat{C_{d i f f}} S) = 0 \end{matrix}$ $\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

$\begin{matrix} u \in H_{D}^{1} (Ω), \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u\in \mathbf H^1_D(\Omega), \end{array}$$

$\begin{matrix} (ρ u, v)_{Ω} + ((C_{0} + C_{d i f f} m_{1 + a}^{- 1}) ε (u), ε v)_{Ω} = (ρ f, v)_{Ω} - (C_{0} F, ε v)_{Ω} \\ - (\hat{C_{d i f f}} {\hat{m_{1 + a}}}^{- 1} \hat{m_{a}} G, ε v)_{Ω} \forall v \in H_{D}^{1} (Ω), \end{matrix}$ $$\begin{array}{} \displaystyle (\rho\,\mathbf u,\mathbf v)_\Omega +((\mathrm C_0+\mathrm {C_{diff}} m_{1+a}^{-1})\boldsymbol{\varepsilon}({\mathbf u}),\boldsymbol{\varepsilon}{\mathbf v})_\Omega = (\rho\,\mathbf f,\mathbf v)_\Omega -(\mathrm C_0\mathrm F,\boldsymbol{\varepsilon}{\mathbf v})_\Omega\\ \displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad~ -(\widehat{\mathrm {C_{diff}}}\widehat{m_{1+a}}^{-1}\widehat{m_a}\mathrm G, \boldsymbol{\varepsilon}{\mathbf v})_\Omega \quad \forall\mathbf v\in \mathbf H^1_D(\Omega), \end{array}$$(38b)

and define (see Proposition 21(a))

$\begin{matrix} E := ε (u) + F, \\ S := {\hat{m_{1 + a}}}^{- 1} (R ε (u) + \hat{m_{a}} G) = R m_{1 + a}^{- 1} ε (u) + {\hat{m_{1 + a}}}^{- 1} \hat{m_{a}} G, \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm E:= \boldsymbol{\varepsilon}({\mathbf u})+\mathrm F,\\\displaystyle \mathrm S:= \widehat{m_{1+a}}^{-1}(\mathrm R\boldsymbol{\varepsilon}({\mathbf u})+\widehat{m_a}\mathrm G) =\mathrm R m_{1+a}^{-1} \boldsymbol{\varepsilon}({\mathbf u}) +\widehat{m_{1+a}}^{-1}\widehat{m_a}\mathrm G, \end{array}$$

so that $\begin{matrix} S + \hat{m_{a}} S = R ε (u) + \hat{m_{a}} G \end{matrix}$ $\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{matrix} S = {\hat{m_{a}}}^{- 1} (R ε (u) - S) + G . \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm S=\widehat{m_a}^{-1}(\mathrm R\boldsymbol{\varepsilon}({\mathbf u})-\mathrm S)+\mathrm G. \end{array}$$

At the same time, by (34) and Proposition 21(a), we have

$\begin{matrix} C_{d i f f} m_{1 + a}^{- 1} ε (u) + \hat{C_{d i f f}} {\hat{m_{1 + a}}}^{- 1} \hat{m_{a}} G = \hat{C_{d i f f}} (R m_{1 + a}^{- 1} ε (u) + {\hat{m_{1 + a}}}^{- 1} \hat{m_{a}} G) = \hat{C_{d i f f}} S, \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm {C_{diff}} m_{1+a}^{-1}\boldsymbol{\varepsilon}({\mathbf u}) +\widehat{\mathrm {C_{diff}}} \widehat{m_{1+a}}^{-1}\widehat{m_a}\mathrm G =\widehat{\mathrm {C_{diff}}} (\mathrm R m_{1+a}^{-1}\boldsymbol{\varepsilon}({\mathbf u}) +\widehat{m_{1+a}}^{-1}\widehat{m_a}\mathrm G) =\widehat{\mathrm {C_{diff}}}\mathrm S, \end{array}$$

so that (38) implies

$\begin{matrix} (ρ u, v)_{Ω} + (C_{0} E + \hat{C_{d i f f}} S, ε v)_{Ω} = (ρ f, v)_{Ω} \forall v \in H_{D}^{1} (Ω), \end{matrix}$ $$\begin{array}{} \displaystyle (\rho\,\mathbf u,\mathbf v)_\Omega +(\mathrm C_0\mathrm E+\widehat{\mathrm {C_{diff}}}\mathrm S,\boldsymbol{\varepsilon}{\mathbf v})_\Omega = (\rho\,\mathbf f,\mathbf v)_\Omega \qquad\forall\mathbf v\in \mathbf H^1_D(\Omega), \end{array}$$

which is equivalent to

$\begin{matrix} ρ u - d i v (C_{0} E + \hat{C_{d i f f}} S) = ρ f, \\ γ_{N} (C_{0} E + \hat{C_{d i f f}} S) = 0. \end{matrix}$ $$\begin{array}{} \displaystyle \rho\,\mathbf u -\mathrm{div}(\mathrm C_0\mathrm E+\widehat{\mathrm {C_{diff}}}\mathrm S) = \rho\,\mathbf f,\\\displaystyle\qquad~~~ \gamma_N(\mathrm C_0\mathrm E+\widehat{\mathrm {C_{diff}}}\mathrm S) =0. \end{array}$$

Summing up, we have (u, E, S) ∈ D(𝓐) and (u, E, S) = 𝓐(u, E, S) + (f, F, G).

Theorem 23

For $\begin{matrix} α \in W_{+}^{2, 1} (0, \infty; H^{1 / 2} (Γ_{D})), β \in W_{+}^{1, 1} (0, \infty; H^{- 1 / 2} (Γ_{N})) \end{matrix}$ $\begin{array}{} \displaystyle \boldsymbol\alpha\in W^{2,1}_+(0,\infty;\mathbf H^{1/2}(\Gamma_D)), ~\boldsymbol\beta\in W^{1,1}_+(0,\infty;\mathbf H^{-1/2}(\Gamma_N)) \end{array}$ , andf ∈ L¹(0, ∞; L²(Ω)), there exists a unique

$\begin{matrix} (u, E, S) \in C^{1} ([0, \infty); H), \end{matrix}$ $$\begin{array}{} \displaystyle (\mathbf u,\mathrm E,\mathrm S)\in\mathcal C^1([0,\infty); \mathcal H), \end{array}$$

such that

$\begin{matrix} ρ \dot{u} (t) = d i v (C_{0} E (t) + \hat{C_{d i f f}} S (t)) + f^{(- 1)} (t) t \geq 0, \end{matrix}$ $$\begin{array}{} \displaystyle \rho \dot{\mathbf u}(t) =\mathrm{div}(\mathrm C_0\mathrm E(t)+\widehat{\mathrm {C_{diff}}}\mathrm S(t)) +\mathbf f^{(-1)}(t) \qquad t\ge 0, \end{array}$$(39a)

$\begin{matrix} \dot{E} (t) = ε (u (t)) & t \geq 0, \end{matrix}$ $$\begin{array}{} \displaystyle \dot{\mathrm E}(t) =\boldsymbol{\varepsilon}({\mathbf u(t))} &&&&t\ge0, \end{array}$$(39b)

$\begin{matrix} S (t) + \hat{m_{a}} \dot{S} (t) = R ε (u (t)) & t \geq 0, \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm S(t)+\widehat{m_a}\dot{\mathrm S}(t) =\mathrm R\boldsymbol{\varepsilon}({\mathbf u(t))} &&&& t\ge 0, \end{array}$$(39c)

$\begin{matrix} γ_{D} u (t) & = α (t) & t \geq 0, \end{matrix}$ $$\begin{array}{} \displaystyle \gamma_D \mathbf u(t) &=\boldsymbol\alpha(t) &&&&t\ge 0, \end{array}$$(39d)

$\begin{matrix} γ_{N} (C_{0} E (t) + \hat{C_{d i f f}} S (t)) = β^{(- 1)} (t) & t \geq 0 \end{matrix}$ $$\begin{array}{} \displaystyle \gamma_N(\mathrm C_0\mathrm E(t)+\widehat{\mathrm {C_{diff}}}\mathrm S(t)) =\boldsymbol\beta^{(-1)}(t) &&&&t\ge 0 \end{array}$$(39e)

and

$\begin{matrix} u (0) = 0, E (0) = 0, S (0) = 0. \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u(0)=0, \qquad \mathrm E(0)=0, \qquad \mathrm S(0)=0. \end{array}$$(39f)

Proof

For each t ≥ 0 we solve the elliptic problem

$\begin{matrix} u_{n h} (t) \in H^{1} (Ω), \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u_{\mathrm{nh}}(t)\in \mathbf H^1(\Omega), \end{array}$$(40a)

$\begin{matrix} ρ u_{n h} (t) = d i v (C_{0} ε u_{n h} (t), \end{matrix}$ $$\begin{array}{} \displaystyle \rho\,\mathbf u_{\mathrm{nh}}(t) =\mathrm{div}(\mathrm C_0\boldsymbol{\varepsilon}{\mathbf u_{\mathrm{nh}}(t)}, \end{array}$$(40b)

$\begin{matrix} γ_{D} u_{n h} (t) = α (t), \end{matrix}$ $$\begin{array}{} \displaystyle \gamma_D\mathbf u_{\mathrm{nh}}(t)=\boldsymbol\alpha(t), \end{array}$$(40c)

$\begin{matrix} γ_{N} C_{0} ε u_{n h} (t) = β^{(- 1)} (t), \end{matrix}$ $$\begin{array}{} \displaystyle \gamma_N\mathrm C_0\boldsymbol{\varepsilon}{\mathbf u_{\mathrm{nh}}(t)} =\boldsymbol\beta^{(-1)}(t), \end{array}$$(40d)

and note that

$\begin{matrix} ∥ u_{n h}^{(ℓ)} (t) ∥_{Ω} + ∥ ε u_{n h}^{(ℓ)} (t) ∥_{Ω} ≲ ∥ α^{(ℓ)} (t) ∥_{1 / 2, Γ_{D}} + ∥ β^{(ℓ - 1)} (t) ∥_{- 1 / 2, Γ_{N}} ℓ = 0, 1, 2. \end{matrix}$ $$\begin{array}{} \displaystyle \| \mathbf u_{\mathrm{nh}}^{(\ell)}(t)\|_\Omega +\|\boldsymbol{\varepsilon}{\mathbf u_{\mathrm{nh}}^{(\ell)}(t)}\|_\Omega \lesssim \|\boldsymbol\alpha^{(\ell)}(t)\|_{1/2,\Gamma_D} +\|\boldsymbol\beta^{(\ell-1)}(t)\|_{-1/2,\Gamma_N} \qquad\ell=0,1,2. \end{array}$$(41)

In a second step, we define the function F : [0, ∞) → 𝓗

$\begin{matrix} F (t) := (u_{n h} (t) - {\dot{u}}_{n h} (t) + ρ^{- 1} f^{(- 1)} (t), ε u_{n h} (t) - {\dot{u}}_{n h} (t), {\hat{m_{a}}}^{- 1} R ε u_{n h} (t)) \end{matrix}$ $$\begin{array}{} \displaystyle F(t):=(\mathbf u_{\mathrm{nh}}(t)-\dot{\mathbf u}_{\mathrm{nh}}(t) +\rho^{-1}\mathbf f^{(-1)}(t),\quad \boldsymbol{\varepsilon}{\mathbf u_{\mathrm{nh}}(t)-\dot{\mathbf u}_{\mathrm{nh}}(t)},\quad \widehat{m_a}^{-1}\mathrm R\boldsymbol{\varepsilon}{\mathbf u_{\mathrm{nh}}(t)}) \end{array}$$

and note that $\begin{matrix} F \in W_{+}^{1, 1} (0, \infty; H) \end{matrix}$ $\begin{array}{} \displaystyle F\in W^{1,1}_+(0,\infty;\mathcal H) \end{array}$ and that

$\begin{matrix} ∥ F^{(ℓ)} (t) ∥_{H} ≲ \sum_{k = ℓ}^{ℓ + 1} (∥ α^{(k)} (t) ∥_{1 / 2, Γ_{D}} + ∥ β^{(k - 1)} (t) ∥_{- 1 / 2, Γ_{N}}) + ∥ f^{(ℓ - 1)} (t) ∥_{Ω} ℓ = 0, 1. \end{matrix}$ $$\begin{array}{} \displaystyle \|F^{(\ell)}(t)\|_{\mathcal H} \lesssim \sum\limits_{k=\ell}^{\ell+1}\left( \|\boldsymbol\alpha^{(k)}(t)\|_{1/2,\Gamma_D} +\|\boldsymbol\beta^{(k-1)}(t)\|_{-1/2,\Gamma_N}\right) +\|\mathbf f^{(\ell-1)}(t)\|_\Omega\qquad \ell=0,1. \end{array}$$(42)

We use this function to solve the non-homogeneous initial value problem

$\begin{matrix} {\dot{U}}_{0} (t) = A U_{0} (t) + F (t) U_{0} (0) = 0, \end{matrix}$ $$\begin{array}{} \displaystyle \dot U_0(t)=\mathcal AU_0(t)+F(t) \qquad U_0(0)=0, \end{array}$$(43)

which has a unique solution U₀ = (u₀, E₀, S₀) ∈ 𝓒¹([0, ∞);𝓗), admitting the bounds

$\begin{matrix} ∥ U_{0} (t) ∥_{H} ≲ \int_{0}^{t} ∥ F (τ) ∥_{H} d τ, ∥ {\dot{U}}_{0} (t) ∥_{H} ≲ \int_{0}^{t} ∥ \dot{F} (τ) ∥_{H} d τ . \end{matrix}$ $$\begin{array}{} \displaystyle \|U_0(t)\|_{\mathcal H}\lesssim \int_0^t \|F(\tau)\|_{\mathcal H}\mathrm d\tau, \qquad \|\dot U_0(t)\|_{\mathcal H}\lesssim \int_0^t \|\dot F(\tau)\|_{\mathcal H}\mathrm d\tau. \end{array}$$(44)

In the next step, we define the triple (u, E, S) ∈ 𝓒¹([0, ∞);𝓗) by

$\begin{matrix} u (t) := u_{0} (t) + u_{n h} (t), \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u(t) :=\mathbf u_0(t)+\mathbf u_{\mathrm{nh}}(t), \end{array}$$(45a)

$\begin{matrix} E (t) := E_{0} (t) + ε (u_{n h} (t)), \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm E(t) :=\mathrm E_0(t)+\boldsymbol{\varepsilon}({\mathbf u_{\mathrm{nh}}(t)}), \end{array}$$(45b)

$\begin{matrix} S (t) := S_{0} (t) \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm S(t) :=\mathrm S_0(t) \end{array}$$(45c)

Taking into account the definition of F, equations (40) and (43) show that (u, E, S) satisfy (39).

Corollary 24

If α, β, andfsatisfy the conditions of Theorem 23, (u, E, S) is the solution of (39), and we define

$\begin{matrix} σ (t) := C_{0} \dot{E} (t) + \hat{C_{diff}} \dot{S} (t), \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol{\sigma}(t):={\mathrm C}_0\dot{\mathrm E}(t)+\widehat{\mathrm C_{\text{diff}}}\dot{\mathrm S}(t), \end{array}$$(46)

the pair (u, σ) satisfies the bounds for allt ≥ 0

$\begin{matrix} ∥ u (t) ∥_{Ω} ≲ \sum_{k = 0}^{1} \int_{0}^{t} (∥ α^{(k)} (τ) ∥_{1 / 2, Γ_{D}} + ∥ β^{(k - 1)} (τ) ∥_{- 1 / 2, Γ_{N}}) d τ \\ + \int_{0}^{t} ∥ f^{(- 1)} (τ) ∥_{Ω} d τ, \\ ∥ ε (u (t)) ∥_{Ω} + ∥ σ (t) ∥_{Ω} ≲ \sum_{k = 1}^{2} \int_{0}^{t} (∥ α^{(k)} (τ) ∥_{1 / 2, Γ_{D}} + ∥ β^{(k - 1)} (τ) ∥_{- 1 / 2, Γ_{N}}) d τ \\ + \int_{0}^{t} ∥ f (τ) ∥_{Ω} d τ . \end{matrix}$ $$\begin{array}{} \displaystyle \qquad\qquad\qquad~\|\mathbf u(t)\|_\Omega \lesssim \sum\limits_{k=0}^{1}\int_0^t \left( \|\boldsymbol\alpha^{(k)}(\tau)\|_{1/2,\Gamma_D} +\|\boldsymbol\beta^{(k-1)}(\tau)\|_{-1/2,\Gamma_N}\right)\mathrm d\tau \\ \displaystyle\qquad\qquad\qquad\qquad\qquad\quad+\int_0^t \|\mathbf f^{(-1)}(\tau)\|_\Omega\mathrm d\tau,\\ \displaystyle~\|\boldsymbol{\varepsilon}({\mathbf u(t))}\|_\Omega + \|\boldsymbol\sigma(t)\|_\Omega \lesssim \sum\limits_{k=1}^{2}\int_0^t \left( \|\boldsymbol\alpha^{(k)}(\tau)\|_{1/2,\Gamma_D} +\|\boldsymbol\beta^{(k-1)}(\tau)\|_{-1/2,\Gamma_N}\right)\mathrm d\tau \\ \displaystyle\qquad\qquad\qquad\qquad\qquad\quad +\int_0^t\|\mathbf f(\tau)\|_\Omega\mathrm d\tau. \end{array}$$

As a consequenceσ(0) = 0.

Proof

Use the decomposition (45) and the estimates (41), (42), and (44).

Corollary 25

Ifα, β, andfsatisfy the conditions of Theorem 23, (u, E, S) is the solution of (39), and we define

$\begin{matrix} σ (t) := C_{0} \dot{E} (t) + \hat{C_{diff}} \dot{S} (t), \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol\sigma(t):=\mathrm C_0\dot{\mathrm E}(t)+\widehat{\mathrm C_\text{diff}}\dot{\mathrm S}(t), \end{array}$$

the pair

$\begin{matrix} (u, σ) \in (C^{1} ([0, \infty); L^{2} (Ω)) \cap C ([0, \infty); H^{1} (Ω))) \times C ([0, \infty); L^{2} (Ω)) \end{matrix}$ $$\begin{array}{} \displaystyle (\mathbf u,\boldsymbol\sigma)\in \left(\mathcal C^1([0,\infty);\mathbf L^2(\Omega)) \cap \mathcal C([0,\infty);\mathbf H^1(\Omega))\right) \times \mathcal C([0,\infty);\mathbb{L}^2(\Omega)) \end{array}$$(47)

satisfies the equations

$\begin{matrix} ρ \ddot{u} (t) = d i v σ (t) + f (t) & a . e . - t, \end{matrix}$ $$\begin{array}{} \displaystyle \rho\,\ddot{\mathbf u}(t) =\mathrm{div}\,\boldsymbol\sigma(t)+\mathbf f(t) && a.e.-t, \end{array}$$(48a)

$\begin{matrix} σ (t) + a \dot{σ} (t) = C_{0} ε (u (t)) + (a C_{0} + C_{diff}) ε (\dot{u} (t)) & a . e . - t, \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol\sigma(t)+a\,\dot{\boldsymbol\sigma}(t) =\mathrm C_0\boldsymbol{\varepsilon}({\mathbf u(t))}+(a\mathrm C_0+{\mathrm C_{\text{diff}}})\boldsymbol{\varepsilon}({\dot{\mathbf u}(t))} & & a.e.-t, \end{array}$$(48b)

$\begin{matrix} γ_{D} u (t) = α (t) & t \geq 0, \end{matrix}$ $$\begin{array}{} \displaystyle \gamma_D\mathbf u(t) =\boldsymbol\alpha(t) && t\ge 0, \end{array}$$(48c)

$\begin{matrix} γ_{N} σ (t) = β (t) & a . e . - t, \end{matrix}$ $$\begin{array}{} \displaystyle \gamma_N\boldsymbol\sigma(t) =\boldsymbol\beta(t) && a.e.-t, \end{array}$$(48d)

with initial conditions

$\begin{matrix} u (0) = 0, \dot{u} (0) = 0, σ (0) = 0. \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u(0)=0, \quad \dot{\mathbf u}(0)=0, \quad \boldsymbol\sigma(0)=0. \end{array}$$(48e)

Proof

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, ∞) → H¹(Ω) follows from (39b). Note that (48c) is (39d). By (39a), we have that div(C₀E+ $\begin{matrix} \hat{C_{diff}} \end{matrix}$ $\begin{array}{} \displaystyle \widehat{\mathrm C_\text{diff}} \end{array}$S) ∈ 𝓒([0, ∞);L²(Ω)) and u̇(0) = 0, which was the missing initial condition (recall Corollary 24). We also have

$\begin{matrix} ρ \dot{u} (t) = d i v σ^{(- 1)} (t) + f^{(- 1)} (t), γ_{N} σ^{(- 1)} (t) = β^{(- 1)} (t), \end{matrix}$ $$\begin{array}{} \displaystyle \rho\dot{\mathbf u}(t)=\mathrm{div} \,\boldsymbol\sigma^{(-1)}(t) + \mathbf f^{(-1)}(t), \qquad \gamma_N\boldsymbol\sigma^{(-1)}(t)=\boldsymbol\beta^{(-1)}(t), \end{array}$$(49)

which are integrated versions of (48a) and (48d). Note next that by Proposition 21(d) and (34), we have

$\begin{matrix} \hat{C_{diff}} S (t) + a \hat{C_{diff}} \dot{S} (t) = \hat{C_{diff}} (S (t) + \hat{m_{a}} \dot{S} (t)) = \hat{C_{diff}} R ε (u (t)) = C_{diff} ε (u (t)) . \end{matrix}$ $$\begin{array}{} \displaystyle \widehat{\mathrm C_\text{diff}}\mathrm S(t)+a\widehat{\mathrm C_\text{diff}}\dot{\mathrm S}(t) =\widehat{\mathrm C_\text{diff}}(\mathrm S(t)+\widehat{m_a}\dot{\mathrm S}(t)) =\widehat{\mathrm C_\text{diff}}\mathrm R\boldsymbol{\varepsilon}({\mathbf u(t))} ={\mathrm C_{\text{diff}}}\boldsymbol{\varepsilon}({\mathbf u(t))}. \end{array}$$

We also have

$\begin{matrix} C_{0} E (t) + a C_{0} \dot{E} (t) = C_{0} ε (u)^{(- 1)} (t) + a C_{0} ε (u) (t), \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm C_0{\mathrm E}(t)+a\mathrm C_0\dot{\mathrm E}(t) =\mathrm C_0\boldsymbol{\varepsilon}({\mathbf u})^{(-1)}(t)+a\mathrm C_0\boldsymbol{\varepsilon}({{\mathbf u}})(t), \end{array}$$

and therefore

$\begin{matrix} σ^{(- 1)} (t) + a σ (t) = C_{0} ε (u)^{(- 1)} (t) + (a C_{0} + C_{diff}) ε (u (t)) . \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol\sigma^{(-1)}(t)+a\boldsymbol\sigma(t)= \mathrm C_0\boldsymbol{\varepsilon}({\mathbf u})^{(-1)}(t)+(a\mathrm C_0+{\mathrm C_{\text{diff}}})\boldsymbol{\varepsilon}({\mathbf u(t))}. \end{array}$$(50)

Equations (49) and (50) identify continuous functions of t taking values in L²(Ω), H^–1/2(Γ_N), and 𝕃²(Ω) 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

$\begin{matrix} ρ \ddot{u} - d i v σ \in L^{1} (0, \infty; L^{2} (Ω)), \\ γ_{N} σ \in L^{1} (0, \infty; H^{- 1 / 2} (Γ_{N})), \\ a \dot{σ} - (a C_{0} + C_{diff}) ε (\dot{u}) \in L^{1} (0, \infty; L^{2} (Ω)) . \end{matrix}$ $$\begin{array}{} \displaystyle \qquad\qquad\quad\rho\ddot{\mathbf u}-\mathrm{div}\,\boldsymbol\sigma \in L^1(0,\infty;\mathbf L^2(\Omega)),\\ \displaystyle\qquad\qquad\qquad\quad~~~\gamma_N\boldsymbol\sigma \in L^1(0,\infty;\mathbf H^{-1/2}(\Gamma_N)),\\ \displaystyle a\dot{\boldsymbol\sigma}-(a\mathrm C_0+{\mathrm C_{\text{diff}}})\boldsymbol{\varepsilon}({\dot{\mathbf u}}) \in L^1(0,\infty;\mathbb{L}^2(\Omega)). \end{array}$$

This finishes the proof.

Corollary 26

Let $\begin{matrix} α \in W_{+}^{3, 1} (0, \infty; H^{1 / 2} (Γ_{D})), β \in W_{+}^{2, 1} (0, \infty; H^{- 1 / 2} (Γ_{N})), \end{matrix}$ $\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 ∈ W^1,1(0, ∞;L²(Ω)). Additionally, let (u, E, S) solve (39) andσbe defined by(46). We have

$\begin{matrix} u \in C^{2} ([0, \infty); L^{2} (Ω)) \cap C^{1} ([0, \infty); H^{1} (Ω)), \\ σ \in C^{1} ([0, \infty); L^{2} (Ω)) \cap C ([0, \infty); H (d i v, Ω)) \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u \in \mathcal C^2([0,\infty);\mathbf L^2(\Omega)) \cap \mathcal C^1([0,\infty);\mathbf H^1(\Omega)),\\ \boldsymbol\sigma \in \mathcal C^1([0,\infty);\mathbb{L}^2(\Omega)) \cap \mathcal C([0,\infty); \mathbb H(\mathrm{div},\Omega)) \end{array}$$

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 𝓒²([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 C_diff = {0}, then M = 𝕃²(Ω) and R = I is just the canonical inclusion of 𝕃²(Ω) into its completion with respect to the norm $\begin{matrix} (a C_{diff} \cdot, \cdot)_{Ω}^{1 / 2} . \end{matrix}$ $\begin{array}{} \displaystyle (a{\mathrm C_{\text{diff}}}\,\cdot\,,\,\cdot\,)_\Omega^{1/2}. \end{array}$

When there exists c_diff > 0 such that

$\begin{matrix} (C_{diff} M) : M \geq c_{d i f f} ∥ M ∥^{2} \forall M \in R_{sym}^{d \times d}, \end{matrix}$ $$\begin{array}{} \displaystyle ({\mathrm C_{\text{diff}}}\,\mathrm M):\mathrm M\ge c_{\mathrm{diff}} \|\mathrm M\|^2 \qquad\forall\mathrm M\in \mathbb R^{d\times d}_\text{sym}, \end{array}$$(51)

there is no need to use the completion process since M = M̂ = 𝕃²(Ω) and then $\begin{matrix} \hat{C_{diff}} \end{matrix}$ $\begin{array}{} \displaystyle \widehat{\mathrm C_\text{diff}} \end{array}$ = C_diff, $\begin{matrix} \hat{m_{a}} = m_{a}, \end{matrix}$ $\begin{array}{} \displaystyle \widehat{m_a}=m_a, \end{array}$, and R is the identity operator. The space

$\begin{matrix} H := L^{2} (Ω; R^{d}) \times L^{2} (Ω) \times L^{2} (Ω), \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal H:=L^2(\Omega;\mathbb R^d)\times \mathbb{L}^2(\Omega) \times \mathbb{L}^2(\Omega), \end{array}$$

is now endowed with the norm

$\begin{matrix} ∥ (u, E, S) ∥_{H}^{2} := (ρ u, u)_{Ω} + (C_{0} E, E)_{Ω} + (a C_{diff} S, S)_{Ω}, \end{matrix}$ $$\begin{array}{} \displaystyle \|(\mathbf u,\mathrm E,\mathrm S)\|_{\mathcal H}^2 :=(\rho\,\mathbf u,\mathbf u)_\Omega +(\mathrm C_0\mathrm E,\mathrm E)_\Omega +(a\,{\mathrm C_{\text{diff}}}\mathrm S,\mathrm S)_\Omega, \end{array}$$

which is equivalent to the usual norm. This makes the analysis of this strictly diffusive viscoelastic problem much simpler.

When C_diff = 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

$\begin{matrix} ρ \dot{u} = d i v C_{0} E + f^{(- 1)}, \dot{E} = ε (u) \end{matrix}$ $$\begin{array}{} \displaystyle \rho\,\dot{\mathbf u}= \mathrm{div}\, \mathrm C_0 \mathrm E+\mathbf f^{(-1)}, \qquad \dot{\mathrm E}=\boldsymbol{\varepsilon}({\mathbf u}) \end{array}$$

of the classical linear elastic wave equation.

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

$\begin{matrix} u (0) = u_{0}, \dot{u} (0) = v_{0}, σ (0) = σ_{0}, \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u(0)=\mathbf u_0, \qquad \dot{\mathbf u}(0)=\mathbf v_0, \qquad \boldsymbol\sigma(0)=\boldsymbol\sigma_0, \end{array}$$

which would be the natural ones for the formulation (48). This is done by modifying the system (39), using initial conditions u(0) = u₀, E(0) = 0, S(0) = 0, and adding ρv₀ to the right-hand side of (39a) and G₀, with

$\begin{matrix} C_{diff} G_{0} = a σ_{0} - C_{1} ε (u_{0}) = a σ_{0} - (C_{diff} + a C_{0}) ε (u_{0}), \end{matrix}$ $$\begin{array}{} \displaystyle {\mathrm C_{\text{diff}}}\mathrm G_0=a\boldsymbol\sigma_0-\mathbf C_1\boldsymbol{\varepsilon}({\mathbf u_0})= a\boldsymbol\sigma_0-({\mathrm C_{\text{diff}}}+a\mathbf C_0)\boldsymbol{\varepsilon}({\mathbf u_0}), \end{array}$$

to the right-hand-side of (39c). The general case is much more complicated, given that the initial conditions for σ(0) must match C₀ε(u₀) 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.

10.1

Maxwell’s model

Maxwell’s model can be understood as the particular case of Zener’s model when C₀ = 0 and C_diff is strictly positive. However, its analysis is not included in the treatment given in Section 9, due to the fact that C₀ was used to define the norm of the space 𝓗. We thus start again, with a new space

$\begin{matrix} H := L^{2} (Ω) \times L^{2} (Ω), \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal H:=\mathbf L^2(\Omega)\times \mathbb{L}^2(\Omega), \end{array}$$

endowed with the norm

$\begin{matrix} ∥ (u, S) ∥_{H}^{2} := (ρ u, u)_{Ω} + (a C_{diff} S, S)_{Ω}, \end{matrix}$ $$\begin{array}{} \displaystyle \|(\mathbf u,\mathrm S)\|_{\mathcal H}^2 :=(\rho\mathbf u,\mathbf u)_\Omega +(a\,{\mathrm C_{\text{diff}}}\mathrm S,\mathrm S)_\Omega, \end{array}$$

which is equivalent to the usual norm. The domain of the operator

$\begin{matrix} A (u, S) := (ρ^{- 1} d i v C_{diff} S, m_{a}^{- 1} (ε (u) - S)) \end{matrix}$ $$\begin{array}{} \displaystyle \mathcal A(\mathbf u,S):= (\rho^{-1}\mathrm{div}\, {\mathrm C_{\text{diff}}}\mathrm S, m_a^{-1}(\boldsymbol{\varepsilon}({\mathbf u})-\mathrm S)) \end{array}$$

is now

$\begin{matrix} D (A) := H_{D}^{1} (Ω) \times (S \in L^{2} (Ω) : d i v C_{diff} S \in L^{2} (Ω), γ_{N} C_{diff} S = 0) . \end{matrix}$ $$\begin{array}{} \displaystyle D(\mathcal A):=\mathbf H^1_D(\Omega) \times \left\{ \mathrm S\in \mathbb{L}^2(\Omega)\,: \mathrm{div}\,{\mathrm C_{\text{diff}}}\mathrm S\in \mathbf L^2(\Omega), \gamma_N {\mathrm C_{\text{diff}}}\mathrm S=0\right\}. \end{array}$$

The operator 𝓐 is maximal dissipative. The proof of surjectivity of U ↦ U – 𝓐U starts with the solution of the coercive problem

$\begin{matrix} u \in H_{D}^{1} (Ω), \\ (ρ u, v)_{Ω} + (C_{diff} (1 + a)^{- 1} ε (u), ε (v))_{Ω} = (ρ f, v)_{Ω} \\ - ((a / (1 + a)) C_{diff} G, ε (v))_{Ω} \forall v \in H_{D}^{1} (Ω), \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u\in \mathbf H^1_D(\Omega),\\ \displaystyle(\rho\,\mathbf u,\mathbf v)_\Omega +({\mathrm C_{\text{diff}}} (1+a)^{-1}\boldsymbol{\varepsilon}({\mathbf u}),\boldsymbol{\varepsilon}({\mathbf v}))_\Omega = (\rho\,\mathbf f,\mathbf v)_\Omega\\ \displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad~ -((a/(1+a)) {\mathrm C_{\text{diff}}} \mathrm G,\boldsymbol{\varepsilon}({\mathbf v}))_\Omega \quad\forall\mathbf v\in \mathbf H^1_D(\Omega), \end{array}$$

for given (f, G) ∈ 𝓗. (Note that the strict positivity of C_diff is key for this argument to hold.) This is followed by the definition of

$\begin{matrix} S = (1 + a)^{- 1} (ε (u) + a G) . \end{matrix}$ $$\begin{array}{} \displaystyle \mathrm S=(1+a)^{-1}(\boldsymbol{\varepsilon}({\mathbf u})+a \mathrm G). \end{array}$$

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) ∈ 𝓒¹([0, ∞);𝓗) of the problem

$\begin{matrix} ρ \dot{u} (t) = d i v C_{diff} S (t) + f^{(- 1)} (t) t \geq 0, \\ S (t) + a \dot{S} (t) = ε (u (t)) t \geq 0, \\ γ_{D} u (t) = α (t) t \geq 0, \\ γ_{N} C_{diff} S (t) = β^{(- 1)} (t) t \geq 0, \end{matrix}$ $$\begin{array}{} \displaystyle \,\qquad\quad\rho \dot{\mathbf u}(t) =\mathrm{div} {\mathrm C_{\text{diff}}}\mathrm S(t) +\mathbf f^{(-1)}(t) \qquad t\ge 0,\\ \,\,\displaystyle\mathrm S(t)+a\,\dot{\mathrm S}(t) =\boldsymbol{\varepsilon}({\mathbf u(t))} \qquad\quad\qquad\qquad\quad t\ge 0,\\ \displaystyle\quad\quad\,\,\,\,\gamma_D \mathbf u(t) =\boldsymbol\alpha(t) \,\,\qquad\qquad\qquad\quad~~~~~~~ t\ge 0,\\ \displaystyle\quad\gamma_N {\mathrm C_{\text{diff}}}\mathrm S(t) =\boldsymbol\beta^{(-1)}(t)\qquad\qquad\qquad\quad~~ t\ge 0, \end{array}$$

with vanishing initial conditions. Introducing the stress tensor

$\begin{matrix} σ (t) := C_{diff} \dot{S} (t), \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol\sigma(t):={\mathrm C_{\text{diff}}}\dot{\mathrm S}(t), \end{array}$$

we have a solution of

$\begin{matrix} ρ \ddot{u} (t) = d i v σ (t) + f (t) a.e. - t, \\ σ (t) + a \dot{σ} (t) = C_{diff} ε (\dot{u} (t)) a.e. - t, \\ γ_{D} u (t) = α (t) t \geq 0, \\ γ_{N} σ (t) = β (t) a.e. - t, \end{matrix}$ $$\begin{array}{} \displaystyle \qquad\quad\rho\,\,\,\ddot{\mathbf u}(t) =\mathrm{div}\,\boldsymbol\sigma(t)+\mathbf f(t) \,\,\,\,\, \mbox{a.e.}-t,\\ \,\,\displaystyle\boldsymbol\sigma(t)+a\,\dot{\boldsymbol\sigma}(t) ={\mathrm C_{\text{diff}}}\boldsymbol{\varepsilon}({\dot{\mathbf u}(t))} \qquad\mbox{a.e.}-t,\\ \displaystyle\qquad\quad\gamma_D\mathbf u(t) =\boldsymbol\alpha(t) \qquad\qquad\quad \!t\ge 0,\\ \displaystyle\qquad\quad\gamma_N\boldsymbol\sigma(t) =\boldsymbol\beta(t) \qquad\qquad\quad\mbox{a.e.}-t, \end{array}$$

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 C₀ and C_diff 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{matrix} C_{0} + m_{1 + a}^{- 1} C_{diff} \end{matrix}$ $\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{matrix} (C_{0} \cdot, \cdot)_{Ω}^{1 / 2} \end{matrix}$ $\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.

10.2

Voigt’s model

The analysis of Voigt’s viscoelastic model (Zener’s model with a = 0, C₀ strictly positive and C_diff ≥ 0), including areas transitioning to classical linear elasticity, follows from a simple modification of the ideas of Section 9. In Voigt’s model C_diff = C₁ 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:

$\begin{matrix} H := L^{2} (Ω) \times L^{2} (Ω), \\ ∥ (u, E) ∥_{H}^{2} := (ρ u, u)_{Ω} + (C_{0} E, E)_{Ω}, \\ D (A) := ((u, E) \in H_{D}^{1} (Ω) \times L^{2} (Ω) : \begin{matrix} d i v (C_{0} E + C_{diff} ε (u)) \in L^{2} (Ω), \\ γ_{N} (C_{0} E + C_{diff} ε (u)) = 0 \end{matrix}), \\ A (u, E) := (ρ^{- 1} d i v (C_{0} E + C_{diff} ε (u)), ε (u)) . \end{matrix}$ $$\begin{array}{} \displaystyle \qquad\quad~\mathcal H := \mathbf L^2(\Omega)\times \mathbb L^2(\Omega),\\ \displaystyle\| (\mathbf u,\mathrm E)\|^2_{\mathcal H} :=(\rho\mathbf u,\mathbf u)_\Omega+(\mathbf C_0\mathrm E,\mathrm E)_\Omega,\\ \displaystyle\!\!\qquad D(\mathcal A) := \left\{ (\mathbf u,\mathrm E)\in \mathbf H^1_D(\Omega)\times \mathbb L^2(\Omega)\,:\, \begin{array}{l} \mathrm{div}(\mathbf C_0\mathrm E+{\mathrm C_{\text{diff}}} \boldsymbol{\varepsilon}({\mathbf u}))\in \mathbf L^2(\Omega),\\ \displaystyle\gamma_N(\mathbf C_0\mathrm E+{\mathrm C_{\text{diff}}} \boldsymbol{\varepsilon}({\mathbf u}))=0 \end{array}\right\},\\ \displaystyle\quad\!\!\mathcal A(\mathbf u,\mathrm E) :=(\rho^{-1} \mathrm{div}(\mathbf C_0\mathrm E+{\mathrm C_{\text{diff}}} \boldsymbol{\varepsilon}({\mathbf u})),\boldsymbol{\varepsilon}({\mathbf u})). \end{array}$$

A simple argument shows that

$\begin{matrix} (A (u, E), (u, E))_{H} = - (C_{diff} ε (u), ε (u))_{Ω} \leq 0 \forall (u, E) \in D (A) . \end{matrix}$ $$\begin{array}{} \displaystyle (\mathcal A(\mathbf u,\mathrm E),(\mathbf u,\mathrm E))_{\mathcal H} =-({\mathrm C_{\text{diff}}}\boldsymbol{\varepsilon}({\mathbf u}),\boldsymbol{\varepsilon}({\mathbf u}))_\Omega\le 0 \qquad\forall (\mathbf u,\mathrm E)\in D(\mathcal A). \end{array}$$

If we take (f, F) ∈ 𝓗, solve the coercive problem

$\begin{matrix} u \in H_{D}^{1} (Ω), \\ (ρ u, v)_{Ω} + ((C_{0} + C_{diff}) ε (u), ε (v))_{Ω} = (ρ f, v)_{Ω} - (C_{0} F, ε (v))_{Ω} \forall v \in H_{D}^{1} (Ω), \end{matrix}$ $$\begin{array}{} \displaystyle \mathbf u\in \mathbf H^1_D(\Omega),\\ \displaystyle(\rho\,\mathbf u,\mathbf v)_\Omega +((\mathrm C_0+{\mathrm C_{\text{diff}}} )\boldsymbol{\varepsilon}({\mathbf u}),\boldsymbol{\varepsilon}({\mathbf v}))_\Omega = (\rho\,\mathbf f,\mathbf v)_\Omega -(\mathrm C_0\mathrm F,\boldsymbol{\varepsilon}({\mathbf v}))_\Omega \quad \forall\mathbf v\in \mathbf H^1_D(\Omega), \end{array}$$

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) ∈ 𝓒¹([0, ∞);𝓗), vanishing at zero, and solving

$\begin{matrix} ρ \dot{u} (t) = d i v (C_{0} E (t) + C_{diff} ε (u (t))) + f^{(- 1)} (t) t \geq 0, \\ \dot{E} (t) = ε (u (t)) t \geq 0, \\ γ_{D} u (t) = α (t) t \geq 0, \\ γ_{N} (C_{0} E (t) + C_{diff} ε (u (t))) = β^{(- 1)} (t) t \geq 0. \end{matrix}$ $$\begin{array}{} \displaystyle \qquad\qquad\qquad\qquad~~\rho \dot{\mathbf u}(t) =\mathrm{div}(\mathrm C_0\mathrm E(t)+{\mathrm C_{\text{diff}}}\boldsymbol{\varepsilon}({\mathbf u(t))}) +\mathbf f^{(-1)}(t) \qquad\! t\ge 0,\\ \displaystyle\qquad\qquad\qquad\qquad\quad\dot{\mathrm E}(t) =\boldsymbol{\varepsilon}({\mathbf u(t))} \qquad\qquad\qquad\qquad\qquad\qquad\qquad~ t\ge0,\\ \displaystyle\qquad\qquad\qquad\qquad\gamma_D \mathbf u(t) =\boldsymbol\alpha(t)\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad~ t\ge 0,\\ \displaystyle\gamma_N(\mathrm C_0\mathrm E(t)+{\mathrm C_{\text{diff}}}\boldsymbol{\varepsilon}({\mathbf u(t))}) =\boldsymbol\beta^{(-1)}(t)\qquad\qquad\qquad\qquad\qquad\qquad\quad~~~ t\ge 0. \end{array}$$

The associated stress tensor is defined by

$\begin{matrix} σ (t) := C_{0} \dot{E} (t) + C_{diff} ε (\dot{u} (t)) = C_{0} ε (u (t)) + C_{diff} ε (\dot{u} (t)) \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol\sigma(t):=\mathbf C_0\dot{\mathrm E}(t)+{\mathrm C_{\text{diff}}}\boldsymbol{\varepsilon}({\dot{\mathbf u}(t))}= \mathbf C_0\boldsymbol{\varepsilon}({\mathbf u(t))}+{\mathrm C_{\text{diff}}}\boldsymbol{\varepsilon}({\dot{\mathbf u}(t))} \end{array}$$

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.

11.1

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

$\begin{matrix} ρ u_{t t} = σ_{x} & x \in [0, 1], t \in [0, 40], \end{matrix}$ $$\begin{array}{} \displaystyle \rho u_{tt} = \sigma_x \qquad & x\in [0,1], \quad t\in [0,40], \end{array}$$(52a)

$\begin{matrix} u (0, t) = g (t) & t \in [0, 40], \end{matrix}$ $$\begin{array}{} \displaystyle u(0,t) = g(t) && t\in [0,40], \end{array}$$(52b)

$\begin{matrix} σ (1, t) = 0 & t \in [0, 40], \end{matrix}$ $$\begin{array}{} \displaystyle \sigma(1,t) = 0&& t\in [0,40], \end{array}$$(52c)

$\begin{matrix} u (x, 0) = u_{t} (x, 0) = 0 & x \in [0, 1], \end{matrix}$ $$\begin{array}{} \displaystyle u(x,0) = u_t(x,0) = 0 &\qquad & x\in [0,1], \end{array}$$(52d)

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

The window function g(t) is smoothly changing in (0:5;1:5)∪(2:5;3:5), and constant on the rest of the domain. The function h(t) has a similar shape, although the upper plateau (forced normal stress) is longer.

$\begin{matrix} σ + a \partial_{t}^{ν} σ = C_{0} u_{x} + C_{1} \partial_{t}^{ν} u_{x}, \end{matrix}$ $$\begin{array}{} \displaystyle \sigma+ a\,\partial^\nu_t \sigma=\mathrm C_0 u_x+\mathrm C_1 \partial^\nu_t u_x, \end{array}$$(52e)

for parameters C₀, C₁, a and ν to be determined. For the discretization in space we use 𝓟₄ 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 C₁, 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 C₁ are chosen so that the parameter C_diff, which controls the diffusion, takes same values for the three different models, except for Maxwell, where C_diff = 0 results in a model that is identically zero. To avoid this while still getting comparable results, we use C₁ = 0.05 for the first value in the Maxwell model. In Figure 2, as we increase C₁, 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.

Effect of changing ν using the parameters given in Table 2.

Table 1

The parameters used to create the plots given in Figure 2. The first choice of C₁ for the Zener and Voigt models reduce the model to linear elasticity

	Zener	Maxwell	Voigt
C_o	1.5	0	1.5
C_l	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

Table 2

The parameters used to create the plots given in Figure 3.

	Zener	Maxwell	Voigt
C_o	1.5	0	1.5
C_l	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

$\begin{matrix} u (0, t) = 0 σ (1, t) = h (t), \end{matrix}$ $$\begin{array}{} \displaystyle u(0,t)=0 \qquad \sigma(1,t)=h(t), \end{array}$$

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.

Each plot shows the displacement of a 1D viscoelastic rod evaluated at the right endpoint, x = 1 for different values of ν. From left to right are the Zener, Maxwell, and Voigt models.

Table 3

The parameters used to create the plots given in Figure 4.

	Zener	Maxwell	Voigt
C_o	1	0	1
C_l	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.

11.2

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.

The Dirichlet boundary conditions used for the 2D spacetime simulations: single pulse and a periodic train of pulses, which will make the solution transition to time-harmonicity.

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.

Table 4

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
C_o	1.75	1.5	1.5	1.75 (x < 1/2), 1.5 (x ≥ 1/2)
C₁	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.

Table 5

The parameters used in the Maxwell model simulations to create the 2D spacetime plots shown in Figure 7.

	Elastic	Maxwell	Fractional Maxwell	Heterogeneous Domain
C_o	1.75	0	0	1.75(x < 1/2), 0 (x ≥ 1/2)
C₁	1.75	1.75	1.75	1.75
a	1	1	1	1
ν	1	1	0.3	1
ρ	10	10	10	10

Table 6

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
C_o	1.75	1.75	1.75	1.75
C₁	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

11.3

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

$\begin{matrix} \ddot{u} (t) = d i v σ (t) Ω \times [0, 50], \\ γ_{D} u (t) = 0.25 {(w (t), 0, 0)}^{⊤} Γ_{D} \times [0, 50], \\ γ_{N} σ (t) = 0 Γ_{N} \times [0, 50], \end{matrix}$ $$\begin{array}{} \displaystyle \quad\,\ddot{{\mathbf u}}(t) = \mathrm{div}\,\boldsymbol\sigma(t) \qquad\qquad\qquad\qquad\qquad\Omega\times[0,50],\\ \displaystyle\gamma_D{\mathbf u}(t) = 0.25 \left( w(t), 0, 0 \right)^\top \qquad\qquad\qquad~~\Gamma_D\times[0,50],\\ \displaystyle\gamma_N\boldsymbol \sigma(t) = 0 \qquad\qquad\qquad\qquad\qquad\qquad~~~\Gamma_N\times[0,50], \end{array}$$

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

$\begin{matrix} σ + \partial^{ν} σ = 2 ε (u) + (\nabla \cdot u) I + 5 (2 ε (\partial^{ν} u) + (\nabla \cdot \partial^{ν} u) I), \end{matrix}$ $$\begin{array}{} \displaystyle \boldsymbol\sigma + \partial^\nu\boldsymbol\sigma= 2\boldsymbol{\varepsilon}({\mathbf u})+(\nabla\cdot\mathbf u) {\mathbf I} +5 \left( 2\boldsymbol{\varepsilon}({\partial^\nu \mathbf u})+(\nabla\cdot\partial^\nu\mathbf u)\,\mathbf I\right), \end{array}$$

where

$\begin{matrix} ν := (\begin{matrix} 0, y \in [0, 5), \\ 1, y \in [5, 10] . \end{matrix}) \end{matrix}$ $$\begin{array}{} \displaystyle \nu := \left\{\begin{array}{c} 0,\quad y \in [0,5),\\[5pt] 1,\quad y\in [5,10]. \end{array}\right. \end{array}$$

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 𝓟₂ 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.

Snapshots for the 3D simulation showing the norm of the stress. From left to right, then from top to bottom, time-step = 9, 30, 55, 70, 100, 200, 350, 500.

The results of the same simulation as Figure 9, without the colormap. From left to right, then from top to bottom, time-step = 9, 30, 55, 70, 100, 200, 350, 500.

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.

eISSN:: 2444-8656
Language:: English

Publication timeframe:: 2 times per year
Journal Subjects:: Life Sciences, other, Mathematics, Applied Mathematics, General Mathematics, Physics

Journal RSS Feed

Analysis of models for viscoelastic wave propagation

Published Online: Oct 03, 2018

Page range: 55 - 96

Received: Dec 01, 2017

Accepted: Mar 09, 2018

DOI: https://doi.org/10.21042/AMNS.2018.1.00006

Keywords
Hyperbolic PDE, viscoelasticity, Laplace transforms, semigroups of operators, stability analysis, fractional derivatives

© 2018 Thomas S. Brown et al., published by Sciendo

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 License.

Fig. 1

Fig. 2

Fig. 3

Fig. 4

Fig. 5

Fig. 6

Fig. 7

Fig. 8

Fig. 9

Fig. 10

Analysis of models for viscoelastic wave propagation

Published Online: Oct 03, 2018

Page range: 55 - 96

Received: Dec 01, 2017

Accepted: Mar 09, 2018

DOI: https://doi.org/10.21042/AMNS.2018.1.00006

KeywordsHyperbolic PDE, viscoelasticity, Laplace transforms, semigroups of operators, stability analysis, fractional derivatives

© 2018 Thomas S. Brown et al., published by Sciendo

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 License.

Fig. 1

Fig. 2

Fig. 3

Fig. 4

Fig. 5

Fig. 6

Fig. 7

Fig. 8

Fig. 9

Fig. 10

Keywords
Hyperbolic PDE, viscoelasticity, Laplace transforms, semigroups of operators, stability analysis, fractional derivatives