A method to accurately define arbitrary algorithmic damping character as viscous damping

Undesired oscillations often emerge in numerical solutions, especially in the case of dynamic problems. These are mainly spurious oscillations which must be eliminated or reduced at least to provide accurate results. Numerical methods with damping effect are especially useful to achieve this goal. However, the concrete shape of the damping characteristics has a great impact on the effectiveness. Dissipative numerical methods mostly have a specific damping character with very limited alteration possibility. In this article, a novel numerical method has been introduced where the dissipative effect is exerted via viscous damping. Using the proposed method, a great variety of damping curves can be defined accurately, straight through the determination of the algorithmic damping ratio. The newly developed technique is mainly useful for applications where the shape of the damping characteristics significantly affects the accuracy.

for calculations with hydrodynamic shockwaves. Metzger [11] tested his iterative technique successfully on dynamic relaxation problems.
Numerical methods can be effectively described using certain accuracy measures [12] which are the spectral radius, the relative period error and the algorithmic damping ratio. In this paper, the lastly mentioned indicator is going to be examined in detail which characterizes the damping effect of the applied method. For that purpose, the amplitude decay [13] is also used oftentimes. However, using this measure is less advantageous, as it cannot describe the damping characteristics. In the literature, many articles have been used an approximate scheme [14,15] for the calculation of the algorithmic damping. Nevertheless, it is beneficial to apply the accurate scheme [16] of the algorithmic damping ratio in all cases, as it only means a minor additional complexity in the calculation. It is important to emphasize that there exists a similarly named, but different measure which is the Physical Damping Ratio (PDR) [17]. This usually characterizes the viscous or material damping in the model; the difference between algorithmic damping and PDR is thoroughly described in [18].
As it has been mentioned earlier, one possibility to exert a damping effect is the numerical damping which is contained in the algorithm. In case of transient problems solved with FEM, the time integration is mostly realized by time stepping methods [19]. There exist both forward and backward increment, single-stage and multistage schemes, while the damping effect may be alterable with one or more parameters. The shape of the damping character is also an important aspect; it can be both progressive, degressive and linear. The Newmark method [5] is a widespread backward increment (implicit) scheme whose damping can be controlled with one parameter. Due to its simple formulation and unconditional stability, this method is built in many FE simulation software like Abaqus or Ansys. The HHT-α method [20] is the improved version of Newmark's scheme where the degressive damping character can be set with an additional parameter. In contrast with that, Bathe developed an implicit method [21] whose damping cannot be set by parameters; it only depends on the actual time step size. The Wilson-Θ method [22] is also a well-known implicit scheme where damping levels can be set substantially higher compared to similar methods. Another widely applied implicit method is the generalized-α scheme [23] with alterable, slightly progressive damping character. The formerly mentioned method from Hulbert and Chung [7] is the forward increment (explicit) modification of the generalized-α scheme with more possible adjustment in damping. Among explicit techniques, Kim has provided a novel two stage method [8] with adjustable progressive characteristics. In Chang's method [14], an interesting approach is applied, as the control parameter is generated from a finite sum containing material parameters. Thus, the shape of the damping character can be altered in a wide range, from slightly degressive to heavily progressive. Besides, there exist methods with decreasing damping character like Heppler and Hansen's scheme [24] and even methods with partially negative damping character like Wang and Au's scheme [25]. However, the shape of the characteristics and the maximal damping value is very limited by most of the above-mentioned methods. Thus, it is often difficult to find a scheme that possesses the proper features for the actual application.
Another approach of exerting a damping effect is the application of viscous damping in the FEM model or the addition of a viscous pressure term to the dynamic equations [26]. A solution of this kind was firstly described by VonNeumann and Richtmeyer [6] in 1950; therefore, it is mostly known as VonNeumann's artificial viscosity or bulk viscosity method. This scheme is based on the Central Difference Method (CDM) [27] which does not have numerical damping, so the amount of dissipation can be controlled through only one parameter. The damping characteristics are linear, so the effect is equivalent with a stiffness matrix proportional Rayleigh damping [28]. This technique is built in several industrial software like Abaqus [29] or LS-Dyna [30]. A similar approach has been described by Munjiza et al. [31] who applied a certain kind of Caughey damping [32]. By using this method, both linear and progressive characteristics can be achieved with remarkably simple methodology. Metzger has provided a novel adaptive technique [11] for the determination of the damping matrix combined with CDM. This is an iterative method based on the eigenmodes of the model and is mainly practical for contact problems and large deformations. Lee elaborated a proportional viscous damping model [33] based on bell-shaped basis functions. Using this method, the numerical damping character can be fitted on known data points with a very good approximation. However, all of the above-mentioned methods can only offer an approximate determination for most damping character shapes.
The proposed method is based on the CDM, so the exerted damping effect can be derived straightly from the viscous damping term in the FE model. Using the newly developed method, a great variety of shapes are achievable as damping characteristics. As it has been mentioned above, most of the existing methods are only able to produce a certain kind of damping character and the available damping levels are also limited. Finding the appropriate damping curve is critical in many applications what makes the above-mentioned methods difficult to use in certain cases. The main goal of this paper is to provide a technique for the determination of the aimed damping character for applications where the shape of the damping curve is especially important.
In this article, a novel technique is presented to specify the algorithmic damping precisely in a wide range. The motivation of the research is to introduce a method by which arbitrary damping characteristics can be exerted accurately in a FE model as viscous damping. To achieve this, the differences between the PDR and the algorithmic damping ratio must be cleared for viscous damping. Deriving the relationship between these quantities is especially important to define an aimed damping character accurately. The main improvement in the research is the exposition of the distortion in the numerical solution when the aimed character of the viscous damping is given as equivalent PDR compared to when it is defined as equivalent algorithmic damping ratio.

Discretization in space and time
In the main part of the article, the formulation of the equation of motion shall be considered as first. The proposed method is based on the CDM which is responsible for the discretization in time. The spatial discretization is realized using the FEM. Methods with damping effect are mainly used for dynamic problems that can be characterized by a semi-discretized equation of motion written as a second-order ordinary differential equation according to where M, C and K denote the mass, damping and stiffness matrices, respectively. On the right hand side, f t means the load vector, u t is the nodal displacement vector,u t means the nodal velocity vector whileü t denotes the nodal acceleration vector at time instant t. Besides, the solution of equation (1), requires both boundary and initial conditions. The former depends on the applied model while the latter can be written as where u 0 means the initial displacement vector while v 0 denotes the initial velocity vector at t = 0.

The central difference method
The CDM is a forward increment time stepping method with a simple formulation and second-order accuracy.
On the one hand, calculations can be run very fast due to its explicitness and the lack of numerical damping is also a beneficial feature in our case. On the other hand, the CDM possesses only conditional stability limiting the available time step size. The formulation of this method can be given by the formula of the velocity and the acceleration at the current time step t the following way.
where t means the time step size. By substituting (3) and (4) into (1) and reordering it to express u t+ t , the following linear system of equations can be determined.
However, there is a serious problem with this formulation, as the inversion of the term tC + 2M is a computationally expensive procedure. Hence, it is worthwhile to replace equation (4) with the following formula to calculateu t .u After the substitution of (3) and (6) into (1) and the deduction for u t+ t , the equation below can be obtained.
This formulation is more advantageous, as it only requires the inversion of M. If the mass matrix is diagonalized [34], then the inversion of the resulting lumped mass matrix is computationally cheap. Besides, considering equation (3), (6) and (7) leads to the conclusion that the calculation ofü t ,u t and u t+ t only requires displacements at different time moments. This feature is especially important in applications (e.g., contact problems, simulation of impacts, etc.) where the acceleration has a steep shift in its time evolution. In order to analyze the accuracy indicators which are presented in the next section, it is practical to write equation (7) into the following form.
where A CDM means the amplification matrix. For the central difference method, this matrix can be given as Knowing the inhomogeneous component (b CDM ) of equation (8) is not substantial for the later examinations.

Accuracy indicators
In order to analyze the accuracy and the stability of the proposed method, three indicators are used which are the spectral radius, the relative period error and the algorithmic damping ratio. To apply these characteristic features, equation (1) must be transformed into the frequency domain, so it can be rewritten as a series of uncoupled Single-Degree-of-Freedom (SDoF) equations.
where ω i is the circular eigenfrequency, ξ i denotes the PDR and i means the index number whose maximal value is the Degree of Freedom (DoF) of the FE model. To determine these parameters based on the FEM model, the M, C and K matrices are transformed into the frequency domain using the transformation matrix denoted with T which can be composed from the eigenvectors.
By this transformation, the parameters ω i and ξ i are determined using equations (11)- (13). Based on the amplification matrix, the spectral radius can be computed according to where λ 1 and λ 2 mean the first and second eigenvalues of A CDM . The numerical solution remains stable only if spectral radius is less than or equals one. Belonging to the maximal permissible value of ρ, a stability limit can be determined in the time step size denoted with t c . This critical time step can be computed according to the following equation.
The bifurcation point is also an important measure that is frequently used to specify the maximal applicable damping level. If the stability condition is satisfied, the eigenvalues of A CDM can be written as according to Kolay and Ricles [16] where σ and ε mean the real and the imaginary part of the eigenvalues, respectively, ξ denotes the algorithmic damping ratio, while Ω = ω t where ω means the algorithmic circular eigenfrequency. The bifurcation point is defined as the time step size value where both λ 1 and λ 2 become real values (ε = 0). This particular value of time step is denoted with t b ; the spectral radius associated with this value is defined by The spectral radii of the CDM for various numerical damping values are shown in Fig. 1 where the stability limit and the bifurcation point are also marked on a particular curve. Knowing λ 1,2 given by (16), Ω can be computed as Using (18), ξ can be given by The ω quantity is the discrete equivalent of the exact circular eigenfrequency (algorithmic circular eigenfrequency), while ξ means the algorithmic damping ratio that has a similar connection to the PDR, as it is the discrete equivalent of ξ . This measure describes the numerical dissipation character of the examined time stepping method. This characteristics may be both computed as the function of the time step size for a given circular eigenfrequency and vice versa. Another important indicator by the examination of the accuracy is the relative period error which characterizes the numerical dispersion of the examined method in the function of the time step size. This indicator is defined as where the exact period is T = 2π/ω and the numerically computed period is T = 2π/ω. The relative period errors of the CDM for various numerical damping values are shown in Fig. 2.

Determination of the damping characteristics
In the proposed method, the damping characteristics are defined by the appropriate specification of the damping matrix in equation (1). Knowing the values of ξ in equation (10) for every single eigenmode, the C matrix can be created using the transformation matrix T. However, the connection between ξ and ξ must be examined to be able to accurately specify the aimed damping character. In order to establish the formula of ξ for the CDM based on (19), the SDoF version of A CDM must be constituted based on (10)- (13).
Knowing A CDM SDOF given by the equation (21), ξ can be calculated using equation (19) based on equation (18) and (16). As the algorithmic damping ratio traditionally means the fictitious numerical damping introduced by the numerical method, a different designation must be applied for ξ which includes the contribution from the physical viscous damping. Hence, the combined effect of the sheer algorithmic and viscous damping shall be mentioned as Apparent Damping Ratio (ADR). Thus, the ADR of the CDM can be given by the following two-stage definition.ξ Comparing the PDR with the matching ADR at identical parameters, a significant difference occurs between the resulting damping curves. In Fig. 3 a progressive, a degressive and a linear PDR character is plotted in the function of circular eigenfrequencies with identical initial tangent. As it can be observed, there appears a substantial amount of distortion in the ADR compared to the original curves. Thus, the aimed damping characteristics cannot be given as the PDR in the function of ω; it must be specified by the ADR. Hence, the corresponding ξ values must be calculated for every given ξ based on equation (22)- (23). As these formulae cannot be treated analytically, the belonging ξ values can be solved with the use of an appropriate root finding algorithm.
The proposed method provides a very versatile determination for the damping characteristics. According to the shape of damping curve, there are no constraints, as the damping values can be specified for each circular eigenfrequency individually. However, there exists a maximal damping limit for each ω which separates the convergent region from the divergence. To ascertain these particular damping values, t c and t b must be deducted for the CDM according to (15) and (17), respectively.
for each ω and the belonging PDR. At the bifurcation point, the imaginary part of the eigenvalues becomes zero (ε = 0). This condition results the following equation as the formula for t b .
Considering equation (24) and (25), it can be easily deducted that inequality t b ≤ t c is fulfilled for every ω and ξ . Hence, the maximal applicable parameters are ascertained at the bifurcation point. Below this limit value, the applicable PDR region can be calculated based on within the region of 0 ≤ ξ < 1. Examining the ADR at the maximal possible PDR determined by equation (26), it can be analytically deducted that the left side limit of ξ at this particular value gives the following result.
This limit value is also called critical damping. Thus, the maximal applicable ADR value is independent from the time step size and the circular eigenfrequency. This result provides a great freedom for the determination of the damping characteristics, as any value is possible within the region of 0 ≤ ξ < 1 for every ω, by any t < t max where t max is determined for the undamped case as The distortion between the PDR and the ADR is especially significant by ω values close to ω max (see Fig. 4). The applied time step size also affects the amount of distortion substantially as it can be concluded by comparing Figs. 4 and 5. Thus, it can be inferred that by time steps near t max and by frequencies close to ω max equation (22) and (23) must be applied to determine the damping characteristics accurately.

Model
In order to demonstrate the practical benefits of the proposed method, a simple numerical example is considered. In Fig. 6, the mechanical model of a wave propagation problem is shown. The right side endpoint of the rod of length l is fixed while a v 0 initial velocity is defined for every other point of the rod. Because of the fixed endpoint, an elastic deformation occurs alongside the x direction which can be described as a propagating wave. Similar examples are used oftentimes to test experimental numerical methods [2,8,35]. The above presented test example is solved using the FEM with both one-dimensional (1D) and twodimensional (2D) models. The FE model and the boundary conditions for the 1D case are shown in Fig. 7. This model is compiled using 1D truss elements with 1-Degree-of-Freedom (DoF) per node which is the longitudinal displacement. It contains uniform elements with E = 90 MPa and ρ = 7.85 · 10 −9 kg/mm 3 The FE model and the boundary conditions for the 2D case are shown in Fig. 8. This model is compiled using 2D quadrilateral elements with 2 Degrees of Freedom (DoF) per node which are the longitudinal and the transversal displacement. It contains uniform elements with the same material parameters and initial velocity like by the 1D model. The parameters of the mesh are the following: l = 100 mm, h = 10 mm, n = 100, m = 10, l e = h e = 1 mm.

Results
As it is proved in section 4, the definition of an aimed damping character as equivalent PDR results a significant distortion in the ADR (see Fig. 3). In order to test this effect in the above-mentioned 1D and 2D test examples, three different characteristics are prescribed according to Fig. 9. The progressive, linear and degressive characters are given in the form of f (x) = 2 x − 1, f (x) = x and f (x) = log 2 (x + 1), respectively. Considering the damping curves at t = 0.4 t max , the ADR equivalent to the aimed character offers a perfect match to the aimed character, while the equivalent PDR results a severe distortion in ξ by all of the considered characters.
Considering the curves in Fig. 9, it is easy to recognize that the distortion mainly occurs at the higher circular eigenfrequencies (which is consistent with Figs. 4 and 5). The amount of distortion shall be quantified using the relative error between the results of the equivalent PDR and ADR. This can be given in a comprehensive form according to where N means the number of time steps, while PDR i and AD i are the considered quantities at the ith time step by equivalent PDR and ADR, respectively. The relative error calculated for for ξ shows a growing difference between the equivalent PDR and ADR as the damping curve turns from progressive to linear and degressive. This phenomenon emerges because at higher ξ values, the discrepancy between the equivalent PDR and ADR will also be higher (see Figs. 4 and 5). Thus, in case of equivalent PDR, the shape of the damping curves remains close to the aimed character at the lower frequency range, but undergoes significant changes at higher frequencies.
The 1D representation of the considered wave propagation problem has an exact solution that can be deducted using the wave equation with the corresponding boundary conditions. This feature provides a very useful reference for the evaluation of the numerical results. Applying the damping characters shown in Fig.  9 for the 1D model, numerical results appear to have a significant difference compared to the exact solution (see Fig. 10). Comparing the definitions of the damping character as viscous damping using equivalent PDR and ADR, it can be concluded that the different ξ characters in Fig. 9 induces a noticeable discrepancy in the resulting velocity functions in Fig. 10. As it has been expected based on Fig. 9, the velocity curves run farther from the exact solution due to the higher amount of damping when the aimed damping character is defined as equivalent PDR.
The relative errors computed for the velocity show an increasing tendency as the damping character turns from progressive to linear and degressive. These results are in accordance with the distortions between the damping characters in Fig. 9. The higher discrepancy in the damping curves results a higher relative error in the velocity functions.
The 2D representation of the test example cannot be treated analytically; hence, there is no reference for the numerical results. However, comparing the two different definitions of damping, a significant distortion occurs in the resulting velocity functions by equivalent PDR. Besides, the 2D case results the same tendency in the relative errors as in the 1D case, but with slightly higher values. The latter is reasoned by the more complex FE model of the 2D case with higher eigenfrequencies (Fig. 11).
In order to verify the accuracy the numerical results, comparative simulations were made in ANSYS Mechanical APDL software for both the 1D and 2D cases. The FE models were defined using the same parameters as above. The element type LINK180 was defined in the 1D case, while PLANE182 was applied for the 2D model. The CDM method is not implemented in Mechanical APDL in general, but with module LS-DYNA, it can be applied. The appropriate damping matrix was pre-composed and imported. The comparative results showed a very precise math for both the 1D and 2D models, proving that our own FE code produces accurate calculations.

The effect of time step size and element number
In the sections above, the differences between the equivalent PDR and ADR were determined for one particular time step size and element number. However, E rel (PDR) must be determined for different temporal and spatial discretizations to examine the effect of the time step size and the element number. In the following, the damping character is defined in Fig. 9b for both the 1D and 2D cases.
Firstly, the 1D case shall be considered in the function of these two parameters (Fig. 12). The contour plot shows that the relative error is roughly independent from the time step size. However, the number of elements has a significant impact on E rel (PDR). It decreases as the spatial discretization becomes finer, but the gradient of decline also diminishes.
As computations are more resource intensive for the 2D model, the relative errors are examined in the function of the time step size and the element number separately. In Fig. 13a, the time step size is set as t = 8 · 10 −7 [s], while in Fig. 13b, n = 640 element number is considered. The results show a similar connection between the relative error and the number of element and the time step size, respectively. E rel (PDR) does not improve significantly with the reduction in the time step size. However, increasing the number of elements reduces the relative error with a diminishing gradient. An other possible application of the CDM with viscous damping is the reproduction of the damping effect of methods using numerical damping. To put this feature to the test, the well-known Newmark method [5] is considered. The basic formulation of this method is determined by the following two equations. where In order to accurately reproduce the numerical results of the Newmark method using the CDM with viscous damping, the best possible match is aimed in the numerical features of the two methods. Thus, not only a matching damping characteristics must be defined, but the PE also must be aligned to the Newmark method. Based on Fig. 14, it can be concluded that a perfect match cannot be achieved, as the PE curves for the Newmark method run in the positive, while for the CDM in the negative range. However, through the appropriate setting of the time step size, the PE values can be taken very close to each other, as it is suggested in Fig. 14. Thus, time step size is set to be t = 0.05 t max in this test example , while the control parameter is given as δ = 0.8. After setting the time step size appropriately, a matching damping characteristics must be defined for the CDM (Fig. 15). As the applied time step is a very small proportion of t max , the definition in PDR does not result a significant distortion in ξ (see Figs. 4 and 5). Thus, the Newmark method's damping character can be  given directly as the PDR for the chosen time step. In contrast, if the time step size is set to be much higher ( t = 0.3 t max ), then a noticeable distortion occurs between the plotted damping characters.
The accuracy of the Newmark method's reproduction is tested in both of the considered 1D and 2D examples. As it is expected based on Fig. 15, the resulting velocity curves (Figs. 16 and 17) run very close to each other by both models if t = 0.05 t max . On the one hand, the results confirm that the direct definition of PDR is sufficient in this time step range. On the other hand, the reproduction of the Newmark method proved to be accurate too. However, when t = 0.3 t max is chosen, the matching with the Newmark method deteriorates and a significant difference emerges between the equivalent PDR and ADR. Thus, it can be stated that the close match of the period error is also required to reproduce numerical damping accurately with the use of viscous damping.

Conclusion
In this paper, a versatile method has been introduced to determine damping characteristics with arbitrary shape. The proposed technique is based on the central difference method, as it has no numerical damping. The damping effect is straightly exerted via viscous damping in the FE model. As it has been proved, the proposed method provides a great accuracy in the composition of the damping matrix, as the damping characteristics are straightly defined through the ADR. Besides, this approach offers a great variety in the determination of the aimed damping character. In the shape of the damping curve, there is no limitation, as any value below the critical damping can be given for the ADR at any eigenfrequency and t < t max . The proposed method was tested in both 1D and 2D test examples where it has been proved that the damping characteristics at time steps close to t max must be defined as the ADR. Another useful feature of the presented method is the reproduction of numerical damping with the use of viscous damping. Nevertheless, this can only be fulfilled if a close match is accessible in the period error compared to the aimed numerical damping method. Further inspections suggest that the application of the proposed method with specific damping characters can be especially useful for contact problems.