论坛 › J01：有限元数值方法 › Smith的程序与ANSYS程序计算结果不一致

grandpollux 发表于 2012-3-10 22:25:05

Smith的程序与ANSYS程序计算结果不一致

最近，我研究了一下《有限元方法编程》第二版（Smith与Griffiths合编）中的例程5.2。对于书中的问题，我采用了ANSYS重新进行了分析，以下是分析使用的命令流：

finish
/clear

/prep7
et,1,plane42,,,2 !平面应变问题
mp,ex,1,1.E6
mp,prxy,1,0.3

n, 1,0.0000E+00,0.0000E+00
n, 2,0.0000E+00,-0.5000E+01
n, 3,0.0000E+00,-0.1000E+02
n, 4,0.1000E+02,0.0000E+00
n, 5,0.1000E+02,-0.5000E+01
n, 6,0.1000E+02,-0.1000E+02
n, 7,20.0,0.0
n, 8,0.2000E+02,-0.5000E+01
n, 9,0.2000E+02,-0.1000E+02
n, 10,0.3000E+02,0.0000E+00
n, 11,0.3000E+02,-0.5000E+01
n, 12,0.3000E+02,-0.1000E+02

e,2,5,4,1
e,3,6,5,2
e,5,8,7,4
e,6,9,8,5
e,8,11,10,7
e,9,12,11,8

/solu
antype,static
f,1,fy,-5.0
f,4,fy,-5.0

d,1,ux,0
d,2,ux,0
d,3,all,0
d,6,all,0
d,9,all,0
d,12,all,0
d,10,ux,0
d,11,ux,0
iswrite,on
solve
finish

/post1
presol,s,x
presol,s,y
presol,s,xy

其计算的结点位移结果为：
NODE      UX   
       1   0.0000   
       2   0.0000   
       3   0.0000   
       4 -0.27011E-06
       50.15050E-05
       6   0.0000   
       70.31518E-06
       80.56193E-06
       9   0.0000   
      10   0.0000   
      11   0.0000   
      12   0.0000
NODE      UY   
       1 -0.90221E-05
       2 -0.42485E-05
       3   0.0000   
       4 -0.37089E-05
       5 -0.18870E-05
       6   0.0000   
       70.79676E-06
       80.26713E-06
       9   0.0000   
      10 -0.10849E-07
      110.59749E-07
      12   0.0000
单元积分点应力为：
!INITIAL STRESS RECORD FOR ELEMENT       2
!   SX            SY            SZ          SXY            SYZ         SXZ
eis,       2,1
-0.300421   -0.932356   -0.369833      0.103297       0.00000       0.00000   
-0.300421   -0.632702   -0.279937      0.103297       0.00000       0.00000   
-0.204939   -0.632702   -0.251292      0.103297       0.00000       0.00000   
-0.204939   -0.932356   -0.341189      0.103297       0.00000       0.00000   
!
!INITIAL STRESS RECORD FOR ELEMENT       1
!   SX            SY            SZ          SXY            SYZ         SXZ
eis,       1,1
-0.241079      -1.03949   -0.384171      0.793198E-01   0.00000       0.00000   
-0.241079   -0.664949   -0.271808      0.793198E-01   0.00000       0.00000   
-0.353698   -0.664949   -0.305594      0.793198E-01   0.00000       0.00000   
-0.353698      -1.03949   -0.417956      0.793198E-01   0.00000       0.00000   
!
!INITIAL STRESS RECORD FOR ELEMENT       4
!   SX            SY            SZ          SXY            SYZ         SXZ
eis,       4,1
-0.127013   -0.381934   -0.152684      0.120921       0.00000       0.00000   
-0.127013   -0.108594   -0.706821E-010.120921       0.00000       0.00000   
-0.186843   -0.108594   -0.886311E-010.120921       0.00000       0.00000   
-0.186843   -0.381934   -0.170633      0.120921       0.00000       0.00000   
!
!INITIAL STRESS RECORD FOR ELEMENT       3
!   SX            SY            SZ          SXY            SYZ         SXZ
eis,       3,1
-0.147111   -0.333461   -0.144171      0.503106E-01   0.00000       0.00000   
-0.147111   -0.350822E-01 -0.546579E-010.503106E-01   0.00000       0.00000   
-0.501469E-01 -0.350822E-01 -0.255687E-010.503106E-01   0.00000       0.00000   
-0.501469E-01 -0.333461   -0.115082      0.503106E-01   0.00000       0.00000   
!
!INITIAL STRESS RECORD FOR ELEMENT       6
!   SX            SY            SZ          SXY            SYZ         SXZ
eis,       6,1
-0.113808E-020.409506E-010.119437E-010.176247E-01   0.00000       0.00000   
-0.113808E-020.146360E-010.404937E-020.176247E-01   0.00000       0.00000   
-0.367899E-010.146360E-01 -0.664618E-020.176247E-01   0.00000       0.00000   
-0.367899E-010.409506E-010.124819E-020.176247E-01   0.00000       0.00000   
!
!INITIAL STRESS RECORD FOR ELEMENT       5
!   SX            SY            SZ          SXY            SYZ         SXZ
eis,       5,1
-0.403818E-010.745721E-010.102571E-01 -0.290092E-01   0.00000       0.00000   
-0.403818E-01 -0.159042E-02 -0.125917E-01 -0.290092E-01   0.00000       0.00000   
-0.247269E-01 -0.159042E-02 -0.789518E-02 -0.290092E-01   0.00000       0.00000   
-0.247269E-010.745721E-010.149536E-01 -0.290092E-01   0.00000       0.00000
而Smith的程序计算得到的结果为：
Node         Displacement
    1   0.0000E+00 -0.8601E-05
    2   0.0000E+00 -0.4056E-05
    3   0.0000E+000.0000E+00
    4    -0.5472E-07 -0.3771E-05
    5   0.1225E-05 -0.1910E-05
    6   0.0000E+000.0000E+00
    7   0.1792E-060.5860E-06
    8   0.6280E-060.1707E-06
    9   0.0000E+000.0000E+00
   10   0.0000E+000.1129E-06
   11   0.0000E+000.1066E-06
   12   0.0000E+000.0000E+00
The Gauss Point stresses for element    1   are :
Point    1
-0.4299E+00 -0.1058E+010.1431E+00
Point    2
-0.2511E+00 -0.6411E+000.8629E-01
Point    3
-0.3304E+00 -0.1016E+010.8352E-01
Point    4
-0.1516E+00 -0.5985E+000.2667E-01
The Gauss Point stresses for element    2   are :
Point    1
-0.2856E+00 -0.9141E+000.8499E-01
Point    2
-0.1427E+00 -0.5806E+000.1394E+00
Point    3
-0.3808E+00 -0.9549E+000.3735E-01
Point    4
-0.2379E+00 -0.6215E+000.9177E-01
The Gauss Point stresses for element    3   are :
Point    1
-0.1513E+00 -0.3680E+000.6412E-01
Point    2
0.2766E-03 -0.1429E-010.1010E+00
Point    3
-0.2159E+00 -0.3957E+000.1359E-01
Point    4
-0.6432E-01 -0.4197E-010.5050E-01
The Gauss Point stresses for element    4   are :
Point    1
-0.2331E+00 -0.4231E+000.1477E+00
Point    2
-0.9447E-01 -0.9964E-010.1211E+00
Point    3
-0.1867E+00 -0.4032E+000.1015E+00
Point    4
-0.4805E-01 -0.7974E-010.7493E-01
The Gauss Point stresses for element    5   are :
Point    1
0.1061E-020.7274E-01 -0.4210E-01
Point    2
-0.2619E-010.9144E-02 -0.2217E-01
Point    3
-0.3382E-010.5779E-01 -0.3302E-01
Point    4
-0.6108E-01 -0.5806E-02 -0.1308E-01
The Gauss Point stresses for element    6   are :
Point    1
-0.4855E-010.1373E-010.3616E-01
Point    2
-0.5281E-010.3775E-020.8267E-02
Point    3
0.2640E-030.3465E-010.3758E-01
Point    4
-0.4001E-020.2469E-010.9688E-02
显然，上述两程序结果是不一样的，请问这是为什么？

grandpollux 发表于 2012-3-11 13:31:45

自已顶一个。
哪位高手能给点线索？？？

cobalt 发表于 2012-3-11 18:35:01

从位移的对比来看，两个软件给出的结果数据都很小，都是E-6或者更小，我感觉这种对比说明不了什么，起码我觉得不能证明Smith的程序有错误。
结果太小了，很有可能是计算机存储误差造成的。
我遇到过这样的例子，6年前某个程序的计算结果与现在重新计算的结果对比存在差异，这话差异就是E-6数量级的，我感觉是不同的计算机或者操作系统的差异造成的。另外，我也遇到过使用不同版本编译器编译的程序计算结果也存在差异的情况。
所以，结果楼主再做一个例子验证一下。

grandpollux 发表于 2012-3-12 15:56:05

本帖最后由 grandpollux 于 2012-3-12 15:56 编辑

谢谢楼上的回复。
但我不认为是计算机造成的精度问题。
刚用ANSYS对《有限元方法编程》（Smith与Griffiths合编）第四版中的图5.11中的问题进行了分析，发现ANSYS的计算结果仍与书中程序的计算结果不同。
个人认为：ANSYS分析时所使用的Plane42缺省状态下使用了非协调元思想，不是书中一般意义上的四节点单元。所以两者位移解答是不一样的。但我不确定这个判断，也不懂得如何关闭或打开ANSYS中的非协调元选项。

jay5109 发表于 2012-3-14 18:43:45

受到存储精度和求解起差别等影响，这样的对比出来的结果不可能完全一致，楼主可以试这对比下两者的分布规律

unipx 发表于 2012-3-20 22:19:37

我用ABAQUS平面应变4点积分单元（不是非协调元也不是一点积分单元）算了一下这个例子，结果如下：
*Node
      1,          30.,         0.
      2,         0.,         0.
      3,          30.,          10.
      4,         0.,          10.
      5,          20.,         0.
      6,          10.,         0.
      7,          30.,         5.
      8,          10.,          10.
      9,          20.,          10.
   10,         0.,         5.
   11,   20.0359497,   4.99692011
   12,   10.0041866,   4.99927092

            结点            U.U1            U.U2
            标签          @Loc 1          @Loc 1
-------------------------------------------------
               1   68.3168E-39   136.642E-39
               2   -238.85E-39    -3.96193E-36
               3   262.878E-39    -229.254E-09
               4    -1.40646E-36    -9.64107E-06
               5   1.01298E-36    -1.03976E-36
               6   1.62306E-36    -5.13496E-36
               7    -118.964E-39   59.3522E-09
               8    -159.678E-09    -3.60043E-06
               9   256.926E-09   1.10702E-06
            10    -1.20296E-36    -4.42389E-06
            11   558.743E-09   357.068E-09
            12       1.523E-06    -1.88429E-06
发现ABAQUS的结果与ANSYS的类似（还略小一点），而我用自己写的程序算一下，发现与smith的结果几乎一样。
曾经比较过ABAQUS的单元刚阵，发现SMITH程序的平面应变单元刚阵与ABAQUS是一样的，然而结果就是不一样。
我赞同楼主，也不认为这是求解器的问题，不信做一个平面应力的简单算例，发现自己写的程序可以与商业软件的结果几乎完全一样。

也许商业软件的奇妙之处就在于它总有一些你搞不懂的地方。

tonnyw 发表于 2012-3-21 04:32:13

In Ansys, by default the plane strain element is incompatible type. It would be better if you can get the results organized. Any time, you make a comparison between different softwares, you have to be sure that all the input are the same.

pasuka 发表于 2012-3-21 09:36:11

翻翻弹性力学书籍，找一个有精确解的例子，不同的网格密度再用二种程序（Q4协调元）计算一下，看看哪个收敛性更好

grandpollux 发表于 2012-4-18 21:36:06

感谢楼上各位的回复。特别是tonnyw的那张图片，让人受教了。

grandpollux 发表于 2012-4-27 20:43:58

我刚才在ANSYS中重算了前面的算例，但在设置单元的命令中增加了一个参数，
et,1,plane42,,1,2   ！参数1代表不使用内部自由度，即单元是传统教科书的四节点单元，而不是ANSYS缺省的非协调元！
结果表明：ANSYS计算的结果与Smith程序的计算结果吻合。

再次感谢tonnyw的帮助。

查看完整版本: Smith的程序与ANSYS程序计算结果不一致