在按照上述方法对复合材料进行数值模拟,命令流如下
/PREP7 !进入前处理器
MU0=0.25
tang0=100
normal0=100
dis=0.001
FKN0=1E6
max_step=4
fibre_damage=0
interface_damage=0
ET,1,42!平面单元
KEYOPT,1,3,0 !设置平面应力分析
KEYOPT,1,5,2 !设置输出结点应力
MP,EX,1,400e3 !纤维材料参数
MP,PRXY,1,0.17
MP,EX,2,115e3 !基体材料参数,进入塑性后用21 个应力-应变点拟合
MP,PRXY,2,0.3
TB, MISO,2,1,21,0
TBPT,,0.008260869565,950
TBPT,,0.01030599278,966.25
TBPT,,0.01148279762,982.5
TBPT,,0.01319060099,998.75
TBPT,,0.01568580449,1015
TBPT,,0.01934123137,1031.25
TBPT,,0.02469575862,1047.5
TBPT,,0.03252378039,1063.75
TBPT,,0.04393191531,1080
TBPT,,0.06049295633,1096.25
TBPT,,0.08443048782,1112.5
TBPT,,0.1188721142,1128.75
TBPT,,0.1681951823,1145
TBPT,,0.2384966466,1161.25
TBPT,,0.3382288502,1177.5
TBPT,,0.4790561384,1193.75
TBPT,,0.677004215,1210
TBPT,,0.9539960686,1226.25
TBPT,,1.339896413,1242.5
TBPT,,1.875222619,1258.75
TBPT,,2.614726016,1275
MP, MU,3,MU0 !界面摩擦系数0.25
RECTNG,8.875/100,8.875,0,0.071 ! 创建模型
RECTNG,0,8.875,0.071,0.284
RECTNG,0,8.875,0.284,0.426
RECTNG,0,8.875,0.426,0.639
RECTNG,0,8.875,0.639,0.781
RECTNG,0,8.875,0.781,0.8875
ASEL,S,,,1,5,2 !分配材料属性
AATT,1,1,1
ASEL,S,,,2,6,2
AATT,2,1,1
ASEL,ALL
LSEL,S,,,1,3,2
LESIZE,all,,,99
LSEL,S,,,5,23,2
LESIZE, ALL,,,100
LSEL,S,,,2,4,2
LESIZE, ALL,,,2
LSEL,S,,,6,8,2
LESIZE, ALL,,,6
LSEL,S,,,14,16,2
LESIZE, ALL,,,6
LSEL,S,,,10,12,2
LESIZE, ALL,,,4
LSEL,S,,,18,20,2
LESIZE, ALL,,,4
LSEL,S,,,22,24,2
LESIZE, ALL,,,3
LSEL,ALL
AMESH,1,5,2 !网格化模型
AMESH,2,6,2
*do,i,1,99 !耦合界面重合位置结点
cp,i,ux,2726+i,1011-i
*enddo
*do,i,100,198
cp,i,uy,2627+i,1110-i
*enddo
*do,i,199,297
cp,i,ux,609+i,2423-i
*enddo
*do,i,298,396
cp,i,uy,510+i,2522-i
*enddo
*do,i,397,495
cp,i,ux,902-i,1623+i
*enddo
*do,i,496,594
cp,i,uy,1001-i,1524+i
*enddo
*do,i,595,693
cp,i,ux,-292+i,2112-i
*enddo
*do,i,694,792
cp,i,uy,-391+i,2211-i
*enddo
*do,i,793,890
cp,i,ux,521+i,994-i
*enddo
*do,i,891,988
cp,i,uy,423+i,1092-i
*enddo
cp,next,all,2725,911
cp,next,all,2726,907
cp,next,all,806,2125
cp,next,all,807,2119
cp,next,all,2018,406
cp,next,all,2019,402
cp,next,all,301,1418
cp,next,all,302,1412
cp,next,all,1313,103
cp,next,all,1312,101
!下面是创建5 对界面接触单元
ET,2,169 !定义目标单元
ET,3,172 !定义接触单元
R,3,,,FKN0,0.1,0, !定义接触刚度为1E6
MAT,3
REAL,3
LSEL,S,,,3
TYPE,2
NSLL,S,1
ESLN,S,0
ESURF !生成目标单元,下面类似
LSEL,S,,,5
TYPE,3
NSLL,S,1
ESLN,S,0
ESURF !生成接触单元,下面类似
ET,4,169 !定义目标单元
ET,5,172 !定义接触单元
R,4,,,FKN0,0.1,0, !定义接触刚度为1E6
MAT,3
REAL,4
LSEL,S,,,9
TYPE,4
NSLL,S,1
ESLN,S,0
ESURF
LSEL,S,,,7
TYPE,5
NSLL,S,1
ESLN,S,0
ESURF
ET,6,169 !定义目标单元
ET,7,172 !定义接触单元
R,5,,,FKN0,0.1,0, !定义接触刚度为1E6
MAT,3
REAL,5
LSEL,S,,,11
TYPE,6
NSLL,S,1
ESLN,S,0
ESURF
LSEL,S,,,13
TYPE,7
NSLL,S,1
ESLN,S,0
ESURF
ET,8,169 !定义目标单元
ET,9,172 !定义接触单元
R,6,,,1FKN0,0.1,0, !定义接触刚度为1E6
MAT,3
REAL,6
LSEL,S,,,17
TYPE,8
NSLL,S,1
ESLN,S,0
ESURF
LSEL,S,,,15
TYPE,9
NSLL,S,1
ESLN,S,0
ESURF
ET,10,169 !定义目标单元
ET,11,172 !定义接触单元
R,7,,,FKN0,0.1,0, !定义接触刚度为1E6
MAT,3
REAL,7
LSEL,S,,,19
TYPE,10
NSLL,S,1
ESLN,S,0
ESURF
LSEL,S,,,21
TYPE,11
NSLL,S,1
ESLN,S,0
ESURF
LSEL,ALL
ALLSEL,ALL
*DIM, FS,,998,1
*VREAD,FS(1,1),G:\cell\fiberstrength,csv,,jik,1,998
(1f6.1)
/SOLU
*DIM,comp,array,494 !定义界面结点等效应力数组
*DIM,interface,array,494 !
*DIM,normal,array,494 !定义界面结点法向应力数组
*DIM,tang,array,494 !定义界面结点剪切应力数组
*DIM,fes,array,998 !定义纤维单元拉伸应力数组
Rescontrol,define,all,1,1 !将每个载荷步和子步结果读出到文件parameter.txt 中
! 下面是施加第一次求解的位移载荷
/SOLU !进入中间处理器
LSEL,S,,,1
LSEL,a,,,8
LSEL,a,,,12
LSEL,a,,,16
LSEL,a,,,20
LSEL,a,,,24
DL,ALL, ,SYMM
LSEL,ALL
dl,2,1,ux,dis/20
dl,6,2,ux,dis/20
dl,10,3,ux,dis/20
dl,14,4,ux,dis/20
dl,18,5,ux,dis/20
dl,22,6,ux,dis/20
nlgeom,on !设置大变形分析
nsubst,10 !设置载荷子步数
cnvtol,f,,0.1,,0.1 !设置力收敛准则
ALLS
solve
/POST1
save,job%first%,db
*do,i,2,max_step!进入重启动循环求解
/POST1 !进入后处理器
prnsol !列出所有结点的所有方向的应力结果
presol !列出所有单元的所有方向的应力结果
etable,sx,s,x !获取所有单元x 方向应力解
*do,j,1,998
*get,fes(j),elem,j,etab,sx !获取每个纤维单元的拉伸应力值
*enddo
*do,j,1,99
*get,normal(j),node,1011-j,s,y !获取界面结点法向应力值,下面类似
*get,tang(j),node,1011-j,s,xy !获取界面结点剪切应力值,下面类似
comp(j)=(normal(j)/normal0)*(normal(j)/normal0)+(tang(j)/tang0)*(tang(j)/tang0)
!获取判定界面失效的结点应力值,下面类似
*enddo
*do,j,100,198
*get,normal(j),node,708+j,s,y
*get,tang(j),node,708+j,s,xy
comp(j)=(normal(j)/normal0)*(normal(j)/normal0)+(tang(j)/tang0)*(tang(j)/tang0)
*enddo
*do,j,199,297
*get,normal(j),node,704-j,s,y
*get,tang(j),node,704-j,s,xy
comp(j)=(normal(j)/normal0)*(normal(j)/normal0)+(tang(j)/tang0)*(tang(j)/tang0)
*enddo
*do,j,298,396
*get,normal(j),node,5+j,s,y
*get,tang(j),node,5+j,s,xy
comp(j)=(normal(j)/normal0)*(normal(j)/normal0)+(tang(j)/tang0)*(tang(j)/tang0)
*enddo
*do,j,397,494
*get,normal(j),node,598-j,s,y
*get,tang(j),node,598-j,s,xy
comp(j)=(normal(j)/normal0)*(normal(j)/normal0)+(tang(j)/tang0)*(tang(j)/tang0)
*enddo
parsav,all,parameter,txt !保存所有参数到文件parameter.txt 中
/SOLU
antype,,rest !激发重启动分析
parres,new,parameter,txt !从文件parameter.txt 中读出所有参数数据
*do,j,1,998
*if,fes(j),gt,fs(j),then
ekill,j
fibre_damage=fibre_damage+1
*endif
*enddo !判断纤维单元是否失效
*do,j,1,99
*if,comp(j),gt,1.0,then
CPDELE,j,j,any
CPDELE,j+99,j+99,any
interface_damage=interface_damage+1
interface(j)=1
*endif
*enddo !
基于界面应力失效准则,判断结点是否分离,下面类似
*do,j,100,198
*if,comp(j),gt,1.0,then
cpdele,j+99,j+99,any
cpdele,j+198,j+198,any
interface_damage=interface_damage+1
interface(j)=1
*endif
*enddo
*do,j,199,297
*if,comp(j),gt,1.0,then
cpdele,j+198,j+198,any
cpdele,j+297,j+297,any
interface_damage=interface_damage+1
interface(j)=1
*endif
*enddo
*do,j,298,396
*if,comp(j),gt,1.0,then
cpdele,j+297,j+297,any
cpdele,j+396,j+396,any
interface_damage=interface_damage+1
interface(j)=1
*endif
*enddo
*do,j,397,494
*if,comp(j),gt,1.0,then
cpdele,j+396,j+396,any
cpdele,j+494,j+494,any
interface_damage=interface_damage+1
interface(j)=1
*endif
*enddo
dl,2,1,ux,i*0.1/20
dl,6,2,ux,i*0.1/20
dl,10,3,ux,i*0.1/20
dl,14,4,ux,i*0.1/20
dl,18,5,ux,i*0.1/20
dl,22,6,ux,i*0.1/20
nlgeom,on !设置大变形分析
autos on !设置自动步长
cnvtol,f,,0.1,,0.1 !设置力收敛准则
!cnvtol,f,,1,,1 !设置力收敛准则
NROPT, FULL, , OFF
ncnv,0 !遇到求解不收敛时不退出
solve !求解
tiwebull=i-1
save,job%tiwebull%,db
*enddo
以上程序将纤维和基体界面处的位移进行耦合,表示粘接良好状态,两者只将通过接触单元描述界面脱粘后两者之间的应力传递,脱粘后,脱粘处的耦合删除。
然而,在执行以上程序时,载荷导致界面失效后,一旦相应的耦合被删除后,在以后计算中出现以下错误,导致程序执行不下去。
出现的错误提示如下:
One or more elements have become highly distorted. Excessive
distortion of elements is usually a symptom indicating the need for
corrective action elsewhere. Try incrementing the load more slowly
(increase the number of substeps or decrease the time step size). You
may need to improve your mesh to obtain elements with better aspect
ratios. Also consider the behavior of materials, contact pairs,
and/or constraint equations. If this message appears in the first
iteration of first substep, be sure to perform element shape checking.
实际上在最开始加载过程中,若纤维和基体界面上不加任何耦合,光有接触单元时,在第一步计算就不收敛,除非将接触刚度系数FKN设置为10以下,但此时计算结果明显不对。
也就是若纤维和基体界面上如果有节点没耦合上,则程序就运行不下去。
出现上述问题的原因是什么呢? |