断裂力学作业解剖.docx
《断裂力学作业解剖.docx》由会员分享，可在线阅读，更多相关《断裂力学作业解剖.docx（17页珍藏版）》请在冰豆网上搜索。
断裂力学作业解剖
理论与应用断裂力学作业
1、解析法求裂尖应力强度因子
如图所示,无限大板,中心穿透裂纹,位于裂纹表面处作用集中法向力P和集中切向力Q,已知：
用极限公式求解裂纹尖端应力强度因子。
解：
已知：
将题目中的带入上式得：
所以，得：
2、Green函数法求裂尖应力强度因子
如图,已知
（平面应力）
（平面应变）
P为作用在裂纹中心处y向对称集中力，L为集中力作用点到裂纹中点的y向距离。
当a＝5mm、15mm、25mm，L＝5mm、15mm、25mm时，求裂纹尖端应力强度因子，并总结规律。
解：
取P＝1000N，泊松比为0.3。
分别计算平面应力和平面应变两种情况下的应力强度因子。
由得
采用复化辛普森数值积分公式。
分析：
复化辛普森数值积分公式需要用到被积函数在积分上界的函数值，由K1的表达式可以看出被积函数在积分上界处的函数值是无穷大。
但是通过被积函数的导数得知当x→a时导数为无穷大，说明被积函数存在垂直渐近线，且渐近线是y＝a。
根据函数的特性（根号下的式子互为倒数），使得在x→a的很小的左临域（a-δ，a）（δ>0）内函数值急剧趋于无穷大。
也就导致了被积函数在左临域（a-δ，a）内函数值趋于无穷大。
因此，可以用f（a-δ）代替f（a）,所得的积分误差应该不会太大。
C＋＋程序如下：
#include
#include
#defineA00.0//积分初始点
#defineP1000//外载
#defineN5000//积分区间等分数！
！
！
N必须为偶数。
#definePSB0.3
doubleA[3]={0.005,0.015,0.025};
doubleL[3]={0.005,0.015,0.025};
intmain（intargc,char*argv[]）
{
doublesun=0.0,y,x,a1,a2,a3,a4,t,buchang,terta;
inti,j=1,k,m;
FILE*fp1,*fp2;
if（（fp1=fopen（"平面应力状态.txt","w"））==NULL）
printf（"无法打开这个文件!
\n"）;
if（（fp2=fopen（"平面应变状态.txt","w"））==NULL）
printf（"无法打开这个文件!
\n"）;
loop1:
printf（"请输入1或2：
1表示平面应力状态，2表示平面应变状态\n"）;
scanf（"%d",&i）;
if（!
（i==1||i==2））
{
printf（"您的输入有误，请重新输入\n"）;
gotoloop1;
}
else
{
if（i==1）
t=（3.0-PSB）/（1.0+PSB）;
if（i==2）
t=3.0-4.0*PSB;
for（k=0;k<3;k++）
{
for（m=0;m<3;m++）
{
buchang=（A[k]-A0）/N;
x=A0+buchang;
terta=0.07*buchang;//临域半径
sun=0.0;
j=1;
while
（1）
{
if（j<=N-1）//进行N－1次循环
{
if（j%2）//奇数项求和
{
a1=sqrt（double（A[k]+x）/double（A[k]-x））+sqrt（double（A[k]-x）/double（A[k]+x））;
a2=（t-1）/（t+1）+（4.0/（t+1.0））*（L[m]*L[m]/double（x*x+L[m]*L[m]））;
a3=3.14159265*（x*x+L[m]*L[m]）;
a4=sqrt（3.14159265*A[k]）;
y=（a1*（double（P））*L[m]*a2/a3）/double（a4）;//被积函数表达式
sun+=4*y;
x+=buchang;
}
if（!
（j%2））//偶数项求和
{
a1=sqrt（double（A[k]+x）/double（A[k]-x））+sqrt（double（A[k]-x）/double（A[k]+x））;
a2=（t-1）/（t+1）+（4.0/（t+1.0））*（L[m]*L[m]/double（x*x+L[m]*L[m]））;
a3=3.14159265*（x*x+L[m]*L[m]）;
a4=sqrt（3.14159265*A[k]）;
y=（a1*（double（P））*L[m]*a2/a3）/double（a4）;//被积函数表达式
sun+=2*y;
x+=buchang;
}
j++;
}
elsebreak;//控制结束循环
}
//求积分下界处被积函数值
a1=sqrt（double（A[k]+A0）/double（A[k]-A0））+sqrt（double（A[k]-A0）/double（A[k]+A0））;
a2=（t-1）/（t+1）+（4.0/（t+1.0））*（L[m]*L[m]/double（A0*A0+L[m]*L[m]））;
a3=3.14159265*（A0*A0+L[m]*L[m]）;
a4=sqrt（3.14159265*A[k]）;
y=（a1*（double（P））*L[m]*a2/a3）/double（a4）;
sun+=y;
//求积分上界处被积函数值
a1=sqrt（double（A[k]+（A[k]-terta））/double（A[k]-（A[k]-terta）））+sqrt（double（A[k]-（A[k]-terta））/double（A[k]+（A[k]-terta）））;
a2=（t-1）/（t+1）+（4.0/（t+1.0））*（L[m]*L[m]/double（（A[k]-terta）*（A[k]-terta）+L[m]*L[m]））;
a3=3.14159265*（（A[k]-terta）*（A[k]-terta）+L[m]*L[m]）;
a4=sqrt（3.14159265*A[k]）;
y=（a1*（double（P））*L[m]*a2/a3）/double（a4）;
sun+=y;
sun*=buchang/3.0;
printf（"%f\n",sun）;
if（i==1）
fprintf（fp1,"%10.5f\n",sun）;
if（i==2）
fprintf（fp2,"%10.5f\n",sun）;
}
}
}
fclose（fp1）;
fclose（fp2）;
return0;
}
将程序算得数值解和将A,L，k值代入
得到的精确解比较得到下面的结果。
（1）平面应力状态下
平
面
应
力
状
态
a值（mm）
L值（mm）
数值解（N/m^3/2）
精确解（N/m^3/2）
5
5
7475.54152
7475.512
15
3999.19194
3999.165
25
2542.78623
2542.768
15
5
4654.26124
4654.256
15
4316.00591
4315.989
25
3502.84414
3502.825
25
5
3586.43118
3586.429
15
3586.21229
3586.204
25
3051.87705
3051.151
表一
结果分析：
由表一数据分析可知：
1、当a固定时，裂纹尖端应力强度因子随L的增加而减小,并且裂纹长度越小，影响越明显。
2、当L固定时，裂纹尖端应力强度因子随a的增加而减小。
L越小，影响越明显。
（2）平面应变状态下
平
面
应
变
状
态
a值（mm）
L值（mm）
数值解（N/m^3/2）
精确解（N/m^3/2）
5
5
7656.88818
7656.859
15
4145.17409
4145.146
25
2639.51097
2639.492
15
5
4682.35481
4682.350
15
4420.70645
4420.689
25
3614.87523
3614.855
25
5
3595.08214
3595.081
15
3638.27907
3638.271
30
3424.26449
3424.251
表二
结果分析：
由表二数据分析得知：
1、当a固定时，裂纹尖端应力强度因子随L的增加而减小,并且裂纹长度越小，影响越明显。
2、当L固定时，裂纹尖端应力强度因子随a的增加而减小。
L越小，影响越明显。
3、在相同情况下，平面应力状态下的应力强度因子小于平面应变状态。
原因是，平面应变屈服区尺寸远比平面应力屈服区尺寸小，平面应变下进入塑性的应力，高于平面应力下进入塑性的应力。
4、利用复化辛普森公式计算得到的数值解和精确解比较接近，具有较高的计算精度。
3、边界配置法编程计算裂纹尖端应力强度因子
如图所示,三点弯曲试件,已知，，，a/W=0.1,0.2,0.3,，，，用边界配置法编程计算裂纹尖端应力强度因子。
习题图
解：
已知Williams应力函数为：
其中：
由于以上三点弯曲试件相对于裂纹面是对称的，故Williams应力函数只需取偶数部分，即：
取2m项，截尾后得：
在试件边界上取m个点，建立2m阶线性方程组，以求解。
取应力函数和其在边界上的外法向导数建立边界条件。
即：
其中，为边界配置点极坐标。
将应力函数改写成如下形式：
是配置点的编号。
其中：
取m＝20，其边界配置如图：
C’D上9个点，B’A上8个点，B’C’上3个点。
其中：
Δ＝＝0.2，
Δ1＝0.25
配置点的坐标为：
上表面段xi＝W–a，，
（=1,3,5,7,9,11,13,15,17）
下表面段xi＝–a，，
（=2,4,6,8,10,12,14,16）
端面段，，
（=18,19,20）
最后求解20阶线性方程组，得到，再由公式
求出应力强度因子。
程序使用LU分解法解线性方程组。
程序如下：
#include
#include
#defineM20//边界点数
#defineB1.0//厚度cm
#defineW20.0//高度cm
#defineP1000.0//力
#defineS14.0*double（W）//试件长
#defineS23.2*double（W）//试件有效长S一撇
#defineN3//需要计算a/W的个数
doubleABW[N]={0.1,0.2,0.3};//a/w值
intmain（intargc,char*argv[]）
{
doublezuobiao[4][M],buchang1,buchang2;
doubleaa[2*M][2*M]
升级会员 