calcul de concentration de carbone.m

碳性材料在高温情况下

clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcul de concentration de carbone 
% d'apres l'article "Predicting IGSCC/IGA susceptibility of Ni-Cr-Fe alloy
% by modelling of grain boundary chromium depletion" William E. Mayo
%
% Cette version necessite la Symbolic Math Toolbox
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Definition des constantes
%%%%%%%%%%%%%%%%%%%%%%%%%%%

C0 = 16 ;           % Concentration initiale en w/o
Beta = 0.1365 ;     % Coefficient de la formule de Zener
D = 4.91e-18 ;      % Diffusivité en micrometre/s^2
R = 12e-6 ;         % rayon du grain en m
tmax = 1000 ;       % secondes
pas_temps = 10 ;
nb_pas_espace = 100 ;

% % integrale de x*erf(x) + exp(-x^2)/(sqrt(pi)) entre 0 et R
% e = R*erf(R) + exp(-R^2)/(sqrt(pi)) - 0*erf(0) + exp(-0^2)/(sqrt(pi))


% Calcul des concentrations a l'interface
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = sym('x','real') ;

nb_pas_temps = tmax / pas_temps +1 ;

% initialisation de C_interface
C_interface = zeros(nb_pas_temps,1);
C_interface(1) = C0 ;

% initialisation de r_interface
r_interface = zeros(nb_pas_temps,1);
r_interface(1) = 0 ;

% initialisation de C
C(1).temps = 0 ;
C(1).C = C0 * ones(100,1) ;


t = 0 ;

for k = 2:nb_pas_temps
    t = t + pas_temps ;
    
    % Calcul de C_interface
    e(k) = double(int(1-erf(x/(2*sqrt(D*t))),x,0,inf)); 
    C_interface(k) = C0 - 0.91*Beta*sqrt(D*t)/e(k) ;
    
    % Calcul de r_interface
    r_interface(k) = Beta*sqrt(D*t)*10^6 ; 
    
    % calcul de C(x)
%     pas_espace = (R-r_interface(k))/(nb_pas_espace-1) ;
%     
%     C(k).temps = t ;
%     
%     j=0;
%     for i=r_interface(k):pas_espace:R
%         j = j+1;
%         C(k).C(j) = C_interface(k) + (C0-C_interface(k))*erf((i-r_interface(k))/(2*sqrt(D*t))) ;
%     end
%     
end
    
% Affichage des courbes
%%%%%%%%%%%%%%%%%%%%%%%




% for k = 2:nb_pas_temps
%     pas_espace = (R-r_interface(k))/(nb_pas_espace-1) ;
%     plot(r_interface(k):pas_espace:R,C(k).C)
% end