The Caughey absorbing layer method – implementation and validation in Ansys software

This paper addresses the Caughey Absorbing Layer Method (CALM) performance in the one-dimensional problem and its implementation in commercial software, with possibility of direct extension to two-dimensions. The adequacy and numerical efficiency is evaluated using three different error measures and five different variations of the damping profile. Other parameters that are subjected to evaluation are the length of the absorbing layer in relation to the wavelength to absorb, the value of the loss factor at the end of the absorbing layer, and the ratio of the load to layer frequency. The problem is firstly analysed theoretically, resulting in estimates for the wave reflection due to transition and truncation of the model. In order to confirm that no spurious waves will be present in the finite element solution, the numerical implementation is validated by comparison with the analytical solution. The analysis of the error measures on the numerical results obtained for various combinations of the model’s parameters lead to the conclusion that CALM is effective at mitigating waves reflected from the boundaries. The optimum loss factor as a function of the ratio of the length of the absorbing layer to the wavelength to absorb is determined through parametric analyses. Although the optimal damping is frequency dependent, it was shown that the CALM’s effectiveness can be extended to a wider range of frequencies by increasing the smoothness of the damping profile.

The Caughey absorbing layer method -implementation and validation in Ansys software

INTRODUCTION
One of the most significant difficulties in the numerical study of elastic wave propagation in solids, particularly when using the finite element (FE) method (see Hughes (1987)), is the difficulty to simulate semi-infinite or unbounded domains.
For example, in the analysis of soil vibrations, the soil region is neither confined to a closed space nor isolated from the surroundings.This means that modelling the area of interest without special considerations will result in reflections of the elastic waves at the boundaries of the model.These reflections will become superimposed with the actual solution, making it inaccurate.As the waves are locked inside the finite domain, the energy is not released as it happens in reality.
In an analytical analysis, it is common to admit the soil as a semi-infinite medium.This approach was employed by Boussinesq, who studied the stresses on the soil due to a static load (see Karol (1960)).Since the domain of most numerical methods must be itself finite, various truncation techniques have been proposed over the last decades, such as: (i) local absorbing boundary conditions (Lysmer and Kuhlemeyer, 1969;Lindman, 1975;Engquist andMajda, 1977, 1979); (ii) the boundary element method (Banerjee and Butterfield, 1981); (iii) the infinite element method (Bettess, 1977); (iv) absorbing layers, including Perfectly Matched Layers (PML, Bérenger (1994)) and the Caughey Absorbing Layer Method (CALM, Semblat et al. (2011)).
The local absorbing boundary conditions are among the simpler methods.They are implemented as dampers at the boundary, whose properties are derived from changing the wave equation to admit only outgoing waves.The Lysmer-Kuhlemeyer boundary in particular uses damping coefficients that are proportional to the mass density of the medium and the speed of the waves to absorb.The rate of absorption depends on the angle of incidence of the wave, and is usually tuned to perfectly absorb only at normal incidence.Although effective for simple problems, absorbing boundaries may lead to instabilities at discontinuities (such as layers with different mechanical properties).Furthermore, in order to avoid rigid body motion, these boundary conditions must be complemented with additional elements, like springs, whose stiffness is difficult to estimate.
The boundary element method changes the nature of the numerical problem, from a volume discretisation to a boundary discretisation.This requires a fundamental solution that satisfies the appropriate conditions at infinity.For dynamical problems, the Sommerfeld's radiation condition is used, which states that "the energy … radiated from the sources must scatter to infinity; no energy may be radiated from infinity into ... the field" (Sommerfeld, 1949).Although the BEM is very robust, the computational cost is much higher than traditional FE -for many problems where the surface to volume ratio is high, the boundary method may be less efficient than volumediscretisation methods (see Katsikadelis (2002)).
The infinite element method is closer to the traditional FE approach.It consists in modelling the interior domain with conventional finite elements and, at the model's boundaries, introducing elements with special shape functions.These special shape functions grow without bound as the coordinate approaches infinity, therefore simulating an infinite element.Unfortunately, only some commercial FE programs include this formulation.
Absorbing layer methods have been widely used since the introduction of the PML by Bérenger (1994), based on previous absorbing layer formulations (e.g.Holland and Williams (1983)).The absorbing layer method applies a layer of material with some damping capability at the boundaries of the medium of interest.Waves behave normally inside the original medium, but decay as they travel inside the absorbing layer, attenuating or preventing reflections at the boundaries of the model.However, some reflection is expected to occur at the interface between the normal medium and the absorbing layers, due to the material discontinuity.Structures 12 (2015) 1540-1564 The PML in particular can be implemented with a complex coordinate stretching (see Chew and Weedon (1994)).The analytical formulation does not introduce reflections at the interface between the two materials (hence perfectly matched), but this property is partially lost after discretisation.

Latin American Journal of Solids and
The main drawback of the PML is that its implementation is not straightforward, particularly in the time domain -it requires a split-field formulation or convolution operations.This makes it difficult to use in FE commercial software.
The CALM, proposed by Semblat et al. (2011), is not perfectly matched, but much simpler to implement than the PML.The absorbing layer has the mechanical properties of the medium of interest, but exhibits Caughey (or Rayleigh) damping tuned to ensure that the rate of absorption for the desired frequency is above an arbitrary value.It has the advantage of being intrinsically multi-directional, unlike local absorbing boundaries and the PML.Since it only requires manipulation of the FE damping matrix, it can be easily implemented in FE commercial software.The authors have found promising results in the implementation of the CALM in the commercial FE software Ansys (see Kohnke (1999)) using an implicit dynamics formulation for one-and twodimensional plane-stress problems (presented in Rodrigues and Dimitrovová (2014)).
In the present paper, the authors implement the Rayleigh formulation of the CALM in the commercial FE software Ansys, using implicit time integration, for a one-dimensional plane-strain problem.The model is subjected to a displacement in one of its extremities, which follows a wavelet function with a well-defined dominant frequency, and therefore the wavelength inside the medium is also well-defined.The CALM's performance is analysed using three different error measures: the maximum amplitude of the reflected waves, the maximum L 2 -norm of the reflected waves, and the time integration of the L 2 -norm of the reflected waves.
It is observed that the optimum properties for the absorbing layer require a compromise between reducing round-trip and transition reflections.In this context round-trip reflection refers to the waves reflected at the end of the absorbing layer and transition reflection represent the waves that reflect at the interface between the medium of interest and the absorbing layer.
Round-trip reflection can be reduced simply by increasing the damping in the layer, but this increases transition reflection.To mitigate this problem, the damping is introduced smoothly along the layer.Different variations of the damping profile are tested (constant, linear, quadratic, cubic and exponential) and their efficiency is compared.It is concluded that a variable damping profile leads to less reflections than constant damping across the absorbing layer.The linear variation is the most efficient for short layers, while the quadratic is best for long layers.
For each damping profile, a parametric optimization is performed to obtain the optimum loss factor at the end of the layer as a function of the dimensionless ratio of the length of the absorbing layer to the wavelength to absorb.From this process, empirical estimates for the optimum loss factor are proposed and shown to fit the numerical results almost perfectly for most cases.These formulas can be used to tune the CALM for different problems, as long as the wavelength to absorb can be estimated.
Lastly, a sensitivity analysis is performed to determine how a difference between the frequency of the load and the frequency that the CALM is tuned to absorb affects the results.Unlike Semblat et al. (2011) suggested, the efficiency of the CALM does not increase with this difference, due to a higher transition reflection coefficient.However, it is seen that the higher order damping profiles are Latin American Journal of Solids and Structures 12 (2015) 1540-1564 less sensitive to this discrepancy, and may be a better choice for problems were a wide range of frequencies is present, or were the predominant frequency is not known.
In brief, the new contributions of this paper are: (i) previously unpublished analytical estimates for the effectiveness of the CALM, which support the empirical observation that a compromise between round-trip and transition reflections must be found; (ii) implementation of the CALM in Ansys implicit dynamics module (in particular the variable Rayleigh damping); (iii) empirical estimates for the optimum loss factor as a function of the ratio of the layer's length to the wavelength for various damping profiles, which can be used to tune the CALM for different problems; (iv) demonstration of how the CALM's effectiveness can be extended over a wider range of frequencies by employing a smoother damping profile.
The paper is organized in the following way: in section 2, the theoretical principle behind the CALM is presented, theoretical values for the reflection coefficients are proposed, and the solution by modal superposition is discussed.Section 3 details the implementation of the CALM in the FE commercial software Ansys, which is then validated by comparison with the solution by modal superposition; the error measures are presented and then used to perform a parametrical optimization of the CALM parameters; lastly, a sensitivity analysis is performed.Section 4 summarizes the results and discusses future work.

THE CAUGHEY ABSORBING LAYER METHOD
The Caughey Absorbing Layer Method was proposed by Semblat et al. (2011) as a simple and reliable alternative to the Perfectly Matched Layer (PML).
Like other absorbing layer methods, the CALM consists in defining an absorbing layer at the boundaries of the elastic medium under consideration.This absorbing layer is modelled with the same elastic properties as the interior medium, but damping is added to attenuate all waves that leave the interior domain.Unlike the classical formulation of the PML, the CALM presents multidirectional attenuating properties.
To ensure that the absorbing layer exhibits the desired behaviour, the damping properties must be carefully chosen to apply an adequate level of attenuation to the relevant frequency.

The Rayleigh damping formulation
One of the simplest ways to define damping in FE analysis is to consider the Rayleigh formulation (see Liu and Gorman (1995)), where the damping matrix (C) is assumed to be a linear combination of the stiffness (K) and mass (M) matrixes:  and  are known as the Rayleigh coefficients.This approach is convenient because FE methods already require the assembly of the mass and stiffness matrices.
It is well known (see, for example, Clough and Penzien (1993)) that the loss factor (  , the ratio of energy dissipated to the energy stored in the system for every oscillation) is approximately dou-Latin American Journal of Solids and Structures 12 (2015) 1540-1564 ble the damping ratio (  , the ratio of the damping coefficient to the critical damping coefficient), which, for Rayleigh damping, relates to the frequency of excitation (  ) and to the Rayleigh coefficients: Following Semblat et al. (2011), the loss factor will reach minimum when the frequency of excitation is Since this is the minimum absorption possible, if the Rayleigh coefficients is defined to get a desired loss factor min  for a certain frequency r  , all excitations will be damped at least as much as the value defined.As the frequency of excitation moves away from r  , the more the damping will be felt (see Figure 1).The Rayleigh coefficients that correspond to min  are obtained from Equations (2) and (3): Conversely, given a particular frequency i  , the loss factor associated with it can be obtained from Equations ( 2) and (4): Latin American Journal of Solids and Structures 12 (2015) 1540-1564 Semblat's approach is to define the frequency of minimum attenuation, r  , as being equal to the predominant frequency of the waves to absorb.In theory, this method ensures that all waves will be attenuated at least by a factor proportional to min  .The Rayleigh damping coefficients to adopt are then calculated using Equation (4) for the predominant frequency and the desired minimum attenuation.However, this does not take into account that when the waves travel from the medium into the absorbing layer, the sudden introduction of damping causes reflections, the amplitude of which increases with the intensity of damping.This phenomenon will be discussed over the next sections.

Absorbing layer properties
The CALM formulation presented above allows the application of an arbitrary minimum factor of attenuation to all vibrations in the selected medium.The minimum value of attenuation to apply ( min  ) and the frequency to which apply it ( r  ) are choices that depend on the particular problem being studied and how that problem is modelled.
The case study presented in this paper is that of a semi-infinite rod: a one-dimensional medium with one free end and that extends to infinity.An excitation is applied at the free end, producing waves that travel along the semi-infinite medium.Since this excitation has a well-defined dominant frequency, f  , the wavelength  can be obtained from the frequency and elastic properties of the medium.As usual, small displacement theory is adopted.
The infinite nature of this problem can be avoided by modelling a finite rod with an absorbing layer ending in a Dirichlet boundary, as shown in Figure 2. The rod is assumed to be homogeneous, with Young modulus E and mass density  .The length of the medium of interest is h , and the length of the absorbing layer ( L ) is defined as a multiple (a ) of the wavelength (  ) produced by the excitation.This wavelength is a function of the pressure wave speed ( P c ) in the medium and the frequency of the excitation ( f  ): where the pressure wave speed is defined as: It is not necessary for a to be an integer, it can be any positive real number, so the length of the absorbing layer, L , can also be any positive real number.and Structures 12 (2015) 1540-1564 The loss factor of the absorbing layer,   y  , is a function of the distance from the interface between both media, y .This loss factor can either be constant (   monotonically along the layer, from no damping to maximum damping (   . The Rayleigh damping is then applied according to Equation (4).
It has been noted (Oskooi, 2008;Oskooi et al., 2008;Oskooi and Johnson, 2011) that absorbing layers, like PMLs, perform better when a variable absorbing profile is defined.The smoothness of the variation was also found to influence their effectiveness.This is a result of the presence of two sources of reflection in the model: round-trip reflections and transition reflections.
Round-trip reflection occurs because the absorbing layer never absorbs completely the incoming waves.As the waves travel along the layer, their amplitude decreases exponentially according to the loss factor, but they still reflect of the boundary condition and travel back into the medium of interest -hence the "round-trip" designation.
Transition reflection occurs wherever there is a change of the properties of the medium of propagation.The steeper the variation, the higher the reflection will be.This is the case at the interface between the medium of interest and the absorbing layer when the loss factor is constant.This problem is mitigated when the loss factor increases gradually along the layer.
Although Semblat et al. (2011) tested the effectiveness of the CALM numerically and found that maximizing its effectiveness required finding a compromise between minimization of both types of reflections, no theoretical justification for this observation was provided.

Analytical solution by modal superposition
The problem of wave propagation in a single dimension due to an applied force at 0 x  can be expressed using the following differential equation: where   , u x t is the displacement of the medium at the position x at time t ,   b x is the viscous damping,  is the mass per length of the medium,   P t is the magnitude of the force over time and   x  is the Dirac delta function.This problem can be solved by modal superposition, expressing the solution as the sum of an infinite number of functions: where   i v x is known as the i -th undamped vibration mode and   i q t is the corresponding generalized coordinate.For the boundary conditions shown in Figure 2 (free end at where i  is the i -th undamped natural frequency.The generalized coordinates are the solutions to the system of equations If the damping is constant over the medium (   b x b  ), then all equations are independent: In the case of the CALM, the damping ratio i  for each frequency must be obtained from Equation (5).
Equation ( 12) allows one to obtain exact generalized coordinates for any desired mode.The system of equations ( 11), on the other hand, has to be solved numerically for a finite number of modes.In both cases an approximate solution is obtained by modal superposition.In the first case, if further precision is needed, one can solve Eq. ( 12) for higher modes and update the solution.In the second case, if higher modes are needed, the entire system of equations must be solved again.
In the case of the CALM the damping is not constant over the length of the medium, and therefore the system of equations ( 11) must be solved numerically.

Theoretical reflection coefficients
In the study of the phenomenon of reflection it is useful to introduce the concept of the reflection coefficient: the ratio between the amplitude of the reflected waves and the amplitude of the incoming waves (see Chapman (2004)).Although analytical expressions for the reflection coefficients associated with the PML have been published by Oskooi (2008), no such formulas exist yet for CALM.
The round-trip reflection coefficient can be approximated by considering how much damping the waves are subjected to while traveling through the absorbing layer.By analogy with a single degree of freedom oscillator, it is assumed that the amplitude of the waves decays exponentially due to damping.This assumption was supported by the analytical solution of a rod with constant damping subjected to a sinusoidal force on the free-end (solving Equation ( 12) for     ).The logarithmic decrement (  ) expresses the decrease in amplitude per cycle of oscillation, and is a function of the damping ratio  , which is itself a function of the loss factor (Equation (2)): If the characteristics of the wave are known, it is possible to express the variation of the logarithmic decrement over the length of the absorbing layer as a function of   y  , since a cycle of oscillation occurs over the distance of a wavelength: where is the modified wavelength due to the presence of damping.This value can be obtained by considering the theoretical expression for the damped frequency: Latin American Journal of Solids and Structures 12 (2015) 1540-1564 which leads to the damped wavelength: Since the waves must travel twice the length of the layer before returning to the medium of interest, the total decrease in amplitude after the round trip can be obtained by integrating (16) over  and multiplying by two.If one expresses the variable loss factor as the product of  by a dimensionless damping profile   the loss factor becomes: Equation ( 18) shows that the loss factor is proportional to the area of the damping profile.
A typical damping profile for absorbing layer techniques is a simple monomial equation of degree , for which the loss factor is: Since the amplitude of the waves decays exponentially, the round trip reflection coefficient is It is evident from Equation (20) that increasing both the loss factor (  ) and the relative length of the layer ( a ) reduces the amplitude of the reflected waves (by increasing the integral of the damping profile).Increasing the degree d of the damping profile increases the amplitude of the reflected waves (it reduces the integral of the damping profile).However, increasing d leads to a smoother damping profile, which reduces the transition reflection (see Oskooi (2008)).
Transition reflection can also be estimated theoretically.For a P-wave at normal incidence, it is well known (see Lowrie (2007)) that the transition reflection coefficient is: For the calculation of the damped wavelength, the wave velocity was assumed to be constant.However, for purposes of calculating the transition reflection, the wavelength can be assumed to be fixed and the wave velocity expressed as a function of the loss factor instead: For low damping, T R increases slowly with  , but as  gets closer to 2 (critical damping), T R grows very rapidly to 1 (total reflection of the incident waves).For 2   , oscillatory systems present super-critical damping, where no oscillation occurs.Therefore, Equation ( 21) is not valid for 2   , and the validity of Eq. ( 19) cannot be guaranteed as well.
Unlike the round-trip reflection, the transition reflection coefficient for a damping profile with 0 d  is quite difficult to derive.Although some work has been done for linear variation (see Wolf (1937) and Liner and Bodmann ( 2010)), it is not directly applicable to the problem at hand.However, Oskooi (2008) has shown numerically that T R decreases as both d and L increase.Particularly, for one-dimensional waves, Oskooi observed the following relation According to Equations ( 20), ( 22) and ( 23 2011): higher loss factor increases reflection at the interface, longer absorbing layers reduce round-trip reflection and a linear damping profile performs better than a constant one.

Limitations
In theory this methodology should work with any kind of dynamic analysis, both for time and frequency domain, since the use of Rayleigh damping ensures that the modes of vibration are uncoupled.
Preliminary tests have shown that the explicit central difference time integration method used by the LS-Dyna software (see Hallquist (2006)) requires an unreasonable computational cost to solve this problem -the time-step is five to six orders of magnitude lower than what is recommended for the same problem without the desired damping.
It is known that the stable time-step for explicit central difference time integration in the presence of damping is: where max  is the frequency of the highest mode and max  is the corresponding damping ratio.As the damping ratio increases, the stable time-step decreases.Since Rayleigh damping is felt more strongly in the higher modes, it is inevitable that it will lead to very small time-steps and therefore to high computation times.All models are therefore implemented in Ansys' implicit dynamics module.

NUMERICAL IMPLEMENTATION
To test the effectiveness of the CALM in the FEM, the one-dimensional problem presented in Figure 2 is implemented in Ansys (see Kohnke (1999)) as a mesh of quadrilateral plane-strain elements.This allows generalization to two-dimensional models for future work.A plane stress analysis was presented in previous work by the authors (Rodrigues and Dimitrovová, 2014).
Figure 3 shows the FE model: the rod is discretised over its length as a regular mesh and the length of the model is equal to four times the wavelength to absorb, plus the length of the absorbing layer, which is also a multiple of the wavelength, as defined before.The Young modulus ( E ) is 200 MPa and the mass density (  ) is 2000 kg/m 3 .Since the objective is to model a uniaxial problem, the Poisson ratio is defined as zero.Both in plane strain and plane stress analysis, if 0   the pressure wave speed is still given by Equation ( 7).This leads to P c  316 m/s.
The excitation on the free end of the rod is an impulse displacement, with its time history equal to a second-order Ricker Wavelet (see Hosken (1988)): where 0 U is the maximum amplitude of the wave, t is the time coordinate, p t is the fundamental period of the wavelet, s t is the time shift and t is the time coordinate.Assuming a fundamental frequency f   500 rad/s, the period of the wavelet is 2 0.012566 The time shift was assumed to be equal to the fundamental period, and the maximum amplitude equal to 10 -3 m.
Knowing the fundamental frequency of the excitation, the wavelength of the pressure waves can be estimated as 2 3.974 m The size of the elements was chosen to be exactly 24  , to reduce numerical wave dispersion.This also means that the results obtained are independent of the wave-speed of the material and the frequency of the applied displacement.The time instant when the wave-front of the Ricker wavelet reaches the interface between the medium of interest and the absorbing layer ( w t ) will be necessary to analyse the efficiency of the CALM: The duration of the analysis is defined as the time needed for the waves to travel from the origin to the end of the absorbing boundary and back to the origin again: The duration of each time-step was assumed to be proportional to the time it takes for a wave point to traverse a single finite element:   24 24 This time step is usually considered to be the critical time-step size for explicit integration methods in undamped problems.Although the implicit method converges independently of the time-step size, Equation (30) leads to a good resolution of the solution, which will be important when analysing the results.Ansys' implicit integration method automatically determines if multiple sub-steps are needed to ensure accuracy of the results.This definition of the time-step is also useful for comparison, since it depends only on the frequency of the load.
The discretisation of the absorbing layer is straight-forward when   y    , but for a nonconstant damping profile, there are two options: assume a constant loss factor for each element, but different from element to element (using the value of   y  at the middle point of the element), or include the variation of the loss factor in the definition of the damping matrices.Both options were initially implemented, but for the level of mesh refinement considered ( 24  ), the results are practically indistinguishable.Since the latter option is more complex to implement and brings no benefits, all further tests were performed using constant damping inside each element.
It is worth noting that Ansys' pre-packaged Rayleigh damping implementation applies the same a coefficient to the whole model.To define the desired damping, six discrete damping elements have to be added to each quadrilateral element of the absorbing layer, connecting each pair of nodes in the element, as illustrated in Figure 4.The  -damping is included on the element formulation, so the damping coefficients to be added do not depend on the type of analysis (since the mass matrix is the same for plane-stress and plane-strain analysis).x y e e  Latin American Journal of Solids and Structures 12 (2015) 1540-1564

Model validation
Before testing the effectiveness of the CALM implementation on FE, it is important to validate the model by comparing the results with an analytical solution, to ensure no spurious reflections occur.
To that effect, a model with relative length of the absorbing layer of 4 a  is considered.The damping is constant inside the layer ( A force was applied on the free end, with variation in time defined by the Ricker wavelet (Equation ( 25)) and maximum intensity of 1 N.All other parameters are as above.
To solve the problem analytically, the damping profile over the whole medium must be defined.In this case, since , the damping profile is: where   H x is the Heaviside step function (0 for 11) can be written in matrix form for a finite number of modes: where i  is the damping ratio for the frequency i  , which for the CALM is obtained from ( 5), and 1 is the unit vector.Each element in the matrices C , W and V is computed analytically according to: This system of equations is then solved using Matlab's (see Moler (2008)) explicit integration tool for the first 300 modes.The solution is then assembled by modal superposition.
The results are compared with the solution obtained for the FE model in Ansys.Figure 5 shows both solutions at instants (just before and shortly after the pulse enters the absorbing layer, respectively).Figure 6 shows the solution in points over the full analysis time (inside the medium of interest and the absorbing layer, respectively).
It is clear from Figures 5 and 6 that the traveling pulse is well represented by the FE model, particularly before it enters the absorbing layer, since until that instant the problem is linear and elastic, and the mesh is sufficiently discretized.The major difference is that the amplitude of the transition reflection is considerably lower in the FE model (less than half of what is observed in the analytical solution).However, since these reflections are an undesired effect of the presence of the absorbing layer, this difference does not compromise the accuracy of the solution.
It is also visible that the wavelength increases when the pulse travels inside the absorbing layer, as postulated in the derivation of the round-trip reflection coefficient (Equation ( 16)).

Error measures for the evaluation of CALM's effectiveness
Since the results of the numerical model were found satisfactory, the next step is to evaluate the effectiveness of the CALM, which can be done by measuring the error introduced in the solution due to both round-trip and transition reflection.To assess the effectiveness of each implementation, objective measures of this error will be defined.
Latin American Journal of Solids and Structures 12 (2015) 1540-1564 Three different error measures are considered: the maximum amplitude of the reflected waves, the maximum 2 L -norm of the reflected waves and the time integration of the 2 L -norm of the reflected waves.To define any of these error measures, it is necessary to separate the vibrations that are due to the presence of the absorbing layer from the ones that would occur if the medium was truly semi-infinite.
Consider the numerical solution inside the medium of interest for the problem with the absorbing layer for a given time instant i t to be i u .The vector i u contains the displacement at the timestep i for each node inside the medium of interest, but not inside the absorbing layer.
The solution for a discretised semi-infinite medium can be approximated if one considers a sufficiently long model with no absorbing layer.For the case-study, a model with a total length of 16 is enough to ensure that the waves do not reflect of the Dirichlet boundary for the duration of the analysis ( f t ).If the solution for the long model is i  u , then the displacements due to reflection for the CALM model can be obtained by simple subtraction: where m is the time step corresponding to f t .Assuming 1 0 The first error measure is the maximum displacement of the reflected waves, and is analogous to the reflection coefficient discussed in Section 2.2.It can be defined as The maximum L 2 -norm of the reflected waves is normalized by dividing it by the maximum L 2norm of the waves in the long model: where the L 2 -norm of a vector is defined as Lastly, the time integration of the L 2 -norm of the reflected waves is similar to Equation (37), but replaces the maximum of the L 2 -norm by an approximation to the integral of the L 2 -norm over the time interval from the instant when the wave-front of the Ricker wavelet reaches the interface ( w t ) to the end of the simulation ( f t ), divided by the time range to preserve dimensional consistency.
The approximation is simply the sum of the L 2 -norm for all the time-steps in the range, divided by the number of time steps: Latin American Journal of Solids and Structures 12 (2015) 1540-1564 where w is the time-step corresponding to the instant w t .

Parametric optimization of the loss factor
A parametric test is then performed to find the optimum loss factor as a function of a , which varies from a  0.25 to a  4 with increment a   0.25.The loss factor at the end of the absorbing layer initially varied from The damping profile also varied: the monomial profile (Equation ( 19)) was used with d  {0,1,2,3} (constant loss factor and linear, quadratic and cubic variation).A fifth profile, proposed by Oskooi et al. (2008), was tested: This profile has the property of being C  and therefore is smoother than any monomial function, which should reduce transition reflections.The theoretical round-trip reflection coefficient is which is lower than that of a monomial profile with d  2.
The three error measures are computed for each combination of loss factor, damping profile and absorbing layer length.The optimum loss factor (   ) is the one that minimizes the selected error measure for a particular combination of the absorbing layer length and the damping profile.
It should be noted that, for very short absorbing layers, the optimum loss factor may not be in the range 0 <  < 2. In those cases, the range must be extended to higher values.
After determining   using , which leads to a better approximation of   .The results for 0 d  (constant loss factor) are summarized in Table 1.
The optimum loss factors obtained for each error measure are close to each other, and in particular, the values obtained for max u and 2 t L are practically the same, save for a few exceptions.
Table 2 presents the same information for a linear profile (d  1). 2 that the error for the linear profile is much lower than that of the constant loss factor for any length of the absorbing layer with the exception of a  0.25.It is also evident that the optimum loss factor is considerably higher.This agrees with Equation (20), that shows that higher values of d require higher  to keep round-trip reflections low.• The constant loss factor approach improves very poorly with the increase of a .It is only the best choice for a very short layer (a  0.25), and even in that case the amplitude of the reflected waves is not negligible ( max u  0.14); • For 0.5 1 a   , the linear damping profile presents the best results;

It is clear from Table
• For a > 1, the quadratic profile outperformed all other options; • The cubic profile behaves poorly for short absorbing layers, and even for long layers, it rarely outperforms the linear profile; • The exponential profile is very close to the quadratic one for short layers (a < 1), but does not keep up as a increases.It is possible that the cubic damping profile performs better for longer absorbing layers, as was observed by Oskooi et al. (2008), where, for some problems, the advantage of higher order profiles only becomes apparent for very long layers (a > 10).However, very long absorbing layers usually defeat the purpose of their implementation: reducing the size of the model so its solution becomes computationally expedite (or even feasible).

Optimum loss factor as a function of the layer's length
Having the optimum values of the loss factor as a function of the layer's length in proportion to the wavelength allows generalization of the obtained results to other material properties.Since the properties of the material are already taken into account in the calculation of the wave speed and wavelength, it seems reasonable to expect the optimum loss factor to be independent of the material parameters.
For that purpose, empirical formulas may be devised from statistical analysis of the results.Various tests show that the best fit is a power function, with a as the independent variable.The power function is also a good theoretical approach: as the length of the layer tends to zero, the loss factor needed to absorb the incident elastic waves grow to infinity; as the length grows to infinity, the loss factor diminishes until no damping is needed at all.Table 3 presents the approximation for all the variations and the three error measures, along with the coefficient of determination ( 2 R ) of the approximation.
Table 3 clearly shows that the higher order damping profiles require a higher loss factor at the end of the layer.The exponential variation (Equation ( 40)) is between d  1 and d  2 in that regard.The consistency between the optimum loss factor obtained using the different error measures is evident, which gives confidence in the results.Lastly, it is clear that there is a better approximation for the loss factor obtained using 2 t L than the other two measures, particularly max u .This suggests that 2 t L may be a more reliable measure, since it is a function of the solution over the relevant time-range, instead of simply taking the maximum value.

Latin
The formulas presented in Table 3 were tested for different combinations of material parameters ( E ,  and  , the Poisson ratio), with no noticeable effect in any of the error measures, as long as all dimensional relations are respected.

Sensitivity analysis
Since it is not always possible to clearly define the frequency content of the loads, or it may happen that a wide range of frequencies are relevant, it is important to test if misjudging the prevailing frequency does not lead to a drastic drop in the efficiency of the CALM.
To that effect, a model where the frequency of the absorbing layer ( CALM


) was changed to take different values from the Ricker wavelet ( R   500 rad/s), from CALM   5 rad/s to 50 000 rad/s.The absorbing layer length is equal to the wavelength, and the five different damping profiles were tested.The loss factor used was the numerical optimum obtained in section 3.2 for a  1 and the various values of d .Figures 10 to 14 show the error measures as functions of the ratio of the load frequency to the absorbing layer frequency.
From the analysis of the results, one can see that the damping parameters defined in Equation (4) lead to maximum absorption of the selected frequency.As CALM  moves away from R  , all three error measures increase.This is in contrast with the behaviour intended by Semblat et al. (2011) when they proposed the damping parameters in Equation ( 4).According to Equation (2), the waves with a frequency different from the one the CALM are subjected to a greater loss factor η, which decreases the round-trip reflection but increases the transition reflection.Essentially, the compromise between round-trip and transition reflection found in section 3.3 is no longer valid for the new frequency.Finding a new compromise is only possible if one knows the ratio , in which case it would be possible to simply tune the CALM to the load frequency.
A very important trend that can be observed is that the increase in error due to an inaccurate estimation of the load frequency is less pronounced for smoother damping profiles.This is effect is particularly visible in Figure 13 (d  3), where the error measures barely change in the range 1 < CALM R   < 10.This results suggest that using smooth damping profiles leads to a much more robust implementation of the CALM, which is particularly useful for complex problems where there is no dominant frequency of the waves.Structures 12 (2015) 1540-1564 For d  1, the efficiency of the absorbing layer is less affected when the frequency of the load is higher than what the absorbing layer is prepared for.Previous work (Rodrigues and Dimitrovová, 2014) has shown that the absorbing layer is less efficient at absorbing low frequency elastic waves, which explains these results.It is therefore preferable to underestimate the dominating frequency of the loads than overestimating it.

CONCLUSIONS
The Caughey Absorbing Layer Method has been shown to work effectively to mitigate the problem of wave reflections at the boundaries for one-dimensional models.
Three different error measures were used to analyse the influence of various parameters in the results, with good agreement between all three.The tested parameters were: the absorbing layer length in relation to the wavelength to be absorbed (a ), the loss factor at the end of the layer (  ), the damping profile variation (constant, linear, quadratic, cubic and exponential) and the ratio between the frequency of the load and that of the absorbing layer ( CALM R

 
).It was shown that the linear damping profile led to better results when the absorbing layer is relatively short ( a ≤ 1), but for longer layers the quadratic damping profile performed better (1 4 a   ).It is possible that the cubic variation performs better for longer layers.The exponential damping profile did not perform as well for the CALM as Oskooi (2008) has reported for the PML.It was as effective as the quadratic profile for short layers (a < 1), where the linear profile was already a better choice, but for longer layers it was outperformed by all the tested profiles except for the constant loss factor.

Latin
For situations where the frequency of the load is incorrectly estimated, it was seen that increasing the smoothness of the damping profile increased the range of frequencies for which it was effective.Although the cubic profile performed worse than the linear and the quadratic ones for tuned absorbing layers ( CALM R    ), it showed a higher tolerance for detuned damping coefficients ( R

  CALM 
).This may make the cubic (or higher order) damping profile a better choice for problems with a wide range of frequencies to absorb, or when it is difficult to estimate the predominant frequency.
Empirical formulas for the optimum loss factor (   ) as a function of the length of the layer in relation to the wavelength to absorb (a ) were proposed.Each combination of error measure and damping profile variation resulted in a different formula, with an average coefficient of determination 2 0.98 R  .

Future work
The CALM was tested for two-dimensional problems with success by Semblat et al. (2011).The authors of this paper have already published results for the plain-stress two-dimensional problem (Rodrigues and Dimitrovová, 2014), where the method was shown to be more efficient than the Lysmer viscous boundaries.The plane-strain two-dimensional problem is currently being studied, with promising results.
The empirical approximations of the optimum loss factor will be tested in two-dimensional problems, to validate their applicability and verify if they are independent of the properties of the problem (except for the damping profile and length of the absorbing layer in relation to the wavelength).
An important development to pursue would be to implement the finite elements with independent Rayleigh damping in the Ansys software, so the method could be applied to more complex problems, including three-dimensional models.
Lastly, the theoretical analysis of the CALM must be further developed.In first place, the assumption that the decay in amplitude of the reflected waves is minimal for the desired frequency was not verified for the case study.Second, the theoretical prediction of the coefficients of reflection provided a reasonable approximation, but without the accuracy required to propose a theoretical optimum loss factor.
21) where 1  and 2  are the mass densities and ,1 P c and ,2 P care the pressure wave velocities for the original medium and the new medium, respectively.This formula assumes a steep transition, and therefore can be only applied for   a longer layer results in a smoother transition); increasing d increases RT R but reduces T R .This is in accordance with the numerical results observed bySemblat et al. (

Figure 4 :
Figure 4: Discrete damping elements for the absorbing layer.ex

Figure 5 :
Figure 5: Displacement across the medium at (a) 4.0 P t c   and (b) 6.6 P t c   .

Figure 7 :
Figure 7: max u for the optimum loss factor as a function of the layer's length.

FigureLL
Figure 8: 2 max L for the optimum loss factor as a function of the layer's length.Latin American Journal of Solids and Structures 12 (2015) 1540-1564 Optimum loss factor as a function of the layer's length and damping profile variation.

Figure 14 :
Figure 14: Error measures as function of the ratio of the load to layer frequency (  

Table 1 :
Numerical optimum loss factor as a function of the layer's length (d = 0).

Table 2 :
Numerical optimum loss factor as a function of the layer's length (d = 1).Figures 7 to 9 show the value of each error measure as a function of a for the corresponding optimum loss factor, for all the tested damping profiles (both axes are in logarithmic scale).All figures show some overall trends: