Deprecated: The each() function is deprecated. This message will be suppressed on further calls in /home/zhenxiangba/zhenxiangba.com/public_html/phproxy-improved-master/index.php on line 456
c hippon.f を改造中．
c
c
c 集成材梁の上下面に軸方向に溝を掘って鋼リブを挿入して補強した剛性梁に
c 温度（と含水による膨張）による荷重を与えて（右下だけを）tyokuon.f に
c 1/4解析させるための入力データを作るプログラム。
c 計算例は
c 後藤文彦，高橋明洋，薄木征三，佐々木貴信：
c 熱膨張・乾燥収縮を受ける鋼板挿入集成 材梁の有限要素解析，
c 第2回木橋技術に関するシンポジウム論文報告集， pp. 97-100, (2003).
c とかだけど、鋼板・集成材境界面の要素の変形が大きすぎるので、
c あまりいい計算例ではない。
c
c O___________x(上下対称面）
c LLLLLL
c LLLLLL
c LLLLLL
c #}LLLL Lが集成材要素（場所によってサイズは3種類）
c #}LLLL #が鋼材要素（サイズは1種類）
c | }が接着剤要素（サイズは1種類）
c y（左右対称面）
c
c 直方体要素の有限要素プログラム tyokuon.f で
c z軸に沿って横たわる片持ち梁の自由端に集中荷重（y方向）を与える問題を
c yz面を対称面とする半解析で解かせるための入力データを作るプログラム。
c このプログラムを実行させると、hippon.d というデータができあがるから
c tyokuonhoge )
c と実行させればよい。
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 03/10/9版
c プログラムを作った人：後藤文彦
c このプログラムは無保証です。
c このプログラムは改造／再配布して構いません。
c 但し、私的利用に留まらない場合は、もとのプログラムが
c http://www.str.ce.akita-u.ac.jp/~gotou/programoj/
c から取ってきたものと分かるように明示して下さい。
c 改造／再配布後も上記の条件を保持して下さい。
c バグの報告は、
c http://www.str.ce.akita-u.ac.jp/~gotou/jpg/meeru.jpg
c にお願いします。
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
program hippou
implicit real*8 (a-h, o-z)
dimension ijkab(9400),a(3,3),b(3,3),c(3,3),
& e66a(6,6),e66b(6,6),e66c(6,6),
& f(9400*3),fab(6,24),n8(8),
& bt1(24,6),bt2(24,6),bt3(24,6),bt4(24,6),bt5(24,6),bt6(24,6)
c
c スカイラインの1次元配列のサイズは、同じ自由度数でも要素分割
c に依存するので、自由度数から一意には決められない。
c mage.f の要素分割のやり方だと、軸長（z軸）の分割数がxy断面
c の分割数に比べて多いとスパース部分が増える。
c 仮にスパース部分が皆無だったとすると
c およそ自由度×自由度／２のサイズになるが、こんなサイズは想定外。
c 取り敢えず、9400を最大節点数、その３倍を最大自由度数にしておく。
c 具体的な要素分割例としては、(縦+1)×(横+1)×(軸長+1)が、
c (6+1)*(12+1)*(99+1)=9191,(20+1)*(20+1)*(20+1)=9261 など。
c mage.f で幅（半解析だから幅の半分）×高さ×軸長が
c 1*(8〜20)*(50〜100)ぐらいの範囲のものを解く場合、有効数字３桁
c の収束解が得られる要素分割数は 6分割×10分割×100分割ぐらい。
c 図心が節点上に来るようにするには高さ方向の分割数を偶数にする。
c 幅（の半分）に対して桁高が10倍以上高い場合は、幅の分割数は
c 1 でもいいぐらいかも。
c
c 断面横方向の要素分割数：nx
c 断面縦方向の要素分割数：ny
c 軸方向の要素分割数：nz
write(*,*) '要素分割数を入力して下さい'
write(*,*) '（横,縦,軸長の順にコンマ区切りで）'
write(*,*) '例えば、20,4,2'
read(*,*) nx,ny,nz
c 接着剤要素の断面横方向の要素分割数：nxc
c 接着剤要素の断面縦方向の要素分割数：nyc
write(*,*) ' 接着剤要素断面の分割数を入力して下さい'
write(*,*) ' （横,縦の順にコンマ区切りで）'
c write(*,*) ' 例えば、',2,',',ny/2
write(*,*) '鋼板を上下にぶっとおす時の縦の分割数は',ny
read(*,*) nxc,nyc
c 鋼材要素の断面横方向の要素分割数：nxa
c 鋼材要素の断面縦方向の要素分割数：nya
nya=nyc
write(*,*) '鋼材要素の幅方向の分割数を入力して下さい'
write(*,*) ' 標準は 1'
read(*,*) nxa
c 荷重増分：df
c write(*,*) '荷重増分を入力して下さい'
c read(*,*) df
df=1.d0
c 荷重増分ステップ数：nstep
c write(*,*) '荷重増分ステップ数を入力して下さい'
c read(*,*) nstep
nstep=1
c 温度変化：dt
write(*,*) '温度変化を（℃で）入力して下さい'
read(*,*) dt
c 木材含水率変化：dh
write(*,*) '含水率変化を（％で）入力して下さい'
read(*,*) dh
c 梁断面の幅：yoko, 桁高：tate
yoko=100.d-3
tate=400.d-3
c 鋼板の断面の横の長さ:yokoa0, 縦の長さ:tatea0
c 接着剤層の横の長さ:yokoc0, 縦の長さ:tatec0
if( (nxa.eq.nx).and.(nya.eq.ny) ) then
c もし全部鋼材だったら
yokoa0=yoko
tatea0=tate
else
c 鋼材が全部ではないとき
yokoa0=5.d-3
tatea0=400.d-3
yokoc0=1.d-3
tatec0=tatea0
end if
c 鋼材要素の断面の横の長さ:yokoa, 縦の長さ:tatea
c 接着剤要素の断面の横の長さ:yokoc, 縦の長さ:tatec
if( (nxa.eq.0).and.(nya.eq.0) ) then
c もし木材だけだったら
yokoa0=0.d0
tatea0=0.d0
yokoa=0.d0
tatea=0.d0
yokoc0=0.d0
tatec0=0.d0
yokoc=0.d0
tatec=0.d0
else
c 鋼材があるときは
yokoa=yokoa0/nxa
tatea=tatea0/nya
yokoc=yokoc0/nxc
tatec=tatec0/nyc
end if
c 木材要素(分類番号4）の断面の横の長さ:yokob, 縦の長さ:tateb
if( (nxa.eq.nx).and.(nya.eq.ny) ) then
c もし全部鋼材だったら
yokob=0.d0
tateb=0.d0
yokoc=0.d0
tatec=0.d0
else
c 鋼材が全部ではないとき
if(nya.eq.ny) then
c もし鋼板が上から下までぶっとおっていたら
yokob=(yoko-yokoa0-yokoc0)/(nx-nxa-nxc)
tateb=0.d0
else
yokob=(yoko-yokoa0-yokoc0)/(nx-nxa-nxc)
tateb=(tate-tatea0)/(ny-nya)
end if
end if
c 鋼板が底面に貼ってある（左右にぶっとおっている）ようなものは、
c aを木材,bを鋼材に置き換えて90゜回転させて解く。
c
c 梁の軸長：ziku
ziku=10.d0
ziku1=ziku/float(nz)
write(*,*) 'yokoa=',yokoa
write(*,*) 'tatea=',tatea
write(*,*) 'yokob=',yokob
write(*,*) 'tateb=',tateb
write(*,*) 'yokoc=',yokoc
write(*,*) 'tatec=',tatec
write(*,*) 'ziku1=',ziku1
call bunrui(nx,ny,nz,nxa,nya,nxc,ijkab)
c
c ********** 境界条件 *************
c
c 拘束箇所の数：nkasyo
c 同じ節点を定義しなおした場合は、後から定義した方が有効
c 一箇所目で定義した場所を二箇所目で定義しなおすこともできる
c 固定端の面（z=0面）が1箇所
c 固定端面の上辺（x軸）が2箇所目
c 左右対称面（y=0面）が3〜3+nz箇所目
c 固定端面の左端（y軸）が3+nz〜3+nz+ny箇所目
c 上下対称面（x=0面）が3+nz+ny+1〜3+nz+ny+(ny+1)*nz箇所目
c 対称軸（z軸）が3+nz+ny+(ny+1)*nz+1〜3+nz+ny+(ny+1)*nz+nz箇所目
c 対称点（原点）が 3+nz+ny+(ny+1)*nz+nz+1箇所目
kkasyo=3+nz+ny+(ny+1)*nz+nz+1
c 1箇所目の拘束箇所における拘束開始節点番号:kkais1
kkais1=1
c 1箇所目の拘束箇所における拘束終了節点番号:kryou1
kryou1=(nx+1)*(ny+1)
c tyoku.fでは拘束開始点から拘束終了点まで節点番号順に全て
c 同じ拘束条件を与える。要素番号の順番は、最前面最上段左端の節点を
c 1として、右(x)方向に進み、右端まできたらy方向へ1段下の段の左端から
c 右へすすみ、真下の右端までいったら、z方向へ1面奥の面の左上から、
c また同じように1段づつ右へ進み、最下段の右端までいったら、もう一つ
c 奥の面へと進んでいき、一番奥の面の最下段の右端が最終要素番号
c 1箇所目の拘束箇所における拘束条件
c u,v,w 方向の変位をそれぞれ拘束するか開放するか(0/1)
c まず、固定端面のz方向変位のみをすべて拘束する
bu1=1.d0
bv1=1.d0
bw1=0.d0
c 固定端の上下対称線（上辺）のyz方向を拘束する
kkais2=1
kryou2=nx+1
bu2=1.d0
bv2=0.d0
bw2=0.d0
c 上下対称面、左右対称面の拘束箇所はいっぱいあるので、
c 下のデータ出力のところで、doループで直接 書き出す
c
c 因みに、固定端の全ての節点でxy方向の変位を拘束してしまうと、
c （つまり、固定端で材料がちぢめないような境界条件だと）
c w=LP/EAの解からは、ちょっとずれる（要素数を増やしても）。
c
c
c
c ******* 材料定数 *********
c 鋼材の材料定数
c ヤング率：exxa,eyya,ezza
exxa=200.d9
eyya=exxa
ezza=exxa
c ポアソン比：poixya,poixza,poiyza
poixya=1.d0/3.d0
poixza=poixya
poiyza=poixya
c せん断弾性係数：gxya,gxza,gyz
gxya=ezza/2.d0/(1.d0+poixya)
gxza=gxya
gyza=gxya
c
c aに柔性（コンプライアス）行列を入れる
a(1,1)=1.d0/exxa
a(1,2)=-poixya/exxa
a(1,3)=-poixza/exxa
a(2,1)=a(1,2)
a(2,2)=1.d0/eyya
a(2,3)=-poiyza/eyya
a(3,1)=a(1,3)
a(3,2)=a(2,3)
a(3,3)=1.d0/ezza
c
c 鋼材の線膨張率：aa
aa=12.d-6
c 鋼材の温度歪：epa
epa=aa*dt
c
c
c 木材の材料定数（今は確認のため変な値を入れてる）
c ヤング率：exxb,eyyb,ezzb
exxb=0.5d9
eyyb=1.d9
ezzb=10.d9
c ポアソン比：poixyb,poixzb,poiyzb
poixyb=0.2d0
poixzb=0.025d0
poiyzb=0.03d0
c せん断弾性係数：gxyb,gxzb,gyzb
gxyb=0.03d0
gxzb=0.6d0
gyzb=0.7d0
c bに柔性（コンプライアス）行列を入れる
b(1,1)=1.d0/exxb
b(1,2)=-poixyb/exxb
b(1,3)=-poixzb/exxb
b(2,1)=b(1,2)
b(2,2)=1.d0/eyyb
b(2,3)=-poiyzb/eyyb
b(3,1)=b(1,3)
b(3,2)=b(2,3)
b(3,3)=1.d0/ezzb
c
c 木材の軸（繊維）方向の線膨張率：tbz
tbz=3.5d-6
c 木材の軸直角方向の線膨張率：tbx,tby
tbx=30.d-6
tby=30.d-6
c 参考：産業調査会『木材活用事典』（事典出版センター）
c
c 木材の軸（繊維）方向の収縮率（含水率1%当たり）：hbz
hbz=0.015d-2
c 木材の軸直角方向の収縮率（含水率1%当たり）：hbx,hby
hbx=0.2d-2
hby=0.2d-2
c 参考：木材活用辞典,(株)産業調査会辞典出版センター(1994)
c http://www.cfqlcs.go.jp/technical_information/tree_and_life/wl38.htm
c
c 温度と湿度の変化による木材の歪：epbx,epby,epbz
epbx=tbx*dt + hbx*dh
epby=tby*dt + hby*dh
epbz=tbz*dt + hbz*dh
c
c write(*,*) 'epa,epbx,epby,epbz=',epa,epbx,epby,epbz
c
c 接着剤の材料定数
c 木質複合建築構造技術の開発平成１３年度報告書−構造分科会−
c 国土交通省国土技術政策総合研究所、独立行政法人建築研究所、
c （財）日本建築センター, p.101-p.108
c ２．６ 接着耐久性に関する検討
c ２．６．３ 実験および解析
c （１）ペンギンセメント#1031の木材／ブラスと処理鋼板接着性
c および樹脂硬化物の物性について
c
c ヤング率：exxc,eyyc,ezzc
exxc=2.4d9
eyyc=exxc
ezzc=exxc
c ポアソン比：poixya,poixza,poiyza
poixyc=0.41
poixzc=poixyc
poiyzc=poixyc
c せん断弾性係数：gxya,gxza,gyz
gxyc=ezzc/2.d0/(1.d0+poixyc)
gxzc=gxyc
gyzc=gxyc
c
c aに柔性（コンプライアス）行列を入れる
c(1,1)=1.d0/exxc
c(1,2)=-poixyc/exxc
c(1,3)=-poixzc/exxc
c(2,1)=c(1,2)
c(2,2)=1.d0/eyyc
c(2,3)=-poiyzc/eyyc
c(3,1)=c(1,3)
c(3,2)=c(2,3)
c(3,3)=1.d0/ezzc
c
c サンスター技研株式会社
c ケミカル事業部新規開発グループマーケティング部
c 【＃1090 の特徴】
c http://www.sunstar-engineering.com/sunstar-engineering/
c chemical/pdf/secondary.pdf
c
c 接着剤の線膨張率：ac
ac=100.d-6
c 接着剤の温度歪：epc
epc=ac*dt
c
c 鋼材、接着剤、木材の柔性行列の逆行列と、せん断弾性係数をそれぞれ
c 66a,e66c,e66bに代入し、それぞれの応力−ひずみ行列を作る。
call gyaku(a,gxya,gxza,gyza,e66a)
call gyaku(b,gxyb,gxzb,gyzb,e66b)
call gyaku(c,gxyc,gxzc,gyzc,e66c)
c
c 要素ごとの節点温度（湿度）荷重ベクトル：fab(n,24)を求める
c nは以下の分類番号
c yokoa*tateaの鋼材： 1
c yokoc*tateaの接着剤：2
c yokob*tateaの集成材：3
c yokoa*tatebの集成材：4
c yokoc*tatebの集成材：5
c yokob*tatebの集成材：6
call ondo(yokoa,tatea,ziku1,e66a, epa, epa, epa,fab,1,bt1)
call ondo(yokoc,tatea,ziku1,e66c, epc, epc, epc,fab,2,bt2)
call ondo(yokob,tatea,ziku1,e66b,epbx,epby,epbz,fab,3,bt3)
call ondo(yokoa,tateb,ziku1,e66b,epbx,epby,epbz,fab,4,bt4)
call ondo(yokoc,tateb,ziku1,e66b,epbx,epby,epbz,fab,5,bt4)
call ondo(yokob,tateb,ziku1,e66b,epbx,epby,epbz,fab,6,bt6)
c
open(9,file='ondo.d')
write(9,*) '分類番号1〜4の要素の節点温度荷重'
write(9,*) 'fx'
do i=1,22,3
write(9,9) fab(1,i),fab(2,i),fab(3,i),fab(4,i),fab(5,i),fab(6,i)
end do
write(9,*) 'fy'
do i=2,23,3
write(9,9) fab(1,i),fab(2,i),fab(3,i),fab(4,i),fab(5,i),fab(6,i)
end do
write(9,*) 'fz'
do i=3,24,3
write(9,9) fab(1,i),fab(2,i),fab(3,i),fab(4,i),fab(5,i),fab(6,i)
end do
close(9)
9 format(1p6d11.3)
c
c tyokuon.fで出力された節点変位（fort.12に出力）に、ouryogu.fで
c 歪-変位行列と応力-歪行列をかけて各要素の応力を求める時のために、
c 要素分割数、要素ごとの分類番号、分類番号ごとの歪-変位行列を
c フォーマットなしでfort.11に書き出しておく。
write(11) nx,ny,nz,ijkab,e66a,e66b,e66c,bt1,bt2,bt3,bt4,bt5,bt6
rewind(11)
c
c 要素ごとの節点温度荷重を（要素の8節点の節点番号ごとに割り振って）
c 重ね合わせて全体の温度荷重ベクトル：f
c を求める
c 全節点数：nset
nset=(nx+1)*(ny+1)*(nz+1)
c 全自由度数：nziyuu
nziyuu=nset*3
do 71 i=1,nziyuu
f(i)=0.d0
71 continue
c
do 1001 k=1,nz
do 1001 j=1,ny
do 1001 i=1,nx
n8(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
n8(2)=n8(1)+1
n8(3)=n8(1)+(nx+1)
n8(4)=n8(3)+1
n8(5)=n8(1)+(nx+1)*(ny+1)
n8(6)=n8(2)+(nx+1)*(ny+1)
n8(7)=n8(3)+(nx+1)*(ny+1)
n8(8)=n8(4)+(nx+1)*(ny+1)
c
ijk=nx*ny*(k-1)+nx*(j-1)+i
nab=ijkab(ijk)
c
do 1001 m=1,3
do 1001 jj=1,8
j24=(n8(jj)-1)*3
c s(i24+l,j24+m)=s(i24+l,j24+m)+sm(l,m,ii,jj)
jj1=jj-1
f(j24+m)=f(j24+m)+fab(nab,3*jj1+m)
1001 continue
c
c
c ***** 出力の設定 *****
c 比較のために出力したい梁−柱理論値：nhikak
c nhikak=1: w 方向に引っ張った時の先端の伸びを出力
c nhikak=2: y 方向に曲げた時の先端のたわみを出力
c nhikak=3: x 方向に曲げた時の先端のたわみを出力
nhikak=2
c
c
c 変位を画面出力したい節点の数：ngamen
ngamen=1
c 1つ目の出力したい節点番号：ngame1
ngame1=(nx+1)*(ny+1)*nz
& +1
c
c 変位をファイル出力したい節点の数：nfairu
c 例えば、先端の節点すべてだったら：
nfairu=(nx+1)*(ny+1)
c 1つ目の出力したい節点番号：nfair1
nfair1=(nx+1)*(ny+1)*nz+1
c
c
c hippon.f 用の入力データを hippon.d というファイルに出力する
c
open(7,file='hippon3.d')
write(7,1000) nx
write(7,1000) ny
write(7,1000) nz
c 全要素数：nyou
nyou=nx*ny*nz
do 70 i=1,nyou
write(7,1000) ijkab(i)
70 continue
write(7,500) yokoa
write(7,500) tatea
write(7,500) yokob
write(7,500) tateb
write(7,500) yokoc
write(7,500) tatec
write(7,500) ziku
c
c 境界条件
write(7,1000) kkasyo
write(7,1000) kkais1
write(7,1000) kryou1
write(7,500) bu1
write(7,500) bv1
write(7,500) bw1
write(7,1000) kkais2
write(7,1000) kryou2
write(7,500) bu2
write(7,500) bv2
write(7,500) bw2
c 上下対称面（zx面）のy方向変位を拘束
do k=1,nz
write(7,1000) (nx+1)*(ny+1)*k+1
write(7,1000) (nx+1)*(ny+1)*k+nx+1
write(7,500) 1.d0
write(7,500) 0.d0
write(7,500) 1.d0
end do
c 固定端の左右対称面(y軸)のx,z方向変位を固定
do j=0,ny
write(7,1000) (nx+1)*j+1
write(7,1000) (nx+1)*j+1
write(7,500) 0.d0
write(7,500) 1.d0
write(7,500) 0.d0
end do
c 左右対称面(yz面)のx方向変位を固定
do k=1,nz
do j=0,ny
write(7,1000) (nx+1)*(ny+1)*k+(nx+1)*j+1
write(7,1000) (nx+1)*(ny+1)*k+(nx+1)*j+1
write(7,500) 0.d0
write(7,500) 1.d0
write(7,500) 1.d0
end do
end do
c 対称線（z軸）のxy方向変位を拘束
do k=1,nz
write(7,1000) (nx+1)*(ny+1)*k+1
write(7,1000) (nx+1)*(ny+1)*k+1
write(7,500) 0.d0
write(7,500) 0.d0
write(7,500) 1.d0
end do
c 対称点（原点）を固定
write(7,1000) 1
write(7,1000) 1
write(7,500) 0.d0
write(7,500) 0.d0
write(7,500) 0.d0
c
c
c *********** 載荷条件 ************
c
c 境界条件に合わせて温度荷重を与えるのは面倒なので
c 全節点に温度荷重を与えてしまって、mage2.f で
c 境界条件を入れる時に、拘束節点の荷重ベクトルの成分を
c 0にする。
c 全荷重ベクトルfそのものをtyokuon.fに渡すので、
c 載荷箇所とか載荷開始点終了点は要らない。
do 20 i=0,nset-1
write(7,500) f(i*3+1)
write(7,500) f(i*3+2)
write(7,500) f(i*3+3)
20 continue
c
do i=1,3
do j=1,3
write(7,900) a(i,j)
end do
end do
write(7,500) gxya
write(7,500) gxza
write(7,500) gyza
c
do i=1,3
do j=1,3
write(7,900) b(i,j)
end do
end do
write(7,500) gxyb
write(7,500) gxzb
write(7,500) gyzb
c
do i=1,3
do j=1,3
write(7,900) c(i,j)
end do
end do
write(7,500) gxyc
write(7,500) gxzc
write(7,500) gyzc
c
c
write(7,500) df
write(7,1000) nhikak
c 端部の図心変位のみを出す場合
write(7,1000) ngamen
write(7,1000) ngame1
c 端部の全節点の変位をファイル出力
write(7,1000) nfairu
do i=nfair1,nset
write(7,1000) i
end do
write(7,1000) nstep
close(7)
c 小数点以下5桁の実数
500 format(1pd13.5)
c 小数点以下9桁の実数
900 format(1pd17.9)
c 10桁までの整数
1000 format(i10)
stop
end
c
c
c 要素番号(i,j,k)の要素の分類番号を決めるサブルーチン
c yokoa*tateaの鋼材： 1
c yokoc*tateaの接着剤：2
c yokob*tateaの集成材：3
c yokoa*tatebの集成材：4
c yokoc*tatebの集成材：5
c yokob*tatebの集成材：6
c
subroutine bunrui(nx,ny,nz,nxa,nya,nxc,ijkab)
implicit real*8 (a-h, o-z)
dimension ijkab(9400)
c もしnxa=0,nya=0（鋼材なし）だったら100に飛ぶ
if( (nxa.eq.0).and.(nya.eq.0)) goto 100
c
c 左下の鋼板要素：横nxa*縦nya分割
do i=1,nxa
ja=ny-nya+1
do j=ja, ny
do k=1,nz
ijk=nx*ny*(k-1)+nx*(j-1)+i
ijkab(ijk)=1
end do
end do
end do
c もしnxa=nx,nya=ny(全部鋼材)だったら、ここで戻る
if( (nxa.eq.nx).and.(nya.eq.ny) ) return
c
c 鋼板右隣の接着剤要素：横(nxc)*縦(nya)分割
do i=nxa+1,nxa+nxc
ja=ny-nya+1
do j=ja, ny
do k=1,nz
ijk=nx*ny*(k-1)+nx*(j-1)+i
ijkab(ijk)=2
end do
end do
end do
c
c 鋼板右隣の木材要素：横(nx-nxc-nxa)*縦(nya)分割
iac=nxa+nxc
do i=iac+1,nx
ja=ny-nya+1
do j=ja, ny
do k=1,nz
ijk=nx*ny*(k-1)+nx*(j-1)+i
ijkab(ijk)=3
end do
end do
end do
c もし鋼板が上から下までぶっとおっていたら、ここで戻る
if(nya.eq.ny) return
c
c 鋼板上の左端の木材要素：横nx*縦(ny-nya)分割
do i=1,nxa
do j=1,ny-nya
do k=1,nz
ijk=nx*ny*(k-1)+nx*(j-1)+i
ijkab(ijk)=4
end do
end do
end do
c
c 接着剤要素の上の木材要素：横(nxc)*縦(ny-nya)分割
do i=nxa+1,nxa+nxc
do j=1,ny-nya
do k=1,nz
ijk=nx*ny*(k-1)+nx*(j-1)+i
ijkab(ijk)=5
end do
end do
end do
c
100 continue
c 残りの木材要素：横(nx-nxa-nxc)*縦(ny-nya)分割
iac=nxa+nxc
do i=iac+1,nx
do j=1,ny-nya
do k=1,nz
ijk=nx*ny*(k-1)+nx*(j-1)+i
ijkab(ijk)=6
end do
end do
end do
return
end
c
c **********************************************************
c ひずみ−応力行列aの逆数を取って、応力−ひずみ行列を求め、
c それとせん断弾性係数gxy,gxz,gyzとをeに代入して、
c 応力−ひずみ行列を作る。
c
subroutine gyaku(a,gxy,gxz,gyz,e)
implicit real*8(a-h, o-z)
dimension a(3,3),b(3,3),e(6,6)
det=a(1,1)*a(2,2)*a(3,3)+a(2,1)*a(3,2)*a(1,3)
& +a(3,1)*a(1,2)*a(2,3)-a(2,1)*a(1,2)*a(3,3)
& -a(3,1)*a(2,2)*a(1,3)-a(1,1)*a(3,2)*a(2,3)
b(1,1)=(-1.d0)**(1+1)*(a(2,2)*a(3,3)-a(3,2)*a(2,3))
b(1,2)=(-1.d0)**(1+2)*(a(2,1)*a(3,3)-a(3,1)*a(2,3))
b(1,3)=(-1.d0)**(1+3)*(a(2,1)*a(3,2)-a(3,1)*a(2,2))
b(2,1)=(-1.d0)**(2+1)*(a(1,2)*a(3,3)-a(3,2)*a(1,3))
b(2,2)=(-1.d0)**(2+2)*(a(1,1)*a(3,3)-a(3,1)*a(1,3))
b(2,3)=(-1.d0)**(2+3)*(a(1,1)*a(3,2)-a(3,1)*a(1,2))
b(3,1)=(-1.d0)**(3+1)*(a(1,2)*a(2,3)-a(2,2)*a(1,3))
b(3,2)=(-1.d0)**(3+2)*(a(1,1)*a(2,3)-a(2,1)*a(1,3))
b(3,3)=(-1.d0)**(3+3)*(a(1,1)*a(2,2)-a(2,1)*a(1,2))
do i=1,3
do j=1,3
e(i,j)=b(j,i)/det
e(i,3+j)=0.d0
e(3+i,j)=0.d0
e(3+i,3+j)=0.d0
end do
end do
e(4,4)=gxy
e(5,5)=gxz
e(6,6)=gyz
return
end
c
c
c
c 歪−変位行列（いわゆるBマトリクス）の1要素の体積積分を求め、
c それの転置を bt とすると、温度歪に左から、応力-歪行列をかけて
c 更にbtをかけたもの(bt*e66*ep)を温度荷重にする
c 変位関数の定義に用いた変位d1-d8の番号は、Melosh とは違い、
c z 5__6
c ／ ／ ／|
c O →x 1−2 |
c ↓ | | ／8
c y 3−4
c だから、荷重ベクトルの重ね合わせで使っているn8の添え字が、
c tyoku2.f とは違っている。
subroutine ondo(yoko1,tate1,ziku1,e66,epx,epy,epz,fab,nab,bt)
implicit real*8 (a-h, o-z)
dimension fab(6,24),bt(24,6),e66(6,6),ep(6),ee(6)
do i=1,24
fab(nab,i)=0.d0
do j=1,6
bt(i,j)=0.d0
end do
end do
c
a=yoko1/2.d0
b=tate1/2.d0
c=ziku1/2.d0
bt( 1,1)=-b*c
bt( 4,1)= b*c
bt( 7,1)=-b*c
bt(10,1)= b*c
bt(13,1)=-b*c
bt(16,1)= b*c
bt(19,1)=-b*c
bt(22,1)= b*c
c
bt( 2,2)=-a*c
bt( 5,2)=-a*c
bt( 8,2)=+a*c
bt(11,2)=+a*c
bt(14,2)=-a*c
bt(17,2)=-a*c
bt(20,2)=+a*c
bt(23,2)=+a*c
c
bt( 3,3)=-a*b
bt( 6,3)=-a*b
bt( 9,3)=-a*b
bt(12,3)=-a*b
bt(15,3)= a*b
bt(18,3)= a*b
bt(21,3)= a*b
bt(24,3)= a*b
c
c 温度によるせん断歪がないなら、tyokuon.fで節点変位を求める上では、
c bt(i,j)のj>4の成分も必要は必要ないだろうが、
c 鋼と木の熱膨張が異なることによってせん断歪が生じるなら、
c ouryoku.fで応力を求める時のために以下の成分が必要になる。
c
bt( 1,4)= a*c
bt( 2,4)=-b*c
bt( 4,4)= a*c
bt( 5,4)= b*c
bt( 7,4)=-a*c
bt( 8,4)=-b*c
bt(10,4)=-a*c
bt(11,4)= b*c
bt(13,4)= a*c
bt(14,4)=-b*c
bt(16,4)= a*c
bt(17,4)= b*c
bt(19,4)=-a*c
bt(20,4)=-b*c
bt(22,4)=-a*c
bt(23,4)= b*c
c
bt( 1,5)=-a*b
bt( 3,5)=-b*c
bt( 4,5)=-a*b
bt( 6,5)= b*c
bt( 7,5)=-a*b
bt( 9,5)=-b*c
bt(10,5)=-a*b
bt(12,5)= b*c
bt(13,5)= a*b
bt(15,5)=-b*c
bt(16,5)= a*b
bt(18,5)= b*c
bt(19,5)= a*b
bt(21,5)=-b*c
bt(22,5)= a*b
bt(24,5)= b*c
c
bt( 2,6)=-a*b
bt( 3,6)= a*c
bt( 5,6)=-a*b
bt( 6,6)= a*c
bt( 8,6)=-a*b
bt( 9,6)=-a*c
bt(11,6)=-a*b
bt(12,6)=-a*c
bt(14,6)= a*b
bt(15,6)= a*c
bt(17,6)= a*b
bt(18,6)= a*c
bt(20,6)= a*b
bt(21,6)=-a*c
bt(23,6)= a*b
bt(24,6)=-a*c
c
c 温度ひずみ
ep(1)=epx
ep(2)=epy
ep(3)=epz
ep(4)=0.d0
ep(5)=0.d0
ep(6)=0.d0
c
do i=1,6
ee(i)=0.d0
ee(i+3)=0.d0
do j=1,6
ee(i)=ee(i)+e66(i,j)*ep(j)
end do
end do
c
do i=1,24
fab(nab,i)=0.d0
do j=1,6
fab(nab,i)=fab(nab,i)+bt(i,j)*ee(j)
end do
end do
c
do i=1,24
do j=1,6
bt(i,j)=bt(i,j)/a/b/c/8.d0
end do
end do
c
return
end
c
c
c