📄 3dframe.f90
字号:
!===================================================================================
! 空间框架计算程序
! 程序说明：
! 1：编程语言为Fortran 90
! 2：本程序可以求解任意空间杆系结构
! 3：节点位移编码为（X, Y, Z, θx, θy, θz），荷载编码为（Fx, Fy, Fz, Mx, My, Mz）
! 4：单元输入信息为：起点号，终点号，EA, EIy, EIz, GIx, α(α为参考坐标系与整体坐标系
! 夹角
! 程序清单：
!=====================================================================================
!1：Lxz_Tools
module Lxz_Tools
implicit none
integer (kind(1)),parameter ::ikind=(kind(1))
integer (kind(1)),parameter ::rkind=(kind(0.D0))
real (rkind), parameter :: Zero=0.D0,One=1.D0,Two=2.D0,Three=3.D0, &
& Four=4.D0,Five=5.D0,Six=6.D0,Seven=7.D0,Eight=8.D0,Nine=9.D0, &
& Ten=10.D0
end module Lxz_Tools
!2：TypeDef
module TypeDef
use Lxz_Tools
implicit none
integer(ikind),parameter :: NDOF=6,NNode=2 ! 每个节点6个自由度，每个单元2个节点
type::typ_Joint ! 节点类型
real(rkind) :: X,Y,Z ! 节点坐标
integer(ikind):: GDOF(NDOF) ! 节点自由度整体编码
end type typ_Joint
type::typ_Element ! 单元类型
integer(ikind):: JointNo(NNode) ! 节点编号
real(rkind):: EIy,EIz,EA,GIp,Length,CosA,CosB,CosC,A ! 单元属性
real(rkind)::EK(NDOF*NNode,NDOF*NNode),ET(NDOF*NNode,NDOF*NNode) ! 刚度矩阵，坐标转换矩阵
integer(ikind)::GlbDOF(NDOF*NNode) ! 整体自由度编码
real(rkind)::EForce(NDOF*NNode),ELForce(NDOF*NNode) ! 单元内力
end type typ_Element
type::typ_JointLoad ! 节点荷载类型
integer(ikind)::JointNo,LodDOF ! 节点编号，荷载自由度
real(rkind)::LodVal ! 荷载大小
end type typ_JointLoad
type::typ_ElemLoad ! 单元荷载类型
integer (ikind):: ElemNo,Indx ! 单元编号，单元荷载形式编号
real(rkind)::Pos,LodVal ! 荷载位置，荷载大小
end type typ_ElemLoad
contains
subroutine SetElemProp(Elem,Joint) ! 计算单元参数，得到单元刚度矩阵和坐标转换矩阵
type(typ_Element),intent(inout)::Elem(:)
type(typ_Joint),intent(in)::Joint(:)
integer(ikind):: i,EJ1,EJ2
real(rkind)::T(3,3),x,y,z
real(rkind)::cx,cy,cz,cs,ca,sa
do i=1,size(Elem)
Elem(i)%EK=zero;Elem(i)%ET=zero;T=zero
EJ1=Elem(i)%JointNo(1);EJ2=Elem(i)%JointNo(2)
Elem(i)%GlbDOF(1:6)=Joint(EJ1)%GDOF(:)
Elem(i)%GlbDOF(7:12)=Joint(EJ2)%GDOF(:)
x=Joint(EJ2)%X-Joint(EJ1)%X
y=Joint(EJ2)%Y-Joint(EJ1)%Y
z=Joint(EJ2)%Z-Joint(EJ1)%Z
Elem(i)%Length=sqrt(x**2+y**2+z**2)
Elem(i)%CosA=x/Elem(i)%Length
Elem(i)%CosB=y/Elem(i)%Length
Elem(i)%CosC=z/Elem(i)%Length
cx=Elem(i)%CosA;cy=Elem(i)%CosB;cz=Elem(i)%CosC
ca=cos(Elem(i)%A);sa=sin(Elem(i)%A)
cs=sqrt(cx**2+cy**2)
if(cs.NE.zero) then
T(1,1)=cx;T(1,2)=cy;T(1,3)=cz;
T(2,1)=-(ca*cy+sa*cx*cz)/cs;
T(2,2)=(ca*cx-sa*cy*cz)/cs;
T(2,3)=cs*sa;
T(3,1)=(sa*cy-ca*cx*cz)/cs;
T(3,2)=-(sa*cx+ca*cy*cz)/cs;
T(3,3)=cs*ca;
else
T(1,3)=one;
T(2,1)=-sa;T(2,2)=ca;
T(3,1)=-ca;T(3,2)=-sa;
end if
Elem(i)%ET(1:3,1:3)=T;Elem(i)%ET(4:6,4:6)=T;
Elem(i)%ET(7:9,7:9)=T;Elem(i)%ET(10:12,10:12)=T;
T=zero
T(1,1)=Elem(i)%EA/Elem(i)%Length
T(2,2)=12D0*Elem(i)%EIz/(Elem(i)%Length**3)
T(3,3)=12D0*Elem(i)%EIy/(Elem(i)%Length**3)
Elem(i)%EK(1:3,1:3)=T
Elem(i)%EK(7:9,7:9)=T
Elem(i)%EK(1:3,7:9)=-T
Elem(i)%EK(7:9,1:3)=transpose(-T)
T=zero
T(1,1)=Elem(i)%GIp/Elem(i)%Length
T(2,2)=4D0*Elem(i)%EIy/Elem(i)%Length
T(3,3)=4D0*Elem(i)%EIz/Elem(i)%Length
Elem(i)%EK(4:6,4:6)=T;Elem(i)%EK(10:12,10:12)=T
T=zero
T(2,3)=6D0*Elem(i)%EIz/(Elem(i)%Length**2)
T(3,2)=-6D0*Elem(i)%EIy/(Elem(i)%Length**2)
Elem(i)%EK(1:3,4:6)=T;Elem(i)%EK(1:3,10:12)=T
Elem(i)%EK(4:6,1:3)=transpose(T)
Elem(i)%EK(10:12,1:3)=transpose(T)
Elem(i)%EK(7:9,10:12)=-T
Elem(i)%EK(10:12,7:9)=-transpose(T)
T=zero
T(2,3)=6D0*Elem(i)%EIy/(Elem(i)%Length**2)
T(3,2)=-6D0*Elem(i)%EIz/(Elem(i)%Length**2)
Elem(i)%EK(4:6,7:9)=T
Elem(i)%EK(7:9,4:6)=transpose(T)
T=zero
T(1,1)=-Elem(i)%GIp/Elem(i)%Length
T(2,2)=2D0*Elem(i)%EIy/Elem(i)%Length
T(3,3)=2D0*Elem(i)%EIz/Elem(i)%Length
Elem(i)%EK(4:6,10:12)=T
Elem(i)%EK(10:12,4:6)=transpose(T)
end do
end subroutine SetElemProp
end module TypeDef
!3：SolveDisp
module Solve ! 矩阵求解模块
use TypeDef
type :: typ_Kcol
real(rkind),pointer :: row(:)
end type typ_Kcol
contains
!==================================================
subroutine SolveDisp (GDisp, Elem,Joint,GLoad)
!==================================================
real(rkind),intent(out) :: GDisp(:)
type (typ_Element),intent(in) :: Elem(:)
type (typ_Joint),intent(in) :: Joint(:)
real(rkind) :: GLoad(:)
integer(ikind) NElem,NGlbDOF
type (typ_Kcol),allocatable ::Kcol(:)
NElem = size(Elem,dim=1) ! 单元数量
NGlbDOF = size(GDisp,dim=1) ! 总自由度数目
allocate (Kcol(NGlbDOF))
call SetMatBand() ! 设定矩阵带宽
call GStifMat() ! 刚度矩阵求解
call BandSolv() ! 刚度矩阵求解
contains
!--------------------------
subroutine SetMatBand()
!--------------------------
integer (ikind) :: minDOF
integer (ikind),allocatable :: Row1(:)
integer (ikind) :: ie,j
integer (ikind),allocatable::ELocVec(:)
allocate (Row1(NGlbDOF))
allocate(ELocVec(size(Elem(1)%GlbDOF)))
Row1=NGlbDOF
do ie=1,NElem
ELocVec(:)=Elem(ie)%GlbDOF(:)
minDOF=minval(ELocVec,mask=ELocVec>0)
where(ELocVec>0)
Row1(ELocVec)=min(Row1(ELocVec),minDOF)
end where
end do
do j=1,NGlbDOF
allocate (Kcol(j)%row(Row1(j):j))
Kcol(j)%row=Zero
end do
return
end subroutine SetMatBand
!---------------------------------------
subroutine BandSolv()
!---------------------------------------
integer(ikind)::row1,ncol,row,j,ie
real(rkind)::diag(1:NGlbDOF),s
ncol=NGlbDOF
diag(1:ncol)=(/(Kcol(j)%row(j),j=1,ncol)/)
do j=2,ncol
row1=lbound(Kcol(j)%row,1)
do ie=row1,j-1
row=max(row1,lbound(Kcol(ie)%row,1))
s=sum(diag(row:ie-1)*Kcol(ie)%row(row:ie-1)*Kcol(j)%row(row:ie-1))
Kcol(j)%row(ie)=(Kcol(j)%row(ie)-s)/diag(ie)
end do
s=sum(diag(row1:j-1)*Kcol(j)%row(row1:j-1)**2)
diag(j)=diag(j)-s
end do
do ie=2,ncol
row1=lbound(Kcol(ie)%row,dim=1)
GLoad(ie)=GLoad(ie)-sum(Kcol(ie)%row(row1:ie-1)*GLoad(row1:ie-1))
end do
GLoad(:)=GLoad(:)/diag(:)
do j=ncol,2,-1
row1=lbound(Kcol(j)%row,dim=1)
GLoad(row1:j-1)=GLoad(row1:j-1)-GLoad(j)*Kcol(j)%row(row1:j-1)
end do
GDisp(:)=GLoad(:)
return
end subroutine BandSolv
!---------------------
subroutine GStifMat()
!---------------------
integer(ikind)::ie,j,JGDOF
integer (ikind),allocatable::ELocVec(:)
real(rkind),allocatable::EK(:,:),ET(:,:)
integer(ikind)::n
n=size(Elem(1)%GlbDOF)
allocate(ELocVec(n))
allocate(EK(n,n))
allocate(ET(n,n))
do IE=1,NElem
EK=Elem(IE)%EK
ET=Elem(IE)%ET
EK = matmul(transpose(ET),matmul(EK,ET))
ELocVec(:)=Elem(IE)%GlbDOF(:)
do j=1,12
JGDOF=ELocVec(j)
if (JGDOF==0) cycle
where (ELocVec>0.and.ELocVec<=JGDOF)
Kcol(JGDOF)%row(ELocVec)=Kcol(JGDOF)%row(ELocVec)+EK(:,j)
end where
end do
end do
return
end subroutine GStifMat
end subroutine SolveDisp
end module Solve
!4：DispMethod
module DispMethod ! 荷载集成
use Solve
implicit none
contains
subroutine GLoadVec(Elem,Joint,JLoad,ELoad,GLoad)
type (typ_Element),intent(inout)::Elem(:)
type (typ_Joint),intent(in)::Joint(:)
type (typ_JointLoad),intent(in)::JLoad(:)
type (typ_ElemLoad),intent(in)::ELoad(:)
real(rkind)::GLoad(:),EF(NNode*NDOF)
integer (ikind)::i,m,n
real(rkind)::a,b,l,q
do i=1,size(JLoad)
n=JLoad(i)%JointNo
m=JLoad(i)%LodDOF
m=Joint(n)%GDOF(m)
GLoad(m)=GLoad(m)+JLoad(i)%LodVal
end do
do i=1,size(ELoad)
n=ELoad(i)%ElemNo
l=Elem(n)%Length
a=ELoad(i)%Pos
q=ELoad(i)%LodVal
b=l-a
EF=zero
if(ELoad(i)%Indx.eq.1) then ! 均布荷载，x,z平面
EF(5)=q*a*a*(6D0-8D0*a/l+3D0*a*a/(l*l))/12D0
EF(11)=-q*(a**3D0)*(4D0-3D0*a/l)/(12D0*l)
EF(3)=-0.5D0*q*a*(2D0-2D0*a*a/(l*l)+a**3D0/l**3D0)
EF(9)=-0.5D0*q*a**3D0*(2D0-a/l)/l**2D0
end if
if(ELoad(i)%Indx.eq.2) then ! 集中荷载，x,z平面
EF(5)=q*a*b**2D0/l**2D0
EF(11)=-q*a*a*b/l**2D0
EF(3)=-q*b**2D0*(1D0+2D0*a/l)/l**2D0
EF(9)=-q*a**2D0*(1D0+2D0*b/l)/l**2D0
end if
if(ELoad(i)%Indx.eq.-1) then ! 均布荷载，x,y平面
EF(6)=-q*a*a*(6D0-8D0*a/l+3D0*a*a/(l*l))/12D0
EF(12)=q*(a**3D0)*(4D0-3D0*a/l)/(12D0*l)
EF(2)=-0.5D0*q*a*(2D0-2D0*a*a/(l*l)+a**3D0/l**3D0)
EF(8)=-0.5D0*q*a**3D0*(2D0-a/l)/l**2D0
end if
if(ELoad(i)%Indx.eq.-2) then ! 集中荷载，x,y平面
EF(6)=-q*a*b**2D0/l**2D0
EF(12)=q*a*a*b/l**2D0
EF(2)=-q*b**2D0*(1D0+2D0*a/l)/l**2D0
EF(8)=-q*a**2D0*(1D0+2D0*b/l)/l**2D0
end if
if(ELoad(i)%Indx.eq.3) then ! 集中弯矩，x，z平面
EF(5)=q*b*(2D0-3D0*b/l)/l
EF(11)=q*a*(2D0-3D0*a/l)/l
EF(3)=-6D0*q*a*b/l**3D0
EF(9)=6D0*q*a*b/l**3D0
end if
if(ELoad(i)%Indx.eq.-3) then ! 集中弯矩，x，y平面
EF(6)=q*b*(2D0-3D0*b/l)/l
EF(12)=q*a*(2D0-3D0*a/l)/l
EF(2)=6D0*q*a*b/l**3D0
EF(8)=-6D0*q*a*b/l**3D0
end if
Elem(n)%ELForce=Elem(n)%ELForce+EF
EF=matmul(transpose(Elem(n)%ET),EF)
where(Elem(n)%GlbDOF>0)
GLoad(Elem(n)%GlbDOF)=GLoad(Elem(n)%GlbDOF)-EF(:)
end where
end do
i=1
end subroutine GLoadVec
end module DispMethod
program DFrame
use DispMethod
implicit none
integer(ikind):: NElem,NJoint,NGlbDOF,NJLoad,NELoad
type(typ_Element),allocatable::Elem(:)
type(typ_Joint),allocatable::Joint(:)
type(typ_JointLoad),allocatable::JLoad(:)
type(typ_ElemLoad),allocatable::ELoad(:)
real(rkind),allocatable::Disp(:),GLoad(:)
call Input_Data()
call SetElemProp(Elem,Joint)
call GLoadVec(Elem,Joint,JLoad,ELoad,GLoad)
write(16,*) "整体荷载向量"
write(16,*) GLoad
write(16,*) "整体节点位移"
call SolveDisp(Disp,Elem,Joint,GLoad)
write(16,*) Disp
write(16,*)
call Output_Results()
stop
contains
subroutine Input_Data()
integer(ikind)::i
open(5,file='data.ipt',status='old',position='rewind')
open(16,file='data.opt',position='rewind')
read(5,*) NElem,NJoint,NGlbDOF,NJLoad,NELoad
allocate (Elem(NElem))
allocate (Joint(NJoint))
allocate (JLoad(NJLoad))
allocate (ELoad(NELoad))
allocate (Disp(NGlbDOF))
allocate (GLoad(NGlbDOF))
read(5,*) (Joint(i),i=1,NJoint)
read(5,*) (Elem(i)%JointNo,Elem(i)%EA,Elem(i)%EIy,&
Elem(i)%EIz,Elem(i)%GIp,Elem(i)%A,i=1,NElem)
if(NJLoad>0) read(5,*) (JLoad(i),i=1,NJLoad)
if(NELoad>0) read(5,*) (ELoad(i),i=1,NELoad)
end subroutine Input_Data
subroutine Output_Results()
integer(ikind)::i
real(rkind) :: EDisp(12)
do i=1,size(Elem)
EDisp=zero
where(Elem(i)%GlbDOF>0)
EDisp(:)=Disp(Elem(i)%GlbDOF)
end where
write(16,*) "单元 ",i,"位移"
write(16,*) EDisp
EDisp=matmul(Elem(i)%ET,EDisp)
Elem(i)%EForce=matmul(Elem(i)%EK,EDisp)+Elem(i)%ELForce
write(16,*) "单元 ",i,"内力"
write(16,*) Elem(i)%EForce
end do
end subroutine Output_Results
end program DFrame
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -