03-線狀裂紋的虛擬裂紋閉合法

1HUST-XD-B-0111-2014線狀裂紋的虛擬裂紋閉合法233.1虛擬裂紋閉合法的數(shù)學解釋目前有兩種方式對虛擬裂紋閉合法進行數(shù)學解釋：1Raju方法

基本假設：虛擬裂紋尖端后面的張開位移和實際裂紋尖端后面的張開位移近似相等。2Xie-Waas方法基本假設：虛擬子裂紋尖端前沿的應力分布和實際主裂紋尖端前沿的應力分布規(guī)律近似相同。43.1虛擬裂紋閉合法的數(shù)學解釋3.1.1Raju方法圖3.1顯示了一個長度為a的主裂紋，虛擬擴展了一個增量為

a的微小裂紋。Irwin的裂紋閉合積分可以寫為：圖3.1虛擬裂紋閉合法示意圖5將裂紋從a擴展到a+a所需要的功與將裂紋從a+a閉合到a所需的功是相等的。如圖3.2所示，考慮裂紋尖端附近的有限元網(wǎng)格，虛擬裂紋線上分布的應力做的功等于節(jié)點力在節(jié)點位移上做的功。即：3.1.1Raju方法3.1虛擬裂紋閉合法的數(shù)學解釋分布于裂紋尖端前沿且趨近于裂紋尖端區(qū)域的應力可以表示為：6對于線性四邊形單元，位移可以插值為：3.1.1Raju方法將公式3.4和3.3代入公式3.2中，并消去得：依據(jù)基本假設：虛擬裂紋尖端后面的張開位移和初始裂紋尖端后面的張開位移近似相等，公式3.1可重寫為：公式3.6的優(yōu)點是它可以在一步有限元分析中計算出來。如圖3.1所示，裂紋尖端前面的應力分布為：

基于內插值，裂紋尖端后面的張開位移為：

將公式3.7和公式3.8代入公式3.6中，得到：3.1虛擬裂紋閉合法的數(shù)學解釋73.1.1Raju方法3.1虛擬裂紋閉合法的數(shù)學解釋

公式3.9括號中部分簡化為：

公式3.10中的極限表明裂紋尖端的單元需要很小，并且裂紋尖端處的單元越小，所得的GI就越精確，通過數(shù)值實驗可以確定必要的網(wǎng)格細化程度。然而，大量的數(shù)值實驗表明公式3.10的值對有限元網(wǎng)格尺寸并不敏感。因此，公式3.10中的極限可以很好的近似的表達為：

類似的，可以得到II型的計算公式如下：83.1.2Xie-Waas方法3.1虛擬裂紋閉合法的數(shù)學解釋如圖3.3所示，一個主裂紋（parentcrack），其應力強度因子為KI。主裂紋前端有一個虛擬的微小子裂紋（childcrack），其應力強度因子為。根據(jù)線彈性斷裂力學：上述公式對于主裂紋和子裂紋都是成立的。圖3.3主裂紋和虛擬示意圖子裂紋93.1.2Xie-Waas方法3.1虛擬裂紋閉合法的數(shù)學解釋通過虛擬裂紋擴展法（參看公式1.26a），并參考圖3.2，能量釋放率可用下式計算：第一步計算的節(jié)點力第二步計算的節(jié)點張開位移對于主裂紋而言，上平面（）和下平面（）的垂直位移分別為：

主裂紋后部的張開位移很容易得到，為：

同樣地，對虛擬子裂紋而言，有15IVIIIIII34ΔaΔcXY1534XY1’IVIIIIII注意上標不同！103.1虛擬裂紋閉合法的數(shù)學解釋3.1.2Xie-Waas方法

接下來建立和之間的關系：

對于主裂紋而言，沿x軸方向，裂紋前端應力分布為：

依據(jù)基本假設：虛擬子裂紋尖端前沿的應力分布和初始裂紋尖端前沿的應力分布規(guī)律近似相同。式3.16所描述的應力分布同樣適用于裂紋。因此，根據(jù)應力強度因子的定義，有注意上標不同！

當裂紋尖端使用相同的網(wǎng)格尺寸時，即，公式3.19退化為：

由于F和v均只是在第一步中計算，因此可以省去上標，這也是公式3.20的優(yōu)點！113.1虛擬裂紋閉合法的數(shù)學解釋3.1.3節(jié)點力的計算使用虛擬裂紋閉合法時涉及節(jié)點力和節(jié)點位移。由于節(jié)點位移是有限元分析中的基本變量，因而可以直接獲得，相對簡單，但在計算節(jié)點力時則應該注意。注意到計算中使用的節(jié)點力是內力，內力將裂紋尖端的單元共享節(jié)點約束在一起。因此，裂紋尖端周圍的單元對這些節(jié)點力都有貢獻。圖3.4顯示了在裂紋尖端周圍有四個單元的情況下的自由體受力圖。如果沒有外力作用在裂紋面上，節(jié)點力在每個坐標方向上的合力為0，即：上標代表圖中的單元編號。因此，節(jié)點力既可以通過上面的兩個單元（單元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虛擬裂紋閉合法的數(shù)學解釋3.1.3節(jié)點力的計算133.1虛擬裂紋閉合法的數(shù)學解釋3.1.3節(jié)點力的計算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注意，此時需要對裂紋尖端節(jié)點進行處理！裂紋尖端兩個節(jié)點，節(jié)點用彈簧單元連接！裂紋尖端一個節(jié)點！*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虛擬裂紋閉合法的數(shù)學解釋3.1.3節(jié)點力的計算15THEFOLLOWINGTABLEISPRINTEDATTHEINTEGRATIONPOINTSFORELEMENTTYPESPRING2ANDELEMENTSETSPXELEMENTPTFOOT-S11NOTEALLVALUESINTHISTABLEAREZEROTHEFOLLOWINGTABLEISPRINTEDATTHEINTEGRATIONPOINTSFORELEMENTTYPESPRING2ANDELEMENTSETSPYELEMENTPTFOOT-S11NOTE5000214.0232E+043.1虛擬裂紋閉合法的數(shù)學解釋3.1.3節(jié)點力的計算163.2平面內任意線狀裂紋的斷裂單元考慮圖3.6所示的位于（X,Y）平面內的傾斜裂紋。虛擬裂紋擴展量是主節(jié)點1和啞節(jié)點5之間的距離：式中，和分別為節(jié)點1和節(jié)點5在整體坐標系（X,Y）下的坐標。如果它們在每個增量步中都通過位移來更新，那么裂紋的方向也將更新，這在分析大變形問題時很有用。為分離斷裂模式（I型和II型），應變能釋放率必須在裂紋尖端處的局部坐標系下計算，如圖3.6.X軸和軸之間的夾角為：應變能釋放率近似的表示為：173.2平面內任意線狀裂紋的斷裂單元183.2平面內任意線狀裂紋的斷裂單元193.2平面內任意線狀裂紋的斷裂單元203.3經(jīng)典I型問題上的應用由于絕大多數(shù)的斷裂參數(shù)解析解都以應力強度因子的形式給出，因此需要將虛擬裂紋法所計算的能量釋放率轉換成應力強度因子來加以比較。對線彈性材料而言，兩者的關系如下：平面應力狀態(tài)：平面應變狀態(tài)：應力強度因子準則和能量準則是從兩種不同觀點建立起來的準則，它們之間有一定的關系。假設裂紋沿其延長線擴展時，計算裂紋尖端應力松弛的功，可求得能量釋放率和應力強度因子的關系。

設I型裂紋尖端處于平面應變狀態(tài)，在固定施力點的情況下裂紋沿其延長線擴展，裂紋尺寸從a增加到a+da，如圖所示。此時，裂紋擴展段da上的應消失，新表面在法線方向產(chǎn)生位移。比較裂紋長度為a+da和裂紋長度為a時板的能量差異就可以求出能量釋放率。先把裂紋長度看做a+da，計算張開位移（x的原點在裂紋尖端）：21正負號分別對應于上下表面。下角標a+da表示裂紋長度為a+da時的應力強度因子。在da段裂紋表面上加載，載荷的分布和大小按照表面的應力大小取值：加載值由0增加到上式的最終值，da段裂紋逐漸靠攏知道完全閉合，恢復到裂紋為a的狀態(tài)。因此從a+da狀態(tài)恢復到a狀態(tài)的加載過程中，外力對彈性體做了功，提高了應變能。由于任意x處的和之間存在線性關系，做功為：式中，B為物體的厚度。反之，從a狀態(tài)變?yōu)閍+da狀態(tài)，應力從松弛到零，位移從零變化到，松弛過程中釋放的應變能dV=-dW。所以：3.3經(jīng)典I型問題上的應用223.3經(jīng)典I型問題上的應用233.3經(jīng)典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經(jīng)典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