Multi-Directional Viscous Damping Absorbing Boundary in Numerical Simulation of Elastic Wave Dynamic Response

: Numerical seismic wave ﬁeld simulation is essential for studying the dynamic responses in semi-inﬁnite space, and the absorbing boundary se�ing is critical for simulation accuracy. This study addresses spherical waves incident from the free boundary by applying dynamic equations and Rayleigh damping. A new multi-directional viscous damping absorbing boundary (MVDB) method is proposed based on regional a�enuation. An approximate formula for the damping value is established, which can achieve absorbing the boundary se�ing by only solving the mass damping coeﬃcients without increasing the absorbing region grid cells or depending on the spatial and temporal walking distance. The validity and stability of the proposed method are proven through numerical calculations with seismic sources incident from diﬀerent angles. Meanwhile, the key parameters aﬀecting the absorption of the MVDB are analyzed, and the best implementation scheme is provided. In order to meet the requirements of mediums with diﬀerent elastic parameters for boundary absorption and ensure the high eﬃciency of numerical calculations, the damping amplitude control coeﬃcients k can be set between 1.02 and 1.12, the thickness of the absorbing region L is set to 2–3 times of the wavelength of the incident transverse wave, and the thickness of the single absorbing layer is set to the size of the discrete mesh of the model ∆ l .


Introduction
Numerical seismic wave field simulation based on the elastic wave equation is important to study the kinematics and dynamics characteristics of seismic wave propagation in geologic bodies, which is widely used in seismic design [1], geologic survey [2], and non-destructive structure detection [3].Currently, the numerical seismic wave field simulation methods include the finite element method (FEM) [4,5], the finite difference method (FDM) [6,7], and the boundary element method (BDM) [8,9].The FEM has become the most commonly used method for numerical seismic wave field simulation in geotechnical engineering due to its advantages in discretizing the mesh and free boundary treatment of irregular geological structures [10,11], and a variety of commercial software programs have been developed based on the FEM, such as ABAQUS (h ps://www.3ds.com/products/simulia/abaqus),COMSOL (h ps://www.comsol.com/),and ANSYS (h ps://www.ansys.com/)(all accessed on 1 March 2024).During simulation, artificially truncating the numerical model is often necessary due to the limitation of computational and storage capacities, and the artificially truncated boundaries must be effectively dealt with so that the seismic wave does not reflect at the boundaries, which is key to numerical seismic wave field simulation.
Two main ideas have been proposed to deal with artificially truncated boundaries in numerical wave field simulations using FEM software.One is to deal with the seismic wave field reflection at the boundary through free-field stress equilibrium at the nodes of the artificial boundary, such as the viscous boundary (VB) [12], the viscous-spring boundary (VSB) [13], and the equivalent viscous-spring boundary (EVSB) [14].The other is introducing a enuation terms in the wave equations to absorb the stress and displacement fields at the boundary, such as the infinite element boundary (IEM) [15], the perfectly matched layer (PML) [16][17][18][19][20][21], and the incremental damping method (IDM) [22][23][24][25].In geological surveys and structural damage detection, seismic waves are usually input from the free boundary at the ground surface, and the existing boundary processing methods fail to deal effectively with the boundary reflection problem of the surface wave and body wave.
The VSB achieves the purpose of absorbing the reflected wave field at the boundary by se ing up damping units to approximately equalize the stresses on the boundary.Compared with the VB, the VSB consisting of damping and spring units exhibits superior stability and accuracy in effectively simulating elastic recovery within an infinite domain [26].The accuracy of its solution depends on the form of the approximate expression for the incident wave.The approximate expression for the incident wave of the VSB is mainly established by the stress-displacement constitutive relationship of the plane wave and is, therefore, suitable for endogenous seismic wave situations [27].As shown in Figure 1a, the seismic sources of these structures are set deep underground, and the seismic waves in the plane form are transmi ed from the underground to the structure through absorbing boundaries, making this method widely used in the analysis of dynamic soil-structure interaction [28].For example, it is used to study the kinetic energy effects generated by natural earthquakes on structures such as dams, tunnels, and coasts (see e.g., [26,[29][30][31]).However, in transient analysis, as shown in Figure 1b, for the spherical wave incident from the free boundary, due to the lack of reasonable standardization for the spring component parameters in the approximate expression of the spherical wave [32], the absorption effect of the VSB is unsatisfactory in the numerical simulation of geological exploration and structural damage detection in semi-infinite space.With an a enuation term introduced into the wave equation, the infinite element (IEM) uses the shape functions to describe the decay of the displacement amplitude and the propagation of different waves, and the accuracy of the IEM depends on the wave number and decay factor selected [33].They can be accurately determined only for a single type of wave but not for multiple types of waves.As such, spurious reflections will occur for the other types of waves passing through the FEM-IEM interface [34].The PML is one of the most popular boundary schemes in recent years, which includes Classical PML [18,19], M-PML [20,21], C-PML [35,36], and joint boundary schemes framed by the PML, such as IEM-PML [37] and VB-PML [38].The PML defines the a enuation process of the outgoing wave from the near-field boundary to infinity through stretching and a enuation functions.The absorption effect of the PML on outgoing waves depends highly on the value of the stretching function.Although there are many available stretching functions [39][40][41], there is a lack of research on the selection of the best parameters for actual use [37].Furthermore, the PML cannot be directly applied to commercial finite element software and requires programming to implement [42].Even though some scholars have applied it to ABAQUS [5], the process is not convenient as it needs to be realized through ABAQUS-UEL combined with programming in computer languages, such as Fortran and Python.
The IDM scheme currently includes two main forms: one is to derive the mass and stiffness damping coefficients from the frequency and a enuation factor of the seismic source, and the classical method is the Caughey absorbing layer method (CALM) proposed by Semblat et al. [22,23].However, it is not convenient to cycle through two damping coefficients for a continuous absorber layer.The other is to calculate the maximum damping coefficient and determine the damping coefficients for each layer according to the designed a enuation rules, e.g., the non-reflective boundary (NRB) proposed by Wang et al. [24,25].Compared to the CALM, the NRB only requires determining the mass damping coefficient, greatly improving efficiency.However, the NRB requires some numerical experience to determine the maximum damping coefficient and brings about the probability of destabilization for low-frequency pulses.
Grasping the propagation characteristics of seismic waves in a geological body is a prerequisite for geological exploration and structural damage detection.In order to effectively deal with the reflection of seismic waves input from a free boundary at an artificially truncated boundary and ensure the accuracy and efficiency of numerical seismic wave field simulation, this study proposes, based on Rayleigh damping and using the incremental damping method, a multi-directional viscous damping absorbing boundary (MVDB) for the finite element time domain analysis.The MVDB sets up a certain absorption area in the periphery of the calculation area with consistent physical properties and does not require solving the control equations.The boundary seismic wave field is absorbed only by changing the damping value of the absorption area.With a simple implementation process, stable calculation, and effective absorption of the boundary seismic wave field, the MVDB can be well applied to the finite element software numerical seismic wave field simulation in geological exploration and structural damage detection.

MVDB Implementation
Semblat [22,23] proposed the CALM based on Rayleigh damping and second-order Caughey damping.The CALM relates the damping to the generalized Maxwell's equations by determining the quality factor (Q) through a linear viscoelastic rheological model [43] and then calculating the minimum damping frequency ( ) and the mass damping coefficient (α) and stiffness damping coefficient (β) based on the relationship (see e.g., [23,44]) between Q 1 and the damping ratio ( ).With high solution stability and absorption effect, the CALM can serve as a simple alternative to the PML [45,46] and can be implemented in FEM software.Wang [24,25] examined the Lamb wave dispersion characteristics in thin plate structures and introduced the NRB suitable for FEM transient dynamic analysis.The NRB determines the maximum mass damping coefficient ( ) in the absorbing layer through numerical experience and then calculates the damping coefficient (α) in each absorbing layer according to the damping control equation (see e.g., [24]).The NRB can effectively a enuate the Lamb wave reflection at the boundary and is easily implemented in ABAQUS, like the CALM.As its implementation process is more convenient than that of the CALM, the numerical simulation cost is greatly shortened (i.e., the storage capacity of computers).The CALM and NRB effectively suppress outgoing wave reflection at the boundary in practical engineering, but they are sometimes not ideal.For example, the a enuation mechanism of the CALM is relatively complicated, and the computational accuracy depends on the quality factor (Q).The of the NRB depends on some numerical experience and may fail to absorb low-frequency pulses.This section introduces the basic ideas and implementation routes of the MVDB from the basic dynamics equations and demonstrates the computational results.
In the time integration scheme, the central difference method is used in this study to solve the dynamics equations.Under the basic assumptions of elasticity theory (homogeneity and linear elasticity), the force vectors at any node b in space at time t can be derived from the basic dynamics equations as follows [47,48]: where [M] is the mass matrix, [K] is the stiffness matrix, [D] is the damping matrix, and { ( )}, ̇( ) , and ̈( ) are the displacement components and the first-order and second-order derivatives of the displacement components, respectively.[ ] ̈( ) , [ ] ̇( ) , and [ ]{ ( )} are the nodal loads due to acceleration-induced mass, nodal loads due to velocity-induced damping, and nodal loads due to displacement-induced stiffness, respectively.In the center difference method, the velocity and acceleration can be approximated as follows using equal time steps ∆ = ∆ : Therefore, the equations of motion for any node on the absorbing layer at the moment can be obtained by bringing Equations ( 2) and (3) into (1), taking into account the damping as follows: The load on the b node can be minimized when the damping force at the b node in the end absorbing layer is approximately balanced with the acceleration and the residual nodal load generated by the stiffness.
The Rayleigh damping model in dynamic calculations is a classical approach for structural dynamic response analysis, which is easy to implement and has been implemented in most FEM software because the decoupled governing equations are controlled by only two constants (α, β).This study uses Rayleigh damping to describe the outgoing wave energy a enuation in the absorbing layer, and the Rayleigh damping model in the explicit direct integration method can be expressed as follows [49,50]: where In the above equations, [ ] is the mass damping matrix, [ ] is the stiffness damping matrix, [ ] and [ ] are the mass and stiffness matrices of the whole model, respectively, and α and β are the mass damping coefficient and stiffness damping coefficient, respectively.Although second-order Rayleigh damping is generally more concise in mathematical form, the physical meaning is not clear.Thus, reverse calculations are often adopted by the damping ratio ( ).Therefore, Equation ( 5) can be expressed in the form of [51,52].
where is the mass damping ratio, is the stiffness damping ratio, Rayleigh is the Rayleigh damping ratio, and ω is the circular frequency corresponding to and , and the variation curve of ξ with ω is shown in Figure 2.

Figure 2. Damping ratio versus frequency curve.
As shown in Figure 2, the damping ratios corresponding to different frequency bands vary greatly, which can cause serious distortions in the relative responses of different vibration modes.Generally, the interval ( , ) is chosen for the value of in the seismic calculation of structures, as shown in the dashed frame.Since the Rayleigh damping ratio in this band is smaller than the actual damping ratio ( ), the calculated damping in this band becomes smaller than the actual damping.Therefore, the response of the structure becomes greater than the actual response to ensure a safe design (Equation ( 8)).However, the damping should be maximized so that the outgoing wave is not reflected at the boundary in the wave field simulation.Therefore, the choice of should be set in the frequency band of interest outside the ( , ) interval.
According to the variation curve of with (Figure 2), the value of is linearly proportional to and inversely proportional to .The circular frequency is in the ( , ) interval, and is less than or equal to the actual damping ratio.Therefore, there exists a minimum Rayleigh damping ratio ( ) as follows: where is the damping coefficient, is the critical damping coefficient, and is the circular frequency corresponding to .When is outside this frequency band, its damping ratio is all greater than the actual damping ratio, so the damping is also all greater than the actual damping.In order to maximize the damping effect treatment, the limit is taken for at this time as follows: Since ω is a number tending to 0, Equation ( 7) can be approximated as follows: To ensure the accuracy and convergence of the analysis area, we only apply damping in the absorption layer.Furthermore, to maximize the damping efficiency in the MVDB, a Seismic calculations  n mass damping mode with β equal to 0 is adopted for the absorption layer.Therefore, for the proposed method, α is the key parameter to realize the simulation of the seismic wave field in infinite space.This parameter can be estimated experimentally and depends on the source frequency, the number of absorbing layers, the physical properties of the propagation medium, etc.In this study, an empirical formula for MVDB-based viscous damping (Equation ( 11)) is derived after extensive numerical argumentation, and the damping scheme is easy to implement and more efficient under the same conditions. where is the primary function, (2 ) is the damping amplitude control function, is the excitation frequency, k is the damping amplitude control coefficient, n is the number of absorbing layers, ∆ is the thickness of a single absorbing layer, and L is the range of absorbing layers ( = ∆ × ).
The control parameters of the MVDB determine the absorption effect and the computational cost.The CALM divides the absorbing region into 1~5 layers and continuous layers equally in this study, and the NRB sets the absorbing region into 10 layers.For the reader's ease of reference, the proposed method defines individual absorbing material thicknesses as grid cell dimensions and standardizes them in subsequent comparisons of the calculations.In order to fully consider the "shock" problem when the step size is below the lower limit and ensure the accuracy of the solution [53,54], the grid size is set as follows: where is the minimum wavelength (i.e., transverse wave wavelength).The working mechanism of the MVDB is shown in Figure 3.This study defines the MVDB as a locally a enuated artificial boundary, which gradually a enuates the outgoing sca ered wave energy through viscous damping in the MVDB to ensure that the pulsed wave is not reflected into the analysis area after entering the absorbing layer.In this way, the transmission process of the pulsed wave from the finite domain to the infinite space is simulated.The core idea of the MVDB is to change only the damping value of each absorbing material layer on the boundary to a enuate the outgoing wave energy at the boundary without changing the physical properties of the absorbing layer itself, thus avoiding reflections due to physical differences between the absorbing layer and the analysis area.

Example Verification
In order to verify the absorption effect of the proposed method under the conditions of various seismic wave fields, the numerical models of the NRB, CALM, and MVDB based on the damping absorbing layer method were established to compare the examples.Additionally, the criteria were unified, with the absorbing region size L being taken as 2 ( = 2m).In Equation ( 11 In this study, the model was discretized using four-node plane stress elements, reduced integral (CPS4R) cells based on ABAQUS/explicit dynamics analysis, and hourglass stiffness correction using the artificial damping method (see e.g., [55]), with a torsion control of 0.1 and an hourglass control of 0.45.The cell dimension was ∆ = ∆ = ∆ = 20 cm.The seismic source was a sinusoidal function modulated by a Hanning window as follows: where the source frequency f adopts 100 Hz, the number of Hanning periods n is three, the amplitude scale factor A is 1 × 10 −3 , the effective duration of source excitation t is 0.03 s, the calculation time is 0.7 s, and the waveform is shown in Figure 5.

Endogenous Seismic Wave Field
Figure 6 shows the snapshots of the wave fields for the artificially truncated boundaries and three damping absorbing boundaries.It can be seen that all these methods effectively absorb the wave energy radiated outward.When based on the damping calculation rule and the working mechanism of the NRB, this method may lead to instability in the absorption of low-frequency pulses due to excessive damping increment.Applying excessive damping increases the impedance change largely between the absorbing layers and the initial damping, resulting in outgoing wave reflection in the absorbing layer.In contrast, the proposed method and the CALM are superior in absorbing the outgoing wave energy at the boundary.
For a more accurate analysis of the a enuation effect of our approach on body waves, the vertical displacement amplitudes of sensors A (Figure 7a) and B (Figure 7b) were computed.Their reflection intervals (A: 0.23-0.35s, B: 0.18-0.28s) were enlarged by 100 times, and it can be seen that the developed method and the CALM have an absorption effect on transmi ed longitudinal (P) wave absorption, as shown in Figure 7a.With increasing energetic S waves (Figure 7b), the application of the CALM fails to a enuate sufficient reflected wave energy.This is probably due to the underdamped problem of the continuous damping layer arising from the defective working mechanism of the quality factor (Q) in the CALM.

Exogenous Seismic Wave Field
In this paper, the wave source on the free boundary was categorized as the exogenous seismic wave field incident on the free boundary for research, such as in mineral prospecting and elastic wave damage detection.Figure 8 presents the snapshots of the wave fields of the artificially truncated boundaries and three damping absorbing boundaries, and it can be seen that the P wave energy is effectively a enuated at the boundaries.A small number of reflected surface (RR) waves and RS waves were generated at the side boundaries when using the NRB method.In contrast, the proposed method and the CALM were relatively stable, and the expected objectives were met.
In addition, extracting the vertical displacement amplitudes collected by sensors C and D, as given in Figure 9a,b, can be er illustrate the energy changes on the boundaries.In Figure 9a, the C sensor records the direct wave signals at t ≈ 0.16 s, and a faint energy reflection is detected in the interval of 0.25-0.5 s.The displacement amplitudes in this interval were amplified using a magnification factor of 100, and it can be found that the absorption scheme adopting the NRB does not effectively reduce the transmi ed surface (R) wave reflection.This is a ributed to the excessive impedance of the NRB, reflecting the R wave in the absorbing layer.The wave superposition occurs when the limited wave energy reflected by the absorbing layer and the R wave overlap.The absorbing layer was employed as a new source to radiate energy to the surrounding area, and it was retransmi ed back to the C sensor.The proposed method and the CALM can favorably a enuate the R wave energy on the boundary.Figure 9b shows the vertical displacement amplitudes collected by the D sensor at the bo om boundary.It is indicated that all three absorbing boundaries effectively reduce the P wave.However, the NRB had relatively low effective a enuation of the S wave, and the RS wave was collected at t ≈ 0.3 s and 0.4 s.In addition, the phase of the displacement amplitude curve in this absorption scheme changes at t ≈ 0.19 s and is presumably a ributed to the transient reflection from the absorbing layer.When adopting the MVDB, no prominent reflected wave field is found at the boundary, and only faint reflected signals are observed when the displacement amplitude of the reflection interval (0.3-0.4 s) is enlarged by 100 times.It is indicated that the proposed method can effectively suppress the energy reflection of the R and body waves at the boundary.This phenomenon reflects the propagation law of the pulse wave in the semi-infinite domain and is consistent with the results obtained by the CALM.In numerical simulation validation, the MVDB and CALM perform acceptably in absorbing efficiency and stability, and the former is more convenient to implement and has a lower computation cost.The CALM requires the inverse inference of and based on a predetermined a enuation factor Q 1 coupled with a circular frequency .When applying continuous damping layers, their numbers can be obtained by solving the system of two-element equations cyclically.For the MVDB, after determining L and the thickness of individual absorbing layers ∆l, these parameters can be directly introduced into the equation to solve .Compared to the CALM, the MVDB only requires solving for one unknown quantity.For a detailed explanation, the key workflow of the CALM and MVDB is summarized in Figure 10.

Applicability of Obliquely Incident Wave Fields
In seismic response studies under near-field conditions and damage detection of welded seams, the seismic wave incidence angle is not always perpendicular to the free boundary.For this reason, the effect of obliquely incident waves should be highlighted in the simulation of seismic wave fields under near-field conditions [56], in rock layer fracture damage investigation [57], and in structural non-destructive testing [58].To verify the a enuation effect of the MVDB on obliquely incident pulse waves, a numerical model of homogeneous linear elastic and isotropic soil oblique incidence angles was established.Typically, three source wave fields are involved in the problem of obliquely incident P waves in a single-layer medium at the free boundary.To avoid the interference caused by the wave field superposition and to facilitate the observation of the wave energy a enuation, the source was set at the junction of the absorbing boundaries to eliminate the transmissive wave field in the opposite direction.∆ is taken as the maximum grid size (Equation ( 12)), k is set at 1.08, and L is determined as two times the minimum wavelength ( 2), and the excitation scheme is shown in Figure 11.For example, at an incident angle of 30° and 45°, in this paper, snapshots of the wave fields at 0.063 s, 0.112 s, 0.189 s, 0.28 s, 0.36 s, and 0.45 s are obtained.As shown in Figure 12, the P wave reaches the bo om and side boundaries at t ≈ 0.189 s and traverses the bo om boundary at t ≈ 0.28 s, and no reflection occurs in this duration.The S wave is about to reach the lateral boundary at t ≈ 0.189 s and continues to propagate towards infinity.Similarly, no reflections are generated, which is identical to expectations.

Effect of Damping Control Parameters on Wave Energy A enuation
In Section 2, it has been pointed out that adjusting the damping control parameters of the MVDB approximation formula (Equation ( 11)) may allow intervening in the absorption mechanism, controlling the absorption process, and achieving artificial boundaries in different media and spatial problems.To differentiate the intrinsic a enuation properties of the uniform soil and the a enuation efficiency of the MVDB, the sensors were set up 2 m away from the bo om boundary until the outermost absorbing layer (Sensor t), and the sensor spacing was 0.2 m.For L = 0.5 , , 2 , and 3 , the sensor numbers corresponding to the outermost absorbing layer are Sensor 16, Sensor 21, Sensor 31, and Sensor 41, respectively.The geological conditions are referenced in Section 3, and the model is shown in Figure 13.

The Effect of the Thickness of a Specific Single Absorbing Layer ∆
Regarding the effect of ∆ on energy absorption, the wavelength occupies ten grids (Equation ( 12)), and ∆ is taken to be one and two times the grid size (20 and 40 cm); with k = 1.08,L is taken to be 2 and n = L/∆ .In this paper, the number of absorbing layers, as well as the damping growth curves, were compiled for individual absorbing layer thicknesses of ∆ and 2∆ , respectively, as shown in Figure 14a.It can be seen that the number of absorbing layers n is ten when the absorbing layer thickness is taken to be two times the grid size.By contrast, the quantity of absorbing layers doubles the original one when the thickness equals the grid size.With more absorbing layers, the increasing curve of α becomes smoother, and the damping of the absorbing layers increases in an orderly manner.In this case, the differences in impedance and wave velocity between the absorbing layers have li le variation, and reflection of the pulsed wave between absorbing layers is less likely to occur.In addition, the instantaneous energy peaks of 31 sensors were extracted, and it can be seen that the a enuation curve of the outgoing wave becomes smoother when the thickness of a single absorbing layer is taken as ∆ , and the a enuation ability for pulse waves is enhanced by about 1.4-1.8times that of a two-fold mesh size, as shown in Figure 14b.

The Influence of Absorption Region Thickness L
In exploring the effect of L on wave energy absorption, L values are determined at 0.5 , , e., 1 m, 2 m, 4 m, and 6 m, respectively), k = 1.08, ∆ is the maximum mesh size (∆ = 20 cm), and n = L/∆ .The model was calculated to have a size of 50 m × 50 m, and the thicknesses of the unilateral absorbing region are 1 m and 4 m when L is 0.5 and 2 .In this paper, the number of a single absorbing layer n and the corresponding increase in are responsible for different L. As shown in Figure 15a, the rapid increase in and the relatively steep curves at L = 0.5 indicate that considerable damping is forced to be arranged in a reduced number of absorbing layers, leading to a surge in damping and a huge difference in wave impedance, which triggers interlayer reflections at the junction of the computational region and the absorbing layer and between the layer interiors.When L is equivalent to , the increase in slightly slows down.In the case of L = 2 and 3 , the curve increase is relatively smooth, indicating that the damping required to consume the wave energy matches the number of absorbing layers, and a larger L indicates a more continuous rise in .As shown in Figure 15b, from Sensor 16 to Sensor t, the instantaneous energy amplitudes of four groups of models (L = 0.5 , ) were counted, and the wave energy under damping can be observed in the a enuation amplitude.When L is 0.5 and 1 , the energy decay behavior is weak despite the rapid decrease in the wave energy, and cusps are generated, indicating that interlayer reflections may occur at this position due to excessive damping increase.For this reason, the snapshots of the wave fields at the two L values were specifically extracted, as shown in Figure 16.It can be seen that the outgoing wave is reflected in the absorbing layer when L is 0.5 , and the MVDB works in an orderly manner when L is 2 , verifying the conjecture that the interlayer reflection is generated at the cusp.In addition, the total damping values (Equation ( 5)) of the MVDB under different L and the residual wave energy of the pulse signal at the outermost absorbing layer were statistically investigated, as seen in Figure 17.As shown in Figure 17, the thickness of L determines the total damping values in the absorbing region.A larger L induces a greater total damping value (Figure 17a) and a be er absorption effect (Figure 17b).When L is 0.5 , the wave energy at the bo om boundary is a enuated by about one third; for L = , the wave energy decays to about 1/15 of the original.It reduces to about 1/120 of the original in the case of L = 2 , and when increasing L to 3 , the wave energy decreases to about 1/230 of the original.Combined with Figures 15-17, L is positively correlated with the stability of the MVDB.A larger L indicates a higher accuracy of the MVDB.
Notably, different materials have various physical properties, leading to varying wavelengths.In the present work, the determination of the MVDB action range depends on the wavelength.In this sense, the applicability of the MVDB viscous damping approximation (Equation ( 11)) to other materials should be considered.For a clear interpretation, other materials were introduced for arithmetic examples, and the meaning of L is explained through numerical results.
L refers to the extent of the entire absorbing region, and its core value lies in the design of the number of absorbing material layers when determining the thickness of a single absorbing material.Due to the different physical properties of materials, a research gap remains on the rigorous definition of the number of absorbing layers for studied media and the desirable calculations.Seismic waves at a consistent frequency exhibit varying wavelengths when propagating in different media, and wavelength is only related to the wave velocity and source frequency.For this consideration, it was selected as a medium to design a suitable absorbing layer scheme for different materials according to Equation (11).The given ∆ and the source frequency were used to convert the wavelength into the thickness of the absorbing region L, which is material specific, and then the number of suitable absorbing layers n was derived.In this sense, the significance of L is determining a reasonable number of absorbing layers n.Due to limited space, this paper only cites the numerical case of concrete materials for explanation.
C30 concrete has Young's modulus E = 20 GPa, Poisson's ratio = 0.2, and density ρ = 2400 kg/m 3 , and the effective computational size of the model is 80 cm × 80 cm.The thickness of the unilateral absorbing region is 10 cm for Model 1 and 8 cm for Model 2. The mesh of Model 2 was refined (∆ = 0.4 cm), and the number of absorbing layers is kept consistent with that of Model 1, where k = 1.08, n = L/∆ , calculation time = 1 ms, and time step = 5 × 10 −4 ms.MVDB parameters were obtained from Equations ( 11) and ( 12), as shown in Table 1, and the model is schematically shown in Figure 18. Figure 19 shows the snapshots of the wave fields of Models 1 and 2. It can be seen that the wave energy a enuation of the two models is highly consistent.Despite a small difference in L, the approximate results can still be obtained due to the consistent numbers of absorbing layers, and such a slight difference may be due to the a enuation properties of the remaining thick concrete medium (0.4 ).Based on the above analysis, the conclusion of the 0.5-3 calculation is only feasible for the current mesh accuracy (Equation (12)).Nonetheless, it still reflects that the MVDB is equally applicable to other materials.

Damping Amplitude Control Coefficient k
In order to verify the effect of k on wave energy absorption, the MVDB was numerically calculated for k = 1.02, 1.07, 1.09, 1.10, and 1.12, L is taken as 2 and n = L/∆ , and ∆ is determined as the singular mesh size (20 cm).As shown in Figure 20a,b, k is positively correlated with the total damping of the absorbing layers and the peak value of damping.With the increase in k, the slope of the mass damping coefficient curve becomes larger, and the total amount of damping in the absorption layers increases exponentially.
In addition, the correlation between k and wave energy decay can be intuitively reflected by monitoring the residual wave energy of the outermost absorbing layer (Sensor 31) of the four models.As shown in Figure 21, when k increases from 1.02 to 1.07, the residual wave energy of the outermost absorbing layer is reduced by about 19.8%.When k = 1.09, 1.10, and 1.12, the residual wave energy is reduced by about 31%, 34%, and 36%, respectively.

Discussion
This is a semi-infinite spatial dynamic characterization study based on a commercial finite element framework to provide convenient and efficient absorbing boundaries for seismic wave field simulations.Here, the a enuation efficiency of the absorbing material is determined using the mass damping coefficient to set the absorption boundaries, and this method is also applicable to other finite element software.Although the MVDB improves work efficiency, it still cannot overcome the limitation of depending on the thickness of the absorption region like the artificial boundary of other absorption layer methods.
In addition, the case analysis demonstrates that the k value correlates to the seismic source frequency.As an exponential term of the equation, the damping amplitude control coefficient k significantly affects the growth rate of .As k increases, both the maximum and total values of damping increase significantly.Nevertheless, when the seismic source frequency is high or L is thin, an excessive k value may cause the impedance of the absorbing layer to change more rapidly, resulting in the reflection of the pulse wave between the absorbing layers.This work seeks to enhance the proposed method's applicability to various conditions.For this objective, seismic wave fields at multiple seismic source frequencies were simulated, followed by the calculation of the damping growth rates and the residual wave energy in the outermost absorption layer.By comparing the minimum values of residual wave energy, we obtained the approximate range of k values, which vary with the frequencies.The approximate range and frequency of k values are summarized in Table 2, which is recommended for the selection of k.

Conclusions
This study introduces a new boundary-reflected wave processing scheme (MVDB) for structural damage detection and geological exploration.Compared with the traditional absorbing boundary, the MVDB does not need to define the wave equation of the absorbing region or derive the expression of the unilateral outgoing wave, and a be er absorbing effect can be achieved only by se ing up a certain thickness of the damping incremental absorbing layer at the periphery of the calculation region.Compared with the CALM, the MVDB can set boundary damping without increasing the grid cells of the absorbing region and only needs to solve the mass damping coefficients.Therefore, it is more advantageous for treating reflected waves on artificially truncated boundaries.With the MVDB, the geometrically centered explosion source and the input spherical wave on the free boundary are simulated, respectively, both of which verify the effectiveness of the proposed method.The main conclusions are as follows: (1) The MVDB implementation process is simple, and the physical meaning is clear.Independent of the space-time discrete step, its calculation process is stable.According to the numerical results of a two-dimensional uniform soil, the MVDB effectively attenuated the wave energy of the body wave and R wave in the endogenous and exogenous seismic wave field under the vertical incidence of the P wave, showing high accuracy and stability.The MVDB is also applicable to the oblique incidence problem on the free boundary.(2) The thickness of the absorbing region L determines the total damping value of the absorbing region and the smoothness of the damping growth, making L the key factor affecting the wave energy a enuation, and the energy a enuates to about 1/4, 1/140, and 1/230 of the original when L is , 2 , and 3 .Based on the numerical results, the 2~3 thickness of the absorption region can meet the general engineering requirements for accuracy.
(3) When the thickness of the absorption region L is constant, and the thickness of a single absorption layer is set as a discrete grid size ∆ , the a enuation trend of the boundary wave field energy is relatively stable without interlayer reflection, thus achieving a be er absorption effect.(4) The damping amplitude control coefficient k in the MVDB is key in determining the initial damping value and smoothing the damping growth.When the thickness of the absorbing region is thinner, an excessively large k triggers interlayer reflections due to the uneven damping growth in the absorbing layer.A be er absorption effect is ensured for incident waves of different frequencies when k is between 1.02 and 1.12.
As a new boundary scheme, the MVDB provides a promising way to study seismic wave fields in semi-infinite space.In future work, specific research tasks will be carried out based on the MVDB, such as advanced prediction of tunneling and non-destructive testing of anchor rods.Additionally, the efficiency of the MVDB in 3D mode also needs to be evaluated.

Figure 1 .
Figure 1.Schematic of the seismic wave field.(a) Unbounded input waves.(b) Free boundary input waves.

Figure 4 .
Figure 4.The uniform soil calculation model.(a) Endogenous seismic wave field.(b) Exogenous seismic wave field.

Figure 6 .Figure 7 .
Figure 6.Wave field snapshots in the case of the endogenous seismic wave field.

Figure 8 .Figure 9 .
Figure 8. Snapshots of the wave fields for the exogenous seismic wave field.
frequency based on the relationship between the main frequency f of the seismic source and the period Reverse calculations and based on the circular frequency and the damping ratio Implement Continuous CALM

Figure 11 .
Figure 11.Numerical model of the obliquely incident wave field.

Figure 12 .
Figure 12.Snapshots of the wave fields for an obliquely incident uniform soil.

Figure 13 .
Figure 13.Numerical model of uniform soil.

Figure 19 .
Figure 19.Snapshots of the wave fields.

Figure 20 .
Figure 20.The relationship between the damping amplitude control coefficient k and wave energy absorption.(a) growth amplitude.(b) Total damping and peak damping.

Figure 21 .
Figure 21.Residual wave energy amplitude detected by Sensor 31.The thickness of the absorbing region is 2 .

Table 2 .
Reference values of damping amplitude control coefficient k in the MVDB absorbing layer.