·1514· 工程科学学报，第40卷，第12期 001、002边的中点，节点Q、N分别代表001、002边 由式(4)可得： 的1/4. F1=F6=3.89×103N 设子单元长度为h,三个节点坐标分别为，1= F2=F5=2.83×103N (11) 0,r2=mh,r3=h. F3=F4=-1.06×103N 构建形函数： 进一步求解各圆孔承受应力： W=-0.5n(1-n) F_2.61×103 N2=(1+7)(1-7) (5) 01=06= S0.01×0.007 =5.56×107Pa N3=0.5n(1+7) -1500 式中，N、N2、V3为单元节点形函数，)为局部坐标 02=5= S0.01×0.007=4.04×10'Pa 参数. 13.5 单元坐标变换公式为： =,=3=0.01x0.07=-1.51×10Pa r=N+N52+N53=mh(1-72)+h1+刀)卫 (12) 2 式中，S为内侧圆孔面积的一半 (6) 建立了计算模型后，对模型进行网格划分，由于 可得： 裂纹尖端存在应力集中现象，为保证计算的精确性， 对裂纹尖端进行网格加密处理，如图6所示.施加 -1± 1-8(1-2m) (m-r 约束及载荷前要设立局部坐标系，这里以裂纹尖端 刀=一 2(1-2m) 为圆心，裂纹扩展方向为局部坐标的x轴，垂直裂纹 (7) 扩展方向为局部坐标的y轴，建立局部坐标系.。施 因为-1≤刀≤1，取正号，即： 加载荷及约束，然后得到裂纹尖端有效应力云图，如 -1+1-8(1-2m)(m-方 图7所示. = 2(1-2m) m2 (8) 单元位移函数为： u=N,山1+N242+N343 (9) x=N1+N2x2 N3x3 式中，山1山2山3表示阶段位移，x1、x2x3表示x方向 位移. 单元边长h,将1/2节点向1/4节点位置处移 动，即m=1/4,研究一维问题则x1=0,x2=h/4,x3 =h 图6网格划分图 e器-器2-2+ ax an ax Fig.6 Meshing map (10) 最终输出结果： 由此得应变ε与F成反比，又知应力和应变成 K1.o=5.04×10Pam2, 正比，则应力和F成反比，即将边的中心节点移至1/ Km.o=2.48×10Pam2. 4节点处，奇异单元可反映应力随F的变化，即解决 由式(2)计算得到的解析解为： 了裂纹尖端奇异性的问题 K1.A=5.10×10Pam 因为试件厚度远远小于试件的长和宽，因此数 K1A=2.67×10Pam2. 值模拟计算时可将问题简化为平面问题.建模过程 由此可见K,、K,误差分别为1.18%、7.12%， 边界条件为A、C两点处施加x方向约束，B点施加 计算误差在合理范围内，由此证明用ANSYS计算 y方向约束4（图1），采用PLANE183单元进行建 CTS试件裂纹尖端应力强度因子是可行的. 模计算. 3.2裂纹扩展后应力强度因子的计算 这里以载荷加载角度45°为例进行计算说明， 计算裂纹扩展后的裂纹尖端应力强度因子时将工程科学学报,第 40 卷,第 12 期 OO1 、OO2边的中点,节点 Q、N 分别代表 OO1 、OO2边 的 1 / 4. 设子单元长度为 h,三个节点坐标分别为 r1 = 0,r2 = mh,r3 = h. 构建形函数: N1 = - 0郾 5浊(1 - 浊) N2 = (1 + 浊)(1 - 浊) N3 = 0郾 5浊(1 + 浊 ì î í ïï ïï ) (5) 式中,N1 、N2 、N3 为单元节点形函数,浊 为局部坐标 参数. 单元坐标变换公式为: r = N1 r1 + N2 r2 + N3 r3 = mh(1 - 浊 2 ) + h (1 + 浊)浊 2 (6) 可得: 浊 = - 1 依 1 - 8(1 - 2m) ( m - r ) h 2(1 - 2m) ( m屹 ) 1 2 (7) 因为 - 1臆浊臆1,取正号,即: 浊 = - 1 + 1 - 8(1 - 2m) ( m - r ) h 2(1 - 2m) ( m屹 ) 1 2 (8) 单元位移函数为: u = N1 u1 + N2 u2 + N3 u3 x = N1 x1 + N2 x2 + N3 x { 3 (9) 式中,u1 、u2 、u3表示阶段位移,x1 、x2 、x3 表示 x 方向 位移. 单元边长 h,将 1 / 2 节点向 1 / 4 节点位置处移 动,即 m = 1 / 4,研究一维问题则 x1 = 0,x2 = h / 4,x3 = h 着 = 鄣u 鄣x = 鄣u 鄣浊 鄣浊 鄣x = 1 [ rh - 2浊u2 + ( 1 + 2浊 ) 2 u3 ] (10) 由此得应变 着 与 r成反比,又知应力和应变成 正比,则应力和 r成反比,即将边的中心节点移至 1 / 4 节点处,奇异单元可反映应力随 r的变化,即解决 了裂纹尖端奇异性的问题. 因为试件厚度远远小于试件的长和宽,因此数 值模拟计算时可将问题简化为平面问题. 建模过程 边界条件为 A、C 两点处施加 x 方向约束,B 点施加 y 方向约束[14] (图 1),采用 PLANE183 单元进行建 模计算. 这里以载荷加载角度 45毅为例进行计算说明, 由式(4)可得: F1 = F6 = 3郾 89 伊 10 3 N F2 = F5 = 2郾 83 伊 10 3 N F3 = F4 = - 1郾 06 伊 10 3 ì î í ï ï ï ï N (11) 进一步求解各圆孔承受应力: 滓1 = 滓6 = F1 S = 2郾 61 伊 10 3 0郾 01 伊 0郾 007 = 5郾 56 伊 10 7 Pa 滓2 = 滓5 = F2 S = 1500 0郾 01 伊 0郾 007 = 4郾 04 伊 10 7 Pa 滓3 = 滓4 = F3 S = 13郾 5 0郾 01 伊 0郾 007 = - 1郾 51 伊 10 7 ì î í ï ï ïï ï ï ïï Pa (12) 式中,S 为内侧圆孔面积的一半. 建立了计算模型后,对模型进行网格划分,由于 裂纹尖端存在应力集中现象,为保证计算的精确性, 对裂纹尖端进行网格加密处理,如图 6 所示. 施加 约束及载荷前要设立局部坐标系,这里以裂纹尖端 为圆心,裂纹扩展方向为局部坐标的 x 轴,垂直裂纹 扩展方向为局部坐标的 y 轴,建立局部坐标系. 施 加载荷及约束,然后得到裂纹尖端有效应力云图,如 图 7 所示. 图 6 网格划分图 Fig. 6 Meshing map 最终输出结果: K玉,O = 5郾 04 伊 10 6 Pa·m 1 / 2 , K域,O = 2郾 48 伊 10 6 Pa·m 1 / 2 . 由式(2)计算得到的解析解为: K玉,A = 5郾 10 伊 10 6 Pa·m 1 / 2 , K域,A = 2郾 67 伊 10 6 Pa·m 1 / 2 . 由此可见 K玉、K域 误差分别为 1郾 18% 、7郾 12% , 计算误差在合理范围内,由此证明用 ANSYS 计算 CTS 试件裂纹尖端应力强度因子是可行的. 3郾 2 裂纹扩展后应力强度因子的计算 计算裂纹扩展后的裂纹尖端应力强度因子时将 ·1514·