試しに、FORTRANで一質点の時刻歴応答解析を行うコードを書いてみようと思います

運動方程式

古典力学分野で解析系の研究をしていた人にとっては、馴染みの深い式かと思います。

m \frac{d^2 x}{dt^2} + c \frac{dx}{dt} + kx = -m\frac{d^2 x_0}{dt^2}

という形が一般的でしょうか。
(mは質量、cは粘性減衰定数、kはバネ剛性、右辺は地動加速度による外力項)
これを、

一質点モデル、バネは弾性
入力地震波は阪神大震災の観測波(JMA Kobe NS)
ニューマークのβ法は平均加速度法(β=1/4)

で、解いてみようかと思います。
また、減衰力を合わせた応答スペクトルを用いて、解析の妥当性を確認します。

FORTRANでコードを書く

早速、コードを書いてみようと思います。
まずは文字の定義から。

program equationOfMotion
implicit none
!---文字の定義------
integer :: i, n, step
real(8) :: ydd_EqNS(20000), ydd_EqEW(20000), ydd_EqUD(20000) !地震波(gal)
real(8) :: y(20000), yd(20000), ydd(20000), time(20000) !変位(m)、速度(m/s)、加速度(m/s/s)、時間(s)
real(8) :: beta, dt, omega, T, h, m, k, c, pi, delta_p,epsilon, delta_y

implicit none忘れずに書きましょうね。

平均加速度法なのでβ=1/4、減衰は5%で剛性比例型減衰としました。
ちなみに、一次固有周期は1.0sになるよう、質量とバネ剛性を決めました。

beta = 0.25                       !β法(平均加速度法なので1/4)
h = 0.05                          !減衰定数(5%の減衰を与える)
c = 2.0 * h * k / omega           !粘性係数(剛性比例型)

解析部分は以下のように記述しています。次ステップの変位を仮定して一度運動方程式を解き、その後、誤差が一定の値以下になるまで、収斂計算をしていきます。収斂計算の誤差はε以下として、

 \epsilon = 1.0 \times 10.0^{-5.0}

と設定しました。

!---次ステップの値を求める----------
do n = 1, step, 1
    y(n+1)   = y(n)                                                                                          !変位を仮定
    ydd(n+1) = (y(n+1) - y(n)) / (beta * dt**2) - yd(n) / (beta * dt) - ydd(n) * (1.0 / (2.0 * beta) - 1.0)  !加速度を求める
    yd(n+1)  = yd(n) + 0.50 * (ydd(n) + ydd(n+1)) * dt                                                       !速度を求める

    !---収れん判定（収れんしていない場合にはループの計算を続ける）---
    do
        delta_p = (-1.0) * m * ydd_EqNS(n+1) * 0.01 - (m*ydd(n+1) + c*yd(n+1) + k*y(n+1)) !不釣合力を計算
        
        if (abs(delta_p).lt.epsilon) exit
            delta_y = delta_p / (k + m / (beta * dt**2) + c /(2.0*beta*dt))
            y(n+1) = y(n+1) + delta_y                                                                               !変位修正
            ydd(n+1) = (y(n+1) - y(n)) / (beta * dt**2) - yd(n) / (beta * dt) - ydd(n) * (1.0 / (2.0 * beta) - 1.0) !加速度を求める
            yd(n+1)  = yd(n) + 0.50 * (ydd(n) + ydd(n+1)) * dt                                                      !速度を求める
    end do

解析結果

今回作ったモデルでの解析結果を以下に記します。

FORTRANでの解析結果

青線が入力地震動、赤線が応答値です。

最大応答加速度(絶対値):14.967(m/s/s)
となっており、かなり揺られていることがわかります。

加速度応答スペクトル

今回は一次固有周期1.0sなので、グラフからは 約15(m/s/s) と読み取れます。
これより解析モデルの妥当性が確認できました。

久しぶりにFORTRANでコードを書いたので、割と手こずりましたが、解析できるまでコーディングできたので良かったです。