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 mozi-koodo: EUC
c kaigyou-koodo: LF
c 集成材梁の有限変位・弾塑性解析プログラム。
c 圧縮側の軸方向応力が降伏応力を越えたら塑性化。
c 引張側の軸方向応力が破断応力を越えたら破壊。
c 金属材料は想定してないので、von Mises の降伏条件は使ってない。
c 断面は長方形のみ。
c 載荷履歴に依存する弾塑性問題では、不釣り合い力をトータルな値
c として扱えないので（後述）、荷重増分を十分に小さく与えること。
c その意味で（接線剛性方程式の定式化はトータルラグラジアンだが）
c 疑似的なアップデートラグラジアンになっている。
c harisen.f から改造。
c
c 定式化は、
c 『後藤文彦，小林 裕，斉木 功，岩熊哲夫：
c 空間固定三軸回りの回転自由度を用いた空間梁解析
c 応用力学論文集, 1, 1998, pp.319-327.』
c の式(48)で座標変換をオイラー角で表した接線剛性方程式を用いた。
c 尚、反りに関して、ねじり率の自由度は、
c 『後藤文彦，小林裕，岩熊哲夫：
c オイラー角を用いた簡潔な有限変位解析手法
c 構造工学論文集, 43A, 1997, pp.333-338.
c の定式化のように扱った。
c
c 04/8/26版
c プログラムを主に作った人：後藤文彦
c 増分式を Mathematica で計算してくれた人：小林裕さん
c いくつかの代数計算やスカイライン法のサブルーチンの
c 元プログラムの作者：岩熊先生(http://www.civil.tohoku.ac.jp/~bear/)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c このプログラムは無保証です。 c
c このプログラムは改造／再配布して構いません。 c
c 但し、私的利用に留まらない場合は、もとのプログラムが c
c http://www.str.ce.akita-u.ac.jp/~gotou/programoj/ c
c から取ってきたものと分かるように明示して下さい。 c
c 改造／再配布後も上記の条件を保持して下さい。 c
c バグの報告は、 c
c http://www.str.ce.akita-u.ac.jp/~gotou/jpg/meeru.jpg c
c にお願いします。 c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
program harisori
implicit real*8(a-h,o-z)
dimension nnl(256+1),nnr(256+1),lorg(256*4),lg4(4),
& noutl(256),noutr(256),nuvwl(256),nuvwr(256),
& f(7*256+8),bc(7*256+8),d(14*256),dd(7*256+8),
& ddl(7*256+8),b(7*256+8),alp(256),gmm(256),phi(256),
& d14(14),b14(14),dd14(14),sko(7,7),s(14,14),sk(14,14),
& stuh(84*256+36),stlh(84*256+36),imx(7*256+9),
& msl(7*256+8),nnegtm(7*256+8),anegtm(7*256+8),
& dld(14*256),vec14(14),sgm(256,200,200),q(7*256+7),
& szenbu(256,14,14),dld14(14),ts14(14,14)
c
c 頭にclと付くサブルーチンは、配列類の初期化
call cl4ivs(256+1,nnl,256+1,nnr,256*4,lorg,4,lg)
call cl4ivs(256,noutl,256,noutr,256,nuvwl,256,nuvwr)
call cl4rvs(7*256+8,f,7*256+8,bc,16*256,d,7*256+8,dd)
call cl4rvs(7*256+8,ddl,7*256+8,b,256,alp,256,gmm)
call cl4rvs(256,phi, 14,d14, 14,b14, 14,dd14)
call cl4rms(7,sko, 14,s, 14,sgz, 14,sgx)
call clrmx(14,sgy)
call clrmx(14,sk)
call clrv(84*256+36,stuh)
call clrv(84*256+36,stlh)
call clrv(7*256+7,q)
call cliv(7*256+9,imx)
call cliv(7*256+8,msl)
call cliv(7*256+8,nnegtm)
call clrv(7*256+8,anegtm)
c
c 変数の初期化
ninc=0
ddi2=0
ddl2=0
nlet=0
p1=0.d0
nbcl=0
nacl=0
iyonda=0
c
c 一発目の弧長が大きすぎて座屈してしまった場合に、初期状態から
c 計算しなおす時のために、初期変位（ゼロ）をフォーマットなしで fort.11
c に書き出しておく。
write(11) dd,d,ddl,dp,dpl,ro,p
rewind 11
c 荷重−変位（u,v,w）関係を出力するファイル：uvw.out
open(8,file='uvw.out')
c
c ### まずはデータの入力 ###（変数の説明を兼ねる）
c
c 要素数：nofe
read(*,*) nofe
write(*,*)'要素数:',nofe
c 座屈モード：mode
read(*,*) mode
write(*,*) mode,'番目に低い不安定点を捜します'
c 前に計算したモードに対する倍率：gmode
c 座屈後などを計算する場合は、座屈点付近の接線剛性行列の
c 固有値解析をして得た変位のモードに、倍率 gmode を掛けたものを
c 外乱や（初期）不整として取り込む場合に指定する。
c 不整を与えずに完全形の変位や座屈点などを求める場合は
c gmode=0 を与える。
read(*,*) gmode
c モード計算のための出力が必要：modeo=1
c モード計算のための出力が不要：modeo=0
c 上述のようなモード計算のために、計算終了時（収束点）の
c 接線剛性行列を出力する場合は modeo=1
read(*,*) modeo
nmax=7*256+8
n7=7*nofe+7
n8=7*nofe+8
call maxa(nofe,n8,imx)
imx9=imx(n8+1)
c
c 収束精度：収束時の荷重が、N-R反復における一つ前のステップの
c 荷重と有効数字何桁まで合うようにしたいかに応じて tol の小ささを
c 決める
tol=1.d-10
c 外乱の大きさに応じて収束精度を甘くしたいときとかは、
cc if(gmode.gt.1.d-10) tol=tol*1.d3
c みたいなのを入れたりするとか。cnvr は座屈点探査の精度。
c
c 初期状態における全ての要素の方向をオイラー角で与える
c （このオイラー角は右手系ではないので注意）
do 8 i=1,nofe
read(*,510) alp(i),gmm(i),phi(i)
8 continue
510 format(1p3d17.9)
c
c 断面定数
c ヤング率：e
read(*,*) e
c 降伏後の軸（繊維）方向接線弾性係数：et
read(*,*) et
c 軸圧縮降伏応力：sgmy
read(*,*) sgmy
c 圧縮：負なので sgmy の符号を負にしておく
sgmy=-abs(sgmy)
c 軸引張破断応力 sgmu
read(*,*) sgmu
c せん断弾性係数：g
read(*,*) g
c 長方形断面の幅 haba
read(*,*) haba
c 長方形断面の高さ taka
read(*,*) taka
c 長方形断面の幅の断面要素分割数 nhaba
read(*,*) nhaba
c 長方形断面の高さの断面要素分割数 ntaka
read(*,*) ntaka
c 断面積：a
read(*,*) a
ea=e*a
c x軸回りの断面二次モーメント:xi
read(*,*) xi
eix=e*xi
c y軸回りの断面二次モーメント:yi
read(*,*) yi
eiy=e*yi
c ねじり定数：cj
read(*,*) cj
gj=g*cj
c そりねじり定数:wi
read(*,*) wi
eiw=e*wi
c el:一要素の長さ
read(*,*) el
c 梁の軸長：ell
ell=el*nofe
c y軸回りの曲げによるxz面のせん断変形に対するせん断補正係数:xzk
read(*,*) xzk
c x軸回りの曲げによるyz面のせん断変形に対するせん断補正係数:yzk
read(*,*) yzk
c 構造工学シリーズ7『構造工学における計算力学の基礎と応用』(土木学会)
c のp.94,95 のφxy,φxz に対応するパラメーター：pxz, pyz
pxz=eiy/(g*xzk*a*el**2)
pyz=eix/(g*yzk*a*el**2)
c enko.fでyz面内に横たわるアーチにした場合の中心角：th
th=-alp(nofe)*2.d0
c π(rad)：pi
pi=asin(1.d0)*2.d0
c 比較対照解の種類：nhikak
c nhikak=1:圧縮単純梁のオイラー座屈の解
c nhikak=2：等曲げ円弧アーチの横ねじれ座屈のウラソフの解
c nhikak=3: 片持ち梁の横ねじれ座屈のトラヘアの解（面内変位無視）
read(*,*) nhikak
goto (511,512,513), nhikak
511 continue
c 柱のオイラー座屈公式
eul=(pi/ell)**2*eiy
write(*,*) 'Euler=',eul
goto 599
512 continue
c 等曲げを受ける円弧アーチの横ねじれ座屈のウラソフの解
eix2=eix*1.d6
call vlasov(mode,ell,eix2,eiy,gj,eiw,th,pi,vln,vlp)
write(*,*)'Vlasov（負）=',vln
write(*,*)'Vlasov（正）=',vlp
c と修正ウラソフの解（座屈前の面内変位の影響考慮）
call vlasov(mode,ell,eix,eiy,gj,eiw,th,pi,vln,vlp)
write(*,*)'修正Vlasov（負）=',vln
write(*,*)'修正Vlasov（正）=',vlp
goto 599
513 continue
c 片持ち梁の横ねじれ座屈のトラヘアの解（座屈前面内変位無視）
a1=sqrt(eiy*gj)/ell/ell
a2=sqrt(eiw/gj)*pi/ell
tra=a1*(3.95d0+3.52d0*a2)
write(*,*) 'Trahair=',tra
goto 599
599 continue
c
c
c ### 載荷条件 ###
c 載荷されている節点の数：nofln
read(*,*) nofln
c 載荷節点の番号：nnl
c 載荷節点における変位ベクトル：f(自由度数+1)
do 10 i=1,nofln
read(*,*) nnl(i)
do 20 j=1,7
nf=7*(nnl(i)-1)+j
read(*,*) f(nf)
c 荷重は el と eiy で無次元化
c これに対応して変位も無次元量で求まる。
if(j.lt.4) f(nf)=f(nf)*el*el/eiy
if(j.gt.3) f(nf)=f(nf)*el/eiy
if(j.eq.7) f(nf)=f(nf)/eiy
20 continue
10 continue
c
c ### 境界条件 ####
do 30 i=1,n8
bc(i)=1.d0
30 continue
c 拘束されている節点の数：nofrn
read(*,*) nofrn
c 拘束節点の番号：nnr
c と拘束節点における変位ベクトル：bc(自由度数+1)
do 40 i=1,nofrn
read(*,*) nnr(i)
do 50 j=1,7
nb=7*(nnr(i)-1)+j
read(*,*) bc(nb)
50 continue
40 continue
c 境界条件を局所座標系で与える（節点変位をそこだけ局所系で
c 定義して全体系に変換する座標変換を掛けないでおく）要素の数：noflbc
c （円弧アーチの横倒れ座屈を解く際に、端部の面外回転を半径軸回りに
c 与える場合など）
read(*,*) noflbc
c その要素番号を入れ、次の四項目の境界条件をどう与えるのか
c （局所系1／全体系0）を lorg(要素数*4) に入れる。四項目とは、
c 左の並進、左の回転、右の並進、右の回転
do 60 i=1,noflbc
read(*,*) lbcen
l1=lbcen-1
read(*,500) lorg(l1*4+1),lorg(l1*4+2),lorg(l1*4+3),lorg(l1*4+4)
60 continue
500 format(4i2)
c
c ### 出力条件 ###
c
c 左の節点変位を出力したい要素の要素数：noflo, その要素番号：noutl
read(*,*) noflo
do 70 i=1,noflo
read(*,*) noutl(i)
70 continue
c 右の節点変位を出力したい要素の要素数:nofro, その要素番号:noutr
read(*,*) nofro
do 80 i=1,nofro
read(*,*) noutr(i)
80 continue
c 左の節点並進変位をファイル出力したい要素の要素数:noflof
c その要素番号: nuvwl
read(*,*) noflof
do 75 i=1,noflof
read(*,*) nuvwl(i)
75 continue
c 右の節点並進変位をファイル出力したい要素の要素数:nofrof
c その要素番号: nuvwr
read(*,*) nofrof
do 85 i=1,nofrof
read(*,*) nuvwr(i)
85 continue
c
c 軸応力を（gnuplot作図用に）出力したい要素の要素番号：nsgm
read(*,*) nsgm
c
c 一番最初の一発目の増分ステップの弧長：r
read(*,*) r
c ２ステップ目以降の計算に使う弧長増分：rn
read(*,*) rn
c
c
c ############# ここまでがデータ入力 #############
c
c 座屈点探査の精度：「アキレスと亀」（の系：分割のパラドックス
c 「目的地に達するにはその半分の点に達しなければならず、その半分の
c 点に達するには更にその半分の点に達しなければならず……」）
c 方式で座屈点を探査する際に、前の荷重ステップのN-R収束荷重と、
c 今の荷重ステップのN-R収束荷重とが、
c 有効数字何桁まで合うようにしたいかに応じて cnvr の小ささを
c 決める
cnvr=5.d-8
c 各N-R収束過程における収束精度は tol で与える。
c
c N-R収束荷重pが前のステップの収束荷重p1とある精度内で同じ(p1=p)に
c なった回数：ip1p
ip1p=0
c
c ### 初期条件 ###
c
c 変位増分ベクトル：dd(自由度数+1）
c 変位増分の回転角成分は空間固定3軸回りの微少回転角
c 変位ベクトル：d(節点自由度数*要素数）
c トータルな変位は、要素どうしが繋がる２節点のオイラー角を
c 別々に覚えておく。
c ddl は dd を各増分ステップの繰り返しで総和したもの
c dp は荷重増分（弧長法なので、増分変位と同様に解かれる未知変数）
do 100 i=1,n7
dd(i)=0.d0
d(i)=0.d0
ddl(i)=0.d0
100 continue
dp=0.d0
p=0.d0
c 各断面要素の軸応力 sgm(要素番号,x座標,y座標)の初期化
do i=1,nofe
do l=1,200
do m=1,200
sgm(i,l,m)=0.d0
end do
end do
end do
c 弧長の初期値
ro=r
dpl=r
c まだ一度も収束値が求まっていないことを表す:nccl（後述）
nccl=0
c ### 各荷重ステップの一発目の弧長増分過程 ###
1000 continue
nofcs=0
c 計算開始からの増分繰り返し数：ninc
ninc=ninc+1
c if(ninc.gt.6000) then
c write(*,*)'(増分繰り返しが 6000 を越えました)'
c close(8)
c stop
c endif
c
c 全変位増分ddと全変位dの中から、1要素ぶんずつをそれぞれ
c dd14とd14に抜き出す
do 1100 i=0,nofe-1
do 120 j=1,14
dd14(j)=dd(i*7+j)
d14(j)=d(i*14+j)
c
c dldは各増分ステップで（1発目の弧長+N-R過程で）収束した荷重増分
c 要素ごとに前のステップのdld14をdldから読み出す。
dld14(j)=dld(i*14+j)
c dld14を読み出したらdldはゼロに戻して良い。
dld(i*14+j)=0.d0
c
120 continue
do 130 j=1,4
lg4(j)=lorg(i*4+j)
130 continue
c
c 前のステップのトータルな d14 で要素ごとに座標変換 ts14 を作る
call trns(ts14,d14(4),d14(5),d14(6))
c （直線梁に限定するので初期形状の座標変換は作らない）
c ts14 の転置を掛けて居所系の dld14 を求める
call mxtv(14,ts14,dld14,vec14)
c （直線梁に限定するので初期形状の座標変換はなし）
c
c dd14とd14を使って1要素ずつ接線剛性行列を計算していく
c 求まった1要素ぶんの接線剛性行列はskに入る
c
c ----------------------------------------
c 要素番号：nyou
nyou=i+1
c 前のステップのdld14を使って剛性行列を作る
c まずは断面要素の降伏判定をして断面定数を断面積分で求める
c 弾性なら、c11=EA,c22=EIy,c33=EIx,c44=EIw,c55=GJ
c
call danmen(nyou,nofe,nhaba,ntaka,el,haba,taka,
& e,et,g,sgmy,sgmu,sgm,vec14,pi,c11,c22,c33,c44)
c
c 求まった断面定数で線形剛性行列（反りを含む1節点7自由度）を
c s(14,14) s に入れる
c
call smx(el,eiy,pxz,pyz,c11,c22,c33,c44,gj,s)
c
c 要素ごとの剛性行列 s を全要素について szenbu に覚えさせる
do ii=1,14
do ij=1,14
szenbu(nyou,ii,ij)=s(ii,ij)
end do
end do
c
c stmxに行く前にdld14はゼロにする必要がある
call clrv(14,dld14)
c
c 接線剛性行列を計算する。
call stmx(nofe,i,alp,gmm,dd14,d14,sk,b14,lg4,s,dld14)
c
do 140 j=1,14
d(i*14+j)=d14(j)
140 continue
c
c stmx から返ってきたdld14は局所系だから全体系に変換する。
c ts14 を掛けて全体系の vec14 を求める
call mxv(14,ts14,dld14,vec14)
cc vec14 を固定座標まで回転させる ヤメ（理由は忘れだけど：
cc call mxv(14,ti14,dld14,vec14)
cc 初期座標方向の局所系で dld14 を定義すべきと判断したのかな）
c
c このステップの一発めの dld14 を dld の中に要素全部に
c ついて覚えさせる．
do j=1,14
dld(i*14+j)=vec14(j)
end do
c
c sk をスカイライン法の1次元ベクトルに書き込む
c stuh に上三角、stlh に下三角が入る
call skyln(nofe,i,sk,sko,stuh,stlh,imx)
1100 continue
c
cccccccccccccccccccccccccccccccccccccccccccccccc
c
do 150 i=1,n7
i8=imx(n8)+n8-i
stuh(i8)=-f(i)
stlh(i8)=ddl(i)/ro
ddl(i)=0.d0
b(i)=0.d0
150 continue
do 160 i=1,n8
stlh(imx(i))=b(i)
160 continue
stuh(imx(n8))=dpl/ro
stlh(imx(n8))=r
ro=r
dpl=0.d0
c 境界条件を入れる
call bncn(nofe,bc,stuh,stlh,imx)
c 代数方程式をガウスの消去法で解いて変位ベクトルを求め、
c その際に対角項の符号から、行列の負の固有値の数 jdt を求める
call solve(stuh,stlh,dp,imx,msl,nnegtm,anegtm,
& det,idt,jdt,n8)
do 170 i=1,n8
dd(i)=stlh(imx(i))
170 continue
c
c ### ニュートン−ラフソン 収束過程 ###
c（基本的に上の「一発目の弧長増分過程」の計算と同じ）
2000 continue
c 各増分ステップでの収束までの繰り返し数：nofcs
nofcs=nofcs+1
do 180 i=1,n7
ddl(i)=ddl(i)+dd(i)
ddi2=ddi2+dd(i)**2
ddl2=ddl2+ddl(i)**2
180 continue
dpl=dpl+dp
p=p+dp
arc=ddi2/ddl2
ddi2=0
ddl2=0
if(nofcs.gt.6) then
write(*,*)'反復数が6を越えても収束しないので弧長を半分にする'
c 反復数が6を越えて辛うじて収束するようなことが何回も重なると、
c 答えがおかしくなってたりするので、その回数を nlet で数えて、
c 適当な回数に達したら強制的に終了
nlet=nlet+1
if(nlet.gt.10) then
write(*,*) '収束ステップが10を越えたので強制終了'
stop
endif
goto 4001
endif
if(arc.lt.tol) goto 3000
c 不釣合力 b の初期化（最終的な不釣合力を出すのは1200の下）
do 190 i=1,n7
b(i)=0.d0
cc トータルラグラジアンの場合は、トータルな不釣合力を用いる
cc b(i)=p*f(i)
190 continue
c
c N-R 収束過程の要素ごとのループのはじまり
c
do 1200 i=0,nofe-1
do 200 j=1,14
dd14(j)=dd(i*7+j)
d14(j)=d(i*14+j)
c 前の段階でのdld14をdld から読み出してくる
dld14(j)=dld(i*14+j)
b14(j)=b(i*7+j)
200 continue
c
do 210 j=1,4
lg4(j)=lorg(i*4+j)
210 continue
c
c stmx へ局所系のdldを送るため dld14 を局所系まで回転させる
cc dld14 は初期形状座標で使うので（上述）、ここでdld14を
cc 初期形状座標まで回転させる必要はない。
cc call mxtv(14,ti14,dld14,vec14) ヤメ
c dld14 に ts14 の転置を掛けて居所系の dld14 を求める
call mxtv(14,ts14,dld14,vec14)
c
c szenbu から要素ごとの剛性行列 s を読み出す
c 要素番号：nyou
nyou=i+1
do ii=1,14
do ij=1,14
s(ii,ij)=szenbu(nyou,ii,ij)
end do
end do
c
c 接線剛性行列を１要素ずつ計算してskに入れる
call stmx(nofe,i,alp,gmm,dd14,d14,sk,b14,lg4,s,dld14)
c
c print *, 'sk'
c print '(1p7d11.3)', ( (sk(i,j),j=8,14),i=8,14 )
c print '(1p7d11.3)', (b14(j),j=1,7)
c print '(1p7d11.3)', (b14(j),j=8,14)
c stop
c
c
do 220 j=1,14
d(i*14+j)=d14(j)
220 continue
c
c stmx から返ってきたdld14は局所系だから全体系に変換する。
c ts14 を掛けて全体系の vec14 を求める
call mxv(14,ts14,dld14,vec14)
cc vec14 を固定座標まで回転させる ヤメ（理由は忘れだけど：
cc call mxv(14,ti14,dld14,vec14)
cc 初期座標方向の局所系で dld14 を定義すべきと判断したのかな）
c
c このステップの一発めの dld14 を dld の中に要素全部に
c ついて覚えさせる．
do j=1,14
dld(i*14+j)=vec14(j)
end do
c
c skをスカイライン法の1次元配列に書き込む
c stuhが上三角行列、stlhが下三角行列
call skyln(nofe,i,sk,sko,stuh,stlh,imx)
do 230 j=1,14
b(i*7+j)=b14(j)
230 continue
1200 continue
ccccccccc 要素ごとのループの終わり cccccccccccccccc
c
c 毎ステップで求まる増分変位Δdで求めた KΔd を累積した
c b=b+KΔd を目標荷重 pf から引いて不釣合力を求める。
c
do i=1,n7
q(i)=q(i)+b(i)
b(i)=p*f(i)-q(i)
end do
c
do 240 i=1,n7
i8=imx(n8)+n8-i
stuh(i8)=-f(i)
stlh(i8)=ddl(i)
240 continue
do 250 i=1,n8
stlh(imx(i))=b(i)
250 continue
stuh(imx(n8))=p
stlh(imx(n8))=0.d0
c 境界条件を入れる
call bncn(nofe,bc,stuh,stlh,imx)
c 座屈点の探査のとき、荷重が前の増分ステップと誤差範囲に
c になったら終了する。モード計算用の出力が必要な場合(modeo=1）は
c fort.7 に収束時の接線剛性行列をフォーマットなしで書き出す。
if( (abs((p1-p)/p).lt.cnvr).and.(gmode.lt.1.d-20) ) then
c 弾塑性解析で増分荷重を小さくしすぎると座屈点に達する前の
c 何でもないとこでp1=pとなって終了してしまうことがあるので、
c p1=pが2回ぐらい連続しないと終了しないようにする。
ip1p=ip1p+1
print *,'前ステップと収束荷重が',ip1p,'回同じになりました'
if(ip1p.gt.1) then
rn=rn*2.d0
call sgmgnu(nsgm,nhaba,ntaka,sgm)
if(modeo.eq.1) then
write(7) nofe+1,stuh,stlh,imx,imx9
rewind 7
endif
close(7)
close(8)
c
c 出力した接線行列の固有値解析をせずに、最後のN-R反復で求まる
c （微少な）変位増分を固有ベクトル（座屈モード）の近似として
c 利用する場合のために、一応、近似モードもここで出力しておく。
c こういう近似ができるということは、
c 久田俊明・野口裕久 著『非線形有限要素法の基礎と応用』
c（丸善株式会社）p.317とかに書いてある。
c
call modeout(nofe,dd)
c
c 但し、接線剛性方程式の回転自由度は空間固定3軸回りの成分なので、
c 座屈モードの回転成分もその成分に対応して求まる。
c 座屈後へ分岐させる外乱としては、3軸回り成分をオイラー角に与えても
c 十分に用を足すとは思うが。
c オイラー角成分がほしいときは、変換行列をかけるとかして。
c
stop
endif
endif
c
c
c 接線剛性方程式を解く
call solve(stuh,stlh,dp,imx,msl,nnegtm,anegtm,
& det,idt,jdt,n8)
do 260 i=1,n8
dd(i)=stlh(imx(i))
260 continue
c
c print *, 'dd'
c print '(1p7d11.3)', (dd(j),j=8,14)
c
goto 2000
c
c ### 出力 ###
c
3000 continue
arc=0.d0
c 収束しないために弧長が小さくなり過ぎたのを、また大きくしたい
c ときは、nlet の大きさに応じて弧長を大きくしてやるとか
c if(nlet.gt.1) rn=rn*(2.d0)**nlet
c ここを通過しているということは、増分ステップが収束したという
c ことなので増分ステップ内の繰り返し数 nlet を0にする。
nlet=0
c 座屈点の探査は、座屈点を超えたら、越える前の荷重ステップまで
c 戻り、弧長を半分にして進み、また行き過ぎたらまた戻って弧長を
c 半分にして次の荷重ステップを計算して座屈点に近づいていく。
c この座屈点を跨いだいったりきたりの計算において、同じ点を重複して計算
c しないように、「アキレスと亀」方式で弧長を毎回半分にしていく。
c
c 接線剛性行列の負の固有値数 jdt と固有ベクトルとして得られる
c 座屈モードの次数とはだいたい対応していることから、jdt が
c 指定した座屈モードの次数(mode)を越えたかどうかで座屈点を
c 越えたかどうかを判断する。座屈点を超えたかどうかにも
c 様々な場合があり得るので、その状況をnacl,nbcl,ncclを使って判断する。
c
c 一発目の弧長を大きく取りすぎて、最初から座屈点を超えてしまった場合
c は、正しい答えがまだ求まっていないので、最初から計算をやり直す
c 必要がある。そのため、まだ一つも収束値が求まっていない時は
c nccl=0, 一つでも収束値が（座屈点を超える前に）求まったら
c nccl=1 とする。
if(jdt.lt.mode) nccl=1
write(*,*) 'N-R反復数=',nofcs
write(*,570) det, idt,jdt
570 format(' ','det=',f11.3,'d',i4,' , ','負の固有値数：',i4)
goto (571,572,573), nhikak
571 write(*,581) ' p=',p,' PL/EA=',p*ell/ea,' Euler=',eul
goto 579
572 write(*,581) ' p=',p,' ML/EIy=',p*ell/eiy
c write(*,581) ' p=',p,' ML^2/8EI=',p*ell**2/8.d0/eix
write(*,580) ' 修Vla正=',vlp*ell/eiy,' 修Vla負=',vln*ell/eiy
goto 579
573 continue
rmd=ell/sqrt(eiy/ea)
write(*,581)' p=',p,' PL^2/EIy=',p*ell**2/eiy,' λ=',rmd
goto 579
579 continue
c
c write(*,580)' p=',p,' Ml/EIy=',p*ell/e/yi
580 format(' ',a,1pd13.5,a,1pd13.5)
581 format(' ',a,1pd13.5,a,1pd13.5,a,1pd13.5)
c 節点変位は、u, v, w, α, γ, φ, λの順に出力。
c α, γ, φ はオイラー角、λはねじれ率
c d は無次元化されているので、出力するときは有次元化する
c
c 指定要素の左節点を出力する場合
if(noflo.ne.0) then
do 460 i=1,noflo
j=noutl(i)-1
write(*,590) noutl(i),'左',d(j*14+1)*el,d(j*14+2)*el,
& d(j*14+3)*el,d(j*14+4),d(j*14+5),d(j*14+6),d(j*14+7)/el
460 continue
endif
c 指定要素の右節点を出力する場合
if(nofro.ne.0) then
do 470 i=1,nofro
j=noutr(i)-1
write(*,590) noutr(i),'右',d(j*14+8)*el,d(j*14+9)*el,
& d(j*14+10)*el,d(j*14+11),d(j*14+12),d(j*14+13),d(j*14+14)/el
470 continue
endif
c 指定要素の左節点のuvwをファイル出力
if(noflof.ne.0) then
do 480 i=1,noflof
j=nuvwl(i)-1
write(8,560) p,d(j*14+1)*el,d(j*14+2)*el,d(j*14+3)*el
480 continue
endif
c 指定要素の右節点のuvwをファイル出力
if(nofrof.ne.0) then
do 490 i=1,nofrof
j=nuvwr(i)-1
write(8,560) p,d(j*14+8)*el,d(j*14+9)*el,d(j*14+10)*el
490 continue
endif
560 format(1p4d11.3)
590 format(i3,a,1p6d11.3,1pd10.2)
if(mode.eq.0) goto 4003
c 接線剛性行列の負の固有値数jdtが指定した座屈モード次数を
c 越えたら、座屈したと見なす。
if(jdt.ge.mode) then
write(*,*) ' zakutu'
write(*,*) ' 座屈したみでだよ................'
call sgmgnu(nsgm,nhaba,ntaka,sgm)
stop
c 最初に座屈点を超えたら、まずnbcl=1,nacl=0としておく。
nbcl=1
nacl=0
endif
if(p.lt.p1) then
write(*,*) '除荷したみでだよ：p=',p
stop
endif
4004 continue
c
c ### 外乱の入れ方 ###
c 一発目の弧長増分ステップが収束した後に入れる。
c 予め座屈点を計算しておいて、一発目の弧長を座屈点の直前の
c 分岐しやすい荷重で外乱が入るように調整する。
if (nccl.eq.1) then
c
c 「アキレスと亀」方式で座屈点に達したら、そこで外乱を読み込む
c 場合は、
c if(abs((p1-p)/p).lt.1.d-1) then
c とか。
if(abs(gmode).lt.1.d-20) goto 5002
c 但し、外乱は一度しか読み込まない
if(iyonda.eq.1) goto 5002
write(*,*)'外乱読み込み中'
c 外乱を読み込んだ後に収束精度を甘くしたいときは、
c cnvr=1.d-3
c とか。外乱を1度 読んだらiyonda=1
iyonda=1
c
c
call modein(nofe,d,gmode)
c
c
rn=rn*1.5d0
endif
5002 continue
p1=p
c 座屈点を超えたために座屈前に戻って1ステップ以上、座屈点を超えない
c 点を解いているとき（nbcl=0, nacl=1）は
c 弧長は半分にするけど、座屈してない状態だから前のステップには
c 戻らなくていい。rn=rn/2の後、4002まで飛んで1000に戻る。
if(nacl.eq.1) rn=rn/2.d0
c 座屈点を超えたばかりの状態(nbcl=1, nacl=0)のときは、
c 座屈前に戻って弧長を半分にする（それでも、また座屈した時は
c 上でnbcl=1になってるので、またここを通過する）。
c 次にここを通ときに座屈前に戻っていることを教えるために
c nbcl=0, nacl=1 とする（nbcl=nacl=1の組み合わせはない）
if(nbcl.eq.1) then
nbcl=0
nacl=1
goto 4001
endif
c こういうふうに goto が交差した書き方をしたのは、
c 一旦ここでif文を閉じておいて、上から4001へ飛べるようにするため。
goto 4002
4001 continue
read(11) dd,d,ddl,dp,dpl,ro,p
rewind 11
rn=rn/2.d0
c 一発目の弧長が大きすぎて座屈してしまった場合（nccl=0）は、
c 弧長を一割小さくして初期状態から計算しなおす。
if(nccl.eq.0) then
r=r*0.9d0
ro=r
dpl=ro
endif
4002 continue
c 一度でも収束値が求まっている(nccl=1)なら、収束するたびに、
c 変位をfort.11にフォーマットなしで書き出す。但し、
c 座屈点を超えた場合は、上記のread(11)で、座屈前の変位を読み込んだ
c ものが、ここで書き出される。つまり、ここでは常に座屈していない
c 状態での最も後に計算した変位が出力される。そして、次に座屈点を
c 越えたときに、この変位の点まで戻れるようにしておく。
if(nccl.eq.1) then
write(11) dd,d,ddl,dp,dpl,ro,p
rewind 11
r=rn
endif
4003 continue
write(*,*)'----------------------------------------------------'
c ここまでで一つの弧長増分ステップの終わり。次の弧長増分ステップ
c に戻る。
goto 1000
end
c
c ものすごく長いけど、ここまでがメインプログラム。
c 以下がサブルーチン
c
c 等曲げを受ける円弧アーチの横ねじれ座屈に対するウラソフの解析解
subroutine vlasov(n,blen,eix,eiy,gj,eiw,th,pi,vln,vlp)
implicit real*8(a-h,o-z)
if(th.lt.1.d-6) th=1.d-6
r=blen/th
f=(n*pi/blen)**2*eiw+gj
a=eiy/eix-1.d0+(1.d0/eix-eiy/eix/eix)*f
b=(eiy+(1.d0-2.d0*eiy/eix)*f)/r
c=((n*pi/blen)**2-1.d0/r/r)*eiy*f
d=b**2-4.d0*a*c
vln=(-b+sqrt(d))/2.d0/a
vlp=(-b-sqrt(d))/2.d0/a
return
end
c
c 断面要素の降伏判定をしながら断面定数を断面積分する
subroutine danmen(nyou,nofe,nhaba,ntaka,el,haba,taka,
& e,et,g,sgmy,sgmu,sgm,dld14,pi,c11,c22,c33,c44)
c
implicit real*8 (a-h,o-z)
dimension x(200),y(200),depsz1(200,200),depsz2(200,200),
& dsigz(200,200),eet(200,200),sgm(256,200,200),dld14(14)
c
c まず、この要素の各断面要素の et を求める
do l=1,nhaba
do m=1,ntaka
c
c nyou番目の要素の断面要素（座標：l,m）の軸応力：sigma
sigma=sgm(nyou,l,m)
c
c 集成材の場合は、圧縮のみ弾塑性、引張は弾性域で破断なので
c 降伏判定は圧縮側つまりsigma<0のときのみ行う。
c （引張側は破断判定のみ行う）圧縮：負なのでsgmy<0で与えてある
if(sigma.lt.sgmy) then
eet(l,m)=et
if((nyou.eq.2).and.(l.eq.1).and.(m.eq.1)) then
write(*,*) ' sosei'
write(*,*) '要素,断面座標,et,sgm=',nyou,l,m,eet(l,m),sigma
endif
c 降伏判定された sigma は sgmy に置き換えて、次の dsgm で
c 除荷されたら、eet=e になるようにしておく（除荷は e で戻るように）。
sgm(nyou,l,m)=sgmy
else
if(sigma.gt.sgmu) then
write(*,*) ' hadan'
write(*,*) ' 要素番号：',nyou,' 断面要素座標：',l,m,
& 'で破断が生じました'
write(*,*) 'sigma=',sigma
call sgmgnu(nyou,nhaba,ntaka,sgm)
stop
end if
eet(l,m)=e
end if
end do
end do
c
an=haba/float(nhaba)
bn=taka/float(ntaka)
f=an*bn
c 断面図心を原点とした場合の断面要素の座標値を求める
nhaba2=nhaba/2
ntaka2=ntaka/2
do i=1,nhaba2
x(nhaba2-i+1)= -(an*i -an/2.d0)
x(nhaba2+i)= an*i -an/2.d0
end do
do j=1,ntaka2
y(ntaka2-j+1)= bn*j -bn/2.d0
y(ntaka2+j)= -(bn*j -bn/2.d0)
end do
c
c 節点変位を以下の変数に書き写す
u1= dld14(1) *el
v1= dld14(2) *el
w1= dld14(3) *el
thx1=dld14(4)
thy1=dld14(5)
thz1=dld14(6)
q1= dld14(7) /el
u2= dld14(8) *el
v2= dld14(9) *el
w2= dld14(10)*el
thx2=dld14(11)
thy2=dld14(12)
thz2=dld14(13)
q2= dld14(14)/el
c
c 断面定数を断面積分する前に初期化
c11=0.d0
c22=0.d0
c33=0.d0
c44=0.d0
c c55=0.d0
c
c
c 要素両端の各断面要素の軸方向歪成分を計算
c
do l=1,nhaba
do m=1,ntaka
xl=x(l)
ym=y(m)
call sori(7,haba,taka,xl,ym,pi,w)
c
depsz1(l,m)=(w2-w1)/el
& -(6.d0/(el**2)*(u2-u1)+2.d0/el*(2.d0*thy1+thy2))
& *(xl-ym*thz1)
& -(6.d0/(el**2)*(v2-v1)+2.d0/el*(2.d0*thx1+thx2))
& *(ym+xl*thz1)
& -w*(6.d0/(el**2)*(thz2-thz1)+2.d0/el*(2.d0*q1+q2))
& +1.d0/2.d0*((thy1**2)+(thx1**2)
& +((xl**2)+(ym**2))*(q1**2))
depsz2(l,m)=(w2-w1)/el
& -(6.d0/(el**2)*(u1-u2)-2.d0/el*(thy1+2.d0*thy2))
& *(xl-ym*thz2)
& -(6.d0/(el**2)*(v1-v2)-2.d0/el*(thx1+2.d0*thx2))
& *(ym+xl*thz2)
& -w*(6.d0/(el**2)*(thz1-thz2)-2.d0/el*(q1+2.d0*q2))
& +1.d0/2.d0*((thy2**2)+(thx2**2)
& +((xl**2)+(ym**2))*(q2**2))
c
c 各断面要素（ファイバー要素）の軸方向応力成分を求める
dsigz(l,m)=eet(l,m)*( depsz1(l,m) + depsz2(l,m) )/2.d0
sgm(nyou,l,m)=sgm(nyou,l,m)+dsigz(l,m)
c
c
c 断面定数の断面積分
c 初期不整のある問題や座屈後など、断面の降伏が非対称に生じる問題では、
c c12,c13などの断面一次モーメントも計算する必要がある。
c
c11=c11+eet(l,m)*f
c22=c22+eet(l,m)*xl**2*f
c33=c33+eet(l,m)*ym**2*f
c44=c44+eet(l,m)* w**2 *f
end do
end do
c
haba2=haba/2.d0
taka2=taka/2.d0
c call timj(7,haba2,taka2,pi,tj)
c c55=g*tj * taka*haba**3/3.d0
c
c write(*,*) 'eet='
c write(*,*) (eet(i,i),i=1,5)
c write(*,*) 'a=',c11/eet(1,1)
c write(*,*) 'xi=',c33/eet(1,1)
c write(*,*) 'yi=',c22/eet(1,1)
c write(*,*) 'wi=',c44/eet(1,1)
c write(*,*) 'ej=', c55/g
c stop
c
return
end
c
c 初期形状が直線梁でないものや大きな回転を伴う問題を解くために
c 断面要素の軸応力を3*3の応力成分の33成分に入れて座標変換するには、
c 以下のような応力の座標変換が必要となる。これを全断面要素について
c やるには、sgm(256,200,200,3,3)のようなサイズの配列が必要になるし、
c 計算時間も増える。その割には、直線梁の座屈点の計算程度には、殆ど
c 影響しないと考えて今回は使っていない。曲がり梁でも座屈点までなら
c 軸応力成分のみを座標変換なしで扱ってもまあ大丈夫かも。
c
c 全体系の応力を局所系に座標変換 Sl=T^T Sg T
c subroutine ttsgt(t,sg,sl)
c implicit real*8(a-h,o-z)
c dimension t(14,14),sg(3,3),sl(3,3),c(3,3)
c do 10 i=1,3
c do 10 k=1,3
c c(i,k)=0.d0
c do 10 j=1,3
c c(i,k)=c(i,k)+sg(i,j)*t(j,k)
c10 continue
c do 20 i=1,3
c do 20 k=1,3
c sl(i,k)=0.d0
c do 20 j=1,3
c sl(i,k)=sl(i,k)+t(j,i)*c(j,k)
c20 continue
c return
c end
c
c 局所系の応力を全体系に座標変換 Sg=T Sl T^T
c subroutine tsltt(t,sl,sg)
c implicit real*8(a-h,o-z)
c dimension t(14,14),sl(3,3),sg(3,3),c(3,3)
c do 10 i=1,3
c do 10 k=1,3
c c(i,k)=0.d0
c do 10 j=1,3
c c(i,k)=c(i,k)+sl(i,j)*t(k,j)
c10 continue
c do 20 i=1,3
c do 20 k=1,3
c sg(i,k)=0.d0
c do 20 j=1,3
c sg(i,k)=sg(i,k)+t(i,j)*c(j,k)
c20 continue
c return
c end
c
c
c そりねじり定数
c 正規化そり関数の級数について和を取る
subroutine sori(kyuu,a,b,x,y,pi,w)
implicit real*8 (a-h, o-z)
w=0.d0
do 10 n=1,kyuu,2
w=w+(-1.d0)**( (n+1)/2 )/ n**3 *
& sinh(n*pi*y/2.d0/a) / cosh(n*pi*b/2.d0/a)
& *sin(n*pi*x/2.d0/a)
10 continue
w=32.d0*a**2/pi**3* w + x*y
return
end
c
cc ねじり定数も断面積分する場合
cc ねじり定数
cc Timoshenko が計算した正解
cc 手で積分して、級数も極限を取れる部分は取ってある
cc subroutine timj(kyuu,a,b,pi,tj)
cc implicit real*8 (a-h, o-z)
cc tj=0.d0
cc do 10 n=1,kyuu,2
cc tj=tj+1.d0/n**5*tanh(n*pi*b/2.d0/a)
cc10 continue
cc tj=1.d0-192.d0/pi**5*a/b*tj
cc return
cc end
c
c
c 梁の線形剛性行列（無次元化してある）
c 2節点14自由度
subroutine smx(el,eiy,pxz,pyz,c11,c22,c33,c44,c55,s)
implicit real*8(a-h,o-z)
dimension s(14,14)
s(3,3)= c11*el*el/eiy
s(10,3)= -c11*el*el/eiy
s(10,10)= c11*el*el/eiy
s(2,2)= 12.d0*c33/eiy /(1.d0+12.d0*pyz)
s(4,2)= -6.d0*c33/eiy /(1.d0+12.d0*pyz)
s(4,4)= (4.d0+12.d0*pyz)*c33/eiy /(1.d0+12.d0*pyz)
s(11,2)= -6.d0*c33/eiy /(1.d0+12.d0*pyz)
s(11,4)=(2.d0-12.d0*pyz)*c33/eiy /(1.d0+12.d0*pyz)
s(1,1)= 12.d0 *c22/eiy /(1.d0+12.d0*pxz)
s(5,1)= 6.d0 *c22/eiy /(1.d0+12.d0*pxz)
s(5,5)= (4.d0+12.d0*pxz)*c22/eiy /(1.d0+12.d0*pxz)
s(12,1)= 6.d0 *c22/eiy /(1.d0+12.d0*pxz)
s(12,5)=(2.d0-12.d0*pxz)*c22/eiy /(1.d0+12.d0*pxz)
s(6,6)= 12.d0*c44/eiy/el/el+c55/eiy*6.d0/5.d0
s(13,6)=-12.d0*c44/eiy/el/el-c55/eiy*6.d0/5.d0
s(7,6)= -6.d0*c44/eiy/el/el-c55/eiy/10.d0
s(14,6)= -6.d0*c44/eiy/el/el-c55/eiy/10.d0
s(7,7)= 4.d0*c44/eiy/el/el+c55/eiy*2.d0/15.d0
s(14,7)= 2.d0*c44/eiy/el/el-c55/eiy/30.d0
s(14,14)= 4.d0*c44/eiy/el/el+c55/eiy*2.d0/15.d0
s(9,2)=-s(2,2)
s(9,9)=s(2,2)
s(9,4)=-s(4,2)
s(11,9)=s(9,4)
s(11,11)=s(4,4)
s(8,1)=-s(1,1)
s(8,5)=-s(5,1)
s(8,8)=s(1,1)
s(12,8)=-s(5,1)
s(12,12)=s(5,5)
s(13,7)=-s(7,6)
s(13,13)=s(6,6)
s(14,13)=s(13,7)
c
do 100 i=1,13
do 100 j=i+1,14
s(i,j)=s(j,i)
100 continue
c
return
end
c
c 接線剛性行列（変位に関する増分式）の計算
c
SUBROUTINE STMX(nofe,is,alp,gmm,dd,d,sl,b,lg,ss,dld14)
implicit real*8 (a-h,o-z)
DIMENSION ss(14,14),sk(14,14),sl(14,14),d(14),dl(14),dd(14),
& t(14,14),ti(14,14),t14(14,14),eti(14,14),alp(256),gmm(256),
& cf1(7,7),cf2(7,7),t1(7,7),ta(7,7),tg(7,7),tp(7,7),b(15),lg(4),
& deti(14,14,14),r7x14(7,14),dr7x14(14,7,14),dss(14,14,14),
& dt0(7,14,14),dmy1(14),dmy2(14),dd3(14),dld14(14),
& tngnt1(14,14),tngnt3(14,14),tngnt4(14,14),sngnt4(14,14)
call cl4rms(14,t14, 14,sl, 14,sk, 14,ti)
call cl4rms(7,ta, 7,tg, 7,tp, 7,cf1)
call cl4rms(7,cf2, 14,eti, 7,t1, 14,t)
call clrv(14,dl)
call cl4rms(14,tngnt1,14,tngnt3,14,tngnt4,14,sngnt4)
call cl5rms(deti,r7x14,dr7x14,dss,dt0)
call clrmx(14,sngnt4)
c
c 曲がり梁の要素の初期方向を与える座標変換：ti
call trns(ti,alp(is+1),gmm(is+1),0.d0)
c 3軸回りの回転成分をオイラー角成分に変換する行列：eti
call euler(eti,d(4),d(5),d(11),d(12))
c 端部の境界条件が局所系で与えられていた時は、その変位に
c 掛けられるtiを外す（単位行列に変える）
c 左端の並進変位
if(lg(1).eq.1) then
do 1 i=1,3
do 1 j=1,3
ti(i,j)=0.d0
ti(i,i)=1.d0
1 continue
endif
c 左端の回転変位
if(lg(2).eq.1) then
do 2 i=4,6
do 2 j=4,6
ti(i,j)=0.d0
ti(i,i)=1.d0
2 continue
endif
c 右端の並進変位
if(lg(3).eq.1) then
do 3 i=8,10
do 3 j=8,10
ti(i,j)=0.d0
ti(i,i)=1.d0
3 continue
endif
c 右端の回転変位
if(lg(4).eq.1) then
do 4 i=11,13
do 4 j=11,13
ti(i,j)=0.d0
ti(i,i)=1.d0
4 continue
endif
c 局所系の変位増分ddを初期形状座標ti系に変換する
c dd3=ti^T * dd （^T は転置）
call mxtv(14,ti,dd,dd3)
c この局所系 dd3 をオイラー化する前に、各ステップの
c 増分歪を計算するための増分変位（ステップごとの1発目の弧長+NR収束）
c dld14 を計算する。
do i=1,14
dld14(i)=dld14(i)+dd3(i)
end do
c dd3 をオイラー角成分に変換
c dd=eti * dd3 = eti * ti^T * dd
call mxv(14,eti,dd3,dd)
c オイラー角成分にしてから足し合わせてトータルな d を求める
do 80 i=1,14
d(i)=d(i)+dd(i)
80 continue
c オイラー角が更新されたところで、座標変換行列を再定義:t
call trns(t,d(4),d(5),d(6))
c 3軸成分をオイラー角成分に変換する行列も再定義：eti
call euler(eti,d(4),d(5),d(11),d(12))
c 変位ベクトル
u1=d(1)
v1=d(2)
w1=d(3)
a1=d(4)
g1=d(5)
p1=d(6)
q1=d(7)
u2=d(8)
v2=d(9)
w2=d(10)
a2=d(11)
g2=d(12)
p2=d(13)
q2=d(14)
ca1=cos(a1)
cg1=cos(g1)
cp1=cos(p1)
ca2=cos(a2)
cg2=cos(g2)
cp2=cos(p2)
sa1=sin(a1)
sg1=sin(g1)
sp1=sin(p1)
sa2=sin(a2)
sg2=sin(g2)
sp2=sin(p2)
c
c 相対変位
dl(7)=q1
dl(8)=u2-u1-ca1*sg1
dl(9)=v2-v1-sa1
dl(10)=w2-w1-(ca1*cg1-1.d0)
gg=(g1+g2)/2.d0
aa=(a1+a2)/2.d0
pp=(p1+p2)/2.d0
cpp=cos(pp)
caa=cos(aa)
cgg=cos(gg)
spp=sin(pp)
saa=sin(aa)
sgg=sin(gg)
dl(11)=-cgg*(a2-a1)+sgg*caa*(p2-p1)
dl(12)=(g2-g1)+saa*(p2-p1)
dl(13)=sgg*(a2-a1)+cgg*caa*(p2-p1)
dl(14)=q2
c
c t の増分
ta(1,1)=-sp1*ca1*sg1
ta(1,2)=-cp1*ca1*sg1
ta(1,3)=-sa1*sg1
ta(2,1)=-sp1*sa1
ta(2,2)=-cp1*sa1
ta(2,3)=ca1
ta(3,1)=-sp1*ca1*cg1
ta(3,2)=-cp1*ca1*cg1
ta(3,3)=-sa1*cg1
tg(1,1)=-cp1*sg1-sp1*sa1*cg1
tg(1,2)=sp1*sg1-cp1*sa1*cg1
tg(1,3)=ca1*cg1
tg(3,1)=-cp1*cg1+sp1*sa1*sg1
tg(3,2)=sp1*cg1+cp1*sa1*sg1
tg(3,3)=-ca1*sg1
tp(1,1)=-sp1*cg1-cp1*sa1*sg1
tp(1,2)=-cp1*cg1+sp1*sa1*sg1
tp(2,1)=cp1*ca1
tp(2,2)=-sp1*ca1
tp(3,1)=sp1*sg1-cp1*sa1*cg1
tp(3,2)=cp1*sg1+sp1*sa1*cg1
do 162 i=1,3
do 162 j=1,3
ta(3+i,3+j)=ta(i,j)
tg(3+i,3+j)=tg(i,j)
tp(3+i,3+j)=tp(i,j)
162 continue
c
do 100 i=1,6
do 100 j=1,6
dt0(4,i,j)=ta(i,j)
dt0(4,i+7,j+7)=ta(i,j)
dt0(5,i,j)=tg(i,j)
dt0(5,i+7,j+7)=tg(i,j)
dt0(6,i,j)=tp(i,j)
dt0(6,i+7,j+7)=tp(i,j)
100 continue
c
c 相対変位の増分
cf1(1,1)=-1.d0
cf1(1,4)=sa1*sg1
cf1(1,5)=-ca1*cg1
cf1(2,2)=-1.d0
cf1(2,4)=-ca1
cf1(3,3)=-1.d0
cf1(3,4)=sa1*cg1
cf1(3,5)=ca1*sg1
cf1(4,4)=cgg-(p2-p1)/2.d0*sgg*saa
cf1(4,5)=(a2-a1)/2.d0*sgg+(p2-p1)/2.d0*cgg*caa
cf1(4,6)=-sgg*caa
cf1(5,4)=(p2-p1)/2.d0*caa
cf1(5,5)=-1.d0
cf1(5,6)=-saa
cf1(6,4)=-sgg-(p2-p1)/2.d0*cgg*saa
cf1(6,5)=(a2-a1)/2.d0*cgg-(p2-p1)/2.d0*sgg*caa
cf1(6,6)=-cgg*caa
cf2(1,1)=1.d0
cf2(2,2)=1.d0
cf2(3,3)=1.d0
cf2(4,4)=-cgg-(p2-p1)/2.d0*sgg*saa
cf2(4,5)=cf1(4,5)
cf2(4,6)=sgg*caa
cf2(5,4)=cf1(5,4)
cf2(5,5)=1.d0
cf2(5,6)=saa
cf2(6,4)=cf1(6,4)+2.d0*sgg
cf2(6,5)=cf1(6,5)
cf2(6,6)=cgg*caa
cf2(7,7)=1.d0
c
c その他、煩雑な増分式は、farejo の中で作る
call farejo(d,dl,t,cf1,cf2,dt0,r7x14,dr7x14,deti)
c
call mxtv8(t,dl,dmy1)
call mxv8(ss,dmy1,dmy2)
do 140 i=1,14
do 140 j=1,14
tngnt1(i,j)=0.d0
do 140 k=1,7
c dr7x17 no naka no sori no zoubun ha 0
tngnt1(i,j)=tngnt1(i,j)+dr7x14(j,k,i)*dmy2(k+7)
140 continue
c
call mxm14(ss,r7x14,sngnt4)
call mtmx14(r7x14,sngnt4,tngnt3)
c
do 260 i=1,14
do 260 j=1,14
sl(i,j)=tngnt1(i,j)+tngnt3(i,j)
260 continue
c
call mxmx(14,sl,eti,sk)
call mxtmx(14,eti,sk,sl)
c
call mtv14(r7x14,dmy2,dmy1)
c
do 300 i=1,14
do 300 j=1,14
tngnt4(i,j)=0.d0
do 300 k=1,14
tngnt4(i,j)=tngnt4(i,j)+deti(j,k,i)*dmy1(k)
300 continue
c
call mxmx(14,tngnt4,eti,sngnt4)
c
do 290 i=1,14
do 290 j=1,14
sk(i,j)=sngnt4(i,j)+sl(i,j)
290 continue
c
c stmx 内の上のddやdd3は既にtiの転置がかかった局所座標だけど、
c 最終的に解くddを全体系にするには、このskの右からもtiを掛ける
call mxmxt(14,sk,ti,sngnt4)
call mxmx(14,ti,sngnt4,sl)
c
c 不釣合力 b の計算（なぜトータルに求めてはいけないか）
c
c 載荷履歴に依存する弾塑性を解くには、トータルなdに対して
c 不釣り合い力を求めることはできない（f=K(d)d を材料非線形に対しても
c dの関数として与えることができてdの増分が取れるならできるかも
c 知れない？けど）ので、ここではΣK(Δd)Δd をbに累積しておいて、
c 行番号1200 の下辺りで、目標荷重pfからbを引いて不釣合力を求める。
c K(d)d≒ΣK(Δd)Δd が成り立つには、Δd が十分に小さい必要があり、
c そのためには、増分荷重Δfを、そのためには増分荷重Δrをじゅうぶんに
c 小さく与える必要がある。つまり、この部分でアップデートラグラジアン的
c な解き方となる。
c
c もし、不釣合力をトータルなまま計算したとすると、例えば
c f
c ↑ ／
c / →d のようなバイリニアなモデルを与えても、
c
c f
c ↑ __／
c / →d みたいな答えになってしまう。これは傾きがK1からK2に
c 変わったときに、pf-K2(d)で不釣合力を計算しているために、答えが、
c 原点を通る f=k2(d) の解に修正されてしまうからである。
c こうならないようにするには、
c pf-(K1Δd+K1Δd+K1Δd+...)+(K2Δd+K2Δd+...) のように履歴が
c 累積される形で不釣合力を求めなければならない。
c
c
c 上記のslには右から既にtiの転置が掛かってしまっているので、
c dd3 にtiを掛けて全体系にする
call mxv(14,ti,dd3,dmy1)
call mxv(14,sl,dmy1,dmy2)
c
do i=1,14
b(i)=b(i)+dmy2(i)
end do
c
c print *, 'ti'
c print '(1p7d11.3)', ( (ti(i,j),j=8,14),i=8,14 )
c print *, 'sl'
c print '(1p7d11.3)', ( (sl(i,j),j=8,14),i=8,14 )
c print *, 'dl'
c print '(1p7d11.3)', (dl(j),j=8,14)
c print *, 'b'
c print '(1p7d11.3)', (b(j),j=1,7)
c print '(1p7d11.3)', (b(j),j=8,14)
c
c
cc 弾性問題をトータルラグラジアンとして解く場合
cc call mxtv(14,t,dl,dmy2)
cc call mxv(14,ss,dmy2,dmy1)
cc call mtv14(r7x14,dmy1,dmy2)
cc call mxtv(14,eti,dmy2,dmy1)
cc call mxv(14,ti,dmy1,dmy2)
cc do 550 i=1,14
cc トータルに求める場合は+dmy2ではなく-dmy2となる
cc 不釣り合い力の初期化のところも参照
cc b(i)=b(i)-dmy2(i)
cc550 continue
c print *, 'b-total'
c print '(1p7d11.3)', (b(j),j=1,7)
c print '(1p7d11.3)', (b(j),j=8,14)
c iprint=iprint+1
c if(iprint.eq.2) stop
return
c888 format(14f5.1)
end
c
c オイラー角で座標変換行列を定義する
subroutine trns(t,al,gm,ph)
IMPLICIT REAL*8 (A-H,O-Z)
dimension t(14,14)
c write(*,*)'I am in sub trns'
t(1,1)=cos(ph)*cos(gm)-sin(ph)*sin(al)*sin(gm)
t(1,2)=-sin(ph)*cos(gm)-cos(ph)*sin(al)*sin(gm)
t(1,3)=cos(al)*sin(gm)
t(2,1)=sin(ph)*cos(al)
t(2,2)=cos(ph)*cos(al)
t(2,3)=sin(al)
t(3,1)=-cos(ph)*sin(gm)-sin(ph)*sin(al)*cos(gm)
t(3,2)=sin(ph)*sin(gm)-cos(ph)*sin(al)*cos(gm)
t(3,3)=cos(al)*cos(gm)
do 10 i=1,3
do 10 j=1,3
t(i+3,j+3)=t(i,j)
t(i+7,j+7)=t(i,j)
t(i+10,j+10)=t(i,j)
10 continue
t(7,7)=1.d0
t(14,14)=1.d0
return
end
c
c 3軸回りの回転角成分をオイラー角成分に変換する
subroutine euler(t,a1,g1,a2,g2)
IMPLICIT REAL*8 (A-H,O-Z)
dimension t(14,14)
do 10 i=1,14
t(i,i)=1.d0
10 continue
t(4,4)=-cos(g1)
t(4,6)=sin(g1)
t(5,4)=-sin(g1)*tan(a1)
t(5,6)=-cos(g1)*tan(a1)
t(6,4)=sin(g1)/cos(a1)
t(6,6)=cos(g1)/cos(a1)
t(11,11)=-cos(g2)
t(11,13)=sin(g2)
t(12,11)=-sin(g2)*tan(a2)
t(12,13)=-cos(g2)*tan(a2)
t(13,11)=sin(g2)/cos(a2)
t(13,13)=cos(g2)/cos(a2)
return
end
c
c
c 以下の mx...とかmt...という一連のサブルーチンは行列の
c かけ算を行うためのもの。
c
subroutine mtmx14(a,b,c)
implicit real*8(a-h,o-z)
dimension a(7,14),b(14,14),c(14,14)
do 10 i=1,7
do 10 k=1,7
c(i,k)=0.d0
do 10 j=1,7
c(i,k)=c(i,k)+a(j,i)*b(j+7,k)
c(i,k+7)=c(i,k+7)+a(j,i)*b(j+7,k+7)
c(i+7,k)=c(i+7,k)+a(j,i+7)*b(j+7,k)
c(i+7,k+7)=c(i+7,k+7)+a(j,i+7)*b(j+7,k+7)
10 continue
do 20 j=1,14
c(7,j)=c(7,j)+b(7,j)
20 continue
return
end
c
c
subroutine mxm14(a,b,c)
implicit real*8(a-h,o-z)
dimension a(14,14),b(7,14),c(14,14)
do 10 i=1,7
do 10 k=1,7
c(i,k)=0.d0
do 10 j=1,7
c(i,k)=c(i,k)+a(i,j+7)*b(j,k)
c(i,k+7)=c(i,k+7)+a(i,j+7)*b(j,k+7)
c(i+7,k)=c(i+7,k)+a(i+7,j+7)*b(j,k)
c(i+7,k+7)=c(i+7,k+7)+a(i+7,j+7)*b(j,k+7)
10 continue
do 20 i=1,14
c(i,7)=c(i,7)+a(i,7)
20 continue
return
end
c
subroutine mxmx(n,a,b,c)
implicit real*8(a-h,o-z)
dimension a(n,n),b(n,n),c(n,n)
do 10 i=1,n
do 10 k=1,n
c(i,k)=0.d0
do 10 j=1,n
c(i,k)=c(i,k)+a(i,j)*b(j,k)
10 continue
return
end
c
subroutine mxtmx(n,a,b,c)
implicit real*8(a-h,o-z)
dimension a(n,n),b(n,n),c(n,n)
do 10 i=1,n
do 10 k=1,n
c(i,k)=0.d0
do 10 j=1,n
c(i,k)=c(i,k)+a(j,i)*b(j,k)
10 continue
return
end
c
subroutine mxmxt(n,a,b,c)
implicit real*8(a-h,o-z)
dimension a(n,n),b(n,n),c(n,n)
do 10 i=1,n
do 10 k=1,n
c(i,k)=0.d0
do 10 j=1,n
c(i,k)=c(i,k)+a(i,j)*b(k,j)
10 continue
return
end
c
subroutine mxv(n,a,b,c)
implicit real*8(a-h,o-z)
dimension a(n,n),b(n),c(n)
do 10 i=1,n
c(i)=0.d0
do 10 j=1,n
c(i)=c(i)+a(i,j)*b(j)
10 continue
return
end
c
subroutine mxtv(n,a,b,c)
implicit real*8(a-h,o-z)
dimension a(n,n),b(n),c(n)
do 10 i=1,n
c(i)=0.d0
do 10 j=1,n
c(i)=c(i)+a(j,i)*b(j)
10 continue
return
end
c
subroutine mxv8(a,b,c)
implicit real*8(a-h,o-z)
dimension a(14,14),b(14),c(14)
do 10 i=1,14
c(i)=0.d0
do 10 j=1,7
c(i)=c(i)+a(i,j+7)*b(j+7)
10 continue
do 20 i=1,14
c(i)=c(i)+a(i,7)*b(7)
20 continue
return
end
c
subroutine mxtv8(a,b,c)
implicit real*8(a-h,o-z)
dimension a(14,14),b(14),c(14)
do 10 i=1,14
c(i)=0.d0
do 10 j=1,7
c(i)=c(i)+a(j+7,i)*b(j+7)
10 continue
do 20 i=1,14
c(i)=c(i)+a(7,i)*b(7)
20 continue
return
end
c
subroutine mtv14(a,b,c)
implicit real*8(a-h,o-z)
dimension a(7,14),b(14),c(14)
do 10 i=1,14
c(i)=0.d0
do 10 j=1,7
c(i)=c(i)+a(j,i)*b(j+7)
10 continue
c(7)=c(7)+b(7)
return
end
c
subroutine maxa(nofe,n8,imx)
dimension imx(n8+1)
c
imx(1)=1
do 10 j=2,15
imx(j)=imx(j-1)+j-1
10 continue
do 20 nj=1,nofe-1
do 20 j=9,15
imx(7*nj+j)=imx(7*nj+j-1)+j-1
20 continue
imx(n8)=imx(n8-1)+14
imx(n8+1)=imx(n8)+n8
return
end
c
c 各要素ごとの接線剛性行列を、スカイライン法の上三角行列と
c 下三角行列にわけながら全要素について重ね合わせる
subroutine skyln(nofe,i,s,so,stuh,stlh,imx)
implicit real*8(a-h,o-z)
dimension s(14,14),so(7,7),
& stuh(84*nofe+36),stlh(84*nofe+36),imx(7*nofe+9)
c
do 10 j=1,7
do 10 k=j,7
ij=7*i+j
ik=7*i+k
im=imx(ik)+ik-ij
stuh(im)=s(j,k)+so(j,k)
stlh(im)=s(k,j)+so(k,j)
so(j,k)=s(j+7,k+7)
so(k,j)=s(k+7,j+7)
10 continue
do 15 j=1,7
do 15 k=1,7
ij=7*i+j
ik=7*i+k
im7=imx(ik+7)+ik+7-ij
stuh(im7)=s(j,k+7)
stlh(im7)=s(k+7,j)
15 continue
if(i.eq.nofe-1) then
do 20 j=1,7
do 20 k=j,7
ij=7*i+j
ik=7*i+k
i7=imx(ik+7)+ik-ij
stuh(i7)=s(j+7,k+7)
stlh(i7)=s(k+7,j+7)
so(j,k)=0.d0
so(k,j)=0.d0
20 continue
endif
return
end
c
c 境界条件を入れる
subroutine bncn(nofe,bc,stuh,stlh,imx)
implicit real*8(a-h,o-z)
dimension bc(7*nofe+8),
& stuh(84*nofe+36),stlh(84*nofe+36),imx(7*nofe+9)
c
do 10 i=0,nofe-1
do 10 j=1,7
do 10 k=j,14
ij=7*i+j
ik=7*i+k
im=imx(ik)+ik-ij
stuh(im)=bc(ij)*bc(ik)*stuh(im)
stlh(im)=bc(ij)*bc(ik)*stlh(im)
10 continue
do 15 j=7*nofe+1,7*nofe+7
do 15 k=j,7*nofe+7
im=imx(k)+k-j
stuh(im)=bc(j)*bc(k)*stuh(im)
stlh(im)=bc(j)*bc(k)*stlh(im)
15 continue
n8=7*nofe+8
do 20 i=1,n8
if(bc(i).gt.1.d-6) goto 20
stuh(imx(i))=1.d0
stlh(imx(i))=0.d0
im=imx(n8)+n8-i
stuh(im)=bc(i)*stuh(im)
stlh(im)=bc(i)*stlh(im)
20 continue
return
end
c
c
c 接線剛性方程式を解いて変位増分を求めるとともに、接線剛性行列の
c 負の固有値数を求める。
c 岩熊先生が、東京大学の大計のライブラリーかどこかの「ガウスの消去法」の
c サブルーチンを2次元骨組をスカイライン法で解くために書き換えたものを
c 拝借した。
c
subroutine solve(stuh,stlh,dp,imx,msl,nnegtm,anegtm,
& det,idt,jdt,n01)
c
c 方程式は掃き出し法によりここで解かれる
c
implicit real*8 (a-h,o-z)
dimension stuh(12*n01-60),stlh(12*n01-60)
dimension imx(n01+1),msl(n01),nnegtm(n01),anegtm(n01)
c
c 行列式 = det D idt
c jdt : 負の固有値の数
c
do 100 i=1,n01
msl(i)=1+i+imx(i)-imx(i+1)
100 continue
det=stuh(1)
idt=0
jdt=0
if( det.gt.0.d0 ) go to 80
jdt=1
nnegtm(jdt)=1
anegtm(jdt)=stuh(1)
80 continue
c
c 前進消去
c
do 10 j=2,n01
ist=msl(j)+1
ied=j-1
if( ist.gt.ied ) go to 20
do 30 i=ist,ied
kst=msl(j)
if( msl(i).gt.kst ) kst=msl(i)
ked=i-1
ij=imx(j)+j-i
do 30 k=kst,ked
kj=imx(j)+j-k
ki=imx(i)+i-k
stuh(ij)=stuh(ij)-stuh(kj)*stlh(ki)
stlh(ij)=stlh(ij)-stlh(kj)*stuh(ki)
30 continue
20 continue
do 40 i=msl(j),ied
ij=imx(j)+j-i
stuh(ij)=stuh(ij)/stuh(imx(i))
stlh(ij)=stlh(ij)/stuh(imx(i))
40 continue
do 50 k=msl(j),ied
kj=imx(j)+j-k
stuh(imx(j))=stuh(imx(j))-stuh(kj)*stlh(kj)*stuh(imx(k))
stlh(imx(j))=stlh(imx(j))-stlh(kj)*stlh(imx(k))
50 continue
c
if( j.eq.n01 ) go to 10
if( stuh(imx(j)).gt.0.d0 ) go to 90
jdt=jdt+1
nnegtm(jdt)=j
anegtm(jdt)=stuh(imx(j))
90 continue
det=det*stuh(imx(j))
jjjdt=idint(dlog10(dabs(det)))
idt=idt+jjjdt
det=det/10.d0**(jjjdt)
c
10 continue
c
c 後進消去
c
do 60 j=1,n01
stlh(imx(j))=stlh(imx(j))/stuh(imx(j))
60 continue
do 70 j=2,n01
i=n01+2-j
ked=i-1
do 70 k=msl(i),ked
ii=imx(i)+i-k
stlh(imx(k))=stlh(imx(k))-stuh(ii)*stlh(imx(i))
70 continue
dp=stlh(imx(n01))
c
return
end
c
c 以下のcl..というサブルーチン群は配列の中に 0 を書き込んで
c 初期化するためのもの
c
subroutine clrv(n,a)
implicit real*8(a-h,o-z)
dimension a(n)
do 10 i=1,n
a(i)=0.d0
10 continue
return
end
c
subroutine cliv(n,l)
implicit real*8(a-h,o-z)
dimension l(n)
do 10 i=1,n
l(i)=0
10 continue
return
end
c
subroutine cl4rvs(n1,a1,n2,a2,n3,a3,n4,a4)
implicit real*8(a-h,o-z)
dimension a1(n1),a2(n2),a3(n3),a4(n4)
do 10 i=1,n1
a1(i)=0.d0
10 continue
do 20 i=1,n2
a2(i)=0.d0
20 continue
do 30 i=1,n3
a3(i)=0.d0
30 continue
do 40 i=1,n4
a4(i)=0.d0
40 continue
return
end
c
subroutine cl4ivs(n1,l1,n2,l2,n3,l3,n4,l4)
implicit real*8(a-h,o-z)
dimension l1(n1),l2(n2),l3(n3),l4(n4)
do 10 i=1,n1
l1(i)=0
10 continue
do 20 i=1,n2
l2(i)=0
20 continue
do 30 i=1,n3
l3(i)=0
30 continue
do 40 i=1,n4
l4(i)=0
40 continue
return
end
c
subroutine clrmx(n,a)
implicit real*8(a-h,o-z)
dimension a(n,n)
do 10 i=1,n
do 10 j=1,n
a(i,j)=0.d0
10 continue
return
end
c
subroutine cl4rms(n1,a1,n2,a2,n3,a3,n4,a4)
implicit real*8(a-h,o-z)
dimension a1(n1,n1),a2(n2,n2),a3(n3,n3),a4(n4,n4)
do 10 i=1,n1
do 10 j=1,n1
a1(i,j)=0.d0
10 continue
do 20 i=1,n2
do 20 j=1,n2
a2(i,j)=0.d0
20 continue
do 30 i=1,n3
do 30 j=1,n3
a3(i,j)=0.d0
30 continue
do 40 i=1,n4
do 40 j=1,n4
a4(i,j)=0.d0
40 continue
return
end
c
subroutine cl5rms(deti,r7x14,dr7x14,dss,dt0)
implicit real*8(a-h,o-z)
dimension deti(14,14,14),r7x14(7,14),dr7x14(14,7,14),
& dss(14,14,14),dt0(7,14,14)
c
do 10 i=1,14
do 10 j=1,14
do 10 k=1,14
deti(i,j,k)=0.d0
dss(i,j,k)=0.d0
10 continue
c
do 20 i=1,7
do 20 j=1,14
r7x14(i,j)=0.d0
do 20 k=1,14
dr7x14(j,i,k)=0.d0
dt0(i,j,k)=0.d0
20 continue
return
end
c
c
c
c 増分行列の幾つかがここで作られる
c 尚、このサブルーチン内の増分式は、
c 小林裕さんが Mathematica で計算してくれました。
c
subroutine farejo(d,dl,t,cf1,cf2,dt0,r7x14,dr7x14,deti)
IMPLICIT REAL*8 (A-H,O-Z)
dimension d(14),dl(14),t(14,14),cf1(7,7),cf2(7,7),
& dt0(7,14,14),r7x14(7,14),dr7x14(14,7,14),deti(14,14,14),
& a(6,3),da(14,6,3),dcf1(14,7,7),dcf2(14,7,7)
dimension tdb1(14,7,7),tdb2(14,7,7),dtb1(14,7,7),
& dtb2(14,7,7),dt1(14,14,14),t1(7,7)
do 1007 i=1,6
do 1007 j=1,3
a(i,j)=0.d0
do 1007 k=1,14
da(k,i,j)=0.d0
1007 continue
do 1013 k=1,14
do 1013 i=1,14
do 1013 j=1,14
deti(k,i,j)=0.d0
1013 continue
do 1018 i=1,7
do 1018 j=1,14
r7x14(i,j)=0.d0
do 1018 k=1,14
dr7x14(k,i,j)=0.d0
1018 continue
do 1024 k=1,14
do 1024 i=1,7
do 1024 j=1,7
dcf1(k,i,j)=0.d0
dcf2(k,i,j)=0.d0
1024 continue
do 460 k=1,7
do 460 i=1,14
do 460 j=1,14
dt1(k,i,j)=dt0(k,i,j)
dt1(7+k,i,j)=0.d0
460 continue
do 470 i=1,7
do 470 j=1,7
t1(i,j)=t(i,j)
470 continue
u1=d(1)
v1=d(2)
w1=d(3)
a1=d(4)
g1=d(5)
p1=d(6)
q1=d(7)
u2=d(8)
v2=d(9)
w2=d(10)
a2=d(11)
g2=d(12)
p2=d(13)
q2=d(14)
c
ca1=cos(a1)
cg1=cos(g1)
cp1=cos(p1)
ca2=cos(a2)
cg2=cos(g2)
cp2=cos(p2)
sa1=sin(a1)
sg1=sin(g1)
sp1=sin(p1)
sa2=sin(a2)
sg2=sin(g2)
sp2=sin(p2)
ta1=tan(a1)
ta2=tan(a2)
tg1=tan(g1)
tg2=tan(g2)
tp1=tan(p1)
tp2=tan(p2)
c
gmm=(g1+g2)/2.d0
all=(a1+a2)/2.d0
phh=(p2-p1)/2.d0
allp=(a2-a1)/2.d0
c
cgg=cos(gmm)
caa=cos(all)
cph=cos(phh)
callp=cos(allp)
sgg=sin(gmm)
saa=sin(all)
sph=sin(phh)
sallp=sin(allp)
c
ru=dl(8)
rv=dl(9)
rw=dl(10)
rx=dl(11)
ry=dl(12)
rz=dl(13)
c
dcf1(4,1,4)=ca1*sg1
dcf1(5,1,4)=sa1*cg1
dcf1(4,1,5)=sa1*cg1
dcf1(5,1,5)=ca1*sg1
dcf1(4,2,4)=sa1
dcf1(4,3,4)=ca1*cg1
dcf1(5,3,4)=-sa1*sg1
dcf1(4,3,5)=-sa1*sg1
dcf1(5,3,5)=ca1*cg1
dcf1(4,4,4)=-phh*sgg*caa/2.d0
dcf1(11,4,4)=dcf1(4,4,4)
dcf1(5,4,4)=(-sgg-phh*cgg*saa)/2.d0
dcf1(12,4,4)=dcf1(5,4,4)
dcf1(6,4,4)=sgg*saa/2.d0
dcf1(13,4,4)=-dcf1(6,4,4)
dcf1(4,4,5)=(-sgg-phh*cgg*saa)/2.d0
dcf1(5,4,5)=(allp*cgg-phh*sgg*caa)/2.d0
dcf1(6,4,5)=-cgg*caa/2.d0
dcf1(11,4,5)=(sgg-phh*cgg*saa)/2.d0
dcf1(12,4,5)=dcf1(5,4,5)
dcf1(13,4,5)=-dcf1(6,4,5)
dcf1(4,4,6)=sgg*saa/2.d0
dcf1(5,4,6)=-cgg*caa/2.d0
dcf1(11,4,6)=dcf1(4,4,6)
dcf1(12,4,6)=dcf1(5,4,6)
dcf2(4,4,4)=-phh*sgg*caa/2.d0
dcf2(5,4,4)=(sgg-phh*cgg*saa)/2.d0
dcf2(6,4,4)=sgg*saa/2.d0
dcf2(11,4,4)=dcf2(4,4,4)
dcf2(12,4,4)=dcf2(5,4,4)
dcf2(13,4,4)=-dcf2(6,4,4)
dcf1(4,5,4)=-saa*phh/2.d0
dcf1(11,5,4)=dcf1(4,5,4)
dcf1(4,5,6)=-caa/2.d0
dcf1(11,5,6)=dcf1(4,5,6)
dcf1(6,5,4)=-caa/2.d0
dcf1(13,5,4)=-dcf1(6,5,4)
dcf1(4,6,4)=-phh*cgg*caa/2.d0
dcf1(5,6,4)=(-cgg+phh*sgg*saa)/2.d0
dcf1(6,6,4)=cgg*saa/2.d0
dcf1(11,6,4)=dcf1(4,6,4)
dcf1(12,6,4)=dcf1(5,6,4)
dcf1(13,6,4)=-dcf1(6,6,4)
dcf1(4,6,5)=(-cgg+phh*sgg*saa)/2.d0
dcf1(5,6,5)=(-allp*sgg-phh*cgg*caa)/2.d0
dcf1(6,6,5)=sgg*caa/2.d0
dcf1(11,6,5)=(cgg+phh*sgg*saa)/2.d0
dcf1(12,6,5)=dcf1(5,6,5)
dcf1(13,6,5)=-dcf1(6,6,5)
dcf1(4,6,6)=cgg*saa/2.d0
dcf1(5,6,6)=sgg*caa/2.d0
dcf1(11,6,6)=dcf1(4,6,6)
dcf1(12,6,6)=dcf1(5,6,6)
dcf2(4,6,4)=-phh*cgg*caa/2.d0
dcf2(5,6,4)=(cgg+phh*sgg*saa)/2.d0
dcf2(6,6,4)=cgg*saa/2.d0
dcf2(11,6,4)=dcf2(4,6,4)
dcf2(12,6,4)=dcf2(5,6,4)
dcf2(13,6,4)=-dcf2(6,6,4)
dcf2(4,6,5)=dcf1(4,6,5)
dcf2(5,6,5)=dcf1(5,6,5)
dcf2(6,6,5)=dcf1(6,6,5)
dcf2(11,6,5)=dcf1(11,6,5)
dcf2(12,6,5)=dcf1(12,6,5)
dcf2(13,6,5)=dcf1(13,6,5)
dcf2(4,4,6)=-sgg*saa/2.d0
dcf2(5,4,6)=cgg*caa/2.d0
dcf2(11,4,6)=dcf2(4,4,6)
dcf2(12,4,6)=dcf2(5,4,6)
dcf2(4,4,5)=dcf1(4,4,5)
dcf2(5,4,5)=dcf1(5,4,5)
dcf2(6,4,5)=dcf1(6,4,5)
dcf2(11,4,5)=dcf1(11,4,5)
dcf2(12,4,5)=dcf1(12,4,5)
dcf2(13,4,5)=dcf1(13,4,5)
dcf2(4,5,4)=dcf1(4,5,4)
dcf2(6,5,4)=dcf1(6,5,4)
dcf2(11,5,4)=dcf1(11,5,4)
dcf2(13,5,4)=dcf1(13,5,4)
dcf2(4,5,6)=caa/2.d0
dcf2(11,5,6)=dcf2(4,5,6)
dcf2(4,6,6)=-cgg*saa/2.d0
dcf2(11,6,6)=-cgg*saa/2.d0
dcf2(5,6,6)=-sgg*caa/2.d0
dcf2(13,5,5)=-dcf2(6,5,5)
dcf2(12,6,6)=-sgg*caa/2.d0
c
deti(5,4,4)=sg1
deti(5,5,4)=-cg1*ta1
deti(4,5,4)=-sg1/ca1/ca1
deti(5,6,4)=cg1/ca1
deti(4,6,4)=sg1*sa1/ca1/ca1
deti(5,4,6)=cg1
deti(5,5,6)=sg1*ta1
deti(4,5,6)=-cg1/ca1/ca1
deti(5,6,6)=-sg1/ca1
deti(4,6,6)=cg1*sa1/ca1/ca1
deti(12,11,11)=sg2
deti(12,12,11)=-cg2*ta2
deti(11,12,11)=-sg2/ca2/ca2
deti(12,13,11)=cg2/ca2
deti(11,13,11)=sg2*sa2/ca2/ca2
deti(12,11,13)=cg2
deti(12,12,13)=sg2*ta2
deti(11,12,13)=-cg2/ca2/ca2
deti(12,13,13)=-sg2/ca2
deti(11,13,13)=cg2*sa2/ca2/ca2
c
c
da(1,1,1)=sg1*t(2,1)
da(2,1,1)=ta1*t(2,1)
da(3,1,1)=cg1*t(2,1)
da(4,1,1)=t(2,1)*(ru*sg1*ta1-rv+rw*cg1*ta1)
da(5,1,1)=(-ru*cg1+rw*sg1)*t(2,1)
da(6,1,1)=(-ru*sg1-rv*ta1-rw*cg1)*t(2,2)
da(8,1,1)=t(2,1)*(-sg1)
da(9,1,1)=t(2,1)*(-ta1)
da(10,1,1)=t(2,1)*(-cg1)
c
da(1,1,2)=-t(3,1)
da(3,1,2)=t(1,1)
da(4,1,2)=-t(2,1)*ru*cg1+t(2,1)*rw*sg1-cp1*sa1
da(5,1,2)=-t(1,1)*ru-t(3,1)*rw+ca1*sa1*sp1
da(6,1,2)=-t(1,2)*rw+t(3,2)*ru
da(8,1,2)=t(3,1)
da(10,1,2)=-t(1,1)
c
da(1,1,3)=-t(1,2)
da(2,1,3)=-t(2,2)
da(3,1,3)=-t(3,2)
da(4,1,3)=-t(2,2)*(sg1*ru+ta1*rv+cg1*rw)+t(1,2)*sa1*sg1
& -t(2,2)*ca1+t(3,2)*sa1*cg1
da(5,1,3)=t(3,2)*ru-t(1,2)*rw-t(1,2)*ca1*cg1+t(3,2)*ca1*sg1
da(6,1,3)=-t(1,1)*ru-t(2,1)*rv-t(3,1)*rw
c
c da(4,1,3)=-ru*t(2,1)*sg1+t(1,1)*sa1*sg1-rv*t(2,2)*ta1
c & -t(2,2)*ca1-rw*t(2,2)*cg1+t(3,2)*sa1*cg1
c da(5,1,3)=ru*t(3,1)-t(1,1)*ca1*cg1-rw*t(1,2)+t(3,2)*ca1*sg1
c da(6,1,3)=ru*t(1,2)-rv*t(2,1)-rw*t(3,1)
c
da(8,1,3)=t(1,2)
da(9,1,3)=t(2,2)
da(10,1,3)=t(3,2)
c
da(1,2,1)=t(2,2)*sg1
da(2,2,1)=t(2,2)*ta1
da(3,2,1)=t(2,2)*cg1
da(4,2,1)=t(2,2)*(ru*sg1*ta1-rv+rw*ta1*cg1)
da(5,2,1)=t(2,2)*(-ru*cg1+rw*sg1)
da(6,2,1)=t(2,1)*(ru*sg1+rv*ta1+rw*cg1)
da(8,2,1)=-t(2,2)*sg1
da(9,2,1)=-t(2,2)*ta1
da(10,2,1)=-t(2,2)*cg1
c
da(1,2,2)=-t(3,2)
da(3,2,2)=t(1,2)
da(4,2,2)=t(2,2)*(-ru*cg1+rw*sg1)+t(3,2)*sa1*sg1-t(1,2)*sa1*cg1
da(5,2,2)=-t(1,2)*ru-t(3,2)*rw-t(3,2)*ca1*cg1-t(1,2)*ca1*sg1
da(6,2,2)=-t(3,1)*ru+t(1,1)*rw
c
c da(4,2,2)=sa1*sp1-t(2,2)*cg1+t(1,3)*cp1
c da(5,2,2)=-ca1*sa1*cp1-t(2,2)*cg1+t(1,3)*cp1
c da(6,2,2)=-t(3,1)*ru-t(1,1)*rw
c
da(8,2,2)=t(3,2)
da(10,2,2)=-t(1,2)
c
do 50 i=1,3
da(i,2,3)=t(i,1)
da(7+i,2,3)=-t(i,1)
50 continue
da(4,2,3)=t(2,1)*ru*sg1+t(2,1)*rv*ta1+t(2,1)*rw*cg1
& -t(1,1)*sa1*sg1+t(2,1)*ca1-t(3,1)*sa1*cg1
da(5,2,3)=ca1*cp1+t(1,1)*rw-t(3,1)*ru
da(6,2,3)=-t(1,2)*ru-t(2,2)*rv-t(3,2)*rw
c
da(1,3,1)=sa1*sg1
da(2,3,1)=-ca1
da(3,3,1)=sa1*cg1
da(4,3,1)=-ca1*sg1*ru-sa1*rv-ca1*cg1*rw-1
c
c da(4,3,1)=sa1**2*sg1**2-ca1**2+sa1**2*cg1**2-ru*sg1*ca1-rv*sa1
c & -rw*cg1*ca1
c
da(5,3,1)=-ru*t(2,3)*cg1+rw*sg1*t(2,3)
da(8,3,1)=-t(2,3)*sg1
da(9,3,1)=ca1
da(10,3,1)=-t(2,3)*cg1
c
da(1,3,2)=-ca1*cg1
da(3,3,2)=t(1,3)
da(4,3,2)=-sa1*(ru*cg1-rw*sg1)
da(5,3,2)=-ca1**2-ca1*(ru*sg1+rw*cg1)
da(8,3,2)=ca1*cg1
da(10,3,2)=-t(1,3)
c
da(4,4,1)=-cgg*sg1*t(2,1)+phh*sgg*saa*sg1*t(2,1)
& -phh*caa*ta1*t(2,1)-ry*t(2,1)/ca1/ca1+sgg*cg1*t(2,1)
& +phh*cgg*saa*cg1*t(2,1)+(rx*sg1+ry*ta1+rz*cg1)*ta1*t(2,1)
da(11,4,1)=cgg*sg1*t(2,1)+phh*sgg*saa*sg1*t(2,1)
& -phh*caa*ta1*t(2,1)-sgg*cg1*t(2,1)+phh*cgg*saa*cg1*t(2,1)
da(5,4,1)=-allp*sgg*sg1*t(2,1)-phh*cgg*caa*sg1*t(2,1)
& -rx*cg1*t(2,1)+ta1*t(2,1)-allp*cgg*cg1*t(2,1)
& +phh*sgg*caa*cg1*t(2,1)+rz*sg1*t(2,1)
da(12,4,1)=-allp*sgg*sg1*t(2,1)-phh*cgg*caa*sg1*t(2,1)
& -ta1*t(2,1)-allp*cgg*cg1*t(2,1)+phh*sgg*caa*cg1*t(2,1)
da(6,4,1)=sgg*caa*sg1*t(2,1)+saa*ta1*t(2,1)+cgg*caa*cg1*t(2,1)
& +(-rx*sg1-ry*ta1-rz*cg1)*t(2,2)
da(13,4,1)=-sgg*caa*sg1*t(2,1)-saa*ta1*t(2,1)-cgg*caa*cg1*t(2,1)
c
da(4,4,2)=-t(2,1)*rx*cg1+t(2,1)*rz*sg1+cgg*t(3,1)
& -phh*sgg*saa*t(3,1)+phh*cgg*saa*t(1,1)+sgg*t(1,1)
da(11,4,2)=-cgg*t(3,1)-phh*sgg*saa*t(3,1)-sgg*t(1,1)
& +phh*cgg*saa*t(1,1)
da(5,4,2)=-t(1,1)*rx-t(3,1)*rz+allp*sgg*t(3,1)
& +phh*cgg*caa*t(3,1)-allp*cgg*t(1,1)+phh*sgg*caa*t(1,1)
da(12,4,2)=allp*sgg*t(3,1)+phh*cgg*caa*t(3,1)-allp*cgg*t(1,1)
& +phh*sgg*caa*t(1,1)
da(6,4,2)=t(3,2)*rx-t(1,2)*rz-sgg*caa*t(3,1)+cgg*caa*t(1,1)
da(13,4,2)=sgg*caa*t(3,1)-cgg*caa*t(1,1)
c
da(4,4,3)=-t(2,2)*rx*sg1-t(2,2)*ry*ta1-t(2,2)*rz*cg1+cgg*t(1,2)
& -phh*sgg*saa*t(1,2)+phh*caa*t(2,2)-phh*cgg*saa*t(3,2)
& -sgg*t(3,2)
da(11,4,3)=-cgg*t(1,2)-phh*sgg*saa*t(1,2)+phh*caa*t(2,2)
& +sgg*t(3,2)-phh*cgg*saa*t(3,2)
da(5,4,3)=t(3,2)*rx-t(1,2)*rz+allp*sgg*t(1,2)+phh*cgg*caa*t(1,2)
& -t(2,2)+allp*cgg*t(3,2)-phh*sgg*caa*t(3,2)
da(12,4,3)=allp*sgg*t(1,2)+phh*cgg*caa*t(1,2)+t(2,2)
& +allp*cgg*t(3,2)-phh*sgg*caa*t(3,2)
da(6,4,3)=-t(1,1)*rx-t(2,1)*ry-t(3,1)*rz-sgg*caa*t(1,2)
& -saa*t(2,2)-cgg*caa*t(3,2)
da(13,4,3)=sgg*caa*t(1,2)+saa*t(2,2)+cgg*caa*t(3,2)
c
da(4,5,1)=-cgg*sg1*t(2,2)+phh*sgg*saa*sg1*t(2,2)
& -phh*caa*ta1*t(2,2)-ry*t(2,2)/ca1/ca1-sgg*cg1*t(2,2)
& -phh*cgg*saa*cg1*t(2,2)+ta1*t(2,2)*(rx*sg1+ry*ta1+rz*cg1)
da(11,5,1)=cgg*sg1*t(2,2)+phh*sgg*saa*sg1*t(2,2)
& -phh*caa*ta1*t(2,2)-sgg*cg1*t(2,2)+phh*cgg*saa*cg1*t(2,2)
da(5,5,1)=-allp*sgg*sg1*t(2,2)-phh*cgg*caa*sg1*t(2,2)
& -rx*t(2,2)*cg1-allp*cgg*cg1*t(2,2)+ta1*t(2,2)
& +phh*sgg*caa*cg1*t(2,2)+rz*t(2,2)*sg1
da(12,5,1)=-allp*sgg*sg1*t(2,2)-phh*cgg*caa*sg1*t(2,2)
& -ta1*t(2,2)-allp*cgg*cg1*t(2,2)+phh*sgg*caa*cg1*t(2,2)
da(6,5,1)=sgg*caa*sg1*t(2,2)+saa*ta1*t(2,2)+cgg*caa*cg1*t(2,2)
& +(rx*sg1+ry*ta1+rz*cg1)*t(2,1)
da(13,5,1)=-sgg*caa*sg1*t(2,2)-saa*ta1*t(2,2)-cgg*caa*cg1*t(2,2)
c
da(4,5,2)=-cg1*rx*t(2,2)+t(2,2)*rz*sg1+cgg*t(3,2)
& -phh*sgg*saa*t(3,2)+sgg*t(1,2)+phh*cgg*saa*t(1,2)
da(11,5,2)=-cgg*t(3,2)-phh*sgg*saa*t(3,2)-sgg*t(1,2)
& +phh*cgg*saa*t(1,2)
da(5,5,2)=-t(1,2)*rx-t(3,2)*rz+allp*sgg*t(3,2)
& +phh*cgg*caa*t(3,2)-allp*cgg*t(1,2)+phh*sgg*caa*t(1,2)
da(12,5,2)=+allp*sgg*t(3,2)+phh*cgg*caa*t(3,2)-allp*cgg*t(1,2)
& +phh*sgg*caa*t(1,2)
da(6,5,2)=-t(3,1)*rx+t(1,1)*rz-sgg*caa*t(3,2)+cgg*caa*t(1,2)
da(13,5,2)=sgg*caa*t(3,2)-cgg*caa*t(1,2)
c
da(4,5,3)=t(2,1)*rx*sg1+t(2,1)*ry*ta1+t(2,1)*rz*cg1-cgg*t(1,1)
& +phh*sgg*saa*t(1,1)-phh*caa*t(2,1)+sgg*t(3,1)+phh*cgg*saa*t(3,1)
da(11,5,3)=cgg*t(1,1)+phh*sgg*saa*t(1,1)-phh*caa*t(2,1)
& -sgg*t(3,1)+phh*cgg*saa*t(3,1)
da(5,5,3)=-t(3,1)*rx+t(1,1)*rz-allp*sgg*t(1,1)
& -phh*cgg*caa*t(1,1)+t(2,1)-allp*cgg*t(3,1)+phh*sgg*caa*t(3,1)
da(12,5,3)=-allp*sgg*t(1,1)-phh*cgg*caa*t(1,1)-t(2,1)
& -allp*cgg*t(3,1)+phh*sgg*caa*t(3,1)
da(6,5,3)=-t(1,2)*rx-t(2,2)*ry-t(3,2)*rz+sgg*caa*t(1,1)
& +saa*t(2,1)+cgg*caa*t(3,1)
da(13,5,3)=-sgg*caa*t(1,1)-saa*t(2,1)-cgg*caa*t(3,1)
c
c da(4,6,1)=cgg*sg1*t(2,3)+phh*sgg*saa*sg1*t(2,3)
c & +phh*caa*ca1-ry*sa1+sgg*cg1*t(2,3)
c & +phh*cgg*saa*cg1*t(2,3)+ca1*(-rx*sg1-rz*cg1)
c da(11,6,1)=-cgg*sg1*t(2,3)+phh*sgg*saa*sg1*t(2,3)-sgg*cg1*t(2,3)
c & +phh*caa*ca1+phh*cgg*saa*cg1*t(2,3)
c da(5,6,1)=allp*sgg*sg1*t(2,3)-phh*caa*cgg*sg1*t(2,3)
c & -rx*t(2,3)*cg1-ca1-allp*cgg*cg1*t(2,3)
c & +phh*sgg*caa*cg1*t(2,3)
c da(12,6,1)=allp*sgg*sg1*t(2,3)-phh*cgg*caa*sg1*t(2,3)
c & +ca1-allp*cgg*cg1*t(2,3)+phh*sgg*caa*cg1*t(2,3)
c da(6,6,1)=sgg*caa*sg1*t(2,3)-saa*ca1
c & +cgg*caa*cg1*t(2,3)
c da(13,6,1)=-sgg*caa*sg1*t(2,3)+saa*ca1
c & -cgg*caa*cg1*t(2,3)
c
da(4,6,1)=-ca1*sg1*rx-sa1*ry-ca1*cg1*rw-cgg*sa1*sg1
& +phh*sgg*saa*sa1*sg1+phh*caa*ca1+sgg*sa1*cg1
& +phh*cgg*saa*sa1*cg1
da(11,6,1)=cgg*sa1*sg1+phh*sgg*saa*sa1*sg1+phh*caa*ca1
& -sgg*sa1*cg1+phh*cgg*saa*sa1*cg1
da(5,6,1)=-sa1*cg1*rx+sa1*sg1*rw-allp*sgg*sa1*sg1
& -phh*cgg*caa*sa1*sg1-ca1-allp*cgg*sa1*cg1+phh*saa*caa*sa1*cg1
da(12,6,1)=-allp*sgg*sa1*sg1-phh*cgg*caa*sa1*sg1
& +ca1-allp*cgg*sa1*cg1+phh*saa*caa*sa1*cg1
da(6,6,1)=sgg*caa*sa1*sg1-saa*ca1+cgg*caa*sa1*cg1
da(13,6,1)=-da(6,6,1)
c
da(4,6,2)=-rx*sa1*cg1+rz*sa1*sg1+cgg*cg1*ca1
& -phh*sgg*saa*cg1*ca1+sgg*t(1,3)+phh*cgg*saa*t(1,3)
da(11,6,2)=-cgg*cg1*ca1-phh*sgg*saa*cg1*ca1
& -sgg*t(1,3)+phh*cgg*saa*t(1,3)
da(5,6,2)=-rz*ca1*cg1-rx*ca1*sg1
& +allp*sgg*cg1*ca1+phh*cgg*caa*cg1*ca1
& -allp*cgg*t(1,3)+phh*sgg*caa*t(1,3)
da(12,6,2)=allp*sgg*cg1*ca1+phh*cgg*caa*cg1*ca1
& -allp*cgg*t(1,3)+phh*sgg*caa*t(1,3)
da(6,6,2)=-sgg*caa*cg1*ca1+cgg*caa*t(1,3)
da(13,6,2)=sgg*caa*cg1*ca1-cgg*caa*t(1,3)
c
c
da(4,5,1)=
& -((-p1 + p2)*Cos((a1 + a2)/2)*Cos(p1)*Sin(a1))/2 -
& Cos(a1)*Cos(p1)*(-g1 + g2 + (-p1 + p2)*Sin((a1 + a2)/2)) -
& Cos(a1)*Cos(g1)*Cos(p1)*(-((-p1 + p2)*Cos((g1 + g2)/2)*
& Sin((a1 + a2)/2))/2 - Sin((g1 + g2)/2)) +
& Cos(g1)*Cos(p1)*Sin(a1)*((-p1 + p2)*Cos((a1 + a2)/2)*
& Cos((g1 + g2)/2) + (-a1 + a2)*Sin((g1 + g2)/2)) +
& Cos(p1)*Sin(a1)*Sin(g1)*(-((-a1 + a2)*Cos((g1 + g2)/2)) +
& (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2)) -
& Cos(a1)*Cos(p1)*Sin(g1)*(Cos((g1 + g2)/2) -
& (-p1 + p2)*Sin((a1 + a2)/2)*Sin((g1 + g2)/2)/2)
c
da(4,6,1)=
& (-p1 + p2)*Cos(a1)*Cos((a1 + a2)/2)/2 -
& Sin(a1)*(-g1 + g2 + (-p1 + p2)*Sin((a1 + a2)/2)) -
& Cos(g1)*Sin(a1)*(-((-p1 + p2)*Cos((g1 + g2)/2)*
& Sin((a1 + a2)/2))/2 - Sin((g1 + g2)/2)) -
& Cos(a1)*Cos(g1)*((-p1 + p2)*Cos((a1 + a2)/2)*
& Cos((g1 + g2)/2) + (-a1 + a2)*Sin((g1 + g2)/2)) -
& Cos(a1)*Sin(g1)*(-((-a1 + a2)*Cos((g1 + g2)/2)) +
& (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2)) -
& Sin(a1)*Sin(g1)*(Cos((g1 + g2)/2) -
& (-p1 + p2)*Sin((a1 + a2)/2)*Sin((g1 + g2)/2)/2)
c
da(5,6,1)=
& -Cos(a1) - Sin(a1)*Sin(g1)*((-p1 + p2)*Cos((a1 + a2)/2)*
& Cos((g1 + g2)/2)/2 + (-a1 + a2)*Sin((g1 + g2)/2)/2) +
& Sin(a1)*Sin(g1)*((-p1 + p2)*Cos((a1 + a2)/2)*
& Cos((g1 + g2)/2) + (-a1 + a2)*Sin((g1 + g2)/2)) -
& Cos(g1)*Sin(a1)*((-a1 + a2)*Cos((g1 + g2)/2)/2 -
& (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2)/2) -
& Cos(g1)*Sin(a1)*(-((-a1 + a2)*Cos((g1 + g2)/2)) +
& (-p1 + p2)*Cos((a1 + a2)/2)*Sin((g1 + g2)/2))
c
c
a(1,1)=(-ru*sg1-rv*ta1-rw*cg1)*t(2,1)
a(1,2)=t(3,1)*ru-t(1,1)*rw
a(1,3)=t(1,2)*ru+t(2,2)*rv+t(3,2)*rw
a(2,1)=t(2,2)*(-ru*sg1-rv*ta1-rw*cg1)
a(2,2)=t(3,2)*ru-t(1,2)*rw
a(2,3)=-t(1,1)*ru-t(2,1)*rv-t(3,1)*rw
a(3,1)=(-ru*sg1*sa1+rv*ca1-rw*cg1*sa1)
a(3,2)=ca1*ru*cg1-t(1,3)*rw
a(3,3)=0.d0
a(4,1)=t(2,1)*(-rx*sg1-ry*ta1-rz*cg1)
a(4,2)=t(3,1)*rx-t(1,1)*rz
a(4,3)=t(1,2)*rx+t(2,2)*ry+t(3,2)*rz
a(5,1)=t(2,2)*(-rx*sg1-ry*ta1-rz*cg1)
a(5,2)=t(3,2)*rx-t(1,2)*rz
a(5,3)=-t(1,1)*rx-t(2,1)*ry-t(3,1)*rz
a(6,1)=(-rx*sg1*sa1+ry*ca1-rz*cg1*sa1)
a(6,2)=ca1*rx*cg1-t(1,3)*rz
a(6,3)=0.d0
c
do 100 i=1,3
do 100 j=1,3
r7x14(i,j)=-t(j,i)
r7x14(i,7+j)=t(j,i)
100 continue
do 110 i=1,3
r7x14(3+i,4)=a(3+i,1)+t(1,i)*cf1(4,4)+t(2,i)*cf1(5,4)
& +t(3,i)*cf1(6,4)
r7x14(i,4)=
& a(i,1)+t(1,i)*cf1(1,4)+t(2,i)*cf1(2,4)+t(3,i)*cf1(3,4)
r7x14(i,5)=a(i,2)+t(1,i)*cf1(1,5)+t(3,i)*cf1(3,5)
r7x14(3+i,5)=a(3+i,2)+t(1,i)*cf1(4,5)-t(2,i)+t(3,i)*cf1(6,5)
r7x14(i,6)=a(i,3)
r7x14(3+i,6)=
& a(3+i,3)+t(1,i)*cf1(4,6)+t(2,i)*cf1(5,6)+t(3,i)*cf1(6,6)
r7x14(3+i,11)=t(1,i)*cf2(4,4)+t(2,i)*cf1(5,4)+t(3,i)*cf2(6,4)
r7x14(3+i,12)=t(1,i)*cf2(4,5)+t(2,i)+t(3,i)*cf1(6,5)
r7x14(3+i,13)=t(1,i)*cf2(4,6)-t(2,i)*cf1(5,6)-t(3,i)*cf1(6,6)
110 continue
r7x14(7,14)=1.d0
c
do 433 k=1,14
do 433 i=1,7
do 433 l=1,7
tdb1(k,i,l)=0.d0
tdb2(k,i,l)=0.d0
dtb1(k,i,l)=0.d0
dtb2(k,i,l)=0.d0
do 433 j=1,7
tdb1(k,i,l)=tdb1(k,i,l)+t1(j,i)*dcf1(k,j,l)
tdb2(k,i,l)=tdb2(k,i,l)+t1(j,i)*dcf2(k,j,l)
dtb1(k,i,l)=dtb1(k,i,l)+dt1(k,j,i)*cf1(j,l)
dtb2(k,i,l)=dtb2(k,i,l)+dt1(k,j,i)*cf2(j,l)
433 continue
do 444 k=1,14
do 444 i=1,7
do 444 j=1,7
dr7x14(k,i,j)=tdb1(k,i,j)+dtb1(k,i,j)
dr7x14(k,i,7+j)=tdb2(k,i,j)+dtb2(k,i,j)
444 continue
do 445 k=1,14
do 445 i=1,6
do 445 j=1,3
dr7x14(k,i,j+3)=da(k,i,j)+dr7x14(k,i,j+3)
445 continue
c
return
end
c
c
c 出力した接線行列の固有値解析をせずに、最後のN-R反復で求まる
c （微少な）変位増分を固有ベクトル（座屈モード）の近似として
c 利用する場合のために、近似モードを出力する。
subroutine modeout(nofe,dd)
implicit real*8(a-h,o-z)
dimension dd(7*256+8)
open(9,file='modeu.d')
call modeo7(1,nofe,dd)
close(9)
open(9,file='modev.d')
call modeo7(2,nofe,dd)
close(9)
open(9,file='modew.d')
call modeo7(3,nofe,dd)
close(9)
open(9,file='modex.d')
call modeo7(4,nofe,dd)
close(9)
open(9,file='modey.d')
call modeo7(5,nofe,dd)
close(9)
open(9,file='modez.d')
call modeo7(6,nofe,dd)
close(9)
open(9,file='modep.d')
call modeo7(7,nofe,dd)
close(9)
return
end
c
subroutine modeo7(i7,nofe,dd)
implicit real*8(a-h,o-z)
dimension dd(7*256+8)
c do 10 i=1,nofe-1
do 10 i=0,nofe
dmode=dd(i*7+i7)
write(9,*) i+1, dmode
10 continue
return
end
c
c 変位dに外乱（座屈モードなど）を読み込む
c
subroutine modein(nofe,d,gmode)
implicit real*8(a-h,o-z)
dimension d(14*256)
open(9,file='modeu.d')
call modei7(1,nofe,d,gmode)
close(9)
open(9,file='modev.d')
call modei7(2,nofe,d,gmode)
close(9)
open(9,file='modew.d')
call modei7(3,nofe,d,gmode)
close(9)
open(9,file='modex.d')
call modei7(4,nofe,d,gmode)
close(9)
open(9,file='modey.d')
call modei7(5,nofe,d,gmode)
close(9)
open(9,file='modez.d')
call modei7(6,nofe,d,gmode)
close(9)
open(9,file='modep.d')
call modei7(7,nofe,d,gmode)
close(9)
return
end
c
c 収束時のddの微少量として近似した座屈モードは節点辺り7自由度だけど
c トータルなdはオイラー角の関係で節点辺り14自由度で覚えさせてるから
c そこを重複して書かせる。
subroutine modei7(i7,nofe,d,gmode)
implicit real*8(a-h,o-z)
dimension d(14*256)
read(9,*) idummy,dmode
i=0
d(i*14+i7)=d(i*14+i7)+gmode*dmode
do 10 i=0,nofe-2
i1=i+1
read(9,*) idummy,dmode
d(i*14+7+i7)=d(i*14+7+i7)+gmode*dmode
d(i1*14+i7)=d(i1*14+i7)+gmode*dmode
10 continue
read(9,*) idummy,dmode
n1=nofe-1
d(n1*14+7+i7)=d(n1*14+7+i7)+gmode*dmode
return
end
c
c nsgm番目の要素の軸応力sgmをgnuplot作図用にファイル出力
subroutine sgmgnu(nsgm,nhaba,ntaka,sgm)
implicit real*8(a-h,o-z)
dimension sgm(256,200,200)
open(10,file='gnusgm.out')
c 縦横が200*200とかの分割数をsplotで描くと線で塗りつぶされるので
c 適当に間引くステップ幅を nstep で指定
nstep=9
do i=1,nhaba,nstep
do j=1,ntaka,nstep
write(10,560) i,j,sgm(nsgm,i,j)
end do
write(10,*)
end do
close(10)
return
560 format( 2i4,1pd11.3)
end
c