implicit real*8(a-h, o-z) dimension f(18),& !節点力ベクトル & d(18),& !節点変位ベクトル & x(18),& !境界条件(拘束節点に0, その他に1が入る) & ss(18,18),& !全体剛性行列 & xv(9),xw(9),& !節点ごとの変位境界条件 & fy(9),fz(9),& !節点ごとの荷重条件 & idari(9),migi(9),& !各要素の左節点番号、右節点番号 & s(9,4,4),& !要素剛性行列 & th(9),& !各要素の回転角 & t(4,4),& !座標変換行列(1要素ごとに計算) & s44(4,4),& !要素剛性行列(1要素ごとに移し替える入れ物) & ts(4,4),& !TK(1要素ごとに計算) & tst(4,4) !TKT(1要素ごとに計算) ! pi=asin(1.d0)*2.d0 !nyou=2 !要素数 例題1 !nset=3 !節点数 例題1 ! nyou=5 !要素数 例題2 nset=4 !節点数 例題2 ! ! 初期化 do n=1,nyou do i=1,4 do j=1,4 s(n,i,j)=0.d0 end do end do end do ! do i=1,18 d(i)=0.d0 !; f(i)=0.d0; x(i)=0.d0 do j=1,18 ss(i,j)=0.d0 end do end do ! do i=1,9 xv(i)=1.d0; xw(i)=1.d0 !境界条件は1で初期化 fy(i)=0.d0; fz(i)=0.d0 end do ! ! 例題1 !! 要素1の剛性マトリクス !s(1,2,2)= 1.d0; s(1,2,4)=-1.d0 !s(1,4,2)=-1.d0; s(1,4,4)= 1.d0 !! 要素2の剛性マトリクス !sk2=1.d0/sqrt(2.d0) !s(2,2,2)= sk2; s(2,2,4)=-sk2 !s(2,4,2)=-sk2; s(2,4,4)= sk2 ! ! ! 例題2 ! 要素1の剛性マトリクス s(1,2,2)= 1.d0; s(1,2,4)=-1.d0 s(1,4,2)=-1.d0; s(1,4,4)= 1.d0 ! 要素2の剛性マトリクス s(2,2,2)= 1.d0; s(2,2,4)=-1.d0 s(2,4,2)=-1.d0; s(2,4,4)= 1.d0 ! 要素3の剛性マトリクス s(3,2,2)= 1.d0; s(3,2,4)=-1.d0 s(3,4,2)=-1.d0; s(3,4,4)= 1.d0 ! 要素4の剛性マトリクス s(4,2,2)= 1.d0; s(4,2,4)=-1.d0 s(4,4,2)=-1.d0; s(4,4,4)= 1.d0 ! 要素5の剛性マトリクス sk5=1.d0/sqrt(2.d0) s(5,2,2)= sk5; s(5,2,4)=-sk5 s(5,4,2)=-sk5; s(5,4,4)= sk5 ! !各要素の左節点番号、右節点番号 ! 例題1 !idari(1)=1; migi(1)=2 !idari(2)=2; migi(2)=3 ! 例題2 idari(1)=1; migi(1)=2 idari(2)=2; migi(2)=3 idari(3)=3; migi(3)=4 idari(4)=4; migi(4)=1 idari(5)=1; migi(5)=3 !! !各要素の回転角(上の左右節点番号がz軸に横たわる状態からの) ! 例題1 !th(1)=0.d0 !th(2)=pi+pi/4.d0 ! ! 例題2 th(1)=pi/2.d0 th(2)=0.d0 th(3)=-pi/2.d0 th(4)=pi th(5)=pi/4.d0 ! !境界条件 ! 例題1 !xv(1)=0.d0; xw(1)=0.d0 !xv(3)=0.d0; xw(3)=0.d0 ! ! 例題2 xv(1)=0.d0; xw(1)=0.d0 xv(4)=0.d0 ! do i=1,9 x(2*i-1)=xv(i) x(2*i)=xw(i) end do ! !載荷条件 !fy(2)=1.d0 ! 例題1 ! fz(3)=1.d0 ! 例題2 ! do i=1,9 f(2*i-1)=fy(i) f(2*i)=fz(i) end do ! ! !要素ごとのTKTの計算 do n=1,nyou ! ! 座標変換マトリクス ! call zahyou(t,th(n)) ! do i=1,4 do j=1,4 s44(i,j)=s(n,i,j) end do end do ! call mxtmx(t,s44,ts) ! call mxmx(ts,t,tst) !! !!重ねあわせ i1=2*idari(n)-1 i2=2*idari(n) m1=2*migi(n)-1 m2=2*migi(n) ! ss(i1,i1)=ss(i1,i1)+tst(1,1) ss(i1,i2)=ss(i1,i2)+tst(1,2) ss(i1,m1)=ss(i1,m1)+tst(1,3) ss(i1,m2)=ss(i1,m2)+tst(1,4) ! ! ss(i2,i1)=ss(i2,i1)+tst(2,1) ss(i2,i2)=ss(i2,i2)+tst(2,2) ss(i2,m1)=ss(i2,m1)+tst(2,3) ss(i2,m2)=ss(i2,m2)+tst(2,4) ! ! ss(m1,i1)=ss(m1,i1)+tst(3,1) ss(m1,i2)=ss(m1,i2)+tst(3,2) ss(m1,m1)=ss(m1,m1)+tst(3,3) ss(m1,m2)=ss(m1,m2)+tst(3,4) ! ! ss(m2,i1)=ss(m2,i1)+tst(4,1) ss(m2,i2)=ss(m2,i2)+tst(4,2) ss(m2,m1)=ss(m2,m1)+tst(4,3) ss(m2,m2)=ss(m2,m2)+tst(4,4) ! ! end do !n要素について !! !! !do i=1,8 !print'(8f10.3)', (ss(i,j),j=1,8) !end do !print* !! ! 境界条件を入れる do i=1,nset*2 do j=1,nset*2 ss(i,j)=x(i)*ss(i,j) ss(j,i)=x(i)*ss(j,i) end do end do do i=1,nset*2 if(x(i)<1.d-3) then ss(i,i)=1.d0 end if end do !! !! ! !print*,'f=',(f(j),j=1,18) ! ! call gausu(nset*2,ss,d,f) ! !print* !do i=1,8 !print'(8f10.3)', (ss(i,j),j=1,8) !end do ! do n=1,nset print*,'節点番号:',n,'v=',d(n*2-1),'w=',d(n*2) end do !print*,'d=',(d(j),j=1,18) end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1! ! 以上がメインプログラム !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1! ! subroutine mxmx(a,b,c) implicit real*8(a-h, o-z) dimension a(4,4), b(4,4), c(4,4) ! ! do i=1,4 do j=1,4 c(i,j)=0. do k=1,4 c(i,j)=c(i,j)+a(i,k)*b(k,j) end do end do end do ! return end ! subroutine zahyou(t,theta) implicit real*8(a-h, o-z) dimension t(4,4) do i=1,2 do j=1,2 t(i+2,j)=0.d0 t(i,j+2)=0.d0 end do end do t(1,1)= cos(theta) t(1,2)= sin(theta) t(2,1)=-sin(theta) t(2,2)= cos(theta) t(3,3)=t(1,1) t(3,4)=t(1,2) t(4,3)=t(2,1) t(4,4)=t(2,2) return end ! subroutine mxtmx(a,b,c) implicit real*8(a-h, o-z) dimension a(4,4), b(4,4), c(4,4) ! do i=1,4 do j=1,4 c(i,j)=0. do k=1,4 c(i,j)=c(i,j)+a(k,i)*b(k,j) end do end do end do ! ! return end ! ! subroutine gausu(n,a,x,b) implicit real*8(a-h, o-z) !dimension a(n,n),x(n),b(n) dimension a(18,18),x(18),b(18) ! !ガウスの消去法の参考としたのは、 !名取亮「すうがくぶっくす12 線形計算」(朝倉書店)p.10-15 !! !print* !do i=1,8 !print'(8f10.3)', (a(i,j),j=1,8) !end do ! ! do k=1,n-1 !a(k,k)を消去 do i=k+1,n !k+1行からn行まで do j=k+1,n !k+1列からn列まで a(i,j)=a(i,j)-a(k,j)*a(i,k)/a(k,k) end do b(i)=b(i)-b(k)*a(i,k)/a(k,k) end do end do ! ! 後退代入 x(n)=b(n)/a(n,n) do k=n,1,-1 akjxj=0.d0 do j=k+1,n akjxj=akjxj+a(k,j)*x(j) end do x(k)=(b(k)-akjxj)/a(k,k) end do ! return end ! !