断裂力学作业解剖文档格式.docx
《断裂力学作业解剖文档格式.docx》由会员分享，可在线阅读，更多相关《断裂力学作业解剖文档格式.docx（17页珍藏版）》请在冰豆网上搜索。
根据函数的特性（根号下的式子互为倒数），使得在x→a的很小的左临域（a-δ，a）（δ>
0）内函数值急剧趋于无穷大。
也就导致了被积函数在左临域（a-δ，a）内函数值趋于无穷大。
因此，可以用f（a-δ）代替f（a）,所得的积分误差应该不会太大。
C＋＋程序如下：
#include<
stdio.h>
math.h>
#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"
printf（"
loop1:
请输入1或2：
1表示平面应力状态，2表示平面应变状态\n"
scanf（"
%d"
&
i）;
if（!
（i==1||i==2））
{
您的输入有误，请重新输入\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<
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））//偶数项求和
sun+=2*y;
}
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]）;
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]）;
sun*=buchang/3.0;
%f\n"
sun）;
if（i==1）
fprintf（fp1,"
%10.5f\n"
fprintf（fp2,"
}
}
fclose（fp1）;
fclose（fp2）;
return0;
}
将程序算得数值解和将A,L，k值代入
得到的精确解比较得到下面的结果。
（1）平面应力状态下
平
面
应
力
状
态
a值（mm）
L值（mm）
数值解（N/m^3/2）
精确解（N/m^3/2）
5
7475.54152
7475.512
15
3999.19194
3999.165
25
2542.78623
2542.768
4654.26124
4654.256
4316.00591
4315.989
3502.84414
3502.825
3586.43118
3586.429
3586.21229
3586.204
3051.87705
3051.151
表一
结果分析：
由表一数据分析可知：
1、当a固定时，裂纹尖端应力强度因子随L的增加而减小,并且裂纹长度越小，影响越明显。
2、当L固定时，裂纹尖端应力强度因子随a的增加而减小。
L越小，影响越明显。
（2）平面应变状态下
变
7656.88818
7656.859
4145.17409
4145.146
2639.51097
2639.492
4682.35481
4682.350
4420.70645
4420.689
3614.87523
3614.855
3595.08214
3595.081
3638.27907
3638.271
30
3424.26449
3424.251
表二
由表二数据分析得知：
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分解法解线性方程组。
程序如下：
#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值
doublezuobiao[4][M],buchang1,buchang2;
doubleaa[2*M][2*M]
升级会员 