DOC《程序说明》

程序说明 本程序为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