基于VCCT的裂纹扩展模拟

.基于VCCT的裂纹扩展模拟虚拟裂纹闭合技术(VCCT)最初用于计算裂纹体的能量释放率。因此被广泛用于层合复合材料的界面裂纹扩展模拟，并假定裂纹扩展总是沿着预先定义的路径，特别是在界面处。基于VCCT的裂纹扩展模拟,当前可用的线性单元如plane182和solid185。基于VCCT的裂纹扩展模拟包含下述假定：裂纹扩展沿着预先定义的路径路径通过界面单元来定义分析为准静态分析，不考虑瞬态效应材料为线弹性材料，可以是各向同性，正交各向异性,各向性。裂纹可以位于一种材料或者两种材料的界面。断裂准则基于采用VCCT方法计算的能量释放率。可采用多种断裂准则或自定义的准则。同一分析中可定义多条裂纹。VCCT裂纹扩展模拟使用： 界面单元INTER202（2D）和INTER205（3D） CINT命令计算能量释放率 CGROW命令定义裂纹扩展集，断裂准则，裂纹扩展路径和求解控制参数。12.1.1 VCCT裂纹扩展模拟过程基于VCCT的裂纹扩展模拟假定为准静态模拟。下面为进行模拟的主要步骤： Step1：建立预先定义裂纹路径的有限元模型 Step2：进行能量释放率的计算 Step3：进行裂纹扩展计算裂纹扩展模拟为非线性结构分析，这里详述了一些特点，特别是裂纹扩展的分析细节。12.1.1.1 Step1：建立预先定义裂纹路径的有限元模型标准的非线性求解过程需要建立有限元模型，有正确的求解控制设置，载荷和边界条件。预先定义的裂纹路径离散为界面单元，并分为一个单元组，如下图所示：图12.1采用界面单元离散裂纹路径界面单元可以通过CZMESH命令划分或者能生成界面单元的第三方工具划分。MPC约束单元选项(KEYOPT(2) = 1) 在裂纹扩展前把潜在的裂纹面绑定在一起。当满足断裂准则时，MPC约束随后释放，从而裂纹扩展。在二维问题中，裂纹尖端后的一个界面单元如果在一个指定的子步满足断裂准则则可能张开。在三维问题中，裂纹前沿后的所有界面单元如果满足断裂准则可能张开。裂纹尖端/前沿周围的单元尺寸影响能量释放率的计算精度。当程序采用修正算法,可能不能产生精确的结果。改为沿着裂纹扩展路径使用相同单元尺寸的网格。12.1.1.2 Step2： 进行能量释放率计算基于VCCT的裂纹扩展模拟，必须首先进行能量释放率的计算。计算能量释放率，采用CINT,TYPE,VCCT命令，随后使用CINT命令指定其它选项比如裂纹尖端节点组和裂纹面/边的法向。VCCT计算采用下述假定：当裂纹增加一个小量时，释放的应变能等于裂纹闭合相同的小量所需的能量当裂纹扩展一个小量时，裂纹尖端（前沿）位置的裂尖场（变形）不变。当裂纹扩展接近边界或者两条裂纹彼此接近时假定不在适用。因此，使用VCCT计算要仔细检查分析结果。12.1.1.3.Step 3 进行裂纹扩展计算裂纹扩展计算在应力计算之后，solution步进行。为了进行裂纹扩展计算，必须先定义裂纹扩展集合，然后指定裂纹路径，断裂准则，裂纹扩展求解控制。求解命令CGROW定义裂纹扩展计算所有必需的参数。进行裂纹扩展计算步骤如下： Step 3a：初始裂纹扩展集 Step 3b：指定裂纹路径 Step 3c：指定裂纹计算的ID和断裂准则 Step 3d：指定裂纹扩展的求解控制12.1.1.3.1 Step 3a：初始裂纹扩展集定义裂纹扩展集，使用CGROW,NEW,n命令，其中n是裂纹扩展集的编号12.1.1.3.1 Step 3b：指定裂纹路径定义裂纹路径，采用CGROW,PATH,cmname,其中cmname是界面单元组的名称。12.1.1.3.1 Step 3c：指定裂纹计算的ID和断裂准则指定裂纹计算ID,通过CGROW,CID,n命令，其中n是采用VCCT计算能量释放率的裂纹计算(CINT)的ID。（CINT命令定义的参数和断裂参数计算一致）对于简单的断裂准则，比如临界能量释放率，可以通过命令CGROW,FCOPTION,GC,指定，其中value为临界能量释放率。对于一些更复杂的断裂准则，可以通过材料数据表定义断裂准则。采用CGROW，FCOPTION，MTAB，matid命令，其中matid是材料表的材料ID号。有多种断裂准则可以使用，比如linear,bilinear,B-K,修正B-K，Power Law，和用户自定义的准则。更多信息请参阅TB，CGCR命令和Fracture criterial.对于每个裂纹扩展集，只可以定义一条断裂准则和一个单元组。可以采用不同的断裂准则来定义多个裂纹扩展集。多条裂纹可以同时扩展或者彼此独立。当多条裂纹位于同一界面时，也可以合并为一条裂纹，如下图所示： 图12.2 裂纹扩展和合并 也可以在各自的断裂扩展集中，对同一裂纹定义不同的断裂准则。裂纹可以基于不同的准则（根据哪条准则满足）扩展，并且彼此独立。这种方法对于比较断裂机理很有帮助。使用CGROW 命令 定义求解控制参数如下: 指定求解控制参数 使用CGROW命令:断裂准则系数 (fc) CGROW,FCRAT,value, Value为 该系数初始时间步长 （裂纹扩展开始时）CGROW,DTIME,value, Value为初始时间步长 为避免过度预测承载能力, 指定一个小的初始时间步长随后裂纹扩展的最小时间步长 CGROW,DTMIN,value, value 为最小时间步长大小 随后裂纹扩展的最大时间步长CGROW,DTMAX,value, value 为最大时间步长大小裂纹扩展前沿节点允许的最大裂纹扩展量CGROW,STOP,CEMX,value, value为允许的最大的裂纹扩展量 裂纹扩展模拟很耗时, 当达到感兴趣的指定裂纹扩展量时,使用该命令来终止分析 当裂纹迅速扩展时(比如,裂纹扩展不稳定),使用较小的DTMAX和DTMIN来容许载荷重新平衡。当裂纹不再增长时，指定的时间步长控制被忽略，结果依赖于标准时间步长控制。12.1.1.4.示例:裂纹扩展集定义下面的输入示例定义一个裂纹扩展集：CGROW,NEW,1CGROW,CPATH,cpath1CGROW,FCOPTION,MTAB,5CGROW,DTIME,1.0e-4CGROW,DTMIN,1.0e-5CGROW,DTMAX,2.0e-3.12.1.2裂纹扩展在裂纹扩展模拟中，一个关心的量是裂纹扩展量。VCCT方法测量裂纹扩展是基于已经张开的界面单元的长度，如下述方程和图片所示： a=i图12.3 二维和三维裂纹扩展对于二维问题，裂纹扩展是当前已经张开的界面单元的长度之和（a）。对于三维问题，裂纹扩展在每个裂纹前沿节点进行测量，为沿着裂纹扩展方向的界面单元边长的和（b).裂纹扩展量a（CEXT)是裂纹求解结果的一部分，和裂纹计算的ID号一致，可以和能量释放率一样通过POST1和POST26 后处理命令（比如PRCINT, PLCINT,和CISOL)进行后处理 。12.1.2断裂准则为建立裂纹扩展,必须定义裂纹开始和随后裂纹扩展的断裂准则。对于线弹性断裂力学，断裂准则通常假定为三种断裂模式的临界能量释放率的函数。表达为：f=f(GIC,GIIC,GIIIC,GI,GII,GIII,)对于一些模型可能需要其它的参数。当断裂准则满足时，发生断裂，表述为：ffc其中 fc为断裂准则比率。推荐值为0.95到1.05，默认为1.0。可以使用下列断裂准则：临界能量释放率准则线性断裂准则双线性断裂准则B-K断裂准则修正B-K断裂准则指数断裂准则自定义断裂准则用户自定义选项需要提供子程序来定义你自己的断裂准则。12.1.3.1.临界能量释放率准则临界能量释放率准则使用总的能量释放率（GT）作为断裂准则。总的能量释放率是三种模式的能量释放率的和，表述为： f=GTGTC GT=GI+GII+GIII其中GTC 为临界能量释放率。对于I型断裂模式，断裂准则简化为： f=GIGIC 能量释放率准则是最简单的断裂准则，适用于2D和3D的断裂扩展模拟。 示例12.1临界能量释放率输入 gtcval=10.0 CGROW,FCOPTION,GTC,gtcval12.1.3.2.线性断裂准则线性选项假设断裂准则是三种模式的能量释放率的线性函数。表述为：f=GIGIC+GIIGIIC+GIIIGIIIC其中GIC，GIIC，GIIIC 分别为I型，II型，III型断裂模式的临界能量释放率。这三个值通过命令TBDATA输入，如下：Constant TBDATA Input Comments GIC C1Mode I 的临界能量释放率, GIC 0GIICC2Mode II的临界能量释放率, GIIC 0GIIICC3Mode III的临界能量释放率,GIIIC 0 示例12.2线性准则输入 g1c=10.0 g2c=20.0 g3c=25.0 TB,CGCR,1,LINEAR TBDATA,1,g1c,g2c,g3c三个参数不能同时为零。如果其中一个设为零，相应的项被忽略。当三个临界能量释放率相等，则线性断裂准则简化为临界能量释放率准则。线性断裂准则适用于当三种断裂模式的临界能量释放率明显存在时的三维混合断裂模式。12.1.3.3.双线性断裂准则双线性断裂选项假设断裂准则是I型和II型断裂模式的能量释放率的线性函数，表达式为：其中GIC，GIIC 分别为I型，II型断裂模式的临界能量释放率，和为材料常数。四个值都能通过TBDATA命令定义为温度的函数。如下所示： Example12.3双向性准则输入 g1c=10.0 g2c=20.0 x=2 y=2 TB,CGCR,1,BILINEAR TBDATA,1,g1c,g2c,x,y双线性断裂准则适用于二维混合断裂模式的模拟。12.1.3.4. B-K断裂准则B-K选项表述为：其中GIC，GIIC，分别为I型，II型断裂模式的临界能量释放率,为材料常数。三个值都能够通过TBDATA命令定义为温度的函数。如下所示：B-K准则本来用于复合材料界面断裂，适用于三维混合断裂模式的模拟。 Example12.4B-K 准则输入 g1c=10.0 g2c=20.0 h=2 TB,CGCR,1,BK TBDATA,1,g1c,g2c,h12.1.3.5.修正B-K断裂准则修正B-K选项,表述为:其中GIC，GIIC，GIIIC 分别为I型，II型，III型断裂模式的临界能量释放率,为材料常数。四个值都能够通过TBDATA命令定义为温度的函数。如下所示：GIIC=GIIIC时，修正B-K准则简化为B-K准则。修正B-K准则用于复合材料界面断裂，考虑明显的II型和III型临界能量释放率，适用于三维混合断裂模式的模拟。 Example12.5修正 B-K 准则输入 g1c=10.0 g2c=20.0 g3c=25.0 h=2 TB,CGCR,1,MBK TBDATA,1,g1c,g2c,g3c,h12.1.3.6.指数断裂准则指数准则选项假设断裂准则是三种模式的能量释放率的指数函数。表述为：其中GIC，GIIC，GIIIC 分别为I型，II型，III型断裂模式的临界能量释放率。n1, n2, n3是指数，为常数。六个量都可以通过TBDATA命令定义为温度的函数， 如下所示:三个临界能量释放率不能同时为零。如果其中一个设为零，相应的项被忽略。当指数n1, n2, n3设置为1时，指数准则简化为线性断裂准则。指数准则适用于当三种断裂模式的临界能量释放率明显存在时的三维混合断裂模式。 示例12.6 指数准则输入 g1c=10.0 g2c=20.0 g3c=25.0 n1=2 n2=2 n3=3 TB,CGCR,1,POWERLAWTBDATA,1,g1c,g2c,g3c,n1,n2,n312.1.3.7.用户自定义断裂准则一个自定义的断裂准则可以表述为：f=f(GIC,GIIC,GIIIC,GT,)其中断裂准则是三种断裂模式的能量释放率和材料常数的函数。所有的值通过TBDATA命令输入。你必须提供一个子程序。下面为一个子程序定义线性断裂准则的例子：*deck,user_cgfcrit optimize SUBROUTINE user_cgfcrit (cgi, cid, kct, & nprop, prop, fcscl, & var1, var2, var3, var4)c*cc * primary function:c compute facture criterion for crack growthc user fracture criterion examplec * notice: this routine contains sasi confidential information *c#include impcom.inc#include ansysdef.inccc input argumentsc =c cgi (int,sc , in) CGROW set idc cid (int,sc , in) CINT ID to be usedc kct (int,sc , in) Current crack tip nodec nprop (int,sc , in) number of propertiesc prop (dp ,ar(*), in) property arraycc Output argumentsc =c fcscl (dp, sc , ou) fracture criterionc a return value of one or biggerc indicates fracturecc Misc. argumentsc =c var1 ( , , ) not usedc var2 ( , , ) not usedc var3 ( , , ) not usedc var4 ( , , ) not usedcc*cc * subroutines/functionc * get_cgfpar : API to access fracture datac * wrinqr : ansys standard io functionc * external get_cgfpar external wrinqr integer wrinqrc * argumentc INTEGER cgi, cid, kct, nprop double precision fcscl, & var1, var2, var3, var4 double precision prop(nprop)cc * local variablec integer debugflag, iott integer nn double precision g1c, g2c, g3c, g1, g2, g3 double precision gs(4),da(1)cc * local parameters DOUBLE PRECISION ZERO, ONE parameter (ZERO = 0.0d0, & ONE = 1.0d0)cc*c * initialization fcscl = ZEROc * retrieve energy-release ratesc * for crack cid and crack tip node kctc * gs(1:3) will be returned as G1, G2, G3c * get energy-release rates nn = 3 gs(1:nn) = ZERO call get_cgfpar (GS , cid, kct, 0, nn, gs(1)c * get crack extension nn = 1 da(1) = ZERO call get_cgfpar (DA , cid, kct, 0, nn, da(1)c * energy-release rates g1 =abs(gs(1) g2 =abs(gs(2) g3 =abs(gs(3)c * input property from TBDATA,1,c1,c2,c3 g1c = prop(1) g2c = prop(2) g3c = prop(3)c * linear fracture criterion fcscl = ZERO if (g1c .gt. TINY) fcscl = fcscl + g1/g1c if (g2c .gt. TINY) fcscl = fcscl + g2/g2c if (g3c .gt. TINY) fcscl = fcscl + g3/g3cc * user debug output debugflag = 1 if (debugflag .gt. 0) then iott = wrinqr (WR_OUTPUT) write(iott, 1000) cgi, cid, kct, da(1), fcscl, gs(1:3) 1000 format (5x,user fracture criterion:/ & 5x, crack growth set ID =,i5/ & 5x, crack ID =,i5/ & 5x, crack tip node =,i5/ & 5x, crack extension =,g11.5/ & 5x, calculated fracture parameter =,g11.5/ & 5x, energy-release rates Gs(1:3) =,3g12.5) end if return end12.1.4裂纹扩展模拟例题 该例子采用双悬臂梁,一端有线裂纹。在梁端裂纹上下施加相反方向的相同位移，来张开裂纹，如图所示： 图12.4 双悬臂梁的裂纹扩展下图显示了有限元网格：采用PLANE182单元，打开enhanced strain选项（keyopt（1）=2）来建立实体模型。INTER202单元来建立裂纹路径。假设为平面应变条件。在垂直方向采用6个单元，水平方向200个单元。下图显示了预测的荷载-位移曲线。图12.6 双悬臂梁荷载-位移曲线在裂纹开始扩展前，反力随着位移增大并迅速达到峰值。然后在裂纹扩展的初始阶段，反力迅速下降，随着裂纹增长下降速度变慢。结果和文献结果符合的很好。下面给出最大主应力的云图。图12.7 双悬臂梁的云图下面给出该双悬臂梁的裂纹扩展模拟的输入文件：/BATCH/TITLE, CRACK GROWTH SIMULATION OF A DCB PROBLEM - 2D PLANE STRAIN/PREP7DIS1=0.9DIS2=12.0N1=1000N2=1000N3=10DL=100DH=3A0=30NEL=200NEH=6TOLER=0.1E-5ET,1,182 !* 2D 4-NODE STRUCTURAL SOLID ELEMENTKEYOPT,1,1,2 !* ENHANCE STRAIN FORMULATIONKEYOPT,1,3,2 !* PLANE STRAINET,2,182KEYOPT,2,1,2KEYOPT,2,3,2ET,3,202 !* 2D 4-NODE COHESIVE ZONE ELEMENT!KEYOPT,3,2,2 !* ELEMENT FREE OPTIONKEYOPT,3,3,2 !* PLANE STRAINMP,EX,1,1.353E5 !* E11 = 135.3 GPAMP,EY,1,9.0E3 !* E22 = 9.0 GPAMP,EZ,1,9.0E3 !* E33 = 9.0 GPAMP,GXY,1,5.2E3 !* G12 = 5.2 GPAMP,PRXY,1,0.24MP,PRXZ,1,0.24MP,PRYZ,1,0.46G1C=0.28 !* CRITICAL ENERGY-RELEASE RATEG2C=0.80G3C=0.80TB,CGCR,1,3,LINEAR !* LINEAR FRACTURE CRITERIONTBDATA,1,G1C,G2C,G3C! FE MODELRECTNG,0,DL,DH/2 !* DEFINE AREASRECTNG,0,DL,0,-DH/2LSEL,S,LINE,2,8,2 !* DEFINE LINE DIVISIONLESIZE,ALL,DH/NEHLSEL,INVELESIZE,ALL, , ,NELALLSEL,ALLTYPE,1 !* MESH AREA 2MAT,1LOCAL,11,0,0,0,0ESYS,11AMESH,2CSYS,0TYPE,2 !* MESH AREA 1ESYS,11AMESH,1CSYS,0NSEL,S,LOC,X,A0-TOLER,DLNUMMRG,NODESESLNTYPE,3MAT,5CZMESH,1,Y,0, !* GENERATE INTERFACE ELEMENTSALLSEL,ALLNSEL,S,LOC,X,DL !* APPLY CONSTRAINTSD,ALL,ALLNSEL,ALL!ESEL,S,ENAME,202 !* SELECT INTERFACE ELEMENT TOCM,CPATH,ELEM !* DEFINE CRACK GROWTH PATHNSLENLISTNSEL,S,LOC,X,A0NSEL,R,LOC,Y,0NLISTESLNELISTCM,CRACK1,NODE !* DEFINE CRACK TIP NODE COMPONENTNLISTALLSFINISH/SOLURESC,NONEESEL,S,TYPE,2NSLE,SNSEL,R,LOC,XNSEL,R,LOC,Y,DH/2 !* APPLY DISPLACEMENT LOADING ON TOPD,ALL,UY,DIS1NSEL,ALLESEL,ALLESEL,S,TYPE,1NSLE,SNSEL,R,LOC,XNSEL,R,LOC,Y,-DH/2 !* APPLY DISPLACEMENT LOADING ON BOTTOMD,ALL,UY,-DIS1NSEL,ALLESEL,ALLAUTOTS,ONTIME,1CINT,NEW,1 !* CRACK IDCINT,TYPE,VCCT !* VCCT CALCULATIONCINT