基于Fortran的ABAQUS子程序TJU.Plastic-Explicit |
 |
|||||
Fortran语言的ABAQUS子程序TJU.plastic-explicit
（Fortran 语言 / 配合TJU.SAP2ABAQUS接口导出的inp文件） uploaded to www.liuguohuan.net by 柳 天津大学 建筑工程学院 ； 天津大学 前沿技术研究院； 国家重点实验室-水利工程仿真与安全 更多程序下载
!! -------------------------- 定义公共----------------------------
module BHM implicit none type :: BHM_material ! 定义滞回模型相关的材料参数,派生类名 real*8 :: stiff_0 ! 第一刚度（初始刚度），结构体成员 real*8 :: ratio ! 屈服刚度比（第二刚度/初始刚度） real*8 :: dis_yie ! 屈服位移 real*8 :: forc_yie ! 屈服强度（力） end type BHM_material type :: BHM_spring ! 定义与弹簧编号相关的数组（行数与spring的number对应，列数与弹簧需要输入的属性个数对应） real(kind=8) DIMENSION(200,8):: atrr_spr ! 定义与弹簧编号对应的属性数组【弹簧数（维数）>最大编号，这里取200】，列数足够大即可,这里取8 end type BHM_spring type :: BHM_update ! 定义每一步 需要 更新和保存的变量（行数与spring的number对应，列数与步数相关） real(kind=8) DIMENSION(200,5000):: dis_cur ! 当前步 相对位移 按照40秒的地震波计算,当0.01秒间隔时为4001步,所以列数 5000 已足够 real(kind=8) DIMENSION(200,5000):: dis_inc ! 当前步 相对位移的 增量 real(kind=8) DIMENSION(200,5000):: forc_cur ! 当前步 滞回力 real(kind=8) DIMENSION(200,5000):: dis_save ! 每一步 需要保存的相对位移， real(kind=8) DIMENSION(200,5000):: forc_save ! 每一步 需要保存的滞回力 real(kind=8) DIMENSION(200,5000):: dis_cro ! 每一步 需要保存的交叉点位移，加载时需要 end type BHM_update contains !! -------------------------- 定义公共块 -------------------------------- !! -------------------------- 《《《 BHM 子程序开始 》》》 ----------------- SUBROUTINE USPRNG(RATK,F,DATAK,U,TIME,N,NN,NSPRNG) use BHM ! 将指定module派生类型定义可用于程序 type(BHM_material) :: mat ! 定义相应的结构体名 type(BHM_spring) :: spr type(BHM_update) :: upd IMPLICIT real*8 (A-H,O-Z) real*8 dis_input ! 将每一步骤的U1读进来，自己定义便于辨识 include "../common/concom" ! 使用公共块，使用concom中的INC（增量步数） DIMENSION RATK(*),DATAK(*),U(*),TIME(*),N(*),F(*),NSPRNG(*) !!-------------------- 所需要的输入值,若弹簧参数不一样，需要另行加入（针对不同spring的编号） ---------------- if (INC=1) then mat%stiff_0=100.0 ! 需要输入的 初始刚度值(第一刚度)，赋给成员名数值 mat%ratio=0.3 ! 需要输入的 屈服刚度比，当等于0时，理想弹塑性表示 mat%dis_yie=0.8 ! 需要输入的 屈服位移值 mat%forc_yie=mat%stiff_0*mat%dis_yie ! 计算得到的屈服强度（力） spr%attr_spr(NSPRNG(1),1)=mat%stiff_0 ! 设置第 1 列为初始刚度（第一刚度） spr%attr_spr(NSPRNG(1),2)=mat%ratio ! 设置第 2 列为屈服刚度比 (第二刚度/第一刚度) spr%attr_spr(NSPRNG(1),3)=mat%dis_yie ! 设置第 3 列为 屈服位移值 , 需要 事先确定后 输入 spr%attr_spr(NSPRNG(1),4)=mat%forc_yie ! 设置第 4 列为 屈服强度（力） end if !!-------------------- 所需要的输入值,若弹簧参数不一样，需要另行加入（针对不同spring的编号） ---------------- !!----------------------------- 将当前步的相对位移读入,并与单元编号一一对应上 ---------------------- dis_input=U(1) ! 当前的相对位移赋给 dis_input upd%dis_cur(NSPRNG(1),INC)=dis_input ! 将 marc主程序的 相对位移 交给 subroutine 相应弹簧 当前位移 元素位置 !!----------------------------- 将当前步的相对位移读入,并与单元编号一一对应上 ---------------------- !!----------------------------- 计算每一步的 相对位移的 增量，便于下面通过滞回模型计算恢复力使用 ----------- if (INC<=1) then upd%dis_inc(NSPRNG(1),INC)=upd%dis_cur(NSPRNG(1),INC) ! 第一步的增量等于当前相对位移值 upd%dis_save(NSPRNG(1),INC)=upd%dis_cur(NSPRNG(1),INC) ! 将相对位移值保存起来,计算第二步的相对位移增量时，将使用 else if (INC>=2) then ! 大于等于第二步 upd%dis_inc(NSPRNG(1),INC)=upd%dis_cur(NSPRNG(1),INC) ! 当前步相对位移值增量(当前步值—上一步的值)， & -upd%dis_save(NSPRNG(1),INC-1) ! 便于写滞回规则时判断使用第一或第二刚度 upd%dis_save(NSPRNG(1),INC)=upd%dis_cur(NSPRNG(1),INC) ! 将相对位移值保存起来 end if !!----------------------------- 计算每一步的 相对位移的 增量，便于下面通过滞回模型计算恢复力使用 ----------- !!------------------------------ 如果处于弹性阶段范围 ，开始 ------------------ if (ABS(upd%dis_cur(NSPRNG(1),INC))<=spr%attr_spr(NSPRNG(1),3)) then ! 当前处于弹性范围 if (INC<=1) then ! 如果是第一步 upd%forc_cur(NSPRNG(1),INC)=spr%attr_spr(NSPRNG(1),1) & *upd%dis_cur(NSPRNG(1),INC) ! 直接计算,恢复力=第一刚度*当前步的相对位移 upd%forc_save(NSPRNG(1),INC)=upd%forc_cur(NSPRNG(1),INC) ! 保存起来 else if (INC>=2) then ! 大于等于第二步 upd%forc_cur(NSPRNG(1),INC)=upd%forc_save(NSPRNG(1),INC-1) & +spr%attr_spr(NSPRNG(1),1)*upd%dis_inc(NSPRNG(1),INC) ! 恢复力=前一步恢复力+第一刚度*【相对位移增量=(dis_cur-dis_save(INC-1))】,无论上一步是否处于二次刚度的非线性阶段 upd%forc_save(NSPRNG(1),INC)=upd%forc_cur(NSPRNG(1),INC) ! 保存起来 end if end if !!------------------------------ 如果处于弹性阶段范围 ，结束 ------------------ !!------------------------------ 如果进入二次阶段范围(第一象限),开始 ------------------ if (upd%dis_cur(NSPRNG(1),INC)>spr%attr_spr(NSPRNG(1),3)) then ! 大于屈服位移 if (INC<=1) then ! 如果是 第一步 upd%forc_cur(NSPRNG(1),INC)=spr%attr_spr(NSPRNG(1),4) & +spr%attr_spr(NSPRNG(1),1)*spr%attr_spr(NSPRNG(1),2) & *(upd%dis_cur(NSPRNG(1),INC)-spr%attr_spr(NSPRNG(1),3)) ! 直接计算，屈服力+二次刚度（当前位移-屈服位移） upd%forc_save(NSPRNG(1),INC)=upd%forc_cur(NSPRNG(1),INC) end if if (INC>=2) then ! 如果 大于等于第二步，需要分情况考虑 if (upd%dis_inc(NSPRNG(1),INC)<=0) then ! 当相对位移增量 小于0时 upd%forc_cur(NSPRNG(1),INC)=upd%forc_save(NSPRNG(1),INC-1) & +spr%attr_spr(NSPRNG(1),1)*upd%dis_inc(NSPRNG(1),INC) ! 判断相对位移的增量小于0，这时按照第一刚度卸载 upd%forc_save(NSPRNG(1),INC)=upd%forc_cur(NSPRNG(1),INC) end if if (upd%dis_inc(NSPRNG(1),INC)>0) then ! 当相对位移增量的大于0时，有两种情况：（1）小于交叉点位移值，按第一刚度加载（2）大于交叉点位移值，按第二刚度加载 upd%dis_cro(NSPRNG(1),INC)=(spr%attr_spr(NSPRNG(1),4) ! 计算按照第一刚度加载时，与第二刚度阶段的交叉点位移 & -spr%attr_spr(NSPRNG(1),1)*spr%attr_spr(NSPRNG(1),2)*spr%attr_spr(NSPRNG(1),3) & -upd%forc_save(NSPRNG(1),INC-1)+upd%dis_save(NSPRNG(1),INC-1)) & /(1.0-spr%attr_spr(NSPRNG(1),2))/spr%attr_spr(NSPRNG(1),1) if (upd%dis_cur(NSPRNG(1),INC)<=upd%dis_cro(NSPRNG(1),INC)) then ! 如果相对位移 小于 交叉点位移 upd%forc_cur(NSPRNG(1),INC)=upd%forc_save(NSPRNG(1),INC-1) & +spr%attr_spr(NSPRNG(1),1)*upd%dis_inc(NSPRNG(1),INC) ! 直接按照第一刚度计算 else if (upd%dis_cur(NSPRNG(1),INC)>upd%dis_cro(NSPRNG(1),INC)) then ! 否则，如果相对位移 大于 交叉点位移 upd%forc_cur(NSPRNG(1),INC)=upd%forc_save(NSPRNG(1),INC-1) ! 分阶段计算（交叉点前第一刚度，后为第二刚度） & +spr%attr_spr(NSPRNG(1),1)*(upd%dis_cro(NSPRNG(1),INC)-upd%dis_save(NSPRNG(1),INC-1)) & +spr%attr_spr(NSPRNG(1),2)*spr%attr_spr(NSPRNG(1),1) & *(upd%dis_cur(NSPRNG(1),INC)-upd%dis_cro(NSPRNG(1),INC)) end if upd%forc_save(NSPRNG(1),INC)=upd%forc_cur(NSPRNG(1),INC) ! 保存起来 end if end if end if !!------------------------------ 如果进入二次阶段范围(第一象限), 结束 ------------------ !!------------------------------ 如果进入二次阶段范围(第三象限),开始 ------------------ if (upd%dis_cur(NSPRNG(1),INC)<-spr%attr_spr(NSPRNG(1),3)) then ! 大于屈服位移(负值) if (INC<=1) then ! 如果是 第一步 upd%forc_cur(NSPRNG(1),INC)=-spr%attr_spr(NSPRNG(1),4) & -spr%attr_spr(NSPRNG(1),1)*spr%attr_spr(NSPRNG(1),2) & *(upd%dis_cur(NSPRNG(1),INC)-spr%attr_spr(NSPRNG(1),3)) ! 直接计算，屈服力+二次刚度（当前位移-屈服位移） upd%forc_save(NSPRNG(1),INC)=upd%forc_cur(NSPRNG(1),INC) end if if (INC>=2) then ! 如果 大于等于第二步，需要分情况考虑 if (upd%dis_inc(NSPRNG(1),INC)>=0) then ! 当相对位移增量 小于0时 upd%forc_cur(NSPRNG(1),INC)=upd%forc_save(NSPRNG(1),INC-1) & +spr%attr_spr(NSPRNG(1),1)*upd%dis_inc(NSPRNG(1),INC) ! 判断相对位移的增量小于0，这时按照第一刚度卸载 upd%forc_save(NSPRNG(1),INC)=upd%forc_cur(NSPRNG(1),INC) end if if (upd%dis_inc(NSPRNG(1),INC)<0) then ! 当相对位移增量的大于0时，有两种情况：（1）小于交叉点位移值，按第一刚度加载（2）大于交叉点位移值，按第二刚度加载 upd%dis_cro(NSPRNG(1),INC)=-(spr%attr_spr(NSPRNG(1),4) ! 计算按照第一刚度加载时，与第二刚度阶段的交叉点位移 & -spr%attr_spr(NSPRNG(1),1)*spr%attr_spr(NSPRNG(1),2)*spr%attr_spr(NSPRNG(1),3) & -upd%forc_save(NSPRNG(1),INC-1)+upd%dis_save(NSPRNG(1),INC-1)) & /(1.0-spr%attr_spr(NSPRNG(1),2))/spr%attr_spr(NSPRNG(1),1) if (upd%dis_cur(NSPRNG(1),INC)>=upd%dis_cro(NSPRNG(1),INC)) then ! 如果相对位移 小于 交叉点位移 upd%forc_cur(NSPRNG(1),INC)=upd%forc_save(NSPRNG(1),INC-1) & +spr%attr_spr(NSPRNG(1),1)*upd%dis_inc(NSPRNG(1),INC) ! 直接按照第一刚度计算 else if (upd%dis_cur(NSPRNG(1),INC)<upd%dis_cro(NSPRNG(1),INC)) then ! 否则，如果相对位移 大于 交叉点位移 upd%forc_cur(NSPRNG(1),INC)=upd%forc_save(NSPRNG(1),INC-1) ! 分阶段计算（交叉点前第一刚度，后为第二刚度） & +spr%attr_spr(NSPRNG(1),1)*(upd%dis_cro(NSPRNG(1),INC)-upd%dis_save(NSPRNG(1),INC-1)) & +spr%attr_spr(NSPRNG(1),2)*spr%attr_spr(NSPRNG(1),1) & *(upd%dis_cur(NSPRNG(1),INC)-upd%dis_cro(NSPRNG(1),INC)) end if upd%forc_save(NSPRNG(1),INC)=upd%forc_cur(NSPRNG(1),INC) ! 保存起来 end if end if end if !!------------------------------ 如果进入二次阶段范围(第三象限), 结束 ------------------ RATK(1)=1.0 F(1)=upd%forc_cur(NSPRNG(1),INC) ! 将 每一步的力 返回给F（1），这是required output RETURN END
|
||||||
