C =====================================================
C 清华大学研究生精品课程《钢筋混凝土有限元》教学程序
C 双折线弹塑性梁单元本构程序
C User defined nonlinear beam subroutine for MSC.MARC, Fortran Source code
C 陆新征 江见鲸
C Lu Xinzheng, Jiang Jianjing
C 清华大学土木工程系，北京，100084
C Dept. Civil Engrg. of Tsinghua University
C luxinzheng@sina.com.cn
C last revised: Sep. 2003.
C 主要参考文献：
C 1. MARC Volumn D user subroutine and Special Routines"
C 2. 江见鲸 "钢筋混凝土结构非线性有限元分析"
C =====================================================

C -----------------------------------
C 公共变量模块，中间数据存放与此
C 格式为：单元号(要求大于模型中最大单元编号)，积分点号(大于最大积分点数)，数据块
C ------------------------------------
module My_Var
real*8 My_Data(1000,10,33)
end module

subroutine ubeam(d,fcrp,df,dum1,etot,de,dum2,s,dum3,gs,dum4,temp,
* dtemp,ngens,m,nnn,mat)
use My_Var
implicit real*8 (a-h,o-z) dp
include "../common/concom"
dimension d(ngens,ngens),fcrp(*),df(ngens),etot(ngens),de(ngens),
* s(*),gs(ngens),temp(*),dtemp(*)
dimension statev(33)
statev=My_Data(m,nnn,:)
c
c 基本变量性质说明
C statev(1) : EA0 !初始压拉刚度
C statev(2) : EIxx0 !初始xx方向弯曲刚度
C statev(3) : EIyy0 !初始yy方向刚度
C statev(4) : EIxy0 !初始扭转刚度
C
C statev(5) : F_p_y !轴向抗拉屈服强度
C statev(6) : F_n_y !轴向抗压屈服强度
C statev(7) : Mxx_p_y!xx正方向弯曲强度
C statev(8) : Mxx_n_y!xx负方向弯曲强度
C statev(9) : Myy_p_y!yy正方向弯曲强度
C statev(10): Myy_n_y!yy负方向弯曲强度
C statev(11): Mxy_p_y!xy正方向弯曲强度
C statev(12): Mxy_n_y!xy负方向弯曲强度
C
C statev(13): EA_s_p !轴向抗拉强化刚度
C statev(14): EA_s_y !轴向抗压强化刚度
C statev(15): EIxx_s_p!xx正方向强化刚度
C statev(16): EIxx_s_n!xx负方向强化刚度
C statev(17): EIyy_s_p!yy正方向强化刚度
C statev(18): EIyy_s_n!yy负方向强化刚度
C statev(19): EIxy_s_p!xy正方向强化刚度
C statev(20): EIxy_s_n!xy负方向强化刚度
C
C statev(21): F_p_m !轴向抗拉最大荷载(abs)
C statev(22): F_n_m !轴向抗压最大荷载(abs)
C statev(23): Mxx_p_m!xx正方向最大荷载(abs)
C statev(24): Mxx_n_m!xx负方向最大荷载(abs)
C statev(25): Myy_p_m!yy正方向最大荷载(abs)
C statev(26): Mxx_n_m!yy负方向最大荷载(abs)
C statev(27): Mxy_p_m!xy正方向最大荷载(abs)
C statev(28): Mxy_n_m!xy负方向最大荷载(abs)
C
C statev(29): Former_INC !前一次迭代步数
C statev(30): Former_F !前一次迭代轴力
C statev(31): Former_Mxx !前一次迭代xx弯矩
C statev(32): Former_Myy !前一次迭代yy弯矩
C statev(33): Former_Mxy !前一次迭代xy扭矩
C

if(inc==0.and.ncycle==0.and.incsub==0) then ! 如果是第一次计算，变量赋初值
if(mat==1) then !如果是杆件类型1
B=0.5; H=1.0; Ec=30.d9; fc=30e6
c 截面初始刚度
statev(1)=B*H*Ec !B*H*Ec
statev(2)=B**3*H*Ec/12. !B**3*H*Ec/12
statev(3)=B*H**3*Ec/12. !B*H**3*Ec/12
statev(4)=statev(2)+statev(3)
C 截面屈服强度
statev(5)=B*H*fc/2.
statev(6)=B*H*0.5*fc/2.
statev(7)=7./24.*(0.1*fc)*B**2*H*3.
statev(8)=7./24.*(0.1*fc)*B**2*H*3.
statev(9)=7./24.*(0.1*fc)*B*H**2*3.
statev(10)=7./24.*(0.1*fc)*B*H**2*6.
statev(11)=1.0d10
statev(12)=1.0d10
C 截面屈服后刚度
soft_factor=.01
statev(13)=statev(1)*soft_factor;
statev(14)=statev(1)*soft_factor;
statev(15)=statev(2)*soft_factor;
statev(16)=statev(2)*soft_factor;
statev(17)=statev(3)*soft_factor;
statev(18)=statev(3)*soft_factor;
statev(19)=statev(4)*soft_factor;
statev(20)=statev(4)*soft_factor;
C 截面初始最大荷载假设为屈服荷载
statev(21:28)=statev(5:12);
C 变量清零
statev(29:33)=0.
end if
end if

C
C 在新的加载步开始，判断是否需要更新历史最大截面内力
C
if(inc.gt.statev(29)) then !新的加载步开始
statev(29)=inc
do i=1,ngens
C 如果上一步计算结束的荷载大于历史最大荷载，则更新历史最大荷载
if(statev(29+i)>=0.and.
1 abs(statev(29+i))>
1 statev(20+(I-1)*2+1)) then
statev(20+(I-1)*2+1)=abs(statev(29+i))
end if
if(statev(29+i)<0.and.
1 abs(statev(29+i))>
1 statev(20+(I-1)*2+2)) then
statev(20+(I-1)*2+2)=abs(statev(29+i))
end if
end do
end if

d=0.

do I=1, ngens
if (gs(I)*de(I).le.0.) then !卸载
d(I,I)=statev(I) !卸载刚度等于初始刚度
else !加载或再加载
if(gs(I).ge.0) then !正荷载情况
if(abs(gs(I)+de(i)*statev(i)).ge.
1 statev(20+(I-1)*2+1)) then ! 屈服后骨架线加载
if(abs(gs(i))<statev(20+(I-1)*2+1)) then
ratio=(statev(20+(I-1)*2+1)-abs(gs(i)))/
1 abs(de(i)*statev(i))
d(i,i)=statev(i)*ratio+
1 (1.-ratio)*statev(12+(I-1)*2+1) ! 刚度为屈服后刚度
else
d(i,i)=statev(12+(I-1)*2+1)
end if
else
d(I,I)=statev(I) !刚度为初始刚度
end if
else
if(abs(gs(I)+de(i)*statev(i)).ge.
1 statev(20+(I-1)*2+2)) then ! 屈服后骨架线加载
if(abs(gs(i))<statev(20+(I-1)*2+2)) then
ratio=(statev(20+(I-1)*2+2)-abs(gs(i)))/
1 abs(de(i)*statev(i))
d(i,i)=statev(i)*ratio+
1 (1.-ratio)*statev(12+(I-1)*2+2) ! 刚度为屈服后刚度
else
d(i,i)=statev(12+(I-1)*2+2)
end if
else
d(I,I)=statev(I) !刚度为初始刚度
end if
end if
end if
df(i)=d(i,i)*de(i) !荷载增量等于切线刚度乘以变形增量
gs(i)=gs(i)+df(i) !更新荷载向量
statev(29+i)=(gs(i)) !记录此时的荷载与内存区块
end do

C 为提高收敛性，修正切线刚度矩阵
do I=1,ngens
if(d(I,I)<statev(I))
1 d(i,i)=0.5*d(I,I)+0.5*statev(I)
end do
My_Data(m,nnn,:)=statev
return
end