!6節点要素のダイヤカット円筒を作るプログラム implicit real*8 (a-h, o-z) dimension x(1000000),y(1000000),z(1000000),mesh(1000000,6)& &,mesh2(1000000,6),narabi(10000,10000) pi=2.d0*asin(1.d0)!円周率 r=0.03 !半径 h=128.5d-3 !高さ m=4!周方向分割数 n=4 !高さ方向分割数 m1=m*2 n1=n*2 ny=0 !三角形の分割の細かさ(4**nyだけ分割される) nyouiti=2*m !一段の要素数 th=2.d0*pi/m!一辺が作る中心角 dh=h/n !ダイヤ三角1段ぶんの高さ(斜めになってるから三角形の高さではない) atusa=1.000d-2 !奇数段の節点座標の入力(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 !print'(i6,a,f12.8,a,f12.8,a,f12.8)',& !&nsetki,',',x(i),',',y(i),',',z(1+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.0) y(nset)=r*sin((i-1)*th+th/2.0) z(nset)=dh+dh*j*2 !print'(i6,a,f9.6,a,f9.6,a,f9.6)',& !&nsetgu,',',x(i+m),',',y(i+m),',',z(2+j*2) end do end do !中間節点の節点座標の入力と出力(凹み部奇数段) do j=0,n/2 do i=1,m nset=i+j*(2*m+2*m*2)+m*(n+1) x(nset)=r*(cos((i-1)*th)+cos(i*th))/2.0 y(nset)=r*(sin((i-1)*th)+sin(i*th))/2.0 z(nset)=dh*j*2 !print'(i6,a,f9.6,a,f9.6,a,f9.6)',& !&nsetkiheko,',',x(i+m*2),',',y(i+m*2),',',z(2+j*2) end do end do !中間節点の節点座標の入力と出力(凹み部遇数段) do j=0,n/2-1 do i=1,m nset=i+j*(2*m+2*m*2)+m*(n+1)+m+2*m x(nset)=r*(cos((i-1)*th+th/2.0)+cos(i*th+th/2.0))/2.0 y(nset)=r*(sin((i-1)*th+th/2.0)+sin(i*th+th/2.0))/2.0 z(nset)=dh*j*2+dh !print'(i6,a,f9.6,a,f9.6,a,f9.6)',& !&nsetguheko,',',x(i+m*3),',',y(i+m*3),',',z(2+j*2) end do end do !中間節点の節点座標の入力と出力(段の中間節点) r1=(r*(cos(0.)+cos(th/2.0))/2.0)/cos(th/4.0)!中間節点の半径 do j=0,n-1 do i=1,m*2 nset=i+j*(2*m+m)+m*(n+1)+m x(nset)=r1*(cos((i-1)*th/2.0+th/4.0)) y(nset)=r1*(sin((i-1)*th/2.0+th/4.0)) z(nset)=dh*j+dh/2.0 !print'(i6,a,f9.6,a,f9.6,a,f9.6)',& !&nsetguheko,',',x(i+m*3),',',y(i+m*3),',',z(2+j*2) end do end do !三角形の要素を更に分割するための節点番号に-1を入れる nset=m*2*(n*2+1) do i=1,nset do j=1,nset narabi(i,j)=-1 end do end do !要素番号とその要素の3節点の節点番号の入力と出力 !逆三角形の要素(例外)奇数 !print'("*ELEMENT, TYPE=S6, ELSET=Eall")' 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 mesh(nyou1,6)=m*(n+1)+m+1+(m+m*2)*(j-1)*2 mesh(nyou1,5)=m*(n+1)+m+m*2+(m+m*2)*(j-1)*2 mesh(nyou1,4)=m*(n+1)+m+m*2+m+(m+m*2)*(j-1)*2 !print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& !&nyou1,',',mesh(nyou1,1),',',mesh(nyou1,2),',',& !&mesh(nyou1,3),',',mesh(nyou1,4),',',mesh(nyou1,5)& !&,',',mesh(nyou1,6) 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 mesh(nyou2,5)=m*(n+1)+m+m*2+m+m*2-1+(m+m*2)*(j-1)*2 mesh(nyou2,6)=mesh(nyou2,5)+1 mesh(nyou2,4)=mesh(nyou2,6)+m !print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& !&nyou2,',',mesh(nyou2,1),',',mesh(nyou2,2),',',& !&mesh(nyou2,3),',',mesh(nyou2,4),',',mesh(nyou2,5)& !&,',',mesh(nyou2,6) 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 mesh(nyou3,4)=m*(n+1)+m+(m+m*2)*(j-1)*2 mesh(nyou3,6)=m*(n+1)+m+m*2-1+(m+m*2)*(j-1)*2 mesh(nyou3,5)=mesh(nyou3,6)+1 !print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& !&nyou3,',',mesh(nyou3,1),',',mesh(nyou3,2),',',& !&mesh(nyou3,3),',',mesh(nyou3,4),',',mesh(nyou3,5)& !&,',',mesh(nyou3,6) 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 mesh(nyou4,4)=m*(n+1)+m+m*2+m+(m+m*2)*(j-1)*2 mesh(nyou4,5)=m*(n+1)+m+m*2+m+1+(m+m*2)*(j-1)*2 mesh(nyou4,6)=m*(n+1)+m+m*2+m+m*2+(m+m*2)*(j-1)*2 !print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& !&nyou4,',',mesh(nyou4,1),',',mesh(nyou4,2),',',& !&mesh(nyou4,3),',',mesh(nyou4,4),',',mesh(nyou4,5)& !&,',',mesh(nyou4,6) 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 mesh(nyou5,4)=i+m*(n+1)+(m*2+m*2*2)*(j-1) mesh(nyou5,6)=i*2-1+m*(n+1)+m+(m+m*2)*(j-1)*2 mesh(nyou5,5)=mesh(nyou5,6)+1 !print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& !&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',& !&mesh(nyou5,3),',',mesh(nyou5,4),',',mesh(nyou5,5)& !&,',',mesh(nyou5,6) 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 mesh(nyou5,4)=i+m*(n+1)+m+m*2+(m+m*2)*(j-1)*2 mesh(nyou5,6)=i*2+m*(n+1)+m+m*2+m+(m+m*2)*(j-1)*2 mesh(nyou5,5)=mesh(nyou5,6)+1 !print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& !&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',& !&mesh(nyou5,3),',',mesh(nyou5,4),',',mesh(nyou5,5)& !&,',',mesh(nyou5,6) 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 mesh(nyou5,5)=i*2+m*(n+1)+m+(m+m*2)*(j-1)*2 mesh(nyou5,6)=mesh(nyou5,5)+1 mesh(nyou5,4)=i+m*(n+1)+m+m*2+(m+m*2)*(j-1)*2 !print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& !&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',& !&mesh(nyou5,3),',',mesh(nyou5,4),',',mesh(nyou5,5)& !&,',',mesh(nyou5,6) 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 mesh(nyou5,5)=i*2-1+m*(n+1)+m+m*2+m+(m+m*2)*(j-1)*2 mesh(nyou5,6)=mesh(nyou5,5)+1 mesh(nyou5,4)=i+m*(n+1)+m+m*2+m+m*2+(m+m*2)*(j-1)*2 !print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& !&nyou5,',',mesh(nyou5,1),',',mesh(nyou5,2),',',& !&mesh(nyou5,3),',',mesh(nyou5,4),',',mesh(nyou5,5)& !&,',',mesh(nyou5,6) end do end do !増やす時ためのmesh2(i,j)の設定 iset=nset nyou=2*m*n !mesh2にmeshを入れる do i=1,nyou do j=1,6 mesh2(i,j)=mesh(i,j) end do end do ! do if (ny==0) exit ! do i=1,nyou mesh2(1+(i-1)*4,1)=mesh(i,1) mesh2(1+(i-1)*4,2)=mesh(i,4) mesh2(1+(i-1)*4,3)=mesh(i,6) mesh2(2+(i-1)*4,1)=mesh(i,4) mesh2(2+(i-1)*4,2)=mesh(i,2) mesh2(2+(i-1)*4,3)=mesh(i,5) mesh2(3+(i-1)*4,1)=mesh(i,6) mesh2(3+(i-1)*4,2)=mesh(i,5) mesh2(3+(i-1)*4,3)=mesh(i,3) mesh2(4+(i-1)*4,1)=mesh(i,5) mesh2(4+(i-1)*4,2)=mesh(i,6) mesh2(4+(i-1)*4,3)=mesh(i,4) end do !中間節点を入れて要素を増やす do i=1,nyou !1-4 if(narabi(mesh(i,1),mesh(i,4))==-1) then iset=iset+1 narabi(mesh(i,1),mesh(i,4))=iset !頂点1-4間と narabi(mesh(i,4),mesh(i,1))=iset !頂点4-1間とに節点番号isetをふり mesh2(1+(i-1)*4,4)=iset x(iset)=(x(mesh(i,1))+x(mesh(i,4)))/2.d0 y(iset)=(y(mesh(i,1))+y(mesh(i,4)))/2.d0 z(iset)=(z(mesh(i,1))+z(mesh(i,4)))/2.d0 else mesh2(1+(i-1)*4,4)=narabi(mesh(i,1),mesh(i,4)) end if !4-6 if(narabi(mesh(i,4),mesh(i,6))==-1) then iset=iset+1 narabi(mesh(i,4),mesh(i,6))=iset !頂点1-4間と narabi(mesh(i,6),mesh(i,4))=iset !頂点4-1間とに節点番号isetをふり mesh2(1+(i-1)*4,5)=iset mesh2(4+(i-1)*4,5)=iset x(iset)=(x(mesh(i,4))+x(mesh(i,6)))/2.d0 y(iset)=(y(mesh(i,4))+y(mesh(i,6)))/2.d0 z(iset)=(z(mesh(i,4))+z(mesh(i,6)))/2.d0 else mesh2(1+(i-1)*4,5)=narabi(mesh(i,4),mesh(i,6)) mesh2(4+(i-1)*4,5)=narabi(mesh(i,4),mesh(i,6)) end if !1-6 if(narabi(mesh(i,1),mesh(i,6))==-1) then iset=iset+1 narabi(mesh(i,1),mesh(i,6))=iset !頂点1-6間と narabi(mesh(i,6),mesh(i,1))=iset !頂点6-1間とに節点番号isetをふり mesh2(1+(i-1)*4,6)=iset x(iset)=(x(mesh(i,1))+x(mesh(i,6)))/2.d0 y(iset)=(y(mesh(i,1))+y(mesh(i,6)))/2.d0 z(iset)=(z(mesh(i,1))+z(mesh(i,6)))/2.d0 else mesh2(1+(i-1)*4,6)=narabi(mesh(i,1),mesh(i,6)) end if !2-4 if(narabi(mesh(i,2),mesh(i,4))==-1) then iset=iset+1 narabi(mesh(i,2),mesh(i,4))=iset !頂点1-4間と narabi(mesh(i,4),mesh(i,2))=iset !頂点4-1間とに節点番号isetをふり mesh2(2+(i-1)*4,4)=iset x(iset)=(x(mesh(i,2))+x(mesh(i,4)))/2.d0 y(iset)=(y(mesh(i,2))+y(mesh(i,4)))/2.d0 z(iset)=(z(mesh(i,2))+z(mesh(i,4)))/2.d0 else mesh2(2+(i-1)*4,4)=narabi(mesh(i,2),mesh(i,4)) end if !2-5 if(narabi(mesh(i,2),mesh(i,5))==-1) then iset=iset+1 narabi(mesh(i,2),mesh(i,5))=iset !頂点1-4間と narabi(mesh(i,5),mesh(i,2))=iset !頂点4-1間とに節点番号isetをふり mesh2(2+(i-1)*4,5)=iset x(iset)=(x(mesh(i,2))+x(mesh(i,5)))/2.d0 y(iset)=(y(mesh(i,2))+y(mesh(i,5)))/2.d0 z(iset)=(z(mesh(i,2))+z(mesh(i,5)))/2.d0 else mesh2(2+(i-1)*4,5)=narabi(mesh(i,2),mesh(i,5)) end if !4-5 if(narabi(mesh(i,4),mesh(i,5))==-1) then iset=iset+1 narabi(mesh(i,4),mesh(i,5))=iset !頂点1-4間と narabi(mesh(i,5),mesh(i,4))=iset !頂点4-1間とに節点番号isetをふり mesh2(2+(i-1)*4,6)=iset mesh2(4+(i-1)*4,6)=iset x(iset)=(x(mesh(i,4))+x(mesh(i,5)))/2.d0 y(iset)=(y(mesh(i,4))+y(mesh(i,5)))/2.d0 z(iset)=(z(mesh(i,4))+z(mesh(i,5)))/2.d0 else mesh2(2+(i-1)*4,6)=narabi(mesh(i,4),mesh(i,5)) mesh2(4+(i-1)*4,6)=narabi(mesh(i,4),mesh(i,5)) end if !5-6 if(narabi(mesh(i,5),mesh(i,6))==-1) then iset=iset+1 narabi(mesh(i,5),mesh(i,6))=iset !頂点1-4間と narabi(mesh(i,6),mesh(i,5))=iset !頂点4-1間とに節点番号isetをふり mesh2(3+(i-1)*4,4)=iset mesh2(4+(i-1)*4,4)=iset x(iset)=(x(mesh(i,5))+x(mesh(i,6)))/2.d0 y(iset)=(y(mesh(i,5))+y(mesh(i,6)))/2.d0 z(iset)=(z(mesh(i,5))+z(mesh(i,6)))/2.d0 else mesh2(3+(i-1)*4,4)=narabi(mesh(i,5),mesh(i,6)) mesh2(4+(i-1)*4,4)=narabi(mesh(i,5),mesh(i,6)) end if !3-5 if(narabi(mesh(i,3),mesh(i,5))==-1) then iset=iset+1 narabi(mesh(i,3),mesh(i,5))=iset !頂点1-4間と narabi(mesh(i,5),mesh(i,3))=iset !頂点4-1間とに節点番号isetをふり mesh2(3+(i-1)*4,5)=iset x(iset)=(x(mesh(i,3))+x(mesh(i,5)))/2.d0 y(iset)=(y(mesh(i,3))+y(mesh(i,5)))/2.d0 z(iset)=(z(mesh(i,3))+z(mesh(i,5)))/2.d0 else mesh2(3+(i-1)*4,5)=narabi(mesh(i,3),mesh(i,5)) end if !3-6 if(narabi(mesh(i,3),mesh(i,6))==-1) then iset=iset+1 narabi(mesh(i,3),mesh(i,6))=iset !頂点1-4間と narabi(mesh(i,6),mesh(i,3))=iset !頂点4-1間とに節点番号isetをふり mesh2(3+(i-1)*4,6)=iset x(iset)=(x(mesh(i,3))+x(mesh(i,6)))/2.d0 y(iset)=(y(mesh(i,3))+y(mesh(i,6)))/2.d0 z(iset)=(z(mesh(i,3))+z(mesh(i,6)))/2.d0 else mesh2(3+(i-1)*4,6)=narabi(mesh(i,3),mesh(i,6)) end if end do !mesh(i,j)をmesh2(i,j)に入れる m1=m1*2 n1=n1*2 iset=m1*(n1+1) do i=1,iset do j=1,iset narabi(i,j)=-1 end do end do nyou=nyou*4 do i=1,nyou do j=1,6 mesh(i,j)=mesh2(i,j) end do end do if (nyou==2*n*m*4**ny) exit ! end do !節点番号と座標 print'("*NODE, NSET=Nall")' do i=1,iset print'(i6,a,f12.8,a,f12.8,a,f12.8)',& &i,',',x(i),',',y(i),',',z(i) end do !要素番号と節点番号 print'("*ELEMENT, TYPE=S6, ELSET=Eall")' do j=1,n-1,2 do i=1,m*2*4**(ny)-1,2 print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& &i+(j-1)*m*2*4**(ny),',',& &mesh2(i+(j-1)*m*2*4**(ny),1),',',& &mesh2(i+(j-1)*m*2*4**(ny),2),',',& &mesh2(i+(j-1)*m*2*4**(ny),3),',',& &mesh2(i+(j-1)*m*2*4**(ny),4),',',& &mesh2(i+(j-1)*m*2*4**(ny),5),',',& &mesh2(i+(j-1)*m*2*4**(ny),6) end do end do do j=1,n-1,2 do i=2,m*2*4**(ny),2 print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& &i+(j-1)*m*2*4**(ny),',',& &mesh2(i+(j-1)*m*2*4**(ny),1),',',& &mesh2(i+(j-1)*m*2*4**(ny),3),',',& &mesh2(i+(j-1)*m*2*4**(ny),2),',',& &mesh2(i+(j-1)*m*2*4**(ny),6),',',& &mesh2(i+(j-1)*m*2*4**(ny),5),',',& &mesh2(i+(j-1)*m*2*4**(ny),4) end do end do do j=2,n,2 do i=1,m*2*4**(ny)-1,2 print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& &i+(j-1)*m*2*4**(ny),',',& &mesh2(i+(j-1)*m*2*4**(ny),1),',',& &mesh2(i+(j-1)*m*2*4**(ny),3),',',& &mesh2(i+(j-1)*m*2*4**(ny),2),',',& &mesh2(i+(j-1)*m*2*4**(ny),6),',',& &mesh2(i+(j-1)*m*2*4**(ny),5),',',& &mesh2(i+(j-1)*m*2*4**(ny),4) end do end do do j=2,n,2 do i=2,m*2*4**(ny),2 print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',& &i+(j-1)*m*2*4**(ny),',',& &mesh2(i+(j-1)*m*2*4**(ny),1),',',& &mesh2(i+(j-1)*m*2*4**(ny),2),',',& &mesh2(i+(j-1)*m*2*4**(ny),3),',',& &mesh2(i+(j-1)*m*2*4**(ny),4),',',& &mesh2(i+(j-1)*m*2*4**(ny),5),',',& &mesh2(i+(j-1)*m*2*4**(ny),6) end do 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)h-0.00001) then print '(i6,a)',i,',' end if end do !境界条件の設定 print'("*BOUNDARY")' print'("kotei,1,3")' !材料の設定 print'("*MATERIAL,NAME=EL")' print'("*ELASTIC")' print'("69000.0,0.30000")' !要素の設定,どの材料を使うか等 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 ',',',3,',',-1.00E-003/((m*2)*(2**ny)) !出力のオプション 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 !