Contents here is detailed explanation of our paper¹⁾ on rational complex modulus. It is composed of three collums, which are damping ratio of SHAKE, explained here, and Lysmer's proposal and its discussion issue, and rational complex modulus.

SHAKE²⁾ is a well-known computer program. It is developed in 1972, therefore, it has been used for about 50 years. However, as far as I know, there is no people who has a question in the followings. In addition, error of equation in ref. 2) was not pointed out. It is, of course, a big job to utilize complex modulus in the seismic response analysis of ground. It is also important to point out something not good.

Equation of motion in SHAKE manual¹⁾ is as follows.

\(\rho \cfrac{\partial ^{2}u}{\partial t^2}=G\cfrac{\partial ^{2}u}{\partial x^2}+\eta\cfrac{\partial ^{3}u}{\partial x^2\partial t} \)(1)

where \(u=u(x,t)\) denotes displacement, \(\rho\) denotes density, \(G\) denotes shear modulus, \(\eta\) denotes damping coefficient, \(x\) denotes coordinate axis to the downward, and \(t\) denotes time. Derivation of this equation will be explained later.

A method of separation of variable is employed to solve the equation of motion. The displacement if separated into space and time domain such as \(u(x,t)=U(x)e^{i\omega t}\). Then the solution in space yields as

\(U(x)=Ee^{ikx}+Fe^{-ikx}\)(2)

where \(E\) and \(F\) are amplitude of incident and reflected waves. Here,

\(k^2=\cfrac{\rho \omega ^2}{G+i\omega \eta}=\cfrac{\rho \omega^2}{G^*}\)(3)

is a wave number. The term \(G+i\omega \eta\) in the central part of this equation, \(\omega \eta\) is replaces as

\(\omega\eta=2G\beta\)(4)

Then we can get a famous complex modulus is obtained as

\(G^*=G(1+2i\beta)\)(5)

The problem is \(\beta\) in Eq. (4). In the manual of SHAKE, it is written as critical damping ratio. The critical damping ratio is already explained in Issue related to damping. Moreover, it is written in the all textbook dealing with vibration. The critical damping ratio \(h\) is defined as

\(h=\eta/(2\sqrt{mk})\)(6)

You can find that Eq. (4) and Eq. (6) are quite different; definition itself is different. For example, mass \(m\) is not included in Eq. (4). It is natural because Eq. (1) is derived from the vibration of one-degree-of freedom system, whereas Eq. (4) is a equilibrium equation of the soil element. In other word, where there is no quantity corresponding to the critical damping in Eq. (4); therefore, it is not a ratio

I think this is just a wrong idea by the developers. As explained in the Issues related to damping, when the values of \(\beta\) or \(h\) is sufficiently small, damping constant derived from the hysteretic behavior \(h\) and critical damping ratio \(\beta\) are same value. On the other hand, \(\beta\) defined by Eq. (4) yields damping ratio by using Eq. (5). Therefore, they seems to use the term critical damping ratio for \(\beta\) that is same value with damping ratio.

In the Issued related to damping , hysteretic damping in the one-degree of freedom system yields

\(h=\cfrac{1}{4\pi}\cfrac{S}{W}\),(7)

or \(h=\beta\) only when \(h\) is sufficiently small and the following equation holds.

\(e^{-h\pi}\approx1-h\pi\) (8)

As written there, this equation takes the first two terms of the Maclaurin expansion. The third term is \(\pi h^2/2\). Therefore, the condition that Eq. (7) holds is this term can be neglected because it is sufficiently small compared with unity. The maximum damping ratio of buildings are about 5%. Therefore, \(\pi h^2/2=0.004\) (0.4%) and is sufficiently small. On the other hand, maximum damping ratio becomes about 30% or more in the soil, the third term yields \(\pi h^2/2=0.14\) or 14%, which may not be neglected.

Let's start the discussion from stress-strain relationships because what we are dealing with is soil. Frequently used stress-strain models for viscous material that uses a spring and a dashpot are Maxwell and Voigt models as shown in the following figures.

We use Voigt model where \(C\) denotes viscous coefficient and \(G\) denotes spring constant. The stress-strain relationships of this model yields Eq. (9)

\(\tau=G\gamma + C\dot{\gamma}\)(9)

Steady state condition under harmonic vibration is assumed in SHAKE. Then we define strain as

\(\gamma=\gamma_0\cos{\omega t}\) or \(\cos{\omega t}=\gamma/\gamma_0\)(10)

By substituting this equation into Eq. (9) and by using the relation \(\sin{\omega t}=\pm\sqrt{1-\cos{^2\omega t}}\), we can obtain stress-strain relationships such as

\(\tau=G\gamma_0\cos{\omega t}-C\gamma_0\omega\sin{\omega t}=G\gamma\pm C\omega\sqrt{\gamma_0^2-\gamma^2}\)(11)

Since this equation includes \(\omega\) in the right hand side, it is clear that stress-strain relationships depend on frequency. However, soil is modeled to frequency independent material under the important frequencies in the seismic response analysis (usually several Hz.), therefore Eq. (11) is inconvenient. Then a trick in the mathematic expression is employed; the viscous coefficient is frequency dependent. A new parameter \(\beta\) as defined in the following is introduced.

\(\omega C=2\beta G\) or \(C=2\beta G/\omega\)(12)

A frequency independent stress-strain relationships is obtained by substituting Eq. (11) as

\(\tau=G\gamma\pm 2\beta G\sqrt{\gamma_0^2-\gamma^2}\)(13)

The first term of this equation is a straight line and the second term is a circle. Then sum of the straight line and a circle becomes an elliptic as shown in the following figure.

It is noted that Eq. (12) and Eq. (4) are identical. This means that \(\beta\) used in Eq. (4) is not a critical damping ratio, but is derived from the assumption that stress-strain relationships of the Voigt model is frequency independent. As far as I know, there is no approach of this kind in our field, therefore, no relevant name for \(\beta\) in Eq. (12). Then we temporary named it a damping parameter. I hope that, in the future, a more relevant name will be found.

You may not be satisfied in the above explanation. So let's consider some more.

Let's compute damping constant from Eq. (12). At first, we computes absorbed energy \(\Delta W\) due to hysteretic behavior. The second term of Eq. (12) is a circle with radius \(\gamma_0\). Then area is easily calculated as \(\pi\gamma_0^2\). Then we obtains

\(\Delta W=2\beta G\pi\gamma_0^2\)(14)

On the other hand, strain energy (sometimes called elastic energy) \(W\) is evaluated as

\(W=G\gamma_0^2/2\)(15)

Therefore, damping constant \(h\) yields

\(h=\cfrac{1}{4\pi}\cfrac{\Delta W}{W}=\beta\)(16)

This shows that, as has been used, we can use damping ratio \(h\) obtained from the cyclic shear test. It is noted that approximation shown in Eq. (8) is not used here, this equation can be used at large damping constant.

Next, let's investigate stress-strain relationships in SHAKE in detail. The Voigt model is used in the above example. Use complex modulus to this kind or problem is first proposed by Sorokin³⁾. However, I cannot get the original paper. This paper is referred in many technical papers, but they were written in the old days and I cannot contact the author of the papers. Recently, I found a paper at which e-mail address is written⁴⁾. So I contacted him, but he response was that he saw the paper at the library in his old office and did not have at present. Since he showed the library, I wrote a letter, but there was no response. The method to use complex modulus is called as Sorokin's hypothesis⁵⁾ or Sorokin's damping⁶⁾.

Stress-strain relationships are derived based on the Sorokin model. It looks convenient to distinguish complex and real numbers. Then, a overbar is put on a complex number such as \(\overline{\gamma}\). Steady state under harmonic loading is considered same as in the preceding. Strain is given as

\(\overline{\gamma}=\gamma_0 e^{i\omega t}\)(17)

Stress strain relationships in a complex form is expressed as

\(\overline{\tau}=G\overline{\gamma}+C\dot{\overline{\gamma}}\)(18)

Substitution of Eq. (17) into Eq. (18), and using Eq. (12), we obtains

\(\overline{\tau}=G(1+2i\beta)\gamma_0e^{i\omega t}=\overline{G_S^*}\overline{\gamma}\)(19)

Here,

\(\overline{G_S^*}=G(1+2i\beta)\)(20)

is the complex modulus used in SHAKE. Subscript S indicates Sorokin model. Hysteretic absorbed energy \(\Delta W\) by this Sorokin mode can be evaluated same as the procedure in the preceding and yields

\(\Delta W_S=2G\beta\pi\gamma_0^2\)(21)

The value is same with the one obtained in the laboratory test. In other words, the complex modulus in SHAKE agrees with laboratory test result in both (G\) and \(h\). It is, however, noted that, as can be seen in Eq. (19) or in the above figure, maximum stress is \(\sqrt{1+4\beta^2}G\gamma_0\) and is overestimated compared with laboratory test.

3. Fundament equation of SHAKE

SHAKE¹⁾ deals with one-dimensional seismic response problem of a ground. Let's denote displacement as \(\overline{u}=\overline{u}(z,t)\), \(\rho\) as mass, and \(C\) as damping coefficient, according to the SHAME manual, equation of motion is written as

\(\rho\cfrac{\partial^2\overline{u}}{\partial t^2}=G\cfrac{\partial^2\overline{u}}{\partial z^2}+C\cfrac{\partial^2\overline{u}}{\partial t\partial z^2}\)(22)

where \(z\) denotes coordinate axis toward downward. Displacement is separated into time and space as

\(\overline{u}(z,t)=\overline{U}(z)e^{i\omega t}\)(23)

which result in

\((G+i\omega C)\cfrac{\partial^2\overline{U}}{\partial z^2}=-\rho\omega^2\overline{U}(z)\)(24)

It is noted that minus sign in the right hand side is not exist in the SHAKE manual. This error occurs because mathematic equation that square of imaginary (unit \(i\)) becomes -1 is not considered. Fortunately, they made the same error at another place, which cancel this error; resulting governing equation is correct.

Here, by setting complex modulus \(\overline{G}_S^*\) as Eq. (20) by using Eq. (12), we can obtain frequency independent stress-strain relationships.

Next, Eq. (22) can be obtained by the following procedure. Stress-strain relationships is already shown in Eq. (9). By using Eq. (12) and by substituting strain-displacement relationships \(\gamma=\partial u/\partial z\), we obtains

\(\tau=G\cfrac{\partial t}{\partial z}+\cfrac{2\beta G}{\omega}\cfrac{\partial^2u}{\partial t\partial z}\)(25)

On the other hand, equilibrium equation of infinitesimally small soil element in horizontal direction is

\(\cfrac{\partial u}{\partial z}=\rho\cfrac{\partial^2u}{\partial t^2}\)(26)

Substitution of Eq. (25) after partially differentiating by \(z\) into Eq. (26), we can obtain Eq. (22). Since the Sorokin model is used in Eq. (25), it is clear that SHAKE uses a wave equation with Sorokin model.

It is noted that, at present, not only SHAKE but also other many computer programs based on complex modulus for seismic response analysis of ground do not use Sorokin model. Prof. Lysmer, one of the developer of SHAKE, proposed a new complex modulus⁷⁾ and he suggest to change complex modulus in SHAKE by his proposal by a letter to peoples who bought SHAKE from the University of California. However, this model looks to have problems as mechanical model for the seismic response analysis of ground and damping ratio is evaluated smaller than that in laboratory test or input value. This will be discussed in Lysmer's proposal and its discussion issue.