Out-of-Plane Elastic Waves in 2D Models of Solids: A Case Study for a Nonlocal Discretization Scheme with Reduced Numerical Dispersion

The paper addresses the problem of numerical dispersion in simulations of wave propagation in solids. This characteristic of numerical models results from both spatial discretization and temporal discretization applied to carry out transient analyses. A denser mesh of degrees of freedom could be a straightforward solution to mitigate numerical dispersion, since it provides more advantageous relation between the model length scale and considered wavelengths. However, this approach also leads to higher computational effort. An alternative approach is the application of nonlocal discretization schemes, which employ a relatively sparse spatial distribution of nodes. Numerical analysis carried out to study the propagation of elastic waves in isotropic solid materials is demonstrated. Fourier-based nonlocal discretization for continuum mechanics is introduced for a two-dimensional model undergoing out-of-plane wave propagation. The results show gradual increase of the effectiveness of this approach while expanding the region of nonlocal interactions in the numerical model. A challenging case of high ratio between the model length scale and wavelength is investigated to present capability of the proposed approach. The elaborated discretization method also provides the perspective of accurate representation of any arbitrarily shaped dispersion relation based on physical properties of modelled materials.


Introduction
Numerical dispersion is a well-known phenomenon present in numerical simulations [1,2].This phenomenon appears due to spatial discretization and temporal discretization that are applied to model a body, preferably of a complicated geometry, and to simulate its dynamic response using transient analysis.Spatial discretization results from specific dimensions introduced by an average element in Finite Elements (FEs), a cell size in Cellular Automata (CA), a distance between particles in peridynamics and molecular dynamics, or a distance between degrees of freedom (DOFs) for Finite Difference (FD) schemes [3][4][5][6][7][8].Moreover, an integration time, which is assumed for either explicit or implicit algorithms used to solve problems in the time domain, is responsible for temporal discretization.Hence the introduction of both scales, that is, length and time, into the model inevitably creates limits regarding maximal wavenumbers and frequencies, which may be handled in simulations accurately.The closer the parameters of simulated waves to these limits are the higher the influence of numerical dispersion is observed.As a consequence, characteristic harmonic disturbances are present in time domain responses.These disturbances do not have any physical explanation.A denser mesh of DOFs and smaller scales of length and/or time, preferably with local approaches, are often considered as a straightforward solution to mitigate numerical dispersion.Although this approach often provides more advantageous relation between model properties and considered wavelengths, higher computational effort is the major problem.It is clear that alternative discretization methods are desirable.If these methods were based on a nonlocal problem formulation, efficient and accurate models could be obtained even for relatively sparse spatial distribution of nodes in these models.
There is variety of different approaches that are used to reduce the influence of numerical dispersion in computational mechanics.The majority of these approaches use higher order FD-based discretization schemes [9].Different types of nonlocal FD schemes have been investigated to study the limitation of these approaches with respect to numerical errors.FD schemes with a large stencil have been considered in [10].This stencil is then optimized to maximize resolution characteristics for various spatial truncation orders.An increase of the quality for calculations can be also achieved by consideration of staggered grid-based discretization technique to FD [9,11].A new explicit time integration scheme is proposed in [12] to reduce numerical dispersion and to remove high spurious frequencies.A correction technique for time and space schemes has been proposed in [7] to mitigate numerical dispersion.Other interesting approaches include a non-Cartesian grid-based discretization method [13,14], a nearly-analytic central difference method applied to both spatial and time domains [6,15], a high order (eighth-order) nearly-analytic discrete operator for spatial derivatives, and a second-order scheme for temporal derivatives [16][17][18].A nonlocal and high order homogenization for modelling of wave dispersion and energy dissipation in composite materials are shown in [19][20][21].Higher order partial differential wave equations are also applied to model physical dispersion in periodic media [21].A nonlocal model for simulation of in-plane wave propagation in composites with periodic structures is investigated in [22].Spectral analysis of the dispersive and dissipative characteristics has been carried out to analyze properties of numerical solution for a twodimensional (2D) advection-diffusion equation [23].
This paper proposes an alternative approach for models with reduced numerical dispersion.This approach, based on a nonlocal discretization scheme for a 2D continuum model of solid, employs the Fourier series, allowing for effective approximation of a nondispersive analytical solution of wave equation.The proposed scheme also introduces a characteristic length scale, which by its nature prevents from obtaining a model of an ideally nondispersive material.However, as found in the results, the influence of numerical dispersion is dramatically reduced for a gradually increasing horizon of nonlocal interactions.It is important to note that no mesh rearrangement, for example, increased nodal density, is needed to improve the quality of computations.The definitions of resultant stiffness coefficients for discrete springs, which were used to model elastic properties, are provided.The present study addresses a 2D case of out-ofplane wave propagation.An initial value analytical problem is used to verify the results.
The structure of the paper is as follows.Section 2 introduces the problem solution for a one-dimensional (1D) case, which is then complemented with the results of an initial value problem in Section 3. A nonlocal discretization scheme for a 2D case, used to analyze out-of-plane propagation, is presented in Section 4. Next, Section 5 verifies the results

A E 𝜌 u(x, t) x
Figure 1: 1D model of continuous rod.
using 2D numerical simulations.Finally, Section 6 summarizes the paper and draws major conclusions.

Reduced Numerical Dispersion for a 1D Model of Solid Using Fourier-Based Nonlocal Discretization
For the sake of clarity, a 1D case is first presented to introduce the developed approach for simulation of out-of-plane waves in 2D models.A 1D model of an infinite prismatic rod is considered, as shown in Figure 1.This model represents a continuous and nondispersive medium for propagating longitudinal waves.The material used to build the model is isotropic and homogeneous.The mass is uniformly distributed along the longitudinal direction .The model is characterized using the following geometric and material properties: , the cross-sectional area, , Young's modulus, and , the mass density.The longitudinal deformation along the rod is determined by function (, ), which captures the displacements for a specific coordinate  and a given time .When d' Alembert's principle and the constitutive relation for a linear material are used, the well-known partial differential equation of motion for an infinitesimally thin and transversely cut slice of the rod can be found as [24]  2  (, ) where  is the velocity of longitudinal wave.This velocity can be defined as The introduction of a plane wave solution into (1) leads to the dispersion relation for a nondispersive material as where  0 is the initial longitudinal displacement, j stands for the imaginary unit, and  and  are wavenumber and circular frequency, respectively.Next, a series of discrete massless springs, which link lumped mass elements, can be used to derive a nonlocal 1D discretization scheme, as illustrated in Figure 2.This approach is considered to provide a numerical representation for a continuum model, preferably with reduced numerical dispersion.The entire mechanical domain is discretized with a spatial resolution  that represents the model length scale.Thus, any discrete location along -axis can be found as   = .Similarly, the longitudinal displacement along the rod (, ) can be defined as (  , ) =   ().
Since a nonlocal problem formulation is assumed, the numerical model takes into account both the nearest neighbours forces and long-range interactions.All local and nonlocal reaction forces are taken into account with the stiffness coefficients   attached to the discrete springs.For a mass element localized at the th position (  ), the displacements for all neighbouring masses situated at the coordinates  ± =   ±  are described as  ± .A discrete mass element  is attached to every th DOF, as resulting from uniform distribution of mass in original discretized model.By the introduction of d' Alembert's principle, the equation of motion for th DOF can be found as The plane wave solutions for discrete cases, for both  = 0, that is, at the th DOF, and  ̸ = 0, that is, at the neighbouring locations determined with respect to an actual central node, take the following forms: Herein, it is important to note that perspective for reduction of the influence of spatial discretization on numerical dispersion is arbitrarily chosen to be solely investigated in the present work.Hence, an ideal and analytically determined partial derivative calculated with respect to time can be found as Next, the dispersion relation for a 1D nonlocal discretization scheme can be derived by the introduction of ( 6) and ( 7) into (5) to obtain It is clear that once (4) is considered, (8) describes a nondispersive case if The left-hand side infinite sum in the above equation can be expressed in the form of the Fourier series as As a consequence, for a limited domain  ≤ , the nonperiodic function present at the right-hand side of ( 9) can be approximated using (10) as Equation ( 11) allows for a nondispersive solution for a 1D case until the wavelength  is not greater than the minimum allowed value  MIN = 2.This condition naturally refers to the period captured by the applied Fourier series.Anyway, irrespectively from the above-mentioned criterion, longitudinal displacement, for example, for the initial value problem, is not possible with a better spatial resolution because of the discretization governed by the length scale .In other words, any further decrease of  MIN is not feasible for the actual mesh with given .
The solution of ( 11) leads to the Fourier coefficients which can be finally used to find the stiffness coefficients   by the substitution of ( 12) into (11) and then into (9): As derived from (10), the following condition must be satisfied additionally: Condition ( 14) verifies the correctness of the nondispersive solution and can be proved using the solution of the Basel problem.Equation ( 13) coincides with the findings presented in [25], where the coefficients of a 1D nonlocal FD scheme, in turn, are determined.The normalized values of the stiffness coefficients,   / 1 , are presented in Figure 3.Moreover, (13) confirms the classical two-sign definition of a continuous function defying a nonlocal elastic modulus [26].
At this point, it should be highlighted that (13) remains valid only for the model of an infinite discrete rod with an unlimited number of nonlocal interactions, that is, when  → ∞.In other words, (13) applies only if all possible interactions between all DOFs are allowed.Therefore, this model description leads to an ideal approximation of the function (/2) 2 , which is made feasible with the application of all components of the Fourier series, as determined in (11).Indeed, one can easily find that applying (13) to (8) for a finite range of nonlocal interactions, that is, when  →  for finite , and assuring  → 0, the circular frequency  converges to either 0 or √ 2 for even and odd , respectively.This behaviour cannot be acceptable for numerical solution of wave propagation problems.The lack of convergence for  → 0, irrespective of , which can be either odd or even, is the consequence of limiting the interaction region and cutting off the terms in (8).This happens due to the fact that, for any finite , the last term in the considered finite sum, that is,   , is not correctly calculated according to (13).However, all remaining terms in the limited sum, that is,   for  < , can be still kept unchanged and taken from the infinite series to provide the solution with reduced numerical dispersion, as shown below.The coefficients   for  < , found based on (13), assure better and better approximation of the expression (/2) 2 when  is gradually increased.
The remaining component of a model exhibiting a finite range of nonlocal interactions, that is, the stiffness coefficient   , can be found as follows.For  → 0 and finite  the dispersion relation represented by ( 8) can be rewritten using the form When the first two of the most significant regressors of the Taylor series for cosine, that is, a constant and a quadratic term, are used and the limit  → 0 is assured, (15) can be transferred to

)]
→ 0 (16) and then to Following ( 4), for any finite , the condition must be satisfied to assure the nondispersive case for (17).The left-hand side in ( 18) can be expanded as The introduction of ( 13) and ( 18) into ( 19) allows for the definition of the stiffness coefficient for the last term   as Finally, an approximated nonlocal discretization scheme for a 1D rod with limited interaction region is characterized by the following stiffness coefficients: The increase of the model quality (accuracy) is illustrated in Figure 4, which shows the dispersion curves found from (21) for different values of .
The interpretation of (21) leads to the conclusion that consistency of the provided solution exhibiting reduced numerical dispersion is maintained as long as the half of the values of originally calculated stiffness coefficient   (determined for an infinite rod when  → ∞) are taken into account for DOFs at the boundaries of the interaction region, that is, when  = , for finite .

Initial Value Problem Solution for a 1D Model of Solid Discretized Using a Nonlocal Scheme
The time domain solution for the initial value problem is presented in this section as a reference.This reference will be used to verify the proposed approach.A rod of the length  = 200 m and the cross-sectional area  = 1m 2 is considered.The rod is made of aluminum with the following material properties: Young's modulus  = 70MPa and the mass density  = 2100 kg/m 3 .The distance between discrete masses, that is, the model length scale, equals  = 1 m.Different ranges for nonlocal interactions, defined with the parameter , are considered in the study to show the reduction of numerical dispersion.The analytical solution of the problem is taken into account as the reference.An implicit time integration technique is applied for the 0.01 s with very rigorous conditions regarding wavenumber and wavelength; that is,  = 4 ( = /2) and  = 8 ( = /4).The initial displacements are shown in Figure 5.These displacements represent a 6-cycle series of narrow sin waves; that is,   = 6 in ( 22) with the amplitude  0 = 1m.Relatively high wavenumbers are chosen to show the performance of the application of the Fourierbased discretization scheme.Figure 6 provides the reference between the contributions of wavenumbers for the initial displacements for both cases investigated (i.e.,  = 4 and

𝜆 = 8𝑎) and the dispersion curves for the nonlocal discretization approach.
The dispersion curves shown in Figure 6, presented as the relation between wavenumber  and normalized wave velocity /  , allow for prediction of the model behaviour when subjected to propagation of waves of different lengths.The relations between  and /  are nonlinear, irrespective of the value of  and the width of the region of dominant wavenumbers for initial displacement.However, extending the region for nonlocal interactions increases the quality of numerical simulations since better convergence to the analytical solution is observed.The choice of  = 1 leads to underestimation of the wave velocity for the entire range of considered wavenumbers, which applies for both cases of the initial deformation, that is, for  = 4 and  = 8.Moreover, a significant scatter for the normalized velocities (with the values varying from 0.6 up to 1) and the curve's slope must inevitably result in strong wave dispersion, which is illustrated in Figure 7. Hence, subsequent harmonics contribute to waves that travel with considerably different velocities, which in turn impose the change of the shape of initial displacement along the rod as time passes.The opposite behaviour regarding the wave velocity is found for  = 2. Further increase of  leads to better and better approximation of nondispersive material and allows for more realistic results.Potential readers are kindly encouraged to analyse the results shown in Figure 7 for different wavelengths  and ranges of nonlocal interactions .
The plots in Figure 7 are captured for the time  = 9.87 ms which relates to the time-of-flight for the arbitrarily selected distance 57 m for the longitudinal wave propagating with the theoretical value of velocity   = 5774 m/s.The case of a sixperiod sin wave is considered to show performance of the nonlocal discretization scheme.

Reduced Numerical Dispersion for Out-of-Plane Waves Propagation in a 2D Model of Solid
The application of a nonlocal Fourier-based discretization scheme for a 2D model is used in this section to simulate propagation of out-of-plane waves.A 2D mesh of discretized masses  is introduced to model continuum solid medium, as shown in Figure 8.The mass elements are uniformly localized along -axis and -axis, following the spatial resolution .
Although the model looks similar to a model of a membrane, displacement behaviour is different.The governing equation does not include tension resulting from an external force as for membranes.Only out-of-plane displacements  , () are introduced, allowing for propagation of a transverse wave.Indeed, the assumption is that all lumped mass elements , lying in one plane, are connected with the springs    and    along the directions parallel to -axis and -axis.In fact these springs are modelled as they were localized vertically with respect to the plane, on which all masses are localized.Similarly, the velocity of a transverse wave depends on the mass density of the material  and its elastic properties, that is, the transverse elastic modulus , even though no shear is explicitly considered in the model.However, it should be noted that the model is valid only for infinitesimally small out-of-plane displacements since inplane motion of the discrete masses remains fixed and inplane stresses are neglected.Hence, the model is incapable of simulating longitudinal waves.The proposed nonlocal approach is used to check the applicability and efficiency of the Fourier-based discretization scheme for solving wave propagation problems.The most straightforward description for this 2D model of a solid is introduced.The following paragraph shows that this type of discretization scheme assures correct propagation of transverse waves.
Assuming that all lumped mass elements can move only along -axis, the equation of motion for the out-of-plane displacement  , (), and for limited interaction region, can be written as follows: Next, the plane wave solution given as can be applied to obtain The introduction of ( 25) into ( 23) yields where the theoretical values of wave velocities referring to the analytical solution of wave propagation problem are The area of the vertical cross sections  depends on the thickness ; that is, The analytical solution for a continuous case is found for the wave equation to be after the introduction of a plane wave solution given as  (, , ) =  0  j(x−) .
(30)   As a result, a cone shaped dispersion surface can be obtained as If the elastic moduli   and   differ, the dispersion surface becomes Making a reference to the previous section on a 1D case, the values of    and    stiffness coefficients, for an approximate nondispersive discretization scheme with limited interaction region, can be derived for a 2D structure made of orthotropic material as where (⋅) * relates to the springs aligned according to  direction.
Examples of the dispersion surfaces, with related contour plots found from (33), are presented in Figure 9.The results are shown for different cases of nonlocal interaction regions defined by the values of  = 1, 2, 5, 10.The case of an isotropic material is presented, where   =   .The improvement regarding mitigation of numerical dispersion is especially seen in the contour plots.The case exhibiting local interactions (i.e., for  = 1) features strong mesh anisotropy.This characteristic means that the velocity of wave strongly depends on the angular orientation of the travelling direction with respect to the axes of the mesh of lumped mass elements in a discrete model.The results also show that there are preferential directions for out-of-plane waves propagation.Moreover, an overall significant decrease of the wave velocity is observed irrespective of the direction of propagation, which is also illustrated in Figure 10.Increasing  to 2 leads to partial overestimation of the velocity for some intervals of wavenumbers.The mesh anisotropy is less visible.Similar to the previously studied 1D case, the further increment of the parameter  leads to better quality of modelling a nondispersive medium.The case of  = 5 allows for a good approximation of the nondispersive case for the analytical solution of a continuous 2D model.Examples of the results for transient analyses performed to investigate out-ofplane wave propagation for different values of  are shown in Figures 10,11

Initial Value Problem Solution for Out-of-Plane Waves Propagation in a Nonlocal 2D Model
Following the 1D case, the 2D plate model is used to investigate the reduction of numerical dispersion for different horizons in long-range interactions.A 100 m by 100 m plate with the thickness of 1 m is used in numerical simulations.The model is capable of modelling out-of-plane displacements only, as mentioned in Section 4. Therefore only a transverse wave can be simulated for given shear modulus  = 26.9GPa, which applies for an isotropic material, aluminum, where  =   =   .The mass density is  = 2100 kg/m 3 .The distance between discrete localizations for considered DOFs, that is, the length scale, is  = 1m.The range of nonlocal interactions is defined by the value of  = 1, 2, 5, and 10 for a two-direction, orthogonal discretization scheme, as shown in Figure 8.An implicit time integration scheme is   The results refer to the model displacements determined for the simulation time  = 8.9 ms, which relates to the timeof-flight for the arbitrarily selected distance 32 m for a shear wave propagating with the theoretical value of the velocity   = 3581 m/s.
The results shown in Figures 10,11, and 12 confirm the observations described in the previous section.The wave simulated for the  = 1 case propagates too slow and characterizes significant changes of its shape regarding both the wave front and cross-sectional views.The wavefront takes a rectangular shape with rounded corners instead of ideally circular corners.Similar to the 1D case, higher order harmonics appear, which significantly disturb the initial wave.These harmonics are all visible because they propagate with smaller velocity with respect to the fundamental wave component determined by the given wavelength.An increase of the parameter  from 1 to 5 results in the convergence regarding the theoretical value of velocity and rounded shape of the wavefront.However additional harmonics still remain.A good approximation of a nondispersive case is seen for the reference case with  = 50.The shape of the initial wave is kept and no significant contribution of higher harmonics within the central area of the model is seen any more.The velocity of propagation again converges to the theoretical value.

Figure 2 :
Figure 2: A nonlocal 1D discretization model for continuous and infinite rod.

Figure 3 :
Figure 3: Normalized stiffness coefficients for the nonlocal discretization scheme.

Figure 4 :
Figure 4: Dispersion curves for the Fourier-based approach.The curves displayed are plotted for  equal to 1, 2, 5, and 10, as indicated.The gradual increase of quality is observed for the growing range of nonlocal interactions (increased  value).

Figure 8 :
Figure 8: Discretization scheme for a 2D model used to analyze outof-plane displacements.

Wavenumber
frequency  ( rad/s) Wav enum ber  x ( rad/ m) W av en um be r y ( ra d/ m frequency  ( rad/s) Wa ven um ber  x ( rad /m ) W a v e n u m b e r y ( ra d /m ) frequency  ( rad/s) W av en um be r  x ( ra d/ m ) W av en um be r y ( ra d/ m ) Analytical N  x ( rad/m) Wavenumber  y ( rad/m) Analytical N = 5

Figure 9 :
Figure 9: Dispersion surfaces (a) and their contour representation (b) for circular frequency determined with the nonlocal Fourier-based approach applied for a 2D case for different ranges of nonlocal interactions defined by .

xxFigure 10 :
Figure 10: Out-of-plane wave propagation in a 2D nonlocal model for different interaction regions.

Figure 11 :
Figure 11: Reference case for out-of-plane wave propagation for  = 50.
, and 12.Here Figures 10 and 11 present the deformation of the entire model, whereas Figure 12 provides a cross-sectional view of the displacements along -axis.