计算应力强度因子K

《计算应力强度因子K》由会员分享，可在线阅读，更多相关《计算应力强度因子K（12页珍藏版）》请在金锄头文库上搜索。

1、西北农学院学报用八结点退化三角形等参数奇异性单元计算应力强度因子K,黄松梅刘恩锡(西北农学院西安水利水电科研所)提要本文通过对奇异性的讨论及计算实践指出:如适当地划分网格,用退化三角形奇异单元可以较准确地求得缝顶单元的位移和应力,克 服四边 形奇异单元 贾:一r曲线 在缝顶单元弯曲下降的缺点。推求了用 缝顶单元结点位移计算K:的公式。用这一方法对巳有解析解答的典型问题进行了K;值的计算,计算结果和解析解答良好的吻合。并将这一方法应用于结构形状比较复杂的结构,对于有凹角的裂缝体除缝顶外在凹角处采 用奇异单元可以提高计算K:的精度。在线弹性断裂力学 中计算应力强度因子 是主要任务之一,其中平面张开

2、型的应力强度因 子K:是常 用的。对 于任意结构型式和裂缝 形状的构件,用有限单元法确定应力强度因 子是最为有效的方法。使用 常应变的三角形单元来计算应力强度因子,需要划分尺寸很小的单元,求解上千个联立方程式,且收敛性很慢,常得不到满意的结果。一 近几年来,许多研究者 相继提出了符合缝顶应力、应变分布规律的奇异性单元,用有限单元法计 算K:取得了良好的结果。方法之一是用八结点四边形等参数奇异性单元,即将缝顶单元边界中间结点移至靠近缝顶的边长1/4处,这时缝顶应力、应变将得 到 应有的,一告的奇异性。当四边形的一条边成为一点时就成为退化三角形奇异性单元。这一方法是1973年He nsh eH”等

3、提出的,在这之后一些文献介绍了用四边形奇异单元及退化三角形奇异单元的计算结果,并讨论了奇异性单元的奇异性。但这些文献对四边形和退化三角形两种奇异单元 没有进行深入的比较。本文在讨论四边形 和退化 三角形奇异性的基础上,通过计算实践比较了用四边形和退化三角形奇异性单元计算K;的结果。计算表明 如适当地划分网格用退化三角形奇异性单元,能较准确地求得缝 顶单元的位移和应力,可以克服四边形奇异性单元了:一r曲线在缝顶单元弯曲下降的缺点。这样,就可以不用了:一!曲线外推求K:值,而可以直接用缝顶单元 的结点位移来计算K:。我们推求了用缝顶单元结点位移值求K:的公式。用这一方法计算了紧凑拉伸、双边裂缝、三

4、点弯曲、中心裂缝和单边裂缝等平板裂缝体的K:值,计算结果和已有的解析解答良好的吻合。还讨论了网格尺寸对K:值的影响。并将这一方法用来计算有裂缝的大头坝的应力强 度 因子,除缝顶布置奇异性单元外,凹角处也布置奇异性单元获得了比较好的计算结果。14西北农学院学报1980年第4期一、平面裂缝体缝顶应力奇异性的特点、在弹性力学平面问题中可以用艾雷应力函数必来求解。对于平面裂 缝体1957年Wi l li ams采用下述无穷级数作为应力函数(下面只 列出了 相应于I型裂缝体 的部分)。图必0 O 万Cj音+一”。 sj=1土+1_ 2尸_,J.、。.Tk一COS、卜二一一1夕口- t 乙李+1艺合+一0

5、“根据应力函数可以求得各应力分量 及位移分量。c o a丫,=万 j=1李一1;艺r_;_i 言C1ri言一“+(一1”Co s(言一1)“、1尹、声、沙弋3 2. 盛一了吸、几了, J才、合一,C 。,(音一3,“仃x义= j=1;音一丫止生C;r峭 2人一李+2+(一i),。05(冬一i)夕乙8万古一,c 。,音一3,“合+一,n(音一, “一J一2CJ一20 0 口x丫二刃 j=1一(冬一i)。in(上,3)口2黄松梅、刘恩锡:用八结点退化三角形等参数寄异性单元计算应力强度因子K:1易0 O ZGu=万 Cj rj=1一“+古+一,。0 5音“+音C O s合一2,“(5)土2止22Gv

6、=0 O 刃Cji=1一k一借一一,s,n古0一专,n音一2,0(6)(1)一(6)式中,r、0参见 图1,Cj是常系数与裂缝体的形状和荷载有关,G是剪切应力模量,K=3一4拼(平面应变),或K二(3一拌)/(l+拌)(平面应力),拌是波桑比。对。;j和u,v,当ra(a为裂缝长),取j二1的主项则得:.一一rK一一匹 一, 自口Y=于c o3令(1+s;n令sn等)(7)K-ax x=一灭于井一一 V之万rC Os令(卜S;n令3警) (8)K,口x丫=,二万井 二一-V艺万r03口co ss e厄-“ne s了“o s s e厄一斋杯亨卜2卜1)c os令-=斋必亨卜Zkl)S;令-302

7、(9)(10)382(11)K:和Ci义是的关系是K:=一C:训厄蕊,K:是张开型裂纹体的应力强度因子,它们、的 定K:=1im侧牙丽又a,(x,0)X=+0(1 3)、。Z G一勇 兄爪,二一不一兀二丁斌2万 ,.11i m卫e s二还里2 x二一。训一x(14)从一(7)一(9)可以看出缝顶应力在。为任意值时均具有r一贵的奇异性,从而各应变分量也具有r一于的奇异性。二、退化三角形和四边形奇异单元的奇异性八节点四边形等参数元法 的基本公 式可参见有关著作“,当用等参数元法计算K:了时缝顶单元可划分为矩 形奇异单元如图2,或划分为退化三角形奇异性单元如 图3(a) 助. (时。-一-1召西北农

8、学院学报1,80年第4期76/2 硫图2矩形奇异单元万2日.1- z(a)C产庆混 确一乙一叫卜拓卜乙( b)图3退化三角形奇异单元对退化三角形奇异单元,日uAo落反一=侧丁一刁uBo刁y一侧丁B arsoum证明 了通过缝顶的任一射线都得到+Ar(15)+B,(16)器丹-、。:器兴一+。,(17)(18)(15)一(18)式中A。、A王、B。、B,、C。、C,、D。、D,都是和r无关 的常数。故采用退化三角形奇异性单元 时,应变及应力均具有r一于的奇异性。对于矩形奇异性单元,Bar s。 um”,只证明 了通 过缝顶的二条边界应变具有r一于的奇异性。文献(7讨论过矩形奇异单元内部应变的奇异

9、性,但只涉及应变表达式的 分黄松梅、刘恩锡:用八结点退 化三角形等参数寄异性单元计算应力强度因子K:对母部分,似不完整。下面我们以正方形奇异性单元的对角线为例推导其应变表达式,来讨论单元内的射线是 否具有r一专的奇异性。_彗咭。;了卜士可r1.1一l叠 ,厂!4己l乡/ 多- 。-一-.- 硒-.一,., 口义图4整体坐标(实际单元)图5局部坐标(标准单元)八结点四边形等参数元素由局部座标变换到整体坐标的变换式是(参见图4、5)8 X=刃N(爹,(歹,号)X i二18 Y=刃N(19)专)Y i二1位移的变换式是1二1(若,刀)u:(梦,刀)v;N8刃一一U(20)N8另一一V i二1N;是插

10、值函数,其表达 式是N、(疗,刀)=(i+普普;) (i+刀刀;)一(i一若) (i+,刀.)一(i一刀)(1+苦荟;)荟:,;2/4+(1一若“)(i+刀刀:) (1一彗;“),2/2+(l一刀)(i+若 若*) (i一,*)普;/2(21)以X,=X:=X。=0,X:=L/4,X3=X=X。=L,X。=L/2,Y:二Y:=Y。二0,Y刁二写2,Y。二Y。=Y7=L,Y。二L/4,代入(1)式得X=L /s(i+若)(一刀荟+荟+,+3)(2 2)18西北农学院学报198。年第4期、_、_Y=L/8(1+,)(一普刀+刀+荟+3)在对角线上X二Y,普=刀得x=L/8(i+荟)(3一荟)Y=L

11、/s(i+若),(3一疗)(23)r=记XZ十Y,=了了L8(i+疗)2(3一疗)叼了钊.4 2 侧石(1+幻训不苦一当爹二,时雅可 比 矩阵变换式是(24)、 ,6巧自、 . . . . , .12L8(1一荟)2Y一5分一日J户 刁Y刁叮)二 !专“2一L专(一若 ”L4(i+断)(2一若)阵逆.其11 ,112J一=)8_L(1+荟)(5一3爹) (3一若)I:r122!2(2一萝)(1一若)一(i+参) (i一普)一(1+聋) (1一普2(2一苦)(i一芳)(26)日Uax=11夏刁u 一, ,一二,11, a雪日Ua亏(2 7)乡N;习u,尹b8刃i=1二U一亡S气a一刁(28)了,

12、一叮寸一日六口一l经推导整理得忍x0。84 侧石了丁(5一3扣了百礴3苦(i一荟) (3一荟)u,+2(卜荟)(3彗一s普+i)uZ+3套(i一若)“u,+s(i+普)(i一舀)u。+3普(3一苦)(i+若)u。+2(i+彗)(3彗2一右舀一i)u。+荟(3+4普一3荟2)u,一8(1一若)u。=A。侧丁f(若)(29)费松梅、刘恩锡:用八结 点退 化 三角形等参数奇异性单元计算应力强度因子Kr上 式中A。是与r无关的常数,而f(扣是和r有关的变数,当r趋近于零,彗趋近 于一1,f(荟)铸。,因此。、具有奇异性。但由于f(妇是和r有关的 变数,在对角线方向奇点邻域 的。、并 不与r一专成正比,

13、以数字为例说明如下:当荟=一0.9时,(29 )式第一项是一1。1 05_/宁干,Ul、/T匕当荟=一0.99时,第一项 是一1。244 一Ui 训rL因此就奇异点附近。x是否与r一于成正比的意义 上来说,在对角线方向。x虽有 奇异性,但并不是r一于形式的奇异性。对。,也可作类似的推 导,在除通过 缝顶 的二条边界外单元内 部的其他射线也可得到类似的证明。三、用四边形奇异性单元与退化三角形奇异性单元计算Kl的比较当用有限单元 法计算出裂缝体的应力和位移时,可根据 式(6 )用 位移 法或 据式(2 )用应力法来求K:O=万,如取式。通常是用 裂缝面 的张开位移即口=二的结点位移v来 计 算K,

14、左公当(6)的前四项并利用(11 )则得K,厂丁厂 V:二-二一二丁一,/艺好I万(Zk+2)卜D,r+Zar3/2+D3rZ (3 0)上式中D,=D3=O令、:=老务了平v则得面:=K:+Dr(31)即在r较小 的范围内,可以认为了:和r有近 似的线性关系。因 此可 绘出了:一厂曲线 用外推法求得r=O时的K:值就是了:值。当用0“。时的二Y Y来求K:时可得到类似于式(3 1)的公式面;二K:+Br(32 )一当用矩形奇异性单元 按位移法计算K,时,通常在缝顶单元贾I值 明显弯曲下降如图“所示。即缝顶单元 的位 移与解析解不 相符合,这 是由于矩形的奇异性单元仅在啄条边界上有r一专的奇异

15、性,不保证缝顶的应力和位移在。为任意值时均具有r一十奇异性。这 时就要除去相 应于 缝顶单元 的了,值,用其 他结点位移值求得的面:值 按外推法求K:值。当采 用退化三角形奇异单元 时,如适当划分网格则 可以消除贾,一r曲 线在 攀顶单元处弯 曲下降的现象,如图7所示。即这时 缝顶单元的位移符合解析解,这是脚于退化三角形奇异单元保证了(7)一(9)式在e为任意值时均具有r一于奇异性酵妾求。西北农学院学报19 80年第4期K l 已J叼P悦已喊小,蛇.: 刁。亏权适纷气门几瓜一图6紧凑拉伸 a /w=0.5面,一r关系图用 矩形奇异性单元5 8单元图7紧凑 拉伸 a /w=0.5 面:一r关系

16、图用退化三角形奇异性单元46单元表1给出了a /w=0.5的紧凑拉伸试件用 两种奇异性单元,按位移法和应力法计算K,的结果。位移法 用的是口=时缝顶单元结点的位移按(3 1)计算,应力法用的是e=O缝顶单元的应力按式(3 2)计算。表lw/a=0.5紧凑拉伸试件的K:值奇异性单元边界配置 法结果 单元数K,B侧丽/P用位移法计算结果误 差 %用应力法计算结果误差%退化兰角形10。2010,150。59。8 93,03矩形15 810。2 09。65。912。692 4。2从表1看出用退化三角形缝顶单元(0二0及0二”)的位移及应力和解析 解是基本吻合的,0=O处的应力误 差为3.0 3 %。而用矩 形单元时0=0处的误差高达24.2%,这样用退化三角形奇异性单元的主要优点是可以较准确地求得缝顶单元的应力和位移值。正确地求得缝