📄 2d_8nodes.f90
字号:
!=============================================================
! 平面8节点等参元程序
! 清华大学土木工程系：陆新征
!============================================================
module Elem_Rect8 ! 八节点等参元
implicit none
integer (kind(1)),parameter ::ikind=(kind(1))
integer (kind(1)),parameter ::rkind=(kind(0.d0))
type :: typ_Kcol
real(rkind),pointer :: Row(:)
end type typ_Kcol
type :: typ_GValue !总体控制变量
integer(ikind) :: NNode, NElem, NLoad, NMat, NSupport
integer(ikind) :: NGlbDOF !整体自由度总数
integer(ikind) :: NGENS, NodeDOF,ElemNodeNo
integer(ikind) :: NInt
end type typ_GValue
type Typ_Node !定义节点类型
real(rkind) :: coord(2) !节点坐标
integer(ikind) :: GDOF(2) !整体自由度编码
real(rkind) :: DISP(2) !节点位移
real(rkind) :: dDISP(2) !节点位移增量
real(rkind) :: dForce(2) !节点不平衡力
end type typ_Node
!=============================================================================
Type typ_IntPoint !定义积分点参数
real(rkind) :: EPS(3) !应变
real(rkind) :: SIG(3) !应力
real(rkind) :: D(3,3) !本构矩阵
real(rkind) :: B(3,16) !几何矩阵
real(rkind) :: DETJ !雅克比行列式
end type Typ_IntPoint
type Typ_Rect8 !定义实体单元
integer(ikind) :: NodeNo(8) !节点编号
real(rkind) :: E !弹性模量
real(rkind) :: u !泊松比
real(rkind) :: t !单元厚度
real(rkind) :: EK(16,16) !单元刚度矩阵
type(typ_intpoint) :: IntP(9) !积分点
!........................
end type Typ_Rect8
type Typ_Load ! 定义荷载类型
integer(ikind) :: NodeNo !荷载作用节点编号
real(rkind) :: Value(2) !荷载大小
end type Typ_Load
type Typ_Support ! 定义支座类型
integer(ikind) :: NodeNo !支座节点编号
integer(ikind) :: DOF !支座约束自由度
real(rkind) :: Value !支座位移大小
end type Typ_Support
contains
subroutine Get_Elem_K(GValue,Elem,Node) ! 核心程序，得到单元的刚度矩阵
implicit none
type(typ_GValue) :: GValue
type(typ_Node) :: Node(:)
type(typ_Rect8) :: Elem(:)
real*8 :: InvJ(2,2) ! 雅克比矩阵及其逆矩阵
real*8 :: Coord(8,2),IntPCoord(9,2),IntPwt(9)
real*8 :: dN(2,8) !节点坐标，积分点坐标，权函数，形函数导数
integer :: I,ElemNo
call Get_IntP_Prop(IntPCoord,IntPWt) ! 计算积分点信息
do ElemNo=1,size(Elem)
Elem(ElemNo)%EK=0.0
do I=1, GValue%ElemNodeNo; !得到节点坐标
coord(I,:)=Node(Elem(ElemNo)%NodeNo(i))%coord;
end do
Elem(ElemNo)%EK=0d0
DO I=1,size(Elem(ElemNo)%IntP)
! 计算本构矩阵
call Get_D(Elem(ElemNo)%IntP(I)%D,Elem(ElemNo)%E,Elem(ElemNo)%u)
! 计算形函数对局部坐标导数
call Get_dN_dxi(dN,IntPCoord(i,1),IntPCoord(i,2))
InvJ=matmul(dN, Coord)
! 得到雅克比行列式
Elem(ElemNo)%IntP(I)%DETJ=InvJ(1,1)*InvJ(2,2)-InvJ(1,2)*InvJ(2,1)
InvJ=matinv(InvJ); ! 对雅克比行列式求逆
dN = matmul(InvJ,dN) ! 得到形函数对整体坐标的导数
call Get_B (Elem(ElemNo)%IntP(I)%B,dN) ! 得到几何矩阵
Elem(ElemNo)%EK=Elem(ElemNo)%EK+&
matmul(matmul(transpose(Elem(ElemNo)%IntP(I)%B),&
Elem(ElemNo)%IntP(I)%D),Elem(ElemNo)%IntP(I)%B)*&
Elem(ElemNo)%IntP(I)%DetJ*IntPWt(i)*Elem(ElemNo)%t
end do
end do
return
end subroutine
subroutine Get_D(D,e,v) ! 得到本构矩阵，平面应力
implicit none
real*8 :: e,v
real*8 :: D (:,:)
D=0.d0
D(1,1)=1D0; D(2,2)=1d0
D(1,2)=v; D(2,1)=v
D(3,1:2)=0d0; D(1:2,3)=0d0;
D(3,3)=(1.d0-v)/2.d0;
D=E/(1.d0-v*v)*D;
return
end subroutine Get_D
subroutine Get_IntP_Prop(IntPCoord,IntPWt) ! 得到高斯积分点的参数
implicit none
real*8 ::IntPCoord(:,:),IntPWt(:) ! 高斯积分点坐标，高斯积分点权重
real*8 :: z,x(3),y(3),w(3)
integer :: i,j
z=.2d0*sqrt(15.d0)
x=(/-1.d0,0.d0,1.d0/); y=(/1.d0,0.d0,-1.d0/)
w=(/5.d0/9.d0,8.d0/9.d0,5.d0/9.d0/)
do i=1,3
do j=1,3
IntPCoord((i-1)*3+j,1)=x(j)*z
IntPCoord((i-1)*3+j,2)=y(i)*z
IntPWt((i-1)*3+j)=w(i)*w(j)
end do
end do
return
end subroutine Get_IntP_Prop
subroutine Get_dN_dxi(dN_dxi,xi,eta) ! 得到形函数对局部坐标系的导数
implicit none
real*8 :: dN_dxi(:,:), xi,eta
real*8:: x1, x2, x3, x4, x5 ,x6 ,x7 ,x8
x1=.25*(1.-eta);
x2=.25*(1.+eta);
x3=.25*(1.-xi);
x4=.25*(1.+xi)
dN_dxi(1,1)=x1*(2.*xi+eta)
dN_dxi(1,2)=-8.*x1*x2
dN_dxi(1,3)=x2*(2.*xi-eta)
dN_dxi(1,4)=-4.*x2*xi
dN_dxi(1,5)=x2*(2.*xi+eta)
dN_dxi(1,6)=8.*x2*x1
dN_dxi(1,7)=x1*(2.*xi-eta)
dN_dxi(1,8)=-4.*x1*xi
dN_dxi(2,1)=x3*(xi+2.*eta)
dN_dxi(2,2)=-4.*x3*eta
dN_dxi(2,3)=x3*(2.*eta-xi)
dN_dxi(2,4)=8.*x3*x4
dN_dxi(2,5)=x4*(xi+2.*eta)
dN_dxi(2,6)=-4.*x4*eta
dN_dxi(2,7)=x4*(2.*eta-xi)
dN_dxi(2,8)=-8.*x3*x4
return
end subroutine Get_dN_dxi
subroutine Get_B(B,dN_dx) ! 得到几何矩阵
implicit none
real*8 :: B(:,:), dN_dx(:,:)
integer::k,l,m,n , ih,nod;
real*8 :: x,y,z
B=0.
do m=1,8
k=2*m; l=k-1; x=dN_dx(1,m); y=dN_dx(2,m)
B(1,l)=x; B(3,k)=x; B(2,k)=y; B(3,l)=y
end do
return
end subroutine Get_B
function matinv(A) result (B) ! 计算2x2逆矩阵
real(rkind) ,intent (in):: A(:,:)
real(rkind) , pointer::B(:,:)
real(rkind) :: x
integer :: N
N=size(A,dim=2)
allocate(B(N,N))
B(1,1)=A(2,2)
B(2,2)=A(1,1)
B(1,2)=-A(1,2)
B(2,1)=-A(2,1)
x=A(1,1)*A(2,2)-A(1,2)*A(2,1)
B=B/x
return
end function matinv
end module
module Data_Input ! 数据输入模块
use Elem_Rect8
implicit none
contains
subroutine DataInput(GValue,Node,Elem,Load,Support)
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
type(typ_Rect8),pointer :: Elem(:)
type(typ_Load),pointer :: Load(:)
type(typ_Support),pointer :: Support(:)
call ReadGValue(GValue) !读入整体控制变量
call ReadNode(GValue,Node) !读入节点坐标信息
call ReadElem(GValue,Elem) !读入单元原始信息
call ReadMaterial(GValue,Elem)
call ReadLoad(GValue,Load) !读入荷载信息
call ReadSupport(GValue,Support) !读入支座信息
return
end subroutine DataInput
subroutine ReadGValue(GValue)
type(typ_GValue) :: GValue
GValue%NGENS=3
GValue%NodeDOF=2
GValue%ElemNodeNo=8
GValue%NInt=3
return
end subroutine ReadGValue
subroutine ReadNode(GValue,Node) ! 读入结点坐标
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
integer(ikind) :: I,J
open(55,file='node.txt')
read(55,*) GValue%NNode
GValue%NGlbDOF=2*GValue%NNode
allocate(Node(GValue%NNode))
do I=1, GValue%NNode
read(55,*) J,Node(I)%Coord(1:2)
end do
close(55)
write(*,*) " 结点信息读入完毕"
return
end subroutine ReadNode
subroutine ReadElem(GValue,Elem) ! 读入单元信息
type(typ_GValue) :: GValue
type(typ_Rect8),pointer :: Elem(:)
integer(ikind) :: I,J,K
integer(ikind) :: TransNode(8),N(8),N1(8) ! 结点坐标变换
TransNode=(/1,7,5,3,8, 6,4, 2/) ! 注意单元节点编号匹配
open(55,file='Element.txt')
read(55,*) GValue%NElem
allocate(Elem(GValue%NElem))
do I=1,GValue%NElem
read(55,*) J, N
do J=1,size(N)
K=TransNode(J)
N1(K)=N(J)
end do
Elem(I)%NodeNo=N1
end do
close(55)
write(*,*) " 单元信息读入完毕"
return
end subroutine ReadElem
subroutine ReadMaterial(GValue,Elem) ! 读入材料信息
type(typ_GValue) :: GValue
type(typ_Rect8),pointer :: Elem(:)
integer (ikind) :: NElem
real(rkind) :: E, mu
integer(ikind) :: I,J
integer(ikind),allocatable :: K(:)
open(55,file='Material.txt')
read(55,*) E, mu
Elem(:)%E=E
Elem(:)%u=mu
Elem(:)%t=0.5
close(55)
write(*,*) " 材料信息读入完毕"
return
end subroutine ReadMaterial
subroutine ReadLoad(GValue,Load) ! 读入荷载信息
type(typ_GValue) :: GValue
type(typ_Load),pointer :: Load(:)
integer(ikind) :: I,J
open(55,file='Load.txt')
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
close(55)
write(*,*) " 荷载信息读入完毕"
return
end subroutine ReadLoad
subroutine ReadSupport(GValue,Support) ! 读入支座信息
type(typ_GValue) :: GValue
type(typ_Support),pointer :: Support(:)
integer(ikind) :: I,J
open(55,file='Support.txt')
read(55,*) GValue%NSupport
allocate(Support(GValue%NSupport))
do I=1, GValue%NSupport
read(55,*) J, Support(I)%NodeNo, Support(I)%DOF, Support(I)%Value
end do
close(55)
write(*,*) " 支座信息读入完毕"
return
end subroutine ReadSupport
end module
module Mat_solve ! 矩阵求解模块
use Elem_Rect8
implicit none
integer :: NGlbDOF
contains
subroutine Matsolve(GValue,Elem,Node,Load, Support)
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
type(typ_Rect8),pointer :: Elem(:)
type(typ_Load),pointer :: Load(:)
type(typ_Support),pointer :: Support(:)
type(typ_Kcol),allocatable :: Kcol(:)
real(rkind),allocatable :: GLoad(:), GDisp(:)
real(rkind) :: Penatly
integer(ikind) :: BandWidth(2*GValue%NNode)
integer(ikind) :: ELocVec(16)
integer(ikind) :: I,J,K
integer(ikind) :: MinDOFNum
NGlbDOF=2*GValue%NNode
Penatly=1.0 ! 罚函数
allocate(GLoad(NGlbDOF))
allocate(GDisp(NGlbDOF))
!查找带宽
do I=1,NGlbDOF
BandWidth(I)=I
end do
do I=1,GValue%NElem
do J=1,size(Elem(I)%NodeNo)
ELocVec(J*2-1)=2*Elem(I)%NodeNo(J)-1
ELocVec(J*2 )=2*Elem(I)%NodeNo(J)
end do
MinDOFNum=minval(ELocVec)
do J=1,size(ElocVec)
BandWidth(ELocVec(J))=min(MinDOFNum,BandWidth(ELocVec(J)))
end do
end do
write(*,*) " 完成带宽查找"
allocate(Kcol(NGlbDOF))
do I=1, NGlbDOF
allocate(Kcol(I)%Row(BandWidth(I):I))
Kcol(I)%Row=0.d0
end do
do I=1, GValue%NElem
do J=1,size(Elem(I)%NodeNo)
ELocVec(J*2-1)=2*Elem(I)%NodeNo(J)-1
ELocVec(J*2 )=2*Elem(I)%NodeNo(J)
end do
do J=1,size(ElocVec)
do K=1,J
if(ELocVec(J)>ELocVec(K)) then
Kcol(ELocVec(J))%row(ELocVec(K))=&
Kcol(ELocVec(J))%row(ELocVec(K))+Elem(I)%EK(J,K)
else
Kcol(ELocVec(K))%row(ELocVec(J))=&
Kcol(ELocVec(K))%row(ELocVec(J))+Elem(I)%EK(J,K)
end if
end do
end do
Penatly=max(Penatly,maxval(abs(Elem(I)%EK)))
end do
write(*,*) " 完成总刚集成"
GLoad=0.d0
do I=1, size(Load)
J=Load(I)%NodeNo
GLoad(J*2-1:J*2)=GLoad(J*2-1:J*2)+Load(I)%Value(:)
end do
do I=1,GValue%NSupport
J=2*(Support(I)%NodeNo-1)+Support(I)%DOF
GLoad(J)=GLoad(J)+Support(I)%Value*Penatly*1.0D8
Kcol(J)%Row(J)=KCol(J)%Row(J)+Penatly*1.0d8
end do
write(*,*) " 完成支座集成"
call BandSolv(Kcol,Gload,GDisp)
call Get_Node_dDisp(GValue,Node,GDisp)
do I=1,NGlbDOF
deallocate(Kcol(I)%Row)
end do
deallocate(Kcol)
deallocate(GDisp)
deallocate(GLoad)
return
end subroutine
subroutine Get_Node_dDisp(GValue,Node,GDisp) ! 从整体位移向量得到结点位移向量
type(typ_GValue) :: GValue
type(typ_Node) :: Node(:)
real(rkind) :: GDisp(:)
integer(ikind) :: I,J
do I=1, GValue%NNode
Node(I)%dDisp(1:2)=GDisp(I*2-1:I*2)
Node(I)%Disp=Node(I)%Disp+Node(I)%dDisp
end do
return
end subroutine Get_Node_dDisp
!---------------------------------------
subroutine BandSolv(Kcol,GLoad,diag)
!---------------------------------------
type (typ_Kcol) :: Kcol(:);
real(rkind) :: GLoad(:),diag(:);
integer :: row1,ncol,row,j,ie
real(rkind) :: 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
diag(:)=GLoad(:)
return;
end subroutine
end module
module Data_Output ! 数据输出模块
use Elem_Rect8
implicit none
contains
subroutine DataOutput(GValue,Elem,Node,Load, Support) ! 输出节点位移
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
type(typ_Rect8),pointer :: Elem(:)
type(typ_Load),pointer :: Load(:)
type(typ_Support),pointer :: Support(:)
integer i
open (55,file='out.txt')
do i=1,size(Node)
write(55,'(I5, 2G14.5)') i, Node(i)%Disp(1:2)
end do
close(55)
return
end subroutine
end module
Program Main ! 主程序
use Elem_Rect8
use Data_Input
use Mat_Solve
use Data_Output
implicit none
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
type(typ_Rect8),pointer :: Elem(:)
type(typ_Load),pointer :: Load(:)
type(typ_Support),pointer :: Support(:)
write(*,*) "开始读取数据"
call DataInput(GValue,Node,Elem,Load,Support) ! 读入数据
write(*,*) "开始计算单元刚度矩阵"
call Get_Elem_K(GValue,Elem,Node) ! 得到单元刚度矩阵
write(*,*) "开始求解整体方程"
call MatSolve(GValue,Elem,Node,Load, Support) ! 求解整体平衡方程
write(*,*) "开始输出结果"
call DataOutput(GValue,Elem,Node,Load, Support) ! 输出结果
stop
end program
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -