发布ABAQUS 的 UEL算例 3， 4 ，9节点INP和FOR文件

amani · 发表于 2010-5-23 11:02:14

本帖最后由 amani 于 2010-5-23 11:19 编辑

1# 蓝牙

你好，我看了你的平面四节点UEL程序，有两个问题，
1。为什么要用 DO Kintk=1,NNODE来控制整个循环呢？一般是用高斯积分点
   来控制的。你好像每个方向用了两个高斯积分点。
2。在最后，SSEMBLE RHS AND LHS
         DO K1=1,8
            DO K2=1,3
               RHS(K1,1)=RHS(K1,1)-B(K2,K1)*SSTRESS(K2)*THICK
            ENDDO
应该是 “RHS(K1,1)=RHS(K1,1)-B(K1,K2)*SSTRESS(K2)*THICK”吧？

            DO K2=1,8
               DO KK=1,3
                  AMATRX(K1,K2)=AMATRX(K1,K2)+S(K1,KK)*B(KK,K2)*THICK
               ENDDO
            ENDDO
         ENDDO
单元刚度矩阵和残余力RHS矩阵中为什么没有雅可比矩阵的行列式呢？
我看到书上是有的。

我认为：在ABAQUS里面，刚度矩阵K=-偏导（ RHS）/偏导（u）

浙江大学化工系  贝壳航母

蓝牙 · 发表于 2010-5-23 21:46:28

我的工作地点不允许用qq类软件
我留个邮箱给你，希望能和你多交流
caodongf@163.com

蓝牙 · 发表于 2010-5-23 22:08:27

本帖最后由蓝牙于 2010-5-23 22:28 编辑

1# 蓝牙

你好，我看了你的平面四节点UEL程序，有两个问题，
1。为什么要用 DO Kintk=1,NNODE来控制整个循环呢？一般是用高斯积分点
来控制的。你好像每个方向用了两个高斯积分点。
2。在最后，SSEMB ...
amani 发表于 2010-5-23 11:02

1.NNODE指的是节点数（不是NDOF自由度数），采用等参单元，节点数和高斯积分点数相同，所以用NNODE没有错啊？

2. RHS(K1,1)=RHS(K1,1)-B(K2,K1)*SSTRESS(K2)*THICK
本来是B的转置乘以SSTRESS的，我在这里偷了懒，没有求转置就直接乘了

3 B矩阵已经用到了雅克比矩阵了，理解为x和y是xi和eta的函数，等参元实际就是个黑盒子，输入的是x y，不管盒子里参与了多少运算，输出的依然是x y的函数，只是积分区间变成了以xi和eta为边界，变量也变成了xi和Eta

蓝牙 · 发表于 2010-5-23 22:19:16

本帖最后由蓝牙于 2010-5-23 22:52 编辑

不知道我讲清楚了没有，如果你依然有困惑，我可以发公式推导给你！

关于刚度矩阵K=-偏导（ RHS）/偏导（u）9 J* L/ h" e7 \9 P8 F% G* X; H, ~4 g
用户手册uel中Usage with general nonlinear procedures章节有关于刚度矩阵的描述
应该就是你表述的那样（至少在静态和修正的risk准静态下如此）

我的RHS就是通过刚度矩阵得到的，按照上面的公式反推，我的刚度矩阵也是对的！

amani · 发表于 2010-5-24 09:39:31

本帖最后由 amani 于 2010-5-24 09:42 编辑

其他所有我看懂了，还有一个问题：

我知道应变矩阵B的求解过程中用雅克比矩阵作了换算，但是在求单元刚度矩阵的积分的上下限是【-1.，1】，我觉得应该将DX,DY转化成D(XI)和D(eta), 因此应该多乘一个雅可比矩阵的行列式。这个在书上的公式很清楚的。

amani · 发表于 2010-5-24 10:28:31

另外你那个UEL程序能和RIKS方法结合吗？

蓝牙 · 发表于 2010-5-24 20:46:24

本帖最后由蓝牙于 2010-5-24 21:04 编辑

其他所有我看懂了，还有一个问题：
我知道应变矩阵B的求解过程中用雅克比矩阵作了换算，但是在求单元刚度矩阵的积分的上下限是【-1.，1】，我觉得应该将DX,DY转化成D(XI)和D(eta), 因此应该多乘一个雅可比矩阵的行列 ...
amani 发表于 2010-5-24 09:39

这是我的疏忽，我上传的for不是最终版，确实要乘以行列式的值，除此RHS表达式也要乘以行列式的值，我在下载文件中进行了更正，谢谢你的提醒

蓝牙 · 发表于 2010-5-24 21:22:32

本帖最后由蓝牙于 2010-5-24 21:25 编辑

另外你那个UEL程序能和RIKS方法结合吗？
amani 发表于 2010-5-24 10:28

RISKS方法在帮助手册中有说明
除了定义单元刚度矩阵外还要定义RHS（1），RHS（2），状态变量的更新，这样你就可以采用RIKS方法了，我现在在实现UEl的材料非线性，还在啃径向返回方面的书籍，暂时还没有想过RISK方法！但是我相信多花点时间也应该不是特别困难！

amani · 发表于 2010-5-24 21:25:48

本帖最后由 amani 于 2010-7-8 21:57 编辑

呵呵，我觉得你能写出来确实不错的，只是我对有限元很感兴趣，一直做有限元研究的，
如果你的UEL能够和RIKS方法结合的话，就更牛了，希望能和你进一步讨论。
我看了用户手册uel中Usage with general nonlinear procedures章节，有关于刚度矩阵的描述
全部是基于牛顿-拉夫森方法，遇到收敛性问题时，一般只能通过控制迭代时间不长或者引入粘性稳定性方法，这和米国田纳西大学高燕飞博士那篇改善cohesive的收敛性的论文思想是一致的。
最近，根据我多年的有限元分析经验，写了一篇评论性的复合材料有限元分析方面的论文，发表在英国的MATERIALS AND DESIGN 上，对有限元数值收敛性的问题的解决方法作了肤浅的综述评论，供大家欣赏。
Liu PF, Zheng JY. Recent developments on damage modeling and finite element analysis for composite laminates: a review. Materials and Design, 2010(31): 3825-3834.
浙江大学化工系贝壳航母

蓝牙 · 发表于 2010-5-24 21:36:09

本帖最后由蓝牙于 2010-5-24 21:40 编辑

呵呵，我觉得你能写出来确实不错的，只是我对有限元很感兴趣，一直做有限元研究的，
如果你的UEL能够和RIKS方法结合的话，就更牛了，希望能和你进一步讨论。
我看了用户手册uel中Usage with general nonlinear pro ...
amani 发表于 2010-5-24 21:25

amani真是牛人啊，一般写综述性文章的人，都是大牛，等文章出来后一定拜读！
在求解方法和稳定性方面，risk（弧长法确实有它长处，对此也是一知半解。
现在一边编程序一边恶补理论！

amani · 发表于 2010-5-24 21:54:27

30# 蓝牙

这个是我发表的评论性论文的ONLINE版本的，供欣赏

蓝牙 · 发表于 2010-5-25 07:48:06

文章已经下载
文章已经in press了，恭喜了
先大概的看了下，还没来得及细看，问个问题，在纤维复合材料进行多尺度的有限元仿真中，有人考虑微结构的尺寸效应吗？我知道在颗粒复合复合材料中，因为经典的弹塑性理论中是不包含尺度项的，经典的弹塑性理论也就不在适用！

amani · 发表于 2010-5-25 13:11:02

尺寸效应很多的，我是多尺度，跟尺寸效应不是一回事。

amani · 发表于 2010-5-25 13:12:36

本帖最后由 amani 于 2010-5-26 19:39 编辑

28# 蓝牙
我还有一个问题，你的FOR程序感觉还是有些问题

      Gauss(1,1)=SQRT(1/3.0D0)*-1.0D0
      Gauss(1,2)=SQRT(1/3.0D0)*-1.0D0
      Gauss(2,1)=SQRT(1/3.0D0)
      Gauss(2,2)=SQRT(1/3.0D0)*-1.0D0
      Gauss(3,1)=SQRT(1/3.0D0)
      Gauss(3,2)=SQRT(1/3.0D0)
      Gauss(4,1)=SQRT(1/3.0D0)*-1.0D0
      Gauss(4,2)=SQRT(1/3.0D0)*1.0D0

   DO Kintk=1,NNODE
         xi=Gauss(Kintk,1)
         eta=Gauss(Kintk,2)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       Make B Matrix  in Natural Coordinate System                C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
         DO I=1,2
         DO J=1,NNODE
            BN(I,J)=ZERO
         ENDDO
         ENDDO
C
         BN(1,1)=0.25D0*(-1+eta)
         BN(1,2)=0.25D0*(1-eta)
         BN(1,3)=0.25D0*(1+eta)
         BN(1,4)=0.25D0*(-1-eta)
         BN(2,1)=0.25D0*(-1+xi)
         BN(2,2)=0.25D0*(-1-xi)
         BN(2,3)=0.25D0*(1+xi)
         BN(2,4)=0.25D0*(1-xi)

应该是下面更清晰一些，每个方向只有2个积分点。4个形状函数，每个形状函数只有4种组合。8个形状函数对局部坐标的导数函数，每个只有2种组合。
      Gauss(1,1)=SQRT(1/3.0D0)
      Gauss(1,2)=-SQRT(1/3.0D0)
      Gauss(2,1)=SQRT(1/3.0D0)
      Gauss(2,2)=-SQRT(1/3.0D0)

   DO Kk1=1,2
         DO,kk2=1,2
         xi=Gauss(kk1,kk2)
      eta=Gauss(kk1,kk2)
C
         DO I=1,2
         DO J=1,NNODE
            BN(I,J)=ZERO
         ENDDO
         ENDDO
C
         BN(1,1)=0.25D0*(-1+eta)
         BN(1,2)=0.25D0*(1-eta)
         BN(1,3)=0.25D0*(1+eta)
         BN(1,4)=0.25D0*(-1-eta)
         BN(2,1)=0.25D0*(-1+xi)
         BN(2,2)=0.25D0*(-1-xi)
         BN(2,3)=0.25D0*(1+xi)
         BN(2,4)=0.25D0*(1-xi)

蓝牙 · 发表于 2010-5-28 20:17:34

当然可以，目的是一样的啊？

rmshao · 发表于 2010-5-28 23:43:34

看不了附件，那位好心的给我发一个吧，多谢
345853811@qq.com

amani · 发表于 2010-5-29 07:34:48

我运行了你的程序，好像有点问题嘛。没有结果出来

蓝牙 · 发表于 2010-5-29 08:07:15

你的abaqus安装是否通过了verification？

amani · 发表于 2010-5-29 11:16:40

本帖最后由 amani 于 2010-5-29 11:44 编辑

您好！我的安装没有问题的。我运行了您的UEL程序，四节点的有问题，输出的单元刚度矩阵全部是零，请你检查一下把。9节点的可以算出结果。另外我觉得高斯积分点的数目并没有必要和节点数目相同的。再问一个问题，uel的后处理有没有方便一点的软件？

洪亮 · 发表于 2010-5-29 13:44:20

请问各位“uel”的具体含义是什么啊？

账号		自动登录	找回密码
密码			注册

账号		自动登录	找回密码
密码			立即注册

[子程序] 发布ABAQUS 的 UEL算例 3， 4 ，9节点INP和FOR文件

本帖子中包含更多资源

本帖子中包含更多资源