ABAQUS用户子程序学习小结

1、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为13；3 子程序（*.for）文件中如何输出调试信息：write(6,*)coords(1),coor

2、ds(1)，在*.dat文件中可看到输出，如果希望write输出到msg文件中，则写为write(7,*)coords.；4 用户子程序dload中coords数组的含义：coords(1)也是一个数组，存贮单元集合中所有单元积分点的x坐标，coords(2)存贮y坐标，相应inp文件中的写法为：*dloadpy,pynu其中py为单元集合名称，定义方法为：*elset, elset=beam, generate1, 5, 1.*elset,elset=pybeam5 dload中f的定义方法：f只有定义在单元积分点上才有效，例如：f=1.0*coords (1)附一个简单实例：beam.in

3、p文件：*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=b311, 1,

4、22, 2, 33, 3, 44, 4, 55, 5, 6*elset, elset=beam, generate1, 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=rect0.2, 5.0.,0.,-1.*end instance*nset, nset=ends, instance=pa

5、rt-1-11, 6*nset, nset=_m4, internal, instance=part-1-16,*nset, nset=_m5, internal, instance=part-1-11,*end assembly* * materials* *material, name=steel*elastic210000., 0.3*elset,elset=pybeam* * boundary conditions* * name: disp-bc-1 type: symmetry/antisymmetry/encastre*boundary_m4, encastre* -* * st

6、ep: step-1* *step, name=step-1*static* * loads* * name: cforce-1 type: concentrated force*dloadpy,pynu* * output requests* * * field output: f-output-1* *output, field, variable=preselect* * field output: f-output-2* *output, field*element outputsf, * * history output: h-output-1* *output, history*n

7、ode output, nset=endscf1, cf2, cf3, cm1, cm2, cm3, rf1, rf2rf3, rm1, rm2, rm3, u1, u2, u3, ur1ur2, ur3*el print, freq=999999*node print, freq=999999*end stepbbb.for文件 subroutine dload(f,kstep,kinc,time,noel,npt,layer,kspt,coords, 1 jltyp,sname)c include aba_param.incc dimension time(2), coords (3) c

8、haracter*80 sname write(6,*)coords(3),coords(3) f=1.0*coords (1) return end运行方法：在abaqus command提示符后输入：abaqus job=beam user=bbb interactive子程序如下！如何编写，调研子程序？子程序usdfld中给固化度赋了一个初值110-4。subroutine usdfld(field,statev,pnewdt,direct,t,ceient,1 time,dtime,omname,orname,neield,nstatv,noel,npt,layer,2 kspt,ks

9、tep,kinc,ndi,nshr,coord,jmac,jmatyp,matlayo,laccfla)c includeaba_param.incc 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

10、field(1)=statev(1)c return end(2)disp在本例中用户子程序disp定义温度边界条件。 subroutine disp(u,kstep,kinc,time,node,noel,jdof,coords)c includeaba_param.incc 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)=3

11、63.+(time(2)-3600.)/30. else u(1)=398. end ifc 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 includeaba_param,incc dimension coords(3),time(2),field(nfield),h(2

12、) charactrr*80 samec h(1)=15 if(time(2).if.1800.) thensink=303.+time(2)/30.else if(time(2).if.3600.) thensink=363.else if(time(2).if.4650.) thensink=363.+(time(2)-3600.)/30elsesink=398.end ifcreturnend（4）hetval本例中用户子程序hetval定义复合材料内部的反应热。statev（1）即为从用户子程序usdfld中传递过来的固化度；statev（2）为固化率（固化度的变化率），由与固化度、温

13、度、升温速率等变量相关的固化动力学方程确定，决定了反应热的放热率；flux(1)表示反应热产生的单位时间、单位体积的热量。本例中，用euler法对固化动力学方程进行数值积分，得到的固化度再通过statev（1）传递回用户子程序usdfld，读者可以采用其他的数值积分方法。由于有些内容涉及相关部门的内部技术，所以未给出函数的具体形式。subroutine hetval(cmname,temp,time,dtime,statev,flux,1 fredef,dpred)cinclude aba_param,incccharacter*80 cmnamedimension temp(2),state

14、v(3),predef(1),time(2),flux(2),1 dpred(1)cif(temp(1).lt.363.15) thenstatev(2)=0.0elsestatev(2)=statev(1)=statev(1)+statev(2)*dtimeend ifflux(1)=242971740*statev(2)creturnend 我的多个子程序就是简单拼装到一个fortran文件的，一个子程序完毕，紧跟另一个子程序，实践证明：简单拼装是可行的，你之所以出错，我怀疑最大可能是源程序本身就有问题？还有，你的两个子程序都是dflux，如果像我这样第一个uflux后再接另一个uflux

15、，我觉得不可行，因为两个子程序同名，而简单拼装是将不同子程序放到一个fortran文件里，建议：要么只使用一个uflux，但根据不同情况调用不同热源，要么使用三个uflux，第一个用于调用另两个uflux用，可以参阅使用umat调用不同的材料本构第一个umat里只有如下语句：if (cmname(1:4) .eq. mat1) thencall umat_mat1(argument_list) else if(cmname(1:4) .eq. mat2) then call umat_mat2(argument_list) end if如材料名的前四字为mat1，则调用子程序umat_mat1

16、；否则，如为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

17、 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是指的什么点的坐标呢？

18、set是指的一组受热载的单元。 搞不太清楚，子程序是新手，哪位大侠做过，帮忙解释一下！ 關於this point指的就是你定義的節點， 每個節點都有座標定義，若是三維的幾何模型，那麼節點就會有x,y,z三個座標， 所以相對應於副程式中的coords array，那麼它就是1x3的array。 至於它其中說的geometry nonlinear需注意的事項， 是說一般若我們不考慮幾何非線性，那麼整個分析是屬於小變形，因為小變形所以忽略其位移量， 所以，在整個分析中導入副程式中的節點座標為原始定義座標，自始至終都是一樣。 但是，若整個分析為幾何非線性，那麼它會是大變形，所以其節點位移量就不可忽略。

19、 如此一來，在分析時導入副程式中的節點座標，會隨分析過程中的變化將位移量考慮進去，不斷更新。 對於你的問題，它傳遞的就是set這單元組的所有節點座標。 btw，chen兄上次提到的問題，我最近太忙還無法深入研究，請見諒! 感谢seansheu老大，没关系，我说过， 以你们的工作为重，小弟的问题现已解决了一部分，就是计算速度还需要提高，也就是要改默认的收敛准则，可对这方面还不太清楚！ 对于这个问题，我有点不太清楚的，那么你说是结点的坐标，那么三维单元有8个结点是不是全放进coords里了，那么在子程序调用的时候，我怎么知道是调用的哪个单元哪个结点的坐标呢？ 我这有个高斯热源的dflux子程序，你

20、看看，通过这个，实现了热源的移动，其中的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=90k=t*u*i/(pi*a*b) f=1.2 v=0.0026 hx=v*time(2) x=coor

21、ds(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)是指的积分点的坐标吧,其实你要想看子程序是怎么调用的,你可以多定义些解相关

22、的变量,用它们来判断某一时刻你调用的是哪个单元的,哪个积分点等等,而且这些都可以输出到后处理文件里面.你所说的载荷移动,实际上是定义的同一区域不同时刻的载荷分布来实现的.载荷的移动实现,也可以用noel,npt来判断的,比如,当到了某时刻的时候,你可以判断不同的单元定义不同的载荷. 哇，遇到大牛了！smarthere兄，有没有更详细一点的说明呢？例子也行！就重在看其他子程序变量调用的子程序，如何做呢？谢谢！另，你说用noel、npt来判断移动载荷，这句话如何解呢？恕小弟愚笨！ to use more than one user-defined mechanical material model

23、, 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 constitutiv