斷裂力學(xué)作業(yè)

1、理論與應(yīng)用斷裂力學(xué)作業(yè)1、解析法求裂尖應(yīng)力強(qiáng)度因子如圖所示, 無限大板, 中心穿透裂紋, 位于裂紋表面處作用集中法向力P和集中切向力Q, 已知：用極限公式求解裂紋尖端應(yīng)力強(qiáng)度因子。解：已知： 將題目中的帶入上式得： 所以，得：2、Green函數(shù)法求裂尖應(yīng)力強(qiáng)度因子如圖, 已知 （平面應(yīng)力） （平面應(yīng)變）P為作用在裂紋中心處y向?qū)ΨQ集中力，L為集中力作用點到裂紋中點的y向距離。當(dāng)a5 mm、15 mm、25mm，L5 mm、15 mm、25mm 時，求裂紋尖端應(yīng)力強(qiáng)度因子，并總結(jié)規(guī)律。解：取P1000N，泊松比為0.3。分別計算平面應(yīng)力和平面應(yīng)變兩種情況下的應(yīng)力強(qiáng)度因子。由 得 采用復(fù)化辛普森數(shù)

2、值積分公式。分析：復(fù)化辛普森數(shù)值積分公式需要用到被積函數(shù)在積分上界的函數(shù)值，由K1的表達(dá)式可以看出被積函數(shù)在積分上界處的函數(shù)值是無窮大。但是通過被積函數(shù)的導(dǎo)數(shù)得知當(dāng)xa時導(dǎo)數(shù)為無窮大，說明被積函數(shù)存在垂直漸近線，且漸近線是ya。根據(jù)函數(shù)的特性（根號下的式子互為倒數(shù)），使得在xa的很小的左臨域（a-，a）（>0）內(nèi)函數(shù)值急劇趨于無窮大。也就導(dǎo)致了被積函數(shù)在左臨域（a-，a）內(nèi)函數(shù)值趨于無窮大。因此，可以用f(a-) 代替f(a),所得的積分誤差應(yīng)該不會太大。C程序如下：#include <stdio.h>#include <math.h>#define A0 0.0

3、 /積分初始點#define P 1000 /外載#define N 5000 /積分區(qū)間等分?jǐn)?shù) ！N 必須為偶數(shù)。 #define PSB 0.3double A3=0.005,0.015,0.025;double L3=0.005,0.015,0.025;int main(int argc, char* argv)double sun=0.0,y,x,a1,a2,a3,a4,t,buchang,terta;int i,j=1,k,m;FILE *fp1,*fp2;if(fp1=fopen("平面應(yīng)力狀態(tài).txt","w")=NULL)printf(&

4、quot;無法打開這個文件!n");if(fp2=fopen("平面應(yīng)變狀態(tài).txt","w")=NULL) printf("無法打開這個文件!n");loop1: printf("請輸入1或2：1表示平面應(yīng)力狀態(tài)，2表示平面應(yīng)變狀態(tài)n");scanf("%d",&i); if(!(i=1|i=2)printf("您的輸入有誤，請重新輸入n");goto loop1; else if(i=1) t=(3.0-PSB)/(1.0+PSB); if(i=2) t

5、=3.0-4.0*PSB; for(k=0;k<3;k+) for(m=0;m<3;m+) buchang=(Ak-A0)/N; x=A0+buchang; terta=0.07*buchang; /臨域半徑 sun=0.0; j=1; while(1) if(j<=N-1) /進(jìn)行N1次循環(huán) if(j%2) /奇數(shù)項求和 a1=sqrt(double(Ak+x)/double(Ak-x)+sqrt(double(Ak-x)/double(Ak+x); a2=(t-1)/(t+1)+(4.0/(t+1.0)*(Lm*Lm/double(x*x+Lm*Lm); a3=3.141

6、59265*(x*x+Lm*Lm); a4=sqrt(3.14159265*Ak); y=(a1*(double(P)*Lm*a2/a3)/double(a4); /被積函數(shù)表達(dá)式 sun+=4*y; x+=buchang; if(!(j%2) /偶數(shù)項求和 a1=sqrt(double(Ak+x)/double(Ak-x)+sqrt(double(Ak-x)/double(Ak+x); a2=(t-1)/(t+1)+(4.0/(t+1.0)*(Lm*Lm/double(x*x+Lm*Lm); a3=3.14159265*(x*x+Lm*Lm); a4=sqrt(3.14159265*Ak);

7、 y=(a1*(double(P)*Lm*a2/a3)/double(a4); /被積函數(shù)表達(dá)式 sun+=2*y; x+=buchang; j+; else break; /控制結(jié)束循環(huán) /求積分下界處被積函數(shù)值 a1=sqrt(double(Ak+A0)/double(Ak-A0)+sqrt(double(Ak-A0)/double(Ak+A0); a2=(t-1)/(t+1)+(4.0/(t+1.0)*(Lm*Lm/double(A0*A0+Lm*Lm); a3=3.14159265*(A0*A0+Lm*Lm); a4=sqrt(3.14159265*Ak); y=(a1*(double

8、(P)*Lm*a2/a3)/double(a4); sun+=y; /求積分上界處被積函數(shù)值 a1=sqrt(double(Ak+(Ak-terta)/double(Ak-(Ak-terta)+sqrt(double(Ak-(Ak-terta)/double(Ak+(Ak-terta); a2=(t-1)/(t+1)+(4.0/(t+1.0)*(Lm*Lm/double(Ak-terta)*(Ak-terta)+Lm*Lm); a3=3.14159265*(Ak-terta)*(Ak-terta)+Lm*Lm); a4=sqrt(3.14159265*Ak); y=(a1*(double(P)

9、*Lm*a2/a3)/double(a4); sun+=y; sun*=buchang/3.0; printf("%fn",sun); if(i=1) fprintf(fp1,"%10.5fn ",sun); if(i=2) fprintf(fp2,"%10.5fn ",sun); fclose(fp1);fclose(fp2);return 0;將程序算得數(shù)值解和將A,L，k值代入得到的精確解比較得到下面的結(jié)果。（1）平面應(yīng)力狀態(tài)下平面應(yīng)力狀態(tài)a值(mm)L值(mm)數(shù)值解(N/m3/2)精確解(N/m3/2)557475.5415

10、27475.512153999.191943999.165252542.786232542.7681554654.261244654.256154316.005914315.989253502.844143502.8252553586.431183586.429153586.212293586.204253051.877053051.151表一結(jié)果分析：由表一數(shù)據(jù)分析可知：1、當(dāng)a固定時，裂紋尖端應(yīng)力強(qiáng)度因子隨L的增加而減小,并且裂紋長度越小，影響越明顯。2、當(dāng)L固定時，裂紋尖端應(yīng)力強(qiáng)度因子隨a的增加而減小。L越小，影響越明顯。（2）平面應(yīng)變狀態(tài)下 平面應(yīng)變狀態(tài)a值(mm)L值(mm)數(shù)值解(

11、N/m3/2)精確解(N/m3/2)557656.888187656.859154145.174094145.146252639.510972639.4921554682.354814682.350154420.706454420.689253614.875233614.8552553595.082143595.081153638.279073638.271303424.264493424.251表二結(jié)果分析：由表二數(shù)據(jù)分析得知：1、當(dāng)a固定時，裂紋尖端應(yīng)力強(qiáng)度因子隨L的增加而減小,并且裂紋長度越小，影響越明顯。2、當(dāng)L固定時，裂紋尖端應(yīng)力強(qiáng)度因子隨a的增加而減小。L越小，影響越明顯。3、在相

12、同情況下，平面應(yīng)力狀態(tài)下的應(yīng)力強(qiáng)度因子小于平面應(yīng)變狀態(tài)。原因是，平面應(yīng)變屈服區(qū)尺寸遠(yuǎn)比平面應(yīng)力屈服區(qū)尺寸小，平面應(yīng)變下進(jìn)入塑性的應(yīng)力，高于平面應(yīng)力下進(jìn)入塑性的應(yīng)力。4、利用復(fù)化辛普森公式計算得到的數(shù)值解和精確解比較接近，具有較高的計算精度。3、邊界配置法編程計算裂紋尖端應(yīng)力強(qiáng)度因子如圖所示,三點彎曲試件,已知，a/W =0.1,0.2,0.3,，用邊界配置法編程計算裂紋尖端應(yīng)力強(qiáng)度因子。習(xí)題圖解：已知Williams應(yīng)力函數(shù)為： 其中： 由于以上三點彎曲試件相對于裂紋面是對稱的，故Williams應(yīng)力函數(shù)只需取偶數(shù)部分，即：取2m項，截尾后得：在試件邊界上取m個點，建立2m階線性方程組，以求解

13、。取應(yīng)力函數(shù)和其在邊界上的外法向?qū)?shù)建立邊界條件。即： 其中，為邊界配置點極坐標(biāo)。將應(yīng)力函數(shù)改寫成如下形式： 是配置點的編號。其中：取m20，其邊界配置如圖：CD上9個點，BA上8個點，BC上3個點。其中：0.2，10.25配置點的坐標(biāo)為：上表面段 xiW a，(=1,3,5,7,9,11,13,15,17)下表面段 xi a，(=2,4,6,8,10,12,14,16)端面段 ，(=18,19,20)最后求解20階線性方程組，得到，再由公式求出應(yīng)力強(qiáng)度因子。程序使用LU分解法解線性方程組。程序如下：#include <stdio.h>#include <math.h>

14、#define M 20 /邊界點數(shù)#define B 1.0 /厚度 cm#define W 20.0 /高度 cm#define P 1000.0 /力#define S1 4.0*double(W) /試件長#define S2 3.2*double(W) /試件有效長 S一撇#define N 3 / 需要計算a/W的個數(shù)double ABWN=0.1,0.2,0.3; / a/w值int main(int argc, char* argv) double zuobiao4M,buchang1,buchang2; double aa2*M2*M,bb2*M,lu2*M2*M; /定義系

15、數(shù)矩陣 double x2*M,y2*M; /存儲結(jié)果 double aN,sun; /定義裂紋長度數(shù)組 double oumikaM; /邊界外法線矢量與x軸夾角?；《?int i,j,k,m,t=0; FILE *fp1,*fp2,*fp3,*fp4; buchang1=S2/(M-4.0); / 上下配置點間距 buchang2=0.25*W; /左斷面配置點間距 for(i=0;i<N;i+) /計算裂紋長度 ai=W*ABWi; if(fp1=fopen("k.txt","w")=NULL) printf("不能打開這個文件！n

16、"); for(t=0;t<N;t+) for(j=0,i=0;i<M/2-1;j+=2,i+) /計算上邊界配置點x坐標(biāo) zuobiao0j=W-at; for(j=0,i=0;i<M/2-1;j+=2,i+) /計算上邊界配置點y坐標(biāo) zuobiao1j=(double(j)/2.0)*buchang1; oumikaj=0; for(j=1,i=0;i<M/2-2;j+=2,i+) /計算下邊界配置點x坐標(biāo) zuobiao0j=-at; for(j=1,i=0;i<M/2-2;j+=2,i+) /計算下邊界配置點y坐標(biāo) zuobiao1j=(dou

17、ble(j+1)/2.0)*buchang1; oumikaj=3.14159265; for(j=M-3;j<M;j+) /計算左邊界配置點x坐標(biāo) zuobiao0j=(double(W)/2.0-at)+(j+1)-(M-1)*buchang2; for(j=M-3;j<M;j+) /計算左邊界配置點y坐標(biāo) zuobiao1j=S2/2.0; oumikaj=3.14159265/2.0; for(j=0;j<M;j+) /轉(zhuǎn)化為極坐標(biāo) zuobiao2j=sqrt(pow(double(zuobiao0j),2.0)+pow(double(zuobiao1j),2.0)

18、; / r坐標(biāo) for(j=0,i=0;i<M/2-1;j+=2,i+) /計算上邊界配置點cita坐標(biāo) zuobiao3j=atan(double(zuobiao1j)/double(zuobiao0j); for(j=1,i=0;i<M/2-2;j+=2,i+) /計算下邊界配置點cita坐標(biāo) zuobiao3j=atan(double(zuobiao1j)/double(zuobiao0j)+3.14159265; for(j=M-3;j<M;j+) /計算左邊界配置點cita坐標(biāo) zuobiao3j=atan(double(zuobiao1j)/double(zuob

19、iao0j); for(i=0;i<M;i+) /求方程的系數(shù) for(j=0;j<2*M;j+) aaij=pow(double(zuobiao2i),1.0+double(j+1.0)/2.0)*(double(j+1.0)/2.0+pow(-1.0,j+1.0)*cos(1.0+double(j+1.0)/2.0)*zuobiao3i)/(1.0+double(j+1.0)/2.0)-cos(double(j+1.0)/2.0-1.0)*zuobiao3i); for(i=0;i<M;i+) /求方程的系數(shù) for(j=0;j<2*M;j+) aaM+ij=pow

20、(double(zuobiao2i),double(j+1.0)/2.0)*(-cos(double(j+1.0)/2.0*zuobiao3i-oumikai)-double(j+1.0)/2.0*cos(double(j+1.0)/2.0-2.0)*zuobiao3i+oumikai)+(double(j+1.0)/2.0+pow(-1.0,j+1.0)*cos(double(j+1.0)/2.0*zuobiao3i+oumikai); for(i=0,j=0;j<M/2-1;i+=2,j+) /求方程常數(shù)向量 bbi=double(P)*(S1-(i+1-1)*S2/(double(

21、M)-4.0)/(4.0*double(B); for(i=1,j=0;j<M/2-2;i+=2,j+) /求方程常數(shù)向量 bbi=0.0; for(i=M-3;i<M;i+) /求方程常數(shù)向量 bbi=-double(P)/double(B)*(double(S1)/2.0-double(S2)/2.0)*(pow(double(i+1.0-M+3.0)/4.0,3.0)-3.0/2.0*pow(double(i+1.0-M+3.0)/4.0,2.0); for(i=M;i<(2*M-3);i+) /求方程常數(shù)向量bbi=0.0;for(i=2*M-3;i<2*M;i

22、+) /求方程常數(shù)向量 bbi=double(P)/double(B)*(pow(double(i+1.0-M+3.0)/4.0,3.0)-3.0/2.0*pow(double(i+1.0-M+3.0)/4.0,2.0);/* for(i=0;i<2*M;i+)for(j=0;j<2*M;j+)aaij*=P*double(W)/B;*/for(j=0;j<2*M;j+) /LU分解，確定L ,U 矩陣 lu0j=aa0j; for(j=1;j<2*M;j+) luj0=aaj0/lu00; for(k=1;k<2*M;k+) for(j=k;j<2*M;j

23、+) sun=0.0; for(m=0;m<k;m+) sun+=lukm*lumj; lukj=aakj-sun; for(i=k+1;i<2*M;i+) sun=0.0; for(m=0;m<k;m+) sun+=luim*lumk; luik=(aaik-sun)/lukk; y0=bb0; for(k=1;k<2*M;k+) /確定y sun=0.0; for(j=0;j<k;j+) sun+=lukj*yj; yk=bbk-sun; x2*M-1=y2*M-1/lu2*M-12*M-1; /求解最終結(jié)果 for(k=2*M-2;k>=0;k-) s

24、un=0.0; for(j=k+1;j<2*M;j+) sun+=lukj*xj; xk=(yk-sun)/lukk; if(t=0) printf("a/W=0.1時， D1=%f ,K1=%f n",x0,-sqrt(2*3.14159265)*x0); if(t=1) printf("a/W=0.2時， D1=%f ,K1=%f n",x0,-sqrt(2*3.14159265)*x0); if(t=2) printf("a/W=0.3時， D1=%f ,K1=%f n",x0,-sqrt(2*3.14159265)*x0); fprintf(fp1,"D1=%5.5f ,K1=%