2023年7月3日发(作者：)
clear
clc
ratio = 1.624;% for Mg
crss_basal = 2.1;% CRSS for basal in Mg
crss_pyr = 40;% CRSS for
crss_ext = 11.5;% CRSS for extension twin in Mg
crss_comp = 125;% CRSS for compression twin in Mg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
basal = [0,0,0,1];%basal plane
pyr1 = [-1,0,1,-1];%pyramidal I plane
pyr2 = [-1,-1,2,-2];%pyramidal II plane
ext = [1,0,-1,2];%extension twin
comp = [-1,0,1,-1];%compression twin
basal_dir = [1,1,-2,0];%basal Burgers vector
pyr1_dir = [1,1,-2,-3];%pyramidal
pyr2_dir = [1,1,-2,-3];%pyramidal
ext_dir = [1,0,-1,-1];%extension twinning direction (yita 1)
comp_dir = [1,0,-1,-2];%compression twinning direction (yita 1)
%four index
alpha = linspace(0,pi/6,31);
beta = linspace(0,pi/2,91);
load_dir = zeros(length(alpha),length(beta),4);
for i = 1: length(alpha)
for j = 1:length (beta)
load_dir(i,j,:) = [1,2/(sqrt(3)*cot(alpha(i))-1), ...
-(sqrt(3)*cot(alpha(i))+1)/(sqrt(3)*cot(alpha(i))-1), ...
3*tan(beta(j))/ratio/(sqrt(3)*cos(alpha(i))-sin(alpha(i)))];
end
end
m = zeros(length(alpha),length(beta),5);
for i = 1: length(alpha)
for j = 1:length (beta)
m(i,j,1) = Ang_d(load_dir(i,j,:),basal_dir,ratio)* ...
Ang_d(load_dir(i,j,:),Norm(basal),ratio);
m(i,j,2) = Ang_d(load_dir(i,j,:),pyr1_dir,ratio)*...
Ang_d(load_dir(i,j,:),Norm(pyr1),ratio);
m(i,j,3) = Ang_d(load_dir(i,j,:),pyr2_dir,ratio)*...
Ang_d(load_dir(i,j,:),Norm(pyr2),ratio);
m(i,j,4) = Ang_d(load_dir(i,j,:),ext_dir,ratio)*...
Ang_d(load_dir(i,j,:),Norm(ext),ratio);
m(i,j,5) = Ang_d(load_dir(i,j,:),comp_dir,ratio)*...
Ang_d(load_dir(i,j,:),Norm(comp),ratio);
if m(i,j,1)<0
m(i,j,1) = 0;
end
if m(i,j,2)<0
m(i,j,2) = 0; end
if m(i,j,3)<0
m(i,j,3) = 0;
end
if m(i,j,4)<0
m(i,j,4) = 0;
end
if m(i,j,5)<0
m(i,j,5) = 0;
end
end
end
figure
title('Schmid Factor');
xlabel('Angle between load direction and basal [°]');
ylabel('Schmid factor,m');
beta = linspace(0,length(beta),length(beta));
hold on
plot(beta',m(:,:,1),'k','DisplayName','basal ');
plot(beta',m(:,:,2),'b');
plot(beta',m(:,:,3),'c');
plot(beta',m(:,:,4),'m');
plot(beta',m(:,:,5),'r');
hold off
xlswrite('SF_.xls', [m(1,:,1);m(31,:,1)]','basal');
xlswrite('SF_.xls', [m(1,:,2);m(31,:,2)]','pyr1');
xlswrite('SF_.xls', [m(1,:,3);m(31,:,3)]','pyr2');
xlswrite('SF_.xls', [m(1,:,4);m(31,:,4)]','ext');
xlswrite('SF_.xls', [m(1,:,5);m(31,:,5)]','comp');
CRSS = ones(length(alpha),length(beta),5)*1e3;
for i = 1: length(alpha)
for j = 1:length (beta)
if crss_basal/m(i,j,1)<1e3
CRSS(i,j,1) = crss_basal/m(i,j,1);
end
if crss_pyr/m(i,j,2)<1e3
CRSS(i,j,2) = crss_pyr/m(i,j,2);
end
if crss_pyr/m(i,j,3)<1e3
CRSS(i,j,3) = crss_pyr/m(i,j,3);
end
if crss_ext/m(i,j,4)<1e3
CRSS(i,j,4) = crss_ext/m(i,j,4);
end
if crss_comp/m(i,j,5)<1e3
CRSS(i,j,5) = crss_comp/m(i,j,5);
end
end
end
figure title('Critical Resolved Shear Stress');
xlabel('Angle between load direction and basal [°]');
ylabel('CRSS/m [MPa]');
hold on
plot(beta',CRSS(:,:,1),'k');
plot(beta',CRSS(:,:,2),'b');
plot(beta',CRSS(:,:,3),'c');
plot(beta',CRSS(:,:,4),'m');
plot(beta',CRSS(:,:,5),'r');
hold off
xlswrite('CRSS_.xls', [CRSS(1,:,1);CRSS(31,:,1)]','basal');
xlswrite('CRSS_.xls', [CRSS(1,:,2);CRSS(31,:,2)]','pyr1');
xlswrite('CRSS_.xls', [CRSS(1,:,3);CRSS(31,:,3)]','pyr2');
xlswrite('CRSS_.xls', [CRSS(1,:,4);CRSS(31,:,4)]','ext');
xlswrite('CRSS_.xls', [CRSS(1,:,5);CRSS(31,:,5)]','comp');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function angle_cos = Ang_d(A,B,C)
angle_cos = zeros(length(A(:,1)));
for i = 1:length(A(:,1))
angle_cos(i) = (A(i,1)*B(i,1)+A(i,2)*B(i,2)+(A(i,1)*B(i,2)+...
A(i,2)*B(i,1))/2+1/3*C^2*A(i,4)*B(i,4))/...
sqrt(A(i,1)^2+A(i,2)^2+A(i,1)*A(i,2)+1/3*A(i,4)^2*C^2)/...
sqrt(B(i,1)^2+B(i,2)^2+B(i,1)*B(i,2)+1/3*B(i,4)^2*C^2);
end
end
发布者：admin，转转请注明出处：http://www.yc00.com/xiaochengxu/1688383702a129871.html
评论列表（0条）