!mentat用の3節点要素のダイヤカット円筒を作るプログラム !一枚の板から作る implicit real*8 (a-h, o-z) dimension x(2000),y(2000),z(2000),& &mesh(2000,3),mesh2(2000,3),mesh3(4000,4) pi=2.d0*asin(1.d0)!円周率 b=53d-3*pi h=104.5d-3 !高さ m=12 !周方向分割数 n=10 !高さ方向分割数 r=b/(2*m*sin(pi/m))!外半径 nyouiti=2*m !一段の要素数 th=2.d0*pi/m!一辺が作る中心角 d=(1-cos(pi/m))*r dh=sqrt((h/n)**2-d**2)!一段あたりの鉛直高さ dh2=dh*n atusa=0.8d-3 r2=r-(atusa/cos(pi/m))!内側の半径 nsettensuu=m*(n+1) nyousosuu=2*m*n !奇数段の節点座標の入力(xとy) !遇数段の節点座標の入力(xとy) !節点座標の入力(z) !節点番号と節点座標の入力と出力(奇数段) do j=0,n/2 do i=1,m nset=i+j*2*m x(nset)=r*cos((i-1)*th) y(nset)=r*sin((i-1)*th) z(nset)=dh*j*2 x(nset+nsettensuu)=r2*cos((i-1)*th) y(nset+nsettensuu)=r2*sin((i-1)*th) z(nset+nsettensuu)=dh*j*2 end do end do !節点番号と節点座標の入力と出力(偶数段) do j=0,n/2-1 do i=1,m nset=i+j*2*m+m x(nset)=r*cos((i-1)*th+th/2) y(nset)=r*sin((i-1)*th+th/2) z(nset)=dh+dh*j*2 x(nset+nsettensuu)=r2*cos((i-1)*th+th/2) y(nset+nsettensuu)=r2*sin((i-1)*th+th/2) z(nset+nsettensuu)=dh+dh*j*2 end do end do !nset=0 !do j=1,n+1 !do i=1,m !nset=nset+1 !print'(i6,a,f12.8,a,f12.8,a,f12.8)',& !&nset,',',x(nset),',',y(nset),',',z(nset) !end do !end do !要素番号とその要素の3節点の節点番号の入力と出力 !逆三角形の要素(例外)奇数 do j=1,n/2 nyou1=nyouiti+(j-1)*2*nyouiti mesh(nyou1,3)=1+m*(j-1)*2 mesh(nyou1,1)=mesh(nyou1,3)+m mesh(nyou1,2)=mesh(nyou1,1)+m-1 mesh2(nyou1,3)=1+m*(j-1)*2+nsettensuu mesh2(nyou1,1)=mesh(nyou1,3)+m+nsettensuu mesh2(nyou1,2)=mesh(nyou1,1)+m-1+nsettensuu !print'(i6,a,i6,a,i6,a,i6)',& !&nyou1,',',mesh(nyou1,1),',',mesh(nyou1,2),',',mesh(nyou1,3) end do !逆三角形の要素(例外)偶数 do j=1,n/2 nyou2=-1+j*2*nyouiti mesh(nyou2,3)=m*2+m*(j-1)*2 mesh(nyou2,1)=mesh(nyou2,3)+1 mesh(nyou2,2)=mesh(nyou2,1)+m-1 mesh2(nyou2,3)=m*2+m*(j-1)*2+nsettensuu mesh2(nyou2,1)=mesh(nyou2,3)+1+nsettensuu mesh2(nyou2,2)=mesh(nyou2,1)+m-1+nsettensuu !print'(i6,a,i6,a,i6,a,i6)',& !&nyou2,',',mesh(nyou2,1),',',mesh(nyou2,2),',',mesh(nyou2,3) end do !三角形の要素(例外)奇数 do j=1,n/2 nyou3=nyouiti-1+(j-1)*2*nyouiti mesh(nyou3,2)=1+m*(j-1)*2 mesh(nyou3,3)=mesh(nyou3,2)+m*2-1 mesh(nyou3,1)=mesh(nyou3,2)+m-1 mesh2(nyou3,2)=1+m*(j-1)*2+nsettensuu mesh2(nyou3,3)=mesh(nyou3,2)+m*2-1+nsettensuu mesh2(nyou3,1)=mesh(nyou3,2)+m-1+nsettensuu !print'(i6,a,i6,a,i6,a,i6)',& !&nyou3,',',mesh(nyou3,1),',',mesh(nyou3,2),',',mesh(nyou3,3) end do !三角形の要素(例外)偶数 do j=1,n/2 nyou4=j*2*nyouiti mesh(nyou4,2)=m+1+m*(j-1)*2 mesh(nyou4,3)=mesh(nyou4,2)+m mesh(nyou4,1)=mesh(nyou4,2)+m-1 mesh2(nyou4,2)=m+1+m*(j-1)*2+nsettensuu mesh2(nyou4,3)=mesh(nyou4,2)+m+nsettensuu mesh2(nyou4,1)=mesh(nyou4,2)+m-1+nsettensuu !print'(i6,a,i6,a,i6,a,i6)',& !&nyou4,',',mesh(nyou4,1),',',mesh(nyou4,2),',',mesh(nyou4,3) end do !三角形要素奇数 do j=1,n/2 do i=1,m-1 nyou5=-1+i*2+(j-1)*2*nyouiti mesh(nyou5,1)=i+m*(j-1)*2 mesh(nyou5,2)=mesh(nyou5,1)+1 mesh(nyou5,3)=mesh(nyou5,1)+m mesh2(nyou5,1)=i+m*(j-1)*2+nsettensuu mesh2(nyou5,2)=mesh(nyou5,1)+1+nsettensuu mesh2(nyou5,3)=mesh(nyou5,1)+m+nsettensuu !print'(i6,a,i6,a,i6,a,i6)',& !&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',mesh(nyou5,3) end do end do !三角形要素遇数 do j=1,n/2 do i=1,m-1 nyou5=nyouiti+2+(i-1)*2+(j-1)*2*nyouiti mesh(nyou5,1)=i+m*(j-1)*2+m mesh(nyou5,2)=mesh(nyou5,1)+1 mesh(nyou5,3)=mesh(nyou5,1)+m+1 mesh2(nyou5,1)=i+m*(j-1)*2+m+nsettensuu mesh2(nyou5,2)=mesh(nyou5,1)+1+nsettensuu mesh2(nyou5,3)=mesh(nyou5,1)+m+1+nsettensuu !print'(i6,a,i6,a,i6,a,i6)',& !&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',mesh(nyou5,3) end do end do !逆三角形要素奇数 do j=1,n/2 do i=1,m-1 nyou5=2+(i-1)*2+(j-1)*2*nyouiti mesh(nyou5,3)=i+m*(j-1)*2+1 mesh(nyou5,1)=mesh(nyou5,3)+m mesh(nyou5,2)=mesh(nyou5,3)+m-1 mesh2(nyou5,3)=i+m*(j-1)*2+1+nsettensuu mesh2(nyou5,1)=mesh(nyou5,3)+m+nsettensuu mesh2(nyou5,2)=mesh(nyou5,3)+m-1+nsettensuu !print'(i6,a,i6,a,i6,a,i6)',& !&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',mesh(nyou5,3) end do end do !逆三角形要素遇数 do j=1,n/2 do i=1,m-1 nyou5=nyouiti+1+(i-1)*2+(j-1)*2*nyouiti mesh(nyou5,3)=i+m*(j-1)*2+m mesh(nyou5,1)=mesh(nyou5,3)+1+m mesh(nyou5,2)=mesh(nyou5,3)+m mesh2(nyou5,3)=i+m*(j-1)*2+m+nsettensuu mesh2(nyou5,1)=mesh(nyou5,3)+1+m+nsettensuu mesh2(nyou5,2)=mesh(nyou5,3)+m+nsettensuu !print'(i6,a,i6,a,i6,a,i6)',& !&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',mesh(nyou5,3) end do end do !c3d4要素の設定 !一つ目の要素 do i=1,nyousosuu mesh3((i-1)*3+1,1)=mesh(i,3) mesh3((i-1)*3+1,3)=mesh(i,2) mesh3((i-1)*3+1,4)=mesh(i,1) mesh3((i-1)*3+1,2)=mesh2(i,1) end do !2つ目の要素 do i=1,nyousosuu mesh3((i-1)*3+2,1)=mesh(i,3) mesh3((i-1)*3+2,3)=mesh(i,2) mesh3((i-1)*3+2,4)=mesh2(i,1) mesh3((i-1)*3+2,2)=mesh2(i,2) end do !3つ目の要素 do i=1,nyousosuu mesh3((i-1)*3+3,1)=mesh(i,3) mesh3((i-1)*3+3,3)=mesh2(i,2) mesh3((i-1)*3+3,4)=mesh2(i,1) mesh3((i-1)*3+3,2)=mesh2(i,3) end do !!inpファイルを書く print'("*NODE, NSET=Nall")' do i=1,nsettensuu*2 print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',i,',',x(i),',',y(i),',',z(i) end do print'("*ELEMENT,TYPE=C3D4, ELSET=Eall")' do i=1,m*2*n*3 print'(i6,a,i6,a,i6,a,i6,a,i6)',i,',',mesh3(i,1),',',& &mesh3(i,2),',',mesh3(i,3),',',mesh3(i,4) end do print'("*NSET, NSET=kotei")' do i=1,nsettensuu*2 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,nsettensuu*2 if (z(i)dh2-0.00001) then print '(i6,a)',i,',' end if end do !境界条件の設定 print'("*BOUNDARY")' print'("kotei,1,3")' !材料の設定 print'("*MATERIAL,NAME=EL")' print'("*ELASTIC")' print'("2843.4,0.312896")' !要素の設定,どの材料を使うか等 print'("*SOLID SECTION,ELSET=Eall,MATERIAL=EL")' !計算のオプション print'("*STEP")' print'("*STATIC,SOLVER=SPOOLES")' !荷重 print'("*CLOAD")' print'(A,a,i6,a,1pe13.5)','saika ',',',3,',',-1000.0E-006/(m*2) !出力のオプション 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