implicit real*8(a-h, o-z) dimension f(27),& !節点力ベクトル & d(27),& !節点変位ベクトル & x(27),& !境界条件(拘束節点に0, その他に1が入る) & ss(27,27),& !全体剛性行列 & xv(9),xw(9),xt(9),& !節点ごとの変位境界条件:v,w,θ & fy(9),fz(9),cx(9),& !節点ごとの荷重条件:S,N,M & idari(9),migi(9),& !各要素の左節点番号、右節点番号 & s(9,6,6),& !要素剛性行列 & th(9),& !各要素の回転角 & t(6,6),& !座標変換行列(1要素ごとに計算) & s66(6,6),& !要素剛性行列(1要素ごとに移し替える入れ物) & ts(6,6),& !TK(1要素ごとに計算) & tst(6,6) !TKT(1要素ごとに計算) ! pi=asin(1.d0)*2.d0 ! nyou=2 !要素数 例題1 nset=3 !節点数 例題1 ! ! 初期化 do n=1,nyou do i=1,6 do j=1,6 s(n,i,j)=0.d0 end do end do end do ! do i=1,27 d(i)=0.d0 do j=1,27 ss(i,j)=0.d0 end do end do ! do i=1,9 xv(i)=1.d0; xw(i)=1.d0; xt(i)=1.d0 !境界条件は1で初期化 fy(i)=0.d0; fz(i)=0.d0; cx(i)=0.d0 end do ! ! 例題1 ! 要素1の剛性マトリクス el=1.d0 !1本のL ea=1.d0 !EA ei=1.d0 !EI s(1,2,2)= ea/el ; s(1,2,5)= -ea/el s(1,5,2)= -ea/el ; s(1,5,5)= ea/el s(1,1,1)= 12.d0*ei/el**3; s(1,1,3)= -6.d0*ei/el**2 s(1,1,4)=-12.d0*ei/el**3; s(1,1,6)= -6.d0*ei/el**2 s(1,3,1)= -6.d0*ei/el**2; s(1,3,3)= 4.d0*ei/el s(1,3,4)= 6.d0*ei/el**2; s(1,3,6)= 2.d0*ei/el s(1,4,1)=-12.d0*ei/el**3; s(1,4,3)= 6.d0*ei/el**2 s(1,4,4)= 12.d0*ei/el**3; s(1,4,6)= 6.d0*ei/el**2 s(1,6,1)= -6.d0*ei/el**2; s(1,6,3)= 2.d0*ei/el s(1,6,4)= 6.d0*ei/el**2; s(1,6,6)= 4.d0*ei/el ! 要素2の剛性マトリクス do i=1,6 do j=1,6 s(2,i,j)=s(1,i,j) end do end do ! !do i=1,6 !print'(9f10.2)',(s(1,i,j),j=1,9) !end do !print* !! ! ! !各要素の左節点番号、右節点番号 ! 例題1 idari(1)=1; migi(1)=2 idari(2)=2; migi(2)=3 ! !各要素の回転角(上の左右節点番号がz軸に横たわる状態からの) ! 例題1 th(1)=0.d0 th(2)=pi+pi/2.d0 ! ! !境界条件 ! 例題1 xv(1)=0.d0; xw(1)=0.d0; xt(1)=0.d0 xv(3)=0.d0; xw(3)=0.d0; xt(3)=0.d0 ! ! do i=1,9 x(3*i-2)=xv(i) x(3*i-1)=xw(i) x(3*i)=xt(i) end do ! !載荷条件 cx(2)=1.d0 ! 例題1 ! do i=1,9 f(3*i-2)=fy(i) f(3*i-1)=fz(i) f(3*i)=cx(i) end do ! ! !要素ごとのTKTの計算 do n=1,nyou ! ! 座標変換マトリクス ! call zahyou(t,th(n)) ! do i=1,6 do j=1,6 s66(i,j)=s(n,i,j) end do end do ! call mxtmx(t,s66,ts) ! call mxmx(ts,t,tst) ! !do i=1,6 !print'(9f10.2)',(tst(i,j),j=1,9) !end do !print* !!! !!重ねあわせ i1=3*idari(n)-2 i2=3*idari(n)-1 i3=3*idari(n) m1=3*migi(n)-2 m2=3*migi(n)-1 m3=3*migi(n) ! ss(i1,i1)=ss(i1,i1)+tst(1,1) ss(i1,i2)=ss(i1,i2)+tst(1,2) ss(i1,i3)=ss(i1,i3)+tst(1,3) ss(i1,m1)=ss(i1,m1)+tst(1,4) ss(i1,m2)=ss(i1,m2)+tst(1,5) ss(i1,m3)=ss(i1,m3)+tst(1,6) ! ss(i2,i1)=ss(i2,i1)+tst(2,1) ss(i2,i2)=ss(i2,i2)+tst(2,2) ss(i2,i3)=ss(i2,i3)+tst(2,3) ss(i2,m1)=ss(i2,m1)+tst(2,4) ss(i2,m2)=ss(i2,m2)+tst(2,5) ss(i2,m3)=ss(i2,m3)+tst(2,6) ! ss(i3,i1)=ss(i3,i1)+tst(3,1) ss(i3,i2)=ss(i3,i2)+tst(3,2) ss(i3,i3)=ss(i3,i3)+tst(3,3) ss(i3,m1)=ss(i3,m1)+tst(3,4) ss(i3,m2)=ss(i3,m2)+tst(3,5) ss(i3,m3)=ss(i3,m3)+tst(3,6) ! ! ! ss(m1,i1)=ss(m1,i1)+tst(4,1) ss(m1,i2)=ss(m1,i2)+tst(4,2) ss(m1,i3)=ss(m1,i3)+tst(4,3) ss(m1,m1)=ss(m1,m1)+tst(4,4) ss(m1,m2)=ss(m1,m2)+tst(4,5) ss(m1,m3)=ss(m1,m3)+tst(4,6) ! ss(m2,i1)=ss(m2,i1)+tst(5,1) ss(m2,i2)=ss(m2,i2)+tst(5,2) ss(m2,i3)=ss(m2,i3)+tst(5,3) ss(m2,m1)=ss(m2,m1)+tst(5,4) ss(m2,m2)=ss(m2,m2)+tst(5,5) ss(m2,m3)=ss(m2,m3)+tst(5,6) ! ss(m3,i1)=ss(m3,i1)+tst(6,1) ss(m3,i2)=ss(m3,i2)+tst(6,2) ss(m3,i3)=ss(m3,i3)+tst(6,3) ss(m3,m1)=ss(m3,m1)+tst(6,4) ss(m3,m2)=ss(m3,m2)+tst(6,5) ss(m3,m3)=ss(m3,m3)+tst(6,6) ! ! end do !n要素について !! !! !do i=1,8 !print'(8f10.3)', (ss(i,j),j=1,8) !end do !print* !! ! 境界条件を入れる do i=1,nset*3 do j=1,nset*3 ss(i,j)=x(i)*ss(i,j) ss(j,i)=x(i)*ss(j,i) end do end do do i=1,nset*3 if(x(i)<1.d-3) then ss(i,i)=1.d0 end if end do !! !! !do i=1,9 !print'(9f10.2)',(ss(i,j),j=1,9) !end do ! !print*,'f=',(f(j),j=1,18) ! ! call gausu(nset*3,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*3-2),'w=',d(n*3-1), 'th=',d(n*3) 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(6,6), b(6,6), c(6,6) ! ! do i=1,6 do j=1,6 c(i,j)=0. do k=1,6 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(6,6) do i=1,3 do j=1,3 t(i+3,j)=0.d0 t(i,j+3)=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)=1.d0 t(4,4)=t(1,1) t(4,5)=t(1,2) t(5,4)=t(2,1) t(5,5)=t(2,2) t(6,6)=1.d0 return end ! subroutine mxtmx(a,b,c) implicit real*8(a-h, o-z) dimension a(6,6), b(6,6), c(6,6) ! do i=1,6 do j=1,6 c(i,j)=0. do k=1,6 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(27,27),x(27),b(27) ! !ガウスの消去法の参考としたのは、 !名取亮「すうがくぶっくす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 ! !