数据处理--文章《按照朱老师修改边框_不对称共线裂纹的应力强度

myEcllipic2.m
for i=1:10
? ? a=2;
? ? b=1;
? ? c=2+0.1*(i-1);
%left crack line%
sqr_alpha_1=-(c-b)/(b+c);
sqr_alpha=(c-b)/(b+c);
sqr_k=(a-b)*(c-b)/((a+b)*(b+c));
g=2/sqrt((a+b)*(b+c));
v0=mfun('EllipticK',sqrt(sqr_k));
v1=mfun('EllipticPi',sqr_alpha,sqrt(sqr_k));
v2=(0.5/((sqr_alpha-1)*(sqr_k-sqr_alpha)))*(sqr_alpha*mfun('EllipticE',sqrt(sqr_k))+(sqr_k-sqr_alpha)*v0+(2*sqr_alpha*sqr_k+2*sqr_alpha-sqr_alpha^2-3*sqr_k)*v1);
J0=g*v0;
J1=-g*b*(v0-2*v1);
J2=g*b^2*(v0-4*v1+4*v2);



%right crack line%



sqr_alpha_1=-c*(b-a)/(a*(b+c));
sqr_alpha=(b-a)/(b+c);
sqr_k=(a-b)*(c-b)/((a+b)*(b+c));
ratio=sqr_alpha/sqr_alpha_1;
g=2/sqrt((a+b)*(b+c));
v0=mfun('EllipticK',sqrt(sqr_k))+(sqr_k-sqr_alpha)*v0+(2*sqr_alpha*sqr_k+2*sqr_alpha-sqr_alpha^2-3*sqr_k)*v1);
S0=g*v0;
S1=g*a*(sqr_alpha_1/sqr_alpha)*(v0-(ratio^2-1)*v1);
S2=g*a^2*(sqr_alpha_1/sqr_alpha)^2*(v0+(ratio^2-1)*v1+(ratio^2-1)^2*v2);


A=[J0 ? J1;S0 S1]
C=-[J2;S2]
inA=inv(A)*C
Ka(i)=(inA(1)+inA(2)*a+a^2)/sqrt((a^2-b^2)*(a+c))%sqrt(2*pi)*
Kb(i)=(inA(1)+inA(2)*b+b^2)/sqrt((a-b)*(a+c)*(b+b))%sqrt(2*pi)*
nKb(i)=(inA(1)+inA(2)*(-b)+b^2)/sqrt((a-(-b))*(-b+c)*(b+b))%sqrt(2*pi)*
Kc(i)=(inA(1)+inA(2)*(-c)+c^2)/sqrt((a-(-c))*(-b+c)*(c+b))%sqrt(2*pi)*
end
j=1:10
plot(j,Ka)
hold on
plot(j,-Kb)
plot(j,nKb)
plot(j,-Kc)

————————————————————————————————————————————
画图设置的备份：myplot.m ?
set(gcf,'color','none') % 去掉灰色背景
xlabel('Length of crack A2B2(alpha=0 b=1 a=2)')
ylabel('SIF/(5.0120D)')
hold on

plot(k1(:,1),k1(:,2),'r-x',3),'b-o');grid on

legend('K1(A1)','K1(B1)');

edit-> copy figure