- #1
William White
- 256
- 79
Hi,
I am modelling the stress field of a dislocation (an extra half line of atoms in a crystal) in a 2D plane using MatLab.
The dislocation can line at different angles. I can rotate the dislocation for any angle; but in relaity the angle is limited by the structure of the crystal: theta is either 0°, 54.7° or 125.2°. i.e. the dislocation must lie on a crystal plane and these are the angles allowed.
My code sets up material properties u; nu and b
They sets up a grid. This gives values of x and y in the 2D plane at which the stress is calculated
The next lines define the magnitude of the stresses Qxx, Qty and Qxy at all values of x and y.
Now if I want to rotate the dislocation, by say 54.7°; x and y are redefined as x2 and y2 and the stresses recalculated as QxxPRIME etc.
The next step of the problem is to move the dislocation by a distance, s, along the plane it sits on. If theta, is say, 0°, which aligns conveniently with the x-axis, then the dislocation should move along the x-axis by a distance s.
If theta is 90°, it moves along the y-axis by a distance s
If theta is 54.7°, it moves along a plane inclined at 54.7, a distance s.
My problem is getting the dislocation to move correctly. The origin 0,0 must stay where it is. The dislocation is centered on the origin. So when x=y=0 then the stresses are infinite - that is the stress at the centre of the dislocation is infinite.
So I am guessing what I cannot do is just add values to x and y, because then x and y will not be zero at the centre of the dislocation.(The goal of the exercise is to have several (in the end, tens of thousands) of dislocations moving along 1D planes in a 2D area and to calculate the entire stress at any point at any time.)If one runs the code twice, with angle at 54.7 and again at 0, one can see the effect of the rotation, which I am happy with.Thank youu = 80000; % Shear modulus, measure in Mpa
nu = 0.3; % Poission’s ratio of material
b = .00025; % Burgers Vector measured in microns
[x,y]=meshgrid(-100:1:100,-100:1:100); % area of grid to which stress field is mapped
angle = 54.7; % angle of slip plane in degrees
theta = angle*pi/180; % angle converted to radians
Qxx = -(u .* b / (2.*pi .* (1-nu))) .* ((y .* (3.*x.^2 +y.^2)) ./ ((x.^2 + y.^2).^2)); % Qxx stress
Qyy = (u .* b / (2.*pi .* (1-nu))) .* ((y .* (x.^2 -y.^2)) ./ ((x.^2 + y.^2).^2)); % Qyy stress
Qxy = (u .* b / (2.*pi .* (1-nu))) .* ((x .* (x.^2 -y.^2)) ./ ((x.^2 + y.^2).^2)); % Qxy stressx2 = x.*cos(theta) + y.*sin(theta); % rotates dislocation
y2 = y.*cos(theta) - x.*sin(theta); % rotates dislocationQxxPRIME = -(u .* b / (2.*pi .* (1-nu))) .* ((y2 .* (3.*x2.^2 +y2.^2)) ./ ((x2.^2 + y2.^2).^2)); % Qyy stress PRIME
QyyPRIME = (u .* b / (2.*pi .* (1-nu))) .* ((y2 .* (x2.^2 -y2.^2)) ./ ((x2.^2 + y2.^2).^2)); % Qyy stress PRIME
QxyPRIME = (u .* b / (2.*pi .* (1-nu))) .* ((x2 .* (x2.^2 -y2.^2)) ./ ((x2.^2 + y2.^2).^2)); % Qxy stress PRIME
contourf(x(1,:),y(:,1)',QxxPRIME)
axis equal
I am modelling the stress field of a dislocation (an extra half line of atoms in a crystal) in a 2D plane using MatLab.
The dislocation can line at different angles. I can rotate the dislocation for any angle; but in relaity the angle is limited by the structure of the crystal: theta is either 0°, 54.7° or 125.2°. i.e. the dislocation must lie on a crystal plane and these are the angles allowed.
My code sets up material properties u; nu and b
They sets up a grid. This gives values of x and y in the 2D plane at which the stress is calculated
The next lines define the magnitude of the stresses Qxx, Qty and Qxy at all values of x and y.
Now if I want to rotate the dislocation, by say 54.7°; x and y are redefined as x2 and y2 and the stresses recalculated as QxxPRIME etc.
The next step of the problem is to move the dislocation by a distance, s, along the plane it sits on. If theta, is say, 0°, which aligns conveniently with the x-axis, then the dislocation should move along the x-axis by a distance s.
If theta is 90°, it moves along the y-axis by a distance s
If theta is 54.7°, it moves along a plane inclined at 54.7, a distance s.
My problem is getting the dislocation to move correctly. The origin 0,0 must stay where it is. The dislocation is centered on the origin. So when x=y=0 then the stresses are infinite - that is the stress at the centre of the dislocation is infinite.
So I am guessing what I cannot do is just add values to x and y, because then x and y will not be zero at the centre of the dislocation.(The goal of the exercise is to have several (in the end, tens of thousands) of dislocations moving along 1D planes in a 2D area and to calculate the entire stress at any point at any time.)If one runs the code twice, with angle at 54.7 and again at 0, one can see the effect of the rotation, which I am happy with.Thank youu = 80000; % Shear modulus, measure in Mpa
nu = 0.3; % Poission’s ratio of material
b = .00025; % Burgers Vector measured in microns
[x,y]=meshgrid(-100:1:100,-100:1:100); % area of grid to which stress field is mapped
angle = 54.7; % angle of slip plane in degrees
theta = angle*pi/180; % angle converted to radians
Qxx = -(u .* b / (2.*pi .* (1-nu))) .* ((y .* (3.*x.^2 +y.^2)) ./ ((x.^2 + y.^2).^2)); % Qxx stress
Qyy = (u .* b / (2.*pi .* (1-nu))) .* ((y .* (x.^2 -y.^2)) ./ ((x.^2 + y.^2).^2)); % Qyy stress
Qxy = (u .* b / (2.*pi .* (1-nu))) .* ((x .* (x.^2 -y.^2)) ./ ((x.^2 + y.^2).^2)); % Qxy stressx2 = x.*cos(theta) + y.*sin(theta); % rotates dislocation
y2 = y.*cos(theta) - x.*sin(theta); % rotates dislocationQxxPRIME = -(u .* b / (2.*pi .* (1-nu))) .* ((y2 .* (3.*x2.^2 +y2.^2)) ./ ((x2.^2 + y2.^2).^2)); % Qyy stress PRIME
QyyPRIME = (u .* b / (2.*pi .* (1-nu))) .* ((y2 .* (x2.^2 -y2.^2)) ./ ((x2.^2 + y2.^2).^2)); % Qyy stress PRIME
QxyPRIME = (u .* b / (2.*pi .* (1-nu))) .* ((x2 .* (x2.^2 -y2.^2)) ./ ((x2.^2 + y2.^2).^2)); % Qxy stress PRIME
contourf(x(1,:),y(:,1)',QxxPRIME)
axis equal