ABAQUS用户子程序学习小结

ABAQUS用户子程序学习小结 1 FORTRAN语言中的“I-N规则”:I、J、K、L、M、N开头的为整型变量，其他开头为实型变量; 1.1 Fortran书写格式:1-5列:标号区;6列:续写标号区，一般就写"1";7-72列:语句区，本区内空格无效;注释行以C开头，本行内书写格式无要求;参考周煦《FORTRAN77结构化程序设计》，中国科学技术出版社，1995，38-40页 2 DIMENSION COORDS(3)表示声明一个含3个元素的数组，下标分别为1、2、3，访问形式为COORDS(n)，n为1,3; 3 子程序(*.for)文件中如何输出调试信息: WRITE(6,*)'COORDS(1)',COORDS(1)，在*.dat文件中可看到输出，如果希望WRITE输出到msg文件中，则写为WRITE(7,*)'COORDS...; 4 用户子程序DLOAD中COORDS数组的含义:COORDS(1)也是一个数组，存贮单元集合中所有单元积分点的X坐标，COORDS(2)存贮Y坐标，相应INP文件中的写法为: *DLOAD PY,PYNU 其中PY为单元集合名称，定义方法为: *Elset, elset=BEAM, generate 1, 5, 1 ... *ELSET,ELSET=PY BEAM 5 DLOAD中F的定义方法: F只有定义在单元积分点上才有效，例如: F=1.0*COORDS (1) 附一个简单实例: beam.inp文件: *Heading ** Job name: Job-1 Model name: beam *Preprint, echo=NO, model=NO, history=NO, contact=NO ** ** PARTS ** *Part, name=PART-1 *End Part ** ** ASSEMBLY ** *Assembly, name=Assembly ** *Instance, name=PART-1-1, part=PART-1 *Node 1, 0., 0. 2, 20., 0. 3, 40., 0. 4, 60., 0. 5, 80., 0. 6, 100., 0. *Element, type=B31 1, 1, 2 2, 2, 3 3, 3, 4 4, 4, 5 5, 5, 6 *Elset, elset=BEAM, generate 1, 5, 1 ** Region: (Section-1-BEAM:BEAM), (Beam Orientation:BEAM) ** Section: Section-1-BEAM Profile: Profile-1 *Beam Section, elset=BEAM, material=STEEL, temperature=GRADIENTS, section=RECT 0.2, 5. 0.,0.,-1. *End Instance *Nset, nset=ENDS, instance=PART-1-1 1, 6 *Nset, nset=_M4, internal, instance=PART-1-1 6, *Nset, nset=_M5, internal, instance=PART-1-1 1, *End Assembly ** ** MATERIALS ** *Material, name=STEEL *Elastic 210000., 0.3 *ELSET,ELSET=PY BEAM ** ** BOUNDARY CONDITIONS ** ** Name: Disp-BC-1 Type: Symmetry/Antisymmetry/Encastre *Boundary _M4, ENCASTRE ** ---------------------------------------------------------------- ** ** STEP: Step-1 ** *Step, name=Step-1 *Static ** ** LOADS ** ** Name: CFORCE-1 Type: Concentrated force *DLOAD PY,PYNU ** ** OUTPUT REQUESTS ** ** ** FIELD OUTPUT: F-Output-1 ** *Output, field, variable=PRESELECT ** ** FIELD OUTPUT: F-Output-2 ** *Output, field *Element Output SF, ** ** HISTORY OUTPUT: H-Output-1 ** *Output, history *Node Output, nset=ENDS CF1, CF2, CF3, CM1, CM2, CM3, RF1, RF2 RF3, RM1, RM2, RM3, U1, U2, U3, UR1 UR2, UR3 *El Print, freq=999999 *Node Print, freq=999999 *End Step bbb.for文件 SUBROUTINE DLOAD(F,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,COORDS, 1 JLTYP,SNAME) C INCLUDE 'ABA_PARAM.INC' C DIMENSION TIME(2), COORDS (3) CHARACTER*80 SNAME WRITE(6,*)'COORDS(3)',COORDS(3) F=1.0*COORDS (1) RETURN END 运行方法: 在Abaqus Command提示符后输入: abaqus job=beam user=bbb interactive 子程序如下～如何编写，调研子程序, 子程序USDFLD中给固化度赋了一个初值1×10-4。 SUBROUTINE USDFLD(FIELD,STATEV,PNEWDT,DIRECT,T,CEIENT, 1 TIME,DTIME,OMNAME,ORNAME,NEIELD,NSTATV,NOEL,NPT,LAYER, 2 KSPT,KSTEP,KINC,NDI,NSHR,COORD,JMAC,JMATYP,MATLAYO,LACCFLA) C INCLUDE'ABA_PARAM.INC' C CHARACTER*80 CMNAME,ORNAME CHARACTER*3 FLGRAY(15) DIMENSION FIELD(NFIELD),STATEV(NSTATV),DIRECT(3,3), 1 T(3,3),TIME(2) DIMENSION ARRAY(15),JARRAY(15),JMAC(*),JMATYP(*),COORD(*) C IF(KINC.EO,1) THEN STATEV(1)=1E-4 ELSE END IF FIELD(1)=STATEV(1) C RETURN END (2)DISP 在本例中用户子程序DISP定义温度边界条件。 SUBROUTINE DISP(U,KSTEP,KINC,TIME,NODE,NOEL,JDOF,COORDS) C INCLUDE’ABA_PARAM.INC’ C DIMENSION U(3),TIME(2),COORDS(3) C IF(TIME(2),LE,1800.) THEN U(1)=303.+TIME(2)/30. ELSE IF(TIME(2),LE,3600.) THEN U(1)=363, ELSE IF(TIME(2),LE,4650,) THEN U(1)=363.+(TIME(2)-3600.)/30. ELSE U(1)=398. END IF C RETURN END (3)FILM 在本例中用户子程序FILM定义对流热边界条件，H(1)给定表面对流换热系数，SINK给定与图16-21的工艺曲线一致的环境温度。 SUBROUTINE FILM(H,SINK,TIMP,JSTEP,JINC,TIME,NOEL,NPT,COORDS, 1 JLTYP,FIELD,NFIELD,SNAME,JUSERNODE,AREA) C INCLUDE’ABA_PARAM,INC’ C DIMENSION COORDS(3),TIME(2),FIELD(NFIELD),H(2) CHARACTRR*80 SAME C H(1)=15 IF(TIME(2).IF.1800.) THEN SINK=303.+TIME(2)/30. ELSE IF(TIME(2).IF.3600.) THEN SINK=363. ELSE IF(TIME(2).IF.4650.) THEN SINK=363.+(TIME(2)-3600.)/30 ELSE SINK=398. END IF C RETURN END (4)HETVAL 本例中用户子程序HETVAL定义复合材料内部的反应热。STATEV(1)即为从用户子程序USDFLD中传递过来的固化度;STATEV(2)为固化率(固化度的变化率)，由与固化度、温度、升温速率等变量相关的固化动力学方程确定，决定了反应热的放热率;FLUX(1)表示反应热产生的单位时间、单位体积的热量。本例中，用Euler法对固化动力学方程进行数值 积分，得到的固化度再通过STATEV(1)传递回用户子程序USDFLD，读者可以采用其他的数值积分方法。由于有些内容涉及相关部门的内部技术，所以未给出函数的具体形式。 SUBROUTINE HETVAL(CMNAME,TEMP,TIME,DTIME,STATEV,FLUX, 1 FREDEF,DPRED) C INCLUDE ABA_PARAM,INC’ C CHARACTER*80 CMNAME DIMENSION TEMP(2),STATEV(3),PREDEF(1),TIME(2),FLUX(2), 1 DPRED(1) C IF(TEMP(1).LT.363.15) THEN STATEV(2)=0.0 ELSE STATEV(2)=<此处函数有固化度、温度、升温速率等变量相关的固化动力学方程确定> STATEV(1)=STATEV(1)+STATEV(2)*DTIME END IF FLUX(1)=242971740*STATEV(2) C RETURN END 我的多个子程序就是简单拼装到一个fortran文件的，一个子程序完毕，紧跟另一个子程序，实践证明:简单拼装是可行的，你之所以出错，我怀疑最大可能是源程序本身就有问题,还有，你的两个子程序都是DFLUX，如果像我这样第一个UFLUX后再接另一个UFLUX，我觉得不可行，因为两个子程序同名，而简单拼装是将不同子程序放到一个fortran文件里，建议:要么只使用一个UFLUX，但根据不同情况调用不同热源，要么使用三个UFLUX，第一个用于调用另两个UFLUX用，可以参阅使用UMAT调用不同的材料本构 第一个UMAT里只有如下语句: IF (CMNAME(1:4) .EQ. 'MAT1') THEN CALL UMAT_MAT1(argument_list) ELSE IF(CMNAME(1:4) .EQ. 'MAT2') THEN CALL UMAT_MAT2(argument_list) END IF 如材料名的前四字为MAT1，则调用子程序UMAT_MAT1;否则，如为MAT2，则调用子程序UMAT_MAT2 第二、三个UMAT的子程序为UMAT_MAT1，UMAT_MAT2，并要分别对CMNAME赋值以上为个人意见，仅供参考，请谨慎使用 【讨论】初学子程序:coords(3)的问题～ 在standard中，许多子程序中的参数都设计到从inp中提供的参数，比方说time参数是分析步的total time和current step time，这些还能看懂，不太清楚的就是coords这个参 数，它在help中指的是An array containing the coordinates of this point. These are the current coordinates if geometric nonlinearity is accounted for during the step (see “Procedures: overview,” Section 6.1.1); otherwise, the array contains the original coordinates of the point. 其中的this point是指的哪个点, 比方说，我在dflux加载时，可以这样定义不均匀分布热载 *dflux， set,s3nu， 这样由于nu的存在，它需要调用我提供的子程序，那么它传递的coords是指的什么点的坐标呢,set是指的一组受热载的单元。 搞不太清楚，子程序是新手，哪位大侠做过，帮忙解释一下～ 關於this point指的就是你定義的節點， 每個節點都有座標定義，若是三維的幾何模型，那麼節點就會有x,y,z三個座標， 所以相對應於副程式中的Coords Array，那麼它就是1x3的Array。 至於它其中說的geometry nonlinear需注意的事項， 是說一般若我們不考慮幾何非線性，那麼整個分析是屬於小變形，因為小變形所以忽略其位移量， 所以，在整個分析中導入副程式中的節點座標為原始定義座標，自始至終都是一樣。 但是，若整個分析為幾何非線性，那麼它會是大變形，所以其節點位移量就不可忽略。 如此一來，在分析時導入副程式中的節點座標，會隨分析過程中的變化將位移量考慮進去，不斷更新。 對於你的問題，它傳遞的就是set這單元組的所有節點座標。 BTW，chen兄上次提到的問題，我最近太忙還無法深入研究，請見諒! 感谢seansheu老大，没关系，我说过， 以你们的工作为重，小弟的问题现已解决了一部分，就是计算速度还需要提高，也就是要改默认的收敛准则，可对这方面还不太清楚～ 对于这个问题，我有点不太清楚的，那么你说是结点的坐标，那么三维单元有8个结点是不是全放进coords里了，那么在子程序调用的时候，我怎么知道是调用的哪个单元哪个结点的坐标呢, 我这有个高斯热源的DFLUX子程序，你看看，通过这个，实现了热源的移动，其中的a,b为特征半径，v为焊接速度你可以看看～ SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS, 1 JLTYP,TEMP,PRESS,SNAME) C INCLUDE 'ABA_PARAM.INC' C DIMENSION FLUX(2),COORDS(3),time(2) CHARACTER*80 SNAME I=3.1415 a=0.001 b=0.001 t=0.85 U=21 I=90 k=t*U*I/(PI*a*b) F=1.2 v=0.0026 hx=v*time(2) x=coords(1)-hx xx=x**2 yy=(coords(2))**2 flux(1)=k*f*exp(-f*((xx/a**2)+(yy/b**2))) RETURN END 子程序功能就是实现在已定义的单元集里从first point开始(,)，实现载荷移动，这里 有一段原文，是讲这个程序的实现的，烦请老大再帮忙看看，这阵子就为这个犯愁了～ dflux定义的热源是定义在单元积分点上的,this point指的是高斯积分点.你可以再用一个子程序,定义一个解相关变量,来输出定义的热源大小,这样再后处理中你就可以很清楚的看到.coords(3)是指的积分点的坐标吧,其实你要想看子程序是怎么调用的,你可以多定义些解相关的变量,用它们来判断某一时刻你调用的是哪个单元的,哪个积分点等等,而且这些都可以输出到后处理文件里面.你所说的载荷移动,实际上是定义的同一区域不同时刻的载荷分布来实现的.载荷的移动实现,也可以用noel,npt来判断的,比如,当到了某时刻的时候,你可以判断不同的单元定义不同的载荷. 哇，遇到大牛了～smarthere兄，有没有更详细一点的说明呢,例子也行～就重在看其他子程序变量调用的子程序，如何做呢,谢谢～另，你说用noel、npt来判断移动载荷，这句话如何解呢,恕小弟愚笨～ To use more than one user-defined mechanical material model, the variable CMNAME can be tested for different material names inside user subroutine UMAT as illustrated below: IF (CMNAME(1:4) .EQ. 'MAT1') THEN CALL UMAT_MAT1(argument_list) ELSE IF(CMNAME(1:4) .EQ. 'MAT2') THEN CALL UMAT_MAT2(argument_list) END IF UMAT_MAT1 and UMAT_MAT2 are the actual user material subroutines containing the constitutive material models for each material MAT1 and MAT2, respectively. Subroutine UMAT merely acts as a directory here. The argument list may be the same as that used in subroutine UMAT. 不知道这段对你有没有用,如果你需要用noel,npt来判断,我想这个也是可以实现的.当然还有其他的办法. 比如,你可以在你的子程序里面直接判断noel等,因为含有noel变量的子程序都是按照单元,积分点的序号来运行的,直到所有的执行完.就好比你定义了一个单元集,但是你只是需要里面某几个,或者有一定规律的几个单元在子程序里面有作用(这个在inp文件里面把单元选成相应集合也可做到),那么你在子程序里面就可以用noel来判断找出所需要的有规律的单元. 不知道我说的有没有道理,大家共同进步吧.有些时候你自己定义几个简单的子程序,运行简单的分析过程,你就会发现很多问题都明白了. 元旦快乐! 谢谢smarthere兄(姐)，你的意思小弟明白了～这就回去看看noel/npt方面的说明去～，它们的解释偶明白，noel是指的调用此子程序时的单元号、npt则是此单元的积分点号，现在不太清楚的是怎么能够在子程序里去调用它们，以及如何就你前面所说的将我这个当前子程序调用时单元号及积分点号能输出来, 还需要你来指点～ 多谢～ 元旦快乐～ 仿照上面的提示: sub aaa(***) if (noel.eq.'needed'.and.npt.eq.'needed') then call aaa1(***) write(6/7,*),noel,npt, endif end sub aaa1(***) **** *** end