truss.f90

完整有限元程序范例,内含C++核心算法的WORD文档及实例应用程序

!===========================================================
! 非线性分析程序设计范例
! 利用Newton-Raphson法求解平面弹塑性桁架
! 程序作者：陆新征
! 北京清华大学土木工程系
! 2006.5.
!===========================================================

module TypDef
	implicit none
	integer (kind(1)),parameter ::ikind=(kind(1))
	integer (kind(1)),parameter ::rkind=(kind(0.D0))
	
	type :: Typ_GValue  ! 整体控制参数
		integer (ikind) :: NNode, NElem, NLoad, NSupport ! 节点，单元，荷载，支座数目
		integer (ikind) :: NStep, CurrentStep ! 总加载步数, 当前加载步数
		real (rkind) :: Tol,Error ! 误差容限, 当前误差
	end type Typ_GValue

	type :: Typ_Node ! 定义节点
		real(rkind) :: Coord(2) ! 节点坐标
		real(rkind) :: Disp(2),dDisp(2)  ! 节点总位移，节点位移增量
	end type Typ_Node

	type :: Typ_Truss ! 定义桁架截面
		integer(ikind) :: NodeNo(2) ! 节点编号
		real(rkind) :: E, A, L ! 弹性模量，截面积，长度
		real(rkind) :: fy, fy_temp, Eh ! 屈服强度，，屈服强度临时变量，强化模量
		real(rkind) :: Et     ! 迭代刚度
		real(rkind) :: T(4,4) ! 坐标转化矩阵
		real(rkind) :: B(1,2) ! 几何矩阵
		real(rkind) :: K(4,4) ! 单元刚度矩阵
		real(rkind) :: EPS, dEPS, SIG ! 迭代开始时应变，应变增量，应力
		real(rkind) :: STRESS ! 迭代结束时应力
		real(rkind) :: R(4) ! 单元节点反力
	end type Typ_Truss

	type :: Typ_Load ! 定义荷载
		integer (ikind) :: NodeNo ! 作用节点
		real(rkind) :: Value(2)   ! 荷载大小
	end type Typ_Load

	type :: Typ_Support ! 定义支座
		integer (ikind) :: NodeNo ! 作用节点
		integer (ikind) :: DOF    ! 所约束的自由度
	end type Typ_Support

	contains

	subroutine Get_dDisp(GValue, Node, GDisp) ! 计算节点位移增量
		type(typ_GValue) :: GValue
		type(Typ_Node) :: Node(:)
		real(rkind) :: GDisp(:)
		integer (ikind) :: i
		do i=1, size(Node)
			Node(i)%dDisp(1:2)=Node(i)%dDisp(1:2)+GDisp(i*2-1:i*2)
		end do
		return
	end subroutine Get_dDisp

	subroutine Get_Stress(Elem) ! 计算应力增量和迭代弹性模量
		type(Typ_Truss) :: Elem
		real(rkind) :: S_trial, dE_pla, x
		Elem%Et=Elem%E
		S_trial=Elem%SIG+Elem%dEPS*Elem%E ! 得到弹性应力增量
		if(abs(S_trial)>Elem%fy) then ! 如果超过屈服应力
			dE_pla=(abs(S_trial)-Elem%fy)*sign(1.d0,S_trial)/Elem%E ! 得到塑性应变增量
			Elem%Stress=Elem%fy*sign(1.d0,S_trial)+dE_pla*Elem%Eh
			Elem%fy_temp=abs(Elem%Stress) ! 更新屈服应力
			if(Elem%dEPS.ne.0.d0) then
				Elem%Et=Elem%Eh
			end if
		else
			Elem%Stress=S_trial
		end if		
		return
	end subroutine Get_Stress

	subroutine Get_EK(Elem) ! 组装单元刚度矩阵
		type(Typ_Truss) :: Elem
		real(rkind) :: K(2,2)
		K=matmul(transpose(Elem%B),Elem%B)*Elem%Et*Elem%L*Elem%A
		Elem%K=0.d0
		Elem%K(1,1)=K(1,1); 
		Elem%K(3,3)=K(2,2); 
		Elem%K(1,3)=K(1,2); 
		Elem%K(3,1)=K(2,1)
		Elem%K=matmul(transpose(Elem%T),matmul(Elem%K,Elem%T)) ! K=T'KT
		return
	end subroutine Get_EK

	subroutine Get_R(Elem) ! 得到单元反力
		type(Typ_Truss) :: Elem
		real(rkind) :: R0(2,1), R1(4) ! 单元坐标系下的反力
		R0=transpose(Elem%B)*Elem%Stress*Elem%L*Elem%A ! R=B*SIG*L*A
		R1=0.d0;
		R1(1)=R0(1,1); R1(3)=R0(2,1);
		Elem%R=matmul(transpose(Elem%T),R1) ! R=T'R'
		return
	end subroutine Get_R
 
	subroutine Get_EPS(Node,Elem) ! 得到应变增量
		type(typ_Node) :: Node(:)
		type(typ_Truss) :: Elem
		real(rkind) :: D(4)
		D(1:2)=Node(Elem%NodeNo(1))%dDisp(1:2)
		D(3:4)=Node(Elem%NodeNo(2))%dDisp(1:2)
		D=matmul(Elem%T,D) ! D'=TD
		Elem%dEPS=D(1)*Elem%B(1,1)+D(3)*Elem%B(1,2)
		return
	end subroutine Get_EPS
	
	subroutine Intial_Elem(Node,Elem) ! 初始化单元，得到单元坐标变化矩阵和几何矩阵
		type(typ_Node) :: Node(:)
		type(typ_Truss) :: Elem(:)
		integer(ikind) :: ElemNo
		real(rkind) :: x1,x2,y1,y2,cosA,sinA
		do ElemNo=1,size(Elem)
			x1=Node(Elem(ElemNo)%NodeNo(1))%Coord(1)
			x2=Node(Elem(ElemNo)%NodeNo(2))%Coord(1)
			y1=Node(Elem(ElemNo)%NodeNo(1))%Coord(2)
			y2=Node(Elem(ElemNo)%NodeNo(2))%Coord(2)
			Elem(ElemNo)%L=sqrt((x1-x2)**2+(y1-y2)**2) ! 得到单元长度
			cosA=(x2-x1)/Elem(ElemNo)%L
			sinA=(y2-y1)/Elem(ElemNo)%L
			Elem(ElemNo)%T=0.d0
			Elem(ElemNo)%T(1,1)= cosA;	Elem(ElemNo)%T(1,2)=sinA; ! 得到单元坐标转换矩阵
			Elem(ElemNo)%T(2,1)=-sinA;	Elem(ElemNo)%T(2,2)=cosA;
			Elem(ElemNo)%T(3:4,3:4)=Elem(ElemNo)%T(1:2,1:2)
			Elem(ElemNo)%B(1,1)=-1.d0/Elem(ElemNo)%L ! 得到单元几何矩阵
			Elem(ElemNo)%B(1,2)=+1.d0/Elem(ElemNo)%L
			Elem(ElemNo)%EPS=0.d0; Elem(ElemNo)%SIG=0.d0;
			Elem(ElemNo)%dEPS=0.d0; Elem(ElemNo)%STRESS=0.d0;
			Elem(ElemNo)%fy_temp=Elem(ElemNo)%fy; Elem(ElemNo)%Et=Elem(ElemNo)%E
			Elem(ElemNo)%R=0.d0;
		end do
		return
	end subroutine Intial_Elem
end module TypDef

module Main_Iteration ! 主迭代模块
	use TypDef
	implicit none

	contains
	
	subroutine Iteration(GValue, Node, Elem, Load, Support) !主迭代程序
		type(typ_GValue) :: GValue
		type(typ_Node) :: Node(:)
		type(typ_Truss) :: Elem(:)
		type(typ_Load) :: Load(:)
		type(typ_Support) :: Support(:)
		real(rkind) :: GLoad(2*GValue%NNode), GDisp(2*GValue%NNode) ! 总荷载向量，总位移向量
		real(rkind) :: GK(2*GValue%NNode,2*GValue%NNode) ! 总刚度矩阵
		real(rkind) :: RLoad(2*GValue%NNode) ! 节点反力
		integer (ikind) :: NCycle ! 迭代次数
		integer (ikind) :: i
		NCycle=1
		do
			print *, "    第",NCycle,"次迭代"
		
			call Get_GLoad(GValue,Load, GLoad) ! 得到当前荷载
			call Get_RLoad(GValue, Elem, RLoad) ! 得到节点反力
			GLoad=GLoad-RLoad ! 得到节点不平衡力
			do i=1, GValue%NElem;	call Get_EK(Elem(I));		end do ! 得到单元刚度矩阵
			call Get_GK (GValue, Elem, GK)  ! 得到整体刚度矩阵
			call Get_Support(GValue, GK, GLoad, Support) ! 组装支座信息
			call Check_Error(GValue, Load, GLoad) ! 检查是否收敛
			print *, "    当前误差为",GValue%Error, "误差容限为", GValue%Tol
			if(GValue%Error<=GValue%Tol) then
				print *, "    当前步迭代收敛"
				Elem(:)%EPS=Elem(:)%EPS+Elem(:)%dEPS ! 更新单元应变
				Elem(:)%SIG=Elem(:)%STRESS           ! 更新单元应力
				Elem(:)%fy=Elem(:)%fy_temp           ! 更新屈服应力
				do i=1,size(Node)
					Node(i)%Disp=Node(i)%Disp+Node(i)%dDisp ! 更新节点位移
					Node(i)%dDisp=0.d0                      ! 节点位移增量清零
				end do
				exit !  退出迭代
			end if
			GK=matinv(GK) ! 刚度矩阵求逆
			GDisp=matmul(GK,GLoad)  ! 得到位移增量
			call Get_dDisp(GValue, Node, GDisp) ! 得到位移增量
			do i=1, GValue%NElem;	call Get_EPS(Node,Elem(I));	end do ! 得到应变增量
			do i=1, GValue%NElem;	call Get_Stress(Elem(I));	end do ! 得到当前应力
			do i=1, GValue%NElem;	call Get_R(Elem(I));		end do ! 得到节点反力
			NCycle=NCycle+1
		end do
		return
	end subroutine Iteration
	
	subroutine Check_Error(GValue, Load, GLoad) ! 检查是否收敛，采用无穷范数 e<|R|/|Load|
		type(typ_GValue) :: GValue
		type(typ_Load) :: Load(:)
		real(rkind) :: GLoad(:)
		real(rkind) :: x,y
		integer (ikind) :: i
		x=maxval(abs(GLoad)) ! 得到不平衡力的最大值
		y=0.d0
		do i=1,size(Load) 
			y=max(y,maxval(abs(Load(i)%Value))/real(GValue%NStep))	 
		end do
		GValue%Error=x/y ! 得到相对误差
		return
	end subroutine

	subroutine Get_Support(GValue, GK, GLoad, Support) ! 采用化零置一法集成支座信息
		type(typ_GValue) :: GValue
		type(typ_Support) :: Support(:)
		real(rkind) :: GK(:,:),GLoad(:)
		integer(ikind) :: i,j
		do i=1,size(Support)
			j=Support(i)%NodeNo*2-2+Support(i)%DOF ! 得到支座对应的自由度编号
			GK(j,:)=0.d0
			GK(:,j)=0.d0
			GK(j,j)=1;
			GLoad(j)=0.d0
		end do
		return
	end subroutine Get_Support

	subroutine Get_GK (GValue, Elem, GK) ! 组装整体刚度矩阵
		type(typ_GValue) :: GValue
		type(typ_Truss) :: Elem(:)
		real(rkind) :: GK(:,:)
		integer(ikind) :: i,j,k,N(4)
		GK=0.d0;
		do i=1, size(Elem)
			N=(/Elem(i)%NodeNo(1)*2-1, Elem(i)%NodeNo(1)*2, Elem(i)%NodeNo(2)*2-1, &
				Elem(i)%NodeNo(2)*2/)
			do j=1,4
				do k=1,4
					GK(N(j),N(k))=GK(N(j),N(k))+Elem(i)%K(j,k)
				end do
			end do
		end do
		return
	end subroutine Get_GK
	
	subroutine Get_RLoad(GValue, Elem, RLoad) ! 得到结构反力
		type(typ_GValue) :: GValue 
		type(typ_Truss) :: Elem(:)
		real(rkind) :: RLoad(:)
		real(rkind) :: R(4)
		integer (ikind) :: i
		RLoad=0.d0
		do i=1,size(Elem)
			R=Elem(i)%R
			RLoad(Elem(i)%NodeNo(1)*2-1)=RLoad(Elem(i)%NodeNo(1)*2-1)+R(1)
			RLoad(Elem(i)%NodeNo(1)*2  )=RLoad(Elem(i)%NodeNo(1)*2  )+R(2)
			RLoad(Elem(i)%NodeNo(2)*2-1)=RLoad(Elem(i)%NodeNo(2)*2-1)+R(3)
			RLoad(Elem(i)%NodeNo(2)*2  )=RLoad(Elem(i)%NodeNo(2)*2  )+R(4)
		end do		
		return
	end subroutine Get_RLoad

	subroutine Get_GLoad(GValue,Load, GLoad) ! 得到总体外荷载
		type(typ_GValue) :: GValue
		type(typ_Load) :: Load(:)
		real(rkind) :: GLoad(:)
		integer (ikind) :: i
		GLoad=0.d0
		do i=1, size(Load)
			GLoad(Load(i)%NodeNo*2-1)=GLoad(Load(i)%NodeNo*2-1)+Load(i)%Value(1)
			GLoad(Load(i)%NodeNo*2  )=GLoad(Load(i)%NodeNo*2  )+Load(i)%Value(2)
		end do
		GLoad=GLoad*real(GValue%CurrentStep)/real(GValue%NStep)
		return
	end subroutine Get_GLoad

	function matinv(A) result (B) ! 对矩阵求逆
		real(rkind) ,intent (in)::A(:,:)
		real(rkind) , pointer ::B(:,:)
		integer(ikind):: N,I,J,K;	real(rkind)::D,T
		integer(ikind), allocatable::IS(:),JS(:)
		N=size(A,dim=2); allocate(B(N,N))
		allocate(IS(N));allocate(JS(N))
		B=A
		do 	K=1,N
			D=0.0D0
			do I=K,N
				do J=K,N
					if(abs(B(I,J))>D) then
						D=abs(B(I,J)); IS(K)=I; JS(K)=J
					end if
				end do
			end do
			do J=1,N
				T=B(K,J); B(K,J)=B(IS(K),J); B(IS(K),J)=T
			end do
			do I=1,N
				T=B(I,K);	B(I,K)=B(I,JS(K));	B(I,JS(K))=T
			end do
			B(K,K)=1.d0/B(K,K);
			do J=1,N; if(J.NE.K) B(K,J)=B(K,J)*B(K,K);  end do
			do I=1,N
				if(I.NE.K) then
					do J=1,N; if(J.NE.K) B(I,J)=B(I,J)-B(I,K)*B(K,J); end do
				end if
			end do
			do I=1,N; if(I.NE.K) B(I,K)=-B(I,K)*B(K,K);	end do
		end do
		do K=N,1,-1
			do J=1,N
				T=B(K,J); B(K,J)=B(JS(K),J); B(JS(K),J)=T
			end do
			do I=1,N
				T=B(I,K); B(I,K)=B(I,IS(K)); B(I,IS(K))=T
			end do
		end do
		return
	end function matinv

end module Main_Iteration


module Data_Input
	use TypDef
	implicit none

	contains

	subroutine DataInput(GValue, Node, Elem, Load, Support) ! 读入数据
		type(typ_GValue) :: GValue
		type(typ_Node),pointer :: Node(:)
		type(typ_Truss),pointer :: Elem(:)
		type(typ_Load),pointer :: Load(:)
		type(typ_Support),pointer :: Support(:)
		integer (ikind) :: i,j
		open(55,file='input.txt')
			read(55,*)  GValue%NStep, GValue%Tol ! 读入荷载步数，误差容限
			read(55,*)  GValue%NNode ! 读入节点数目
			allocate(Node(GValue%NNode))
			do i=1, GValue%NNode
				read(55,*) j, Node(i)%Coord ! 读入节点坐标
				Node(i)%Disp=0.d0; Node(i)%dDisp=0.d0
			end do
			
			read(55,*) GValue%NElem ! 读入单元数目
			allocate(Elem(GValue%NElem))
			do i=1, GValue%NElem
				read(55,*) j, Elem(i)%NodeNo(1), Elem(i)%NodeNo(2), Elem(i)%A, &
					Elem(i)%E, Elem(i)%fy, Elem(i)%Eh 
					! 读入单元节点编号，单元截面积，初始弹性模量，屈服强度，硬化模量。
			end do
			call Intial_Elem(Node,Elem) !赋予单元初始条件

			read(55,*) GValue%NLoad ! 读入节点荷载数
			allocate(Load(GValue%NLoad)) 
			do i=1, GValue%NLoad
				read(55,*) j, Load(i)%NodeNo, Load(i)%Value(1:2) ! 读入节点荷载
			end do

			read(55,*) GValue%NSupport ! 读入支座数
			allocate(Support(GValue%NSupport))
			do i=1, GValue%NSupport
				read(55,*) j, Support(i)%NodeNo, Support(i)%DOF ! 读入节点，支座约束自由度
			end do

		close(55)
		return
	end subroutine DataInput

end module Data_Input

module Data_Output
	use TypDef
	implicit none

	contains
	subroutine DataOutput(GValue, Node, Elem, Load, Support) ! 计算结果输出模块
		type(typ_GValue) :: GValue
		type(typ_Node),pointer :: Node(:)
		type(typ_Truss),pointer :: Elem(:)
		type(typ_Load),pointer :: Load(:)
		type(typ_Support),pointer :: Support(:)
		integer (ikind) :: i
		open(55,file="dispout.txt")
		do i=1, size(Node)
			write(55,'(I3,2G13.5)') i, Node(i)%Disp(1:2) ! 输出最终节点位移
		end do
		close(55)
		open(55,file="Elemout.txt")
		do i=1,size(Elem)
			write(55,'(I3,2G13.5)') i, Elem(i)%EPS, Elem(i)%SIG ! 输出最终单元应变，应力
		end do
		close(55)
		return
	end subroutine DataOutput
end module Data_Output

program Main
	use TypDef
	use Main_Iteration
	use Data_Input
	use Data_Output
	type(typ_GValue) :: GValue
	type(typ_Node),pointer :: Node(:)
	type(typ_Truss),pointer :: Elem(:)
	type(typ_Load),pointer :: Load(:)
	type(typ_Support),pointer :: Support(:)
	integer (ikind) :: Inc ! 增量步
	
	call DataInput(GValue, Node, Elem, Load, Support) ! 读入数据
	do Inc=1, GValue%NStep
		print *,"****************************************************"
		print *,"当前增量步", inc
		GValue%CurrentStep=Inc
		call Iteration(GValue, Node, Elem, Load, Support) ! 进行迭代
	end do
	call DataOutput(GValue, Node, Elem, Load, Support) ! 输出结果
	stop
end program