《含圆孔双边裂纹板的有限元计算》由会员分享，可在线阅读，更多相关《含圆孔双边裂纹板的有限元计算（10页珍藏版）》请在金锄头文库上搜索。
1、 级 专 业 学 号 学生姓名 指导教师 含圆孔双边裂纹板的有限元计算含圆孔双边裂纹板的有限元计算完成日期 摘要摘要：通过MATLAB有限元的分析方法,针对在平面中的受力模型,分析了在静态荷载作用下的应力分布。关键词关键词：MATLAB 有限元 孔边应力 应力强度因子 裂纹引言引言：有限元分析利用数学近似的方法对真实物理系统（几何和载荷工况）进行模拟。利用简单而又相互作用的元素，即单元，就可以用有限数量的未知量去逼近无限未知量的真实系统。有限元法最初被称为矩阵近似方法，应用于航空器的结构强度计算，并由于其方便性、实用性和有效性而引起从事力学研究的科学家的浓厚兴趣。经过短短数十年的努力，随着计算
2、机技术的快速发展和普及，有限元方法迅速从结构工程强度分析计算扩展到几乎所有的科学技术领域，成为一种丰富多彩、应用广泛并且实用高效的数值分析方法。 一、 问题描述含圆孔不对称双边裂纹板受拉伸荷载的解析解有一个圆孔，半径为 1m，弹性模量为 80.851Gpa，泊松比为 0.32，承受均布拉伸荷载，其大小为q=1,x轴两端孔边有长度不等的裂纹，右边裂纹长度L1=4m，左边裂纹长度L2=3m。分析孔边应力分布情况并画出孔边应力分布图，并求解应力强度因子。并给出孔边应力及强度因子表。如图 1 所示：1q1q41L32 L 1R图 1.二、公式推导含圆孔不对称双边裂纹板受拉伸荷载的解析解对于较复杂的带裂
3、纹的孔口问题, 下面利用复变函数方法讨论带不对称双边裂纹的圆孔口问题的解析解。假设在无限大平面中含有一个带裂纹的圆孔洞, 如图 1, 半径为R， 左边裂纹长为c-R，右边裂纹长为d-R。以圆孔中心为原点, 裂纹所在直线为x 轴建立直角坐标系。如图 2、3 所示： 1AC1ii1C1AD1Dyx )(1AAdiRRRiR)(1CCD1DOc图 2 图 3作保角映射如下： )1)(1()1)(1(2)1)(1()1 ()1 (421 42224222gfgfgfRz（1）其中RccRgRddRf2,22222映射后如图 3 所示（d）1, （Ri）D, （-c）-1, （-Ri）D1,1111则该
4、映射将物理平面上带裂纹形孔洞的无限大平面保角映射到数学平面上的单位圆内部。显然，只有一个奇点=，经计算得 （2） 2/ )(,RegfRs其中Res表示留数.又（1）式对求导得（3） )1)(1()1)(1(2)1)(1/()1)(1() 1)(1(2)1)(1(4) 1(21 42224222222 22gfgfgfgfgfR（1）式的共轭为（4） )1)(1()1 () 1(2)1)(1()1 ()1 (4/121 42224222 1gfgfgfR当孔洞边缘及其所带裂纹面上不受力,仅在受到沿y 轴方向的单向拉应力q令时X=Y=0, B=, . YX4q 2qCiB（5） 04q（6） 0
5、2q（7） 00022qqf经计算可知：（8） df idi 00 021 21考虑到（9 0222222220)1)(1() 1)(1(2)1)(1()()1 ()1 ( )1 ( gfgfmgfmgfmH）（10 L4 33 22 10321aaa）422242)1)(1()1 () 1(2)1)(1(gfgfm（11） 应用柯西积分， （7）式左边（12） 0)( 210Hdi 将公式（6）代入（7）式得（13） dqqidf i 22 21 210 0因为=是和的奇点，根据留数理论得 1（14） ),(Re1 21 sdi（15） ),(Re1 21 211 sdidi将（14） ，
6、（15）和（2）代入（13）得（16） 2)(),(Re0gfqRqsqq将（16）代入（5）有（17） 2)( 43gfqRq又（17）对式对求导得（18） 22)( 43 gfqR-q则， wgfqRq w22)( 43（19），在孔边=0，所以孔边应力)()(23Re)(Re42gfqRq（20）)()(23Re)(Re42gfqRq公式中 Re 为取其实部。 应力强度因子（21）)(/ )(2)()()()(lim2200210 0 CiKKKI三、数值计算结果及分析1、在 MATLAB 中建立问题模型，如图 4 所示：-4-3-2-1012345-101图 42、孔边应力分布图如图
7、5 所示：11.522.533.544.555.5-7.5-7-6.5-6-5.5-5-4.5-4-3.5-3-2.5图 5表 1：精确孔边应力分布表-坐标系中A1D1C1CDA1.24221.67152.12214.16114.61175.0410- 3.01 86-7.0000-3.0182-3.0182-7.0000-3.01862、应力强度因子分布如图 6 所示00.511.522.533.5400.511.522.533.54图 6表 2 左边裂缝 L2 保持不变，右边裂缝 L1 变化时，右边裂缝的应力强度因子L100.0010.010.10.30.511.52 KNaN0.140.
8、43981.30632.01312.36632.823.07153.2540表 3 右边裂缝 L1 保持不变，左边裂缝 L2 变化时，右边裂缝的应力强度因子L200.511.522.53K3.29503.33303.40753.49493.58723.68113.7750四、结论对于无限大板带裂纹的圆孔口问题, 本文利用复变函数方法将复杂的孔洞问题简单化以求其解, 并求出了裂纹端点的应力强度因强度因子表达式。该问题的有效解决进一步证明了复变函数方法的巧妙。本文应用等效裂纹的观念计算孔边裂纹的应力，其结果与传统方法的结果一致，因此该方法是可靠的。除此之外，该方法是建立在传统方法之上的，而且在工程
9、断裂问题中操作简便。五、附录程序代码如下： 1、% 原始图像 q=1; L1=4; L2=3; R=1; c=L2+R; d=L1+R;syms a p z w %即为 ,u 为 ,U 为 ,v 为 ,V 为 ,p 为 ,z 为 f=(R2+d2)/2/R/d; g=(R2+c2)/2/R/c; z=0:pi/2000:2*pi; a=cos(z)+i*sin(z); m=sqrt(f2-1)*(1+a).4+2*(f*g+1)*(1-a.2).2+(g2-1)*(1-a).4); n=(f2-1)*(1+a).2+2*(f*g+1)*(a.2+1)+(g2-1)*(1-a).2; s=n./
10、m; w=R/4./a.*(f*(1+a).2+g*(1-a).2+m); plot(w)2、%孔边应力计算 q=1; L1=4; L2=3; R=1; c=L2+R; d=L1+R; f=(R2+d2)/2/R/d; g=(R2+c2)/2/R/c;%求精确角度 w=1; w1=1/189*(40+40*w2-19*w+4*(100-2010*w2-95*w+100*w4- 95*w3)(1/2)/w; w2=1/189*(40+40*w2-19*w-4*(100-2010*w2-95*w+100*w4- 95*w3)(1/2)/w; z1=angle(w1),z2=angle(w2) w=
11、-1; w1=1/189*(40+40*w2-19*w+4*(100-2010*w2-95*w+100*w4- 95*w3)(1/2)/w; w2=1/189*(40+40*w2-19*w-4*(100-2010*w2-95*w+100*w4- 95*w3)(1/2)/w; z3=angle(w1),z4=angle(w2) z=z1:pi/2000:z4; a=cos(z)+i*sin(z); m=sqrt(f2-1)*(1+a).4+2*(f*g+1)*(1-a.2).2+(g2-1)*(1-a).4);n=(f2-1)*(1+a).2+2*(f*g+1)*(a.2+1)+(g2-1)*(
12、1-a).2; wds=R/4*(a.2-1)./a.2).*(f+g+n./m); U=-3*q/4+q*R*(f+g)/2./a.2./wds; xx=4*real(U); figure plot(z,xx) hold on; z=z3+2*pi:pi/2000:z2+2*pi; plot(z,xx)3、%应力强度因子计算 %左边裂缝长度变 clc clear L2=0.001:0.05:3 q=1;L1=4;R=1;c=L2+R;d=L1+R; z=0; a=cos(z)+i*sin(z);f=(R2+d2)/2/R/d; g=(R2+c.2)./(c.*2*R); m=sqrt(f2-
13、1)*(1+a)4+2*(f*g+1)*(1-a2)2+(g.2-1)*(1-a)4); n=(f2-1)*(1+a)2+2*(f*g+1)*(a2+1)+(g.2-1)*(1-a)2; w=R/4/a*(f*(1+a)2+g*(1-a)2+m); wds=R/4*(a2-1)/a2)*(f+g+n./m); uds=-3/4*q*wds+q*R*(f+g)/2/a2; werds =1/2/L1*(L13*L2+L13+L12*L2.2+5*L12+5*L12*L2+2*L1*L2.2+8*L1+8 *L1*L2+4*L2+L2.2+4)./(L2+1)/(L1+2); a=1; K=2*sqrt(pi)*(uds./sqrt(werds) plot(L2,K) hold on%右边裂缝长度变 L1=0.001:0.05:4; L2=3; c=L2+R;d=L1+R; z=0; a=cos(z)+i*sin(z); f=(R2+d.2)./(d.*2*R); g=(R2+c2)/(c*2*R); m=sqrt(f.2-1)*(1+a)4+2*(f*g+1)*(1-a2)2+(g2-1)*(1-a)4); n=(f.2-1)*(1+a)2+2*(f*g+1)*(a2+1)+(g2-1)*(1-a)2; w=
