一个例题引发的困惑

yft117

最近，我在书上看到过一个例题，是关于轴的扭转分析，换了另外一种建摸方式，结果就相差一倍，问题也没检查出来，哪位大侠能否指点下。
原题条件为对一个直径为50mm，长度为120mm得实心轴在轴的两端施加扭矩1.5e3 N.m。求轴的最大剪应力。
例题程序为
/PREP7  
ET,1,PLANE42
ET,2,SOLID45
MP,EX,1,2.08e11
MP,PRXY,1,0.3  
rectng,0,0.025,0,0.12   
lesize,1,,,4
lesize,2,,,8
mshape,0
mshkey,1
amesh,1
extopt,esize,6  
extopt,aclear,1
vrotat,1,,,,,,1,4,360   
wprot,0,90  
cswpla,11,1,1,1
nsel,s,loc,x,0.025  
nrotat,all  
/solu   
d,all,ux
nsel,r,loc,z,-0.12  
f,all,fy,2500   
allsel,all  
da,2,all
da,6,all
da,10,all   
da,14,all   
SOLVE   
FINISH  
/POST1  
!*  
RSYS,11
AVPRIN,0,0  
AVRES,2,
/EFACET,1   
LAYER,0
PLESOL, S,XY, 0,1.0
所得结果为60.15
而我换种模型后所得的结果为他的两倍，
源程序如下：
/prep7
et,1,solid45
mp,ex,1,2.08e11
mp,prxy,1,0.3
CYL4,0,0,0.025, , , ,0.12   
wprota,,,90
vsbw,1
wprota,,90
vsbw,all
lsel,s,,,1,8,1
LESIZE,all, , ,5, , , , ,1
lsel,s,,,9,10,1
lsel,a,,,13,14,1
lesize,all,,,20,,,,,1
mshape,0
mshkey,2
vmesh,all
nummrg,all
numcmp,all
WPCSYS,-1,0
cswpla,11,1,1,1
nsel,s,loc,x,0.025
nrotat,all
d,all,ux
nsel,r,loc,z,0.12
f,all,fy,2500
allsel,all
nsel,s,loc,z,0,0
d,all,all
/solu
allsel,all
solve
/post1
RSYS,11
AVPRIN,0,0  
AVRES,2,
/EFACET,1   
LAYER,0
PLESOL, S,XY, 0,1.0

yan198202

由于你采用的45号单元退化导致的
如果你将45号单元改为95号单元，结果应该就差不多了

yan198202

对不起
弄错了
我试了一下
应该不是这个问题

yan198202

原来是你沿柱体方向划分数量不同造成的

yft117

数目是一样的，都是24个节点，我都数过几次了，估计不是这方面的原因

imlini

首先这个问题材料力学有理论解，是61Mpa。

   其次我把扭矩加载在整个圆柱面而不是一端的边上，这样得出的结果两种网格都是55Mpa左右，差别不大。

   接着尝试把单元换成有中间节点的solid186，使用第二种mesh，结果最大值变成了400Mpa左右。

   对比两种网格划分，在z方向上单元数量基本是1：2的关系，我觉得这个问题应该是应力集中造成的。因为LZ将这个大的力(2500N)加载在节点上，而且查看element result，所以在加载的地方是应力集中的。第一种情况因为单元大，所以应力集中影响不算大，在细分以后（第二种）单元变小，应力相差一倍。而将扭矩加载在整个圆柱面上时，应力集中会得到很大改善，所以两种网格计算结果差不多。

   结论：1 不要将力加载在线上或过少节点上， 以避免应力集中。
             2 不要盲目相信软件。软件没有错，就是看用户怎样使用了。
             3 约束圆柱面径向位移没有任何意义。
             4 尽信书不如无书。看帮助文件才是王道

[ 本帖最后由 imlini 于 2008-2-21 14:45 编辑 ]

yft117

非常感谢您的指导和教育，今后一定好好向您学习
对于您所说的结论我对此有不同观点。
1.体单元一般没有扭转的自由度，在施加扭矩时，容易出现约束方程不够而产生大量的警告，对结果精度有影响，因此，在施加扭矩时，一般需要对扭转自由度进行约束。
2，对于采用第二种模型，在轴径节点上施加扭矩，我曾经尝试过增加线段画分的单元数，结果单元数越多，结果越大，与理论解相差更远，不是你所说的由于节点过少造成应力集中而导致结果的错误。
3.将扭矩施加在整个面上，我尝试过，结果是55MPa左右，与你说的一致，但从剪应力云图上发现，云图不对称，与理论解有较大差距。
不知道您是否有时间，我的qq为156831846，有时间我们在网上详细讨论下这个问题。

[ 本帖最后由 yft117 于 2008-2-22 09:30 编辑 ]

imlini

2，对于采用第二种模型，在轴径节点上施加扭矩，我曾经尝试过增加线段画分的单元数，结果单元数越多，结果越大，与理论解相差更远，不是你所说的由于节点过少造成应力集中而导致结果的错误。
==================================
可能我表达的不清楚，我想说的是，由于应力＝力/面积，网格划分越细，单元面积越小，应力集中情况越严重。个人认为例题那样做不一定正确，只是凑巧和理论解接近。

不好意思我上班不能上qq，下班坚决不看专业相关

yft117

刚对命令流进行调试和修改，没出现警告，求得的结果为62.6，我所说的由于体模型自由度约束不够而造成大量的警告是不正确的，但不知道为啥，使用后处理方式，在剪应力云图中显示的最大剪应力很小，通过等效应力强度才得到62.6，接近理论解，但扭转角为1.876e-3，与理论解3.67e-3相差一半，不知道为哈，能解释下吗
关于例题上面说的方法，在论坛中那些功勋会员已经作了很详细的解释，他们也认为不正确的，但我看到几本书上的例题一般采用这种方法，肯定有它的道理，毕竟这些编书的估计在这方面肯定有很深入的研究，其应力应变云图都很接近理论情况。

hg_boy

闲着无聊，也做了一下：
前处理全部在hypermesh7.0里面完成。

单元类型：solid45;

把圆柱体的一个端面上节点的3个自由度都约束住。
扭矩的施加：

把另外一个端面上的节点都连接到一个independent node上，形成刚性区域。（Cerig）
扭矩直接施加在该independent node上。

[ 本帖最后由 hg_boy 于 2008-2-22 14:28 编辑 ]

hg_boy

shear yz:

hg_boy

von misee and  shear xz:

[ 本帖最后由 hg_boy 于 2008-2-22 14:06 编辑 ]

hg_boy

最大扭转角度：0.003698
（不晓得是度还是弧度）

hg_boy

剪应力结果不呈圆周对称，是因为结果是以直角坐标系显示的。

imlini

扭矩的施加：

把另外一个端面上的节点都连接到一个independent node上，形成刚性区域。（Cerig）
扭矩直接施加在该independent node上。
========================
请问是把1500Nm直接加载在一个这个节点上吗？需要对这个节点加位移约束吗？

hg_boy

原帖由 imlini 于 2008-2-22 15:54 发表
扭矩的施加：

把另外一个端面上的节点都连接到一个independent node上，形成刚性区域。（Cerig）
扭矩直接施加在该independent node上。
========================
请问是把1500Nm直接加载在一个这个节点上吗？ ...

扭矩直接施加在那个点上，你从我发的那个图也可以看的出来的。


不需要加约束，也不应该加约束，你吧他给约束死了，都动不了，力矩怎么传递啊。

ansys里面的CErig，作用类同与nastran里面的rbe2.

imlini

彻底的晕了，有人知道为什么要加刚性区域吗？为什么不这样做就不对呢？

hg_boy

所谓刚性区域，就是nastran里面的rbe2。
rbe2就是一个depnedent node,连接到很多independent node.
在程序内部，independent node 的位移被强制恒等于dependent node的位移。

这样，就相当于用刚度无穷大的一根杆或beam，把dependent node 和independent node连接起来了。

这样，你把扭矩作用在independent node上的时候，因为它是刚性连接到其他dependent node上的，扭矩就会牵动所有的independent node，这样，扭矩就施加到了整个圆柱的端面上了。

cerig不是唯一的办法。
用mpc单元也可以做到的。

yft117

非常感谢您的意见，您能给出您的命令流吗，我的经过修改后的命令流为
/prep7
et,1,45
et,2,21
r,1,0
mp,ex,1,2.08e11
mp,prxy,1,0.3
cylind,0.025,,0,0.12,0,360
wprot,,,90
vsbw,all
wprot,,90
vsbw,all
vatt,1,1
vmesh,all
n,100000,,,0.12
type,2
e,100000
nsel,s,loc,z,0.12
cerig,100000,all,uxyz
nsel,all
finish
/solu
f,100000,mz,1.5e3
allsel
nsel,s,loc,z,0
d,all,all
allsel
solve
finish

/post1
plnsol,s,eqv
set,1,last,0.5
plnsol,s,int
finish
不知道为啥得出的结果与您的差距如此大呢，我的syz为31.3，最小为-31.3，按照第三强度理论获得的最大应力为62.6，
最大位移为0.469e-4

imlini

set,1,last,0.5
＝＝＝＝＝＝＝＝
这句话的问题。改成set,1,last,1就好了，不知道为什么会加上这句话