基于Newton-Raphson法的平面桁架弹塑性分析有限元程序.docVIP

程序说明 本程序为Newton-Raphson法求解非线性问题的范例。为了便于读者掌握非线性程序编写的基本思路，本程序设计时着重以下特点： 计算模型尽可能简单，采用最简单的单元——平面桁架单元，和最简单的非线性本构关系——双线性弹塑性本构关系。以减少读者学习难度，充分提高程序的可读性和实用性。 尽可能多些注释，尽量符合矩阵运算法则，减少程序长度，便于读者理解。 本程序执行的基本流程为 上述核心迭代过程相应程序代码如下: 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 算 例 分析下图所示平面桁架，构件截面积A均为1，弹性模量E为200E3，屈服强度fy为200，硬化模量Eh为2E3。 程序数据输入格式 第一行，输入荷载增量步数，收敛标准 例如：10, 0.001 第二行开始，输入节点数目和节点编号和坐标 8 1 ,0.0000, 0.0000 2 ,1