基于python二次开发的围线积分(contour integral)+网格重划分(remeshing)+结果映射

1、前面帖子已经详细介绍了如何使用围线积分(contourintegral)+网格重划分(remeshing)来模拟裂纹的任意路径扩展，并提取裂纹扩展路径上的应力强度因子。那么在这个帖子上我们进一步升级原理算法，使得模拟与实际实验过程更加符合，具体内容如下：首先，我们对之前的裂纹扩展原理进行大体的回顾，程序设计的总体原理图如下，旧裂纹：裂纹扩展一个增量步：这个裂纹扩展分析过程与看到的几乎所有文献是一样的，是一个常规的裂纹扩展分析过程：提取分析步末（或者取第一次应力强度因子大于临界值的那个时刻）的应力强度因子，根据最大周向拉应力开裂准则计算复合I型应力强度因子，如果裂纹扩展，进行新裂纹的计算，在原来

2、模型的基础之上，根据新的裂纹角度及裂尖坐标，sketch出新的裂纹，然后重新进行网格划分，但是，分析步、载荷都并进行任何改变，而是将整个载荷又重新加载一遍，再次得到载荷分析步末（或者取第一次应力强度因子大于临界值的那个时刻）的应力强度因子，再进行裂纹是否扩展的判断，这相当于把离散的不同裂纹长度单独计算模型整合成了考虑裂纹扩展角度的连续性裂纹扩展。这样的裂纹扩展过程只能得到裂纹扩展路径与实验大体一致，在考虑时间的分析中，不能得到裂纹随时间的扩展过程，反而，（例如薄板中心裂纹）因为裂纹长度越大，应力强度因子越大，所以得到更长的裂纹对应更短的开裂时间，这与实验过程是相悖的。那么，问题来了：实际实验中

3、裂纹是在加载过程中连续扩展的，而并不是裂纹扩展后载荷归零又重新从零加载了，而是裂纹在某时刻开裂后，在剩余载荷的进一步的加载过程中裂纹进一步扩展。通过前面的分析我们知道，利用前面的裂纹扩展原理，通常对于一些静力结构分析来说，可以得到相对合理的裂纹扩展路径，但是，每次裂纹扩展分析都是全部载荷重新施加一遍，并不能取到裂纹长度随时间的变化曲线，而且对于与时间相关的分析而言，裂纹长度随时间的变化就尤为重要了，那么，我们的优化就是解决这一个问题，解决办法：在前面的理论基础上，利用ABAQUS的inp文件中的关键词*mapsolution来将旧网格的结果映射到新网格，施加剩余载荷，进行下一次的裂纹扩展。这个

4、裂纹分析过程似乎看起来很容易，但是实际实施起来困难重重，我们经过长时间的探索终于克服了重重困难，实现了整个裂纹动态扩展过程。下面介绍裂纹扩展分析的整体原理，II否舟折步.建立Ccracl（历史输出巫启动输出）修改Inp.潘加*mmpsolutionM;TpythOn_fe开发的井CooniMrlHtegMiMJM*製纹随机扩展程俘原理團计算新裂纹尖端，裂纹扩展一个增量.建宜分析歩（銅余时闸段，历史输出，重启动输出:ASAUUScommand握交穆改后的伯p丈件夙Mb提取开裂时刻的几何变形模理裂飯血阖拣触筑立（:防止单元穿透建立job.谕出inp曙寵黔15匝新建立的束亟新建立接触重新网恪划廿几何

5、模型建立重摘装配网辂划分joh堤交|边界棗件施加|裁荷施加初始模型建立分析步后要设置重启动输出（*mapsolution的计算需要重启动文件），如果裂纹扩展，计算新的裂纹尖端，提取odb开裂时刻的几何变形模型，裂纹扩展一个增量，再次重新分配材料属性，建立独立实体，创建新的分析步，时间长度为上次开裂后的剩余时间段，重新建立接触，重新施加约束，重新网格划分，重新施加载荷（载荷是对应上次开裂后的剩余载荷段），输出inp文件，在i叩文件中添加*mapsolution命令,然后使用命令行提交新的inp文件。需要注意的是：从odb提取几何变形模型后，原来完整的模型边缘会变成由很多小的线段组成（小线段的多少

6、与提取特征角度featureangle和网格密度有关），小线段的多少也影响了后续的载荷及约束的施加；由于1中小线段的生成，所以载荷和约束的施加变的尤为困难，主要是因为边的选取通过findat命令也很难选择到，对于不同的模型，我们要视情况而定；每次裂纹扩展的增量不要太大，因为，如果增量太大，新的裂纹尖端生成后，旧的网格结果映射到新的网格上，会导致裂纹尖端应力出现较大偏差，导致不合理的结果。介绍了整个分析原理后，我们通过编制的二次开发程序分析一个例子来验证程序的正确性，采用薄板边缘裂纹，施加了位移载荷。初始几何模型，约束，载荷，网格：第一次裂纹扩展后，提取的裂纹扩展时刻的几何变形图及旧裂纹尖端与新

7、生成的seam裂纹图示如下：分析步的重新建立（第一次开裂后的剩余时间段）斗；ElepMari.gerNameProcedLireNlgeomTame19InitialN/AN；A1VStep-1Stac;GeneralOFF0.7d41507279S72B9载荷及约束的施加（剩余载荷段，总体载荷位移为0.05mm）lEdrt.BoundjryConditionName：T0P21Typre;DisplaLement/RotationStep:&tep-lStatic.GeneralRegion:T21RUR3:Amplitude:CSYi(Global)Disthbutian:Ns-tg1：T

8、hrdispla匚emuntvalu?willbemHiniainedinSLibsquentstops.OK重新划分后的网格模型：修改的inp文件:*mapsolutiQiirsteplfinc=36STE?TStep-1Step/naitie=33ep-l/nlgeom=NOjinc=L0000*StaticD.0001f0.744151,-le-lz,0.01*+BOUNDAEYCONDITIONS命令提交job：iTtabaigLisjob=cpaclc2o丄djoh=ciacklcpus=6interactiueIitelParal丄匕丄oluuiu矗工迅丄Gopyriglit198

9、52014Intc1Covporationfi11rigihtsreserved-IntelCRiConpoeeiHE2013GP1Update3(paclagie262J裂纹扩展过程：裂纹经过40次扩展后即将扩展到模型边界，从图中可以看出，在裂纹扩展过程中，已经张开的裂纹面附近可能出现较大的应力值，这是由于mapsolution映射结果导致的，这并不影响裂纹的扩展，因此可以忽略。每次裂纹扩展时的K1与裂纹扩展对应的总时间的关系：每次裂纹扩展时的KL与裂纹扩展对应的总片间的关系从图中可以看出，裂纹扩展过程中只要裂纹尖端应力强度因子达到临界值，裂纹便发生扩展值得注意的是裂纹扩展过程中裂纹尖端的K1并不会逐渐递增，这是因为裂纹尖端的k1值并不是提取的每次裂纹扩展加载末的K1,下面我们给出裂纹扩展过程中每次夹杂