今天木木给大家带来的是梁单元的UEL子程序,作为杆系单元重要的组成成员,该部分是得好好讲讲了,而且要全面的讲,尽木木的知识能力所及讲的透彻,让刚入手UEL的小白们看的“起劲”!参考文档和源文件可在后台回复:梁单元UEL即可。
今天分享的内容主要有:梁单元介绍、等效节点荷载处理、状态变量的自定义输出、理论解对比等等,精彩超多,请进一步往下翻阅~
梁单元的讲解在前几期的有限元基础编程的推文中有详细讲述,感兴趣的小伙伴可以跳转相关页面查阅。
案例简介
我们今天讲的案例如图二所示,梁长
这里我们简单讲述一下变形协调条件,图3所示(图2分解)为一承受均布荷载的悬臂梁,在B端点产生的位移
由以上公式可得到右端点节点力、弯矩以及左端点节点力。
为了程序的易读性,我们可以将经常使用的的计算部分编写成一个子程序,使用时直接call
一下即可,比如在UMAT
中,我们会常常用到Mises应力公式,所以我们可以将Mises的求解编写为一个子程序,在UMAT中调用即可,大大减少了代码量,更能提升易读性。
本程序中,我编写了矩阵置零函数,用于数组的初始化,大家在自己的程序中编制矩阵相乘函数和矩阵求逆函数,大大提升编写代码的效率。
call K_MATRIX_ZERO(AMATRX,ndofel,ndofel)
call K_MATRIX_ZERO(RHS,ndofel,nrhs)
call K_MATRIX_ZERO(SVARS,4,1)
C---------------------矩阵置零------------------------
subroutine K_MATRIX_ZERO (A,n,m)
INCLUDE 'ABA_PARAM.INC'
dimension A(n,m)
do i=1,n
do j=1,m
A(i,j) = 0.d0
end do
end do
return
end
在本次案例中,我们注意到,梁体受到向下的均布荷载,以往的UEL案例分析中,我们没考虑过这种面力、体力的作用,可以参考一下势能泛函:
具体含义可见附带文档。本案例的面力处理方式:将其转化为单元节点荷载,转化的原理可参考附带文档,也可查阅有限元相关教材,推荐曾攀老师的《有限元基础教程》。
H(1)=-P*DL/TWO
H(2)=-P*DL2/TWELVE
H(3)=-P*DL/TWO
H(4)=P*DL2/TWELVE
!---注意RHS的变化------------
do I=1,NDOFEL
RHS(I,1)=RHS(I,1)-H(I)
end do
为了与理论解析解做对比,可以以状态变量的形式将右端点的节点力和弯矩以及左端点的节点力输出到指定位置,这里输出到指定文件中。涉及到Fortran语言中文件的使用,编制了计算节点力的子程序,如下所示。
CALL Node_Force(SVARS,AMATRX,U,H)
!-----------------------------
SUBROUTINE Node_Force(SVARS,AMATRX,U,H)
INCLUDE 'ABA_PARAM.INC'
DOUBLE PRECISION,DIMENSION(4):: H,SVARS,U
DIMENSION AMATRX(4,4)
SVARS(1:4)=SVARS(1:4)+MATMUL(AMATRX,U(1:4))
SVARS(1)=SVARS(1)+H(1)
SVARS(2)=SVARS(2)+H(2)
SVARS(3)=SVARS(3)+H(3)
SVARS(4)=SVARS(4)+H(4)
RETURN
END
状态变量输出
OPEN(UNIT=16,FILE='D:\TEMP\STRESS.DAT',STATUS='UNKNOWN',
1 FORM='FORMATTED')
IF(JELEM .EQ. 1) THEN
! JELEM为单元编号
WRITE(16,*) '-----------NEXT STEP-------------'
END IF
WRITE(16,*) 'ELEMENT NUMBER',JELEM
WRITE(16,*) 'NODE',JELEM,'FORCE',SVARS(1)
WRITE(16,*) 'NODE',JELEM,'MOMENT',SVARS(2)
WRITE(16,*) 'NODE',JELEM+1,'FORCE',SVARS(3)
WRITE(16,*) 'NODE',JELEM+1,'MOMENT',SVARS(4)
WRITE(16,*) ' '
即可在'D:\TEMP'位置处创建STRESS.DAT文件,按照以上语句进行写出。
与解析解一致。最后来看一下位移云图吧!
今天的分享就止于此了,下期我们再见啦~~喜欢本次案例的小伙伴可以点点小赞和在看哦,欢迎转发至朋友圈,让更多的对二次开发感兴趣的伙伴看到。