%% FATIGUE CRACK GROWTH VIDEO CREATION
% This script creates an educational video demonstrating fatigue crack
growth
% analysis using Paris' law with animations and progressive visualization.
%
% The video covers:
% 1. Introduction to fracture mechanics concepts
% 2. Paris Law visualization
% 3. Real-time crack growth simulation
% 4. Parametric study visualization
% 5. Design implications and conclusions
clear; clc; close all;
%% VIDEO SETUP
videoFileName = 'Fatigue_Crack_Growth_Analysis.mp4';
fps = 20; % Frames per second for smooth animation
videoDuration = 45; % Total video duration in seconds
totalFrames = videoDuration * fps;
% Create video writer object
videoWriter = VideoWriter(videoFileName, 'MPEG-4');
[Link] = fps;
open(videoWriter);
fprintf('Creating fatigue analysis video...\n');
fprintf('Video will be saved as: %s\n', videoFileName);
fprintf('Duration: %.1f seconds at %d fps\n\n', videoDuration, fps);
%% MATERIAL AND GEOMETRY PARAMETERS
% Same as original script
C = 6.9e-12;
m = 3.0;
K_Ic = 24;
W = 0.1;
t = 0.005;
sigma_max = 150;
R_ratio = 0.1;
sigma_min = R_ratio * sigma_max;
Delta_sigma = sigma_max - sigma_min;
a0 = 0.001;
%% PRELIMINARY CALCULATIONS
% Calculate critical crack length
a_test = linspace(a0, W/2, 1000);
K_max_test = zeros(size(a_test));
for i = 1:length(a_test)
[K_max_test(i), ~] = calculateStressIntensity(a_test(i), sigma_max,
sigma_min, W);
end
idx_crit = find(K_max_test >= K_Ic, 1, 'first');
if isempty(idx_crit)
ac = W/2;
else
ac = a_test(idx_crit);
end
%% PART 1: INTRODUCTION FRAME
fig = figure('Position', [100, 100, 1200, 700], 'Color', 'white');
Page 1 of 11
% Title and introduction
subplot('Position', [0.1, 0.7, 0.8, 0.25]);
axis off;
text(0.5, 0.7, 'FATIGUE CRACK GROWTH ANALYSIS', ...
'HorizontalAlignment', 'center', 'FontSize', 28, 'FontWeight', 'bold');
text(0.5, 0.4, 'Using Paris Law and Fracture Mechanics', ...
'HorizontalAlignment', 'center', 'FontSize', 18, 'Color', [0.2, 0.2,
0.6]);
text(0.5, 0.2, sprintf('Material: Aluminum 7075-T6 | Stress Range: %.0f
MPa', Delta_sigma), ...
'HorizontalAlignment', 'center', 'FontSize', 14);
% Paris Law equation
subplot('Position', [0.1, 0.4, 0.8, 0.2]);
axis off;
text(0.5, 0.8, 'Paris Law Governing Equation:', ...
'HorizontalAlignment', 'center', 'FontSize', 16);
text(0.5, 0.5, '$$\frac{da}{dN} = C (\Delta K)^m$$', ...
'HorizontalAlignment', 'center', 'FontSize', 24, 'Interpreter',
'latex');
text(0.5, 0.2, sprintf('where: C = %.1e, m = %.1f', C, m), ...
'HorizontalAlignment', 'center', 'FontSize', 14);
% Plate schematic
subplot('Position', [0.15, 0.1, 0.7, 0.25]);
hold on; box on;
rectangle('Position', [0, -0.5, W*10, 1], 'FaceColor', [0.8, 0.8, 1],
'EdgeColor', 'b', 'LineWidth', 2);
plot([W*5-0.5, W*5+0.5], [0, 0], 'r-', 'LineWidth', 3);
text(W*5, 0.6, '2a', 'HorizontalAlignment', 'center', 'FontSize', 12);
text(W*5, -0.8, 'σ_max', 'HorizontalAlignment', 'center', 'FontSize', 12);
% Draw arrow using custom function
drawArrow([W*5, -0.7], [W*5, -0.3], 'Color', [0, 0, 0.7], 'LineWidth', 2);
title('Center-Cracked Plate Configuration', 'FontSize', 12);
axis equal; axis([0, W*10, -1, 1]);
% Write frame
frame = getframe(fig);
writeVideo(videoWriter, frame);
pause(1); % Pause for emphasis
%% PART 2: PARIS LAW ANIMATION
% Animate the Paris law curve
fig = figure('Position', [100, 100, 1200, 700], 'Color', 'white');
for frameNum = 1:50
clf;
% Calculate Delta K range
DeltaK_range = logspace(0, log10(40), 100);
dadN_range = C * DeltaK_range.^m;
% Main plot
subplot('Position', [0.1, 0.1, 0.55, 0.8]);
loglog(DeltaK_range, dadN_range*1e6, 'b-', 'LineWidth', 3);
hold on; grid on;
Page 2 of 11
% Progressively plot the curve
idx = round(frameNum/50 * length(DeltaK_range));
if idx > 0
loglog(DeltaK_range(1:idx), dadN_range(1:idx)*1e6, 'r-',
'LineWidth', 4);
end
xlabel('Stress Intensity Range, \DeltaK (MPa√m)', 'FontSize', 12);
ylabel('Crack Growth Rate, da/dN (μm/cycle)', 'FontSize', 12);
title('Paris Law: da/dN vs \DeltaK', 'FontSize', 14, 'FontWeight',
'bold');
legend('Paris Law Curve', 'Progress', 'Location', 'northwest');
% Add equation
annotation('textbox', [0.7, 0.7, 0.25, 0.15], ...
'String', {'Paris Law:', '$$ \frac{da}{dN} = C (\Delta K)^m $$',
...
sprintf('C = %.1e', C), sprintf('m = %.1f', m)}, ...
'FontSize', 14, 'Interpreter', 'latex', ...
'BackgroundColor', [0.95, 0.95, 1], 'EdgeColor', 'blue');
% Progress indicator
progress = frameNum/50;
annotation('textbox', [0.7, 0.5, 0.25, 0.1], ...
'String', sprintf('Plotting Paris Law Curve\nProgress: %.0f%%',
progress*100), ...
'FontSize', 12, 'BackgroundColor', [0.9, 0.9, 0.9]);
frame = getframe(fig);
writeVideo(videoWriter, frame);
end
%% PART 3: CRACK GROWTH SIMULATION ANIMATION
fig = figure('Position', [100, 100, 1200, 700], 'Color', 'white');
% Pre-calculate the crack growth for animation
dN_anim = 100;
a_current = a0;
N_current = 0;
a_history = a0;
N_history = 0;
DeltaK_history = [];
while a_current < ac
[~, Delta_K] = calculateStressIntensity(a_current, sigma_max,
sigma_min, W);
if Delta_K <= 0
break;
end
dadN = C * (Delta_K)^m;
da = dadN * dN_anim;
a_current = a_current + da;
N_current = N_current + dN_anim;
a_history = [a_history; a_current];
N_history = [N_history; N_current];
DeltaK_history = [DeltaK_history; Delta_K];
Page 3 of 11
if a_current >= ac || a_current >= W/2
break;
end
end
% Animate the crack growth
numFramesAnimation = min(100, length(a_history));
frameStep = max(1, floor(length(a_history)/numFramesAnimation));
for frameNum = 1:frameStep:length(a_history)
clf;
% Left: Crack growth curve
subplot('Position', [0.08, 0.1, 0.4, 0.8]);
plot(N_history(1:frameNum)/1000, a_history(1:frameNum)*1000, 'b-',
'LineWidth', 3);
hold on; grid on;
% Mark critical point
plot(N_history(frameNum)/1000, a_history(frameNum)*1000, 'ro', ...
'MarkerSize', 10, 'MarkerFaceColor', 'r');
% Add critical crack length line
yline(ac*1000, 'r--', 'LineWidth', 2, 'Label', 'Critical Crack
Length');
xlabel('Number of Cycles, N (×10³)', 'FontSize', 12);
ylabel('Crack Length, a (mm)', 'FontSize', 12);
title('Crack Growth vs Cycles', 'FontSize', 14, 'FontWeight', 'bold');
legend('Crack Growth', 'Current State', 'Location', 'northwest');
xlim([0, max(N_history)/1000]);
ylim([0, ac*1000*1.1]);
% Right: Paris law with current point
subplot('Position', [0.58, 0.1, 0.4, 0.8]);
DeltaK_plot = logspace(0, log10(max(DeltaK_history)*1.2), 100);
dadN_plot = C * DeltaK_plot.^m;
loglog(DeltaK_plot, dadN_plot*1e6, 'b-', 'LineWidth', 2);
hold on; grid on;
% Plot history
if frameNum > 1
loglog(DeltaK_history(1:frameNum), ...
diff(a_history(1:frameNum+1))/dN_anim*1e6, 'ro', ...
'MarkerSize', 6, 'MarkerFaceColor', 'r');
end
% Current point
if frameNum > 1
current_dadN = diff(a_history(frameNum-1:frameNum))/dN_anim*1e6;
loglog(DeltaK_history(frameNum-1), current_dadN, 'go', ...
'MarkerSize', 12, 'MarkerFaceColor', 'g');
end
xlabel('\DeltaK (MPa√m)', 'FontSize', 12);
ylabel('da/dN (μm/cycle)', 'FontSize', 12);
title('Paris Law Verification', 'FontSize', 14, 'FontWeight', 'bold');
legend('Theoretical', 'History', 'Current', 'Location', 'northwest');
Page 4 of 11
% Information box
annotation('textbox', [0.1, 0.85, 0.3, 0.1], ...
'String', {sprintf('Current Crack Length: %.2f mm',
a_history(frameNum)*1000), ...
sprintf('Current Cycles: %.0f', N_history(frameNum)), ...
sprintf('Current \\DeltaK: %.2f MPa√m',
DeltaK_history(min(frameNum, length(DeltaK_history))))}, ...
'FontSize', 12, 'BackgroundColor', [0.95, 0.95, 1], 'EdgeColor',
'blue');
frame = getframe(fig);
writeVideo(videoWriter, frame);
end
%% PART 4: PARAMETRIC STUDY ANIMATION
fig = figure('Position', [100, 100, 1200, 700], 'Color', 'white');
% Perform parametric study
Delta_sigma_range = linspace(80, 180, 20);
a0_range = linspace(0.0005, 0.005, 20);
Nf_matrix = zeros(length(a0_range), length(Delta_sigma_range));
fprintf('Computing parametric study for animation...\n');
for i = 1:length(a0_range)
for j = 1:length(Delta_sigma_range)
sig_max = Delta_sigma_range(j);
sig_min = R_ratio * sig_max;
% Find critical length
a_test_local = linspace(a0_range(i), W/2, 200);
K_max_test_local = zeros(size(a_test_local));
for k = 1:length(a_test_local)
[K_max_test_local(k), ~] =
calculateStressIntensity(a_test_local(k), sig_max, sig_min, W);
end
idx_crit_local = find(K_max_test_local >= K_Ic, 1, 'first');
if isempty(idx_crit_local)
ac_local = W/2;
else
ac_local = a_test_local(idx_crit_local);
end
% Integrate
a_curr = a0_range(i);
N_curr = 0;
max_iter = 1e5;
iter = 0;
while a_curr < ac_local && a_curr < W/2 && iter < max_iter
[~, DK] = calculateStressIntensity(a_curr, sig_max, sig_min,
W);
if DK <= 0, break; end
dadN = C * (DK)^m;
da = dadN * dN_anim;
a_curr = a_curr + da;
N_curr = N_curr + dN_anim;
iter = iter + 1;
Page 5 of 11
end
Nf_matrix(i, j) = N_curr;
end
% Update progress bar
if mod(i, 5) == 0
fprintf('Progress: %.0f%%\n', i/length(a0_range)*100);
end
end
% Animate the contour plot buildup
[X, Y] = meshgrid(Delta_sigma_range, a0_range*1000);
for frameNum = 1:10
clf;
% Progressively increase contour levels
subplot('Position', [0.1, 0.1, 0.8, 0.8]);
% Start with fewer contours, end with more
contour_levels = round(linspace(5, 20, 10));
contourf(X, Y, Nf_matrix/1000, contour_levels(frameNum), 'LineStyle',
'none');
hold on;
% Add contour lines
contour(X, Y, Nf_matrix/1000, contour_levels(frameNum), 'k',
'LineWidth', 0.5);
colorbar;
xlabel('Stress Range, \Delta\sigma (MPa)', 'FontSize', 12);
ylabel('Initial Crack Size, a_0 (mm)', 'FontSize', 12);
title('Fatigue Life Design Map (×10³ cycles)', 'FontSize', 14,
'FontWeight', 'bold');
grid on;
% Add progress annotation
annotation('textbox', [0.7, 0.85, 0.25, 0.1], ...
'String', {sprintf('Parametric Study Animation'), ...
sprintf('Frame %d/10', frameNum)}, ...
'FontSize', 12, 'BackgroundColor', [0.95, 0.95, 1]);
frame = getframe(fig);
writeVideo(videoWriter, frame);
end
% Final detailed contour plot
for i = 1:5 % Hold final frame
clf;
subplot('Position', [0.1, 0.1, 0.8, 0.8]);
contourf(X, Y, Nf_matrix/1000, 20, 'LineStyle', 'none');
hold on;
contour(X, Y, Nf_matrix/1000, 15, 'k', 'LineWidth', 0.5);
% Highlight specific contours
[C2, h2] = contour(X, Y, Nf_matrix/1000, [10, 50, 100, 200], 'r',
'LineWidth', 2);
Page 6 of 11
clabel(C2, h2, 'FontSize', 10, 'Color', 'r');
colorbar;
xlabel('Stress Range, \Delta\sigma (MPa)', 'FontSize', 12);
ylabel('Initial Crack Size, a_0 (mm)', 'FontSize', 12);
title('Fatigue Life Design Map with Critical Contours', 'FontSize', 14,
'FontWeight', 'bold');
grid on;
% Add design implications
annotation('textbox', [0.15, 0.85, 0.3, 0.1], ...
'String', {'Design Implications:', ...
'• Lower stress = Longer life', ...
'• Smaller flaws = Safer design', ...
'• Red contours show critical zones'}, ...
'FontSize', 11, 'BackgroundColor', [1, 0.95, 0.95]);
frame = getframe(fig);
writeVideo(videoWriter, frame);
end
%% PART 5: COMPARISON OF MATERIALS ANIMATION
fig = figure('Position', [100, 100, 1200, 700], 'Color', 'white');
% Different material parameters
materials = {
'Aluminum 7075-T6', 6.9e-12, 3.0, [0, 0, 1];
'Titanium 6Al-4V', 3.5e-11, 2.5, [1, 0, 0];
'Steel A533B', 1.0e-11, 3.5, [0, 0.7, 0];
'Inconel 718', 2.0e-12, 3.8, [0.8, 0, 0.8]
};
% Animate each material
for matIdx = 1:size(materials, 1)
for frameNum = 1:15
clf;
% Calculate crack growth for this material
C_mat = materials{matIdx, 2};
m_mat = materials{matIdx, 3};
color_mat = materials{matIdx, 4};
a_curr = a0;
N_curr = 0;
a_hist = a0;
N_hist = 0;
while a_curr < ac && a_curr < W/2
[~, DK] = calculateStressIntensity(a_curr, sigma_max,
sigma_min, W);
if DK <= 0, break; end
dadN = C_mat * (DK)^m_mat;
da = dadN * dN_anim;
a_curr = a_curr + da;
N_curr = N_curr + dN_anim;
a_hist = [a_hist; a_curr];
N_hist = [N_hist; N_curr];
end
Page 7 of 11
% Plot all materials up to current one
subplot('Position', [0.1, 0.1, 0.55, 0.8]);
hold on; grid on; box on;
for i = 1:matIdx
% Recalculate for each material
C_temp = materials{i, 2};
m_temp = materials{i, 3};
color_temp = materials{i, 4};
a_temp = a0;
N_temp = 0;
a_temp_hist = a0;
N_temp_hist = 0;
while a_temp < ac && a_temp < W/2
[~, DK] = calculateStressIntensity(a_temp, sigma_max,
sigma_min, W);
if DK <= 0, break; end
dadN = C_temp * (DK)^m_temp;
da = dadN * dN_anim;
a_temp = a_temp + da;
N_temp = N_temp + dN_anim;
a_temp_hist = [a_temp_hist; a_temp];
N_temp_hist = [N_temp_hist; N_temp];
end
plot(N_temp_hist/1000, a_temp_hist*1000, 'Color', color_temp,
...
'LineWidth', 2 + (i == matIdx)*1, ...
'LineStyle', '-');
end
xlabel('Number of Cycles, N (×10³)', 'FontSize', 12);
ylabel('Crack Length, a (mm)', 'FontSize', 12);
title('Crack Growth in Different Materials', 'FontSize', 14,
'FontWeight', 'bold');
% Create legend
legend_strs = materials(1:matIdx, 1);
legend(legend_strs, 'Location', 'northwest', 'FontSize', 10);
yline(ac*1000, 'r--', 'Critical Crack Length', 'LineWidth', 2);
% Material properties display
subplot('Position', [0.7, 0.1, 0.25, 0.8]);
axis off;
y_pos = 0.8;
for i = 1:matIdx
text(0.1, y_pos, materials{i, 1}, ...
'Color', materials{i, 4}, 'FontSize', 12, 'FontWeight',
'bold');
text(0.1, y_pos-0.05, sprintf('C = %.1e, m = %.1f', ...
materials{i, 2}, materials{i, 3}), 'FontSize', 10);
y_pos = y_pos - 0.1;
end
% Highlight current material
if matIdx > 0
Page 8 of 11
annotation('rectangle', [0.68, 0.82-(matIdx-1)*0.1, 0.3, 0.08],
...
'EdgeColor', materials{matIdx, 4}, 'LineWidth', 2);
end
frame = getframe(fig);
writeVideo(videoWriter, frame);
end
end
%% PART 6: CONCLUSION AND SUMMARY
fig = figure('Position', [100, 100, 1200, 700], 'Color', 'white');
% Summary statistics
[min_Nf, min_idx] = min(Nf_matrix(:));
[min_i, min_j] = ind2sub(size(Nf_matrix), min_idx);
[max_Nf, max_idx] = max(Nf_matrix(:));
[max_i, max_j] = ind2sub(size(Nf_matrix), max_idx);
% Title
subplot('Position', [0.1, 0.7, 0.8, 0.2]);
axis off;
text(0.5, 0.7, 'ANALYSIS COMPLETE', ...
'HorizontalAlignment', 'center', 'FontSize', 28, 'FontWeight', 'bold');
text(0.5, 0.4, 'Key Findings and Design Recommendations', ...
'HorizontalAlignment', 'center', 'FontSize', 18, 'Color', [0.2, 0.2,
0.6]);
% Key findings
subplot('Position', [0.1, 0.4, 0.45, 0.25]);
axis off;
text(0, 0.9, 'KEY FINDINGS:', 'FontSize', 16, 'FontWeight', 'bold',
'Color', 'b');
text(0, 0.7, sprintf('1. Paris Law accurately predicts crack growth'),
'FontSize', 12);
text(0, 0.6, sprintf('2. Critical crack length: %.2f mm', ac*1000),
'FontSize', 12);
text(0, 0.5, sprintf('3. Base case fatigue life: %.0f cycles',
N_history(end)), 'FontSize', 12);
text(0, 0.4, '4. Crack growth accelerates near failure', 'FontSize', 12);
text(0, 0.3, '5. Material properties significantly affect life',
'FontSize', 12);
% Design recommendations
subplot('Position', [0.6, 0.4, 0.35, 0.25]);
axis off;
text(0, 0.9, 'DESIGN RECOMMENDATIONS:', 'FontSize', 16, 'FontWeight',
'bold', 'Color', 'r');
text(0, 0.7, '1. Minimize stress concentrations', 'FontSize', 12);
text(0, 0.6, '2. Implement regular inspections', 'FontSize', 12);
text(0, 0.5, '3. Use damage-tolerant design', 'FontSize', 12);
text(0, 0.4, '4. Consider material selection carefully', 'FontSize', 12);
text(0, 0.3, '5. Account for worst-case scenarios', 'FontSize', 12);
% Extreme cases
subplot('Position', [0.1, 0.1, 0.8, 0.2]);
axis off;
text(0, 0.8, 'EXTREME CASES IN PARAMETRIC STUDY:', 'FontSize', 14,
'FontWeight', 'bold');
Page 9 of 11
text(0, 0.5, sprintf('Worst Case: %.0f cycles (a₀=%.2f mm, Δσ=%.1f MPa)',
...
min_Nf, a0_range(min_i)*1000, Delta_sigma_range(min_j)), ...
'FontSize', 12, 'Color', 'r');
text(0, 0.2, sprintf('Best Case: %.0f cycles (a₀=%.2f mm, Δσ=%.1f MPa)',
...
max_Nf, a0_range(max_i)*1000, Delta_sigma_range(max_j)), ...
'FontSize', 12, 'Color', 'g');
% Credits
annotation('textbox', [0.7, 0.05, 0.25, 0.05], ...
'String', 'Created with MATLAB | Fracture Mechanics Simulation', ...
'FontSize', 10, 'HorizontalAlignment', 'right', ...
'BackgroundColor', [0.9, 0.9, 0.9]);
for i = 1:fps*3 % Hold final frame for 3 seconds
frame = getframe(fig);
writeVideo(videoWriter, frame);
end
%% CLOSE VIDEO
close(videoWriter);
fprintf('\nVideo creation complete!\n');
fprintf('Video saved as: %s\n', videoFileName);
fprintf('Total duration: %.1f seconds\n', videoDuration);
%% LOCAL FUNCTIONS (MUST BE AT THE END OF THE FILE)
function [K_max, Delta_K] = calculateStressIntensity(a, sigma_max,
sigma_min, W, geometry_factor)
Delta_sigma = sigma_max - sigma_min;
if nargin < 5
geometry_factor = sqrt(sec(pi * a / W));
end
K_max = sigma_max * geometry_factor * sqrt(pi * a);
Delta_K = Delta_sigma * geometry_factor * sqrt(pi * a);
end
function drawArrow(startPoint, endPoint, varargin)
% Custom function to draw an arrow
% startPoint: [x, y] starting point
% endPoint: [x, y] ending point
% Optional: 'Color', [r,g,b], 'LineWidth', width
% Parse optional arguments
color = [0, 0, 1]; % Default blue
lineWidth = 2;
for i = 1:2:length(varargin)
if strcmpi(varargin{i}, 'Color')
color = varargin{i+1};
elseif strcmpi(varargin{i}, 'LineWidth')
lineWidth = varargin{i+1};
elseif strcmpi(varargin{i}, 'Length')
% Calculate end point based on length
direction = endPoint - startPoint;
direction = direction / norm(direction);
endPoint = startPoint + direction * varargin{i+1};
Page 10 of 11
end
end
% Draw the arrow line
plot([startPoint(1), endPoint(1)], [startPoint(2), endPoint(2)], ...
'Color', color, 'LineWidth', lineWidth);
% Calculate arrowhead
arrowLength = norm(endPoint - startPoint);
headLength = min(0.2, arrowLength * 0.2);
headWidth = headLength * 0.6;
direction = (endPoint - startPoint) / arrowLength;
perpendicular = [-direction(2), direction(1)];
% Arrowhead points
p1 = endPoint;
p2 = endPoint - direction * headLength + perpendicular * headWidth;
p3 = endPoint - direction * headLength - perpendicular * headWidth;
% Draw arrowhead
patch([p1(1), p2(1), p3(1)], [p1(2), p2(2), p3(2)], ...
color, 'EdgeColor', color);
end
Page 11 of 11