热应力－结构耦合分析问题

wliustc

[请教]热应力－结构耦合分析问题

小弟为ansys新手，在学习时遇到问题。
这是别人文章里的一个结果，我在重复实现时，总是无法出现一致结果，请大家帮忙看看问题出在什么地方，谢谢！
激光辐照充内压圆柱壳体（壳有一定厚度）

                               
登录/注册后可看大图

                               
登录/注册后可看大图

物理模型如上图所示
P为内压，F为等效拉应力
我的问题是，热分析的温度分布基本和文献一致，但是相应的热应力以及结构分析的结果，远远弱于文献中的结论。
下面是用间接法分析的命令流（激光采用gui方式加载）
FINISH
/CLEAR
/BATCH
/FILNAME, thermal,1
/TITLE, thermal
!SI单位制下：
!************basic parameters*************************
Ro=0.027        !outer radius(m)
ri=0.023        !inner radius(m)
Ra=0.040        !laser as(m)
d=0.004         !thickness(m)
L=0.080         !length(m)
!****parameters of laser****
rl=0.020        !radius of laser(m)
I=2.313*1.5e6    !average power(W/mm^2) [150W/cm^2=1.50e6W/m^2]
!****parameters of press****
p=5e6         !inner press(Pa)  [5MPa]
p0=31.63e6      !limit press(Pa)
/PREP7
!柱壳材料为Ly12铝合金
ET,1,SOLID70                               !element type
MPTEMP,1,20,100,200,250,300,350                     !temp table(`c) 6个点
!材料密度
MP,DENS,1,2.8e3                             !density(Kg/m^3) [2800Kg/m^3]
!比热    c = 770.045 + 1.671*t - 0.000152*t^2
MPDATA,c,1,1,0.900e3,0.921e3,1.047e3,1.092e3,1.130e3,1.172e3      !specific heat(J/(Kg.`c)) [900J/(Kg.`c)]
!热传导系数
MPDATA,KXX,1,1,0.121e3,0.121e3,0.121e3,0.121e3,0.121e3,0.121e3      !thermal conductivity(W/(mm.`c)) [121W/m.c]
TOFFST,273    !参考温度，初始温度
!************mode mesh*************************
K,200,0,0,0
wprot,0,0,90
PCIRC,Ro,ri,0,180,
VOFFST,1,-L, ,
wprot,0,0,-90
VSBW,    1
PCIRC,RL, ,90,270,
VOFFST,1,RA, ,
VSBA,       2,       3
VDELE,      1, , ,1
K,100,-RA,0,0
KWPAVE,     100
wprot,0,0,90
VSBW,       5
VSBW,       3
KWPAVE,     200
wprot,0,0,-90
SAVE
!划分单元
Lesize,2,,,4     !厚度划分
Lesize,4,,,4
Lesize,6,,,4
Lesize,8,,,4
Lesize,14,,,4
Lesize,15,,,4
Lesize,38,,,4
LESIZE,30, , ,20,4, , , ,1   !光斑区不等分划分
LESIZE,31, , ,20,4, , , ,1
LESIZE,33, , ,20,4, , , ,1
LESIZE,34, , ,20,4, , , ,1
Lesize,16,,,10
Lesize,17,,,10
Lesize,21,,,10
Lesize,22,,,10
Lesize,23,,,10
Lesize,24,,,10
MSHAPE,0,3d
MSHKEY,1
VMESH,4   
MSHKEY,0
VSWEEP,2,7,23         !过渡区不等分划分

VSWEEP,1            !其他部分
VSWEEP,5
VSWEEP,6
ALLSEL,ALL
SAVE
FINISH
/SOL
!*  
!分析设置
ANTYPE,4
TRNOPT,FULL
LUMPM,0   
NLGEOM,1
OUTPR,BASIC,ALL,
OUTRES,ALL,ALL,
!加载载荷
!初始温度
TUNIF,30,
!加载激光热流密度
*DEL,_FNCNAME   
*DEL,_FNCMTID   
*DEL,_FNC_C1
*DEL,_FNC_C2
*DEL,_FNCCSYS   
*SET,_FNCNAME,'myload'  
*DIM,_FNC_C1,,1
*DIM,_FNC_C2,,1
*SET,_FNC_C1(1),I   
*SET,_FNC_C2(1),RL  
*SET,_FNCCSYS,0
! /INPUT,guass.func,,,1
*DIM,%_FNCNAME%,TABLE,6,17,1,,,,%_FNCCSYS%  
!   
! Begin of equation: I*exp(-2*({X}^2+{Y}^2)/RL^2)   
*SET,%_FNCNAME%(0,0,1), 0.0, -999   
*SET,%_FNCNAME%(2,0,1), 0.0
*SET,%_FNCNAME%(3,0,1), %_FNC_C1(1)%
*SET,%_FNCNAME%(4,0,1), %_FNC_C2(1)%
*SET,%_FNCNAME%(5,0,1), 0.0
*SET,%_FNCNAME%(6,0,1), 0.0
*SET,%_FNCNAME%(0,1,1), 1.0, -1, 0, 0, 0, 0, 0  
*SET,%_FNCNAME%(0,2,1), 0.0, -2, 0, 1, 0, 0, -1
*SET,%_FNCNAME%(0,3,1),   0, -3, 0, 1, -1, 2, -2
*SET,%_FNCNAME%(0,4,1), 0.0, -1, 0, 2, 0, 0, -3
*SET,%_FNCNAME%(0,5,1), 0.0, -2, 0, 1, -3, 3, -1
*SET,%_FNCNAME%(0,6,1), 0.0, -1, 0, 2, 0, 0, 2  
*SET,%_FNCNAME%(0,7,1), 0.0, -3, 0, 1, 2, 17, -1
*SET,%_FNCNAME%(0,8,1), 0.0, -1, 0, 2, 0, 0, 3  
*SET,%_FNCNAME%(0,9,1), 0.0, -4, 0, 1, 3, 17, -1
*SET,%_FNCNAME%(0,10,1), 0.0, -1, 0, 1, -3, 1, -4   
*SET,%_FNCNAME%(0,11,1), 0.0, -3, 0, 1, -2, 3, -1   
*SET,%_FNCNAME%(0,12,1), 0.0, -1, 0, 2, 0, 0, 18
*SET,%_FNCNAME%(0,13,1), 0.0, -2, 0, 1, 18, 17, -1  
*SET,%_FNCNAME%(0,14,1), 0.0, -1, 0, 1, -3, 4, -2   
*SET,%_FNCNAME%(0,15,1), 0.0, -1, 7, 1, -1, 0, 0
*SET,%_FNCNAME%(0,16,1), 0.0, -2, 0, 1, 17, 3, -1   
*SET,%_FNCNAME%(0,17,1), 0.0, 99, 0, 1, -2, 0, 0
! End of equation: I*exp(-2*({X}^2+{Y}^2)/RL^2)
!-->
  
SFA,17,1,HFLUX, %MYLOAD%
ALLSEL,ALL
SBCTRAN
TRNOPT,FULL !设置计算方法
*DO,DELTA,0.001,3.5,0.1 !设置循环计算
TIME,DELTA !设置荷载步时间
AUTOTS,0 !关闭自动时间步长
KBC,0 !设置阶梯荷载步
TSRES,ERASE !Defines an array of keytimes at which the time-stepping strategy changes
OUTRES,ALL,ALL !输出所有项
SOLVE !开始计算
*ENDDO !终止循环计算
FINISH
!结构应力分析
/PREP7 !前处理
ETCHG,TTS !转换热分析单元为结构单元
!输入参数
!热胀系数  alpha = 21.079 + 0.0175*t - 0.000018*t^2
MPDATA,ALPX,1,1,21.4e-6,22.7e-6,23.8e-6,24.3e-6,24.71e-6,25.00e-6    !linear expansion coefficient(1/`c)e-6
MPPLOT,KXX,1
!弹性模量  ex = 69.3 - 0.048*t - 0.000144*t^2
MPDATA,EX,1,1,68E9,64E9,54E9,48E9,42E9,34.9E9              !elastic modulus(Pa) [58E9 N/m^2]{}
!泊松比
MPDATA,PRXY,1,1,0.3,0.3,0.3,0.3,0.3,0.3,0.3                !Poisson ratio
MPPLOT,EX,1  
!屈服强度σs和抗拉强度σb和硬化模量
!MPTEMP,1,      47,   87,   143,    197,   243,   289,   343,   403
!yield stress(Pa)   4.11e8,3.983e8,3.794e8,3.303e8,2.742e8,1.613e8,0.811e8,0.672e8
!tension strength(Pa)6.38e8,5.802e8,5.488e8,5.155e8,3.989e8,1.782e8,1.025e8,0.676e8         
!tang mod(Pa)     6.72e8,6.403e8,5.949e8,5.426e8,4.913e8,4.340e8,3.589e8,2.657e8     
!硬化模量EX/10(<300);EX/50(>300)
!假设各材料都是双线性随动硬化弹塑性本构关系
!下面是各材料在各温度点下的屈服应力和屈服后的弹性模量
!Yield stress(Pa)  [411MPa]
!Tang Mod(Pa)
TB,BISO,1,6,2,  
TBTEMP,47   
TBDATA,1,4.110e8,6.720e9
TBTEMP,87  
TBDATA,1,3.983e8,6.403e9
TBTEMP,197  
TBDATA,1,3.303e8,5.426e9
TBTEMP,243
TBDATA,1,2.742e8,4.913e9
TBTEMP,289
TBDATA,1,1.613e8,4.340e9
TBTEMP,343  
TBDATA,1,0.811e8,3.589e9
ALLSEL !选择选择所有
/SOLU !进入求解器
NROPT,FULL !Specifies the Newton-Raphson options in a full transient analysis
NLGEOM,ON !考虑大变形的影响
!对称位移约束
NSEL,S,LOC,X,0  
DSYM,SYMM,X,,
NSEL,S,LOC,Y,0  
DSYM,SYMM,Y,,
!对称位移约束
!柱壳另一端固定？
NSEL,S,LOC,X,-L
D,ALL, ,0, , , ,UX,UY,UZ, , ,
NSEL,S,LOC,X,-L
F,ALL,FX,-P*RI/D

!加载内压  
SFA,3,1,PRES,P
SFA,16,1,PRES,P
SFA,18,1,PRES,P
SFA,24,1,PRES,P
SFA,28,1,PRES,P  
ALLSEL,ALL
SBCTRAN
*DO,DELTA,0.001,3.5,0.1
LDREAD,TEMP,,,DELTA, ,'thermal1','rth',' ' !TEMP LOAD Reads results from the results file and applies them as loads.
TIME,DELTA
AUTOTS,0
KBC,0 !设置阶梯荷载步
TSRES,ERASE !Defines an array of keytimes at which the time-stepping strategy changes
OUTRES,ALL,ALL, !输出所有的结果
SOLVE !计算
*ENDDO
FINISH

                               
登录/注册后可看大图

phymodel.JPG (15.62 KB)
2007-7-31 09:38

mututu_gui

发一下别人文献中的结果呀，不然我们做了，也不知道有无此类问题呀。

venture

!对称位移约束
NSEL,S,LOC,X,0  
DSYM,SYMM,X,,
NSEL,S,LOC,Y,0  
DSYM,SYMM,Y,,
!对称位移约束
!柱壳另一端固定？
NSEL,S,LOC,X,-L
D,ALL, ,0, , , ,UX,UY,UZ, , ,
NSEL,S,LOC,X,-L
F,ALL,FX,-P*RI/D

其他约束看起来没有什么问题，下面的有点问题
NSEL,S,LOC,X,-L
D,ALL, ,0, , , ,UX,UY,UZ, , ,
NSEL,S,LOC,X,-L
F,ALL,FX,-P*RI/D
x=-l的面全约束了，你再加等效拉力有什么用？等效拉力的值好像也不对。改成下面的试试。
NSEL,S,LOC,X,0
F,ALL,FX,-P*RI/(2*D)
KSEL,,LOC,X,-L
KSEL,R,LOC,Y,RI+D
DK,all, , , ,0,UZ, , , , , ,  
要么建1/2模型，筒体上端全约束，下端加等效拉力。

[ 本帖最后由 venture 于 2007-8-15 08:50 编辑 ]

venture

我的问题是，热分析的温度分布基本和文献一致，但是相应的热应力以及结构分析的结果，远远弱于文献中的结论。
那肯定是结构分析出问题了，在结构分析里面找原因。

飘叶和流星

大家好
我用的是间接法模拟焊接温度场跟残余应力场的
我有两个问题亟待解决
（1）我看到好多文献一般都用solid70进行热分析
结构分析用的是solid45单元
但我计算完温度再算应力进行单元转换时solid70转换成了solid185
不知道为什么
（2）我看到一些文献进行应力求解时候，这样做的：施加载荷时，读入热分析的节点温度并制定相应的时间点或载荷步即可。（这个没问题，都会）因为焊接过程中熔池区的金属处于融化状态，即进入零力学性能状态，其所有的应力应变将消失。对此，采用生死单元的方法，在每一步热应力计算时，将对应的温度场计算结果进行选择，超过熔点的单元令其“死掉”，而低于熔点的单元令其“激活”。我不知道这种做法对不对。如果对的话该怎么把节点温度为熔点以上的单元选出来呢，能不能指点一下？谢谢
   补充一点，因为我做的是薄板焊接，所以我只是用高斯移动热源进行加载热源，并没有采用生死单元技术。
  期待大家的帮助，再一次感谢！

lilichunyan

请问，在热应力分析的时候，应力分析部分的程序一定要加位移约束吗？不加可以吗？

jxxansys

我这里有一个热焊接模拟的实例，可供参考。
需要的话，请联系！