Fortranでニューマークのβ法

【管理人おすすめ！】セットで3割もお得！大好評の用語集と図解集のセット⇒　建築構造がわかる基礎用語集＆図解集セット(※既に26人にお申込みいただきました！）

自然現象等の不規則な外乱を受ける場合、微分方程式より解を求めるということは未知数が多くてできません。よって、このような場合、近似解を得ることが工学として一般的です。地震動も不規則な外乱の一つであり、このような外乱が作用した建物の応答（どのように揺れるのか？速度は？加速度はどうなっているのか？）は「ニューマークのβ法」と呼ばれている数値解析より求めています。ニューマークのβ法を理解していない人は当サイトの「ニューマークのβ法について」から勉強してください。

ニューマークのβ法の式

まず、ニューマークのβ法で得られる変位、速度、加速度の式は、

となります。

早速、Fortranで書いてみる。

さて、ダウンロードしたエディタ（僕はTerapad）を開いてプログラムを書いてみます。まだコンパイラしていません。とりあえず、「Terapad」を利用している方は

「名前をつけて保存」→保存形式を「全てのファイル」として→「○○.f90」と保存してください。

そうすると、F90ファイルが生成されます。書いたプログラムはこんな感じ。

下記にニューマークのβ法を計算するソースコードを書いています。コピペしてそのまま使えるものだとは思いますが、式の意味を理解することや自分でプログラムをつくることを必ずやってみないと勉強にはなりませんので、気をつけて下さい。

newmarkb.f90のソースコード

!一自由度系の時刻歴応答解析()

program newmark

implicit none

real gg(10000) ! 入力加速度

real acc(10000) ! 応答加速度

real vel(10000) ! 応答速度

real dis(10000) ! 応答変位

integer i

real mm,dt,tt,Pi,omg,kk,hh,cc,cm,km,mol1,mol2,B

open(1, file='step-load.txt',status='old')

print *, "質量(N/(mm/s/s))を入力してください"

read *, mm

print *, "刻み時間(s)を入力してください"

read *, dt

print *, "固有周期(s)を入力してください"

read *, tt

print *, "減衰定数(%)を入力してください"

read *, hh

print *, "B(0,1/4,1/6,1/8)を入力してください"

read *, B

Pi = 3.14159265358979

omg = (2 * Pi) / tt

kk = mm * omg * omg

cc = hh * (2 * mm * omg)

!'入力加速度の読み込み'

! zisin-data. から gg() に読み込む

do i = 1,2687

read(1,*) gg(i)

end do

!'Nemark's β'

cm = cc / mm

km = kk / mm

do i = 1,10000,1

if(i==1) then

acc(i) = -gg(i)/mm

vel(i) = 0

dis(i) = 0

Else

mol1 = -((gg(i) + cc(i2) * (vel(i-1)) + 0.5 * acc(i-1) * dt)+ kk(i2) * (dis(i-1) + vel(i-1) * dt + (0.5-B) * acc(i-1) * dt * dt))

mol2 = (mm + 0.5 * cc(i2) * dt + B * kk(i2) * dt * dt)

acc(i) = mol1 / mol2

vel(i) = vel(i - 1) + 0.5 * (acc(i - 1) + acc(i)) * dt

dis(i) = dis(i - 1) + vel(i - 1) * dt + (0.5 - B) * acc(i - 1) * dt * dt + B * acc(i) * dt * dt

end if

enddo

open(2, file='acc-kekka.txt',status='replace')

open(3, file='vel-kekka.txt',status='replace')

open(4, file='dis-kekka.txt',status='replace')

open(5, file='dt-kekka.txt',status='replace')

open(6, file='cm-kekka.txt',status='replace')

open(7, file='km-kekka.txt',status='replace')

do i = 1,10000,1

Write(2,*) ' acc+gg =',acc(i) + gg(i)

Write(3,*) ' vel =',vel(i)

Write(4,*) ' dis =',dis(i)

Write(5,*) ' dt =',dt * (i - 1)

Write(6,*) ' cm =',cm

Write(7,*) ' km =',km

end do

close(unit=1)

close(unit=2)

close(unit=3)

close(unit=4)

close(unit=5)

close(unit=6)

close(unit=7)

stop

end program newmark

newmarkb.f90の説明

１つ目のポイントは「open(1, file='step-load.txt',status='old')」で、入力データの読み込みを行っていることです。同じフォルダ内にファイルが無ければ実行エラーが起きるので気をつけて下さい。ここでは、ステップ荷重を表す入力データを読み込みました。

ニューマークのβ法を勉強していくとわかりますが、初期値を設定する必要があります。つまり、変数i=1のときにおける変位や速度、加速度の値とi=2からの状態で条件分岐させる必要があるのです。つまり、「if(i==1) then」とそれ以外の「else」で分岐させています。

あとは、目立って特別なテクニックは使っていないと思います。意外と簡単ですよね。

newmark2.f90をコンパイラしてみる

さて、コンパイラしてみましょう。デスクトップにある「Fortran」フォルダを開いてください。開いたらコマンドプロンプトのショートカットをクリックします。コマンドプロンプトの画面に「g95.exe ○○○.f90」と入力してEnterキ―を押します。フォルダ内に「a.exe」という実行ファイルが生成されます。コマンドプロンプトに「a.exe」と入力すると、プログラムが回り値の設定を行った後、フォルダに「acc-kekka.txt」を初めとする計算結果が生成されたなら成功です。次は→「計算結果の確認」を勉強しましょう。