implicit real*8(a-h,o-z) dimension x(20000),y(20000),z(20000),& &mesh(20000,6),narabi(20000,20000) !open(7,file='setten.inp') !open(8,file='youso.inp') !open(9,file='syoryou.inp') a=43.3015d-3/2.0d0!格子の縦の半分のながさ b=9.5238d-3/2.0d0!格子の横の半分のながさ h=10.0d-3!たかさ atusa=1.0d-3 p=5.0d-6 nz=4 nx=4 ny=10 ntatehan=9!縦板の数 nyokohan=2!横板の数 nyoko=ntatehan*2 ntate=nyokohan*2 pset=p/((ny+1)*ntatehan) ell=a*ntate xn=b/nx zn=a/nz yn=h/ny !節点の入力 nxa=nyoko*nx nza=ntate*nz nitimen=(nxa+1)*nyokohan+(nza+1)*ntatehan-ntatehan*nyokohan nsetten=nitimen*(ny+1) !print'("*NODE, NSET=Nall")' do j=1,ny+1 do l=1,2 do m=1,ntatehan do k=1,nz nset=k+(m-1)*nz+(l-1)*(nitimen-nz*ntatehan)+(j-1)*nitimen x(nset)=b+(m-1)*b*2.0d0 y(nset)=(j-1)*yn z(nset)=(k-1)*zn+(l-1)*(a*(ntate-1)+zn) !print'(i6,a,f12.8,a,f12.8,a,f12.8)',& !&nset,',',x(nset),',',y(nset),',',z(nset) end do end do end do end do ! !print'("*NODE, NSET=Nall")' do j=1,ny+1 do l=1,nyokohan-1 do m=1,ntatehan do k=1,nz*2-1 nset=k+(m-1)*(nz*2-1)+(j-1)*nitimen& &+(l-1)*((nz*2-1)*ntatehan+nxa+1)& &+(nz*ntatehan+nxa+1) x(nset)=b+(m-1)*b*2.0d0 y(nset)=(j-1)*yn z(nset)=k*zn+(l-1)*a*2.0d0+a !print'(i6,a,f12.8,a,f12.8,a,f12.8)',& !&nset,',',x(nset),',',y(nset),',',z(nset) end do end do end do end do ! !print'("*NODE, NSET=Nall")' do j=1,ny+1 do l=1,nyokohan do i=1,nxa+1 nset=i& &+(l-1)*(nxa+1+(nz*2-1)*ntatehan)& &+(j-1)*nitimen+nz*ntatehan x(nset)=(i-1)*xn y(nset)=(j-1)*yn z(nset)=(l-1)*a*2.0d0+a !print'(i6,a,f12.8,a,f12.8,a,f12.8)',& !&nset,',',x(nset),',',y(nset),',',z(nset) end do end do end do !要素の設定 !横板三角要素 do j=1,ny do m=1,nyokohan do i=1,nxa nyou=1+(i-1)*2+(m-1)*nxa*2*ny+(j-1)*nxa*2 mesh(nyou,1)=i+nz*ntatehan+(m-1)*(nxa+1+(nz*2-1)*ntatehan)& &+(j-1)*nitimen mesh(nyou,2)=i+1+nz*ntatehan+(m-1)*(nxa+1+(nz*2-1)*ntatehan)& &+(j-1)*nitimen mesh(nyou,3)=i+1+nitimen+nz*ntatehan+(m-1)*(nxa+1+(nz*2-1)*ntatehan)& &+(j-1)*nitimen end do end do end do !横板三角要素 do j=1,ny do m=1,nyokohan do i=1,nxa nyou=2+(i-1)*2+(m-1)*nxa*2*ny+(j-1)*nxa*2 mesh(nyou,3)=i+nz*ntatehan& &+(m-1)*(nxa+1+(nz*2-1)*ntatehan)& &+(j-1)*nitimen mesh(nyou,2)=i+nitimen+nz*ntatehan& &+(m-1)*(nxa+1+(nz*2-1)*ntatehan)& &+(j-1)*nitimen mesh(nyou,1)=i+1+nitimen+nz*ntatehan& &+(m-1)*(nxa+1+(nz*2-1)*ntatehan)& &+(j-1)*nitimen end do end do end do nyou=nxa*ny*nyokohan*2 !縦板はじ do j=1,ny do l=1,2 do m=1,ntatehan do k=1,nz-1 nyou=nyou+1 mesh(nyou,2)=k+(m-1)*nz+(l-1)*(nitimen-nz*ntatehan)& &+(j-1)*nitimen mesh(nyou+1,1)=mesh(nyou,2)+nitimen mesh(nyou,1)=mesh(nyou,2)+1 mesh(nyou+1,2)=mesh(nyou+1,1)+1 mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou+1,3)=mesh(nyou,1) nyou=nyou+1 end do end do end do end do !縦板交差1-1 do j=1,ny do m=1,ntatehan nyou=nyou+1 mesh(nyou,2)=nz+(m-1)*nz& &+(j-1)*nitimen mesh(nyou+1,1)=mesh(nyou,2)+nitimen mesh(nyou,1)=nz*ntatehan+nx+1+(m-1)*nx*2+(j-1)*nitimen mesh(nyou+1,2)=mesh(nyou,1)+nitimen mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou+1,3)=mesh(nyou,1) nyou=nyou+1 end do end do !縦板交差1-2kara do j=1,ny do l=1,nyokohan-1 do m=1,ntatehan nyou=nyou+1 mesh(nyou,2)=nz*ntatehan+nxa+1& &+nz*2-1+(l-1)*((nz*2-1)*ntatehan+nxa+1)& &+(m-1)*(nz*2-1)& &+(j-1)*nitimen mesh(nyou+1,1)=mesh(nyou,2)+nitimen mesh(nyou,1)=nz*ntatehan+nxa+1+(nz*2-1)*ntatehan& &+(m-1)*nx*2+(l-1)*((nz*2-1)*ntatehan+nxa+1)& &+nx+1+(j-1)*nitimen mesh(nyou+1,2)=mesh(nyou,1)+nitimen mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou+1,3)=mesh(nyou,1) nyou=nyou+1 end do end do end do !縦板交差2 do j=1,ny do l=1,nyokohan-1 do m=1,ntatehan nyou=nyou+1 mesh(nyou,2)=nz*ntatehan+nx+1+(m-1)*nx*2& &+(l-1)*((nz*2-1)*ntatehan+nxa+1)& &+(j-1)*nitimen mesh(nyou+1,1)=mesh(nyou,2)+nitimen mesh(nyou,1)=nz*ntatehan+nxa+1+1+(m-1)*(nz*2-1)& &+(l-1)*((nz*2-1)*ntatehan+nxa+1)+(j-1)*nitimen mesh(nyou+1,2)=mesh(nyou,1)+nitimen mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou+1,3)=mesh(nyou,1) nyou=nyou+1 end do end do end do !縦板交差2-last do j=1,ny do m=1,ntatehan nyou=nyou+1 mesh(nyou,2)=nz*ntatehan+nx+1& &+((nz*2-1)*ntatehan+nxa+1)*(nyokohan-1)+(m-1)*nx*2& &+(j-1)*nitimen mesh(nyou+1,1)=mesh(nyou,2)+nitimen mesh(nyou,1)=nz*ntatehan+nxa+1+1& &+((nz*2-1)*ntatehan+nxa+1)*(nyokohan-1)& &+(m-1)*nz& &+(j-1)*nitimen mesh(nyou+1,2)=mesh(nyou,1)+nitimen mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou+1,3)=mesh(nyou,1) nyou=nyou+1 end do end do !縦板間 do j=1,ny do l=1,nyokohan-1 do m=1,ntatehan do k=1,nz*2-2 nyou=nyou+1 mesh(nyou,2)=k+nz*ntatehan+nxa+1& &+(m-1)*(nz*2-1)+(l-1)*((nz*2-1)*ntatehan+nxa+1)& &+(j-1)*nitimen mesh(nyou+1,1)=mesh(nyou,2)+nitimen mesh(nyou,1)=k+1+nz*ntatehan+nxa+1& &+(m-1)*(nz*2-1)+(l-1)*((nz*2-1)*ntatehan+nxa+1)+(j-1)*nitimen mesh(nyou+1,2)=mesh(nyou,1)+nitimen mesh(nyou,3)=mesh(nyou,2)+nitimen mesh(nyou+1,3)=mesh(nyou,1) nyou=nyou+1 end do end do end do end do !write(7,*) nsetten !do i=1,nsetten !write(7,'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)')& !&i,',',x(i),',',y(i),',',z(i) !end do !write(8,*) nyou !do i=1,nyou !write(8,'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)')& !&i,',',& !&mesh(i,1),',',& !&mesh(i,2),',',& !&mesh(i,3) !end do !write(9,*) atusa !write(9,*) p/(2*(ny+1)*ntatehan) !write(9,*) e !close(7) !close(8) !close(9) !節点を増やしてs6要素にする iset=nsetten !narabiに-1をいれる do j=1,iset do i=1,iset narabi(i,j)=-1 end do end do !中間節点を入れて要素を増やす do i=1,nyou !1-2 if(narabi(mesh(i,1),mesh(i,2))==-1) then iset=iset+1 narabi(mesh(i,1),mesh(i,2))=iset !頂点1-2間と narabi(mesh(i,2),mesh(i,1))=iset !頂点2-1間とに節点番号isetをふり mesh(i,4)=iset x(iset)=(x(mesh(i,1))+x(mesh(i,2)))/2.d0 y(iset)=(y(mesh(i,1))+y(mesh(i,2)))/2.d0 z(iset)=(z(mesh(i,1))+z(mesh(i,2)))/2.d0 else mesh(i,4)=narabi(mesh(i,1),mesh(i,2)) end if !2-3 if(narabi(mesh(i,2),mesh(i,3))==-1) then iset=iset+1 narabi(mesh(i,2),mesh(i,3))=iset !頂点2-3間と narabi(mesh(i,3),mesh(i,2))=iset !頂点3-2間とに節点番号isetをふり mesh(i,5)=iset x(iset)=(x(mesh(i,2))+x(mesh(i,3)))/2.d0 y(iset)=(y(mesh(i,2))+y(mesh(i,3)))/2.d0 z(iset)=(z(mesh(i,2))+z(mesh(i,3)))/2.d0 else mesh(i,5)=narabi(mesh(i,2),mesh(i,3)) end if !1-3 if(narabi(mesh(i,1),mesh(i,3))==-1) then iset=iset+1 narabi(mesh(i,1),mesh(i,3))=iset !頂点1-3間と narabi(mesh(i,3),mesh(i,1))=iset !頂点3-1間とに節点番号isetをふり mesh(i,6)=iset x(iset)=(x(mesh(i,1))+x(mesh(i,3)))/2.d0 y(iset)=(y(mesh(i,1))+y(mesh(i,3)))/2.d0 z(iset)=(z(mesh(i,1))+z(mesh(i,3)))/2.d0 else mesh(i,6)=narabi(mesh(i,1),mesh(i,3)) end if end do ! print'("*NODE, NSET=Nall")' do i=1,iset print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',& &i,',',x(i),',',y(i),',',z(i) end do print'("*ELEMENT, TYPE=S6, ELSET=Eall")' do i=1,nyou print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& &i,',',& &mesh(i,1),',',& &mesh(i,2),',',& &mesh(i,3),',',& &mesh(i,4),',',& &mesh(i,5),',',& &mesh(i,6) end do print'("*NSET, NSET=kotei")' do i=1,iset if (z(i)<0.00001.and.z(i)>-0.00001) then print '(i6,a)',i,',' end if end do print'("*NSET, NSET=saika")' do i=1,iset if (z(i)ell-0.00001) then print '(i6,a)',i,',' end if end do !境界条件の設定 print'("*BOUNDARY")' print'("kotei,1,6")' !材料の設定 print'("*MATERIAL,NAME=EL")' print'("*ELASTIC")' print'("3011.7,0.3290")' !要素の設定,どの材料を使うか等 print'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL")' print'(1pe13.5)',atusa !計算のオプション print'("*STEP")' print'("*STATIC,SOLVER=SPOOLES")' !荷重 print'("*CLOAD")' print'(A,a,i6,a,1pe13.5)','saika ',',',2,',',p !出力のオプション print'("*NODE PRINT,NSET=Nall")' print'("U")' print'("*EL PRINT,ELSET=Eall")' print'("S")' print'("*NODE FILE")' print'("U")' print'("*EL FILE,POSITION=AVERAGED AT NODES")' print'("S")' print'("*END STEP")' ! end