含裂纹的椭圆孔口问题的应力分析.doc

含裂纹的椭圆孔口问题的应力分析级 专 业 学 号 学生姓名 指导教师 完成日期 带裂纹的椭圆孔口问题的应力分析 有一个椭圆孔，长半轴为a=1m，短半轴为b=1m，弹性模量为80.851Gpa，泊松比为0.32，承受均布拉伸荷载，其大小为q=1Mpa，长轴一端的孔边有裂纹，裂纹长度为L=2m，如图1。求解孔边应力分布图及应力强度因子，并给出当L=1,2,3时强度因子图。 （a）含裂纹的椭圆孔口 (b) 单位圆盘图1 含裂纹的椭圆孔口的无限大平面到单位圆盘的映射1 引言断裂现象与材料和结构中的孔洞、缺口或裂纹等缺陷密切相关, 这是因为缺陷附近的应力集中明显。利用复变函数方法，通过保角映射可处理带裂纹的椭圆孔洞的平面弹性问题, 在各种平面问题中, 孔口问题最能显示复变函数解法的优越性, 有些比较复杂的孔口问题, 如果不用这种方法, 几乎就无法求解。 现在从创立平面弹性的复变方法开始,假定保角映射函数 （1） （2）其中R 以及与均为常数. 在解方程时, 常遇到积分 为计算这些积分, 往往要求是简单的有理函数。 这样, 的解析延拓性质简单, 容易处理。相反,不是简单的有理函数, 则上述些函数很复杂, 其解析延拓性质很难讨论。本文解决了带单边裂纹的椭圆孔口的有关问题，求解了孔边应力分布图，并应用应力强度因子的解析解求解应力强度因子，计算了有关两互相垂直的裂纹等问题，同时分析了有关裂纹问题的若干结果。2 带裂纹的椭圆孔口问题的复变方法对于较复杂的带裂纹的孔口问题,下面利用复变函数方法讨论带裂纹的椭圆孔口问题的解析解。假设在无限大平面中含有一个带裂纹的椭圆孔洞,如图1（其程序见附录1）,长半轴为a,短半轴为b,裂纹长为c-a。以椭圆孔中心为原点,裂纹所在直线为x 轴建立直角坐标系。作保角映射如下 （3）其中 则该映射将物理平面上带裂纹形孔洞的无限大平面保角映射到数学平面上的单位圆内部，如图1所示。其中 A= 0.1111+0.9938i B= -0.6000+0.8000i C=-0.6000-0.8000iD= 0.1111-0.9938i对应角 ，（其程序见附录2）当1时，（）a，因此1是（）的可去奇点,该映射函数只可能有两个极点 其中1,1且=1。经计算得其中Res表示留数。由式（3）可得（）的导数 （4）当孔洞边缘及其所带裂纹面上不受力,仅受沿y 轴方向的单向拉应力q时B=， ， X=Y=0 （5） （6）经计算可知 （7）可求得 （8）则 （9）孔边应力表达式为 （10）由式（11）可绘出孔边应力与角度的关系图，如图2所示(程序见附录3)。从图2可以看出随着的增加孔边应力逐渐减小，当增大到一定值时孔边应力变为负值，表明应力由拉应力变为了压应力。在时，即椭圆的（-a,0）点取得正的最大值3.11；在=时，即椭圆的裂缝左端点处取得应力的最小值（最大压应力）为-0.76。 图2 孔边应力图3 关于裂纹问题的若干结果强度因子表达式为 （11）（1） 在公式（12）中,b0时,，则该问题变成了只含一条裂纹的薄板问题，即Griffith的裂纹问题，其强度因子随a、c的变化如图3（程序见附录4（1）。由图3可知，对于给定的c值，强度因子K随a值的增大而增大；对于同一个a值，强度因子K随c值的增大而增大。 图3 含一条裂纹薄板的强度因子随a、c的变化(2) 在式（12）中,当ba时,。则该问题转变为带裂纹的圆形孔洞中裂纹端点应力强度因子的问题。它的大小与水平裂纹的长度c，和圆孔半径a有关。为了进一步研究应力强度因子随参数a、c的变化规律, 令q=1时绘制曲线如图4所示（程序见附录4（2）。由图4可知，对于给定的c值，圆形孔洞裂纹端点的强度因子K随a值的增大而增大；对于同一个a值，圆形孔洞裂纹端点的强度因子K随c值增大而增大。(3) 在式（12）中,当ba0时,。则该问题同样变成了只含一条直线裂纹的薄板问题，即Griffith 9, 10 的裂纹问题，其强度因子随c的变化如图5所示（程序见附表4（3）。由图5可知，当ba0，时,含一条直线裂纹薄板的强度因子K随c值的增大而增大。图4 圆形孔洞裂纹端点强度因子随a、c的变化图5 含一条直线裂纹的薄板的强度因子随c的变化(4) 当a0时，该问题转变为两条相互垂直裂纹的问题。裂纹形状如图6所示图6 两相互垂直的裂纹则强度因子表达式变为 （12）当薄板无穷远处受到沿y轴方向的单向拉应力q时,裂纹端点c处的应力强度因子的大小不仅与水平裂纹的长度c有关,而且与竖直裂纹的长度2b有关。为了进一步研究应力强度因子随参数b、c的变化规律, 令q=1时绘制曲线如图4所示（其程序代码见4（4）。图7 圆形孔洞裂纹端点强度因子随b、c的变化由图7可知，对于给定的c值，强度因子K随a值的增大而增大，当a值增大到一定值之后，K趋于一个常数；对于同一个a值，强度因子K随c值增大而增大。4 结论（1） 通过对带裂纹椭圆孔的应力分析，可以得到，椭圆孔离裂缝较远处的端点比裂缝处端点应力小，而且远端为正，裂缝处为负。（2） 通过对带裂纹椭圆孔端点处强度因子的分析，可以得到：当b0或ba0时，为水平裂缝，裂缝端点强度因子随裂缝长度的增大而变大；当ba时，为带裂缝的圆孔，裂缝端点强度因子随裂缝长度和圆孔半径的增大而变大；当a0时，为带T形裂缝，水平裂缝端点强度因子随竖直裂缝长度增大而变大，最后趋于稳定，水平裂缝长度增大，强度因子也增大。附录：程序1 通过变换函数画出带缝椭圆孔原图像clcclearsyms x la=1;b=1;c=3;k=(c-a-b+sqrt(c2-a2+b2)/(c+a+b+sqrt(c2-a2+b2);for nm=1:101l=(nm-1)*pi/100;x=cos(l)+i*sin(l);wx=(a+b)/2*(1+sqrt(1+x)/(1-x)2+k2)/(1-sqrt(1+x)/(1-x)2+k2)+(a-b)/2*(1-sqrt(1+x)/(1-x)2+k2)/(1+sqrt(1+x)/(1-x)2+k2);n1(nm)=real(wx);n2(nm)=imag(wx);endplot(n1,abs(n2)hold on;plot(n1,-abs(n2)hold off; 2 求A、B点即椭圆上边界和右边界处对应转换后函数图像的角clcclearsyms x a=1;b=1;c=3k=(c-a-b+sqrt(c2-a2+b2)/(c+a+b+sqrt(c2-a2+b2)%solve(z-(a+b)/2*(1+sqrt(1+x)/(1-x)2+k2)/(1-sqrt(1+x)/(1-x)2+k2)+(a-b)/2*(1-sqrt(1+x)/(1-x)2+k2)/(1+sqrt(1+x)/(1-x)2+k2),x)%z=1;%z=i;k1=(-2*z2*k2+2*a2*k2-4*z*a+a2*k4-4*k2*b2+2*z*k4*a+z2*k4+4*(b2*z2-b2*a2+b4)(1/2)/(-4*b2+4*a2+z2*k4+a2*k4-4*k2*b2+2*z*k4*a+4*a2*k2+4*z*k2*a)+2/(-4*b2+4*a2+z2*k4+a2*k4-4*k2*b2+2*z*k4*a+4*a2*k2+4*z*k2*a)*(z4*k4-z4*k6-5*a4*k4-a4*k6+4*z2*a2-8*b4*k2+4*b2*z2+4*b2*a2-8*a4*k2-4*a4+16*z*a*k2*b2-12*z2*k4*a2+4*z3*k2*a-6*z2*k6*a2+6*z2*k4*b2-4*z3*k6*a-12*a3*k2*z+6*a2*k4*b2-4*a3*k6*z-14*z*a3*k4-2*z3*a*k4+16*b2*a2*k2+12*b2*z*k4*a+4*z*k4*a*(b2*z2-b2*a2+b4)(1/2)-4*z2*k2*(b2*z2-b2*a2+b4)(1/2)+2*z2*k4*(b2*z2-b2*a2+b4)(1/2)+2*a2*k4*(b2*z2-b2*a2+b4)(1/2)+4*a2*k2*(b2*z2-b2*a2+b4)(1/2)-8*z*a*(b2*z2-b2*a2+b4)(1/2)-8*k2*b2*(b2*z2-b2*a2+b4)(1/2)(1/2)k2=(-2*z2*k2+2*a2*k2-4*z*a+a2*k4-4*k2*b2+2*z*k4*a+z2*k4+4*(b2*z2-b2*a2+b4)(1/2)/(-4*b2+4*a2+z2*k4+a2*k4-4*k2*b2+2*z*k4*a+4*a2*k2+4*z*k2*a)-2/(-4*b2+4*a2+z2*k4+a2*k4-4*k2*b2+2*z*k4*a+4*a2*k2+4*z*k2*a)*(z4*k4-z4*k6-5*a4*k4-a4*k6+4*z2*a2-8*b4*k2+4*b2*z2+4*b2*a2-8*a4*k2-4*a4+16*z*a*k2*b2-12*z2*k4*a2+4*z3*k2*a-6*z2*k6*a2+6*z2*k4*b2-4*z3*k6*a-12*a3*k2*z+6*a2*k4*b2-4*a3*k6*z-14*z*a3*k4-2*z3*a*k4+16*b2*a2*k2+12*b2*z*k4*a+4*z*k4*a*(b2*z2-b2*a2+b4)(1/2)-4*z2*k2*(b2*z2-b2*a2+b4)(1/2)+2*z2*k4*(b2*z2-b2*a2+b4)(1/2)+2*a2*k4*(b2*z2-b2*a2+b4)(1/2)+4*a2*k2*(b2*z2-b2*a2+b4)(1/2)-8*z*a*(b2*z2-b2*a2+b4)(1/2)-8*k2*b2*(b2*z2-b2*a2+b4)(1/2)(1/2)k3=(-2*z2*k2+2*a2*k2-4*z*a+a2*k4-4*k2*b2+2*z*k4*a+z2*k4-4*(b2*z2-b2*a2+b4)(1/2)/(-4*b2+4*a2+z2*k4+a2*k4-4*k2*b2+2*z*k4*a+4*a2*k2+4*z*k2*a)+2/(-4*b2+4*a2+z2*k4+a2*k4-4*k2*b2+2*z*k4*a+4*a2*k2+4*z*k2*a)*(z4*k4-z4*k6-5*a4*k4-a4*k6+4*z2*a2-8*b4*k2+4*b2*z2+4*b2*a2-8*a4*k2-4*a4+16*z*a*k2*b2-12*z2*k4*a2+4*z3*k2*a-6*z2*k6*a2+6*z2*k4*b2-4*z3*k6*a-12*a3*k2*z+6*a2*k4*b2-4*a3*k6*z-14*z*a3*k4-2*z3*a*k4+16*b2*a2*k2+12*b2*z*k4*a-4*z*k4*a*(b2*z2-b2*a2+b4)(1/2)+4*z2*k2*(b2*z2-b2*a2+b4)(1/2)-2*z2*k4*(b2*z2-b2*a2+b4)(1/2)-2*a2*k4*(b2*z2-b2*a2+b4)(1/2)-4*a2*k2*(b2*z2-b2*a2+b4)(1/2)+8*z*a*(b2*z2-b2*a2+b4)(1/2)+8*k2*b2*(b2*z2-b2*a2+b4)(1/2)(1/2)k4=(-2*z2*k2+2*a2*k2-4*z*a+a2*k4-4*k2*b2+2*z*k4*a+z2*k4-4*(b2*z2-b2*a2+b4)(1/2)/(-4*b2+4*a2+z2*k4+a2*k4-4*k2*b2+2*z*k4*a+4*a2*k2+4*z*k2*a)-2/(-4*b2+4*a2+z2*k4+a2*k4-4*k2*b2+2*z*k4*a+4*a2*k2+4*z*k2*a)*(z4*k4-z4*k6-5*a4*k4-a4*k6+4*z2*a2-8*b4*k2+4*b2*z2+4*b2*a2-8*a4*k2-4*a4+16*z*a*k2*b2-12*z2*k4*a2+4*z3*k2*a-6*z2*k6*a2+6*z2*k4*b2-4*z3*k6*a-12*a3*k2*z+6*a2*k4*b2-4*a3*k6*z-14*z*a3*k4-2*z3*a*k4+16*b2*a2*k2+12*b2*z*k4*a-4*z*k4*a*(b2*z2-b2*a2+b4)(1/2)+4*z2*k2*(b2*z2-b2*a2+b4)(1/2)-2*z2*k4*(b2*z2-b2*a2+b4)(1/2)-2*a2*k4*(b2*z2-b2*a2+b4)(1/2)-4*a2*k2*(b2*z2-b2*a2+b4)(1/2)+8*z*a*(b2*z2-b2*a2+b4)(1/2)+8*k2*b2*(b2*z2-b2*a2+b4)(1/2)(1/2) m=angle(k1) %m =2.2143 对应B点即椭圆左边界 %m =1.4595 对应A点即椭圆上边界3 求孔边应力clcclearsyms x y p la=1;b=1;c=3;q=1;k=(c-a-b+sqrt(c2-a2+b2)/(c+a+b+sqrt(c2-a2+b2);% n1=1 n2=2 为两个极点 Res1=Res（，1） Res2=Res（，2） 为两个极点处的留数n1=-(1+sqrt(1-k2)/(1-sqrt(1-k2);n2=-(1-sqrt(1-k2)/(1+sqrt(1-k2);Res1=2*(a+b)/sqrt(1-k2)/(1-sqrt(1-k2);Res2=-2*(a+b)/sqrt(1-k2)/(1+sqrt(1+k2);l=0.000000000001:0.01: 2.2143x=cos(l)+i*sin(l)nn=(1+x)./(1-x).2+k2n=real(nn)% wdx 为变换函数（）的导数wdx=2*(a+b)*(1+x)./(1-x).3./1-sqrt(n).2./(sqrt(n)-2*(a-b)*(1+x)./(1-x).3./1+sqrt(n).2./(sqrt(n)fdx=-(q/4.*wdx-q*Res1/2./(x-n1).2+q*Res2/2./(x-n2).2) % fdx 为（）的导数Fx=fdx./wdx % Fdx 为（）ff=4*real(Fx) % ff 为孔边应力plot(l,ff)4 有关应力强度因子的计算(1)clcclearsyms a c %a,c符号化 q=1;a=0:1/10:16 %取做a为自变量，范围是0到16，步长为1/10f=q*sqrt(pi*(a+c)./2) %f为强度因子的表达式% 给定一个c值可求出一条f随a变化的曲线,共绘制4条曲线c=1; % 给定一个c值 f1=eval(f) %求出c=1时函数f和自变量a的关系c=2;f2=eval(f)c=3;f3=eval(f)c=4;f4=eval(f)plot(a,f1)hold on;plot(a,f2)hold on;plot(a,f3)hold on;plot(a,f4)hold off;（2）clcclearsyms a b c% ba，k(c-a)/(c+a) 时转化为带裂纹的圆形孔洞问题中裂纹端点的应力强度因子问题syms a b c %a,b,c符号化,a为椭圆长半轴，c为裂纹长度q=1; %取单向拉应力q为1 a=0:1/5:16 %取做a为自变量，范围是0到16，步长为1/5b=c+a %b为裂纹末端的