前回、 SCILABで周波数応答関数を計算する では、大崎順彦（1994）「新・地震動のスペクトル解析入門」を参考にしてSCILABにより地盤の加速度/加速度応答を求めました。
実は大崎（1994）のFESPでは各層中点のひずみの大きさも計算できるようになっています。そこで、同様にSCILABで拵えてみましたので御紹介します。加速度応答の場合ほどに適当な検証方法が実は見当たらないのですが、一応本家のFESPと同じ結果が得られることは確認しています。
#!/bin/bash temp=`mktemp ./XXXX` cat エルセントロ.txt |\ sharp_cancel |\ awk 'NR<=800{print $1,-$2}' > $temp scilab -nw << _SCI_ | sed -e "/^ *$/d" function [z]=strain_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 を返す // 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 ea=exp(P(obj)*omega*thick(obj)/2.); eb=1/ea; zobj=-P(obj)/omega*(A(obj)*ea-B(obj)*eb)*0.01; if b_cond=='inner' then zref=A(ref)+B(ref); // 内部境界面 elseif b_cond=='expose' then zref=2*A(ref); // 解放基盤面 else error("境界条件が不適当です"); exit; end z=zobj/zref; endfunction nlayer=6; //層の数 thick=[5.0 17.0 13.0 31.0 4.0]; // 層厚(m) gamma_w=[1.98 1.95 1.94 1.94 1.94 1.83]; // 単位体積重量(tf/m3) G=[2900. 8000. 13400. 17800 24200 22800]; // せん断弾性係数(tf/m2) h_alpha=[2. 2. 2. 2. 2. 2.]; // マクスウェル型散乱減衰（地盤の減衰） h_beta=[0.020 0.032 0.024 0.014 0.010 0.010]; // 履歴減衰型材料減衰（土の減衰） ACC=9.8; //重力加速度(m/sec2) ref=1; // 基準層 //エルセントロ波の読み込み fd=mopen('$temp','r'); [n,t,a0]=mfscanf(-1,fd,'%f %f'); mclose(fd); ns0=size(a0,'r'); // ns0=800 fs=100; // サンプリング周波数 // 後続のゼロ a=[a0;zeros(1024-ns0,1)]; ns=size(a,'r'); // ns=1024(2の10乗) nspe=ns/2+1; // スペクトルの要素数 // ここまででエルセントロ波の読み込み freq0=[0:fs/ns:fs]; // 周波数軸の設定 x=fft(a,-1); // FFTの実行 for obj=1:nlayer-1 for i=1:nspe // DCからナイキストまで if i==1 then z=0; else f=freq0(i); // 周波数 z=strain_resp(f,'inner',obj,ref); end // if文の終わり st(i)=z*x(i); // strain_respの結果と各周波数の複素フーリェ級数を掛け合わせる end // 周波数ループの終わり for i=nspe-1:-1:2 // 共役複素数の作成 st(ns-i+2)=conj(st(i)); end // 共役作成の終わり y=fft(st,+1); // 時間領域に戻す mprintf('対象層:%d 最大ひずみ%.4f(%%)\n',obj,max(abs(real(y)))*100); //各層の最大ひずみを表示する 単位を(%)にするためには100を掛ける end // 層ループの終わり _SCI_ rm $temp
