帮忙看看单刚计算有没问题，附源程序

ccdier

为验证自己编的单刚程序的有效性和准确性，和abaqus程序进行对比。。。。
abaqus程序：
建模，得到目标单元（4）的节点编号和坐标，输出单刚，如图：

另附，abaqus的inp文件，和单刚输出文件.
fortran程序:
将inp文件中的坐标和单元节点，弄成输入文件：

运行程序计算并进行对比，输出格式虽然不一样，但明显可以看出不同（左abaqus，右fortran）：

找了很久，没找到原因，可能当  局 者迷，希望大家帮忙看看什么原因。。
ps：坐标单位为mm

tonnyw

One part seems wrong to me.
          DO IEVAB=1,24
               DO ISTRE=1,6
                   BTDMATX(IEVAB,ISTRE)=0.0
                   DO JSTRE=1,6
                       BTDMATX(IEVAB,ISTRE)=BTDMATX(IEVAB,ISTRE)+BMATX(JSTRE,IEVAB)*DMATX(JSTRE,ISTRE)         !!!!============!!!
                   END DO
                END DO
          END DO

Is it supposed to be something like: BMATX(IEVAB, JSTRE)*DMATX(JSTRE,ISTRE)

ccdier

tonnyw 发表于 2012-9-1 20:32
One part seems wrong to me.
          DO IEVAB=1,24
               DO ISTRE=1,6

不是的。。。我求的是BT*D，BT是B矩阵的转置矩阵。。。所以是列。。

tonnyw

The comparison has to be an apple-to-apple one. Make sure the brick element in Abaqus is eight-noded conforming brick element using full integration.

ccdier

tonnyw 发表于 2012-9-1 22:40
The comparison has to be an apple-to-apple one. Make sure the brick element in Abaqus is eight-noded ...

abaqus里面是采用完全积分，而非减缩积分。。

tonnyw

What is exactly the element type you are using in Abaqus?

ccdier

tonnyw 发表于 2012-9-2 01:57
What is exactly the element type you are using in Abaqus?

** ELEMENT TYPE  C3D8

pasuka

http://forum.simwe.com/thread-1002748-1-1.html
既然会fortran就和上面帖子里面的代码结果比较一下，相对还简单些

hillyuan

C3D8 is not an element with full integration. It is reaaly a B-bar element

ccdier

hillyuan 发表于 2012-9-2 09:08
C3D8 is not an element with full integration. It is reaaly a B-bar element

那这些种不同的方法，得到的单刚数量级上会差多少？而且，abaqus的单刚第一个节点对角线为2113340,2113340,411...第一第二对角线位置一样，我又修改之后的fortran算出来的确不同。它是第一和第三对角线位置一样。。

tonnyw

3D linear brick element is supposed to be very simple element in terms of implementation. Even simpler than 2D element.

Why not form element stiffness matrix and do some computation to see if you can get correct results?

hillyuan

看了一下你的程序，似无大错。但调试程序是个精细活，当然也有一些技巧。如11楼，8`楼的建议都可考录。

题外话，程序的写法太老套。如

   DO I=1,6
      DO J=1,6
         DMATX(I,J)=0.0               
      END DO
   END DO

写成DMATX=0.0即可。这样程序会快一些，而且可以减少笔误。

Iorifo

一些意见：

（1）涉及积分点坐标诸如0.577350269789626，建议换做双精度（DOUBLE PRECISION or REAL*8）
     否则后面的数字都被截断了;
（2）SHAPEN(1)=(1-z1)*(1-z2)*(1-z3)/8，程序中整型，实型最好不要混合吧，建议1用1.0E/D0，8用
     8.0E/D0;
  (3) KSMATX(IEVAB,JEVAB)=KSMATX(IEVAB,JEVAB)+BTDBMATX(IEVAB,JEVAB) *DVOLU*1.0e3                  !刚度矩阵的单位为N/m,N/mm→N/m
       这里不用乘1.0E3，既然大家都是无量纲化的，这里就不需要特殊处理，商业软件是不会在类似这边
地方做这样的处理的，这就是为什么你的数量级要比ABAQUS高三个量级的原因；
（4）在测试的时候，可以考虑尽量简单一点的单元吧，比如可以用长*宽*高为1.0*1.0*1.0去测试，这样
     一看你的Jacobian行列式的值就知道程序在这段处理有没有处理对，实际上这样的单元的Jacobian行列式值为1.0。

ccdier

Iorifo 发表于 2012-9-3 18:19
一些意见：

（1）涉及积分点坐标诸如0.577350269789626，建议换做双精度（DOUBLE PRECISION or REAL*8）

1）这个我整个项目会全部采用real*8类型，就是浪费了点。。
2）同上，设置之后，常数就默认为双精度了。。
3）实际上不是都有单位的么。。而且之后的荷载也要相应的进行换算，统一就好了。
4）我的程序找到错的原因了。。修正之后和abaqus量级差不多。abaqus采用b-bar方法。
5）谢谢你的回答。。

tonnyw

Iorifo 发表于 2012-9-3 18:19
一些意见：

（1）涉及积分点坐标诸如0.577350269789626，建议换做双精度（DOUBLE PRECISION or REAL*8）

在测试的时候，可以考虑尽量简单一点的单元吧，比如可以用长*宽*高为1.0*1.0*1.0去测试
>>Maybe it should be 长*宽*高为2.0*2.0*2.0.

Iorifo

tonnyw 发表于 2012-9-3 20:48
在测试的时候，可以考虑尽量简单一点的单元吧，比如可以用长*宽*高为1.0*1.0*1.0去测试
>>Maybe it shoul ...

Yeah，你说得对，我弄错了：）

原来一直以为ABAQUS的C3D8是常规的Linear HEXA单元，看来确实用了B-Bar method，
我猜ABAQUS用的应该类似这本书里描述的，
http://solidmechanics.org/text/Chapter8_6/Chapter8_6.htm

lxj502

确实，c3d8单元和ansys中的solid185单元一样，都采用了b-bar方法，来改善体积锁定。二者单刚已对比，一模一样。