Abstract
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.
1. 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–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–18]. A nonlocal and high order homogenization for modelling of wave dispersion and energy dissipation in composite materials are shown in [19–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 two-dimensional (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-of-plane 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 using 2D numerical simulations. Finally, Section 6 summarizes the paper and draws major conclusions.
2. 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]where is the velocity of longitudinal wave. This velocity can be defined as The introduction of a plane wave solutioninto (1) leads to the dispersion relation for a nondispersive material aswhere is the initial longitudinal displacement, 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 asThe plane wave solutions for discrete cases, for both , that is, at the th DOF, and , 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 asNext, the dispersion relation for a 1D nonlocal discretization scheme can be derived by the introduction of (6) and (7) into (5) to obtainIt 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 asAs a consequence, for a limited domain , the nonperiodic function present at the right-hand side of (9) can be approximated using (10) asEquation (11) allows for a nondispersive solution for a 1D case until the wavelength is not greater than the minimum allowed value . 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 is not feasible for the actual mesh with given .
The solution of (11) leads to the Fourier coefficientswhich 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, , 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 , 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 , the circular frequency converges to either or for even and odd , respectively. This behaviour cannot be acceptable for numerical solution of wave propagation problems. The lack of convergence for , 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 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 and finite the dispersion relation represented by (8) can be rewritten using the formWhen 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 is assured, (15) can be transferred to and then toFollowing (4), for any finite , the conditionmust be satisfied to assure the nondispersive case for (17). The left-hand side in (18) can be expanded asThe introduction of (13) and (18) into (19) allows for the definition of the stiffness coefficient for the last term asFinally, 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 .
3. 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 m and the cross-sectional area m2 is considered. The rod is made of aluminum with the following material properties: Young’s modulus MPa and the mass density kg/m3. The distance between discrete masses, that is, the model length scale, equals 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 s long transient analysis with the time step μs. The propagation of longitudinal wave is analyzed, based on the calculated axial displacements for all DOFs along the rod. The initial model deformation is arbitrarily set as with very rigorous conditions regarding wavenumber and wavelength; that is, () and (). The initial displacements are shown in Figure 5. These displacements represent a 6-cycle series of narrow sin waves; that is, in (22) with the amplitude m. Relatively high wavenumbers are chosen to show the performance of the application of the Fourier-based discretization scheme. Figure 6 provides the reference between the contributions of wavenumbers for the initial displacements for both cases investigated (i.e., and ) and the dispersion curves for the nonlocal discretization approach.
(a)
(b)
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 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 and . 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 . 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 ms which relates to the time-of-flight for the arbitrarily selected distance m for the longitudinal wave propagating with the theoretical value of velocity m/s. The case of a six-period sin wave is considered to show performance of the nonlocal discretization scheme.
4. 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 in-plane motion of the discrete masses remains fixed and in-plane 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 ascan be applied to obtainThe introduction of (25) into (23) yieldswhere the theoretical values of wave velocities referring to the analytical solution of wave propagation problem areThe 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 beafter the introduction of a plane wave solution given as As a result, a cone shaped dispersion surface can be obtained as If the elastic moduli and differ, the dispersion surface becomesMaking 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 aswhere 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 . 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 ) 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 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-of-plane wave propagation for different values of are shown in Figures 10, 11, and 12.
5. 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 m by m plate with the thickness of 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 GPa, which applies for an isotropic material, aluminum, where . The mass density is kg/m3. The distance between discrete localizations for considered DOFs, that is, the length scale, is m. The range of nonlocal interactions is defined by the value of , and for a two-direction, orthogonal discretization scheme, as shown in Figure 8. An implicit time integration scheme is also applied to carry out transient analysis. A 5-period planar circularly symmetric sin wave is applied for an initial value problem. The wavelength for the in-plane directions equals . The results of the numerical simulations are presented in Figures 10, 11, 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. The results refer to the model displacements determined for the simulation time ms, which relates to the time-of-flight for the arbitrarily selected distance m for a shear wave propagating with the theoretical value of the velocity m/s.
The results shown in Figures 10, 11, and 12 confirm the observations described in the previous section. The wave simulated for the 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 . 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.
6. Summary and Concluding Remarks
A new approach for a nonlocal discretization scheme is proposed to mitigate numerical dispersion for simulation of wave propagation in solid media. The Fourier-based approach is effectively used to identify the long-range interaction coefficients applied in the mesh of lumped mass elements. The propagation of longitudinal (for the 1D rod case) and transverse (for the 2D plate case) waves is studied based on the developed numerical approach and its applicability and efficiency are discussed. Dispersion curves/surfaces are investigated and the obtained solutions are compared with the solutions of the relevant initial value problems. For both 1D and 2D cases, the dispersion curves and surfaces rapidly converge to their analytical counterparts. Similarly, the solutions in the time domain coincide, making a reference to either analytical expression (in case of 1D problem formulation) or nearly ideal outcomes of numerical simulations when using a very wide range of nonlocal interactions, which unquestionably leads to the desired nondispersive case.
The results show that the increase of the region for nonlocal interactions, represented by the Fourier-based stiffness coefficients, leads to the reduction of numerical dispersion. The results also show that the increase of the quality of numerical simulation, for arbitrarily selected very rigorous wavenumbers, can be achieved without denser meshes. Numerical dispersion reduction is feasible using the nonlocal approach and without any changes made to mesh arrangements.
However, it is important to note that although the approach promotes the increase of the interaction region for better and better approximation of a nondispersive case, any unnecessary extension of the nonlocal neighbourhood leads to a higher computational effort, which could be unjustified. Therefore for each case study a compromise between proper mesh density and the Fourier-based nonlocal interaction terms should be considered to balance both effects, that is, the accuracy required and the computational efficiency acceptable.
The focus of the work presented was entirely on spatial sources of numerical dispersion, assuming the optimal time sampling. The latter was assured by a very dense resolution of samples in the time domain used in the time integration procedures to avoid errors that could originate from incorrect time discretization. The presented approach is used for simulation of out-of-plane wave propagation in 2D models of solids and allows for explicit convergence of a general solution to a nondispersive solution. This approach does not follow the classical approaches of higher order FDs. The solution introduces the stiffness coefficients used to represent elastic properties of modelled physically nondispersive material, based on a nonlocal formulation. Additionally, the present results provide an insight into the physics of the nonlocal theory of elasticity. This theory considers elastic moduli as continuous functions of distance between parts of deformable solid. The presented approach also allows for modelling any type of dispersion determined during real measurements for investigated material. Future work in this area should focus on wave propagation in three-dimensional anisotropic domains and on the effect of time discretization.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The work was supported by the Foundation for Polish Science under Programme FNP-WELCOME/2010-3/2. Massimo Ruzzene and Julian J. Rimoli acknowledge the partial support of the National Science Foundation under Grant CMMI-1130368.