首页/文章/ 详情

abaqus 线弹性各向同性umat开发(附代码)

1月前浏览1143

概述

  帖子介绍了abaqus线弹性各向同性umat的开发过程,目的是科普umat二次开发最基础的工作流程,不涉及复杂高深的推导本构工作,仅仅包括umat子程序接口关键参数介绍、切线刚度矩阵和应力更新计算。并且设计一个立方体表面受动荷载算例,统计了位移、应力等数据,均与abaqus自带的材料模型保持一致。最后给出了子程序计算代码,拿过来就能直接用的那种

如何调用umat子程序

  调用umat子程序的步骤比较简单,在abaqus的cae界面中就可以实现,以下是具体步骤。  打开箭头1的Edit Material界面,在箭头2的General中点击3的User Material,然后下方会出现4的Mechanical Constants,我们需要在这个地方输入自己用到的参数,在这个帖子的线弹性各向同性材料中,第一个参数弹性模量    为    ,第二个参数为泊松比    为    
  上面的操作在inp文件中表现为

*Material, name=Material-1
**密度不会传进umat子程序,但是动力计算需要
*Density
2000.,
**User Material关键字的数据会传进umat中
*User Material, constants=2
 1e+10, 0.25

abaqus umat接口介绍

  umat子程序是abaqus平台提供的用于实现自己本构模型的接口,官方文档上可以找到该子程序的介绍,在介绍的最开始,官方文档先整了个活,给出了一个“Warning”

  Warning: The use of this subroutine generally requires considerable expertise. You are cautioned that the implementation of any realistic constitutive model requires extensive development and testing. Initial testing on a single-element model with prescribed traction loading is strongly recommended.

  我不得不说,这段话起到了开门见山、首尾呼应、吸引读者阅读兴趣、画龙点睛、blabla略略略等的作用,但事实上我第一次看这个子程序的时候并没有注意到这段话,后来才发现这是官方文档在给我们打预防针:小崽子们小心点儿!这个子程序很难!
  但其实大家完全不用担心,umat子程序本身并不难,总共就那么几个参数,但是数学难!难于上青天!推本构是一项十分有挑战性的工作,谁搞谁知道,劝退umat二次开发的不是umat接口本身,而是推公式。
  好了,umat接口的大致情况就是以上,下面是子程序的接口

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,
     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
     4 JSTEP(4)

      user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD
      and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT

      RETURN
      END

  为了让记忆更加的有逻辑,官方文档大致将上面的参数分为了两大类:需要用户更新的参数abaqus传进子程序的参数
  在线弹性各向同性材料中,需要更新的参数有应力项STRESS和切线刚度矩阵DDSDDE
  三维问题的应力项STRESS可表示为

 

  注意这里的应力排序,这是abaqus规定好的,一定要按着这个顺序写程序,不然会报错。我在用umat给用户自定义单元可视化中遇到过这个问题,当时我的应变-位移算子顺序跟abaqus规定的不一样,老是算不准,后来搞了很长时间才校正回来。
  三维问题的切线刚度矩阵DDSDDE表达方式不唯一,这里给出两种,后面的代码中分别将两种表达方式封装成了两个函数,可以直接调用。第一种表达方式为

 

  其中

 

  第二种表达方式为

 

  其中

 
 

  abaqus传进子程序的参数比较多,这里重点讲解线弹性各项同性材料中涉及的参数。
  参数STRAN(NTENS)为当前增量步开始时刻的应变项,括号中NTENS为应变项的个数,三维问题中为6,与应力项对应。
  参数DSTRAN(NTENS)为当前增量步结束时刻的应变增量,是abaqus估算的一个量,我们可以用它来计算当前增量步的应力增量。
  参数NDI是当前材料的正应力或正应变的个数,对应的是切应力或者切应变的个数NSHR,这两个参数的和为NTENS,即

 

  参数PROPS(NPROPS)是材料属性数组,是我们在cae界面中定义的参数,NPROPS是定义的参数个数,在这个帖子中,PROPS(1)是杨氏模量,PROPS(2)是泊松比。

算例

  模型尺寸为    ,杨氏模量为    ,泊松比为    ,采用C3D8单元离散,单元总数为125,节点总数为216,网格图为  模型一端设置固定边界,另一端施加竖直向面力,边界条件和荷载示意图为  荷载采用正弦荷载,计算总时长为10s,增量步长设置为0.1s,采用dynamic implicit分析步,施加的荷载幅值为1e5,幅值曲线为  统计了3.1s竖直向位移云图对比为(左边为abaqus自带材料模型计算结果,右边为umat计算结果  统计了3.1s的    应力云图对比为  统计了3.1s的    应力云图对比为  统计了3.1s的    应力云图对比为  统计了7.9s竖直向位移云图对比为  统计了7.9s的    应力云图对比为  统计了7.9s的    应力云图对比为  统计了7.9s的    应力云图对比为  统计了加载面节点加载向位移时程曲线对比为  统计了加载面节点    时程曲线对比为  统计了加载面节点    时程曲线对比为  统计了加载面节点    时程曲线对比为  可以看到,以上结果吻合得跟造数据一样好!其实就算我复 制一下ODB文件,然后再比较,你们也不知道,这可怎么是好!
  很简单,我把代码扔出来就好了

代码

  inp文件太长了,我就不放在帖子里面了,只把umat代码放出来,如果需要inp文件的可以私聊我,我发给你,但是验证出来数据有问题别吭气儿,人艰不拆,小声儿私信我,我给你解决。

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,
     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS ,DROT,PNEWDT,
     4 CELENT,DFGRDO,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
     
      INCLUDE 'ABA_PARAM.INC'
      
      CHARACTER*80CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD(3,3),DFGRD1(3,3),
     4 JSTEP(4)
     
     
      young=props(1)
      poiss=props(2)    

      DDSDDE=0.D0     
!这里调用的是切线刚度矩阵的第一种写法,第二种也可以
      call elas1(ddsdde,young,poiss)   
            
      STRESS=STRESS+MATMUL(DDSDDE,DSTRAN)
                       
      RETURN
      END
      
!切线刚度矩阵的第一种写法   

      subroutine elas1(ddsdde,young,poiss) 
      
      INCLUDE 'ABA_PARAM.INC'    
      
      double precision ddsdde(6,6)
      double precision young,poiss
      double precision ar   
      
      ar=young/((1+poiss)*(1-2*poiss))
      
      do i=1,3
       do ii=1,3
         DDSDDE(i,ii)=poiss*ar
       enddo
       DDSDDE(i,i)=(1-poiss)*ar       
      enddo
      
      do i=4,6
       DDSDDE(I,I)=((1-2.0D0*poiss)/2.0D0)*ar
      enddo   
      
      return
      end
      
!切线刚度矩阵的第二种写法            
      subroutine elas2(ddsdde,young,poiss) 
      
      INCLUDE 'ABA_PARAM.INC'
      
      DOUBLE PRECISION DDSDDE(6,6)
      DOUBLE PRECISION young,poiss
      DOUBLE PRECISION KK,GG  
        
      kk=young/(3.0D0*(1-2.0D0*poiss))
      gg=young/(2.0D0*(1+poiss))   
      
      do i=1,3
       do ii=1,3
         DDSDDE(i,ii)=kk-(2.0D0/3.0D0)*gg
       enddo
       DDSDDE(i,i)=kk+(4.0D0/3.0D0)*gg       
      enddo
      
      do i=4,6
       DDSDDE(i,i)=gg
      enddo         
      
      return
      end


来源:有限元先生

ACTMechanicalAbaqus二次开发UM材料科普
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-11-02
最近编辑:1月前
外太空土豆儿
硕士 我们穷极一生,究竟在追寻什么?
获赞 20粉丝 1文章 59课程 0
点赞
收藏
作者推荐

基于paraview的非结构化网格可视化

零概述ParaView是一个开源的多平台数据分析与可视化应用,广泛应用于科学计算领域中的大规模数据集可视化。下面是paraview中一些模型的可视化效果,包括液体搅拌、波相关和医学扫描数据。与此同时,ParaView对有限元结果的可视化也有强大的支持,特别是在处理大型数据集和复杂的网格时。有限元分析通常生成包含节点、单元和物理场变量的数据,这些数据可以通过ParaView进行后处理和可视化。在有限元结果可视化中,网格的可视化是重中之重,网格可大致分为结构化网格和非结构化网格,关于结构化网格的可视化,paraview的帮助文档有详细的说明,这里重点讲解非结构化网格的可视化,兼具标量数据和矢量数据的可视化。非结构化网格非结构化网格在处理如汽车、飞机的气动设计、复杂地形的环境模拟、生物医学工程中的器官建模等领域中非常有用,因为这些应用往往涉及到复杂的几何形状和边界条件。尽管非结构化网格的计算成本可能相对较高,但其灵活性和适应性使其成为解决复杂问题的重要工具。非结构化网格的主要特点包括:灵活性:可以很好地适应复杂的几何形状,不需要对计算域进行过多的简化或近似。节点和单元的多样性:网格中的单元可以是三角形、四边形(在二维情况下),或者四面体、六面体、金字塔和楔形等(在三维情况下)。这种多样性使得网格可以针对特定区域进行细化,以提高计算精度。局部细化:在需要更高精度的区域,可以局部增加网格密度,而不必在整个计算域中均匀细化网格,这样可以有效控制计算成本。生成和处理的复杂性:由于非结构化网格的灵活性,其生成和处理通常比结构化网格复杂。需要使用专门的算法来生成网格,并且在进行数值计算时,如有限元分析或流体流动模拟,需要处理不同形状和大小的单元。数据结构:非结构化网格的数据结构通常比结构化网格复杂,需要额外的数据结构来存储单元与节点之间的连接关系。paraview的vtk文件格式VTK(VisualizationToolkit)文件格式是ParaView中使用的一种重要数据格式,该文件记录了数据集的网格信息、节点坐标信息、标量、矢量和张量等信息,vtk文件有两种格式,分别为二进制和字符型(ASCII型),ASCIIVTK文件是一种文本格式的文件,用于存储和描述数据集,以便于可视化和分析。这种格式对于人类可读和编辑,但通常比二进制格式的文件更大,且读写速度较慢。ASCIIVTK文件的结构通常包括以下几个部分:文件头部:以#vtkDataFileVersion3.0开头,指明文件是VTK格式,以及版本号。版本号之后通常会跟随一个描述性的字符串,用于说明数据集的内容或来源。数据类型:接下来会指定数据集的类型,如ASCII或BINARY。尽管文件是ASCII格式,但这里指明的是文件中数据的存储方式,对于ASCII文件,这一行会是ASCII。数据集描述:紧接着是数据集的描述,可以包含任意文本,用于描述数据集的特征或提供其他信息。数据块(DataBlocks):数据部分开始于DATASET关键字,后面跟随数据集的类型,如STRUCTURED_POINTS、UNSTRUCTURED_GRID、POLYDATA等。之后是具体的数据定义,包括点(Points)和单元(Cells)的定义,以及可能的关联数据(如标量、向量等)。点(Points):定义数据集中的几何点,通常以一系列的三元组(x,y,z坐标)表示。单元(Cells):定义网格的拓扑结构,例如,对于非结构化网格,会列出构成每个单元的点的索引。数据字段(DataArrays):可以包含标量、向量或张量数据,这些数据与点或单元相关联。例如,对于一个标量数据,会列出每个点或单元的值。结束标记:数据部分以ENDDATA结束。下面是一个简单的ASCIIVTK文件结构示例#vtkDataFileVersion3.0ASCIIexampleASCIIDATASETUNSTRUCTURED_GRIDPOINTS3float...CELLS210...CELL_TYPES2...CELL_DATA2SCALARScell_scalarsfloat1LOOKUP_TABLEdefault...在这个例子中,定义了一个包含3个点的非结构化网格,由一个三角形单元组成。同时,还定义了一个名为cell_scalars的标量数据字段,包含两个值,分别对应于网格中的两个单元。ASCIIVTK文件的可读性使得它们在需要手动编辑或调试时非常有用,但因为文件较大且读写速度较慢,对于大规模数据集,二进制格式通常是更好的选择。非结构化网格示例下面采用手编的形式讲解非结构化网格在paraview中的显示,该例子包含两个多面体单元,如图所示这两个正方体的边长均为2,边中间的点均为边中点,方便手编坐标。这一些点的坐标为下面逐行讲解这两个非结构化单元的vtk文件(文末附有vtk全文)。首先是vtk文件头#vtkDataFileVersion2.0vtkfrommananulASCIIDATASETPOLYDATA第一行是涉及到paraview的版本。第二行是注释信息,随便写,我这是手动生成的,因此写了:vtkfrommananul,第三行的ASCII表示这是一个字符型的vtk文件,最后一行则代表着数据集的类型,这里表示的是多边形数据。然后是非结构化网格的节点坐标信息。POINTS19float200...points表示下面的数据集是点,19代表非结构化网格的节点总数,后面的float说明坐标的数据类型,后面的每一行依次是19个节点的xyz坐标。下面是非结构化网格的面信息。POLYGONS158240164512876498234049340321645678944518175561013184671110478121141011141341112151444171696131415161718516151289POLYGONS代表下面的数据是非结构化网格的多边形表面信息,一共15个表面,后面的82是所有的表面节点数目总和加上多边形个数,即82=15(多边形个数)+4(第1个多边形节点数目)+5(第2个多边形节点数目)+...+5(最后一个多边形节点数目)后面每一行数据的第一个数字是该多边形的节点数目,后面是具体的节点编号,注意节点标号的索引是从零开始的。为了方便理解,下面给出15个表面的示意图(手指头累断了),三个为一组,一共五组,大家可以跟上面的POLYGONS数据集对比着看。上面是非结构化网格的15个多边形表面,下面开始定义节点上的数据,可以是标量,如密度和泊松比等,可以是矢量,如位移、速度和加速度等,在正式开始给节点赋数据之前,要先声明节点数据,如POINT_DATA1919是我这个例子的节点数目,是可变的。下面首先给节点赋予标量数据,如SCALARSnamefloat1LOOKUP_TABLEmy_table1...第一行的name是数据量的名称,在paraview中会显示,float说明数据的类型是浮点数,后面的1表示我们这个例子中只有一组标量数据。后面的数据行要和节点数目对应。下面是两组矢量数据FIELDFieldData2U319double123...S319double123...这个例子设置了两组矢量数据,分别是:U和S,代表位移和应力。第一行的数据2,即表示有两组数据。第二行和第五行中的3表示该矢量有三个分量,19表示节点数目,double是矢量分量的数据类型。每组矢量后面的数据行行数都是节点数目。非结构化网格可视化效果上面的非结构化网格示例在paraview中显示为上面的示例中,设置了一组标量场,两组矢量场,对于矢量场,paraview会给出自己设计的分量,并且自动计算矢量对应的大小,如U数据,分别给出了xyz分量,和magnitude对于标量场,paraview直接给出数值,如v数据上面就是这个简单非结构化网格示例的可视化效果。非结构化网格vtk文件示例下面给出了上面手编的非结构化网格示例,里面的标量和矢量数据都是随意设置的,仅仅为了展示文件结构和各关键词的含义,搞清楚了这些,就可以手写程序。vtk文件为#vtkDataFileVersion2.0vtkfrommananulASCIIDATASETPOLYDATAPOINTS19float200220020000202212222122022002223123023224124024004204214POLYGONS158240164512876498234049340321645678944518175561013184671110478121141011141341112151444171696131415161718516151289POINT_DATA19FIELDFieldData2U319double12334556778991011111213131415151617171819192021212223232425252627272829293031313233333435353637373839S319double12334556778991011111213131415151617171819192021212223232425252627272829293031313233333435353637373839SCALARSvfloat1LOOKUP_TABLEmy_table12345678910111213141516171819paraview功能强大,上面只是非常简单的非结构化网格示例,仅仅涉及到paraview基本的功能,还有张量可视化等更加复杂的功能,大家可自行探索。点击卡片关注我们来源:有限元先生

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