%use sigmaY which is C*QVol/omega/N, this also effects tc
R=8.3 ;%gas constant in J/mol/K
T_room=270;
Tmelt=1370; %melting of quartz in kelvin, or use Tw=1000C as in Rice 2006?
N=6*10^23 ;% avogadro # in 1/mol %quartz
G=31*10^9; %shear modulus in Pa
Cs=3750; %shear wavespeed in m/s
a_therm0=1*10^(-7);%1.2*10^(-6) m2/s is thermal diffusivity for quartz in m2/s from V&S
%from hanley fig 4. from 28*10^-7 to 6*10^-7
%Gibert and Mainprice show that if
%you assume there is porosity, their fig 8, this reduces diffusivity by 1 OM.
%from vosteen &Schellschmidt fig 5, T=0 alpha=1.2* 10^-6, and this is
%starting point from which we add porosity to reduce alpha by 50%
Q_vol=240000; % Rice2001 says 180KJ, Nakatanis says 90, 200, and 500 at different points in his paper. 200-250 is trypical for quartz
Q_GB=270000; %weakening occurs if QGB>QVol*B
%B=0.82*(4+273/T0)/5; to comply with about 58kJ activation energy as measured in Nakatani and scholz 2004.<1
rho=2650;%
c0=730; %specific heat capacity in J/kg/K
omega_vol=5.0*10^(-29); % activation m^3. for quartzite from Rice 2001.
omega_GB=7.8*10^(-29); % activation m^3. for quartzite from Rice 2001.
tauY=Q_GB/omega_GB/N; %2011 Putlet
%
% for quartzite molar mass =0.0064 kg/mol. divide by denisty 2650kg/m3,
% 2.415*10^-6 m3/mol. so omega*N=2.415*10^-6 m3/mol; divide by N=6*10^23--> omega 0.4*10^-29 ???
sigma=400*10^6; %applied normal stress
sigma0=2300*10^6 % ref normla stress - to be used to calculate
% Armax as function of stress
B0=0.89;
Ar0=0.95;
Ar_max=Ar0*(0.55+0.45*tanh(4*(sigma-sigma0)/sigma0));
%Ar_max is the contcat area that marks transition to ductility.
tc0=2;
imax=22;
Vt=0.5*Cs; % this must be large enough,
%larger than most V slip, or else the physics breaks.
d0=10*10^(-6);% initial contcat size in flash heating theory. proctor uses d=5micron
u=-1;
q=5;
V0=1*10^(-13);
D_therm=q*(sigma/10^6)^u; % Ditoro nature 2010.
b_room=R*T_room/Q_vol/B0;
Vn0=b_room*d0*exp((1-B0)*Q_vol/R/T_room)/tc0; %use measured tc0,e.g.=1s
%or 1000s at room T to calcualte Vn0, the max convergence V.
for i=1:imax
n=i/2+1;
V(i)=V0*10^(n);
t(i)=d0/V(i);
dc(i)=d0;
end
% From Ar(0)=n*pi*(d0/2)^2, and Ar(0)*sigmaY=A*sigma. this gives at
% refernce normla stress sigma0
sigmaY0=B0*Q_vol/omega_vol/N;
k=1*sigma0;
for j=1:15
T0(j)=300+(j-1)*70 % different ambient T
B(j)=B0*exp(-0.0006*(T0(j)-T_room));
% fit from Evans JGR 84, indentation hardness of quartz
% to allow sigmaY to vary with T.
if B(j)>1
B(j)=1;
end
mu0(j)=Q_GB*omega_vol/Q_vol/omega_GB/B(j);
a0(j)=R*T0(j)/Q_GB;
b0(j)=R*T0(j)/Q_vol/B(j);
tau_c=tauY*(1+a0(j)*log(V/Vt));
sigmaY(j)=B(j)*Q_vol/omega_vol/N; %this is the stress at
%time of contact onset .
Tsurf=T0(j)+sigma*tau_c.*sqrt(V.*D_therm/pi./a_therm0)/rho/c0/ ...
sigmaY(j); %surface T needs tau not tauc, thats why it
%multiplied vy sigma and divided by sigmaY
T=Tsurf+tau_c.*sqrt(V*d0/pi./a_therm0)/rho/c0; % contact T , Proctor JGR 2014, eq 2.
for i=1:imax
if T(i)> 0.99*Tmelt
T(i)=Tmelt; %saturate T for T>Tmelt or T_decomposition.
end
end
for k=2:111 %iterate to converege T and tau for obtaining st-st
for i=1:imax
a(i)=R*T(i)/Q_GB;
b(i)=R*T(i)/Q_vol/B(j);
tc(i)=b(i)*(dc(i)/Vn0)*exp((1-B0)*Q_vol/R/T(i));
% the critical time at each T
%I add that at ductility sigma is saturated A which is sigmaY at this level
% i know this from runs that show me that this is
% the value of sigmac for v low strain rate (10^-12)
% this is sigma_c calculated using max Ar
sigma_c(i)=max(sigmaY(j)*(1-b(i)*log(1+dc(i)/V(i)/tc(i))),sigma/Ar_max);
dc(i)=d0*sqrt(sigmaY(1)/sigma_c(i));
tau_c(i)=tauY*(1+a(i)*log(V(i)/Vt));
%recalculate T now with corrected thermal parameters.
a_therm(i)=1*10^(-4)./T(i)-0.5*10^(-7);%fig 4 from hanley
c(i)=c0*(1.7-200/T(i)); % from fitting fig 4, Vosteen & Schellschmidt
ff=sqrt(V(i)/pi./a_therm(i))/rho/c(i);
Tsurf(i)=T0(j)+sigma*tau_c(i)*ff*sqrt(D_therm)/sigma_c(i);
T(i)=Tsurf(i)+tau_c(i)*sqrt(dc(i))*ff;
if T(i)> 0.95*Tmelt
T(i)=Tmelt;
tau_c(i)=0.95*tau_c(i);
a_therm(i)=a_therm(i-1);
c(i)=c(i-1);
elseif T(i)<270
T(i)=500;
tau_c(i)=tauY*(1+a(i)*log(V(i)/Vt));
end
sigmaCon(i,j)=sigma_c(i);
end
end
for i=1:imax
if T(i)> 0.95*Tmelt
T(i)=Tmelt;
a_therm(i)=a_therm(i)*0.7;%fig 4 from hanley
c(i)=c(i-1); % from Vosteen & Schellschmidt
ff=sqrt(V(i)/pi./a_therm(i))/rho/c(i);
sigma_c(i)=sigma_c(i-1);
dc(i)=d0*sqrt(sigmaY(1)/sigma_c(i));
tau_c(i)=(Tmelt-T0(j))/ff/(sqrt(D_therm)*sigma/sigma_c(i)+sqrt(dc(i)));
end
end
%end
TCon(:,j)=T;
dcCon(:,j)=dc;
tauCon(:,j)=tau_c;
sigmaCon(:,j)=sigma_c;
friction(:,j)=tau_c./sigma_c;
tau(:,j)=sigma*friction(:,j);
power(:,j)=tau(:,j).*V(:);
tcCon(:,j)=tc;
aCon(:,j)=a;
bCon(:,j)=b;
dcCon(:,j)=dc;
a_thermCon(:,j)=a_therm;
mu0Con=mu0(j);
mu_approxCon(:,j)=mu0(j)+a.*log(V/Vt)+b.*log(1+dc./tc./V);
mu_fullCon(:,j)=mu0(j)*(1+a.*log(V/Vt))./(1-b.*log(1+dc./tc./V));
end
figure(1)
% semilogx(V,sigma./sigmaCon(:,1),'r+-',V,sigma./sigmaCon(:,3),'+-g',V,sigma./sigmaCon(:,4),'k+-',V,sigma./sigmaCon(:,5),'c+-',V,sigma./sigmaCon(:,6),'b+-',V,sigma./sigmaCon(:,7),'+-y')
% legend(int2str(T0(1)-270),int2str(T0(3)-270),int2str(T0(4)-270),int2str(T0(5)-270),int2str(T0(6)-270),int2str(T0(7)-270))
semilogx(V,sigma./sigmaCon(:,1),'r+-',V,sigma./sigmaCon(:,3),'c+-',V,sigma./sigmaCon(:,4),'k+-',V,Ar_max*ones(size(V)),'-k',V,sigma./sigmaCon(:,6),'k+-',V,sigma./sigmaCon(:,8),'g+-')
ylabel('Ar')
xlabel('V')
hold on
%
%
%
% figure(10)
% plot(T0-270,Ar_max,'g+-')
% legend(num2str(V(8)),num2str(V(10)),num2str(V(1)))
% ylabel('A_r/A')
% xlabel('T_0(C)')
% hold on
figure(3)
plot(T0-270,sigma./sigmaCon(8,:),'r+-',T0-270,sigma./sigmaCon(10,:),'k+-',T0-270,sigma./sigmaCon(12,:),'c+-',T0-270,Ar_max,'g+-')
legend(num2str(V(8)),num2str(V(10)),num2str(V(12)))
ylabel('A_r/A')
xlabel('T_0(C)')
hold on
figure(4)
%plot(T0(4:end)-270,friction(10,4:end))
plot(T0-270,friction(8,:),'r+-',T0-270,friction(10,:),'k+-',T0-270,friction(12,:),'c+-')
%legend(num2str(V(10)))
legend(num2str(V(8)),num2str(V(10)),num2str(V(12)))
ylabel('friction')
xlabel('T_0(C)')
hold on
% figure(2)
% semilogy(T0-270,dhdtCon(1,:),'r+-',T0-270,dhdtCon(6,:),'k+-',T0-270,dhdtCon(10,:),'g+-',T0-270,dhdtCon(12,:),'y+-',T0-270,dhdtCon(14,:),'c+-')
% legend(num2str(V(1)),num2str(V(6)),num2str(V(10)),num2str(V(12)),num2str(V(14)))
% ylabel('dhdt')
% xlabel('T_0(C)')
% hold on
figure(6)
semilogy(T0-270,sigmaCon(1,:),'r+-',T0-270,sigmaCon(6,:),'k+-',T0-270,sigmaCon(12,:),'c+-')
legend(num2str(V(1)),num2str(V(6)),num2str(V(12)))
ylabel('sigmac')
xlabel('T_0(C)')
hold on
figure(5)
loglog(V,sigmaCon(:,1),'r+-',V,sigmaCon(:,5),'k+-',V,sigmaCon(:,8),'c+-')
legend(num2str(T0(1)),num2str(T0(5)),num2str(T0(8)))
ylabel('sigmac')
xlabel('V')
hold on
% figure(9)
% semilogx(V,tauCon(:,1),'r',V,tauCon(:,2),'k',V,tauCon(:,3),'c',V,tauCon(:,4),'b',V,tauCon(:,5),'y',V,tauCon(:,6),'g')
% hold on
% semilogx(V,sigmaCon(:,1),'*r',V,sigmaCon(:,2),'*k',V,sigmaCon(:,3),'*c',V,sigmaCon(:,4),'*b',V,sigmaCon(:,5),'*y',V,sigmaCon(:,6),'*g')
% legend(int2str(T0(1)-270),int2str(T0(2)-270),int2str(T0(3)-270),int2str(T0(4)-270),int2str(T0(5)-270),int2str(T0(6)-270))
% ylabel('shear and normal stress')
% xlabel('V')
% hold on
%
% % figure(1)
% % semilogx(V,sigmaCon(:,1),'r',V,sigmaCon(:,2),'k',V,sigmaCon(:,3),'c',V,sigmaCon(:,4),'b',V,sigmaCon(:,5),'y',V,sigmaCon(:,6),'g')
% % hold on
% % legend(int2str(T0(1)-270),int2str(T0(2)-270),int2str(T0(3)-270),int2str(T0(4)-270),int2str(T0(5)-270),int2str(T0(6)-270))
% % ylabel('normal stress')
% % xlabel('V')
% % hold on
%
%
figure(2)
semilogx(V,friction(:,1),'-r',V,friction(:,3),'-g',V,friction(:,4),'-k',V,friction(:,5),'-c',V,friction(:,6),'-b',V,friction(:,8),'-y')
hold on
% semilogx(V,mu_approxCon(:,1),'-r',V,mu_approxCon(:,3),'-g',V,mu_approxCon(:,4),'-k',V,mu_approxCon(:,5),'-c',V,mu_approxCon(:,6),'-b',V,mu_approxCon(:,7),'-y')
% hold on
% semilogx(V,mu_fullCon(:,1),'+r',V,mu_fullCon(:,3),'+g',V,mu_fullCon(:,4),'+k',V,mu_fullCon(:,5),'+c',V,mu_fullCon(:,6),'+b',V,mu_fullCon(:,7),'+y')
ylabel('friction')
xlabel('V')
legend(int2str(T0(1)-270),int2str(T0(3)-270),int2str(T0(4)-270),int2str(T0(5)-270),int2str(T0(6)-270),int2str(T0(8)-270))
hold on
%
% figure(7)
%
% semilogx(V,TCon(:,1)-273,'*-r',V,TCon(:,3)-273,'*-k',V,TCon(:,5)-273,'*-c',V,TCon(:,7)-273,'o-b',V,TCon(:,9)-273,'o-y')
% ylabel('T (C)')
% xlabel('V (m/s)')
% legend(int2str(T0(1)-270),int2str(T0(3)-270),int2str(T0(5)-270),int2str(T0(7)-270),int2str(T0(9)-270))
% hold on
%
% figure(17)
%
% loglog(V,tcCon(:,1),'*-r',V,tcCon(:,3),'*-k',V,tcCon(:,5),'*-c',V,tcCon(:,7),'o-b',V,tcCon(:,9),'o-y')
% ylabel('tc')
% xlabel('V (m/s)')
% legend(int2str(T0(1)-270),int2str(T0(3)-270),int2str(T0(5)-270),int2str(T0(7)-270),int2str(T0(9)-270))
% hold on
%
% % figure(5)
% % plot(1:k,Ttime)
% % ylabel('T')
% % xlabel('time')
% % hold on
%
% % figure(5)
% % plot(1:k,tautime,'r',1:k,sigmatime,'b')
% % ylabel('tau')
% % xlabel('time')
% % hold on
R=8.3 ;%gas constant in J/mol/K
T_room=270;
Tmelt=1370; %melting of quartz in kelvin, or use Tw=1000C as in Rice 2006?
N=6*10^23 ;% avogadro # in 1/mol %quartz
G=31*10^9; %shear modulus in Pa
Cs=3750; %shear wavespeed in m/s
a_therm0=1*10^(-7);%1.2*10^(-6) m2/s is thermal diffusivity for quartz in m2/s from V&S
%from hanley fig 4. from 28*10^-7 to 6*10^-7
%Gibert and Mainprice show that if
%you assume there is porosity, their fig 8, this reduces diffusivity by 1 OM.
%from vosteen &Schellschmidt fig 5, T=0 alpha=1.2* 10^-6, and this is
%starting point from which we add porosity to reduce alpha by 50%
Q_vol=240000; % Rice2001 says 180KJ, Nakatanis says 90, 200, and 500 at different points in his paper. 200-250 is trypical for quartz
Q_GB=270000; %weakening occurs if QGB>QVol*B
%B=0.82*(4+273/T0)/5; to comply with about 58kJ activation energy as measured in Nakatani and scholz 2004.<1
rho=2650;%
c0=730; %specific heat capacity in J/kg/K
omega_vol=5.0*10^(-29); % activation m^3. for quartzite from Rice 2001.
omega_GB=7.8*10^(-29); % activation m^3. for quartzite from Rice 2001.
tauY=Q_GB/omega_GB/N; %2011 Putlet
%
% for quartzite molar mass =0.0064 kg/mol. divide by denisty 2650kg/m3,
% 2.415*10^-6 m3/mol. so omega*N=2.415*10^-6 m3/mol; divide by N=6*10^23--> omega 0.4*10^-29 ???
sigma=400*10^6; %applied normal stress
sigma0=2300*10^6 % ref normla stress - to be used to calculate
% Armax as function of stress
B0=0.89;
Ar0=0.95;
Ar_max=Ar0*(0.55+0.45*tanh(4*(sigma-sigma0)/sigma0));
%Ar_max is the contcat area that marks transition to ductility.
tc0=2;
imax=22;
Vt=0.5*Cs; % this must be large enough,
%larger than most V slip, or else the physics breaks.
d0=10*10^(-6);% initial contcat size in flash heating theory. proctor uses d=5micron
u=-1;
q=5;
V0=1*10^(-13);
D_therm=q*(sigma/10^6)^u; % Ditoro nature 2010.
b_room=R*T_room/Q_vol/B0;
Vn0=b_room*d0*exp((1-B0)*Q_vol/R/T_room)/tc0; %use measured tc0,e.g.=1s
%or 1000s at room T to calcualte Vn0, the max convergence V.
for i=1:imax
n=i/2+1;
V(i)=V0*10^(n);
t(i)=d0/V(i);
dc(i)=d0;
end
% From Ar(0)=n*pi*(d0/2)^2, and Ar(0)*sigmaY=A*sigma. this gives at
% refernce normla stress sigma0
sigmaY0=B0*Q_vol/omega_vol/N;
k=1*sigma0;
for j=1:15
T0(j)=300+(j-1)*70 % different ambient T
B(j)=B0*exp(-0.0006*(T0(j)-T_room));
% fit from Evans JGR 84, indentation hardness of quartz
% to allow sigmaY to vary with T.
if B(j)>1
B(j)=1;
end
mu0(j)=Q_GB*omega_vol/Q_vol/omega_GB/B(j);
a0(j)=R*T0(j)/Q_GB;
b0(j)=R*T0(j)/Q_vol/B(j);
tau_c=tauY*(1+a0(j)*log(V/Vt));
sigmaY(j)=B(j)*Q_vol/omega_vol/N; %this is the stress at
%time of contact onset .
Tsurf=T0(j)+sigma*tau_c.*sqrt(V.*D_therm/pi./a_therm0)/rho/c0/ ...
sigmaY(j); %surface T needs tau not tauc, thats why it
%multiplied vy sigma and divided by sigmaY
T=Tsurf+tau_c.*sqrt(V*d0/pi./a_therm0)/rho/c0; % contact T , Proctor JGR 2014, eq 2.
for i=1:imax
if T(i)> 0.99*Tmelt
T(i)=Tmelt; %saturate T for T>Tmelt or T_decomposition.
end
end
for k=2:111 %iterate to converege T and tau for obtaining st-st
for i=1:imax
a(i)=R*T(i)/Q_GB;
b(i)=R*T(i)/Q_vol/B(j);
tc(i)=b(i)*(dc(i)/Vn0)*exp((1-B0)*Q_vol/R/T(i));
% the critical time at each T
%I add that at ductility sigma is saturated A which is sigmaY at this level
% i know this from runs that show me that this is
% the value of sigmac for v low strain rate (10^-12)
% this is sigma_c calculated using max Ar
sigma_c(i)=max(sigmaY(j)*(1-b(i)*log(1+dc(i)/V(i)/tc(i))),sigma/Ar_max);
dc(i)=d0*sqrt(sigmaY(1)/sigma_c(i));
tau_c(i)=tauY*(1+a(i)*log(V(i)/Vt));
%recalculate T now with corrected thermal parameters.
a_therm(i)=1*10^(-4)./T(i)-0.5*10^(-7);%fig 4 from hanley
c(i)=c0*(1.7-200/T(i)); % from fitting fig 4, Vosteen & Schellschmidt
ff=sqrt(V(i)/pi./a_therm(i))/rho/c(i);
Tsurf(i)=T0(j)+sigma*tau_c(i)*ff*sqrt(D_therm)/sigma_c(i);
T(i)=Tsurf(i)+tau_c(i)*sqrt(dc(i))*ff;
if T(i)> 0.95*Tmelt
T(i)=Tmelt;
tau_c(i)=0.95*tau_c(i);
a_therm(i)=a_therm(i-1);
c(i)=c(i-1);
elseif T(i)<270
T(i)=500;
tau_c(i)=tauY*(1+a(i)*log(V(i)/Vt));
end
sigmaCon(i,j)=sigma_c(i);
end
end
for i=1:imax
if T(i)> 0.95*Tmelt
T(i)=Tmelt;
a_therm(i)=a_therm(i)*0.7;%fig 4 from hanley
c(i)=c(i-1); % from Vosteen & Schellschmidt
ff=sqrt(V(i)/pi./a_therm(i))/rho/c(i);
sigma_c(i)=sigma_c(i-1);
dc(i)=d0*sqrt(sigmaY(1)/sigma_c(i));
tau_c(i)=(Tmelt-T0(j))/ff/(sqrt(D_therm)*sigma/sigma_c(i)+sqrt(dc(i)));
end
end
%end
TCon(:,j)=T;
dcCon(:,j)=dc;
tauCon(:,j)=tau_c;
sigmaCon(:,j)=sigma_c;
friction(:,j)=tau_c./sigma_c;
tau(:,j)=sigma*friction(:,j);
power(:,j)=tau(:,j).*V(:);
tcCon(:,j)=tc;
aCon(:,j)=a;
bCon(:,j)=b;
dcCon(:,j)=dc;
a_thermCon(:,j)=a_therm;
mu0Con=mu0(j);
mu_approxCon(:,j)=mu0(j)+a.*log(V/Vt)+b.*log(1+dc./tc./V);
mu_fullCon(:,j)=mu0(j)*(1+a.*log(V/Vt))./(1-b.*log(1+dc./tc./V));
end
figure(1)
% semilogx(V,sigma./sigmaCon(:,1),'r+-',V,sigma./sigmaCon(:,3),'+-g',V,sigma./sigmaCon(:,4),'k+-',V,sigma./sigmaCon(:,5),'c+-',V,sigma./sigmaCon(:,6),'b+-',V,sigma./sigmaCon(:,7),'+-y')
% legend(int2str(T0(1)-270),int2str(T0(3)-270),int2str(T0(4)-270),int2str(T0(5)-270),int2str(T0(6)-270),int2str(T0(7)-270))
semilogx(V,sigma./sigmaCon(:,1),'r+-',V,sigma./sigmaCon(:,3),'c+-',V,sigma./sigmaCon(:,4),'k+-',V,Ar_max*ones(size(V)),'-k',V,sigma./sigmaCon(:,6),'k+-',V,sigma./sigmaCon(:,8),'g+-')
ylabel('Ar')
xlabel('V')
hold on
%
%
%
% figure(10)
% plot(T0-270,Ar_max,'g+-')
% legend(num2str(V(8)),num2str(V(10)),num2str(V(1)))
% ylabel('A_r/A')
% xlabel('T_0(C)')
% hold on
figure(3)
plot(T0-270,sigma./sigmaCon(8,:),'r+-',T0-270,sigma./sigmaCon(10,:),'k+-',T0-270,sigma./sigmaCon(12,:),'c+-',T0-270,Ar_max,'g+-')
legend(num2str(V(8)),num2str(V(10)),num2str(V(12)))
ylabel('A_r/A')
xlabel('T_0(C)')
hold on
figure(4)
%plot(T0(4:end)-270,friction(10,4:end))
plot(T0-270,friction(8,:),'r+-',T0-270,friction(10,:),'k+-',T0-270,friction(12,:),'c+-')
%legend(num2str(V(10)))
legend(num2str(V(8)),num2str(V(10)),num2str(V(12)))
ylabel('friction')
xlabel('T_0(C)')
hold on
% figure(2)
% semilogy(T0-270,dhdtCon(1,:),'r+-',T0-270,dhdtCon(6,:),'k+-',T0-270,dhdtCon(10,:),'g+-',T0-270,dhdtCon(12,:),'y+-',T0-270,dhdtCon(14,:),'c+-')
% legend(num2str(V(1)),num2str(V(6)),num2str(V(10)),num2str(V(12)),num2str(V(14)))
% ylabel('dhdt')
% xlabel('T_0(C)')
% hold on
figure(6)
semilogy(T0-270,sigmaCon(1,:),'r+-',T0-270,sigmaCon(6,:),'k+-',T0-270,sigmaCon(12,:),'c+-')
legend(num2str(V(1)),num2str(V(6)),num2str(V(12)))
ylabel('sigmac')
xlabel('T_0(C)')
hold on
figure(5)
loglog(V,sigmaCon(:,1),'r+-',V,sigmaCon(:,5),'k+-',V,sigmaCon(:,8),'c+-')
legend(num2str(T0(1)),num2str(T0(5)),num2str(T0(8)))
ylabel('sigmac')
xlabel('V')
hold on
% figure(9)
% semilogx(V,tauCon(:,1),'r',V,tauCon(:,2),'k',V,tauCon(:,3),'c',V,tauCon(:,4),'b',V,tauCon(:,5),'y',V,tauCon(:,6),'g')
% hold on
% semilogx(V,sigmaCon(:,1),'*r',V,sigmaCon(:,2),'*k',V,sigmaCon(:,3),'*c',V,sigmaCon(:,4),'*b',V,sigmaCon(:,5),'*y',V,sigmaCon(:,6),'*g')
% legend(int2str(T0(1)-270),int2str(T0(2)-270),int2str(T0(3)-270),int2str(T0(4)-270),int2str(T0(5)-270),int2str(T0(6)-270))
% ylabel('shear and normal stress')
% xlabel('V')
% hold on
%
% % figure(1)
% % semilogx(V,sigmaCon(:,1),'r',V,sigmaCon(:,2),'k',V,sigmaCon(:,3),'c',V,sigmaCon(:,4),'b',V,sigmaCon(:,5),'y',V,sigmaCon(:,6),'g')
% % hold on
% % legend(int2str(T0(1)-270),int2str(T0(2)-270),int2str(T0(3)-270),int2str(T0(4)-270),int2str(T0(5)-270),int2str(T0(6)-270))
% % ylabel('normal stress')
% % xlabel('V')
% % hold on
%
%
figure(2)
semilogx(V,friction(:,1),'-r',V,friction(:,3),'-g',V,friction(:,4),'-k',V,friction(:,5),'-c',V,friction(:,6),'-b',V,friction(:,8),'-y')
hold on
% semilogx(V,mu_approxCon(:,1),'-r',V,mu_approxCon(:,3),'-g',V,mu_approxCon(:,4),'-k',V,mu_approxCon(:,5),'-c',V,mu_approxCon(:,6),'-b',V,mu_approxCon(:,7),'-y')
% hold on
% semilogx(V,mu_fullCon(:,1),'+r',V,mu_fullCon(:,3),'+g',V,mu_fullCon(:,4),'+k',V,mu_fullCon(:,5),'+c',V,mu_fullCon(:,6),'+b',V,mu_fullCon(:,7),'+y')
ylabel('friction')
xlabel('V')
legend(int2str(T0(1)-270),int2str(T0(3)-270),int2str(T0(4)-270),int2str(T0(5)-270),int2str(T0(6)-270),int2str(T0(8)-270))
hold on
%
% figure(7)
%
% semilogx(V,TCon(:,1)-273,'*-r',V,TCon(:,3)-273,'*-k',V,TCon(:,5)-273,'*-c',V,TCon(:,7)-273,'o-b',V,TCon(:,9)-273,'o-y')
% ylabel('T (C)')
% xlabel('V (m/s)')
% legend(int2str(T0(1)-270),int2str(T0(3)-270),int2str(T0(5)-270),int2str(T0(7)-270),int2str(T0(9)-270))
% hold on
%
% figure(17)
%
% loglog(V,tcCon(:,1),'*-r',V,tcCon(:,3),'*-k',V,tcCon(:,5),'*-c',V,tcCon(:,7),'o-b',V,tcCon(:,9),'o-y')
% ylabel('tc')
% xlabel('V (m/s)')
% legend(int2str(T0(1)-270),int2str(T0(3)-270),int2str(T0(5)-270),int2str(T0(7)-270),int2str(T0(9)-270))
% hold on
%
% % figure(5)
% % plot(1:k,Ttime)
% % ylabel('T')
% % xlabel('time')
% % hold on
%
% % figure(5)
% % plot(1:k,tautime,'r',1:k,sigmatime,'b')
% % ylabel('tau')
% % xlabel('time')
% % hold on