Energetic boundary element method for accurate solution of damped waves hard scattering problems

The paper deals with the numerical solution of 2D wave propagation exterior problems including viscous and material damping coefficients and equipped by Neumann boundary condition, hence modeling the hard scattering of damped waves. The differential problem, which includes, besides diffusion, advection and reaction terms, is written as a space–time boundary integral equation (BIE) whose kernel is given by the hypersingular fundamental solution of the 2D damped waves operator. The resulting BIE is solved by a modified Energetic Boundary Element Method, where a suitable kernel treatment is introduced for the evaluation of the discretization linear system matrix entries represented by space–time quadruple integrals with hypersingular kernel in space variables. A wide variety of numerical results, obtained varying both damping coefficients and discretization parameters, is presented and shows accuracy and stability of the proposed technique, confirming what was theoretically proved for the simpler undamped case. Post-processing phase is also taken into account, giving the approximate solution of the exterior differential problem involving damped waves propagation around disconnected obstacles and bounded domains.

energy dissipation process, referred to as damping. Usually, the dissipation is generated by the interaction between the waves and the propagation medium. On the other side, for instance in mechanical systems, since damping has the effect of reducing the amplitude of vibrations, it is desirable to have some amount of damping in order to faster achieve an equilibrium configuration. Hence, damping is whether an unavoidable presence in physical reality or a desired characteristic in industrial design (see e.g., the applicative papers [1][2][3][4]).
For these reasons, it is interesting to consider damping terms into the modeling partial differential equations (PDEs): well-known examples are, for instance, the Klein-Gordon or the Telegraph equations (see e.g. [5] and references therein). For what concerns their approximate resolution, while the application of advanced numerical techniques, such as finite element (FEM) and finite difference (FDM) methods, is well established, the Boundary element method (BEM) analysis of dissipation through damped wave equation rewritten as a boundary integral equation (BIE) is a relatively recent research topic [6][7][8] and not yet sufficiently investigated. Note that BEMs are ideally suited for solving problems defined in unbounded domains, such as wave propagation scattering by finite obstacles, because they (i) rewrite the problem over the bounded obstacle and (ii) implicitly fulfill prescribed solution behavior at infinity. Let us further remark that for the numerical solution of this kind of problems, one needs consistent approximations and accurate simulations even on large time domains of analysis.
In principle, both frequency domain [9,10] and time domain [11,12] BEM can be used for hyperbolic initialboundary value problems. Space-time BEM has the advantage that it directly gives the unknown time-dependent quantities; in the case of unbounded domain problems, this is very beneficial. Let us note that Dirichlet absorbing boundary condition methods (see [13][14][15][16][17]) share the same feature. For space-time BEM, the construction of the BIEs, via boundary integral representation formula of the PDE solution, uses the fundamental solution of the hyperbolic differential operator at hand, if available, and jump relations [18]. The mathematical background of time-dependent BIEs is summarized by Costabel in [19].
This paper deals with the numerical solution of wave propagation problems exterior to open arcs in the plane, including viscous and material damping coefficients into the PDE, and equipped by Neumann boundary condition and homogeneous initial conditions. The numerical study of damped waves hard scattering model is a challenging task due to the presence of advection and reaction terms, and the adopted methodology generalizes what has been done in [20] for the much simpler undamped case. The differential problem is written as a space-time BIE whose kernel is given by the hypersingular fundamental solution of the 2D damped waves operator. The integral problem is then solved by a modified energetic BEM, appeared for the first time in its original form in [21]. Here, a suitable kernel treatment is introduced for the efficient evaluation of the discretization linear system matrix entries represented by space-time quadruple integrals with hypersingular kernel in space variables having no analytical time integration available, differently from the undamped case.
The paper is structured as follows: at first, we present the differential model problem on an unbounded 2D domain and its BIE energetic weak formulation. Section 3 briefly summarizes, for reader's convenience, the main features of Energetic BEM discretization which leads to a final linear system having a lower triangular block Toeplitz matrix. Section 4 highlights the new approach introduced for the evaluation of quadruple space-time integrals with hypersingular damped kernel, giving some details about the adopted quadrature schemes, while in Sect. 5, a wide variety of significant results are presented and discussed, showing, from a numerical point of view, stability and accuracy of the proposed technique, theoretically proved for the undamped case in [20]. Post-processing phase is also taken into account, giving the numerical solution of the exterior differential problem involving damped waves propagation around disconnected obstacles and bounded domains. Conclusions and ongoing research are reported in the last Section.

Model problem and its hypersingular BIE energetic weak formulation
In a 2D Cartesian reference system for space variables, we will consider the Neumann problem for the damped wave equation in a bounded time interval [0, T ], exterior to an obstacle given by an open arc Γ ⊂ R 2 : where c is the propagation velocity of a perturbation inside the domain, D ≥ 0 and P ≥ 0 are the viscous and material damping coefficients, respectively. Denoting by Γ − and Γ + the lower and upper faces of the obstacle, n stands for the normal unit vector defined in x ∈ Γ , oriented from Γ − to Γ + ; further, the datumq represents the opposite of the normal derivative of the incident wave along Γ , i.e.,q = −∂u I /∂n. In the acoustic framework, the exterior Neumann problem defines the scattering of a plane wave at a hard obstacle [22]. Remark When D = P = 0, the PDE (1) collapses to the classic wave equation and the considered model problem can be also conceived as the scattering problem by a crack in an unbounded elastic isotropic medium = R 2 \ Γ . As usual, the total displacement field can be represented as the sum of the incident field (the wave propagating without the crack) and the scattered field; this latter, under suitable invariance hypothesis, satisfies the Neumann scalar problem (1)-(3) (see e.g., [23]).
Let us consider the double-layer representation of the solution of (1)-(3) (in the particular case of an open arc, the Neumann boundary condition fixes also the jumps of the solution and therefore the only representation formula is the double-layer potential representation [19]): where the unknown density ϕ = [u] Γ represents the time history of the jump of u along Γ , is the fundamental solution of the 2D damped wave operator, with H [·] the Heaviside distribution, and r := r 2 = x − ξ 2 . Definition (5) is an extension of fundamental solutions found in [24][25][26], where only the viscous damping coefficient was taken into account, and it has been verified with the use of Mathematica software.
The case P > D 2 defines the so-called underdamping configuration, the case P < D 2 the overdamping configuration, while the ideal separation state P = D 2 , referred to the vanishing of both cos(·) and cosh(·) arguments, is called critical damping. Furthermore, if we consider D = P = 0, we recover the fundamental solution of the wave equation in 2D without damping, i.e., Note that the solution of the differential problem can be evaluated at any point outside the obstacle and at any time instant from (4), with a post-processing phase once the density function ϕ(x, t) has been recovered.
To this aim, applying a directional (normal) derivative w.r.t. x in (4), performing a limiting process for x tending to Γ and using the assigned Neumann boundary condition (3), we obtain the hypersingular space-time BIE in the unknown ϕ(x, t), which can be written with the compact notation The equivalence of (8) with the initially given differential model problem (1)-(3) has been stated for general time-dependent problems in [19]. In order to properly manage equation (8), we introduce the following results.

Proposition 1
Defining the auxiliary kernelG(r, t − τ ) as representation formula (4) can be rewritten explicitly for x ∈ R 2 \ Γ and t ∈ (0, T ] as Proof Let us consider the case P ≥ D 2 ; starting from the definition of the double-layer potential in the right-hand side of (4), we observe that Now, considering the integration over Γ × (0, t) of the above expression multiplied by ϕ(ξ , τ ), integrating in the sense of distributions the term containing the derivative with respect to τ and then expressing explicitly one finally deduces the thesis. When P < D 2 we can proceed in the same manner, remembering that cos(i ·) = cosh(·) and sin(i ·) = i sinh(·), being i the imaginary unit.

Proposition 2
The space-time hypersingular integral operator D admits the following explicit expression: Proof Applying ∂/∂n x to the integrand function in (10), after some analytical manipulations, similar to what has been shown in the proof of Proposition 1, (11) is straightforwardly deduced.
Note that the expressions in the right-hand sides of (10) and (11) generalize formulas given for the case of undamped waves (P = D = 0) and for c = 1 in [20]. Problem (8) is then set in a weak form. Its energetic weak formulation can be deduced observing that, multiplying the PDE (1) by u t , integrating over [0, T ] × (R 2 \ Γ ) and using integration by parts in space, one obtains that the energy of the solution u at the final time of analysis T , defined by can be rewritten as Hence, projecting (8) by means of test functions ψ, derived w.r.t. time and belonging to the same functional space of the unknown density ϕ, we can write the energetic weak problem: Remark The theoretical analysis of the quadratic form coming from the left-hand side of (14) was carried out for the undamped case (P = D = 0) in [20]: under suitable hypothesis, coercivity was proved in an infinitedimensional subspace H Δx of H 1 ([0, T ]; H 1/2 0 (Γ )) depending on a piece-wise linear semi-discretization in space variables. Fixing trivial damping parameters in (11)and indicating by D 0 the hypersingular integral operator so obtained for the undamped case, it holds [20] for the great amount of technical details). This property allowed us to deduce stability and convergence of the related space-time Galerkin approximate solution. For the more involved case of non-trivial damping coefficients, these properties are only conjectured, but they will be checked from a numerical point of view with intensive tests on several benchmarks discussed in Sect. 5.

Energetic BEM discretization
In the following, we will briefly summarize the main steps in the discretization of energetic weak problem (14).
For time discretization, we consider a uniform decomposition of the time interval [0, T ] with time step Δt = T /N Δt , N Δt ∈ N + generated by the N Δt+1 time-knots: t k = kΔt, k = 0, . . . , N Δt and we choose piece-wise linear time shape functions defined by where is the ramp function. For the space discretization, we introduce on the obstacle a mesh constituted by M Δx straight elements e 1 , . . . , e M Δx , with 2l i :=length(e i ) ≤ Δx, e i ∩ e j = ∅ if i = j and such that M Δx i=1 e i coincides with Γ if the obstacle is (piece-wise) linear, or is a suitable approximation of Γ , otherwise. The functional background compels one to choose space shape functions belonging to H 1 0 (Γ ); hence we use standard piece-wise linear polynomial boundary element functions w j (x), j = 1, . . . , M Δx − 1, related to the introduced mesh and vanishing at the endpoints of Γ , even if higher degree Lagrangian basis can be employed. Hence, the approximate solution of the problem at hand will be expressed as The Galerkin BEM discretization coming from energetic weak formulation (14) produces the linear system and matrix A has a block lower triangular Toeplitz structure, owing, respectively, to the choice (15) and the fact that involved kernels depend on the difference between t and τ . Each block has dimension M Δ x − 1, which for 2D problems is typically low. The solution of (17) is obtained by a block forward substitution, being the diagonal block non-singular. A FFT-based algorithm for the fast resolution of this kind of structured linear systems can be used to speed up the computation, as suggested in [27].

Handling hypersingular kernel space-time integration in matrix entries evaluation
In order to avoid a heavy notation in the explanation of the fundamental issues involved in the evaluation of the entries of matrix A, we consider, at first, the geometrical case of a straight obstacle Γ . In this setting, the BIE (8) reduces to (18) and, in (17), matrix entries coming from the energetic weak formulation of (18) are quadruple integrals of the form Specifying the choice made for time basis function, quadruple integral (19) results, after some analytical manipulations, in a linear combination of integrals of the form In the case of the undamped wave equation, the analysis of the space hypersingularity O(1/r 2 ) of the integrand function has been performed in detail in [20]; furthermore, the presence in the kernels of the Heaviside distribution , which represents the wavefront, and of the square root c 2 (t − τ ) 2 − r 2 can cause numerical troubles that in [28] have been solved by suitable splitting of the outer integral over Γ and using quadrature schemes which regularize integrand functions with mild singularities. We stress the fact that all the above-recalled space singularities appear after analytical integration in time variables. Unfortunately, these analytical time integrations are no more possible on the damped kernel but we expect for the problem at hand a singular behavior similar to the one just described. Therefore, we consider the first term of the expansion of G(r, t − τ ) with respect to damping parameters, centered in P = D = 0, i.e., the undamped kernel G 0 (r, t − τ ) defined in (6), and we subtract and add back this term to the damped fundamental solution. In this way, we can split (20) as After analytical double time integration in the second quadruple integral of (21), we finally obtain the following expression for (20): where and

Proposition 3
In (22), it holds: thus highlighting the hypersingular nature of space integration.
Proof The limit of F 2 (r, t h , t k ) can be evaluated straightforwardly from (24) and it turns out to be − (t h −t k ) 2 4π . For what concerns the limit of r F 1 (r, t h , t k ), we cannot perform, as already remarked, time integrations in (23) analytically; anyway, since we are interested in the behavior of the integrand function for vanishing r , we observe that in this case the time integration, which for r > 0 must handle a weak singularity due to the presence of [c 2 (t − τ ) 2 − r 2 ] −1/2 , shows a stronger singularity in the neighborhood of τ = t. Hence, if we expand the whole integrand function in (23) w.r.t. r = 0 and τ = t, we obtain, besides regular terms of polynomial type r (t − τ ) ι , ≥ 0, ι ≥ 0, which cannot develop outer space singularities, terms of the type r (t−τ ) ι+1 , ≥ 0, ι = 0, . . . , . These terms, after double time integration, show a space singular behavior at most of type log(r ); therefore, lim r →0 r F 1 (r, t h , t k ) = 0 and the thesis follows.
Remark The above result assures that the difference between G and G 0 in (23) removes outer integration space hypersingularities. These singularities stem only after time integration in the second quadruple integral in (21). Moreover, we will have to deal only with weakly singular numerical time integration in (23) for r > 0, since the case r = 0 has been treated by means of the above analytical considerations.
For the general case of an open arc Γ in the plane, we can proceed in the same manner, starting from the complete expression of the hypersingular integral operator in (11). Hence, we can finally state the following result.

Proposition 4
Using piece-wise linear time basis functions, matrix entries coming from the energetic weak formulation of (8) are of the form where F 1 (r, t h , t k ) and F 2 (r, t h , t k ) are defined in (23) and (24), respectively, and Moreover, , thus highlighting the hypersingular nature of space integration in (26).
Using the standard element by element technique, the evaluation of every double integral (22) (the same treatment is applied to (26)) is reduced, up to the constant −1/(Δt) 2 , to the assembling of local contributions of the type where r = r (s, z) andw i ,w j define one of the Lagrangian basis functions in local space variables over the elements e i and e j , respectively. Due to the hypersingularity O(1/r 2 ) for r → 0, the evaluation of double integrals of type (29) is troublesome when e i ≡ e j and when e i , e j are consecutive, and it has been performed similarly as in [28]. Anyway, for reader's convenience, in the following we describe the double integration over coincident elements of the boundary mesh, i.e., i = j, which is the the most difficult case since hypersingularity appears at any point of e i . Let us start with an analysis of the double integration domain in (29) for i = j. In this case, the distance between the source and the field points can be written as r = |s − z|. Due to the presence of the Heaviside function H [c(t h − t k ) − r ], the double integration domain is represented by the intersection between the square [0, 2 l i ] 2 and the strip defined by |s − z| < c(t h − t k ) where the Heaviside function is not trivial. Having set double integral (29) in this case becomes The numerical quadrature in the variable s has been performed possibly subdividing the outer interval of integration in correspondence to the abscissas  Hence, (30) will be eventually decomposed into the sum of double integrals of the form where Of course when no subdivision is needed, we will have to deal with only one double integral At this stage, exploiting (25), we rewrite double integral (31) as follows: For the numerical evaluation of I 1 we use: Outer integral regularization procedure for integrand functions with endpoints mild singularities, coupled with Gauss-Legendre rule, as introduced and analyzed in [29]. This procedure depends on parameters p, q which suitably push the quadrature nodes towards the lower and upper endpoints of the integration interval, respectively ( p = q = 1 means standard Gauss-Legendre rule); Inner integral interpolatory product rule [30] for |s − z| −1 kernel, with Gauss-Legendre quadrature nodes and modified Gauss-Legendre weights absorbing the strong singularity.
Fixing P = D = 0 and discretization parameters 2l i = 0.1, t h − t k = 0.15, c = 1, F 1 (r, t h , t k ) becomes trivial and it is possible to evaluate I 1 analytically. Therefore, Table 1 reports the relative errors in the numerical evaluation of the integral I 1 , varying p and q and increasing the number of nodes in the outer integration; the inner integration has been performed with 16-nodes product rule, assuring single precision accuracy. Symbol '−' means that the overall single precision accuracy in the evaluation of I 1 has been achieved. Further, Table 2 shows results for an example of evaluation of integral I 1 in the presence of damping. Let us observe that in this case we have to evaluate numerically F 1 (r, t h , t k ): time integration has been done by means of the above-cited regularization procedure [29] coupled with Gauss-Legendre rule, due to the presence of square root function c 2 (t h − t k ) 2 − r 2 . The inner integration of I 1 in z variable has been performed with 16-nodes product rule, as before. For the outer integration we have fixed p = q = 2. Table 2 reports the relative errors (evaluated w.r.t. a reference value obtained fixing 256 nodes in all quadrature rules) in the numerical integration of I 1 , varying the number of outer quadrature nodes and increasing the value of damping parameters. In this case, the stabilization worsens a little bit w.r.t. to the simpler case P = D = 0 numerically analyzed before.
For I 2 , after the inner analytical integration of the hypersingular kernel we have  For consecutive as well as disjoint boundary elements involved as double integration domain, the interest reader can proceed similarly as described in [28]. The overall accuracy of the adopted numerical integration schemes is based on the convergence properties of the above-cited basic quadrature rules and allows to obtain stable and convergent approximate BIE solutions, as we will show in the next Section.

Numerical results
In the following, we present an intensive numerical study of Energetic BEM applied to the analysis of 2D damped waves hard scattering by both connected and disconnected obstacles and considering not only open arcs but also bounded domains.
At first, we consider the model problem (1) We show the results obtained for two different functions f (·), chosen also in [20] for the known asymptotic behavior of the solution, which allow us to validate the space-time approximate solution on Γ . In fact, at the Authors' knowledge, no analytical space-time solution is available for the model problem taken into account. In particular, we will see how differently the solution of a scalar wave propagation problem behaves in the presence of damping terms w.r.t. the undamped case. Firstly, we consider a plane harmonic wave, i.e., a wave which becomes harmonic after a fixed time: where ω represents the frequency. In this case the solution has to become harmonic too, with the same period as the incident wave, i.e., 2π/ ω, where ω = ω/2. The fixed circular frequency ω = 8π is such that the wave length λ = 2π/ω is equal to a quarter the crack length. We choose a uniform decomposition of the crack Γ in 20 subintervals (Δx = 0.05) and we decompose the observation time interval [0, 5] in 100 equal parts (Δt = 0.05).    In Fig. 2, for θ = π 2 and P = 0, D = 0, we show the time harmonic behavior of the Crack Opening Displacement (COD) ϕ at x = 0.25 (left) and at x = 0.35 (right), obtained starting from the energetic weak formulation. Note that the solution becomes immediately not trivial in both points, since the incident wave strikes the whole crack simultaneously.
In order to verify that the period of ϕ coincides with the period of the incident wave, we show in Fig. 3 the approximated COD at time instants t = 1, 3, 5, separated by multiples of the time period; note that, after the first considered time instant, the two successive curves perfectly match each other.
In Figs. 4, 5 (left) we show analogous results for P = 0, D = 4. Note that the presence of the non-trivial viscous damping coefficient changes the profile of the approximate solution, reducing the amplitude of the oscillations.  In Fig. 5 (right), CODs at the above-considered time instants t = 1, 3, 5 are reported for P = 4, D = 0. The profiles are very similar to the undamped case, but, due to the presence of the non-trivial material damping coefficient, the approximate solution presents a slightly more amplified oscillating behavior and at t = 3, 5 the related graphs are not yet completely overlapped, as one can see in Fig. 3. In Table 3 we show the variation of ϕ(0.25, t), t = 1, 3, 5, for D = 0 and different values of parameter P. Now, let us consider an incident plane linear wave, i.e., we fix in (32) f (t) = 0.5 t H [t]. In this case, when t tends to infinity, the Neumann datum (32) tends toq θ (x) = 0.5 n x · (cos θ, sin θ) = 0.5 (0, 1) · (cos θ, sin θ) = 0.5 sin θ =:q θ , independent of time and constant, so we expect that the approximate transient solution ϕ(x, t) of BIE (7) on Γ will tend to the BIE solution related to a simpler elliptic PDE.
When P = 0, D ≥ 0, we can discard in (1) the term depending on P and those derived w.r.t. time, for t → ∞; we can therefore consider the following exterior Boundary Value Problem (BVP) for the Laplace equation: and the related BIE on Γ , whose analytical solution is explicitly known and it reads ϕ ∞ (x) = sin θ √ x(1 − x). Let us remark that the static solution remains the same for every value of D, when P = 0, once θ is fixed.
For an incident angle θ = π/3 and for discretization parameters fixed as Δx = 0.05 and Δt = 0.05, in Fig.  6 we show the approximate solution obtained by Energetic BEM, at the final time instant of analysis T = 5, for P = 0 and different values of the viscous damping parameter D = 0, 0.5, 1, 4 (overdamping configuration), together with the reference Laplace BIE solution. The higher the value of D, the higher the gap between transient and steady-state solutions, meaning that, when we are in the presence of growing viscosity, more time is needed to see the overlapping between the two corresponding plots, i.e., to reach the equilibrium, conceived as the solution of the stationary BIE.   Fig. 8, for θ = π/3 (top) and θ = π/6 (bottom): for short times, if the wave hits the obstacle with a non-perpendicular incident direction, different points of the obstacle are stressed at different time instants generating an asymmetric solution on Γ . When the wavefront has covered the whole crack, for growing times the transient solution converges to the symmetric static one, and therefore the asymmetry tends to disappear.
In Fig. 9, the surface of the approximate solution in space and time variables is plotted for P = 0, D = 1, θ = π/2 (top), and θ = π/6 (bottom), together with the corresponding surface representing the difference between transient and analytic static solutions, on the time interval [0, 6]. Also here there is the evidence that when the wavefront of the incident wave is not perpendicular to the obstacle, different points of Γ are stricken at different time instants at the beginning of the simulation.
In order to observe longtime behavior of the Energetic BEM solutions and to numerically check longtime stability of the energetic formulation (14), we choose a uniform decomposition of Γ in 10 subintervals (Δ x = 0.1) and enlarge the observation time interval, fixing T = 15 and Δt = 0.1. In Fig. 10, for θ = π/3, P = 0 and different values of D ≥ 0, the graphs report the time history of ϕ(·, t) − ϕ ∞ (·) L 1 (Γ ) , i.e., the time history of the difference in L 1 (Γ ) norm between the approximate transient solution and the analytical stationary one related to the Laplace BVP (34). The fastest convergence is related to the critical damping configuration; the error order, which for growing time becomes the same for all values of D, is of course related to the fixed discretization parameters.
In this respect, in Fig. 11, we show the time history of ϕ(·, t) − ϕ ∞ (·) L 1 (Γ ) on the time interval [0, 8] for P = 0, D = 0, 0.125, 0.25 and different discretization parameters, having fixed θ = π/3: the smaller Δx and Δt, the smaller the gap between the approximate transient solution and the equilibrium state, for growing t. The   Fig. 9 Surfaces of approximate solution and its difference with the static one, for P = 0, D = 1, θ = π/2 (top) and θ = π/6 (bottom)  [20] in the 2D undamped case and theoretically proved in [5] for 1D damped problems.   Further, even if the space-time analytical solution is not known, we have studied the space-time energy norm of the approximate solutions. Remembering (13), this norm is defined by In Table 4, we show the values of (35), having fixed [0, T ] = [0, 10], θ = π/3, P = 0, and D = 0, 0.125, 0.25, for diminishing discretization parameters. We observe that the energy norms are converging to a limit value which can be conceived as the energy norm of the exact space-time BIE solution. Let us note that limit values are similar because, apart from an initial different transient phase due to different viscous damping parameters, the limit stationary solution is the same for the three problems here considered. Moreover, they are reasonably lower for higher damping parameters. When P > 0 we can discard in (1) the terms derived w.r.t. time, for t → ∞; hence, for any value of D ≥ 0, we can consider the following exterior BVP for the Helmholtz equation:  Looking at the graphs reporting the time history of ϕ(·, t) − ϕ ∞ (·) L 1 (Γ ) , shown in Fig. 13 for D = 0, P ≥ 0, we observe the expected convergence of each approximate transient solution to the corresponding approximate stationary one obtained using the same space discretization parameter. In particular, curves reveal a growing number of oscillations for growing P > 0 (underdamping configuration), but the same global decay of the undamped case, which presents a monotone convergence after t = 6 time instant. Note that the oscillations are due to intersections between approximate transient solutions and corresponding approximate stationary ones; they can be interpreted as oscillations of diminishing amplitude of each damped solution around its own equilibrium configuration.
In Fig. 14, the surface of the approximate solution in space and time variables is plotted for P = 1, D = 0, θ = π/2 (left), and θ = π/6 (right), on the time interval [0, 6]. Also here there is the evidence that when the  Fig. 15, both for θ = π/3. In both cases, the approximate solutions, plotted on Γ at T = 5, obtained by Energetic BEM with discretization parameters fixed as Δx = 0.05 and Δt = 0.05, reveal optimal accuracy w.r.t. the corresponding stationary solutions, visible directly in Fig. 6 in the case D = 0, P = 0 and comparing Fig. 15 with the Helmholtz BIE solution reported in Fig. 12 for P = 1.
Similarly as before, in Fig. 16, the surface of the approximate solution in space and time variables is plotted for P = 1, D = 1, θ = π/2 (left) and θ = π/6 (right), on the time interval [0, 6]. Also in this critical damping configuration Energetic BEM produces stable solutions.   Fig. 21, we show the transient approximate solutions at the final time instant T of analysis, together with reference stationary solution, for β = 10 −4 (left) and for β = 1 (right). Note that for Γ 10 −4 , which is almost flat, the approximate solution profile is analogous to those obtained for the straight crack, while for β = 1 the approximate solution changes its behavior over the obstacle Γ 1 . At last, we present some results involving the total displacement field u(x, t) obtained by the superposition of the incident wave u I (x, t) and the reflected and diffracted waves caused by the presence of an obstacle in the plane. The temporal profile of the incident wave, from which we have deduced the Neumann datum (32) on Γ , is shown in Fig. 22 and it is similar to that one considered in [32,33]. In the first simulation, Energetic BEM is applied in the case of a disconnected obstacle. In particular, the incident plane wave strikes perpendicularly a breakwater obstacle, made by five disjoint aligned or parallel segments each of length 0.5. The observation time interval is [0, 3]. As uniform temporal discretization step we use Δ t = 0.1 and every segment is uniformly approximated by 5 boundary elements. The total recovered displacement in a square around the obstacle at time instants t = 0.2 + j 0.3, j = 0, . . . , 8 is presented in Fig. 23 for P = 0, D = 0. These results show how the plane wave reaches the obstacle and how Γ degenerates the wavefront. Reflected waves are evident in the lower half of the square. As time increases, the wavefront recovers and the scattering effect caused by the obstacle on the plane wave diminishes. As one can see, the adopted approximation technique furnishes satisfactory results. In Figs. 24, 25, we show analogous results for P = 0, D = 1 and P = 0, D = 2, respectively: the higher the value of D, the lower the effects of reflection and diffraction around the disconnected obstacle. The wavefront is, however, overall weakened.  The proposed methodology can handle also closed boundaries. In the following simulation, the same incident plane wave strikes a unitary circle. The observation time interval is [0, 4]. As uniform temporal discretization step we use Δ t = 0.1 and the boundary of the circle is uniformly approximated by 80 straight boundary elements. Several snapshots related to the total recovered displacement in a square around the plane convex domain are shown in Fig. 26 for P = D = 0. In Fig. 27 for P = 1, D = 0, the positive material damping coefficient amplifies the recovering of the wavefront at the top of the circle, producing circular waves that are not present in the undamped case. Further, in Fig. 28 related to the critical damping configuration P = 1, D = 1, the positive viscosity weakens all phenomena around the circle visible in the previous figure.

Conclusions
We have considered the numerical solution of 2D damped wave propagation exterior problems equipped by Neumann boundary condition, which model hard scattering phenomena. A modified version of Energetic BEM has been applied, in order to take into account the non-analytical time integrability of the hypersingular damped kernel here involved. The original version of this method was already considered in the case of undamped wave equations, revealing, also from theoretical point of view, its accuracy and stability. Numerical results obtained in the damped scenario confirm that these properties are maintained in presence of dissipation. Next steps will involve BEM-FEM coupling for this kind of model problems. and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.