HCP施密特因子Schmid-factor-Mtalab计算.doc
《HCP施密特因子Schmid-factor-Mtalab计算.doc》由会员分享，可在线阅读，更多相关《HCP施密特因子Schmid-factor-Mtalab计算.doc（5页珍藏版）》请在冰点文库上搜索。
clear
clc
ratio=1.624;%forMg
crss_basal=2.1;%CRSSforbasalinMg
crss_pyr=40;%CRSSforinMg
crss_ext=11.5;%CRSSforextensiontwininMg
crss_comp=125;%CRSSforcompressiontwininMg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
basal=[0,0,0,1];%basalplane
pyr1=[-1,0,1,-1];%pyramidalIplane
pyr2=[-1,-1,2,-2];%pyramidalIIplane
ext=[1,0,-1,2];%extensiontwin
comp=[-1,0,1,-1];%compressiontwin
basal_dir=[1,1,-2,0];%basalBurgersvector
pyr1_dir=[1,1,-2,-3];%pyramidalBurgersvector
pyr2_dir=[1,1,-2,-3];%pyramidalBurgersvector
ext_dir=[1,0,-1,-1];%extensiontwinningdirection（yita1）
comp_dir=[1,0,-1,-2];%compressiontwinningdirection（yita1）
%fourindex
alpha=linspace（0,pi/6,31）;
beta=linspace（0,pi/2,91）;
load_dir=zeros（length（alpha）,length（beta）,4）;
fori=1:
length（alpha）
forj=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）;
fori=1:
length（alpha）
forj=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）;
ifm（i,j,1）<0
m（i,j,1）=0;
end
ifm（i,j,2）<0
m（i,j,2）=0;
end
ifm（i,j,3）<0
m（i,j,3）=0;
end
ifm（i,j,4）<0
m（i,j,4）=0;
end
ifm（i,j,5）<0
m（i,j,5）=0;
end
end
end
figure
title（'SchmidFactor'）;
xlabel（'Anglebetweenloaddirectionandbasal[°]'）;
ylabel（'Schmidfactor,m'）;
beta=linspace（0,length（beta）,length（beta））;
holdon
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'）;
holdoff
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;
fori=1:
length（alpha）
forj=1:
length（beta）
ifcrss_basal/m（i,j,1）<1e3
CRSS（i,j,1）=crss_basal/m（i,j,1）;
end
ifcrss_pyr/m（i,j,2）<1e3
CRSS（i,j,2）=crss_pyr/m（i,j,2）;
end
ifcrss_pyr/m（i,j,3）<1e3
CRSS（i,j,3）=crss_pyr/m（i,j,3）;
end
ifcrss_ext/m（i,j,4）<1e3
CRSS（i,j,4）=crss_ext/m（i,j,4）;
end
ifcrss_comp/m（i,j,5）<1e3
CRSS（i,j,5）=crss_comp/m（i,j,5）;
end
end
end
figure
title（'CriticalResolvedShearStress'）;
xlabel（'Anglebetweenloaddirectionandbasal[°]'）;
ylabel（'CRSS/m[MPa]'）;
holdon
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'）;
holdoff
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'）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionangle_cos=Ang_d（A,B,C）
angle_cos=zeros（length（A（:
1）））;
fori=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
升级会员 