Extended Rayleigh Damping Model

In dynamic analysis, frequency domain analysis can be used if the entire structure is linear. However, time history analysis is generally used if nonlinear elements are present. Rayleigh damping has been widely used in time history response analysis. Many articles have reported the problems associated with this damping and suggested remedies. A basic problem is that the frequency area across which the damping ratio is almost constant is too narrow. If the area could be expanded while incurring only a small increase in computational cost, this would provide an appropriate remedy for this problem. In this study, a novel damping model capable of expanding the constant frequency area by more than five times was proposed based on the study of a causal damping model. This model was constructed by adding two terms to the Rayleigh damping model and can be applied to the linear elements in the time history analysis of a nonlinear structure. The accuracy and efficiency of the model were confirmed using example analyses.


INTRODUCTION
With the recent improvements in computer performance and simulation methods, a large number of response analyses using finite element methods (FEM) with a scale of many thousands of elements have been conducted, including Onimaru et al. (2012) and Nakamura et al. (2012).
Various materials for soil and buildings were used in these analyses. This study investigated the damping models used for these materials.
In general, the internal damping of a material does not depend on the frequency. A damping model with frequency independence is called a hysteretic damping (or structural damping) model. In frequency domain analysis, a complex damping model is often used because it satisfies the condition that the damping is independent of the frequency. The complex stiffness of such a model is described as K(1 + 2h·i), where i is an imaginary unit, h is the damping ratio, and K is the stiffness. Therefore, h is constant with ω, as shown in Figure 1.
However, this method can only be used when the entire structure is linear; it cannot be used if even a single non-linear element is present. To consider the effects of the non-linearity of a structure, time history response analyses of complete models must be performed, even if most of the elements remain in a linear condition. The objective of this study is to propose a novel damping model that can be applied to the linear elements, which are used with non-linear elements in a model. The Rayleigh damping model is often used for such cases. This model has a mass proportional part and a stiffness proportional part. In this model, h can be defined for two frequencies, and it is approximately constant between these frequencies. This model is therefore able to express limited frequency independence. Furthermore, the computational cost of the model is small. Other formulations of damping are able to define a greater range of values of h, e.g., Chopra (2001). However, these damping matrices [C] are complete, which greatly increase the computational cost, making them seldom used. Many recent studies have reported the effects of Rayleigh damping on inelastic response history analysis. Hall (2006) studied the problems of Rayleigh damping and noted that the damping force can become unrealistically large and the analysis response too small when during non-linear analysis softening occurs. He suggested that the causes lie in both parts. In the mass proportional part, the damping force becomes too large when the actual resonant frequencies decrease from the initial frequencies due to a reduction in stiffness. In the stiffness proportional part, when the restoring force is limited, the damping force becomes unrealistically large compared with the restoring force.
To address these problems, he proposed a capped viscous damping model that uses only a stiffness proportional damping part with an initial stiffness and a force limiter. Tangent stiffness proportional damping is not recommended because it produces a large change in the damping force and may cause convergence difficulties. Ryan and Polanco (2008) also studied the problems of the Rayleigh damping model for a base-isolated building and recommended the use of only stiffness proportional damping. Charney (2008) reported that excessive damping force can be created when initial stiffness is used in the non-linear condition and suggested the use of tangent stiffness instead. Zareian and Medina (2010) pointed out the same problem. Jehel et al. (2013) compared the effect of using either initial structural stiffness or updated tangent stiffness and studied the control and design of a Rayleigh damping model. Erduran (2012) evaluated the effects of Rayleigh damping on the engineering parameters.
All these studies proposed remedies for the stiffness proportional part of Rayleigh damping. However, no remedies have been proposed for the mass proportional part, other than eliminating it altogether. This is because the frequency area where h is almost constant (hereafter referred to as the constant area) for Rayleigh damping is too narrow. Hall (2006) noted that "A small Δ (tolerance of h) and a large R (constant area where the accuracy of h is under the tolerance) are desirable, but competing goals. " If a new damping model can expand the constant area of Rayleigh damping while imposing only a small increase in the computational cost, it will provide a good solution to this problem (see Figure 2).
In an earlier study, we proposed a range of causal hysteretic damping models (hereafter referred to as a causal model), which can approximately satisfy the requirement of frequency independence (Nakamura, 2007).
In this study and based on these earlier studies, a new damping model is proposed that is able to expand the constant area by more than five times by adding two terms to the Rayleigh damping model. In this paper, the Rayleigh damping model is described, and its accuracy is discussed. Next, an outline is given of causal damping. A new model called extended Rayleigh damping is then proposed, the results of a linear analysis using this new model are reported, and its accuracy and efficiency are assessed.

Outline of the Model
The equation of motion can be given by Eq. 1, where [M s ], [K s ], and {u s (t)} are, respectively, the mass matrix, stiffness matrix, and displacement vector for a structure, and {D s (t)} is the damping vector.
The Rayleigh damping model from Eq. 1 can be expressed by Coefficients α and β are given by Eq. 3 using (ω 1 , h 1 ) and (ω 2 , h 2 ), where ω is the circular frequency, and h is the damping ratio.
As shown in Eq. 2, the Rayleigh damping model has a mass proportional part and a stiffness proportional part. The damping ratio of the former is inversely proportional to ω and that of the latter is proportional to ω. As this paper is concerned with the condition in which the damping is independent of the frequency, only the case of h aim = h 1 = h 2 , is considered, where h aim is the target damping ratio. From this, Eq. 3 is changed to the simple form in Eq. 4.

Accuracy of the Rayleigh Model
The goal of this study is to propose a model that is able to realize a constant damping ratio h aim across a wide frequency area in time history response analysis. However, in the Rayleigh damping model, only two points of ω 1 and ω 2 satisfy the requirement that h(ω) is equal to h aim . We therefore studied the frequency area where the difference of h(ω) is within the tolerance by considering the tolerance of h(ω) against h aim . "Tolerance, " in this case, refers to the allowable upper and lower limits from h aim . Limits of 5 and 10% were used in this study. Figure 3 shows R h (ω) at a tolerance of 10%. The values for ω 1 and ω 2 were set such that the lower limit of R h (ω) was 0.9. This figure gives each circular frequency for ω 1 as a ratio to ω = 1. Table 1 shows the values of the constant area at h aim tolerances of 5, 10, and 20%. Here, W h is the ratio of the maximum circular frequency ω max to the minimum circular frequency ω min . It can be seen that the difference in the damping ratio became less than the tolerance as W h became larger. At tolerances of 5, 10, and 20%, the values for W h were 2.49, 3.71, and 6.88, respectively. These values cannot be said to be satisfactory. In Hall (2006), the tolerance and W h were notated as Δ and R, respectively, and it was observed that a small tolerance and a large W h are desirable. However, they are competing goals for the Rayleigh model.
The resonant frequency of a structure depends on the damping ratio h(ω). The resonant frequency ω ′ 0 with damping is different from the undamped resonant frequency ω 0 and is expressed by Eq. 7. When the accuracy of h(ω) for h aim is low, the accuracy of the resonant frequency also decreases. The accuracy of the resonant frequency is expressed by R res (ω) in Eq. 8.
From the above, although the accuracy of R res (ω) is high, the width W h is small, and the accuracy of the Rayleigh damping model is unsatisfactory. A simple and practical model with a wide constant frequency area for the damping ratio is therefore desirable. This study proposed a new model, by extending the Rayleigh damping model, to expand the width W h largely and to improve effectiveness greatly.

OUTLINE OF THE CAUSAL MODEL
This chapter describes the causal damping model used in this study. It is well known that the frequency-dependent function can be expressed as a series of impulses in the time domain. This Duhamel integral allows the frequency-dependent force to be expressed as the sum of the current displacement and the past displacements with each coefficient (Nakamura, 2012). "Almost frequency independent damping force" is another type of a frequency-dependent function that can be expressed in the same form.

Outline of the Model
In the complex damping model used in the frequency domain, the relationship between the reaction and displacement can be expressed by Eq. 9 in a frequency independent form, where F(ω) is the reaction, u(ω) is the displacement, K 0 is stiffness, h is the damping ratio, and i is an imaginary unit. Equation 10 is obtained by transforming Eq. 9 into the time domain. However, the calculation of I·u(t) is challenging.
Equations 11 and 12 express the unit imaginary function Z(ω) with 0 as the real part and 1 (−1 in the case of ω < 0) as the imaginary part. However, it is well known that this function is not causal and cannot be transformed into the time domain (Inaudi and Kelly, 1995).
A number of studies on a hysteretic damping model that satisfies the causality have been carried out. For example, Inaudi and Kelly (1995) presented a precise hysteretic damping model in the tine domain based on the Hilbert transform and proposed a method for applying it to time history analyses by conducting convergent calculations. Makris (1997) proposed a causal model that can be placed at the limit of viscoelastic model of Biot (1958) and showed that this model is closer to the ideal hysteretic damping than the Biot model. Afterward, Makris and Zhang (2000) compared the model with the Biot model minutely and showed that the Biot model is more practical and reliable than the model for actual time-domain analyses. They also approximated the Biot model simply by the Prony series and showed the efficiency of the model by applying it to the response analyses of a soil structure. Spanos and Tsavachidis (2001) proposed two methods to minimize the computational burden of the Biot model used in a timedomain analysis. One uses a recursive algorithm, and the other uses digital filters. Muscolino et al. (2005) applied the Laguerre polynomial approximation method to the Biot model in order to turn the original integrodifferential equations of motion into a set of differential equations, and the effect of the method was confirmed in both deterministic and random excitation analyses. Chouw (2002) used the causal damping model defined in the Laplace domain for the non-linear problem. Nakamura (2007) also proposed a causal hysteretic damping model using the Hilbert transform pair and showed the high accuracy and the effectiveness of the model. In this paper, this model is used as the causal damping model.
In order to make the function causal, the real part and the imaginary part must be a Hilbert transform pair. Recalculation of the imaginary part using the Hilbert transform allows the real part [Z k (ω)] to be obtained. It is no longer 0, but it instead becomes a frequency-dependent value [hereafter referred to as Z ′ R (ω)]. Z ′ (ω) in Eq. 13 is expressed using Z ′ R (ω). We call this value, the causal unit imaginary function (Figure 4).

FIGURE 4 | Calculated causal unit imaginary function Z ′ (ω).
The area between the upper and the lower limits (−ω lim < ω < ω lim , ω lim = 2πf lim ), where the imaginary part becomes constant, must be set, allowing Eq. 9 to be approximated by Eq. 14. This area can also be expressed as 0 < f < f lim (Hz) using f lim (Hz) by bringing the positive frequency part into focus. Hereafter, this area is called the focused frequency area, and ω lim and f lim are referred to as the upper limit circular frequency and the upper limit frequency, respectively.
The causal function Z ′ (ω) can be transformed into the impulse response function Z ′ (t) in the time domain with good accuracy. This gives an approximate value for the imaginary unit i in the time domain. Equation 14 is rewritten as Eq. 15, which uses a convolution integral in the time domain.
When using the time-domain transform method (Nakamura, 2006), Z ′ (t) is obtained as a series of discrete values. The integral from Eq. 15 is expressed by Eq. 16, whereu(t) is velocity, a 0 and b j are coefficients obtained from the time-domain transform, and Δt is 1/f lim (s).
Equation 17 is obtained by substituting Eq. 16 into Eq. 15. This is the target equation in the time domain. The first term on the right side of this equation corresponds to the ordinary restoring force, and the second term corresponds to the damping force. As a result, the reaction F(t) at the current time t is expressed as the sum of the coefficient products of the current displacement u(t), the current velocityu(t), and past displacements u(t − j·Δt).
In the frequency domain, Z ′ (ω) can be expressed by Eq. 18 using a 0 and b j . F(ω) can be calculated by substituting Z ′ (ω) into Eq. 14.
Equation 17 is rewritten as Eq. 19, where {F E (t)}, {u E (t)}, and [K E ] are the element reaction vector, the element displacement vector, and the element stiffness matrix, respectively. This causal model is clearly independent of the element type, and different types of element can be employed in the same way.
July 2016 | Volume 2 | Article 14 Nakamura (2008) proposed a method for simultaneously calculating the real part and the impulse response in the time domain with only the imaginary part of the complex function. This was done using the Kramers-Kroning relation (Nakamura, 2007). This process also provided a method for obtaining Z ′ (t) by carrying out time-domain transform directly, using only the imaginary part Z I (ω) and without conducting the numerical Hilbert transform. This method yields a model with higher accuracy because the omission of the numerical Hilbert transform reduces the calculation error. This model was used in the current study.
Equation 20 was obtained by expressing Eq. 19 as the damping force vector {D S (t)} of Eq. 1 (the equation of motion). In this equation, the stiffness proportional model with coefficient β of the first term is modified by the sum of the coefficient product terms of the past displacement of the second term (hereafter referred to as the time delay part term). Equation 21 shows the corresponding damping vector {D S (ω)} in the frequency domain.
(21) Figure 5A shows an image of the causal damping model. In the frequency area considered (0~ω lim ), the damping ratio h(ω) is almost constant. However, the characteristics are different at the two ends of the area, as h(ω) changes from 0 to h aim in the vicinity of ω = 0, and from h aim to 2h aim near ω = ω lim . This model is calculated as the sum of the stiffness proportional damping part with h aim at ω lim /2, shown in Figure 5B, and the term of the time delay part, shown in Figure 5C. The former corresponds to the first term Eq. 20, and the latter corresponds to the second term. Figure 6 includes the outside of ω lim . The causal model can be regarded as being made by changing the stiffness proportional part of the Rayleigh damping model to a step form. Table 3 shows the coefficients (a 0 and b j ) of representative cases, where ω lim is the upper limit circular frequency for the 9-term, 4term, and 2-term models. The models are identified by the time delay part b j , and are labeled the C9, C4, and C2 models. Models with a higher number are more accurate. The value of a 0 for all the models is the same as 1/(πf lim ), and Δt is 1/f lim (s).    Table 4 describes the condition of each model in the timedomain transform. The "frequencies used for transform" of each model are expressed as 0~1 (normalized with ω lim ). In this frequency area, the time-domain transform requires that the imaginary part become 1, and the imaginary part is interpolated smoothly.

Coefficients and Characteristics of Representative Models
Although Δt does not need to agree with the time step ΔT for a time history response analysis, Δt should be an integer multiple of ΔT. For example, at values of 0.1 s for Δt and 0.01 s for ΔT, the values at every 10 steps of the results obtained from the response analyses need to be used to derive the past displacement vector {u s (t − jΔt)}. If ΔT is changed in the analysis, it is unnecessary to change Δt. Figure 7 compares the Z ′ (ω) of each model. This figure was obtained by transforming the coefficients in Table 3, which were obtained from the transform to the time domain, back to the frequency domain. The values of the imaginary parts of all models are 1 at the frequencies used for the transform shown in Table 4. Accordingly, the area in which the value of the imaginary part is close to 1 is 0.05ω lim~0 .95ω lim for the C9 model and 0.2ω lim~0 .8ω lim for the C2 model.
In the real part, the maximum values at ω lim /2 are almost the same across all models. However, near the ends (0 or ω lim ), the value for the C9 model decreases sharply, while that of the C2 model decreases more slowly.
The damping ratio of the causal model, h ′ (ω), was obtained by Eq. 22 using the ratio of the imaginary part to the real part of 1 + 2hZ ′ (ω) (Nakamura, 2007).  The accuracy R h (ω) for the target damping ratio h aim was calculated using Eq. 23. Figure 8 shows the R h (ω) of the C9 and C2 models. In this figure, R h (ω) when h aim = 1, 3, 5, and 10 is shown because it changes as h aim changes.
Overall, the shapes of R h (ω) resembled those of Z ′ I (ω). When h aim was small, the accuracy tended to be high, and R h (ω) was therefore close to 1.
Next, a comparison between the C9 and C2 models was carried out, and the accuracy of the C9 model was shown to be higher overall, with a value of R h (ω) close to 1 in the area between 0.1ω lim and 0.9ω lim , whereas that of the C2 model was close to 1 in the area between 0.2ω lim and 0.8ω lim . The accuracy of R h (ω) in the C9 model decreased locally at both ends of the area. This became more notable as the value of h aim increased. In contrast, this tendency was hardly observed in the C2 model. These characteristics are produced by the shapes of Z ′ I (ω) at the ends of the frequency area.
Furthermore, when the value of h aim was large, the accuracy value was less than 1 across a wide frequency area because of the comparatively large value of Z ′ R (ω) in the central region.
The accuracy of the resonant frequency R res (ω) of the causal model was then investigated. R res (ω) was calculated using Eq. 24, obtained by multiplying the value of Eq. 25 by Eq. 8. Figure 9 shows the values of R res (ω) in the C9 and C2 models. The shapes of R res (ω) resembled those of Z ′ R (ω). The accuracy became higher as the value of h aim decreased. Across all models, when h aim was large, R res (ω) was greater than 1 in the central part of the frequency area around ω lim /2. In contrast, R res (ω) tended to fall below 1 at both ends of the frequency area. In the C9 model, this decrease in R res (ω) was significant, whereas the reduction in

Improvement of the Causal Model
This decrease in accuracy R h (ω) area around ω lim /2 when h aim is large was addressed next. Figure 10 gives an outline of the modification method. Figure 10A shows the hollow on the low-frequency side from the center being supplemented by the stiffness proportional damping. As a result, the accuracy decreases slightly on the high-frequency side from the center of the focused frequency area. The damping ratio h ′ (ω), which is obtained by a response analysis for the target damping ratio h aim , can be approximated using Eq. 26. At the central part of the frequency area where the accuracy is lowest, Z ′ I (ω lim /2) becomes a 0 ·(ω lim /2). Then, a 0 is increased ( 1 + 2h aim · Z ′ I (ω lim /2) ) times. A new value for a 0 (assuming a 0a ) is used in Eq. 27, with a 0 = 1/(πf lim ). Table 5 shows the Z ′ R (ω lim /2) value in each model. When h aim is large, an error arises from the difference between Eqs 22 and 26.
From a practical point of view, the use of a 0 (assuming a 0b ) from Eq. 28 is recommended. This equation is obtained by correcting the difference using a quadratic expression for h aim . In the numerical study, cases where h aim was less than or equal to 10% were investigated. From the results, 1.5 and 3.7 were selected as optimum coefficients.

Figure 11
shows R h (ω) calculated using a 0b for the C9 and C2 models. Since R h (ω lim /2) was almost 1, the effect of the correction was confirmed. This correction has little effect on R res (ω).

PROPOSAL FOR A NEW MODEL
As noted above, the Rayleigh model is simple and effective, but the frequency area in which the damping ratio can be regarded as constant is narrow. In this section, we propose a novel model, which has a wide constant frequency area.

Outline of the Proposed Model
Rayleigh damping is a combination of mass proportional damping and stiffness proportional damping, as illustrated in Figure 12. In our new model, the stiffness proportional damping is replaced by causal damping (Figure 13). The damping ratio h ′ (ω) can be approximated by Eq. 29 using Eqs 5 and 26. This allows the constant frequency area to be expanded. Therefore, the new model can be thought of as adding Rayleigh damping to the time delay term.
Frontiers in Built Environment | www.frontiersin.org July 2016 | Volume 2 | Article 14 The new model was constructed as follows: 1) In the constant frequency area, two models with damping ratios of 5 and 10% tolerance were produced. These are referred to as the Extended Rayleigh High accuracy model (ER-H) and Extended Rayleigh Middle accuracy model (ER-M).
2) The indicator of the constant frequency area W h should be more than five times that of Rayleigh damping. 3) In order to prevent the calculation cost from increasing, the C2 model was used, as it has the smallest number of time delay part terms. The damping force vector in the time domain was therefore expressed by Eq. 30, based on Eq. 20. 4) Equation 31 was obtained by setting the coefficients C 0 , C 1 , and C 2 . These values were obtained by numerical trial and error maximization of W h within the damping ratio tolerance. This equation was used in the time history response analysis. a 0 , b 1 , and b 2 are the coefficients of the C2 model in Table 3.
Δt implies 1/f lim (s), and Δt is different from the analysis time step ΔT. 5) Equation 29 was then rewritten as Eq. 32, and Eq. 32 was expressed as Eq. 33 using the coefficients C 0 , C 1 , and C 2 . For the representative damping ratio h aim , four values were selected: 1, 3, 5, and 10%. In each case, the coefficients C 0 to C 2 were established by carrying out numerical trials using Eq. 33 such that the value of W h was maximized. 6) Cases in which h aim was placed between these four values were investigated, and all models were adjusted to obtain favorable results by performing linear interpolation of these coefficient values.
where  Table 6 shows the recommendation values for C 0 , C 1 , and C 2 with regards to each h aim of the ER-H model and the ER-M model. The coefficients in the cases 1~3%, 3~5%, and 5~10% for h aim were obtained using linear interpolation. C 0 exhibited a strong effect in the low-frequency area, C 1 in all frequency areas, and C 2 in the high-frequency area. These are independent values from Δt and ΔT.

Confirmation Using Multiple Mass System Model
To test the proposed model in a time history response analysis, the case of 1~100 Hz, assuming f lim = 100 Hz, was investigated. Figure 14 shows the analysis model, and Table 7 gives its properties. One hundred single degree of vibration models were installed on a rigid beam, which was connected to the fixed end of a rigid spring. The masses were adjusted so that the value for each undamped resonant frequency was 1, 2, 3, . . ., 100 Hz. A time history response analysis was conducted by providing a unit pulse wave as the ground motion acceleration. The Newmark-β method (β = 1/4) was employed for the time integral. The analysis time step ΔT was set at 0.0005 s to allow the response in the high-frequency area to also be evaluated.
The transfer function was calculated by dividing the response acceleration of each point by the ground motion acceleration in the frequency area, to examine the peak height of each transfer function. It is known that the theoretical value p with regards to the peak height of the transfer function for a one degree of freedom (DOF) system can be expressed by Eq. 34 as a function of the damping ratio h. Based on this, the obtained damping ratio h(ω 0 ) was calculated using Eq. 35 as a function of the peak height of each single DOF model p(ω 0 ), where ω 0 is the undamped resonant circular frequency. By comparing the result of this equation with the target damping ratio h aim , the accuracy of the damping ratio R h (ω 0 ) could be calculated.
The theoretical value of the resonant circular frequency of the vibration system with damping ω ′ 0 is given by Eq. 7 as the function of the undamped system resonant circular frequency ω 0 .   The accuracy of the resonant frequency realized for each model R res (ω 0 ) could then be calculated by comparing the circular frequency corresponding to the peak value of each vibration system computed in the response analysis with this theoretical value. Figure 15 shows the results from the analyses carried out using the Rayleigh damping model. The values for h aim were set at 3 and 10% for ω 1 and 10 and 25.5 Hz for ω 2 . This corresponds to the case with 10% tolerance shown in Table 1. White and black circles show the theoretical values at h aim = 3 and 10%, respectively. These values were calculated for R h (ω) using Eq. 5 and for R res (ω) using Eq. 7. The results calculated by the response analyses are shown in the figure as solid and dotted lines. These calculated results corresponded well to the theoretical values. The validity of the analyses in this section was confirmed. R h (ω) did not depend on h aim , and there was a large difference between R h (ω) and h aim in most of the frequency area, with the exception of ω 1 and ω 2 . As a result, R res (ω) was consistent with Table 2.  Figure 16 shows the results of the analysis of the causal model improved using Eq. 28. h aim were set at 3 and 10%. R h (ω) corresponded to Figure 11 and R res (ω) to the C9 model shown in Figure 9. The results for R h (ω) and R res (ω) were equivalent to Figures 9 and 11, respectively.  Figures 17 and 18 show the results of the high accuracy model (ER-H) and the middle accuracy model (ER-M) with h aim at 3 and 10%. In ER-H, all values of R h (ω) were within 0.95~1.05 in the area from 6 Hz to~80 Hz, a tolerance within 5%. In ER-M, all the values of R h (ω) were within 0.90~1.10 from 4 Hz to the vicinity of 85 Hz, giving a tolerance within 10%. The values of R res (ω) in these figures correspond well to those of the C2 model shown in Figure 9. Table 8 shows the constant frequency area confirmed in these analyses. The value of W h from the ER-H model was greater than 13.0. This is 5.3~5.5 times larger than the value from the Rayleigh model with the same tolerance. In the case where h aim = 3% and f lim = 100 Hz, for example, the frequency area from 6~82 Hz was

EXAMPLE ANALYSIS
To confirm the applicability of the proposed model to a practical structure, investigations were carried out using a 3-dimensional finite element model (hereafter referred to as 3D FEM).

Investigations Using a 3D FEM Model
To confirm the effectiveness of the proposed model as a practical multiple DOF model, investigations using a 3D FEM model were carried out. The analysis model was a dome-shaped structure with a diameter of 44 m, a height of 58 m, and a weight of 1400 t, as shown in Figure 19. The model was divided into 672 elements using shell elements, with 32 divisions in the circumferential direction. In the vertical direction, there were 10 divisions in the cylindrical part and 5 divisions in a dome part. The lowest points of the cylindrical part were fixed. Table 9 shows the principal properties of the model. Table 10 shows the principal natural frequency obtained from real eigenfrequencies, and the eigenmodes are given in Figure 20.
The conditions for the time history response analysis were the same as those in the previous section. The time step ΔT was set at 0.0005 s to estimate the high frequency. The response acceleration at the nodal point of the top was obtained by providing a unit impulse in the horizontal direction as ground motion acceleration.
Transfer functions were calculated by dividing this response acceleration by the ground motion acceleration in the frequency area. Values of 3 and 10% were used for the damping ratio h aim of the materials. The transfer functions were horizontal response outputs for horizontal direction inputs.
A frequency domain analysis (FDA) was carried out. This was assumed to give the correct answer to the problem, and the results were compared with those from the proposed model. Constant damping ratios of 3 and 10% were provided by the complex damping. The analysis pitch Δf in the frequency area was set to 0.1 Hz. Figure 21 shows the results from the Rayleigh model, in which Case 1 and Case 2 (corresponding to Rayleigh-1 and -2) FIGURE 19 | Analysis model.   were investigated. In Case 1, a horizontal first-order frequency of 6 Hz was set for ω 1 , and a vertical first-order frequency of 14.6 Hz was set for ω 2 . In Case 2, the horizontal first-order frequency was set for ω 1 , and ω 2 was set to 10 times ω 1 . The results for Case 1 were very close to the FDA results in the area from 0 to nearly 20 Hz. However, in the area of 20 Hz and above, a large difference from the FDA emerged due to excessive damping. In Case 2, the results corresponded closely to FDA in the range 40~80 Hz because the ω 2 of this model was 60 Hz. However, a large difference in the peak values at the horizontal second-order frequency was seen in the vicinity of 15 Hz. These results confirmed that it is difficult for the Rayleigh model to express a constant damping ratio across a wide frequency area. Figure 22 shows the results obtained using the C9 model described in Section "Improvement of the Causal Model. " The value for f lim was set at 100 Hz. For the time delay part, nine terms for the displacement at 0.01, 0.02, 0.03, . . ., 0.09 s before the current time were considered. This model closely matched the response results of FDA in the area from 0 to almost 90 Hz when h aim = 3%, and the area from 0 to almost 80 Hz when h aim = 10%. However, when h aim = 10%, a difference in the peak height of the horizontal first-order mode in the vicinity of 6 Hz was seen. This was caused by the reduction in the damping accuracy at both ends of the constant frequency area. It occurred when the damping ratio was large in the C9 model (cf., the results from C9 in Figure 11). Figures 23 and 24 give the results from the ER-H and ER-M models. The f lim value for ER-H was set at 100 Hz, and the displacement at 0.01 and 0.02 s before the current time were used in the time delay term. The f lim value for ER-M was set at 150 Hz, and displacements at 0.0067 and 0.0133 s were used.
The results from ER-H corresponded quite closely to those from the FDA in the area from 0 to almost 80 Hz. Although the accuracy of ER-M was slightly lower than that of ER-H in deriving the height of the peak value, the results from ER-M corresponded well in general to the results from the FDA in the area from 0 to almost 130 Hz. This confirmed the accuracy of the proposed models.

CONCLUSION
When time history response analysis is carried out using FEM, a Rayleigh damping model is often employed to take account of the constant damping characteristics for the frequency. However, the accuracy of this model is unsatisfactory. In this paper, we proposed a new model in which the stiffness proportional damping term of the Rayleigh model was replaced by that of the causal damping model. In this approach, the term of the coefficient sum for the past two time displacements is added to the standard Rayleigh model. This greatly improved the accuracy of the Rayleigh model. The effectiveness of the new models was confirmed by carrying out a trial analysis.
Future studies will be conducted to investigate the applicability of the new model to a range of model types and to large-scale models, especially those with tens of thousands of elements. We will also investigate its applicability to the element Rayleigh damping model, which is used when a total model has multiple damping ratios. In this study, the investigations were conducted on a linear level. Further studies are needed on the applicability of the model when elements become non-linear and the value of h aim increases.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and approved it for publication.

FUNDING
This work was supported by JSPS KAKENHI Grant Number 26289197.