AEiスクリプトで埋込SRモデルの解析を連続実行する

動的サブストラクチャ解析の一つに、建屋の一部(または全部)が表層地盤に埋め込まれた構造を、地盤の応答解析の結果を引き継ぎ、建屋の応答解析を実施するという方法があります。

この場合、建屋側面地盤の非線形性を等価線形解析で求め、この物性で求めた建屋側面のインピーダンスから地盤剛性と地盤減衰を算出して、建屋に接続して解析を行うことが一般的です。

以下は、その一連の解析を、入力地震動を変更しながら、また地盤物性と建屋減衰を変更してAEiスクリプトで連続実行するデモです。

入力地震動はエルセントロNSと日本建築センター模擬波(基盤波)、八戸NSの3種、地盤物性は標準、+σ、-σの3種、建屋減衰は5%と3%の2種で合計18ケースの解析を実施します。

なお、地盤の応答解析以外は自前のプログラムを使用しております。

地盤の応答解析は吉田望氏が作成されたDYNEQ336を使用させていただきました。
たいへん使いやすく、ドキュメントもしっかりいて素晴らしいものだと思います。

デモ動画のOSはMacOSですが、AEiスクリプトはWindowsでも、Linuxでも動作します。

参考として以下にaeiスクリプトを記載します。

print " aei スクリプトで、入力地震動と地盤物性を変化させ、埋込SRを連続解析して各ケースの解析結果のレポートを作成する"
pause

print "解析定義のExcelから解析ケースのシートを抽出する"
cmd{
	excelcsv -proc ex2csv -excel analysis.xlsx -sheet case -csv case.csv
	cat case.csv
	}
	pause

# 連続解析のループ
for i = 1 to 18
	L = i + 1

	print "解析ケースのCSVから入力波と地盤の情報を変数に代入する"
	csvfile = case.csv
	$NUM  = csv_r[L]_c[1]
	$SOIL = csv_r[L]_c[2]
	$WAVE = csv_r[L]_c[3]
	$WNUM = csv_r[L]_c[4]
	$WSTP = csv_r[L]_c[5]
	$BLDH = csv_r[L]_c[6]
	print "ケースNo   = ",$NUM
	print "地盤条件   = ",$SOIL
	print "入力地震動 = ",$WAVE,"( ",$WSTP,")"
	print "建屋減衰   = ",$BLDH
	pause

	print "各解析データのファイル名を定義する"
	$CASE = CASE<$NUM>
	$INPUT_WAV = <$WAVE>.wav
	$JBNWAV    = <$SOIL>_<$WAVE>
	$SHAKE_DAT = Shake_<$JBNWAV>.dat
	$SHAKE_OUT = Shake_<$JBNWAV>.out
	$SHAKE_WAV = Shake_<$JBNWAV>.wav
	$NOVAK_DAT = Novak_<$JBNWAV>.dat
	$NOVAK_OUT = Novak_<$JBNWAV>.out
	$ADMIT_DAT = Admit_<$SOIL>.dat
	$ADMIT_OUT = Admit_<$SOIL>.out
	$NAMAZ_DAT = <$CASE>_Namaz_<$JBNWAV>_h<$BLDH>.dat
	$NAMAZ_OUT = <$CASE>_Namaz_<$JBNWAV>_h<$BLDH>.out
	$NAMAZ_WAV = <$CASE>_Namaz_<$JBNWAV>_h<$BLDH>.wav
	print "Shake DATA  = ",$SHAKE_DAT
	print "Novak DATA  = ",$NOVAK_DAT
	print "Admit DATA  = ",$ADMIT_DAT
	print "Namaz DATA  = ",$NAMAZ_DAT
	pause

	print "解析定義のExcelから地盤情報を抽出して地盤物性を変数に代入する"
	cmd{
		excelcsv -proc ex2csv -excel analysis.xlsx -sheet <$SOIL> -csv JIBAN.csv
		cat JIBAN.csv
		}
	csvfile = JIBAN.csv
	for j = 1 to 7
		M = j + 1
		$G[j]  = csv_r[M]_c[2]
		$VS[j] = csv_r[M]_c[3]
		$RO[j] = csv_r[M]_c[4]
		$PO[j] = csv_r[M]_c[5]
		$T[j]  = csv_r[M]_c[6]
		$SH[j] = csv_r[M]_c[7]
	next j
	pause

	print "解析定義のExcelから入力波のCSVを抽出してShakeが読み込める書式に変換する"
	cmd{
		excelcsv -proc ex2csv -excel analysis.xlsx -sheet WAVE -csv WAVE.csv
		csv2wav.sh -c WAVE.csv -w <$INPUT_WAV> -r <$WNUM> -s <$WSTP>
		}
	pause

	print "Shakeデータを作成して解析を実行する"
	makedata shake.template $SHAKE_DAT
	print "Shakeのプログラムは 吉田望 氏 が作成された DYNEQ336 を使わせていただいています。"
	cmd{ sub.sh DYNEQ336 $SHAKE_DAT $SHAKE_OUT }
	pause

	print "アドミッタンスばねデータを作成して解析を実行する"
	makedata admit.template $ADMIT_DAT
	cmd{ sub.sh admit $ADMIT_DAT $ADMIT_OUT }
	pause
	
	print "Shake解析結果からNovakばね計算に必要な情報を取得する"
	cmd{
		shake2csv.sh -L $SHAKE_OUT -c SHAKE.csv
		cat SHAKE.csv
		}
	csvfile = SHAKE.csv
	for j = 1 to 7
		$SKG[j]  = csv_r[j]_c[1]
		$SKVS[j] = csv_r[j]_c[2]
	next j
	pause

	print "novakばねデータを作成して解析を実行する"
	makedata novak.template $NOVAK_DAT
	cmd{ sub.sh novak $NOVAK_DAT $NOVAK_OUT }
	pause
	
	print "アドミッタンスばねのKとCを取り出す(Cの計算は暫定的に1.0Hz)"
	cmd{
		admit2csv.sh -L $ADMIT_OUT -c ADMIT.csv -f 1.0
		cat ADMIT.csv
		}
	pause

	print "novakばねのKとCを取り出す(Cの計算は暫定的に1.0Hz)"
	cmd{
		novak2csv.sh -L $NOVAK_OUT -c NOVAK.csv -f 1.0
		cat NOVAK.csv
		}
	pause

	print "アドミッタンスばねとnovakばねの値を変数に代入する(Cの計算は暫定的に1.0Hz)"
	csvfile = ADMIT.csv
	$ADMK[1] = csv_r[1]_c[1]
	$ADMK[2] = csv_r[2]_c[1]
	$ADMC[1] = csv_r[1]_c[3]
	$ADMC[2] = csv_r[2]_c[3]
	print "SWAY K = ",$ADMK[1]," ROCK K = ",$ADMK[2]
	print "SWAY C = ",$ADMC[1]," ROCK C = ",$ADMC[2]
	csvfile = NOVAK.csv
	for j = 1 to 6
		$NVKK[j] = csv_r[j]_c[1]
		$NVKC[j] = csv_r[j]_c[3]
		print "NOVAK K= ",$NVKK[j]," NOVAK C= ",$NVKC[j]
	next j
	pause

	print "固有値解析用データ作成と実行(一度、1.0Hzの減衰で応答解析を実施する)"
	makedata namazu.template $NAMAZ_DAT
	cmd{ sub.sh namazu $NAMAZ_DAT $NAMAZ_OUT }
	pause

	print "解析モデルの地盤-建屋連成1次の振動数を取得する"
	cmd{
		namazu2csv.sh -L $NAMAZ_OUT -c NAMAZU.csv
		cat NAMAZU.csv
		}
	pause

	print "連成1次の減衰を取得する"
	csvfile = NAMAZU.csv
	$FREQ = csv_r[1]_c[1]
	cmd{
		admit2csv.sh -L $ADMIT_OUT -c ADMIT.csv -f $FREQ
		novak2csv.sh -L $NOVAK_OUT -c NOVAK.csv -f $FREQ
		}
	pause

	print "アドミッタンスばねとnovakばねの値を変数に代入する"
	csvfile = ADMIT.csv
	$ADMK[1] = csv_r[1]_c[1]
	$ADMK[2] = csv_r[2]_c[1]
	$ADMC[1] = csv_r[1]_c[3]
	$ADMC[2] = csv_r[2]_c[3]
	print "SWAY K = ",$ADMK[1]," ROCK K = ",$ADMK[2]
	print "SWAY C = ",$ADMC[1]," ROCK C = ",$ADMC[2]
	csvfile = NOVAK.csv
	for j = 1 to 6
		$NVKK[j] = csv_r[j]_c[1]
		$NVKC[j] = csv_r[j]_c[3]
		print "NOVAK K= ",$NVKK[j]," NOVAK C= ",$NVKC[j]
	next j
	pause

	print "応答解析を実施する"
	makedata namazu.template $NAMAZ_DAT
	cmd{ sub.sh namazu $NAMAZ_DAT $NAMAZ_OUT }
	pause
	
	print "応答スペクトル計算"
	pause
	cmd{
		wav2csv -in $INPUT_WAV -csv <$INPUT_WAV>.csv -fmt '(8f10.0)' -step $WSTP -dt 0.01
		wav2csv -in $NAMAZ_WAV -csv <$NAMAZ_WAV>.csv -fmt '(8f10.0)' -step $WSTP -dt 0.01
		duhamel -c <$INPUT_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 1 -b 300
		duhamel -c <$NAMAZ_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 1 -b 300
		duhamel -c <$NAMAZ_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 2 -b 300
		duhamel -c <$NAMAZ_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 3 -b 300
		duhamel -c <$NAMAZ_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 4 -b 300
		duhamel -c <$NAMAZ_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 5 -b 300
		}
	pause

	print "一括後処理(関連するCSVファイルを全てTemplateのExcelに貼り付ける)"
	cmd{
		sr_post.sh <$CASE> <$WAVE> <$SOIL> <$BLDH> <$INPUT_WAV> <$CASE>_<$WAVE>_<$SOIL>_<$BLDH>.xlsx
		open <$CASE>_<$WAVE>_<$SOIL>_<$BLDH>.pdf
		}
	pause

next i

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

    # 第2種変形ベッセル関数
    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)