%=== PROBLEM DEFINITION (Loading and Model Parameters) ====== dN=1; % transition cycle interval dtMeasu=1000; % measurement cycle interval dSmin=5; % min. load sNomin=65; % nominal load sOverl=[125 100 100 125 125]; % overload cNomin=[5 5 5 5 45]; % num. of cycle in a block: nominal load cOverl=[45 45 95 95 5]; % num. of cycle in a block: overload cBlock=[5000 10000 15000 20000 25000];% cycles for each block a0=0.01; % true model param including initi crack, a0 m=3.1; C=log(5.5e-11); dKth=5.2; beta=0.2; n=2.8; sY=580; thick=4e-3; width=150e-3/2; % geometry dimension %============================================================ cycle=[0:dtMeasu:cBlock(end)]'; % % % LOAD RATIO, R nb=length(cBlock); cInter(1)=cBlock(1); if nb>1; for i=2:nb; cInter(i)=cBlock(i)-cBlock(i-1); end;end dSmax=[]; for i=1:nb; B=[sNomin*ones(cNomin(i),1); sOverl(i)*ones(cOverl(i),1)]; dSmax(length(dSmax)+1:cBlock(i),1)= ... repmat(B,cInter(i)/(cNomin(i)+cOverl(i)),1); end; dSmax=dSmax(1:dN:end); R=dSmin./dSmax; % % % HUANG's MODEL aCurre(1)=a0; akOL(1)=a0; rOL(1)=0; for k=1:max(cycle)/dN; ak=aCurre(k); Y=(1-ak./(2*width)+0.326.*(ak./width).^2)...%geometry factor ./sqrt(1-ak./width); loca=ak./width>1; Y(loca)=0; Kmax=Y.*dSmax(k).*sqrt(pi*ak); %SIF from Eq. 6.5 Kmin=Y.*dSmin.*sqrt(pi*ak); % Scaling factor, Mr (Eq. 6.6) Mr=(1-R(k)).^(-beta); % Correction factor, Mp (Eqs. 6.7 and 6.8) alpha=0.35-0.29./(1+(1.08.*Kmax.^2./(thick.*sY.^2)).^2.15); ry=alpha.*(Kmax./sY).^2; if R(k)~=dSmin/sNomin; loca=ry>rOL; rOL(loca)=ry(loca); akOL(loca)=ak(loca); end Mp=(ry./(akOL+rOL-ak)).^n; loca=ak+ry>=akOL+rOL; Mp(loca)=1; % Calculate crack size dKeq=Mr.*Mp.*(Kmax-Kmin); % Eq. 6.5 rate=exp(C).*(dKeq.^m - dKth.^m); % Eq. 6.4 loca=dKeq<=dKth; rate(loca)=0; aCurre(k+1)=aCurre(k)+rate.*dN; end; crackHuang=aCurre(1:dtMeasu/dN:max(cycle)/dN+1); % % % DATA PLOT plot(cycle,crackHuang);