%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % This function numerically computes the stress field for a collection of % arrays of dislocations with the specified Burgers vectors and positions. % % Arguments: % - X, Y: computational grids (as generated by meshgrid()) % - num_dislocation_arrays: number of dislocation arrays % - burgers_vectors: matrix containing Burgers vectors for % dislocation arrays. The i-th row should % contain the Burgers vector for the i-th % dislocation array % - positions: matrix containing positions for dislocation % arrays. The i-th row should contain the (x,y) % coordinates of the i-th dislocation array % - numerical_core_radius: radius to use for numerical core % - delta_function_type: type of delta function to use when smearing % core. Choices are: % 1: first-order trigonometric approximation % to delta function (see Osher & Fedkiw 2003) % 2: second-order trigonometric approximation % to delta function % 3: second-order trigonometric approximation % to delta function % 4: modified Lorentzian approximation % to delta function % - G: shear modulus % - poisson_ratio: Poisson ratio % % Return values: % - sigma_xx: xx-component of stress field % - sigma_yy: yy-component of stress field % - sigma_xy: xy-component of stress field % - d_z: dislocation line field for last % dislocation line array % % % NOTES: % - The line direction of all dislocations is assumed to be in the (0,0,1) % direction. % - Each array is oriented such that the normal to the plane of the % dislocation array is in the (1,0,0) direction and all dislocation % arrays are periodically repeated in the x-direction. % - The stress fields for the dislocation configuration are computed using % the Fourier transform formulation in Xiang et. al. % - The grid is assumed to be uniform in each spatial coordinate direction. % % Kevin T. Chu % MAE, Princeton University % 08/2007 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [sigma_xx, sigma_yy, sigma_xy, d_z] = ... compute_numerical_stress_fields(X, Y, ... num_dislocation_arrays, ... burgers_vectors, positions, ... numerical_core_radius, ... delta_function_type, ... G, poisson_ratio) % check arguments if (num_dislocation_arrays ~= size(burgers_vectors,1)) error('Number of Burgers vectors does not match number of dislocation arrays'); end if (num_dislocation_arrays ~= size(positions,1)) error('Number of positions does not match number of dislocation arrays'); end % compute Nx and Ny [Nx,Ny] = size(X); % compute Lx and Ly Lx = Nx*(X(2,1)-X(1,1)); Ly = Ny*(Y(1,2)-Y(1,1)); % initialize solution for stress sigma_xx = zeros(size(X)); sigma_yy = zeros(size(X)); sigma_xy = zeros(size(X)); % compute stress field by summing up the stress field from % individual dislocations for line_num = 1:num_dislocation_arrays % get Burgers vector and position of current dislocation line b = burgers_vectors(line_num,:); pos = positions(line_num,:); switch (delta_function_type) case 1 % compute delta function using first-order delta function delta_X = zeros(size(X)); x_in_core = find(abs(X - pos(1)) <= numerical_core_radius); delta_X(x_in_core) = 0.5/numerical_core_radius ... * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius)); delta_Y = zeros(size(Y)); y_in_core = find(abs(Y - pos(2)) <= numerical_core_radius); delta_Y(y_in_core) = 0.5/numerical_core_radius ... * (1+cos(pi*(Y(y_in_core)-pos(2))/numerical_core_radius)); case 2 % compute delta functions using second-order delta function delta_X = zeros(size(X)); x_in_core = find(abs(X - pos(1)) <= numerical_core_radius); delta_X(x_in_core) = -1/6/numerical_core_radius ... * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius)); x_in_core = find(abs(X - pos(1)) <= 0.5*numerical_core_radius); delta_X(x_in_core) = 4/3/numerical_core_radius ... * (1+cos(2*pi*(X(x_in_core)-pos(1))/numerical_core_radius))... -1/6/numerical_core_radius ... * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius)); delta_Y = zeros(size(Y)); y_in_core = find(abs(Y - pos(2)) <= numerical_core_radius); delta_Y(y_in_core) = -1/6/numerical_core_radius ... * (1+cos(pi*(Y(y_in_core)-pos(2))/numerical_core_radius)); y_in_core = find(abs(Y - pos(2)) <= 0.5*numerical_core_radius); delta_Y(y_in_core) = 4/3/numerical_core_radius ... * (1+cos(2*pi*(Y(y_in_core)-pos(2))/numerical_core_radius))... -1/6/numerical_core_radius ... * (1+cos(pi*(Y(y_in_core)-pos(2))/numerical_core_radius)); case 3 % compute delta functions using third-order delta function delta_X = zeros(size(X)); x_in_core = find(abs(X - pos(1)) <= numerical_core_radius); delta_X(x_in_core) = 1/48/numerical_core_radius ... * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius)); x_in_core = find(abs(X - pos(1)) <= 0.5*numerical_core_radius); delta_X(x_in_core) = -16/15/numerical_core_radius ... * (1+cos(2*pi*(X(x_in_core)-pos(1))/numerical_core_radius))... +1/48/numerical_core_radius ... * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius)); x_in_core = find(abs(X - pos(1)) <= numerical_core_radius/3); delta_X(x_in_core) = 243/80/numerical_core_radius ... * (1+cos(3*pi*(X(x_in_core)-pos(1))/numerical_core_radius))... -16/15/numerical_core_radius ... * (1+cos(2*pi*(X(x_in_core)-pos(1))/numerical_core_radius))... +1/48/numerical_core_radius ... * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius)); delta_Y = zeros(size(Y)); y_in_core = find(abs(Y - pos(2)) <= numerical_core_radius); delta_Y(y_in_core) = 1/48/numerical_core_radius ... * (1+cos(pi*(Y(y_in_core)-pos(1))/numerical_core_radius)); y_in_core = find(abs(Y - pos(2)) <= 0.5*numerical_core_radius); delta_Y(y_in_core) = -16/15/numerical_core_radius ... * (1+cos(2*pi*(Y(y_in_core)-pos(2))/numerical_core_radius))... +1/48/numerical_core_radius ... * (1+cos(pi*(Y(y_in_core)-pos(2))/numerical_core_radius)); y_in_core = find(abs(Y - pos(2)) <= numerical_core_radius/3); delta_Y(y_in_core) = 243/80/numerical_core_radius ... * (1+cos(3*pi*(Y(y_in_core)-pos(2))/numerical_core_radius))... -16/15/numerical_core_radius ... * (1+cos(2*pi*(Y(y_in_core)-pos(2))/numerical_core_radius))... +1/48/numerical_core_radius ... * (1+cos(pi*(Y(y_in_core)-pos(2))/numerical_core_radius)); case 4 % compute delta functions using a modified Lorentzian function delta_X = 1/(2*(atan(pi)/pi-1/(1+pi^2))) ... * 1/numerical_core_radius ... * ( 1./(1+(pi*(X-pos(1))/numerical_core_radius).^2) ... - 1/(1+pi^2) ); x_outside_core = find(abs(X - pos(1)) > numerical_core_radius); delta_X(x_outside_core) = 0; delta_Y = 1/(2*(atan(pi)/pi-1/(1+pi^2))) ... * 1/numerical_core_radius ... * ( 1./(1+(pi*(Y-pos(2))/numerical_core_radius).^2) ... - 1/(1+pi^2) ); y_outside_core = find(abs(Y - pos(2)) > numerical_core_radius); delta_Y(y_outside_core) = 0; otherwise error('Unsupported order for delta function'); end % compute dislocation line field d_z = delta_X.*delta_Y; % compute FFT of d_z d_z_fft = fft2(d_z); % compute stress field from current dislocation in Fourier space k_x = 0:Nx-1; k_x(ceil((Nx+1)/2):end) = k_x(ceil((Nx+1)/2):end) - Nx; k_x = 2*pi*k_x/Lx; k_y = 0:Ny-1; k_y(ceil((Ny+1)/2):end) = k_y(ceil((Ny+1)/2):end) - Ny; k_y = 2*pi*k_y/Ly; [K_y,K_x] = meshgrid(k_y, k_x); norm_K_sq = K_x.^2 + K_y.^2; sigma_xx_fft = -2*G*i/(1-poisson_ratio)*K_y.^2./norm_K_sq.^2 ... .*(K_x*b(2) - K_y*b(1)).*d_z_fft; sigma_yy_fft = -2*G*i/(1-poisson_ratio)*K_x.^2./norm_K_sq.^2 ... .*(K_x*b(2) - K_y*b(1)).*d_z_fft; sigma_xy_fft = 2*G*i/(1-poisson_ratio)*K_x.*K_y./norm_K_sq.^2 ... .*(K_x*b(2) - K_y*b(1)).*d_z_fft; % correct DC component of stress fields in Fourier space sigma_xx_fft(1,1) = 0; sigma_yy_fft(1,1) = 0; sigma_xy_fft(1,1) = 0; % compute stress contribution in real space sigma_xx_cur = ifft2(sigma_xx_fft); sigma_yy_cur = ifft2(sigma_yy_fft); sigma_xy_cur = ifft2(sigma_xy_fft); % check that imaginary part is small if ( max(max(abs(imag(sigma_xx_cur)))) > 1e-10 ... | max(max(abs(imag(sigma_yy_cur)))) > 1e-10 ... | max(max(abs(imag(sigma_xy_cur)))) > 1e-10 ) max(max(abs(imag(sigma_xx_cur)))) max(max(abs(imag(sigma_yy_cur)))) max(max(abs(imag(sigma_xy_cur)))) warning('Imaginary component of stress nonzero!!'); end sigma_xx = sigma_xx + real(sigma_xx_cur); sigma_yy = sigma_yy + real(sigma_yy_cur); sigma_xy = sigma_xy + real(sigma_xy_cur); end