% This function calculates the dislocation density for a slip zone of size % d in a case detailed by parameter using the bounded bounded quadrature function [Phi,om,W,t,s,Tol]=CalculatePhiSeq(parameter,d) % First eigenvalue N = parameter.N; % Number of Gauss along slip zone ID = parameter.ID; % Logical vector identifying the equation to be removed f = parameter.f; % Coefficient of friction lam1 = parameter.lam1; % Mode I eigenvalue lam2 = parameter.lam2; % Mode II eigenvalue g1 = parameter.g1; % Mode I shear multiplier g2 = parameter.g2; % Mode II shear multiplier % Numerical quadrature for bounded-bounded idx = 1:N+1; % Index up to number of points t = cos(pi.*(2.*idx'-1)/(2*(N+1))); % Observation (collocation) points s = cos(pi.*idx(1:end-1)./(N+1)); % Integration points W = (1-s.^2)./(N+1); % Weights om = (1-s.^2).^(1/2); % Fundamental function [Fxyy,Fxxy,~,~] = WedgeKernels(t(~ID),s,parameter,d); % Kernel functions A = pi*W.*(-f*Fxyy + Fxxy + 1./(t(~ID)-s)); % Matrix A - replace last term with polar kernel? % Right hand side for slip zone at contact edge F = (f-g1).*(d./2.*(t(~ID)+1)).^(lam1-1) ... + (f-g2).*(d./2.*(t(~ID)+1)).^(lam2-1); % right hand side % Calculate phi values Phi = A\F; %% Calculating stresses at removed point [Fxyy,Fxxy,~,~] = WedgeKernels(t(ID),s,parameter,d); % Kernel functions AID = pi*W.*(-f*Fxyy + Fxxy + 1./(t(ID)-s)); % Matrix A - replace last term with polar kernel? % Right hand side for slip zone at contact edge FID = (f-g1).*(d./2.*(t(ID)+1)).^(lam1-1) ... + (f-g2).*(d./2.*(t(ID)+1)).^(lam2-1); % right hand side Tol=AID*Phi-FID; end % Kernelfunctions % % Kernels for cracks evolving from apex of 3/4 plane % -> origin at the surface % -> dislocation at y=0 (hard coded) % Date: 08/10/20 % function [Fxyy,Fxxy,Fyyy,Fyxy] = WedgeKernels(t,s,parameter,d) % Parameter N = parameter.N; theta = parameter.theta; n = parameter.n; facl = parameter.facl; alpha = parameter.alpha; %------------- Prepare numerical integration for free surfaces ------------ idx = 1:n; % Index up to number of points sn = cos(pi.*2.*idx./(2*n+1)); % Integration points tn = cos(pi.*(2.*idx'-1)./(2*n+1)); % Collocation points Wn = 2.*(1-sn)./(2*n+1); % Weights x = d./2.*(t+1); xi = d./2.*(s+1); a = xi; % Radius of dislocation l = facl*d; % Length of free surfaces t = 1-2.^(1-2.*x./l); s = 1-2.^(1-2.*a./l); %-------------------------------------------------------------------------- %------------------------------ Calculation ------------------------------- %-------------------------------------------------------------------------- % Assemble matrix A A = zeros(4*n,4*n); [Krtt,Krrt,Kttt,Ktrt] = PolarKernels_norm(tn,alpha,sn,alpha); A(1:n ,1:4:end) = pi.*Wn.*Krtt; A(1:n ,2:4:end) = pi.*Wn.*Kttt; A(n+1:2*n,1:4:end) = pi.*Wn.*Krrt; A(n+1:2*n,2:4:end) = pi.*Wn.*Ktrt; [Krtt,Krrt,Kttt,Ktrt] = PolarKernels_norm(tn,alpha,sn,-alpha); A(1:n ,3:4:end) = pi.*Wn.*Krtt; A(1:n ,4:4:end) = pi.*Wn.*Kttt; A(n+1:2*n,3:4:end) = pi.*Wn.*Krrt; A(n+1:2*n,4:4:end) = pi.*Wn.*Ktrt; [Krtt,Krrt,Kttt,Ktrt] = PolarKernels_norm(tn,-alpha,sn,alpha); A(2*n+1:3*n,1:4:end) = pi.*Wn.*Krtt; A(2*n+1:3*n,2:4:end) = pi.*Wn.*Kttt; A(3*n+1:4*n,1:4:end) = pi.*Wn.*Krrt; A(3*n+1:4*n,2:4:end) = pi.*Wn.*Ktrt; [Krtt,Krrt,Kttt,Ktrt] = PolarKernels_norm(tn,-alpha,sn,-alpha); A(2*n+1:3*n,3:4:end) = pi.*Wn.*Krtt; A(2*n+1:3*n,4:4:end) = pi.*Wn.*Kttt; A(3*n+1:4*n,3:4:end) = pi.*Wn.*Krrt; A(3*n+1:4*n,4:4:end) = pi.*Wn.*Ktrt; % Right hand side for radial and tangential dislocation Fr = zeros(4*n,N); Ft = zeros(4*n,N); [Krtt,Krrt,Kttt,Ktrt] = PolarKernels_norm(tn,alpha,s,theta); % Kernel functions Fr(1:n,:) = -Krtt./(l./2./((1-s).*log(2))); % Normal traction at pos half-line Fr(n+1:2*n,:) = -Krrt./(l./2./((1-s).*log(2))); % Shear traction at pos half-line Ft(1:n,:) = -Kttt./(l./2./((1-s).*log(2))); % Normal traction at pos half-line Ft(n+1:2*n,:) = -Ktrt./(l./2./((1-s).*log(2))); % Shear traction at pos half-line [Krtt,Krrt,Kttt,Ktrt] = PolarKernels_norm(tn,-alpha,s,theta); % Kernel functions Fr(2*n+1:3*n,:) = -Krtt./(l./2./((1-s).*log(2))); % Normal traction at neg half-line Fr(3*n+1:4*n,:) = -Krrt./(l./2./((1-s).*log(2))); % Shear traction at neg half-line Ft(2*n+1:3*n,:) = -Kttt./(l./2./((1-s).*log(2))); % Normal traction at neg half-line Ft(3*n+1:4*n,:) = -Ktrt./(l./2./((1-s).*log(2))); % Shear traction at neg half-line [Krttp,Krrtp,Ktttp,Ktrtp] = PolarKernels_norm(t,theta,sn,alpha); [Krttn,Krrtn,Ktttn,Ktrtn] = PolarKernels_norm(t,theta,sn,-alpha); %---------------------------- Radial direction ---------------------------- phi = A\Fr; % Extract different distributions of dislocations phirp = phi(1:4:end,:); % Density at pos free surf in rad direction phitp = phi(2:4:end,:); % Density at pos free surf in tan direction phirn = phi(3:4:end,:); % Density at neg free surf in rad direction phitn = phi(4:4:end,:); % Density at neg free surf in tan direction %-------------------------------- Kernels --------------------------------- % Fxxx = pi*(Wn.*Krrrp*phirp + Wn.*Ktrrp*phitp + Wn.*Krrrn*phirn + Wn.*Ktrrn*phitn)*c/2; Fxyy = pi*(Wn.*Krttp*phirp + Wn.*Ktttp*phitp + Wn.*Krttn*phirn + Wn.*Ktttn*phitn)*d/2; Fxxy = pi*(Wn.*Krrtp*phirp + Wn.*Ktrtp*phitp + Wn.*Krrtn*phirn + Wn.*Ktrtn*phitn)*d/2; %-------------------------- Tangential direction -------------------------- phi = A\Ft; % Extract different distributions of dislocations phirp = phi(1:4:end,:); % Density at pos free surf in rad direction phitp = phi(2:4:end,:); % Density at pos free surf in tan direction phirn = phi(3:4:end,:); % Density at neg free surf in rad direction phitn = phi(4:4:end,:); % Density at neg free surf in tan direction %-------------------------------- Kernels --------------------------------- %-------------------------------- Kernels --------------------------------- % Fyxx = pi*(Wn.*Krrrp*phirp + Wn.*Ktrrp*phitp + Wn.*Krrrn*phirn + Wn.*Ktrrn*phitn)*c/2; Fyyy = pi*(Wn.*Krttp*phirp + Wn.*Ktttp*phitp + Wn.*Krttn*phirn + Wn.*Ktttn*phitn)*d/2; Fyxy = pi*(Wn.*Krrtp*phirp + Wn.*Ktrtp*phitp + Wn.*Krrtn*phirn + Wn.*Ktrtn*phitn)*d/2; end % Kernelfunctions % % Kernels for dislocations near a notch % -> origin at the tip of the notch % -> dislocation at r = a and theta = phi % Checked against StressTrans_finalDan.m % % Date: 23/09/20 % % function [Krrr,Krtt,Krrt,Ktrr,Kttt,Ktrt] = PolarKernels_norm(t,theta,s,phi) function [Krtt,Krrt,Kttt,Ktrt] = PolarKernels_norm(t,theta,s,phi) % Krrr = (s-1).^(-1).*(log((-2).*(s-1).^(-1)).^2-2.*cos(theta-1.*phi) ... % .*log((-2).*(s-1).^(-1)).*log((-2).*(t-1).^(-1))+log((-2).*( ... % t-1).^(-1)).^2).^(-2).*(2.*cos(theta-1.*phi).*log((-2).*(s-1) ... % .^(-1)).^3-1.*(2+cos(2.*(theta-1.*phi))).*log((-2).*(s-1).^(-1) ... % ).^2.*log((-2).*(t-1).^(-1))+log((-2).*(t-1).^(-1)).^3).* ... % sin(theta-1.*phi); Krtt = (s-1).^(-1).*(log((-2).*(s-1).^(-1)).^2-2.*cos(theta-1.*phi) ... .*log((-2).*(s-1).^(-1)).*log((-2).*(t-1).^(-1))+log((-2).*( ... t-1).^(-1)).^2).^(-2).*((4+cos(2.*(theta-1.*phi))).*log((-2).*(( ... -1)+s).^(-1)).^2.*log((-2).*(t-1).^(-1))+log((-2).*(t-1).^( ... -1)).^3-2.*cos(theta-1.*phi).*log((-2).*(s-1).^(-1)).*(log((-2) ... .*(s-1).^(-1)).^2+2.*log((-2).*(t-1).^(-1)).^2)).*sin(theta-1 ... .*phi); Krrt = (s-1).^(-1).*(log((-2).*(s-1).^(-1))-1.*cos(theta-phi).* ... log((-2).*(t-1).^(-1))).*(log((-2).*(s-1).^(-1)).^2-2.* ... cos(theta-phi).*log((-2).*(s-1).^(-1)).*log((-2).*(t-1).^( ... -1))+log((-2).*(t-1).^(-1)).^2).^(-2).*(cos(2.*(theta-phi)).* ... log((-2).*(s-1).^(-1)).^2+log((-2).*(t-1).^(-1)).*((-2).* ... cos(theta-phi).*log((-2).*(s-1).^(-1))+log((-2).*(t-1).^(-1) ... ))); % Ktrr = (1/2).*(s-1).^(-1).*(log((-2).*(s-1).^(-1)).^2-2.*cos(theta+( ... % -1).*phi).*log((-2).*(s-1).^(-1)).*log((-2).*(t-1).^(-1))+log( ... % (-2).*(t-1).^(-1)).^2).^(-2).*(2.*log((-2).*(s-1).^(-1)).^3+ ... % ((-7).*cos(theta-1.*phi)+cos(3.*(theta-1.*phi))).*log((-2).*(s-1).^( ... % -1)).^2.*log((-2).*(t-1).^(-1))+6.*log((-2).*(s-1).^(-1)).* ... % log((-2).*(t-1).^(-1)).^2-2.*cos(theta-1.*phi).*log((-2).*((-1) ... % +t).^(-1)).^3); Kttt = (1/2).*(s-1).^(-1).*(log((-2).*(s-1).^(-1)).^2-2.*cos(theta+( ... -1).*phi).*log((-2).*(s-1).^(-1)).*log((-2).*(t-1).^(-1))+log( ... (-2).*(t-1).^(-1)).^2).^(-2).*(2.*log((-2).*(s-1).^(-1)).^3+ ... (-1).*(5.*cos(theta-1.*phi)+cos(3.*(theta-1.*phi))).*log((-2).*(s-1) ... .^(-1)).^2.*log((-2).*(t-1).^(-1))+2.*(1+2.*cos(2.*(theta-1.*phi)) ... ).*log((-2).*(s-1).^(-1)).*log((-2).*(t-1).^(-1)).^2-2.* ... cos(theta-1.*phi).*log((-2).*(t-1).^(-1)).^3); Ktrt = (-1).*(s-1).^(-1).*log((-2).*(t-1).^(-1)).*(log((-2).*((-1)+ ... s).^(-1)).^2-2.*cos(theta-1.*phi).*log((-2).*(s-1).^(-1)).*log( ... (-2).*(t-1).^(-1))+log((-2).*(t-1).^(-1)).^2).^(-2).*(cos( ... 2.*(theta-1.*phi)).*log((-2).*(s-1).^(-1)).^2+log((-2).*(t-1) ... .^(-1)).*((-2).*cos(theta-1.*phi).*log((-2).*(s-1).^(-1))+log((-2) ... .*(t-1).^(-1)))).*sin(theta-1.*phi); end