clear; clc;
n_s = 1e4;
Y2S = @ (x) x * 365 * 24 * 60 * 60; % Transform year into seconds
T_year = 40; % in years
% Generate onset and formulation time
pd_t_form = makedist ('Normal', 'mu', 9.3, 'sigma', 3.162);
t_form_year = random (pd_t_form); % Generate crack formulation time: in year
t_form = Y2S (t_form_year); % in seconds
% Tube parameters
pd_t_p = makedist ('Normal', 'mu', 8.3, 'sigma', 1/3); % Pressure difference: Mpa
P_diff = random (pd_t_p);
pd_t_d = makedist ('Normal', 'mu', 22.23, 'sigma',. 5/3); % Outer diameter: mm
d = random (pd_t_d);
pd_t_t = makedist ('Normal', 'mu', 1.27, 'sigma',. 125 * 1.27 / 3); % inner diameter mm
thickness = random (pd_t_t);
% S_u = 713; % Ultimate tensile strength, MPa
% S_y = 362; % Yeild strength, MPa
T = Y2S (T_year);
% Scotte model parameters: the velocity unit is mm / s, for non-cold workded
% Alloy 600
alpha = 2.8e-9;
K_th = 9;
beta = 1.16;
F = .93;
a_0 = .1; % Initial crack length: mm
a_rup = 62.5; % threshold value mm
t_current = t_form;
sigma_stress = P_diff * d / 2 / thickness; % Stress factor
% Initial values for recording the crack evolution process
a_k = a_0; % a_k: current crack length
t = t_current;
a = a_k;
time_step = 1;
index = 1; % Counter for saving the result
if t <T% If t_current does not exceed $ T $
while t
K = F * sigma_stress * sqrt (pi * a_k / 2); % Stress factor
a_kp = a_k + alpha * (K-K_th) ^ beta * time_step; % Calculate the crack length at the next time point
t_current = t_current + time_step; % Update the time point
a_k = a_kp; % Update current crack length
index = index + 1; % Update counter
time_step = time_step + 1;
if a_kp> a_rup% If the rupture happens
t (index) = t_current;
a (index) = a_kp;
break;
else if t_current> T
t (index) = T;
a (index) = a_kp;
break;
else
t (index) = t_current;
a (index) = a_kp;
end
end
end
else% If t_current exceeds $ T $
t = T; % Censored
end
end
plot (t, a)
xlabel ('Time'), ylabel ('Crack length')
