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 直交異方性材料の直方体構造物を、直方体要素の c 線形／弾性の剛性行列で解くプログラム tyokum.f を c 二種類の材料で解けるようにしたtyoku2.fをgnuplot用 c の出力に対応。 c c 剛性行列は、 c Robert J. Melosh, Strucutural Analysis of Solids, c Jornal of the Structural Division, ASCE, c Vol. 89, No. ST4, 1963, pp. 205-223 c にあるものを用いた。但し、印刷ミスと思われる以下の箇所について、 c p.209: K11(1,1)の+4d'(1,1)を+4d'(1,1)に c p.209: K11(6,5)の+2d'(5,5)を+2d(5,5)に c p.209: K11(8,5)の-4d'(5,5)を-4d(5,5)に c p.209: K22(3,1)の+4d(4,4)を+d(4,4)に c p.210: K33(5,4)の-d'(6,6)を+d'(6,6)に c p.211: K31左(8,1)の2を-2に c p.211: k31左(8,4)の-1を-2に c それぞれ修正した（対応箇所についてはプログラム中にも注意を書いた） c c 05/11/17版 c tyokus.f c のサブルーチン solve をメインに入れて、配列サイズをもう少し c 大きくとれるようにした。 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 にお願いします。 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c program tyoku2ou implicit real*8 (a-h, o-z) dimension bc(9400*3),f(9400*3),e66a(6,6),e66b(6,6), & d(9400*3),ngame(100),nfair(100),ijkab(9400), & kkaisi(9400),kryou(9400),bu(9400),bv(9400),bw(9400), & lkaisi(9400),lryou(9400),fx(9400),fy(9400),fz(9400), & sm1(3,3,8,8),sm2(3,3,8,8),sm3(3,3,8,8),sm4(3,3,8,8), & n8(8),imx(9400*3), & su(45000000), sl(45000000),mksky(45000000),mkx(9400), & msl(9400*3),yokoza(100),tateza(100) c 配列の最大サイズ：nmax nmax=45000000 c c スカイラインの1次元配列のサイズは、同じ自由度数でも要素分割 c に依存するので、自由度数から一意には決められない。 c mage.f の要素分割のやり方だと、軸長（z軸）の分割数がxy断面 c の分割数に比べて多いとスパース部分が増える。 c 仮にスパース部分が皆無だったとすると c およそ自由度×自由度／２のサイズになるが、こんなサイズは想定外。 c 取り敢えず、9400を最大節点数、その３倍を最大自由度数にしておく。 c 具体的な要素分割例としては、(縦+1)×(横+1)×(軸長+1)が、 c (6+1)*(13+1)*(100+1)=9191,(20+1)*(20+1)*(20+1)=9261 など。 c mage.f で幅（半解析だから幅の半分）×高さ×軸長が c 1*(8〜20)*(50〜100)ぐらいの範囲のものを解く場合、有効数字３桁 c の収束解が得られる要素分割数は 6分割×10分割×100分割ぐらい。 c 幅（の半分）に対して桁高が10倍以上高い場合は、幅の分割数は c 1 でもいいぐらいかも。 c c 断面横方向の要素分割数：nx read(*,*) nx c 断面縦方向の要素分割数：ny read(*,*) ny c 軸方向の要素分割数：nz read(*,*) nz c 要素数：nyou nyou=nx*ny*nz c 鋼材の断面横方向の要素分割数：nxa read(*,*) nxa c 鋼材の断面縦方向の要素分割数：nya read(*,*) nya c c 要素番号 i の要素の分類番号 c yokoa*tateaの材料a：1 c yokob*tateaの材料b：2 c yokoa*tatebの材料b：3 c yokob*tatebの材料b：4 do 13 i=1,nyou read(*,*) ijkab(i) 13 continue c c 材料aの要素の断面の横の長さ:yokoa, 縦の長さ:tatea read(*,*) yokoa read(*,*) tatea c それらの半分 yokoa1=yokoa/2.d0 tatea1=tatea/2.d0 c c 材料bの要素の断面の横の長さ:yokob, 縦の長さ:tateb read(*,*) yokob read(*,*) tateb c それらの半分 yokob1=yokob/2.d0 tateb1=tateb/2.d0 c c 梁の軸長：ziku read(*,*) ziku c 一要素の軸長の半分：ziku1 ziku1=ziku/float(nz)/2.d0 c c 節点数：nset nset=(nx+1)*(ny+1)*(nz+1) c 自由度数：nziyuu nziyuu=nset*3 c c ****** 境界条件 ******* do 100 i=1, nziyuu bc(i)=1.d0 100 continue read(*,*) kkasyo do 110 i=1,kkasyo read(*,*) kkaisi(i) read(*,*) kryou(i) read(*,*) bu(i) read(*,*) bv(i) read(*,*) bw(i) do 110 j=kkaisi(i)-1, kryou(i)-1 bc(3*j+1)=bu(i) bc(3*j+2)=bv(i) bc(3*j+3)=bw(i) 110 continue c ****** 載荷条件 ******* do 120 i=1, nziyuu f(i)=0.d0 120 continue read(*,*) lkasyo write(*,*) 'lkasyo=',lkasyo do 130 i=1,lkasyo read(*,*) lkaisi(i) read(*,*) lryou(i) read(*,*) fx(i) read(*,*) fy(i) read(*,*) fz(i) do 130 j=lkaisi(i)-1, lryou(i)-1 f(3*j+1)=fx(i) f(3*j+2)=fy(i) f(3*j+3)=fz(i) 130 continue c c c 材料aの応力−ひずみ行列：e66a c まずはe66aの3*3のところに「ひずみ−応力行列」を読み込む call zeromx(6,e66a) do 140 i=1,3 do 140 j=1,3 read(*,*) e66a(i,j) 140 continue write(*,*) 'ea11=',e66a(1,1) c で、ヤング率やポアソン比を書き写してから exxa=1.d0/e66a(1,1) eyya=1.d0/e66a(2,2) ezza=1.d0/e66a(3,3) poixya=-e66a(1,2)*exxa poixza=-e66a(1,3)*exxa poiyza=-e66a(2,3)*eyya c 逆行列を取って、応力−ひずみ行列に変換する call gyaku(e66a) c せん断弾性係数を読み込む read(*,*) e66a(4,4) read(*,*) e66a(5,5) read(*,*) e66a(6,6) c 一応、せん断弾性係数も書き写しておく gxya=e66a(4,4) gxza=e66a(5,5) gyza=e66a(6,6) c c 材料bの応力−ひずみ行列：e66b c まずはe66aの3*3のところに「ひずみ−応力行列」を読み込む call zeromx(6,e66b) do 141 i=1,3 do 141 j=1,3 read(*,*) e66b(i,j) 141 continue write(*,*) 'eb11=',e66b(1,1) c で、ヤング率やポアソン比を書き写してから exxb=1.d0/e66b(1,1) eyyb=1.d0/e66b(2,2) ezzb=1.d0/e66b(3,3) poixyb=-e66b(1,2)*exxb poixzb=-e66b(1,3)*exxb poiyzb=-e66b(2,3)*eyyb c 逆行列を取って、応力−ひずみ行列に変換する call gyaku(e66b) c せん断弾性係数を読み込む read(*,*) e66b(4,4) read(*,*) e66b(5,5) read(*,*) e66b(6,6) c 一応、せん断弾性係数も書き写しておく gxyb=e66b(4,4) gxzb=e66b(5,5) gyzb=e66b(6,6) c c c 荷重増分：dp read(*,*) dp c c 出力条件： c 比較のために出力したい梁−柱理論値：nhikak c nhikak=1: w 方向の伸びを出力 c read(*,*) nhikak write(*,*) 'nhikak=',nhikak c 変位を画面出力したい節点の数：ngamen read(*,*) ngamen c 変位を画面出力したい節点の節点番号：ngame(i) do 200 i=1,ngamen read(*,*) ngame(i) 200 continue c 変位をファイル出力したい節点の数：nfairu read(*,*) nfairu c 変位をファイル出力したい節点の節点番号：nfair(i) do 210 i=1,nfairu read(*,*) nfair(i) 210 continue c c 増分ステップ数：nstep read(*,*) nstep c c ***********データ入力終わり*********** write(*,*) '材料aの応力-ひずみ行列：' do 211 i=1,3 write(*,212) (e66a(i,j),j=1,3) 211 continue write(*,*) 'Gxy,Gxz,Gyz=' write(*,212) gxya,gxza,gyza c write(*,*) '材料bの応力-ひずみ行列：' do 213 i=1,3 write(*,212) (e66b(i,j),j=1,3) 213 continue write(*,*) 'Gxy,Gxz,Gyz=' write(*,212) gxyb,gxzb,gyzb c 212 format(1p3d11.3) c c c Melesh(1963)の剛性行列をsm1-sm4に書き込む call melesh(yokoa1,tatea1,ziku1,e66a,sm1) call melesh(yokob1,tatea1,ziku1,e66b,sm2) call melesh(yokoa1,tateb1,ziku1,e66b,sm3) call melesh(yokob1,tateb1,ziku1,e66b,sm4) c c まずは、スカイラインの概形をmkskyに書き込む mkx(1)=1 do 980 i=1,nset-1 mkx(i+1)=mkx(i)+i 980 continue lasto=mkx(nset)+nset-1 write(*,*) ' 節点数(<9400)=',nset write(*,*) ' 自由度数(<28200)=',nziyuu write(*,*) ' スカイライン整形配列サイズ(<45000000)=',lasto if(lasto.gt.nmax) then write(*,*) ' スカイライン整形配列が大きすぎます' stop endif c do 990 i=1,nset do 990 j=1,nset ij=mkx(j)+j-i mksky(ij)=0 990 continue c do 1000 k=1,nz do 1000 j=1,ny do 1000 i=1,nx n8(8)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i n8(4)=n8(8)+1 n8(7)=n8(8)+(nx+1) n8(3)=n8(7)+1 n8(5)=n8(8)+(nx+1)*(ny+1) n8(1)=n8(4)+(nx+1)*(ny+1) n8(6)=n8(7)+(nx+1)*(ny+1) n8(2)=n8(3)+(nx+1)*(ny+1) c do 1000 ii=1,8 do 1000 jj=1,8 n8ii=n8(ii) n8jj=n8(jj) if(n8ii.gt.n8jj) goto 1000 ij=mkx(n8jj)+n8jj-n8ii mksky(ij)=1 1000 continue c c スカイラインの対角項番号をimxに書き込む imx(1)=1 do 1001 j=1,nset ifzumi=0 do 1001 i=1,j ij=mkx(j)+j-i if(ifzumi.eq.1) goto 1004 if(mksky(ij).eq.1) then ifzumi=1 do 1002 k=1,3 n3=3*(j-1)+k imx(n3+1)=imx(n3)+(j-i)*3+k 1002 continue endif 1004 continue if(ifzumi.eq.1) mksky(ij)=1 1001 continue c imxyz=imx(nziyuu)+nziyuu-1 write(*,*) ' スカイライン配列サイズ(<45000000)= 約',imxyz if(imxyz.ge.nmax) then write(*,*) 'たぶん配列サイズが大きすぎます' stop endif c c ************************************************************** p=0.d0 open(7,file='gnu.out') open(8,file='tanbu.out') open(9,file='ziku00.out') open(10,file='ziku10.out') open(11,file='ziku01.out') open(12,file='ziku11.out') c 指定された荷重増分ステップまで以下の作業を繰り返す do 2000 i2000=1,nstep p=p+dp c c c 全要素のsmを、それぞれの要素の8節点の節点番号ごとに c 割り振って上三角スカイラインsu,下三角スカイラインslに書き込む do 220 i=1, imxyz su(i)=0.d0 sl(i)=0.d0 220 continue c do 1003 k=1,nz do 1003 j=1,ny do 1003 i=1,nx n8(8)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i n8(4)=n8(8)+1 n8(7)=n8(8)+(nx+1) n8(3)=n8(7)+1 n8(5)=n8(8)+(nx+1)*(ny+1) n8(1)=n8(4)+(nx+1)*(ny+1) n8(6)=n8(7)+(nx+1)*(ny+1) n8(2)=n8(3)+(nx+1)*(ny+1) c do 1003 ii=1,8 do 1003 jj=1,8 n8ii=n8(ii) n8jj=n8(jj) do 1003 l=1,3 do 1003 m=1,3 i24=(n8ii-1)*3 +l j24=(n8jj-1)*3 +m if(i24.le.j24) then ij=imx(j24)+j24-i24 ijk=nx*ny*(k-1)+nx*(j-1)+i n1234=ijkab(ijk) goto (201,202,203,204), n1234 201 continue su(ij)=su(ij)+sm1(l,m,ii,jj) goto 205 202 continue su(ij)=su(ij)+sm2(l,m,ii,jj) goto 205 203 continue su(ij)=su(ij)+sm3(l,m,ii,jj) goto 205 204 continue su(ij)=su(ij)+sm4(l,m,ii,jj) 205 continue else ij=imx(i24)+i24-j24 ijk=nx*ny*(k-1)+nx*(j-1)+i n1234=ijkab(ijk) goto (301,302,303,304), n1234 301 continue sl(ij)=sl(ij)+sm1(l,m,ii,jj) goto 305 302 continue sl(ij)=sl(ij)+sm2(l,m,ii,jj) goto 305 303 continue sl(ij)=sl(ij)+sm3(l,m,ii,jj) goto 305 304 continue sl(ij)=sl(ij)+sm4(l,m,ii,jj) 305 continue endif 1003 continue c c write(*,*) '境界条件' c 境界条件を入れる c su,slの中で定義されていないs(i,j)の空白部分を定義してはいけない c 対角項以外 do 1100 i=1,nset do 1100 j=i+1,nset mij=mkx(j)+j-i if(mksky(mij).eq.1) then do 1103 l=1,3 do 1103 m=1,3 i3=(i-1)*3+l j3=(j-1)*3+m ij=imx(j3)+j3-i3 su(ij)=su(ij)*bc(i3)*bc(j3) sl(ij)=sl(ij)*bc(i3)*bc(j3) 1103 continue endif 1100 continue c 対角項 do 1101 i=1,nset do 1101 l=1,3 i3=(i-1)*3+l c 温度荷重などを全節点に与えてしまった時のため、 c 拘束節点の荷重ベクトルはここで0にする f(i3)=bc(i3)*f(i3) do 1101 m=l,3 j3=(i-1)*3+m ij=imx(j3)+j3-i3 su(ij)=su(ij)*bc(i3)*bc(j3) sl(ij)=sl(ij)*bc(i3)*bc(j3) 1101 continue c c do 1200 i=1,nziyuu if( bc(i).lt.1.d-3) then su(imx(i))=1.d0 f(i)=0.d0 end if 1200 continue c do 2010 i=1,nziyuu c 荷重項はslの対角項に入れる sl(imx(i))=f(i)*p 2010 continue c c write(*,*) '剛性方程式' c 剛性方程式を解く c tyokus.f c (http://gthmhk.virtualave.net/programoj/tyoku/tyokus.html) c では subroutine solve(stuh,stlh,imx,n01) というサブルーチン c に置いていたものを、配列サイズ節約のため、メインに入れた。 c このもとのサブルーチンは、 c 岩熊先生(http://www.civil.tohoku.ac.jp/~bear/）が、 c 東京大学の大計のライブラリーかどこかの「ガウスの消去法」の c サブルーチンを2次元骨組をスカイライン法＋弧長法で解くために c 書き換えたものを私が3次元梁用に書き換えたもの c (http://gthmhk.virtualave.net/programoj/hari/hari.html) c を弧長法なし（＋固有値の計算なし）に書き換えた。 c c c 方程式は掃き出し法によりここで解かれる c c n01=nziyuu c do 101 i=1,n01 msl(i)=1+i+imx(i)-imx(i+1) 101 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 su(ij)=su(ij)-su(kj)*sl(ki) sl(ij)=sl(ij)-sl(kj)*su(ki) 30 continue 20 continue do 40 i=msl(j),ied ij=imx(j)+j-i su(ij)=su(ij)/su(imx(i)) sl(ij)=sl(ij)/su(imx(i)) 40 continue do 10 k=msl(j),ied kj=imx(j)+j-k su(imx(j))=su(imx(j))-su(kj)*sl(kj)*su(imx(k)) sl(imx(j))=sl(imx(j))-sl(kj)*sl(imx(k)) 10 continue c c 後進消去 c do 60 j=1,n01 sl(imx(j))=sl(imx(j))/su(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 sl(imx(k))=sl(imx(k))-su(ii)*sl(imx(i)) 70 continue c return c end c c c c do i=1,nziyuu d(i)=sl(imx(i)) end do c c c ************ 出力 *************** c do 3000 j=1,ngamen nga=ngame(j)*3-3 write(*,*) '--------------------------------------------------' write(*,*) '節点番号：',ngame(j) c yoko=yokoa*nxa+yokob*(nx-nxa) tate2=tateb*(ny-nya*2) tate=tatea*nya*2+tateb*(ny-nya*2) print*,'yoko*tate=',yoko,tate a=yoko*tate aa=(yokoa*nxa)*(tatea*nya) ea=ezza*aa*2+ezzb*(a-aa*2) ga=gyza*aa*2+gyzb*(a-aa*2) print*,'gyza,aa*2,gyzb,a-aa*2=',gyza,aa*2,gyzb,a-aa*2 print*,'gyzb*a=',gyzb*a print*,'ga=',ga xi=yoko*tate**3/12.d0 xia=(yokoa*tate**3-yokoa*tate2**3)/12.d0 xib=xi-xia ezzxi=ezza*xia+ezzb*xib print*,'ezzxi=',ezzxi,ezza*xi yi=tate*yoko**3/12.d0 c c goto 599 goto (511,512,513), nhikak 511 continue c write(*,125) ' f=',p,' PL/EA=',p*ziku/ea goto 599 512 continue tawa=p*ziku**3/3.d0/ezzxi c pk=(12.d0+11.d0*poiyzb)/10.d0/(1.d0+poiyzb) c tawap=tawa+p*ziku/gyzb/a*pk pk0=1.2d0 tawap0=tawa+p*ziku/gyzb/a*pk0 tawapga=tawa+p*ziku/ga*pk0 print*,gyzb,a,ziku,ezzxi vk=gyzb*a/p/ziku*(d(nga+2)-p*ziku**3/3.d0/ezzxi) print*,d(nga+2),vk c tawav=tawa+p*ziku/gyz/a*vk c write(*,115) ' k( 等方 の k )=',1.d0/pk write(*,115) ' k(ν=0の等方性材料)=',1.d0/pk0 write(*,115) ' k(たわみv から逆算)=',1.d0/vk write(*,115) ' 荷重 f=',p write(*,115) ' たわみ v(せん断無視の梁理論)=',tawa write(*,115) ' たわみ v(ν=0で せん断考慮)=',tawap0 write(*,115) ' たわみ v(ν=0で せん断考慮)=',tawapga c write(*,115) ' たわみ v(νから せん断考慮)=',tawap t0=d(nga+2) write(*,114) 'v誤差(梁理論 %)=',(t0-tawa)/tawa*100.d0 write(*,114) 'v誤差(k=5/6 %)=',(t0-tawap0)/tawap0*100.d0 write(*,114) 'v誤差(k=5/6 %)=',(t0-tawapga)/tawapga*100.d0 c write(*,114) 'v誤差(等方k %)=',(t0-tawap)/tawap*100.d0 c 114 format(a, f10.3) write(*,*) ' (h/L)^2, 1/E" ' write(*,112) (tate/ziku)**2, 3.d0*xi*d(nga+2)/p/ziku**3 goto 599 112 format(1p2d13.5) 513 write(*,125) ' f=',p,' PL^3/3EI=',p*ziku**3/3.d0/ezz/yi goto 599 write(*,115) ' f=',p 599 continue write(*,135) ' u,v,w=',d(nga+1),d(nga+2),d(nga+3) 3000 continue c c c 自由端部の全節点変位をgnuplotの３次元プロット用に出力 yokoza(1)=0.d0 j=1 k=nz do i=1,nx ijk=nx*ny*(k-1)+nx*(j-1)+i n1234=ijkab(ijk) if(n1234.eq.1) then yokoza(i+1)=yokoa else yokoza(i+1)=yokob end if end do c do i=2,nx+1 yokoza(i)=yokoza(i)+yokoza(i-1) end do c tateza(1)=0.d0 i=1 k=nz do j=1,ny ijk=nx*ny*(k-1)+nx*(j-1)+i n1234=ijkab(ijk) if(n1234.eq.1) then tateza(j+1)=tatea else tateza(j+1)=tateb end if end do c do j=2,ny+1 tateza(j)=tateza(j)+tateza(j-1) end do c ntan0=(nx+1)*(ny+1)*nz do j=1,ny+1 i1=ntan0+(nx+1)*(j-1) do i=1,nx+1 i2=i1+i c 変位に倍率をかける時は、下の「*1.d0」を適当に変える dx=d(i2*3-2)* 1.d0 dy=d(i2*3-1)* 1.d0 dz=d(i2*3) * 1.d0 write(7,185) dx+yokoza(i),dy+tateza(j),dz end do write(7,*) end do c c 自由端部の全節点変位の出力 nfair1=nfair(1) write(8,*) '-------------------------------------------------' write(8,*) 'u' do i=1,ny+1 j1=nfair(1)+(nx+1)*(i-1) write(8,189) (d(j*3-2),j=j1,j1+nx) end do write(8,*) '-------------------------------------------------' write(8,*) 'v' do i=1,ny+1 j1=nfair(1)+(nx+1)*(i-1) write(8,189) (d(j*3-1),j=j1,j1+nx) end do write(8,*) '-------------------------------------------------' write(8,*) 'w' do i=1,ny+1 j1=nfair(1)+(nx+1)*(i-1) write(8,189) (d(j*3),j=j1,j1+nx) end do write(8,*) '-------------------------------------------------' c c 梁の4角の4辺の全節点変位をそれぞれ軸方向に出力 do k=1,nz+1 za=ziku1*2.d0*float(k-1) n00=(nx+1)*(ny+1)*(k-1)+1 n10=(nx+1)*(ny+1)*(k-1)+nx+1 n01=(nx+1)*(ny+1)*(k-1)+(nx+1)*ny+1 n11=(nx+1)*(ny+1)*k write( 9,195) za,d(n00*3-2),d(n00*3-1),d(n00*3) write(10,195) za,d(n10*3-2),d(n10*3-1),d(n10*3) write(11,195) za,d(n01*3-2),d(n01*3-1),d(n01*3) write(12,195) za,d(n11*3-2),d(n11*3-1),d(n11*3) end do c 2000 continue close(7) close(8) close(9) close(10) close(11) close(12) 115 format(a,1pd13.5) 125 format(a,1pd12.5,a,1pd12.5) 126 format(a,1pd12.5,a,1pd12.5,a,1pd12.5) 135 format(a,1p3d13.5) 145 format(1p4d13.5) 185 format(1p3d13.5) 189 format(1p9d11.3) 175 format(1p7d13.5) 195 format(1p4d13.5) c c ouryoku.fで各要素の応力を求める時のために節点変位dを c formatなしでfort.12に書き出しておく。 write(12) d rewind(12) c end c c c ********************************************************** c c ひずみ−応力行列の逆数を取って、 c 応力−ひずみ行列を求める。 c subroutine gyaku(a) implicit real*8(a-h, o-z) dimension a(6,6), b(3,3) 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 10 i=1,3 do 10 j=1,3 a(i,j)=b(j,i)/det 10 continue return end c c c 剛性行列 subroutine melesh(a,b,c,d,sm) implicit real*8(a-h, o-z) dimension d(6,6),sm(3,3,8,8), & s21a(8,8),s21b(8,8),s31a(8,8),s31b(8,8),s32a(8,8),s32b(8,8) c sm に Melesh(1963) の剛性行列そのままを与えて c メインプログラムで節点順に並び替えた剛性行列を s に入れる。 d11=d(1,1)*b**2*c**2 d22=d(2,2)*a**2*c**2 d33=d(3,3)*a**2*b**2 d44=d(4,4)*b**2*c**2 d55=d(5,5)*a**2*b**2 d66=d(6,6)*a**2*b**2 d44d=d(4,4)*a**2*c**2 d55d=d(5,5)*b**2*c**2 d66d=d(6,6)*a**2*c**2 d44dd=c*d(4,4) d55dd=b*d(5,5) d66dd=a*d(6,6) d12=c*d(1,2) d13=b*d(1,3) d23=a*d(2,3) c sm(1,1,1,1)= 4.d0*d11 +4.d0*d44d+4.d0*d55 sm(1,1,2,1)= 2.d0*d11 -4.d0*d44d+2.d0*d55 sm(1,1,3,1)= d11 -2.d0*d44d-2.d0*d55 sm(1,1,4,1)= 2.d0*d11 +2.d0*d44d-4.d0*d55 sm(1,1,5,1)=-4.d0*d11 +2.d0*d44d+2.d0*d55 sm(1,1,6,1)=-2.d0*d11 -2.d0*d44d +d55 sm(1,1,7,1)= -d11 -d44d -d55 sm(1,1,8,1)=-2.d0*d11 +d44d-2.d0*d55 c Melesh(63)ではsm(1,1,1,1)のd44dがd44になっている。 c sm(1,1,2,2)= 4.d0*d11 +4.d0*d44d+4.d0*d55 sm(1,1,3,2)= 2.d0*d11 +2.d0*d44d-4.d0*d55 sm(1,1,4,2)= d11 -2.d0*d44d-2.d0*d55 sm(1,1,5,2)=-2.d0*d11 -2.d0*d44d +d55 sm(1,1,6,2)=-4.d0*d11 +2.d0*d44d+2.d0*d55 sm(1,1,7,2)=-2.d0*d11 +d44d-2.d0*d55 sm(1,1,8,2)= -d11 -d44d -d55 c sm(1,1,3,3)= 4.d0*d11 +4.d0*d44d+4.d0*d55 sm(1,1,4,3)= 2.d0*d11 -4.d0*d44d+2.d0*d55 sm(1,1,5,3)= -d11 -d44d -d55 sm(1,1,6,3)=-2.d0*d11 +d44d-2.d0*d55 sm(1,1,7,3)=-4.d0*d11 +2.d0*d44d+2.d0*d55 sm(1,1,8,3)=-2.d0*d11 -2.d0*d44d +d55 c sm(1,1,4,4)= 4.d0*d11 +4.d0*d44d+4.d0*d55 sm(1,1,5,4)=-2.d0*d11 +d44d-2.d0*d55 sm(1,1,6,4)= -d11 -d44d -d55 sm(1,1,7,4)=-2.d0*d11 -2.d0*d44d +d55 sm(1,1,8,4)=-4.d0*d11 +2.d0*d44d+2.d0*d55 c sm(1,1,5,5)= 4.d0*d11 +4.d0*d44d+4.d0*d55 sm(1,1,6,5)= 2.d0*d11 -4.d0*d44d+2.d0*d55 sm(1,1,7,5)= d11 -2.d0*d44d-2.d0*d55 sm(1,1,8,5)= 2.d0*d11 +2.d0*d44d-4.d0*d55 c Melesh(1963)ではsm(1,1,6,5)とsm(1,1,8,5)のd55がd55dになっている c sm(1,1,6,6)= 4.d0*d11 +4.d0*d44d+4.d0*d55 sm(1,1,7,6)= 2.d0*d11 +2.d0*d44d-4.d0*d55 sm(1,1,8,6)= d11 -2.d0*d44d-2.d0*d55 c sm(1,1,7,7)= 4.d0*d11 +4.d0*d44d+4.d0*d55 sm(1,1,8,7)= 2.d0*d11 -4.d0*d44d+2.d0*d55 c sm(1,1,8,8)= 4.d0*d11 +4.d0*d44d+4.d0*d55 c c sm(2,2,1,1)= 4.d0*d22 +4.d0*d44 +4.d0*d66 sm(2,2,2,1)=-4.d0*d22 +2.d0*d44 +2.d0*d66 sm(2,2,3,1)=-2.d0*d22 +d44 -2.d0*d66 sm(2,2,4,1)= 2.d0*d22 +2.d0*d44 -4.d0*d66 sm(2,2,5,1)= 2.d0*d22 -4.d0*d44 +2.d0*d66 sm(2,2,6,1)=-2.d0*d22 -2.d0*d44 +d66 sm(2,2,7,1)= -d22 -d44 -d66 sm(2,2,8,1)= d22 -2.d0*d44 -2.d0*d66 c Melesh(63)では、sm(2,2,3,1)のd44が4.d0*d44になっている c sm(2,2,3,1)=-2.d0*d22 +4.d0*d44 -2.d0*d66 c sm(2,2,2,2)= 4.d0*d22 +4.d0*d44 +4.d0*d66 sm(2,2,3,2)= 2.d0*d22 +2.d0*d44 -4.d0*d66 sm(2,2,4,2)=-2.d0*d22 +d44 -2.d0*d66 sm(2,2,5,2)=-2.d0*d22 -2.d0*d44 +d66 sm(2,2,6,2)= 2.d0*d22 -4.d0*d44 +2.d0*d66 sm(2,2,7,2)= d22 -2.d0*d44 -2.d0*d66 sm(2,2,8,2)= -d22 -d44 -d66 c sm(2,2,3,3)= 4.d0*d22 +4.d0*d44 +4.d0*d66 sm(2,2,4,3)=-4.d0*d22 +2.d0*d44 +2.d0*d66 sm(2,2,5,3)= -d22 -d44 -d66 sm(2,2,6,3)= d22 -2.d0*d44 -2.d0*d66 sm(2,2,7,3)= 2.d0*d22 -4.d0*d44 +2.d0*d66 sm(2,2,8,3)=-2.d0*d22 -2.d0*d44 +d66 c sm(2,2,4,4)= 4.d0*d22 +4.d0*d44 +4.d0*d66 sm(2,2,5,4)= d22 -2.d0*d44 -2.d0*d66 sm(2,2,6,4)= -d22 -d44 -d66 sm(2,2,7,4)=-2.d0*d22 -2.d0*d44 +d66 sm(2,2,8,4)= 2.d0*d22 -4.d0*d44 +2.d0*d66 c sm(2,2,5,5)= 4.d0*d22 +4.d0*d44 +4.d0*d66 sm(2,2,6,5)=-4.d0*d22 +2.d0*d44 +2.d0*d66 sm(2,2,7,5)=-2.d0*d22 +d44 -2.d0*d66 sm(2,2,8,5)= 2.d0*d22 +2.d0*d44 -4.d0*d66 c sm(2,2,6,6)= 4.d0*d22 +4.d0*d44 +4.d0*d66 sm(2,2,7,6)= 2.d0*d22 +2.d0*d44 -4.d0*d66 sm(2,2,8,6)=-2.d0*d22 +d44 -2.d0*d66 c sm(2,2,7,7)= 4.d0*d22 +4.d0*d44 +4.d0*d66 sm(2,2,8,7)=-4.d0*d22 +2.d0*d44 +2.d0*d66 c sm(2,2,8,8)= 4.d0*d22 +4.d0*d44 +4.d0*d66 c c sm(3,3,1,1)= 4.d0*d33 +4.d0*d55d +4.d0*d66d sm(3,3,2,1)= 2.d0*d33 +2.d0*d55d -4.d0*d66d sm(3,3,3,1)=-2.d0*d33 +d55d -2.d0*d66d sm(3,3,4,1)=-4.d0*d33 +2.d0*d55d +2.d0*d66d sm(3,3,5,1)= 2.d0*d33 -4.d0*d55d +2.d0*d66d sm(3,3,6,1)= d33 -2.d0*d55d -2.d0*d66d sm(3,3,7,1)= -d33 -d55d -d66d sm(3,3,8,1)=-2.d0*d33 -2.d0*d55d +d66d c sm(3,3,2,2)= 4.d0*d33 +4.d0*d55d +4.d0*d66d sm(3,3,3,2)=-4.d0*d33 +2.d0*d55d +2.d0*d66d sm(3,3,4,2)=-2.d0*d33 +d55d -2.d0*d66d sm(3,3,5,2)= d33 -2.d0*d55d -2.d0*d66d sm(3,3,6,2)= 2.d0*d33 -4.d0*d55d +2.d0*d66d sm(3,3,7,2)=-2.d0*d33 -2.d0*d55d +d66d sm(3,3,8,2)= -d33 -d55d -d66d c sm(3,3,3,3)= 4.d0*d33 +4.d0*d55d +4.d0*d66d sm(3,3,4,3)= 2.d0*d33 +2.d0*d55d -4.d0*d66d sm(3,3,5,3)= -d33 -d55d -d66d sm(3,3,6,3)=-2.d0*d33 -2.d0*d55d +d66d sm(3,3,7,3)= 2.d0*d33 -4.d0*d55d +2.d0*d66d sm(3,3,8,3)= d33 -2.d0*d55d -2.d0*d66d c sm(3,3,4,4)= 4.d0*d33 +4.d0*d55d +4.d0*d66d c sm(3,3,5,4)=-2.d0*d33 -2.d0*d55d -d66d sm(3,3,5,4)=-2.d0*d33 -2.d0*d55d +d66d sm(3,3,6,4)= -d33 -d55d -d66d sm(3,3,7,4)= d33 -2.d0*d55d -2.d0*d66d sm(3,3,8,4)= 2.d0*d33 -4.d0*d55d +2.d0*d66d c Melesh(63)ではsm(3,3,5,4)の+d66dは-d66dになっているが、 c 答えが出ないので+d66dに変えた。 c sm(3,3,5,5)= 4.d0*d33 +4.d0*d55d +4.d0*d66d sm(3,3,6,5)= 2.d0*d33 +2.d0*d55d -4.d0*d66d sm(3,3,7,5)=-2.d0*d33 +d55d -2.d0*d66d sm(3,3,8,5)=-4.d0*d33 +2.d0*d55d +2.d0*d66d c sm(3,3,6,6)= 4.d0*d33 +4.d0*d55d +4.d0*d66d sm(3,3,7,6)=-4.d0*d33 +2.d0*d55d +2.d0*d66d sm(3,3,8,6)=-2.d0*d33 +d55d -2.d0*d66d c sm(3,3,7,7)= 4.d0*d33 +4.d0*d55d +4.d0*d66d sm(3,3,8,7)= 2.d0*d33 +2.d0*d55d -4.d0*d66d c sm(3,3,8,8)= 4.d0*d33 +4.d0*d55d +4.d0*d66d c c do 10 n=1,3 do 10 i=1,8 do 15 j=1,i sm(n,n,i,j)=sm(n,n,i,j)/18.d0/a/b/c 15 continue do 10 j=1,i-1 sm(n,n,j,i)=sm(n,n,i,j) 10 continue c c s21a(1,1)=-2.d0 s21a(2,1)= 2.d0 s21a(3,1)= 1.d0 s21a(4,1)=-1.d0 s21a(5,1)=-2.d0 s21a(6,1)= 2.d0 s21a(7,1)= 1.d0 s21a(8,1)=-1.d0 c s21a(1,3)=-1.d0 s21a(2,3)= 1.d0 s21a(3,3)= 2.d0 s21a(4,3)=-2.d0 s21a(5,3)=-1.d0 s21a(6,3)= 1.d0 s21a(7,3)= 2.d0 s21a(8,3)=-2.d0 c s21b(1,1)=-2.d0 s21b(2,1)=-2.d0 s21b(3,1)=-1.d0 s21b(4,1)=-1.d0 s21b(5,1)= 2.d0 s21b(6,1)= 2.d0 s21b(7,1)= 1.d0 s21b(8,1)= 1.d0 c s21b(1,3)= 1.d0 s21b(2,3)= 1.d0 s21b(3,3)= 2.d0 s21b(4,3)= 2.d0 s21b(5,3)=-1.d0 s21b(6,3)=-1.d0 s21b(7,3)=-2.d0 s21b(8,3)=-2.d0 c s31a(1,1)= 2.d0 s31a(2,1)= 1.d0 s31a(3,1)=-1.d0 s31a(4,1)=-2.d0 s31a(5,1)= 2.d0 s31a(6,1)= 1.d0 s31a(7,1)=-1.d0 s31a(8,1)=-2.d0 c Melesh(63)ではs31a(8,1)は2になっているが、i=1,8まで c 足しても0にならないので -2 の間違いだと思う。 c ということは、Melesh(63)でs31a(8,4)が -1 となっているのも c -2 の間違いだと思う。 s31a(1,2)= 1.d0 s31a(2,2)= 2.d0 s31a(3,2)=-2.d0 s31a(4,2)=-1.d0 s31a(5,2)= 1.d0 s31a(6,2)= 2.d0 s31a(7,2)=-2.d0 s31a(8,2)=-1.d0 c s31b(1,1)= 2.d0 s31b(2,1)= 1.d0 s31b(3,1)= 1.d0 s31b(4,1)= 2.d0 s31b(5,1)=-2.d0 s31b(6,1)=-1.d0 s31b(7,1)=-1.d0 s31b(8,1)=-2.d0 c s31b(1,2)= 1.d0 s31b(2,2)= 2.d0 s31b(3,2)= 2.d0 s31b(4,2)= 1.d0 s31b(5,2)=-1.d0 s31b(6,2)=-2.d0 s31b(7,2)=-2.d0 s31b(8,2)=-1.d0 c s32a(1,1)=-2.d0 s32a(2,1)=-2.d0 s32a(3,1)= 2.d0 s32a(4,1)= 2.d0 s32a(5,1)=-1.d0 s32a(6,1)=-1.d0 s32a(7,1)= 1.d0 s32a(8,1)= 1.d0 c s32a(1,5)=-1.d0 s32a(2,5)=-1.d0 s32a(3,5)= 1.d0 s32a(4,5)= 1.d0 s32a(5,5)=-2.d0 s32a(6,5)=-2.d0 s32a(7,5)= 2.d0 s32a(8,5)= 2.d0 c s32b(1,1)=-2.d0 s32b(2,1)= 2.d0 s32b(3,1)= 2.d0 s32b(4,1)=-2.d0 s32b(5,1)=-1.d0 s32b(6,1)= 1.d0 s32b(7,1)= 1.d0 s32b(8,1)=-1.d0 c s32b(1,5)=-1.d0 s32b(2,5)= 1.d0 s32b(3,5)= 1.d0 s32b(4,5)=-1.d0 s32b(5,5)=-2.d0 s32b(6,5)= 2.d0 s32b(7,5)= 2.d0 s32b(8,5)=-2.d0 c do 20 i=1,8 s21a(i,2)= s21a(i,1) s21a(i,5)=-s21a(i,1) s21a(i,6)=-s21a(i,1) s21a(i,4)= s21a(i,3) s21a(i,7)=-s21a(i,3) s21a(i,8)=-s21a(i,3) c s21b(i,2)=-s21b(i,1) s21b(i,5)= s21b(i,1) s21b(i,6)=-s21b(i,1) s21b(i,4)=-s21b(i,3) s21b(i,7)= s21b(i,3) s21b(i,8)=-s21b(i,3) c s31a(i,4)= s31a(i,1) s31a(i,5)=-s31a(i,1) s31a(i,8)=-s31a(i,1) s31a(i,3)= s31a(i,2) s31a(i,6)=-s31a(i,2) s31a(i,7)=-s31a(i,2) c s31b(i,4)=-s31b(i,1) s31b(i,5)= s31b(i,1) s31b(i,8)=-s31b(i,1) s31b(i,3)=-s31b(i,2) s31b(i,6)= s31b(i,2) s31b(i,7)=-s31b(i,2) c s32a(i,2)=-s32a(i,1) s32a(i,3)=-s32a(i,1) s32a(i,4)= s32a(i,1) s32a(i,6)=-s32a(i,5) s32a(i,7)=-s32a(i,5) s32a(i,8)= s32a(i,5) c s32b(i,2)= s32b(i,1) s32b(i,3)=-s32b(i,1) s32b(i,4)=-s32b(i,1) s32b(i,6)= s32b(i,5) s32b(i,7)=-s32b(i,5) s32b(i,8)=-s32b(i,5) 20 continue c Melesh(63)ではs31a(8,4)=-1.d0となっているが、これだと c (8,j),j=1,8まで総和したときに0にならない。 c do 30 i=1,8 do 30 j=1,8 sm(2,1,i,j)=( d12*s21a(i,j) + d44dd*s21b(i,j) )/12.d0 sm(3,1,i,j)=( d13*s31a(i,j) + d55dd*s31b(i,j) )/12.d0 sm(3,2,i,j)=( d23*s32a(i,j) + d66dd*s32b(i,j) )/12.d0 30 continue do 40 i=1,8 do 40 j=1,8 sm(1,2,i,j)=sm(2,1,j,i) sm(1,3,i,j)=sm(3,1,j,i) sm(2,3,i,j)=sm(3,2,j,i) 40 continue c return end c c subroutine zeromx(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