当前位置:首页 > 程序开发 / Fortran C# VC MATLAB >

    最新论文下载:



















 

基于Fortran的ABAQUS子程序TJU.Plastic-Explicit

 

 

下载次数:

          

Fortran语言的ABAQUS子程序TJU.plastic-explicit

    (Fortran 语言 / 配合TJU.SAP2ABAQUS接口导出的inp文件)         
uploaded to  www.liuguohuan.net by   柳   
天津大学 建筑工程学院 ; 天津大学 前沿技术研究院; 国家重点实验室-水利工程仿真与安全 

更多程序下载
下载链接:

VUMAT详细内容链接http://www.liuguohuan.net/lunwenxz/2013/1212/150.html


!! --------------------------    定义公共----------------------------
        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

 

 

更多程序下载