%int_defects_R6_analysis_ed1.m
%Harry Coules 2018
%
%DESCRIPTION
%This script plot two sets of fracture analysis results in R6-type Failure
%Assessment Diagrams on the same set of axes. The first assessment point is
%the phi=19.7 degree point in a crack in an interacting crack pair. The
%second assessment point is the deepest point (phi=90 degrees) in a single
%crack. The resulting FAD is shown in Figure 3 of the accompanying paper:
%   H. E. Coules, "On predicting the interaction of crack-like defects in
%      ductile fracture", 2018
%
%INPUT
%Before running this script, TUBE_6_R6_analysis_ed1.mat must be loaded into
%the workspace. This contains FEA results (contour integral results,
%limit loads) in addition to material properties etc.

%% Material
%Calculate flow stress
matProps.s_f = (matProps.s_y+matProps.s_uts)/2; %Mean of yield and UTS

%FAC cut-off
matProps.Lr_max = matProps.s_f/matProps.s_y;    %Flow stress/yield stress

%% Calculate Lr and Kr
%Lr
fractSingle.Lr = p./fractSingle.pL;
fractd025.Lr = p./fractd025.pL;

%Kr
fractSingle.Kr_90 = fractSingle.K_90./matProps.KIC;
fractd025.Kr = fractd025.K./matProps.KIC;

%Lrs used for plotting
plotPressures = [0:20:120]';
fractSingle.plotLrs = plotPressures./fractSingle.pL;
fractd025.plotLrs = plotPressures./fractd025.pL;

%Krs used for plotting
fractSingle.plotKrs_90 = interp1(fractSingle.Lr,fractSingle.Kr_90,fractSingle.plotLrs);
fractd025.plotKrs = interp1(fractd025.Lr,fractd025.Kr,fractd025.plotLrs);

%% Calculate FACs
%Option 1 (R6 eq. I.6.3). Note any Lr vector could be used - but must be consistent when plotting.
FACs.f1 = ((1+0.5*(fractSingle.Lr.^2)).^-0.5).*(0.3+0.7*exp(-0.6*(fractSingle.Lr.^6)));
FACs.f1(isnan(FACs.f1))=1;

%Option 2 (R6 eq. I.6.4). Note any Lr vector could be used - but must be consistent when plotting.
FACs.e_ref = interp1(matProps.plasticTable(:,1),matProps.elasticTable(:,2)+matProps.plasticTable(:,2),fractSingle.Lr*matProps.s_y);
FACs.f2 = (((matProps.E*FACs.e_ref)./(fractSingle.Lr*matProps.s_y))+(((fractSingle.Lr.^3)*matProps.s_y)./(2*matProps.E*FACs.e_ref))).^-0.5;
FACs.f2(isnan(FACs.f2))=1;

%Option 3 FAC (Interacting cracks)
FACs.f3 = fractd025.K./fractd025.KJ;
FACs.f3(isnan(FACs.f3))=1;

%Option 3 FAC (Single crack)
FACs.f3S_90 = fractSingle.K_90./fractSingle.KJ_90;
FACs.f3S_90(isnan(FACs.f3S_90))=1;

%% Plot FAD
f1 = figure('pos',[10,10,1200,600]);
hold on

%Assessment lines
h_al_co = line([matProps.Lr_max,matProps.Lr_max],[0 1.2],'linewidth',2,'linestyle','-','color','k');    %Cutoff line
h_al_f1 = plot(fractSingle.Lr,FACs.f1, 'linewidth', 2, 'linestyle', '-.', 'color', 'k');   %Option 1 - note Lr used consistent with calc.
h_al_f2 = plot(fractSingle.Lr,FACs.f2, 'linewidth', 2, 'linestyle', '--', 'color', 'k');   %Option 2 - note Lr used consistent with calc.
h_al_f3S = plot(fractSingle.Lr,FACs.f3S_90, 'linewidth', 2, 'linestyle', '-', 'color', [0,0.4470,0.7410]);  %Option 3 (single)
h_al_f3 = plot(fractd025.Lr,FACs.f3, 'linewidth', 2, 'linestyle', '-', 'color', [0.6350,0.0780,0.1840]);     %Option 3 (interacting)

%Loading lines
h_ll_S = plot(fractSingle.plotLrs,fractSingle.plotKrs_90,'linewidth',1,'color', [0,0.4470,0.7410],'marker','o','markerfacecolor','w');    %Single
h_ll_d025 = plot(fractd025.plotLrs,fractd025.plotKrs, 'linewidth',1, 'color', [0.6350,0.0780,0.1840],'marker','o','markerfacecolor','w');    %Interacting

%Axis properties
axis equal
xlim([0 1.6]);
ylim([0 1.2]);
xlabel('L_r')
ylabel('K_r')

grid on
box on
set(gca,'linewidth',1,'fontsize',14)

%Legend
l1 = legend([h_al_f1,h_al_f2,h_al_f3S,h_al_f3,h_al_co,h_ll_S,h_ll_d025],...
{'Option 1 FAC',...
   'Option 2 FAC',...
   'Option 3 FAC (single)',...
   'Option 3 FAC (twin, d = 0.25t)',...
   'FAC cut-off line',...
   'Pressurisation (single)',...
   'Pressurisation (twin, d = 0.25t)'});
set(l1,'location','eastoutside')