Possible bugs with vector3d and Schmid factor related commands?

[AzdiarGazder]

Objective To calculate the Schmid factor using the angles between an externally imposed load and the slip plane normal and slip direction for an fcc partial slip system.

Example script

close all; clc; clear all; clear hidden;

% define the crystal system
CS = {'notIndexed',...
    crystalSymmetry('m-3m', [3.6 3.6 3.6], 'mineral', 'copper')};

% load an fcc dataset
mtexdata copper
ebsd.CSList = CS;
% consider only indexed fcc data
ebsd = ebsd('indexed');
ebsd = ebsd('copper');
% calculate the grains with a 15° threshold
[grains,ebsd.grainId] = calcGrains(ebsd,'angle',15*degree,'unitCell');

% define the plane & direction of a partial slip system
n = Miller(1,1,1,CS{2},'hkl');
b = Miller(1,1,-2,CS{2},'uvw');
% define the slip system
sS = slipSystem(b,n);
sS = symmetrise(sS)

% define the loading direction
r = diag([1,0,0]); % uniaxial tension
% r = diag([1,0,-1]); % rolling
% define the stress tensor
sigma = stressTensor(r).normalize;

% rotate slip system(s) into specimen coordinates
sSLocal = grains.meanOrientation * sS;
% compute Schmid factor
SF1 = sSLocal.SchmidFactor(sigma);
% sort SF in descending order along the rows
[SF1,sSNumber] = maxk(SF1,size(sS,1),2);

% calculate the angle between an externally imposed load and the plane normal
planeNormal = grains.meanOrientation.*sS(sSNumber).n;
theta = dot(normalize(planeNormal(:,1)),normalize(vector3d(r)),'noSymmetry');

% calculate the angle between an externally imposed load and the slip direction
slipDirN = grains.meanOrientation.*sS(sSNumber).b;
lambda = dot(normalize(slipDirN(:,1)),normalize(vector3d(r)),'noSymmetry');

% re-calculate the Schmid factor
SF2 = theta.*lambda;

Results and issues

(a) When the script is run with the following command:

r = diag([1,0,0]); % uniaxial tension

(a1) The command:

normalize(vector3d(r))

returns

ans = vector3d (show methods, plot)
 size: 1 x 3
    x   y   z
    1   0   0
  NaN NaN NaN
  NaN NaN NaN

Are the NaNs an error?

(a2) There are 3 columns of data returned for the variables theta, lambda and SF2. In all variables, column 1 contains numbers and columns 2 and 3 contain NaNs. For columns 2 and 3, I understand this is because rows 2 and 3 of the vector in point (a1) above contains NaNs. In any case, given that we are calculating one angle for each set of two 3d vectors, my understanding is that there should only be 1 column of data returned. If not, how do you compute one angle?

(a3) SF1 is calculated using MTex's Schmid factor script. SF2 is calculated as the product of theta and lambda. The SF2 values in column 1 matches the SF1 values in column 1. So for the case of uniaxial tension, the Schmid factors calculated using the 2 methods (SF1 & SF2) are equivalent.

(b) When the script is run with the following command:

r = diag([1,0,-1]); % rolling

(b1) The command:

normalize(vector3d(r))

returns

ans = vector3d (show methods, plot)
 size: 1 x 3
    x   y   z
    1   0   0
  NaN NaN NaN
    0   0  -1

Are the NaNs an error?

(b2) There are 3 columns of data returned for the variables theta, lambda and SF2. In all variables, columns 1 and 3 contain numbers and column 2 contains NaNs. For column 2, I understand this is because row 2 of the vector in point (b1) above contains NaNs. In any case, given that we are calculating one angle for each set of two 3d vectors, my understanding is that there should only be 1 column of data returned. If not, how do you compute one angle?

(b3) SF1 is calculated using MTex's Schmid factor script. SF2 is calculated as the product of theta and lambda. The SF2 values in columns 1 or 3 do not match the SF1 values in column 1. So for the case of rolling, the Schmid factors calculated using the 2 methods (SF1 & SF2) are not equivalent. The essential point of difference between (a3) and (b3) is the use of different stress tensors. But the end result should be equivalent - and I do not understand why this is not the case.

In summary, issues needing resolution are (a1), (a2), (b1), (b2) and (b3). I am unclear if they are Mtex issues or coding oversights/errors on my part. In any case, I would be very grateful for some help.

Thank you in advance! :)

MTEX version 5.8.1

[Tijmenvermeij]

In your code, r is a 3x3 matrix, representing a stress tensor, which you later define as 'sigma'. If you apply the vector3d command on the 3x3 matrix, you get 3 vectors out of it (because Mtex assumes you are giving a list of 3 vectors). If you the normalize the 0,0,0 vectors, you get NaNs.

My guess is, if you replace r=diag([1 0 0]) with r=vector3d(1,0,0) you will fix it, altough then you need to use sigma=stressTensor.uniaxial(r); You can also just define another vector representing your loading direction.

[AzdiarGazder]

Hi Tijmen, Thank you so much the comment. However, your suggeston is limited to the uniaxial tension/compression case only. You run into the issues I outlined in Part B of the Results and Issues section when you generalise the Schmid factor implementation for more complex stress states. Warm regards, Azdi

In your code, r is a 3x3 matrix, representing a stress tensor, which you later define as 'sigma'. If you apply the vector3d command on the 3x3 matrix, you get 3 vectors out of it (because Mtex assumes you are giving a list of 3 vectors). If you the normalize the 0,0,0 vectors, you get NaNs.

My guess is, if you replace r=diag([1 0 0]) with r=vector3d(1,0,0) you will fix it, altough then you need to use sigma=stressTensor.uniaxial(r); You can also just define another vector representing your loading direction.

[AzdiarGazder]

Hi everyone,

We've made a bit of progress on the above issue. As defined in Section 3.2.1 of the reference paper: https://www.sciencedirect.com/science/article/pii/S1359645411008792

The SF2 for rolling is calculated as:

SFx = theta(:,1).*lambda(:,1);
SFz = theta(:,3).*lambda(:,3);
SF2 = (SFx-SFz)/2;

In this case, the "2" in the denominator is described as a normalisation factor. As stated in the reference paper, this SF2 comprises "effective" Schmid factor values that are consistent with the above SF1(:,1). However, this only solves one issue.

Two larger questions remain:

1. What is the generalised formulation of the "effective" Schmid factor equation (i.e. - the line where SF2 is calculated) for even more complex (say, 3D) stress states? Is it something along the lines of the sentence in the reference paper that states "...the minus sign in the equation takes care of the fact that a favored (slip/)twinning mode will give extension along (x or)RD but contraction along (z or)ND and hence a large effective Schmid factor."?

2. How do you compute the "effective" angles between an externally imposed load and the plane normal or the slip direction in the case of rolling or even more complex (say, 3D) stress states? Does it make sense to think of "effective" angles or should we only compute the individual angles along the x (or RD), y (or TD) and z (or ND) directions and leave it at that? If this is the case for angles, then does an "effective" Schmid factor like SF2 make sense? Should SFs then also be stated along the 3 principal sample directions only?

Am I tying myself up in knots here? :)

Warm regards, Azdi

[kilir]

Hi, I'm not sure what you mean by computing individual Schmid factors in a 3D stress state but in general, define your stress tensor as stressTensor(m) where m is the appropriate 3x3 matrix.

Wrt 2: If I understand you'd like to compute the angles between principal stress directions and e.g. your slip plane? You get the former as the eigenvectors of your stress tensor, the latter by transforming those crystal directions to specimen coordinates (ior simply the other way around).

Cheers, Rüdiger

[AzdiarGazder]

Hi Rüdiger,

Thank you for replying.

1. I've updated my comment from yesterday to make my first question clearer. Please do have a look at Section 3.2.1 of the reference paper as well.

2. I've updated my comment from yesterday to make my second question clearer. Yes, I'd like to compute the angles between the principal stress directions and the slip plane and slip direction. This is what I mean by the "effective" angle between an externally imposed load and the plane normal or the slip direction Are the lines given below what you are suggesting?

t = dot(planeNormal(1,1).normalize.xyz,eig(sigma))
l = dot(slipDirN(1,1).normalize.xyz,eig(sigma))
SF3 = t*l;

If so, then SF3 does not match SF1(1,1).

Warm regards, Azdi

Hi, I'm not sure what you mean by computing individual Schmid factors in a 3D stress state but in general, define your stress tensor as stressTensor(m) where m is the appropriate 3x3 matrix.

Wrt 2: If I understand you'd like to compute the angles between principal stress directions and e.g. your slip plane? You get the former as the eigenvectors of your stress tensor, the latter by transforming those crystal directions to specimen coordinates (ior simply the other way around).

Cheers, Rüdiger

[Tijmenvermeij]

I think taking the angles between principle stress directions and slip plane/direction might be a bit too simplistic to calculate "partial" (?) Schmid factors. The Schmid factor should just be the ratio between the resolved shear stress on the slip system and some global stress value. The global stress value is just the uniaxial stress (normalized) under uniaxial loading. Under complex loading, you need to project the stress on the slip system, and then get the resolved stress after normalizing the stress tensor (as done in MTex), or do something that ensures the max SF is e.g. 0.5.

I can imagine there are some situations where you can calculate some "partial" Schmid factors based on the principle stresses, but I wouldn't just blindly do this.

Regarding the sign of the SF, this is not important for slip, as it can generally move in two directions (i.e. antipodal symmetry for the slip direction). For Twinning it can be important.

[AzdiarGazder]

Hi Tijmen,

Thank you so much for the reply. Could you please provide code for the solution you're suggesting? Ideally, I'd like to post a working generalised example of Schmid factor -based calculations for complex stress states for everyone to refer to in future.

Warm regards, Azdi

I think taking the angles between principle stress directions and slip plane/direction might be a bit too simplistic to calculate "partial" (?) Schmid factors. The Schmid factor should just be the ratio between the resolved shear stress on the slip system and some global stress value. The global stress value is just the uniaxial stress (normalized) under uniaxial loading. Under complex loading, you need to project the stress on the slip system, and then get the resolved stress after normalizing the stress tensor (as done in MTex), or do something that ensures the max SF is e.g. 0.5.

I can imagine there are some situations where you can calculate some "partial" Schmid factors based on the principle stresses, but I wouldn't just blindly do this.

Regarding the sign of the SF, this is not important for slip, as it can generally move in two directions (i.e. antipodal symmetry for the slip direction). For Twinning it can be important.

[Tijmenvermeij]

I'm pretty sure MTex already does what I suggested, but not totally sure. You can check the slipSystem/schmidFactor code in Mtex to see what it does when you input a stress tensor.