均布扭转荷载如何在壳单元上施加？

graduates

一根箱形简支梁，长L=4000mm，受到均布扭转荷载Mt=2.5KN.m/m的作用，求其截面应力及位移。如果是单纯分析位移或内力的话，无疑用BEAM单元是最好的选择，可是如果要分析截面上应力的分布，就不得不考虑使用SHELL单元了。下面我就用shell63分析了这个问题，可是在加载时，就有点迷惑了，究竟怎么加呢？我想出来三种方法，得到的结果都不同，而且没有一个接近正确答案的，恳请大家指正。
fini
/cle
/filname,gzt-285,1
t1=6
t2=2
a=250+(t1+t2)/2
b=50+t1
L=4000
/prep7
et,1,shell63$keyopt,1,11,2
mp,ex,1,206e3
mp,prxy,1,0.3
R,1,T1
R,2,T2
k,1
k,2,b
k,3,b,a
k,4,,a
k,5,,,L
l,1,2$l,2,3$l,3,4$l,4,1
l,1,5
lsel,,,,1,4
adrag,all,,,,,,5
lsel,all
ldele,5,,,1
numcmp,all
/view,1,1,1,1
/rep
aplo

asel,,,,2,4
aatt,1,1,1
asel,inve
aatt,1,2,1
SAVE

ENUM=50                                !纵向所划分的单元数
lsel,,loc,z,l/2
lesize,all,,,ENUM
lsel,,loc,x,b/2$lesize,all,,,4
lsel,,loc,y,a/2$lesize,all,,,10
alls
mshkey,1$mshape,0,2d
amesh,all
!eplot

nsel,,loc,z$nsel,a,loc,z,l$d,all,ux,,,,,uy
nsel,,loc,z$nsel,r,loc,y,a/2$d,all,uz

!均布扭转荷载为2.5kN.m/m，将其转换为节点荷载——扭矩mz
!nsel,,loc,y,a$nsel,r,loc,x,b/2$nplo,1$f,all,mz,2.5e3*L/ENUM

!将每个节点扭矩mz转化为两个集中力加在腹板中点
!nsel,,loc,y,a/2$nsel,r,loc,x$f,all,fy,2.5e3*L/ENUM/b
!nsel,,loc,y,a/2$nsel,r,loc,x,b$f,all,fy,-2.5e3*L/ENUM/b

!将每个节点扭矩mz转化为两个集中力加在腹板与上翼缘交点处
nsel,,loc,y,a$nsel,r,loc,x$f,all,fy,2.5e3*L/ENUM/b
nsel,,loc,y,a$nsel,r,loc,x,b$f,all,fy,-2.5e3*L/ENUM/b

alls
/solu$solve$fini
/post1
nsel,,loc,z,l/2$prnsol,rot,z

!跨中转角位移理论解为0.01181rad

以跨中腹板中间节点（NODE=1208/545）的转角位移为标准，比较这三种加载方式的结果：

1.均布扭转荷载为2.5kN.m/m，将其转换为节点荷载——扭矩mz
        ROTZ=0.28501E-02

2.将每个节点扭矩mz转化为两个集中力加在腹板中点
        ROTZ=-0.32776E-02

3.将每个节点扭矩mz转化为两个集中力加在腹板与上翼缘交点处
        ROTZ=-0.33165E-02
从结果上看，这三种加载方案都不理想。


还有一个问题，就是求出的理论解是针对整个跨中截面的，而ANSYS分析出来的结果是跨中截面各个
节点的，且各节点不尽相同，取哪个节点为合适呢？上面作比较时我取的腹板中点处节点。