Calculating Schmid Factors For h.c.p

[Schmid h.c.p.]

I am attempting to use MTEX to calculate the Schmid factors for slip and twinning systems in h.c.p materials. The results are not correct but I am unsure what is wrong.  I am new to MTEX but have a small amount of Matlab experience.  

I have tried two approaches, one using the MTEX function calcShearStress and one attempting to calculate the angles between the shear direction and plane normal with angle_outer. Both approaches seem to give incorrect answers so I believe that it is an error with my use of the symmetry.  

I am using MTEX 4.0.9 and I have also tried with MTEX 3.5.0

I want to be able to use Euler angles obtained by Channel 5 hkl EBSD. I would be very grateful if someone can spot what I am doing wrong. The following code is what I have used so far:

%Lattice parameters

a=2.95; b=2.95; c=4.686;


alpha=90*degree; beta=90*degree; gamma=120*degree;

%Grain orientation Euler angles

Phi1=0*degree

Theta=0*degree

Phi2=0*degree

% crystal symmetry

CS = crystalSymmetry('6/mmm', [a,b,c], [alpha, beta, gamma], 'X||a*', 'Y||b', 'Z||c', 'mineral', 'Ti-Hex', 'color', 'light blue');

% specimen symmetry

SS = specimenSymmetry('triclinic');


%Grain orientation


ori = orientation('Euler', Phi1, Theta, Phi2, 'bunge', CS, SS);


%Uniaxial Loading dire tions (orthogonal directions)

xcs = Miller (1,0,-1,0, CS, 'uvw')

ycs = Miller (-1,2,-1,0, CS, 'uvw')

zcs = Miller (0,0,0,1, CS, 'uvw')


%Plot loading directions on a stereographic plot

plot(xcs)

hold all


plot(ycs)

plot(zcs)

hold off

% orientation matrix

g = Euler2mat (Phi1, Theta, Phi2, 'bunge', CS,SS);


gT = transpose(g)

% Obtain loading directions in terms of the crystal orientation


% [1 0 -1 0], [-1 2 -1 0] and [0 0 0 1] in three index array to multiply

% with g matrix

x = [1, 0.5, 0]

y = [0, 1, 0]

z = [0, 0, 1]


xcs = gT * x';

ycs = gT * y';

zcs = gT * z';


xcs = Miller(xcs(1),xcs(2),xcs(3),CS,'uvw')

ycs = Miller(ycs(1),ycs(2),ycs(3),CS,'uvw')

zcs = Miller(zcs(1),zcs(2),zcs(3),CS,'uvw')

xcs = normalize(xcs)

ycs = normalize(ycs)

zcs = normalize(zcs)



% {0001} Slip

% Shear direction

mx = Miller(1,1,-2,0,CS,'uvw');

% Plane

k = Miller(0,0,0,1,CS,'hkl')

k = [0, 0, 1]

 

%calculate plane normal

 

T = [b^2*c^2*sin(alpha)^2, a*b*c^2*(cos(alpha)*cos(beta)-cos(gamma)), a*b^2*c*(cos(alpha)*cos(gamma)-cos(beta));

    a*b*c^2*(cos(alpha)*cos(beta)-cos(gamma)), a^2*c^2*sin(beta)^2, a^2*b*c*(cos(beta)*cos(gamma)-cos(alpha));

    a*b^2*c*(cos(alpha)*cos(gamma)-cos(beta)), a^2*b*c*(cos(beta)*cos(gamma)-cos(alpha)), a^2*b^2*sin(gamma)^2];

 
n = T * k'

n = Miller(n(1),n(2),n(3),CS,'uvw');

 

nx = normalize(n)




rx = xcs


ry = ycs

rz = zcs



% Find symmetrically equivalent slip and plane normal directions

  [m,l] = symmetrise(mx),CS,SS

  [n,l] = symmetrise(nx),CS,SS

 
plot(m,'antipodal')

hold all

plot(n)

hold off



% Find planes and directions that are normal

[r,c] = find(isnull(dot_outer(vector3d(m),vector3d(n))));

 

m = m(r);

n = n(c);

 

plot(m,'antipodal')

hold all

plot(n)

plot(rx)

hold off

 

tau = zeros(numel(m),numel(rx));


%Calculate individual Schmid factors by calculating the angles

 

for i = 1:numel(m)

 

anglem = degree - angle_outer(m(i),rx)

 

ang = anglem / degree

anglen = degree - angle_outer(n(i),rx)

 

ang = anglen / degree

 

R = cos(anglem)*cos(anglen)

 

  tau(i,:) = R;

 

end

 

if numel(m)>1

  [tauMaxx,ind] = max(abs(tau));

 

  m = m(ind);

  n = n(ind);

else

  tauMaxx = tau;

  ind = ones(size(tauMaxx));

end

 

tauMaxx

 

%Clear cache and work out for y direction ------------------

 

clear m

clear n

clear i

clear tau

clear r

clear c

 

  [m,l] = symmetrise(mx),CS,SS

  [n,l] = symmetrise(nx),CS,SS

 

plot(m,'antipodal')

hold all

plot(n)

hold off

 

[r,c] = find(isnull(dot_outer(vector3d(m),vector3d(n))));

 

m = m(r);

n = n(c);

 

plot(m,'antipodal')

hold all

plot(n)

plot(ry)

hold off

 

tau = zeros(numel(m),numel(ry));

 

%Calculate individual Schmid factors by calculating the angles

 

for i = 1:numel(m)

 

anglem = angle_outer(m(i),ry)

%anglem = min(anglem(:))

ang = anglem / degree

anglen = angle_outer(n(i),ry)

%anglen = min(anglen(:))

ang = anglen / degree

 

R = cos(anglem)*cos(anglen)

 

  tau(i,:) = R;

 

end

 

if numel(m)>1

  [tauMaxy,ind] = max(abs(tau));

 

  m = m(ind);

  n = n(ind);

else

  tauMaxy = tau;

  ind = ones(size(tauMaxy));

end

 

tauMaxy

 

%Clear cache and work out for z direction ------------------

 

clear m

clear n

clear i

clear tau

clear r

clear c

 

  [m,l] = symmetrise(mx),CS,SS

  [n,l] = symmetrise(nx),CS,SS

 

plot(m,'antipodal')

hold all

plot(n)

hold off

 

[r,c] = find(isnull(dot_outer(vector3d(m),vector3d(n))));

 

m = m(r);

n = n(c);

 

plot(m,'antipodal')

hold all

plot(n)

plot(rz)

hold off

 

tau = zeros(numel(m),numel(rz));

 

%Calculate individual Schmid factors by calculating the angles

 

for i = 1:numel(m)

 

anglem = angle_outer(m(i),rz)

%anglem = min(anglem(:))

ang = anglem / degree

anglen = angle_outer(n(i),rz)

%anglen = min(anglen(:))

ang = anglen / degree

 

R = cos(anglem)*cos(anglen)

 

  tau(i,:) = R;

 

end

 

if numel(m)>1

  [tauMaxz,ind] = max(abs(tau));

 

  m = m(ind);

  n = n(ind);

else

  tauMaxz = tau;

  ind = ones(size(tauMaxz));

end

tauMaxz




% Display Schmid factors

 

h=msgbox({'x,y,z:', num2str(tauMaxx), num2str(tauMaxy), num2str(tauMaxz)}, 'Schmid Factors')

[Ralf Hielscher]

Hi sl655,

it was not so simply to read you code. Please try to condense it a bit when publishing it here. I rewrote the first half of it and removed some mistakes.

%// crystal symmetry
CS = crystalSymmetry('6/mmm', [2.95,2.95,4.686], 'X||a*', 'Y||b', 'Z||c', 'mineral', 'Ti-Hex', 'color', 'light blue');

%// Grain orientation Euler angles
ori = orientation('Euler', 0*degree, 0*degree, 0*degree, CS);

%%
%// Obtain loading directions in terms of the crystal orientation
%// [1 0 -1 0], [-1 2 -1 0] and [0 0 0 1] in three index array to multiply
%// with g matrix

xcs = round(inv(ori) * xvector);
ycs = round(inv(ori) * yvector);
zcs = round(inv(ori) * zvector);

%%
%// {0001} Slip Shear direction

mx = Miller(1,1,-2,0,CS,'uvw');

%// Plane

k = Miller(0,0,0,1,CS,'hkl')

%// there is not so much a difference between k and nx
nx = k; nx.dispStyle = 'uvw';

%// Find symmetrically equivalent slip and plane normal directions
[m,l] = symmetrise(mx,CS)
[n,l] = symmetrise(nx,CS)

 
plot(m,'antipodal')
hold all
plot(n)
hold off

%%

rx = xcs, ry = ycs, rz = zcs

%// Find planes and directions that are normal

[r,c] = find(isnull(dot_outer(vector3d(m),vector3d(n))));
m = m(r); n = n(c);
 
plot(m,'antipodal')
hold all
plot(n)
plot(rx)
hold off

%%
%// Calculate individual Schmid factors by calculating the angles

%// you wrote here "degree - angle_outer(m(i),rx)" which makes no sense
anglem = angle_outer(m,rx);
anglen = angle_outer(n,rx);
   
tau = cos(anglem).*cos(anglen);
 
if length(m)>1

  [tauMaxx,ind] = max(abs(tau));
  m = m(ind);
  n = n(ind);
else
  tauMaxx = tau;
  ind = ones(size(tauMaxx));
end

Maybe, the following lines of code do a better job and are more understandable

%// stress Tensor
loadingDirection = xvector;

sigma = EinsteinSum(tensor(loadingDirection),1,tensor(loadingDirection),2,'name','stress')

%// {0001} Slip Shear direction

k = Miller(0,0,0,1,CS,'hkl')

m = Miller(1,1,-2,0,CS,'uvw')

%// the maximum Schmid factor this slip system
[tauMax,mActive,nActive] = calcShearStress(sigma,m,k,'symmetrise')

Have you ever read through http://mtex-toolbox.github.io/files/doc/PlasticDeformation.html

I hope this helps,

Ralf.

[Schmid h.c.p.]

Thank you for having a look at it Ralf.

I apologise for the length of the code that I uploaded, I do not have much Matlab or MTEX experience.  But I thought that my problem was worth posting as I could not find an example of h.c.p. Schmid factor calculation.

"round(inv(ori) * xvector)" is a much more efficient way of representing the loading directions with respect to the crystal coordinate system.

In the example that I gave with k plane (0001) the plane normal nx is [0001] but that is a unique case for the (0001) plane in h.c.p. materials.  All other planes will need to have the plane normal calculated.  Is there a more efficient way of computing the plane normal for the crystal coordinate system than with the T matrix in my previous code?  I could write the T matrix into a function if there is not one that already exists.

I thought that a for loop for the angle calculation would be needed to calculate each angle, but anglem = angle_outer(m,rx) is more concise.  The mistake I left in there was from when I was trying to output the angle in degrees for debugging.

And the link that you provided is useful.  I had not found that.

The method of calculating the angles is working fine but the MTEX built in function (calcShearStress) is not correct, suggesting that there is a problem with the function and crystal symmetry other than cubic.  Taking h.c.p. Cadnium as an example, with a = 2.98, b= 2.98 and c = 5.62 the Schmid factor for slip on the (001) plane in the [010] direction with loading direction [021] should be 0.499 (http://www.doitpoms.ac.uk/tlplib/slip/printall.php).  By working through the angles with the code I get a Schmid factor of 0.49914 but the calcShearStress function returns a Schmid factor of 8.3738 which is complete rubbish.

%// crystal symmetry

CS = crystalSymmetry('6/mmm', [2.98,2.98,5.62], 'X||a*', 'Y||b', 'Z||c', 'mineral', 'Cadnium', 'color', 'light blue');

 

%// Grain orientation Euler angles

ori = orientation('Euler', 0*degree, 0*degree, 0*degree, CS);

 

%%

%// Loading direction

 

xcs = Miller(0,1,0.5,CS,'uvw');

 

%%

%// {001} Slip Shear direction

 

mx = Miller(0,1,0,CS,'uvw');

 

%// Plane normal

 

nx = Miller(0,0,1,CS,'uvw');

 

%%

 

%// Find symmetrically equivalent slip and plane normal directions

[m,l] = symmetrise(mx,CS)

[n,l] = symmetrise(nx,CS)

 

 

%%

 

rx = xcs

 

%// Find planes and directions that are normal

 

[r,c] = find(isnull(dot_outer(vector3d(m),vector3d(n))));

m = m(r); n = n(c);

 

%%

%// Calculate individual Schmid factors by calculating the angles

 

anglem = angle_outer(m,rx);

anglen = angle_outer(n,rx);

   

tau = cos(anglem).*cos(anglen);

 

if length(m)>1

 

  [tauMaxx,ind] = max(abs(tau));

  m = m(ind);

  n = n(ind);

else

  tauMaxx = tau;

  ind = ones(size(tauMaxx));

end

 

% Display Schmid factors

h = msgbox({'x,y,z:', num2str(tauMaxx)}, 'Schmid Factors')

%%

%// stress Tensor (Does not work)

loadingDirection = xcs;

 

sigma = EinsteinSum(tensor(loadingDirection),1,tensor(loadingDirection),2,'name','stress')

 

%// Slip Shear direction and plane normal

m=mx

k=nx

 

%// the maximum Schmid factor this slip system

[tauMax,mActive,nActive] = calcShearStress(sigma,m,k,'symmetrise')

 

h = msgbox({'x,y,z:', num2str(tauMax)}, 'Schmid Factors')