03-线状裂纹的虚拟裂纹闭合法

1HUST-XD-B-0111-2014线状裂纹的虚拟裂纹闭合法233.1虚拟裂纹闭合法的数学解释目前有两种方式对虚拟裂纹闭合法进行数学解释：1Raju方法

基本假设：虚拟裂纹尖端后面的张开位移和实际裂纹尖端后面的张开位移近似相等。2Xie-Waas方法基本假设：虚拟子裂纹尖端前沿的应力分布和实际主裂纹尖端前沿的应力分布规律近似相同。43.1虚拟裂纹闭合法的数学解释3.1.1Raju方法图3.1显示了一个长度为a的主裂纹，虚拟扩展了一个增量为

a的微小裂纹。Irwin的裂纹闭合积分可以写为：图3.1虚拟裂纹闭合法示意图5将裂纹从a扩展到a+a所需要的功与将裂纹从a+a闭合到a所需的功是相等的。如图3.2所示，考虑裂纹尖端附近的有限元网格，虚拟裂纹线上分布的应力做的功等于节点力在节点位移上做的功。即：3.1.1Raju方法3.1虚拟裂纹闭合法的数学解释分布于裂纹尖端前沿且趋近于裂纹尖端区域的应力可以表示为：6对于线性四边形单元，位移可以插值为：3.1.1Raju方法将公式3.4和3.3代入公式3.2中，并消去得：依据基本假设：虚拟裂纹尖端后面的张开位移和初始裂纹尖端后面的张开位移近似相等，公式3.1可重写为：公式3.6的优点是它可以在一步有限元分析中计算出来。如图3.1所示，裂纹尖端前面的应力分布为：

基于内插值，裂纹尖端后面的张开位移为：

将公式3.7和公式3.8代入公式3.6中，得到：3.1虚拟裂纹闭合法的数学解释73.1.1Raju方法3.1虚拟裂纹闭合法的数学解释

公式3.9括号中部分简化为：

公式3.10中的极限表明裂纹尖端的单元需要很小，并且裂纹尖端处的单元越小，所得的GI就越精确，通过数值实验可以确定必要的网格细化程度。然而，大量的数值实验表明公式3.10的值对有限元网格尺寸并不敏感。因此，公式3.10中的极限可以很好的近似的表达为：

类似的，可以得到II型的计算公式如下：83.1.2Xie-Waas方法3.1虚拟裂纹闭合法的数学解释如图3.3所示，一个主裂纹（parentcrack），其应力强度因子为KI。主裂纹前端有一个虚拟的微小子裂纹（childcrack），其应力强度因子为。根据线弹性断裂力学：上述公式对于主裂纹和子裂纹都是成立的。图3.3主裂纹和虚拟示意图子裂纹93.1.2Xie-Waas方法3.1虚拟裂纹闭合法的数学解释通过虚拟裂纹扩展法（参看公式1.26a），并参考图3.2，能量释放率可用下式计算：第一步计算的节点力第二步计算的节点张开位移对于主裂纹而言，上平面（）和下平面（）的垂直位移分别为：

主裂纹后部的张开位移很容易得到，为：

同样地，对虚拟子裂纹而言，有15IVIIIIII34ΔaΔcXY1534XY1’IVIIIIII注意上标不同！103.1虚拟裂纹闭合法的数学解释3.1.2Xie-Waas方法

接下来建立和之间的关系：

对于主裂纹而言，沿x轴方向，裂纹前端应力分布为：

依据基本假设：虚拟子裂纹尖端前沿的应力分布和初始裂纹尖端前沿的应力分布规律近似相同。式3.16所描述的应力分布同样适用于裂纹。因此，根据应力强度因子的定义，有注意上标不同！

当裂纹尖端使用相同的网格尺寸时，即，公式3.19退化为：

由于F和v均只是在第一步中计算，因此可以省去上标，这也是公式3.20的优点！113.1虚拟裂纹闭合法的数学解释3.1.3节点力的计算使用虚拟裂纹闭合法时涉及节点力和节点位移。由于节点位移是有限元分析中的基本变量，因而可以直接获得，相对简单，但在计算节点力时则应该注意。注意到计算中使用的节点力是内力，内力将裂纹尖端的单元共享节点约束在一起。因此，裂纹尖端周围的单元对这些节点力都有贡献。图3.4显示了在裂纹尖端周围有四个单元的情况下的自由体受力图。如果没有外力作用在裂纹面上，节点力在每个坐标方向上的合力为0，即：上标代表图中的单元编号。因此，节点力既可以通过上面的两个单元（单元1和单元2）计算，也可以通过下面的两个单元计算：IIVIIIII12*HEADINGABAQUSJOB=Chapter_03_example_3_1INT*PREPRINT,HISTORY=NO,MODEL=NO,CONTACT=NO,ECHO=YES*PARAMETERB=50.0YM=200E3XNU=0.3*NODE1,0.0,0.01001,1000.0,0.010011,0.0,10.011011,1000.0,10.020001,0.0,0.020500,499.0,0.021002,0.0,-1.022002,1000.0,-1.030011,0.0,-10.031011,1000.0,-10.0*NGEN,NSET=botup1,1001,1*NGEN,NSET=topup10011,11011,1*NGEN,NSET=botdn20001,20500,1*NGEN,NSET=middn21002,22002,1*NGEN,NSET=topdn30011,31011,1*NFILLbotup,topup,10,1001middn,topdn,9,1001*ELEMENT,TYPE=CPS4,ELSET=arm1,1,2,1003,100220001,20001,21002,21003,2000220500,20500,21501,21502,50120501,501,21502,21503,50221001,21002,22003,22004,21003*ELGEN,ELSET=arm1,1000,1,1,10,1001,100020001,499,1,120501,500,1,121001,1000,1,1,9,1001,1000*SOLIDSECTION,ELSET=arm,MATERIAL=st4340<B>*MATERIAL,NAME=st4340*ELASTIC<YM>,<XNU>*NSET,NSET=load5006,25006*ELSET,ELSET=cracktip500,501,20500,20501*NSET,NSET=open500,501,20500,20501*NSET,NSET=open500,20500*STEP*STATIC,DIRECT1.0,1.0*BOUNDARY5006,125006,15006,2,2,60.025006,2,2,-60.01001,2*ELPRINT,FREQ=0*NODEPRINT,FREQ=0*ELPRINT,ELSET=cracktipNFORC*NODEPRINT,NSET=openU*NODEPRINT,NSET=loadU,RF*OUTPUT,FIELD*NODEOUTPUTU*ENDSTEPNFORC:Forcesatthenodesofabeamelementfromboththehourglassandtheregulardeformationmodesofthatelement.3.1虚拟裂纹闭合法的数学解释3.1.3节点力的计算133.1虚拟裂纹闭合法的数学解释3.1.3节点力的计算ELEMENTFOOT-NODENFORC1NFORC2NOTE5005001.6253E+04-32.52500501-2.6590E+041.9915E+045001502-1.1234E+04-9578.50015012.1572E+04-1.0305E+045015012.6590E+042.0316E+04501502-1.9036E+046297.5011503-2.0021E+04-1.3851E+0450115021.2466E+04-1.2762E+0420500205001.6253E+0432.5220500215012.1572E+041.0305E+042050021502-1.1234E+049578.20500501-2.6590E+04-1.9915E+04205015012.6590E+04-2.0316E+0420501215021.2466E+041.2762E+042050121503-2.0021E+041.3851E+0420501502-1.9036E+04-6297.通过上面单元（500,501）通过下面单元（20500,20501）14注意，此时需要对裂纹尖端节点进行处理！裂纹尖端两个节点，节点用弹簧单元连接！裂纹尖端一个节点！*ELEMENT,TYPE=SPRING2,ELSET=spx50001,501,20501*SPRING,ELSET=spx1,1<ST>*ELEMENT,TYPE=SPRING2,ELSET=spy50002,501,20501*SPRING,ELSET=spy2,2<ST>3.1虚拟裂纹闭合法的数学解释3.1.3节点力的计算15THEFOLLOWINGTABLEISPRINTEDATTHEINTEGRATIONPOINTSFORELEMENTTYPESPRING2ANDELEMENTSETSPXELEMENTPTFOOT-S11NOTEALLVALUESINTHISTABLEAREZEROTHEFOLLOWINGTABLEISPRINTEDATTHEINTEGRATIONPOINTSFORELEMENTTYPESPRING2ANDELEMENTSETSPYELEMENTPTFOOT-S11NOTE5000214.0232E+043.1虚拟裂纹闭合法的数学解释3.1.3节点力的计算163.2平面内任意线状裂纹的断裂单元考虑图3.6所示的位于（X,Y）平面内的倾斜裂纹。虚拟裂纹扩展量是主节点1和哑节点5之间的距离：式中，和分别为节点1和节点5在整体坐标系（X,Y）下的坐标。如果它们在每个增量步中都通过位移来更新，那么裂纹的方向也将更新，这在分析大变形问题时很有用。为分离断裂模式（I型和II型），应变能释放率必须在裂纹尖端处的局部坐标系下计算，如图3.6.X轴和轴之间的夹角为：应变能释放率近似的表示为：173.2平面内任意线状裂纹的断裂单元183.2平面内任意线状裂纹的断裂单元193.2平面内任意线状裂纹的断裂单元203.3经典I型问题上的应用由于绝大多数的断裂参数解析解都以应力强度因子的形式给出，因此需要将虚拟裂纹法所计算的能量释放率转换成应力强度因子来加以比较。对线弹性材料而言，两者的关系如下：平面应力状态：平面应变状态：应力强度因子准则和能量准则是从两种不同观点建立起来的准则，它们之间有一定的关系。假设裂纹沿其延长线扩展时，计算裂纹尖端应力松弛的功，可求得能量释放率和应力强度因子的关系。

设I型裂纹尖端处于平面应变状态，在固定施力点的情况下裂纹沿其延长线扩展，裂纹尺寸从a增加到a+da，如图所示。此时，裂纹扩展段da上的应消失，新表面在法线方向产生位移。比较裂纹长度为a+da和裂纹长度为a时板的能量差异就可以求出能量释放率。先把裂纹长度看做a+da，计算张开位移（x的原点在裂纹尖端）：21正负号分别对应于上下表面。下角标a+da表示裂纹长度为a+da时的应力强度因子。在da段裂纹表面上加载，载荷的分布和大小按照表面的应力大小取值：加载值由0增加到上式的最终值，da段裂纹逐渐靠拢知道完全闭合，恢复到裂纹为a的状态。因此从a+da状态恢复到a状态的加载过程中，外力对弹性体做了功，提高了应变能。由于任意x处的和之间存在线性关系，做功为：式中，B为物体的厚度。反之，从a状态变为a+da状态，应力从松弛到零，位移从零变化到，松弛过程中释放的应变能dV=-dW。所以：3.3经典I型问题上的应用223.3经典I型问题上的应用233.3经典I型问题上的应用24ELEMENTPTFOOT-SDV1SDV2SDV3SDV4SDV5SDV6SDV7NOTE5000112.9000.0002.4000.0001.0000.0000.2000a=2.4ELEMENTPTFOOT-SDV1SDV2SDV3SDV4SDV5SDV6SDVNOTE

5000114.4470.0003.4000.0001.0000.0000.2000a=3.4ELEMENTPTFOOT-SDV1SDV2SDV3SDV4SDV5SDV6SDVNOTE

5000116.4310.0004.4000.0001.0000.0000.2000a=4.4ELEMENTPTFOOT-SDV1SDV2SDV3SDV4SDV5SDV6SDVNOTE

5000119.2030.0005.4000.0001.0000.0000.2000a=5.4ELEMENTPTFOOT-SDV1SDV2SDV3SDV4SDV5SDV6SDVNOTE

50001113.500.0006.4000.0001.0000.0000.2000a=6.4ELEMENTPTFOOT-SDV1SDV2SDV3SDV4SDV5SDV6SDVNOTE

50001121.190.0007.4000.0001.0000.0000.2000a=7.4ELEMENTPTFOOT-SDV1SDV2SDV3SDV4SDV5SDV6SDVNOTE

50001138.740.0008.4000.0001.0000.0000.2000a=8.4ELEMENTPTFOOT-SDV1SDV2SDV3SDV4SDV5SDV6SDVNOTE

500011113.50.0009.4000.0001.0000.0000.2000a=9.43.3经典I型问题上的应用25clcclear%Units(Nmm)E=200E3v=0.3sig=0.4E3a=2.4:1:9.4;W=10;x=a/W;

k1=sqrt(sec(pi/2.*x)).*(1-0.25.*x.^2+0.06.*x.^4)k2=sqrt(sec(pi/2.*x))k3=sqrt((2./pi./x).*tan(pi/2.*x))plot(x,k1,'ro');holdonplot(x,k2,'bo');holdonplot(x,k3,'ko');holdon

%aSDV1SDV2kkk=[2.42.9000.000

3.44.4470.000

4.46.4310.000

5.49.2030.000

6.4