Deprecated: The each() function is deprecated. This message will be suppressed on further calls in /home/zhenxiangba/zhenxiangba.com/public_html/phproxy-improved-master/index.php on line 456
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
!
!