投稿日:

pythonで作るconeモデル基礎底面ばね

こちらの文献を参考に作成してみました。
http://library.jsce.or.jp/jsce/open/00550/2010/16-01-0030.pdf

元のWolf氏の文献を見つけらず「z0＝近似円錐地盤の頂点高さ」にどのような数値を設定すべきか不明です。

import math
import cmath
import scipy.special

G = 1.8E+07
VS = 100.0
poisson = 0.33
R = 6.0
Z0 = 15.0

for i in range(0,80):

    freq = float(i) * 0.1
    omega = 2.0 * math.pi * freq

    a0 = omega * R / VS

    K = 8.0 * G * R / ( 2.0 - poisson )
    ca0 = math.pi * ( 2.0 -poisson ) / 8.0
    a0cao = a0 * ca0
    Zsway = K * complex(1.0, a0cao) / ( G * R )

    mu = 0.3 * math.pi * ( poisson - 1.0 / 3.0 )
    Kt = 8.0 * G * R**3.0 / ( 3.0 * ( 1.0 - poisson ) )
    a02 = a0**2.0
    kta0 = 1.0 - 1.0 / 3.0 * mu / math.pi * Z0 / R * a02 - 1.0 / 3.0 * a02 / ( ( 2.0 * R / Z0 )**2.0 + a02 )
    cta0 = Z0 / ( 6.0 * R) * a02 / ( ( 2.0 * R / Z0 )**2.0 + a02 )
    Zrock = Kt * complex( kta0, a0 * cta0 ) / ( G * R**3.0  )

    print(a0, ',', Zsway.real, ',', Zsway.imag, ',', Zrock.real, ',', Zrock.imag)
投稿日:

pythonで作るNOVAKばね計算プログラム

novak impedance
novak impedance

建屋側面地盤のインピーダンスを簡易計算する手法としてNOVAKのばねがよく知られています。

こちらのPDFにその理論式が掲載されています。

これをもとにpythonで作成したプログラムが以下です。

import math
import cmath
import scipy.special

G = 1.0E+09
VS = 700.0
poisson = 0.3
R = 10.0
HJ = 1.0
h = 0.00

VP = VS * math.sqrt( ( 2.0 * ( 1.0 - poisson ) ) / ( 1.0 - 2.0 * poisson ) )

for i in range(1,2000):

    freq = float(i) * 0.01
    omega = 2.0 * math.pi * freq
    a0 = omega * R / VS
    b0 = omega * R / VP

    Z = cmath.sqrt( 1.0 + 2.0 * h * complex( 0.0, 1.0 ))
    a0Z = a0 * complex( 0.0, 1.0 ) / Z
    b0Z = b0 * complex( 0.0, 1.0 ) / Z
    # B = math.sqrt( ( 2.0 * ( 1.0 - poisson ) ) / ( 1.0 - 2.0 * poisson ) )
    # b0Z = a0Z / B

    # 第２種変形ベッセル関数
    k0a0 = scipy.special.kv(0, a0Z)
    k0b0 = scipy.special.kv(0, b0Z)
    k1a0 = scipy.special.kv(1, a0Z)
    k1b0 = scipy.special.kv(1, b0Z)

    p1 =  - math.pi * HJ * G * a0**2.0
    p2 = 4.0 * k1b0 * k1a0 + a0Z * k1b0 * k0a0 + b0Z * k0b0 * k1a0
    p3 = b0Z * k0b0 * k1a0 + a0Z * k1b0 * k0a0 + b0Z * a0Z * k0b0 * k0a0
    kuj = p1 * p2 / p3

    print(freq, ',', kuj.real, ',', kuj.imag)