SCILABで周波数応答関数を計算する

地盤の増幅スペクトルを計算します。大崎順彦(1994)のFESPを参考にしました。ですので、大地震時にせん断ひずみが大きくなって剛性率が低下するという「非線形性」にまでは対応できません。

一応、オリジナルのFESPおよび吉田望らのDYNEQ320(弾性材料使用)との比較・検証はやってみました。DYNEQは散乱減衰の設定の仕方が少し異なっているので、近似値を使用しました。

 

f:id:S_E_Hyphen:20180301093901j:plain

しばらくの間、このスクリプトを使って色々試してみようと思っています。

 

#!/bin/bash
temp=`mktemp ./XXXX`
scilab -nw << _SCI_ >$temp

function [z]=acc_resp(f,b_cond,obj,ref)
// 大崎(1994)§9のFESP.forをSCILAB用に
// 書き直したもの
// 加速度-加速度応答を求める
// [引数]
// f : 周波数(実数)
// b_cond : 境界条件(文字列) 'inner' or 'expose'
// obj : 対象層(整数)
// ref : 基準層(整数)
// [その他]
// nlayer(I),thick(ARRAY),gamma_w(ARRAY),G(ARRAY),
// h_alpha(ARRAY),h_beta(ARRAY)はグローバル変数
// として引き渡す
// [返値]
// 周波数応答関数 Zobj/Zref を返す
//                            2018年2月28日 作成
// 
omega=2*%pi*f; // 円振動数
h=h_alpha./omega+h_beta; // 減衰定数の算出
gb=G.*(1+2*h*%i); // 複素剛性の算出
P=sqrt(gamma_w./gb/ACC)*%i; // 伝播定数の算出
for layer=1:nlayer-1
 R(layer)=gb(layer)/gb(layer+1)*P(layer)/P(layer+1); // インピーダンス比の算出
end

A(1)=1.;
B(1)=1.;
for layer=1:nlayer-1
 ec=exp(P(layer)*omega*thick(layer));
 ea=A(layer)*ec;
 eb=B(layer)/ec;
 A(layer+1)=((1.+R(layer))*ea+(1.-R(layer))*eb)/2.;
 B(layer+1)=((1.-R(layer))*ea+(1.+R(layer))*eb)/2.;
end

if b_cond=='inner' then z=(A(obj)+B(obj))/(A(ref)+B(ref)); // 内部境界面
elseif b_cond=='expose' then z=(A(obj)+B(obj))/(2*A(ref)); // 解放基盤面
else error("境界条件が不適当です"); exit; end

endfunction

nlayer=4; //層の数
thick=[3.8 3.2 3.9]; // 層厚(m)
gamma_w=[1.50 1.67 1.85 1.95]; // 単位体積重量(tf/m3)
G=[1200. 2900. 5700. 50000.]; // せん断弾性係数(tf/m2)
h_alpha=[2. 2. 2. 2.]; // マクスウェル型散乱減衰(地盤の減衰)
h_beta=[0.02 0.02 0.02 0.02]; // 履歴減衰型材料減衰(土の減衰)
ACC=9.8; //重力加速度(m/sec2)
obj=1; // 対象層
ref=4; // 基準層

if size(thick,'c')+1 <> nlayer then error("層厚の数が層数と一致しません"); end
if size(gamma_w,'c') <> nlayer then error("単位体積重量の数が層数と一致しません"); end
if size(G,'c') <> nlayer then error("せん断弾性係数の数が層数と一致しません"); end
if size(h_alpha,'c') <> nlayer then error("散乱減衰の数が層数と一致しません"); end
if size(h_beta,'c') <> nlayer then error("材料減衰の数が層数と一致しません"); end
[str,ier]=lasterror();
if ier==10000 then exit; end

df=50/128;
freq0=[0:df:50];
for i=1:size(freq0,'c')
 if i==1 then z=1;
 else
  f=freq0(i); // 周波数
  z=acc_resp(f,'inner',obj,ref);
 end
 mfprintf(6,'%f %f %f\n',freq0(i),real(z),imag(z));
end

_SCI_

cat fesp_out.txt |\
 awk '{print $1,sqrt($2*$2+$3*$3)}' |\
 gmt psxy -JX15c/10c -R0/21/0/8 -BnWeS+t"増幅スペクトルの比較" \
    -Bxa5f1+l"周波数(Hz)" -Bya2f1+l"増幅率" -W2p -K > resp.ps

cat dyneq_out.txt |\
 sed -e "1,2d" |\
 gmt psxy -J -R -B -W1p,red -O -K >> resp.ps

cat $temp |\
 sed -e "/^ *$/d" |\
 awk '{print $1,sqrt($2*$2+$3*$3)}' |\
 gmt psxy -J -R -B -S+0.4c -W1p,blue -O -K >> resp.ps
rm $temp

gmt gmtset FONT_TITLE 16p,GothicBBB-Medium-UniJIS-UTF8-H
gmt gmtset FONT_LABEL 12p,GothicBBB-Medium-UniJIS-UTF8-H
gmt gmtset FONT_ANNOT_PRIMARY 10p,GothicBBB-Medium-UniJIS-UTF8-H

cat << + |\
gmt pslegend -R -J -Dx11c/9c+w3c/1.6c+jTL -F+p+gwhite -O >> resp.ps
N 1
S 0.2c - 0.4c 0/0/0 2p,black 0.8c 大崎FESP
G 0.1c
S 0.2c - 0.4c 255/0/0 1p,red 0.8c DYNEQ320
G 0.1c
S 0.2c + 0.4c 0/0/255 1p,blue 0.8c SCILAB
+