首页/文章/ 详情

开源线弹性静力小变形UEL子程序分享

1月前浏览663
   

概述

       

分享一个适用于线弹性静力小变形的UEL子程序,程序来源于github,地址为(也可私信我获取全部程序)


https://github.com/bibekanandadatta/Abaqus-UEL-Elasticity
 

该UEL可计算的单元类型有


CPE4 C3D8 C3D20 CPE6 CPE8
 

该程序采用UVARM子程序给用户自定义单元进行可视化。在计算之前首先在自定义单元的初始位置嵌套一层abaqus自带的单元,并且将材料属性设置为无限小,因此嵌套的单元对模型的刚度无贡献,在UEL中计算单元应变和应力,并传递参数到UVARM子程序中。

程序结构讲解

     

采用Fortran完成所有的程序设计,主程序为


uel_mech
 

其余的功能采用fortran的module进行封装,供主程序调用,该程序编写的module有









global_parameters.forerror_logging.forlinear_algebra.forlagrange_element.forgauss_quadrature.forsolid_mechanics.forpost_processing.formode.for
 

在主程序中,需要采用include引用上述module,具体的引用代码为









include 'global_parameters.for'     ! global parameters moduleinclude 'error_logging.for'         ! error/ debugging moduleinclude 'linear_algebra.for'        ! linear algebra moduleinclude 'lagrange_element.for'      ! Lagrange element moduleinclude 'gauss_quadrature.for'      ! Guassian quadrature moduleinclude 'solid_mechanics.for'       ! solid mechanics moduleinclude 'post_processing.for'       ! post-processing moduleinclude 'mode.for'
 

程序内部涉及到线性代数的矩阵运算,因此Fortran需要调用mkl计算库,程序提供了env文件,用于程序运行时调用mkl计算库,env文件中的内容为


compile_fortran += ['/Qmkl:sequential']
 

主程序中的uel程序用于向abaqus主程序返回刚度矩阵和右端量等,具体程序为














































































































































! **********************************************************************! ****************** ABAQUS USER ELEMENT SUBROUTINE ********************! **********************************************************************! make sure to have the correct directory      include 'global_parameters.for'     ! global parameters module      include 'error_logging.for'         ! error/ debugging module      include 'linear_algebra.for'        ! linear algebra module      include 'lagrange_element.for'      ! Lagrange element module      include 'gauss_quadrature.for'      ! Guassian quadrature module      include 'solid_mechanics.for'       ! solid mechanics module      include 'post_processing.for'       ! post-processing module      include 'mode.for'      SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,     & PROPS,NPROPS,COORDS,MCRD,NNODE,Uall,DUall,Vel,Accn,JTYPE,TIME,     & DTIME,KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,     & NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROPS,PERIOD)      ! This subroutine is called by Abaqus with above arguments      ! for each user elements defined in an Abaqus model. Users are      ! responsible for programming the element tangent/ stiffness      ! matrix and residual vectors which will be then assembled and      ! solved by Abaqus after applying the boundary conditions.      use global_parameters      use error_logging      use user_element      use post_processing      INCLUDE 'ABA_PARAM.INC'      DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),PROPS(*),     & SVARS(*),ENERGY(8),COORDS(MCRD,NNODE),UAll(NDOFEL),     & DUAll(MLVARX,*),Vel(NDOFEL),Accn(NDOFEL),TIME(2),PARAMS(*),     & JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),     & PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)      ! user coding to define RHS, AMATRX, SVARS, ENERGY, and PNEWDT      integer, intent(in)   :: NDOFEL, NRHS, NSVARS, NPROPS, MCRD      integer, intent(in)   :: NNODE, JTYPE, KSTEP, KINC, JELEM      integer, intent(in)   :: NDLOAD, JDLTYP, NPREDF, LFLAGS      integer, intent(in)   :: MLVARX, MDLOAD, JPROPS, NJPROPS      real(wp), intent(in)  :: PROPS, COORDS, DUall, Uall, Vel, Accn      real(wp), intent(in)  :: TIME, DTIME, PARAMS, ADLMAG, PREDEF      real(wp), intent(in)  :: DDLMAG, PERIOD      real(wp), intent(out)           :: RHS, AMATRX      real(wp), intent(out), optional :: SVARS, ENERGY, PNEWDT      character(len=8)      :: abqProcedure      integer               :: nDim, nStress      character(len=2)      :: analysis      logical               :: nlgeom      integer               :: nInt, nPostVars      integer               :: lenJobName,lenOutDir      character(len=256)    :: outDir      character(len=256)    :: jobName      character(len=512)    :: errFile, dbgFile      type(logger)          :: msg      ! initialize the output variables to be defined for Abaqus subroutine      energy        = zero      amatrx        = zero      rhs(:,nrhs)   = zero      ! Open an error and debug file for the current job.      ! See Abaqus documentation for Fortran unit number.      ! File unit numbers are defined in the error_logging module.      call getJobName(jobName,lenJobName)      call getOutDir(outDir,lenOutDir)      errFile = trim(outDir)//'\aaERR_'//trim(jobName)//'.dat'      dbgFile = trim(outDir)//'\aaDBG_'//trim(jobName)//'.dat'      call msg%fopen( errfile=errFile, dbgfile=dbgFile )      ! change the LFLAGS criteria as needed (check Abaqus UEL manual)      if( (lflags(1)  .eq.  1) .or. (lflags(1)  .eq.  2) ) then        abqProcedure = 'STATIC'      else        call msg%ferror(flag=error, src='UEL',     &       msg='Incorrect Abaqus procedure.', ia=lflags(1))        call xit      end if      ! check if the procedure is linear or nonlinear      if (lflags(2)  .eq.  0) then        nlgeom = .false.      else if (lflags(2) .eq. 1) then        nlgeom = .true.      end if      ! check to see if it's a general step or a linear purturbation step      if(lflags(4)  .eq.  1) then        call msg%ferror(flag=error, src='UEL',     &       msg='The step should be a GENERAL step.', ia=lflags(4))        call xit      end if      ! set parameters specific to analysis and element types      if ((jtype .ge. 1).and.(jtype .le. 4)) then         ! three-dimensional analysis        nDim      = 3        analysis  = '3D'        nStress   = 6      else if ((jtype .ge. 5).and.(jtype .le. 8)) then    ! plane strain analysis        nDim      = 2        analysis  = 'PE'        nStress   = 3      else if ((jtype .ge. 9).and.(jtype .le. 12)) then   ! plane stress analysis        nDim      = 2        analysis  = 'PS'        nStress   = 3      else if ((jtype .ge. 13).and.(jtype .le. 16)) then  ! axisymmetric analysis        ! CAUTION: axisymmetric elements are not validated yet        nDim      = 2        analysis  = 'AX'        nStress   = 4      else        call msg%ferror( flag=error, src='UEL',     &       msg='Element is unavailable.', ia=jtype )        call xit      end if      nInt      = jprops(1)      nPostVars = jprops(2)      ! array containing variables for post-processing      if (.not. allocated(globalPostVars)) then        allocate( globalPostVars(numElem,nInt,nPostVars) )         ! print job-related information the first time        call msg%finfo('---------------------------------------')        call msg%finfo('------- Abaqus SMALL STRAIN UEL -------')        call msg%finfo('--------- SOURCE: uel_mech.for --------')        call msg%finfo('---------------------------------------')        call msg%finfo('--- Abaqus Job: ', ch=trim(jobName))        call msg%finfo('---------------------------------------')        call msg%finfo('------- PROCEDURE       = ', ch=abqProcedure)        call msg%finfo('------- ANALYSIS TYPE   = ', ch=analysis)        call msg%finfo('---------- NLGEOM       = ', la=nlgeom)        call msg%finfo('------- MODEL DIMENSION = ', ia=nDim)        call msg%finfo('------- ELEMENT NODES   = ', ia=nNode)        call msg%finfo('---------------------------------------')        call msg%finfo('-------- INTEGRATION SCHEME -----------')        call msg%finfo('----------- NINT   = ', ia=nInt)        call msg%finfo('---------------------------------------')        call msg%finfo('---------- POST-PROCESSING ------------')        call msg%finfo('--- NO OF ELEMENTS            = ',ia=numElem)        call msg%finfo('--- OVERLAY ELEMENT OFFSET    = ',ia=elemOffset)        call msg%finfo('--- NO OF VARIABLES AT INT PT = ',ia=nPostVars)        call msg%finfo('---------------------------------------')      end if       ! call the element subroutine with extended input arguments       call uelMech(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,     & PROPS,NPROPS,COORDS,MCRD,NNODE,Uall,DUall,Vel,Accn,JTYPE,TIME,     & DTIME,KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,     & NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROPS,PERIOD,     & NDIM,ANALYSIS,NSTRESS,NINT)
     END SUBROUTINE UEL
 

主程序中UVARM用于给用户自定义单元可视化,具体代码为









































! **********************************************************************! ************** ABAQUS USER OUTPUT VARIABLES SUBROUTINE ***************! **********************************************************************
     SUBROUTINE UVARM(UVAR,DIRECT,T,TIME,DTIME,CMNAME,ORNAME,     & NUVARM,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,NDI,NSHR,COORD,     & JMAC,JMATYP,MATLAYO,LACCFLA)
     ! This subroutine is called by Abaqus at each material point (int pt)      ! to obtain the user defined output variables for standard Abaqus      ! elements. We used an additional layer of standard Abaqus elements      ! with same topology (same number of nodes and int pts) on top of      ! the user elements with an offset in the numbering between the user      ! elements and standard elements. This number is defined in the      ! post_processing module and should match with Abaqus input file.
     use global_parameters      use post_processing
     INCLUDE 'ABA_PARAM.INC'
     CHARACTER*80 CMNAME,ORNAME      CHARACTER*3 FLGRAY(15)      DIMENSION UVAR(NUVARM),DIRECT(3,3),T(3,3),TIME(2)      DIMENSION ARRAY(15),JARRAY(15),JMAC(*),JMATYP(*),COORD(*)
     ! the dimensions of the variables FLGRAY, ARRAY and JARRAY      ! must be set equal to or greater than 15.
     ! explicityly define the type for uvar to avoid issues      real(wp)        :: uvar
     ! assign the stored global variables to the UVAR for Abaqus to process      uvar(1:nuvarm)  = globalPostVars(noel-elemOffset,npt,1:nuvarm)
     END SUBROUTINE UVARM
! **********************************************************************! **********************************************************************
 

为验证程序,提供了一系列的inp文件








cantilever_U7_CPE4_elem_80_elastic.inpsingle_elem_uniaxial_U3_C3D8_elastic.inpsingle_elem_uniaxial_U4_C3D20_elastic.inpsingle_elem_uniaxial_U5_CPE3_elastic.inpsingle_elem_uniaxial_U6_CPE6_elastic.inpsingle_elem_uniaxial_U7_CPE4_elastic.inpsingle_elem_uniaxial_U8_CPE8_elastic.inp
 

这里以最后一个inp文件为例,简要讲解inp文件中的修改




















































































*Heading** abaqus/ standard input file for uniaxial test of** single CPE8 (U8) type element for uel_mech.for**节点坐标信息*Node1,           0.,           0....****************************************************用户自定义单元信息*User Element,Type1,2*Element, type=U81, 1, 2, 4, 3, 5, 6, 7, 8***************************************************Nset, nset=All, generate 1,  8,  1*Elset, elset=All 1,*Nset, nset=Bottom 1, 2, 5*Elset, elset=Bottom 1,*Nset, nset=Left 1, 3, 8*Elset, elset=Left 1,*Nset, nset=Right, generate 2,  6,  2*Elset, elset=Right 1,*Nset, nset=Top 3, 4, 7*Elset, elset=Top 1,****************************************************给用户自定义单元赋予属性*uel property,elset=all** E, nu, nInt, nPostVars79E9, 0.3, 9, 6**************************************************** same element connectivity as the user element** uses built-in abaqus element to overlay results** negligible young's modulus won't affect the result**嵌套的一层abaqus自带的单元用于可视化*Element, type=CPE8100001, 1, 2, 4, 3, 5, 6, 7, 8*Elset, elset=elDummy 100001*Solid section, elset=elDummy, material=Dummy**给嵌套的一层单元赋予材料属性,无限小的材料属性*Material, name=Dummy*Elastic1.e-20** strain vector (3) and stress vector (3) as user output**UVARM子程序用于可视化*User output variables6****************************************************给左边界和底边界固定住*BoundaryLeft, 1, 1Bottom, 2, 2****************************************************定义幅值曲线*Amplitude, name=LinearRamp0., 0., 1., 1.***************************************************Step, name=Stretch, nlgeom=NO, inc=1000**静力分析步*Static1, 1, 0.001, 1*给右边界施加力*Boundary, amplitude=LinearRampRight, 1, 1, 0.01*Restart, write, frequency=0 *Output, field, frequency=1, time marks=no*node output, nset=Allu,rf** specify the user-defined output for the dummy element set*element output, elset=elDummyuvarm*End Step**************************************************
 


算例

       

下面以cantilever_U7_CPE4_elem_80_elastic为例运行子程序,直接采用bat文件运行文件,bat文件为



call abaqus job=cantilever_U7_CPE4_elem_80_elastic user=uel_mech cpus=10 intpause
 

程序运行结束后,cmd窗口中会输出





















---------------------------------------------- Abaqus SMALL STRAIN UEL ---------------- SOURCE: uel_mech.for -------------------------------------------------- Abaqus Job: single_elem_uniaxial_U7_CPE4_elastic---------------------------------------------- PROCEDURE       = STATIC------- ANALYSIS TYPE   = PE---------- NLGEOM       = F------- MODEL DIMENSION = 2------- ELEMENT NODES   = 4----------------------------------------------- INTEGRATION SCHEME ---------------------- NINT   = 4------------------------------------------------- POST-PROCESSING --------------- NO OF ELEMENTS            = 50000--- OVERLAY ELEMENT OFFSET    = 100000--- NO OF VARIABLES AT INT PT = 6---------------------------------------
 

这是uel子程序中编写的输出语句,用于监测程序的运行过程

最终未变形的odb文件在cae中显示为

   

可以发现,在正常的节点位置出现了“XXX”形状的图形,这些是用户自定义单元,如果不采用可视化技术,界面中就只有这些“XX”形状的图形。

位移数值(magnitude)云图为

   

UVARM1~UVARM6的数值云图为

   
   
   
   
   
   

总结

       

该程序采用线弹性小变形理论编写了适用于CPE8等单元的UEL程序,采用module功能将不同的代码段封装到不同的文件中,在主程序uel_mech中进行引用。为了给用户自定义单元进行可视化,采用UVARM子程序进行应变应力等数据的传递,通过在自定义单元上嵌套一层属性值无限小的abaqus自带单元完成了自定义单元的可视化。

      

来源:有限元先生
ACTAbaqusUGUM理论材料
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-10-14
最近编辑:1月前
外太空土豆儿
硕士 我们穷极一生,究竟在追寻什么?
获赞 12粉丝 1文章 52课程 0
点赞
收藏
作者推荐

蒙皮单元在UEL二次开发中的作用

零概述为了用户能够在abaqus中实现自己的有限元理论,abaqus提供了二次开发接口uel,uel子程序设计用来完成一个单元的核心计算,包括单元对运动方程左端贡献:刚度矩阵、质量矩阵和阻尼矩阵,以及对方程右端不平衡力的贡献,这里要注意区分“方程右端量”和“右端不平衡力”的区别。uel是一款功能强大的子程序,但同时也是一把双刃剑。用户自定义单元无法直接与abaqus自带的功能完美结合,比如给用户自定义单元施加复杂荷载,面力、体力等,也无法给用户自定义单元定义surface,也就无法与abaqus自带的接触非线性模型相结合。帖子分享了蒙皮单元在uel二次开发中几个常见的作用,包括可视化、施加复杂荷载和定义接触属性,重在讲解思路。壹原理首先分析用户自定义单元为什么无法与abaqus某些功能直接结合。从荷载讲起。在有限元计算中,任何荷载都是需要被转化为节点荷载的,无论是线荷载、面荷载还是体荷载,形函数在荷载转化过程中起到了至关重要的作用,这些节点力最终组装到运动方程的右端,参与迭代计算。但是abaqus主程序是无法得知用户自定义单元的形函数的,因为uel只向主程序输出一些关键矩阵,这些关键矩阵直接被组装到总体模型矩阵中,这么一来,用户自定义单元的形函数只有写程序的人知道,那么复杂荷载也就无法转化到节点上了,因此,用户自定义单元只能定义基于节点的荷载,比如节点位移、节点力和节点速度等等。与此同时,用户自定义单元的可视化也是无法进行的,因为位移、应变等数据的磨平也需要形函数,这就是用户自定义单元的odb文件中只有“X”形状的原因。如下图通过以上分析,主要的问题是用户自定义单元无法提供方程的右端量,也就是节点数据。完成abaqus的所有功能,包括但不限于定义surface、体荷载和接触属性等。那么,只需要在自定义单元原有的节点上面在嵌套一层abaqus内置的单元,被称为蒙皮单元,蒙皮单元的材料属性被设置的无限小,因此蒙皮单元对运动方程的左端量没有贡献,只对右端的荷载向量有贡献,与此同时,也能完成自定义单元的可视化。贰用户自定义单元的可视化设计悬臂梁顶端受节点力算例。悬臂梁尺寸为2×2×8,弹性模量1e10,密度2000,泊松比0.25,悬臂端受垂直轴向的节点力,幅值为1e6,另一端为固定端,边界条件和荷载示意图为网格尺寸为0.5,共计单元数目256,节点数目为425,网格为采用静力通用分析步,计算总时长为1,最小增量步长1e-5。(静力计算中的计算时长没有物理意义,只有数值意义)采用自编的C3D8单元uel与abaqus做对比,设置三种对比工况,分别为:abaqus计算,不带蒙皮单元的uel计算,带蒙皮单元的uel计算,如下图可以看到,三者计算结果一致。且不带蒙皮单元的云图显示效果很差,仅仅显示了“XX”形状,非常不利于结果展示。第三幅图为带有蒙皮单元的uel计算效果,云图同时显示了“XX”和abaqus自带的单元,把自定义单元隐藏之后的效果为隐藏了自定义单元之后,显示的云图与abaqus保持一致。同样地,当隐藏了自定义单元以后,u1方向的位移云图为u2方向的位移云图为初步说明了蒙皮单元在可视化方面的应用是可行的。叁施加面力设计悬臂梁一端受垂直于轴向面力算例。计算模型参数与分析步信息与上述相同。边界条件和荷载示意图为这里为了显示面力,将模型视图调整为透明,可以看到施加的是垂直于轴向的面力。对该计算模型设置两个对比工况,分别为abaqus计算和带有蒙皮单元的uel计算。因为不带有蒙皮单元无法施加面力,所以没有设置不带蒙皮单元的工况。统计了加载向的位移云图,不隐藏自定义单元的时候,显示如下隐藏了自定义单元之后的云图效果为u1方向隐藏了自定义单元之后的云图u2方向隐藏了自定义单元之后的云图以上云图和计算效果说明了蒙皮单元施加面力荷载是可行的。肆施加体力设计悬臂梁一端受轴向体力算例。计算模型参数与分析步信息与上述相同。边界条件和荷载示意图为这里为了显示面力,将模型视图调整为透明,可以看到施加的是垂直于轴向的面力。对该计算模型设置两个对比工况,分别为abaqus计算和带有蒙皮单元的uel计算。因为不带有蒙皮单元无法施加体力,所以没有设置不带蒙皮单元的工况。统计了加载向的位移云图,不隐藏自自定义单元的时候,显示如下隐藏了自定义单元之后的云图效果为u1方向隐藏了自定义单元之后的云图u2方向隐藏了自定义单元之后的云图以上云图和计算效果说明了蒙皮单元施加体力荷载是可行的。伍用于接触非线性问题设计了滑块在斜面滑落的算例,滑块一共有两个运动过程:在斜面滑动和出斜面运动。该算例采用比例边界有限元法计算,涉及到给用户自定义单元定义重力荷载和接触属性,滑块与斜面的几何尺寸如下图:滑块与斜面的网格图如下:该算例中的三种方法均采用动力隐式算法,计算总时长为10s,设置固定增量步长为0.01s,增量步数共计1000。重力加速度为9.81m/s²。具体参数见附件。滑块两个运动过程受力分析见下图:上面左图为斜面内的受力分析,右图为斜面外的受力分析。可知滑块沿斜面的分力与摩擦力大小比值有两种情况:当滑块在斜面内运动的时,遵循如下控制方程:通过上式子,可求得位移和速度的解析解。滑块运动水平向位移时程对比如图误差分析:观察四个时程曲线,发现滑块在斜面上运动的过程中,四种计算结果均吻合良好。而在滑块离开斜面的第二运动过程中,三个数值解吻合较好,但均与解析解有一定的差别。这是因为解析求解是建立在质点的假设上,但是数值解是一个小方块,因此方块第一个角点离开斜面的时候,方块并未完全离开斜面,所以误差会从方块的第一个角点离开斜面开始。在该算例中,完成荷载和接触定义的同时,附加的蒙皮单元兼具用户自定义单元的位移、速度和加速度可视化功能。当蒙皮单元和用户自定义单元共存的时候,位移云图为当只显示用户自定义单元的时候,ABAQUS显示如图当只显示蒙皮单元的时候,ABAQUS显示如图这个例子结合了蒙皮单元的所有作用:可视化、施加荷载和接触属性定义,验证了蒙皮单元在uel二次开发中的作用陆蒙皮单元的限制当然,蒙皮单元也是有限制的。1、蒙皮单元只适用于用户自定义单元和abaqus内置单元单元节点形状相同的时候。为了避免这种限制,可以将几何拓扑信息复杂的自定义单元拆分为最简单的四面体单元即可。2、蒙皮单元只能给自定义单元的位移、速度和加速度进行可视化,不能对应变、应力等数据进行可视化。应力应变等数据的可视化可以采用python处理到第三方软件进行可视化,如sufur、paraview和tecplot等软件。也可以采用UMAT子程序进行可视化,但是采用UMAT同样要求自定义单元的形状与abauqs内置单元相同。注:除用户子程序外,上述算例所有inp文件均来源,如有需要,可留言或私信。点击卡片关注我们来源:有限元先生

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈