! EUC-jp ! 文字コードはEUC-jp ! 自由形式のfortranで書いています。 ! g77でコンパイルする場合は、 !g77 -ffree-form -o smv smv.f ! のようにコンパイルしてください。 ! gfortran でコンパイルする場合は文字コードをutf-8にしてください。 !ble6ccx.fのプログラムを折り畳みの拘束、節点番号を出力させる !ように書き換えたプログラム !折り畳み円筒の座屈荷重を求めるのに必要な !inpファイルを出力させるプログラム ! 09/3/22版 ! プログラムを主に作った人:後藤文彦 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! このプログラムは無保証です。 ! このプログラムは改造/再配布して構いません。 ! 但し、私的利用に留まらない場合は、もとのプログラムが ! http://www.str.ce.akita-u.ac.jp/~kudo/prog/ ! から取ってきたものと分かるように明示して下さい。 ! 改造/再配布後も上記の条件を保持して下さい。 ! バグの報告は、 ! http://www.str.ce.akita-u.ac.jp/~gotou/png/meeru.png ! にお願いします。 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc program ble6ccx implicit real*8 (a-h, o-z) dimension x(3000000),y(3000000),z(3000000),& & mesh(1000000,6),narabi(10000,10000) character mozi*20,mozi0*20 read*, mozi0 !1行目の文字を読み込む read*, nset !節点数の読み込み do i=0,nset-1 read*, x(i),y(i),z(i) !各節点の座標値の読み込み end do ! open(7,file='yousosuu.txt') !yousosuu.fで数えた要素数ファイル read(7,*) nyou !要素数の読み込み read(7,*) nsyuubun,ntakabun !周方向、高さ方向の分割数の読み込み do i=1,nyou !各要素の頂点数?、各頂点の節点番号、カラーコード?の読み込み read*,kazu,mesh(i,1),mesh(i,2),mesh(i,3),mozi end do ! do i=0,nset do j=0,nset !隣り合う節点i,jの中間の節点番号の初期値をすべて-1にしておく narabi(i,j)=-1 end do end do ! ! !中間点に新たに節点番号をふっていくのだが、 !頂点だけの節点数がnset個で、blenderの節点番号は0から始まるから、 !まだ使ってない一番若い節点番号はnset-1番 iset=nset-1 do i=1,nyou !1要素目からnyou要素まで順番に中間節点をふる !頂点1-2間に中間節点番号がまだふられていなければ、 if(narabi(mesh(i,1),mesh(i,2))==-1) then iset=iset+1 !節点番号を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 !その節点番号をi要素目の頂点5に覚えさせる !頂点5のxyz座標を頂点1-2間の中点として与える 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 !頂点1-2間にすでに節点番号がふられていれば、 mesh(i,4)=narabi(mesh(i,1),mesh(i,2)) !それを頂点5の節点番号にする end if ! 頂点2-3,3-1間も以下同様に if(narabi(mesh(i,2),mesh(i,3))==-1) then iset=iset+1 narabi(mesh(i,2),mesh(i,3))=iset narabi(mesh(i,3),mesh(i,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 ! if(narabi(mesh(i,3),mesh(i,1))==-1) then iset=iset+1 narabi(mesh(i,3),mesh(i,1))=iset narabi(mesh(i,1),mesh(i,3))=iset mesh(i,6)=iset x(iset)=(x(mesh(i,3))+x(mesh(i,1)))/2.d0 y(iset)=(y(mesh(i,3))+y(mesh(i,1)))/2.d0 z(iset)=(z(mesh(i,3))+z(mesh(i,1)))/2.d0 else mesh(i,6)=narabi(mesh(i,3),mesh(i,1)) end if nset=iset end do ! ! !!!!!!!!!!!!!!!!!!!!!! !open(9,file='ble6ccx.inp') print'("** blenderのobjファイルをccxのS6要素inp化")' print'("**")' print'("*NODE, NSET=Nall")' do i=0,nset print'(i6,a,f9.6,a,f9.6,a,f9.6)'& &,i+1,',',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)+1,',',mesh(i,2)+1,',',mesh(i,3)+1,',',mesh(i,4)+1& & ,',',mesh(i,5)+1,',',mesh(i,6)+1 end do !同じ内容をkireme.txtの中にコピーする。 open(7,file='kireme.txt') !節点数と各節点 write(7,*) nset+1 do i=0,nset write(7,*) i+1,x(i),y(i),z(i) end do !要素数と各要素 write(7,*) nyou do i=1,nyou write(7,*) i,mesh(i,1)+1,mesh(i,2)+1,mesh(i,3)+1,mesh(i,4)+1,mesh(i,5)+1,mesh(i,6)+1 end do print'("*NSET,NSET=FIX")' print'("** ここに拘束したい節点番号とコンマを1行ずつ並べる")' ! !ここで下端開口部の拘束する節点番号を入力 do i=0,nsyuubun-2,2 !周方向に print*, i*3+1,',' print*, i*3+2,',' print*, i*3+6,',' end do print*,nsyuubun*3/2*(ntakabun+1)+1,',' do i=1,nsyuubun*3/2-1 print*,nsyuubun*3/2*(ntakabun+1)+2+(4*i),',' end do ! !ここで上端開口部の拘束する節点番号を入力 print'("*NSET,NSET=FIX2")' print'("** ここに拘束したい節点番号とコンマを1行ずつ並べる")' do i=1,nsyuubun*3/2 print*, nsyuubun*3*ntakabun/2+i,',' end do do i=0,nsyuubun*3/2-2 print*, (nsyuubun*ntakabun*6)-nsyuubun/2-(nsyuubun-4)+(3*i),',' end do print*, nsyuubun*3*(2*ntakabun+1),',' print*, 30*nsyuubun/2*ntakabun/2-6*nsyuubun/2*(ntakabun/2-1)+ntakabun/2*2+1,',' ! !ここで境界条件を入力 print'("*BOUNDARY")' print'("** 以下のFIX,に続けて拘束変位1=x,2=y,3=zをコンマ区切りで")' print'("FIX,1,3")' print'("FIX2,1,2")' !if(nsyuubun % 4==0) then !print*, nsyuubun*2-(nsyuubun/2-1),',2,3'!yz拘束 !print'("63,2,3")'!yz拘束 !else if(nsyuubun==10) then !print'("105,2,3")'!yz拘束 !else if(nsyuubun==14) then !print'("147,2,3")'!yz拘束 !else !print*, nsyuubun*10+nsyuubun/2,',2,3'!yz拘束 !end if ! !ここで材料定数を定義 print'("*MATERIAL,NAME=EL")' print'("*ELASTIC")' e=210000. ! ヤング率(MPa) poi=0.300000012 ! ポアソン比 print*, e,',',poi t=10.d-4 !板の厚さ(m) print'("*SHELL SECTION,ELSET=Eall,MATERIAL=EL,OFFSET=0.")' print*,t !温度設定 !print'("*INITIAL CONDITIONS,TYPE=TEMPERATURE")' !print'("Nall,273.")' !節点ごとに厚さを変えたい時 !print'("*NSET,NSET=ATUSA")' !do i=1,nset !print'(i6,a)',i,',' !end do !print'("*NODAL THICKNESS")' !print*,'ATUSA,',t ! !ここで載荷する節点番号を入力 print'("*NSET,NSET=SAIKA")' print'("** ここに載荷したい節点番号とコンマを1行ずつ並べる")' do i=1,nsyuubun*3/2 print*, nsyuubun*3*ntakabun/2+i,',' end do do i=0,nsyuubun*3/2-2 print*, (nsyuubun*ntakabun*6)-nsyuubun/2-(nsyuubun-4)+(3*i),',' end do print*, nsyuubun*3*(2*ntakabun+1),',' print*, 30*nsyuubun/2*ntakabun/2-6*nsyuubun/2*(ntakabun/2-1)+ntakabun/2*2+1,',' ! print'("*STEP")' !print'("*STEP,NLGEON")' print'("*STATIC,SOLVER=SPOOLES")' !座屈荷重を求める場合は、以下のBUCKLEを入力 !print'("*BUCKLE")' !print'("5,0.01")' print'("*CLOAD")' p=1.d-4 !載荷荷重(MN) print'("** 以下のSAIKAに続けて載荷方向、載荷荷重を書く")' print'("** 例えば載荷方向が3=zで、載荷荷重がp/nsaikaなら")' !荷重は載荷節点数で割る print*,'SAIKA,3,',p/2/(nsyuubun*3/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")' ! ! !!!!!!!!!!!!!!!!!!!!!! ! open(8,file='ble6ccx.obj') write(8,'("3DG1")') write(8,*) nset+1 do i=0,nset write(8,*) x(i),y(i),z(i) end do ! do i=1,nyou write(8,*) 6,mesh(i,1),mesh(i,2),mesh(i,3),mesh(i,4)& & ,mesh(i,5),mesh(i,6),' ','0xcccccc' end do ! !------------------------------------------------------------------------------ ! open(8,file='suuti.d') write(8,*) ,nset+1 do i=0,nset write(8,*) ,i+1,',',x(i),',',y(i),',',z(i) end do write(8,*) ,nyou do i=1,nyou write(8,*),& i,',',mesh(i,1)+1,',',mesh(i,2)+1,',',mesh(i,3)+1,',',mesh(i,4)+1& & ,',',mesh(i,5)+1,',',mesh(i,6)+1 end do end