2017年5月24日 星期三

查阅相关文献中关于校验刻划实验和模拟结果的方法

1、非应力的比对校验

A model for stresses, crack generation and fracture toughness
calculation in scratched TiN-coated steel surfaces

用abaqus模拟,能模拟出镀层的应力场但是没有裂纹。将拉应力极值出现位置和实验结果比对。

3. With increasing load a new stress-field with a tetra-armed
star shape grows around the contact area. At the tail arms,
stress concentrations are amplified at a distance of 1–2
times the contact length from the edge of the contact at
the border of the scratch groove. The magnitude of these
stress concentrations is of the same level as within the
contact area.
4. The stress concentrations at the tail arms travel with increasing
load to the plane of symmetry and increase in
magnitude considerably faster than within the contact
area, and become the dominating tensile stress as the
contact conditions change to ploughing mode.
5. After about 1mm of sliding in the scratch test, a peak
area of maximum first principal stress is generated in
the back-tail region at the border of the scratch groove,
creating the first visible angular cracks in the coating with
this material-coating combination.

参考:abaqus、应力场分析、与实验校验的模式

2、心不死的JH2论文检索

A numerical model for optical glass cutting based
on SPH method

K9玻璃,JH2本构,2d模型,SPH方法(要继续看???)特点是虽然能看到裂纹产生但是分子会乱飞,和有限元完全不同,类似于md的效果。应力分布很乱。

如果要选这个本构和sph,以后可以参考这一篇进行复现。

K9的JH2本构可以追溯到以下两篇文章:

(甚至于tensile strength也在其中)

Implementation and Validation of the Johnson-Holmquist Ceramic Material Model in LS-Dyna

Silica Float Glass [11]

最早的出处,也是此文中参数的来源,引文11:

1995 High Strain Rate Properties and Constitutive Modeling of Glass

作者正是J和H两人,文中所写的材料是浮法玻璃。

文中明确写道,改材料的tensile strength为0.015Gpa,15Mpa;而强度则为0.30Gpa,300Mpa,这个参数和顾等人的论文中的3.5Gpa相去甚远。

3、在球形压头印压和刻划(后者更重要)方面的深水水域

On the Appropriate Use of Representative Stress Quantities
at Correlation of Spherical Contact Problems

ON THE NUMERICS AND CORRELATION OF SCRATCH TESTING

以上3篇全部讲了弹塑性材料刻划模拟中的一个参数,由Johnson提出的(j-c、j-h、and here again)一个参数

TIM图片20170524232410

并且分区讨论了这个参数对于应力场分布、弹塑性行为的影响。

讨论的用意在于进行“CORRELATION”。

4、明日的任务是继续找相关的实验、方法验证模拟模型的正确性

2017年5月22日 星期一

abaqus 应力场分析

1、统一的图例/标尺

Screenshot from 2017-05-22 14-13-43.png

如图所是右侧一个图例是之前设置好后截图的效果,和中间设置效果对比可知,只要每次在图例的设置中设置好等值线的数量:BASIC-CONTOUR INTERVALS-DISCRETE-24;然后在LIMITS上下限的选项单改上下限为2000、-3000,全部勾选上show location,apply后的效果和右侧截图相同。

注意,超出上限的部分全部以灰色显示,低于下限的全部以黑色显示,因此必须对所有要加入对比的结果进行预先的调查,找到一个适合的上下限,再进行设置。

图例的统一是至关重要的,他将不同模型的不同时刻的应力等值线的值进行了统一,同一种颜色的等值线就是对应唯一的应力数值,该等值线所包围的区域越大,说明应力越大。

2、寻找合适的截面获得信息

首先,将刻划的位移量和刻划深度以整数对应起来是非常便利的设置方法。比如在以下模型中,U3位移为0.01mm的时候,刻划深度U2为0.0005mm,因此0.008mm位移时,深度为0.0004mm。

注意,选取z=0.01的截面进行分析,可以获得设计实验时想获得的“500nm深刻划应力场”的数据,但是相比而言,z=0.008截面经历过刻划前、中和后3个阶段的状态,更能反映刻划过程的应力场变化。因此,这两个面上获取的信息都非常有用。另外,还有Y0Z平面也需要进行仔细分析。

注意2:关于不同刻划深度时应力场分布的大小,实际上是有很多不同之处的。考虑2um的压头从浅到深刻划了500nm;再考虑一个压头的钝圆半径,一边刻划,钝圆半径一边从0增大到500nm,最终刻划深度当然也是500nm。这两种情况下的应力场分布极其不同,后者再深度增加的事后,钝圆半径也在增大,其应力场近似于等比缩放(小球压入深度浅,大球压入深度大,但只要保持压入深度和半径的比率相同,应力场的分布也只是近似放大而已,没有空间分布上的大差异),而前者则是从一小个球面的压入一直过渡到一个半球,压头接触部分的形状不停地发生改变。(如果无法理解这段文字可以去死了)因此研究不同下压深度的应力场分布差异,也是有很大意义的。

3、以动态的方式迅速读懂结果

选好截面,设置好统一标尺以后,可以进行动态的结果输出,输出为avi后可以随时对视频截图。同时,动态的结果更能展示变化的趋势。

5um_0008surf

很明显地,右侧先刻划产生的应力场因受到第二次刻划时候的应力场影响,应力大小都同时降低了。(即等值线包围的区域缩小了)

4、简单分析结果

  • 先后和同时刻划

比如下图为先后刻划0.008截面的应力场:

Screenshot from 2017-05-22 15-20-21Screenshot from 2017-05-22 15-23-13

不管是4um的还是5um的顺序刻划,其中央部分的应力都是+500以下。但是将5um间距的同时刻划和5um的先后刻划结果进行比对,结果却有很大不同:

Screenshot from 2017-05-22 15-20-21Screenshot from 2017-05-22 15-26-17

无法描述应力的分布是否有很大不同,但是至少拉应力都普遍更大。

2017年5月10日

今天的目标是

  • 试着enriched区域进行xfem模拟,看中位裂纹是否产生和扩展
  • 将hard contact改成softened contact,检查应力状态
  • 建一个mm级别模型和brinell的公式进行校核,排除是本构问题还是尺寸效应/接触问题

1、建立如图的一个区域,

Screenshot from 2017-05-10 23-51-12.png图1

结果是中位裂纹萌生了,但是没有扩展现象,扩展条件有问题?针对这个现象有一个想法,就是VCCT的应用最好让裂纹在实体的边缘开始产生,这样才能判断哪里是裂纹尖端的单元。参考:11.4.3中的VCCT tips

For two-dimensional models with contact pairs involving higher-order underlying elements, the initially unbonded portion must extend over complete element faces. In other words, the crack tip in a two-dimensional, higher-order model must start at a corner node on the quadratic slave surfaces. The crack tip must not start at a midside node.

2、softened contact的机理和应用,参考37.1.2中的softened contact,以下两个图可以清晰地解释其机制。

图2ab

  • 硬接触就是指,两个接触面之间距离到0之前无压力,到达0之后可有任意大的压力
  • 软接触是指,距离靠近到C0时(未到0,实际未接触)可以开始有压力,到实际接触时压力到达P0,并且可以穿过从面,压力继续上升

实际结果就是主面能嵌入从面中,考虑到圆弧面在小位移下接触面积的增加是二次的,这个的影响会有,但相比3个数量级的印压力大小仍有很大差异。模拟结果基本符合这个猜想。印压的支反力仍然达到10N以上。

3、本构的理解

fig_3

图3

  • E弹性模量指单位变形量下的应力值,越大则材料越硬
  • sigmay是屈服强度,越大则越难产生塑性变形

在典型的印压实验中,接触应力会快速地到达sigmay,在开始屈服以后,应力的大小将不会增加得太快。更清晰的描述是,对于单个单元,压缩量超过屈服时的应变后,产生应变所需的应力增量小了,其实应该说,硬化系数就是这之中起作用的关键,因为alpha=0.01,因此塑性变形所需应力增量更小。

但是宏观的表现一定为:接触应力迅速上升至屈服应力。

从这一层面看,我们的模拟结果应力大是没有问题的,比如钢的屈服强度是300Mpa,那mises应力必然会快速到达300,然后再慢慢往上升。

4、针对同学提出,可能是参数单位有问题,或者是尺寸效应的问题,我尝试用brinell硬度推倒的一个力和球型压头下压深度的关系对模拟进行校核。

https://en.wikipedia.org/wiki/Brinell_scale

使用wiki附上的材料参数,glass的HB为1500,根据公式Screenshot from 2017-05-11 00-08-29.png

图4

F=46196*R*sigma

  • sigma为下压深度
  • R为压头
  • F为所需下压力

则旧模型中2um压头压入0.5um所需的压力为46mN,而若用1mm压头压入0.5mm则需要23kN的压力。

建模进行校核结果分别为:

Screenshot from 2017-05-11 00-17-50Screenshot from 2017-05-11 00-18-15

12N和9.9kN,后者比较符合计算结果,前者差3个量级以上,说明该理论可能更适合mm级别,在um级别不太符合,而材料参数没有很大问题。

下一步是校核模拟结果和实际um级别印压实验的力。

2017年5月9日 周二

1、自己从零重建了模型

  • 增加了金刚石的材料属性
  • 把压头改成了shell,划分网格,定义材料

换的原因时因为觉得应力过大是接触问题,把压头当成解析刚体的结果,改成具有弹性的金刚石材料应该能缓解问题。这是接触力学里面的知识,在计算的时候,必须考虑两个接触面的弹性,而两个表面变形以后,接触面积也会发生改变。

  • 接触仍然时surface to surface,hard contact和frictionless
  • 加载方式是位移加载,因为不知道设体力是否不正确,而点加载会出错

2、建模的过程有很多小问题

  • 一开始的parts应该选用2D deformable,shell
  • section里面新建setion用SOLID,HOMOGENEOUS,就只需要设材料,且能对应2Dshell进行材料定义
  • sets可以一开始建,但是后续assembly里面会在set前加上parts的名字

3、结果

  • 第一个job test1是-0.5N体力加载的结果,应力e-6量级,太小
  • 第二个job是位移下压500nm,应力和以前没什么区别,非常大

接下来的任务是改hard contact等contact properties;参考下面的网页,试试排错

 

4、daily read

http://imechanica.org/node/7960

https://en.wikipedia.org/wiki/Brinell_scale

借用里面的公式,可以反求出压力F和下压深度d之间的关系:

压头钝圆半径在2um时(0.002)

F/d=47.7

在下压0.0005mm的深度时,1550HB硬度的玻璃,其下压力应大约是0.023N(23mN)

因此模拟还是有很大问题,这个计算结果和实验结果量级相当,模拟要继续改进。

2017年5月8日 周一

1、研究了一下,力加载条件和位移加载之间的差异:

  • 力加载和位移加载不会对应力场分布没有太大影响
  • 力加载时可以用U2测位移量,位移加载时可以用RF2测丫头的支反力作为下压力
  • 力加载和位移加载的条件,对应关系是:500nm深对应10N左右的力
  • 模拟结果和实验结果相差很大,实验时1um深对应刻画力约为100mN左右

2、本构模型正确的情况下,为了尽快排除错误,最直观的做法是,重做出Weibin Gu的结果。如果要重做出Weibin Gu的模拟结果,需要调试方向可能是:

  • 网格类型?
  • 接触形式(切向行为,法向行为)
  • 压头的形状

力和刻划/印压深度之间对应关系的问题必须及时解决。

3、每日一读Documentation

http://50.16.225.63/v2016/books/usb/default.htm?startat=pt04ch11s04aus70.html

线性缩放(?linear scaling)加速VCCT的收敛

  • *CONTROLS, TYPE=VCCT LINEAR SCALING
  • step module: OTHER——general solution control——edit: stepname, VCCT LINEAR SCALING

VCCT使用小贴士

  • 显式算法时,用general contact可以获得更好的鲁棒性
  • 应尽量减少模型中动能对内能(?internal energy)的比例
    • 尽量避免质量缩放
    • 使用尽量平滑的加载函数
  • 用quad网格类型,控制网格的长宽比不要超过2,以正方形为佳
  • 三种形式的断裂条件不应该差异过大
  • 鉴于裂纹扩展过程发生在一瞬间,输出频率最好更大,能更好发现问题

Tips for using the VCCT criterion in Abaqus/Explicit
Crack propagation problems using the VCCT criterion analyzed in Abaqus/Explicit benefit from the robustness of the general contact algorithm in the context of an explicit time integrator. Nevertheless, as is the case in Abaqus/Standard, these analyses remain challenging given the discontinuous nature of the fracture phenomenon. The following tips will help you create a successful Abaqus/Explicit model:

  • Dynamic effects are of utmost relevance when assessing the results from a debonding analysis using the VCCT criterion. In most cases experimental and/or theoretical data are available in quasi-static settings. You must ensure that the Abaqus/Explicit analysis generates low ratios of kinetic energy to internal energy (1% or less). In practical terms this requirement often translates into avoiding the use of mass scaling in the model. Use smooth amplitudes to drive the loading to help reduce the kinetic energy in the model. Running the analysis over a longer period of time will not help in most cases because bond breakage is an inherently fast and localized process.
  • If appropriate, use damping-like behavior in the materials associated with the debonding plates to reduce dynamic vibrations. Unlike Abaqus/Standard, where a pure static equilibrium is achieved at the end of a converged increment, in Abaqus/Explicit the bond breakage at a given location is associated with a dynamic overshoot beyond the static equilibrium position. If the vibrations are significant (kinetic energy is clearly observable), the dynamic overshoot at nodes behind the crack tip may lead to premature debonding of the crack tip.
  • To maximize the accuracy of the debonding simulation, use quad meshes between the slave and master surfaces of the debonding surfaces. Avoid using elements with aspect ratios greater than 2. In most cases mesh refinement will help with obtaining a realistic result.
  • Highly mismatched critical energy values between modes tend to induce crack propagation in continuously changing directions in a manner that may be unstable and unrealistic, particularly for modes II and III. Do not use such values unless experimental data suggest so.
  • Use frequent field output requests to evaluate the debonding evolution as the analysis progresses. In some cases this can point to nontrivial modeling deficiencies that are difficult to identify from a simple data check analysis.
  • Avoid the use of other constraints involving nodes on both surfaces of the debonding interface because the cohesive contact forces will compete with the constraint forces to achieve global equilibrium. Bond breakage might be hard to interpret in these cases.


明天一读,接触。

http://50.16.225.63/v2016/books/usb/default.htm?startat=pt04ch11s04aus70.html

2017年5月5日 周五

1、2dxfem里面的接触尽量简单就行,用surface2surface,无摩擦最好,另外看了很多interaction properties的设置都没有norminal behavior,考虑删除,但是影响可能并不大。

http://shodhganga.inflibnet.ac.in/bitstream/10603/17557/9/09_chapter4.pdf

这一篇里面对两种材料——钢和铝进行了压印实验,文章中用的网格类型是CAX4R,但是在abq2016只能找到CPlainStress4R,还不知道有什么区别。

将塑性部分的应变修改以后,实验结果和之前差异不是太大,应力场的分布似乎变得更合理了,但是仍旧是压头两边快速开裂。

2、主要的时间和经历用在了池沸腾实验的数据处理上。

用MATLAB把数据处理的部分全部做好了,文档的结构还是要和以前一样,一个文件夹内包含若干组实验的数据文件夹,每个实验的数据文件夹内直接包含不同瓦数的测试数据xls文件。然后将所有瓦数中的温度取平均值汇总到一个以该实验文件夹名字为名的xls中,并且将后续的数据计算填好。其中的要点包括:

  • 文件系统的路径、文件列表和文件复制粘贴的操作
  • xls的读取和写入操作(之前是对csv的操作)
  • 矩阵、数组的合并、截取,简单计算函数等

%%本程序用于处理大量池沸腾数据文档的统计和计算%%
% 文件夹结构必须为:
% ...\xxxxxx实验\
% ..............\第1组数据\
% ..............\第2组数据\30w.xls
% 对应的变量为:
% workingdir\
% \alldirs(di).name\
% \temp1(i),name
% 输出为:
% workingdir\第1组数据.xls
%2017年5月5日改
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
workingdir='D:\work\Project boiling\Experiment\1705-3dpoolboiling_by_chen\处理后数据\';%文件夹的路径
alldirs=dir([workingdir]);
for di=3:length(alldirs)
currentdir=[workingdir,alldirs(di).name,'\'];
temp1=dir([currentdir,'*.xls']);
%[num,txt,raw] = xlsread('D:\work\Other projects\MATLAB\mb.xlsx')
writename=[workingdir,alldirs(di).name,'.xls']
copyfile('D:\work\Other projects\MATLAB\mb.xlsx',writename);
for i=1:length(temp1)
filename=[currentdir,temp1(i).name];
%filename='E:\文档 池沸腾数据\QF58\QF58-310.CSV'
FileID = fopen(filename);
%%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%s%s%s%s%f%s
%读取xls文档的方法如下,分为数字矩阵、字符串矩阵和原始矩阵进行读取
[num,txt,raw] = xlsread(filename,'A:G');
%data = str2num(txt{1:7})
fclose(FileID);
zsum=[];
%for t=17:23
%这个循环是用来切分7个通道的数据为7份,然后每个通道数据直接用mean函数求均值,再将均值ysum组合成一行的数据,至此所有原始数据使用完毕
for t=1:7
y=txt(2:length(txt),t);
y=str2double(y);
ysum=mean(y);
%data{t}=[data{t};ysum]
zsum=[zsum ysum];
end
%一堆处理数据的公式
%tave顶上3根热电偶的平均温度
tave=(zsum(1,1)+zsum(1,2)+zsum(1,3))/3;
%dtdz梯度,热电偶距离为0.65cm
dtdz=(zsum(1,6)-zsum(1,5)+zsum(1,4)-zsum(1,3))/2/0.65;
%tw壁面(样品)温度,和tave之间距离0.6cm
tw=tave-0.6*dtdz;
%dtsat/dt饱和过热度和实际过热度,用壁面温度减工质沸点和工质实际温度
dtsat=tw-78;
dt=tw-zsum(1,7);
%qpp热流密度啊大哥
qpp=dtdz*3.77;
%h传热系数
h=qpp/dt*10;
%dlmwrite(filename,[0 zsum],'-append')
w_name={temp1(i).name(6:end-5)};%{}把原来的string或数字放在了一个cell里面,这样输出就不会打散字符串
xlswrite(writename,w_name,1,['A' num2str(i+1)]);
xlswrite(writename,[zsum tave dtdz tw dtsat dt qpp h],1,['B' num2str(i+1)]);
end
end
jobdone= 'JOB DONE!!!!!!'%不想输出太乱所以把大部分都隐藏起来了,只剩下每个输出名会输出,作为进度提示,最后输出一句jobdone提示完成。

2017年5月4日 周四

INDENTATION & SCRATCH SIMULATION WITH ABAQUS

1、确认了本构,新发现如下:

  • 本构塑性部分有问题,塑性为正确值的1/10,希望验证结果塑性流动会更明显。
  • 本构分两部分,一部分是压缩时候的弹-塑性行为,无脆性破坏;另一部分是拉伸时候的弹-迅速破坏的脆性行为。
  • 多裂纹先后产生仍无法实现,正在想办法把激活区域变成小区域的类型,在压头正下方定一个小区域,在两侧考虑定两个小区域,将两者区分开来。如果做刻划的话,要做一个横放的圆锥形区域把刻划划痕包住,并且建立激活区,让裂纹能快速穿过激活区,允许其他裂纹产生。

POOL BOILING AND DATA PROCESSING

2、用MATLAB编写程序处理池沸腾的数据,结果基本定型为如下:

QQ截图20170504225318

首先输入的数据文档和以前的不同,以前是csv,现在是xls。csv的长处是好读,很好处理,坏处是输出函数不能将不同类型数据合并输出。xls的读取一般只能读数,不能读字符串,要读字符串还要用上语句:
 [num,txt,raw] = xlsread(filename,xlsheet,xlrange)
然后自己分取num、txt的矩阵进行使用。由于上位机的作者非常脑残,导致输出数据全是字符串,因此后续还有做个变换str2double(str2num也行,但是温度数据最好double),然而这个变换只能做一列的,不能做一个大矩阵,正好要用一列的求平均值,所以不是太麻烦。

明天还需要做的一点就是把后面的Tave等等数据全部再matlab里面用公式求好并输出到csv文件中,大概还要花一个早上

3、仔细看了一下之前邓师兄写的论文,和师弟讨论了一下后,大概明白师兄做池沸腾分析数据的思路,很遗憾没有很早做这个工作,实际上这类论文真的很好写的。

  • 先算样品表面温度Tw,用Tave+dT/dz*l即可,简单
  • 关键是dT/dz,为了减少焊锡层对数据的影响,这里用的是((T6-T5)+(T4-Tave))/2求的,舍去了5号和4号之间的差值。
  • 数据处理更好的方法是用所有瓦数的温度数据减去0w时候的数据,消除所有误差,这样想来实际上真的可以用类似思想完全消除单个热电偶的系统误差,还有焊锡层的影响。明天重点研究这个。
  • 处理完后还有DTsat和DT两个温度差值要求,一个是表面温度减去工质沸点,一个是表面温度减去工质温度。
  • 再有热流密度q”是瓦数除以样品面积,传热系数h为热流密度除以DT温差
  • 然后就是两个图的套路:热流密度——过热度;传热系数——热流密度

具体的分析明天认真看了才能得出结论。

附件

整理实验数据的matlab程序代码,供以后参考
%UNTITLED1 Summary of this function goes here
% Detailed explanation goes here
workingdir='D:\work\Project boiling\Experiment\1705-3dpoolboiling_by_chen\matlab_test\';%文件夹的路径
alldirs=dir([workingdir]);
for di=3:length(alldirs)
currentdir=[workingdir,alldirs(di).name,'\']
temp1=dir([currentdir,'*.xls']);
for i=1:length(temp1)
filename=[currentdir,temp1(i).name]
%filename='E:\文档 池沸腾数据\QF58\QF58-310.CSV'
FileID = fopen(filename)
%%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%s%s%s%s%f%s
[num,txt,raw] = xlsread(filename,'A:G')
%data = str2num(txt{1:7})
fclose(FileID);
zsum=[]
%for t=17:23
for t=1:7
y=txt(2:length(txt),t)
y=str2double(y);
ysum=mean(y)
%data{t}=[data{t};ysum]
zsum=[zsum ysum]
end
%dlmwrite(filename,[0 zsum],'-append')
if filename(end-6:end-4)=='ONB'
%将ONB和CHF排除在外,如果要视察可以手动打开文档
else
if filename(end-6:end-4)=='CHF'
else
writename=[alldirs(di).name,'.csv']
w_name=str2num(temp1(i).name(6:end-5))
dlmwrite(writename,[w_name zsum],'-append')
end
end
end
end