【讨论】请大家看看我的问题，奇怪的错误，郁闷死了

Mrice

大家好，我现在做一个物理气相沉积过程的温度场和应力场的模拟问题，物理过程是这样的，在一块基板上交替沉积两种材料，形成两种材料的多多层复合涂层，其中基板的温度为900度，两种材料的沉积到基板的初始温度分别为2680度和1600度，我想求解沉积过程的温度场和沉积完毕后的应力场。我现在用的是间接的方法，即先算温度场，在计算温度场过程中使用单元生死来模拟沉积生长过程，然后将最后一个载荷步的温度场作为体载荷加到模型中，进行热应力的计算，我在温度的计算过程中，没什么问题，但是在计算应力的时候总是出现这样的错误：
Error:
      negative radius element 3401
其中3401号单元位于模型的最左上角，这个问题困扰我两个星期了，真的不知道是怎么回事儿，我试着调节网格的密度和材料的属性，但问题都解决不了，始终是最左上角的单元出现错误，在结果中，有位移解，但没有应力，应变解。
     我采用的是轴对称模型，基板上沉积四层材料，下面是我的命令流，真的希望板主和高手能给些意见，谢谢了先
  
/PREP7   
!*   
ET,1,PLANE55
!*   
KEYOPT,1,1,0
KEYOPT,1,3,1
KEYOPT,1,4,0
KEYOPT,1,8,0
KEYOPT,1,9,0
!*   
/REPLOT,RESIZE   
/REPLOT,RESIZE   
!*   
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,400
MPTEMP,3,800
MPDATA,DENS,1,,8220  
MPDATA,DENS,1,,8220  
MPDATA,DENS,1,,8220  
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,400
MPTEMP,3,800
UIMP,1,REFT,,,25
MPDE,ALPX,1  
MPDE,ALPY,1  
MPDE,ALPZ,1  
MPDATA,ALPX,1,,14.4e-6   
MPDATA,ALPX,1,,14.4e-6   
MPDATA,ALPX,1,,14.4e-6   
!*   
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,400
MPTEMP,3,800
MPDE,KXX,1   
MPDE,KYY,1   
MPDE,KZZ,1   
MPDATA,KXX,1,,11.5   
MPDATA,KXX,1,,17.3   
MPDATA,KXX,1,,23.8   
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,400
MPTEMP,3,800
MPTEMP,4,1200   
MPDATA,ENTH,1,,0
MPDATA,ENTH,1,,1328557500   
MPDATA,ENTH,1,,3051469500   
MPDATA,ENTH,1,,5113045500   
!*   
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,400
MPTEMP,3,800
MPTEMP,4,1200   
UIMP,2,REFT,,,1600   
MPDE,ALPX,2  
MPDE,ALPY,2  
MPDE,ALPZ,2  
MPDATA,ALPX,2,,11.6e-6   
MPDATA,ALPX,2,,14.0e-6   
MPDATA,ALPX,2,,16.0e-6   
MPDATA,ALPX,2,,20.8e-6   
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,400
MPTEMP,3,800
MPTEMP,4,1200   
UIMP,2,REFT,,,1600   
MPDE,ALPX,2  
MPDE,ALPY,2  
MPDE,ALPZ,2  
MPDATA,ALPX,2,,11.6e-6   
MPDATA,ALPX,2,,14.0e-6   
MPDATA,ALPX,2,,16.0e-6   
MPDATA,ALPX,2,,20.8e-6   
!*   
!*   
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,400
MPTEMP,3,800
MPTEMP,4,1200   
MPDATA,DENS,2,,7320  
MPDATA,DENS,2,,7320  
MPDATA,DENS,2,,7320  
MPDATA,DENS,2,,7320  
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,400
MPTEMP,3,800
MPTEMP,4,1200   
MPDE,KXX,2   
MPDE,KYY,2   
MPDE,KZZ,2   
MPDATA,KXX,2,,4.3   
MPDATA,KXX,2,,6.4   
MPDATA,KXX,2,,10.2   
MPDATA,KXX,2,,16.1   
!*   
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,400
MPTEMP,3,800
MPTEMP,4,1200   
MPTEMP,5,1550   
MPTEMP,6,1600   
MPTEMP,7,1800   
MPDATA,ENTH,2,,0
MPDATA,ENTH,2,,1375245000   
MPDATA,ENTH,2,,3108621000   
MPDATA,ENTH,2,,5395389000   
MPDATA,ENTH,2,,7352757000   
MPDATA,ENTH,2,,9.5575e9  
MPDATA,ENTH,2,,1.0676e10
!*   
!*   
MPTEMP,,,,,,,,   
MPTEMP,1,0   
MPDATA,DENS,3,,5650  
MPTEMP,,,,,,,,   
MPTEMP,1,0   
UIMP,3,REFT,,,2680   
MPDE,ALPX,3  
MPDE,ALPY,3  
MPDE,ALPZ,3  
MPDATA,ALPX,3,,7.72e-6   
MPTEMP,,,,,,,,   
MPTEMP,1,0   
MPDE,KXX,3   
MPDE,KYY,3   
MPDE,KZZ,3   
MPDATA,KXX,3,,1.05   
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,2630   
MPTEMP,3,2680   
MPTEMP,4,2880   
MPDATA,ENTH,3,,0
MPDATA,ENTH,3,,7108914750   
MPDATA,ENTH,3,,11106289750   
MPDATA,ENTH,3,,11652079750   
!*   
!*   
  
RECTNG,0,0.01,0,0.001,   
RECTNG,0,0.01,0.001,0.00105,
RECTNG,0,0.01,0.00105,0.0011,   
RECTNG,0,0.01,0.0011,0.00115,   
RECTNG,0,0.01,0.00115,0.0012,   
FLST,2,5,5,ORDE,2   
FITEM,2,1   
FITEM,2,-5   
AGLUE,P51X   
  
*SET,_zcf(1),32,0,   
!*   
/SHRINK,0   
/ESHAPE,0.0  
/EFACET,1   
/RATIO,1,1,5
/CFORMAT,32,0   
/REPLOT  
!*   
CM,_Y,AREA   
ASEL, , , ,       1  
CM,_Y1,AREA  
CMSEL,S,_Y   
!*   
CMSEL,S,_Y1  
AATT,       1, ,   1,       0,   
CMSEL,S,_Y   
CMDELE,_Y   
CMDELE,_Y1   
!*   
FLST,5,2,5,ORDE,2   
FITEM,5,6   
FITEM,5,8   
CM,_Y,AREA   
ASEL, , , ,P51X  
CM,_Y1,AREA  
CMSEL,S,_Y   
!*   
CMSEL,S,_Y1  
AATT,       2, ,   1,       0,   
CMSEL,S,_Y   
CMDELE,_Y   
CMDELE,_Y1   
!*   
FLST,5,2,5,ORDE,2   
FITEM,5,7   
FITEM,5,9   
CM,_Y,AREA   
ASEL, , , ,P51X  
CM,_Y1,AREA  
CMSEL,S,_Y   
!*   
CMSEL,S,_Y1  
AATT,       3, ,   1,       0,   
CMSEL,S,_Y   
CMDELE,_Y   
CMDELE,_Y1   
!*   
FLST,5,5,4,ORDE,5   
FITEM,5,3   
FITEM,5,7   
FITEM,5,11   
FITEM,5,15   
FITEM,5,19   
CM,_Y,LINE   
LSEL, , , ,P51X  
CM,_Y1,LINE  
CMSEL,,_Y   
!*   
LESIZE,_Y1, , ,200, , , , ,1
!*   
MSHAPE,0,2D  
MSHKEY,1
!*   
FLST,5,5,5,ORDE,3   
FITEM,5,1   
FITEM,5,6   
FITEM,5,-9   
CM,_Y,AREA   
ASEL, , , ,P51X  
CM,_Y1,AREA  
CHKMSH,'AREA'   
CMSEL,S,_Y   
!*   
AMESH,_Y1   
!*   
CMDELE,_Y   
CMDELE,_Y1   
CMDELE,_Y2   
!*   
  
！以上为材料的热物理性能参数及建模分网过程
  
！以下为模拟层的生长过程引入的一些变量和数组
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!将对称轴处材料2，3的节点Y坐标按照从小到大的顺序
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!存储到数组layer_start_node中
  
nsel,all
*get,nmax,node,,num,max          !定义总的节点数
*get,minx,node,,mnloc,x
*get,maxx,node,,mxloc,x
*get,miny,node,,mnloc,y
*get,maxy,node,,mxloc,y
dep_conv=200
cool_conv=25
dep_temp=900
cool_temp=25
  
mat2=2                            !2号材料为NiCoCrAlY合金
mat3=3                            !3号材料为ZrO2陶瓷
esel,s,mat,,mat2                  !选择2号材料单元
esel,a,mat,,mat3                  !选择3号材料单元
nsle                              
nsel,r,loc,x,minx                 !选择对称轴处材料2，3的节点
          
*get,layer_no,node,,count         !返回对称轴处材料2，3的节点数
!定义每层单元起始节点数组
*dim,layer_start_node,array,layer_no   
jjt=0  
!将对称轴处材料2，3的节点号写入数组                     
*do,i,1,nmax                        
*if,nsel(i),eq,1,then
jjt=jjt+1
layer_start_node(jjt)=i
*endif
*enddo
*do,i,1,layer_no-1               !将layer_start_no数组按照y坐标从小到大排序
*do,j,i+1,layer_no
jjti=layer_start_node(i)
jjtj=layer_start_node(j)
*if,ny(jjti),gt,ny(jjtj),then
jjt=layer_start_node(i)
layer_start_node(i)=layer_start_node(j)
layer_start_node(j)=jjt
*endif
*enddo
*enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!下面程序的作用是将每层单元的单元号存入
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!到二维数组layer_elem_num中
nsel,s,loc,y,maxy
*get,x_direct_node,node,,count    !返回材料2，3的X方向节点数
!定义存储每层单元数目的2维数组
*dim,layer_elem_num,array,x_direct_node-1,layer_no-1
*do,i,1,layer_no-1
jjt=layer_start_node(i)
nsel,s,node,,jjt
elem_no=enextn(jjt,2)
esel,s,elem,,elem_no
*do,j,1,x_direct_node-1
layer_elem_num(j,i)=elem_no
sec_node=nelem(elem_no,2)
nsel,s,node,,sec_node
elem_no=enextn(sec_node,4)
*enddo
*enddo
  
！以下为沉积过程的温度常分析
/solu
antype,trans                         !瞬态分析
trnopt,full
outres,basic,all
kbc,1
  
tintp,0.005,,,1,0.5,0.2
  
nropt,full          !设定牛顿-拉夫森选项
estif,1.0e-8
/auto,1
/replot
  
dt=1e-3            !用于建立初始条件的一个很小的时间段
  
t1=1               !沉积一层单元所需时间
t2=9
  
t=0                 !起始时间
  
t=t+1
!!!!!!!!!!!!!!!!!!!!!!!!!!!该载荷步为基板加热阶段
time,t
deltim,0.1
timint,1,therm
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!下面杀死材料2，3的单元
esel,all
eplot
esel,s,mat,,mat2
esel,a,mat,,mat3
ekill,all
esel,s,live
eplot
nsle,s
nsel,inve
d,all,all,0
  
allsel,all
ic,all,temp,25      !初始温度
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!在模型底部施加温度约束
nsel,s,loc,y,miny
d,all,temp,900            !施加温度约束
  
!基板外侧对流边界条件
nsel,s,loc,x,maxx
nsel,r,loc,y,miny,ny(layer_start_node(1))
sf,all,conv,dep_conv,dep_temp         
!基板上侧对流边界条件
esel,s,live
nsel,s,loc,y,ny(layer_start_node(1))
sf,all,conv,dep_conv,dep_temp           
  
esel,all
nsel,all
  
solve  
esel,s,live
eplot
  
!!!!!!!!!!!!!!!!!!!!!!!!!!去除基板上侧对流边界
esel,s,live
nsel,s,loc,y,ny(layer_start_node(1))
sfdele,all,conv
  
max_tem_mat2=1600         !材料2的初始温度
max_tem_mat3=2680         !材料3的初始温度
  
*do,i,1,layer_no-1
  
!下面的求解用于建立温度的初始条件
t=t+dt
time,t
!!!!!!!!!!!!!!!!!!!!!!!!!deltim,0.001
autots,on
timint,1,therm
  
*do,j,1,x_direct_node-1
jjtelem=layer_elem_num(j,i)
esel,s,elem,,jjtelem
ealive,jjtelem
esel,s,live
eplot
!判断该单元的材料号
*get,mat_ptr,elem,jjtelem,attr,mat
*if,mat_ptr,eq,mat2,then
max_tem=max_tem_mat2
*else
max_tem=max_tem_mat3
*endif
!给新激活单元节点施加初始节点温度
*do,xyz,1,4
  
nj=nelem(jjtelem,xyz)
nsel,s,node,,nj
ddele,nj,all
d,nj,temp,max_tem
  
*enddo
  
*enddo
  
nsel,all
esel,all
  
solve
  
esel,s,live
eplot
  
!下面的步骤用于求解新沉积层对温度场的影响
!计算新沉积层在一秒内温度变化情况，因为在
!此时间内温度变化非常快，所以设置较多子步
t=t+t1-dt
time,t
deltim,0.1,0.01,0.2
autots,on
timint,1,therm
  
*do,j,1,x_direct_node-1
  
jjtelem=layer_elem_num(j,i)
  
*do,xyz,1,4
  
!去除节点温度约束
nj=nelem(jjtelem,xyz)  
nsel,s,node,,nj   
ddele,nj,temp
                       
*enddo
  
!施加上侧对流边界
esel,s,elem,,jjtelem
jjt3=nelem(jjtelem,3)
nsel,s,node,,jjt3
jjt4=nelem(jjtelem,4)
nsel,a,node,,jjt4
sf,all,conv,dep_conv,dep_temp                  
*enddo
!施加右侧对流边界条件
jjtelem=layer_elem_num(x_direct_node-1,i)
esel,s,elem,,jjtelem
jjt2=nelem(jjtelem,2)
nsel,s,node,,jjt2
jjt3=nelem(jjtelem,3)
nsel,a,node,,jjt3
sf,all,conv,dep_conv,dep_temp
  
nsel,all
esel,all
  
!求解开始
solve
计算该沉积层在接下来时间内的温度变化规律
!t=t+t2-dt
!deltim,1,0.5,2
!autots,on
!timint,1,therm
!solve
  
!去除上侧对流边界条件
*do,j,1,x_direct_node-1
  
jjtelem=layer_elem_num(j,i)
esel,s,elem,,jjtelem
jjt3=nelem(jjtelem,3)
nsel,s,node,,jjt3
jjt4=nelem(jjtelem,4)
nsel,a,node,,jjt4
sfdele,all,conv   
  
*enddo
esel,s,live
jjt4=nelem(jjtelem,4)
nsel,s,node,,jjt4
sfdele,all,conv  
  
esel,s,live
eplot
  
*enddo
  
！以下为应力场的分析，就是在这一步要
！！出现上面提到的错误
  
/solu
  
!去除所有边界对流条件
lsclear,fe
  
!去除模型底部施加温度约束
nsel,s,loc,y,miny
ddele,all,temp,900            
finish
  
/PREP7   
ETCHG,TTS   
!*   
KEYOPT,1,1,0
KEYOPT,1,2,0
KEYOPT,1,3,1
KEYOPT,1,5,0
KEYOPT,1,6,0
!*   
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,400
MPTEMP,3,800
MPDE,EX,1   
MPDE,EY,1   
MPDE,EZ,1   
MPDE,NUXY,1  
MPDE,NUYZ,1  
MPDE,NUXZ,1  
MPDE,PRXY,1  
MPDE,PRYZ,1  
MPDE,PRXZ,1  
MPDE,GXY,1   
MPDE,GYZ,1   
MPDE,GXZ,1   
MPDATA,EX,1,,200e9   
MPDATA,EX,1,,179e9   
MPDATA,EX,1,,149e9   
MPDATA,PRXY,1,,  
MPDATA,PRXY,1,,  
MPDATA,PRXY,1,,  
TB,BISO,1,5,2,   
TBTEMP,25   
TBDATA,,1185e6,5e9,,,,   
TBTEMP,540   
TBDATA,,1065e6,5e9,,,,   
TBTEMP,650   
TBDATA,,1020e6,5e9,,,,   
TBTEMP,760   
TBDATA,,740e6,5e9,,,,   
TBTEMP,860   
TBDATA,,330e6,5e9,,,,   
MPTEMP,,,,,,,,   
MPTEMP,1,25  
MPTEMP,2,400
MPTEMP,3,800
MPTEMP,4,1200   
MPDE,EX,2   
MPDE,EY,2   
MPDE,EZ,2   
MPDE,NUXY,2  
MPDE,NUYZ,2  
MPDE,NUXZ,2  
MPDE,PRXY,2  
MPDE,PRYZ,2  
MPDE,PRXZ,2  
MPDE,GXY,2   
MPDE,GYZ,2   
MPDE,GXZ,2   
MPDATA,EX,2,,204e9   
MPDATA,EX,2,,183e9   
MPDATA,EX,2,,154e9   
MPDATA,EX,2,,118e9   
MPDATA,PRXY,2,,  
MPDATA,PRXY,2,,  
MPDATA,PRXY,2,,  
MPDATA,PRXY,2,,  
TB,BISO,2,1,2,   
TBTEMP,0
TBDATA,,270e6,5e9,,,,   
MPTEMP,,,,,,,,   
MPTEMP,1,0   
MPDE,EX,3   
MPDE,EY,3   
MPDE,EZ,3   
MPDE,NUXY,3  
MPDE,NUYZ,3  
MPDE,NUXZ,3  
MPDE,PRXY,3  
MPDE,PRYZ,3  
MPDE,PRXZ,3  
MPDE,GXY,3   
MPDE,GYZ,3   
MPDE,GXZ,3   
MPDATA,EX,3,,150e9   
MPDATA,PRXY,3,,0.25  
TB,BISO,3,1,2,   
TBTEMP,0
TBDATA,,30e6,1e6,,,,
  
finish
  
/solu
t=0
antype,0
nlgeom,on
timint,0,struct
nropt,full
outres,basic,all
t=t+1
time,t
deltim,0.1,0.01,0.5
autots,on
kbc,1
!约束基板Y向位移
nsel,s,loc,y,miny
d,all,uy,0
  
allsel,all
!读取温度分析最后一个载荷步的温度场
ldread,temp,,,13,,'ther-stru','rth'
esel,all
nsel,all
solve
esel,s,live
eplot
  
!finish

ainite

兄弟，先看一下网格，这样的网格能用来计算吗？

Mrice

网格存在什么问题呀，你指的是网格太粗还是单元形状不好亚
我曾经在此基础上将网格加粗了2倍，但是还是出现这样的错误
  
请执教

Mrice

不好意思，打错了
我曾经在此基础上将网格细化了两倍，但还是出现这个问题
  
请执教