止于至善

0%

Schmid因子

发表于 更新于 分类于

在晶体学中，由Schmid定律引出了Schmid因子的概念，该因子可用于描述滑移系中应力张量与分解剪应力之间的关系。

定义

晶体滑移系由两个矢量来定量描述，即滑移方向\(~\boldsymbol{s}~\)与滑移面法线方向\(~\boldsymbol{m}~\)。对于单轴受力情况，记应力矢量为\(~\boldsymbol{p}~\)，则分解剪应力为 \[ \tag{1a} \tau = \|\boldsymbol{p}\| \cos \lang \boldsymbol{p}~,~\boldsymbol{s} \rang \cos \lang \boldsymbol{p}~,~\boldsymbol{m} \rang = \|\boldsymbol{p}\| \cos \theta \cos \varphi \]

将其中的\(~\cos \theta \cos \varphi~\)称为Schmid因子，并记作 \[ \tag{1b} f = \cos \theta \cos \varphi \]

对于多轴受力情况，采用Cauchy应力张量\(~\boldsymbol{\sigma}~\)来描述其应力状态。根据Cauchy公式，可得作用在滑移面（法线方向为\(~\boldsymbol{m}\)）上应力矢量 \[ \tag{2a} \boldsymbol{t} = \boldsymbol{\sigma} \cdot \boldsymbol{m} \]

接着将应力矢量\(~\boldsymbol{t}~\)投影到滑移方向\(~\boldsymbol{s}~\)上，可得分解剪应力 \[ \tag{2b} \tau = \boldsymbol{t} \cdot \boldsymbol{s} = (\boldsymbol{\sigma} \cdot \boldsymbol{m}) \cdot \boldsymbol{s} = \boldsymbol{s} \cdot \boldsymbol{\sigma} \cdot \boldsymbol{m} = \boldsymbol{\sigma} : (\boldsymbol{s} \otimes \boldsymbol{m}) \]

实际上对于多轴应力状态，并没有一个直观的量来描述Schmid因子

算例：面心立方晶体

面心立方晶体通常有12个滑移系，在晶体局部坐标系下的具体方向为：

滑移系编号 滑移面法线方向 m 滑移方向 s
1 (1,1,1) [0,-1,1]
2 [1,0,-1]
3 [-1,1,0]
4 (-1,-1,1) [0,1,1]
5 [-1,0,-1]
6 [1,-1,0]
7 (-1,1,1) [0,-1,1]
8 [-1,0,-1]
9 [1,1,0]
10 (1,-1,1) [0,1,1]
11 [1,0,-1]
12 [-1,-1,0]

需要注意，一般\(~\boldsymbol{s}\)、\(\boldsymbol{m}~\)均取单位矢量，在实际使用时需要将表格中矢量单位化。

式\((1)\)和式\((2)\)提供了两种计算单轴应力下Schmid因子的思路。

Q：为什么采用式\((2)\)的多轴应力公式来计算单轴应力状态下的Schmid因子，似乎没有必要？

A：1、深入理解应力张量；2、式\((1)\)和式\((2)\)相互验证。

在计算面心立方晶体12个滑移系的Schmid因子时，单轴应力一般以全局坐标系下的应力矢量形式（记为\(~\{\sigma\}_g\)）给出，因此还需考虑到晶体局部坐标系与全部坐标系之间的关系，通常采用变换矩阵\(~\boldsymbol{Q}~\)来描述，详见《全局坐标系与局部坐标系下弹性矩阵的变换》。

采用式(1)计算

将全局坐标系下的应力矢量\(~\{\sigma\}_g~\)转换到晶体局部坐标系下，即 \[ \tag{3a} \{\sigma\}_l = \boldsymbol{Q}^\mathrm{T} \cdot \{\sigma\}_g \]

将晶体局部坐标系下的滑移方向\(~\boldsymbol{s}_l~\)、滑移面法线方向\(~\boldsymbol{m}_l~\)转换到全局坐标系下，即 \[ \tag{3b} \boldsymbol{s}_g = \boldsymbol{Q} \cdot \boldsymbol{s}_l \] \[ \tag{3c} \boldsymbol{m}_g = \boldsymbol{Q} \cdot \boldsymbol{m}_l \]

于是可以在晶体局部坐标系与全局坐标系下分别计算第\(~\alpha~\)个滑移系的Schmid因子，即 \[ \tag{4a} f^\alpha = \big(\boldsymbol{s}^\alpha_l \cdot \{\sigma\}_l\big)\big(\boldsymbol{m}^\alpha_l \cdot \{\sigma\}_l\big) \] \[ \tag{4b} f^\alpha = \big(\boldsymbol{s}^\alpha_g \cdot \{\sigma\}_g\big)\big(\boldsymbol{m}^\alpha_g \cdot \{\sigma\}_g\big) \]

不难证明，式\((4\mathrm{a})\)与式\((4\mathrm{b})\)是等价的

采用式(2)计算

首先需要将应力矢量\(~\{\sigma\}_g~\)转化为应力张量\(~\boldsymbol{\sigma}_g~\)。考虑到应力矢量\(~\{\sigma\}_g~\)对应一个单轴应力状态，为了方便计算Schmid因子，可令单轴应力的大小为1。那么，应力张量\(~\boldsymbol{\sigma}_g~\)的三个特征值分别为1、0、0。根据实对称矩阵的相似对角化原理，存在正交矩阵\(~\boldsymbol{A}~\)，使得应力张量\(~\boldsymbol{\sigma}_g~\)满足： \[ \tag{5} \boldsymbol{\sigma}_g = \boldsymbol{A} \cdot \mathrm{diag}([0,0,1]) \cdot \boldsymbol{A}^\mathrm{T} \]

式中，\(\mathrm{diag}([0,0,1])~\)为对角阵，对角阵中元素1的位置会影响正交矩阵\(~\boldsymbol{A}~\)的取值，本文将元素1置于对角阵的第三行第三列。

实际上单轴应力状态对应一个特殊的主应力微体，在该主应力微体上只有一个主应力不为零

接下来需要确定正交矩阵\(~\boldsymbol{A}~\)。应力张量\(~\boldsymbol{\sigma}_g~\)的特征值1对应的特征向量\(~\alpha_3~\)可取为 \[ \tag{6a} \alpha_3 = \frac{\{\sigma\}_g}{\|\{\sigma\}_g\|} \]

特征向量\(~\alpha_2~\)仅需满足垂直于\(~\alpha_3~\)，可取为 \[ \tag{6b} \alpha_2 = \big[-\alpha_3(2)~,~\alpha_3(1)~,~0\big]^\mathrm{T} \]

特征向量\(~\alpha_1~\)则可取为 \[ \tag{6c} \alpha_1 = \alpha_2 \times \alpha_3 \]

三个特征向量均为列向量，且构成正交右手系

\(\alpha_2~\)的取法并不唯一，例如还可取\(~\alpha_2 = \big[0~,~\alpha_3(3)~,~-\alpha_3(2)\big]^\mathrm{T}\)

于是，可得正交矩阵\(~\boldsymbol{A}~\)，即 \[ \tag{6d} \boldsymbol{A} = \big[\alpha_1~,~\alpha_2~,~\alpha_3\big]_{3\times3} \]

由式\((5)\)、\((6)\)计算得到全局坐标系下的应力张量\(~\boldsymbol{\sigma}_g~\)后，还可得到晶体局部坐标系下的应力张量\(~\boldsymbol{\sigma}_l~\)，即 \[ \tag{7} \boldsymbol{\sigma}_l = \boldsymbol{Q}^\mathrm{T} \cdot \boldsymbol{\sigma}_g \cdot \boldsymbol{Q} \]

变换矩阵\(~\boldsymbol{Q}~\)也为正交矩阵

于是可以在晶体局部坐标系与全局坐标系下分别计算第\(~\alpha~\)个滑移系的Schmid因子，即 \[ \tag{8a} f^\alpha = \boldsymbol{s}^\alpha_l \cdot \boldsymbol{\sigma}_l \cdot \boldsymbol{m}^\alpha_l \] \[ \tag{8b} f^\alpha = \boldsymbol{s}^\alpha_g \cdot \boldsymbol{\sigma}_g \cdot \boldsymbol{m}^\alpha_g \]

不难证明，式\((8\mathrm{a})\)与式\((8\mathrm{b})\)也是等价的

MATLAB程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 面心立方晶体的Schmid因子
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT
% sigma_g :  - 全局下载荷方向，列向量 3x1, 即应力矢量
% Q       :  - 变换矩阵, 正交矩阵 3x3, 即Ag=Q*Al, 且Ag与Al中向量按列存储
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUT
% y       :  - 12个Schmid因子，12x4
%            - 第一列: 全局坐标系下由应力矢量计算
%            - 第二列: 局部坐标系下由应力矢量计算
%            - 第三列: 全局坐标系下由应力张量计算
%            - 第四列: 局部坐标系下由应力张量计算
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 调用示例
% 全局坐标系下的随机应力矢量;随机正交变换矩阵;12个Schmid因子
% sigma=rands(3,1);q=orth(rands(3,3));y=Schmid(sigma,q)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y=Schmid(sigma_g,Q)

    [num,~]=size(sigma_g);
    if num==1
        sigma_g=sigma_g';% 转为列向量
    end
    sigma_g=sigma_g/norm(sigma_g);% 单位化

    tm=[1  1  1; 1  1  1; 1  1  1;-1 -1  1;-1 -1  1;-1 -1  1;
       -1  1  1;-1  1  1;-1  1  1; 1 -1  1; 1 -1  1; 1 -1  1];

    ts=[0 -1  1; 1  0 -1;-1  1  0; 0  1  1;-1  0 -1; 1 -1  0;
        0 -1  1;-1  0 -1; 1  1  0; 0  1  1; 1  0 -1;-1 -1  0];

    % 由主应力微体上的应力矢量反推应力张量
    lambda=diag([0,0,1]);% 由于将sigma_g指定为alpha3, 所以是[001]
    alpha3=sigma_g;

    L=(sigma_g==[0;0;0]);
    if L(1) && L(2) && ~L(3)
        A=diag([1,1,1]);
    elseif L(1) && ~L(2) && L(3)
        A=[0,1,0;0,0,1;1,0,0];
    elseif ~L(1) && L(2) && L(3)
        A=[0,0,1;1,0,0;0,1,0];
    else
        alpha2=[-sigma_g(2);sigma_g(1);0]; % 假设的, 只要与alpha3垂直即可
        %alpha2=[0;sigma_g(3);-sigma_g(2)];% 还可假设成这样
        alpha2=alpha2/norm(alpha2);
        alpha1=cross(alpha2,alpha3);
        A=[alpha1,alpha2,alpha3];
    end
    
    sigma_l=Q'*sigma_g;% 局部坐标系下的应力矢量 3x1
    sigma_mat_g=A*lambda*A';     % 全局下应力张量 3x3
    sigma_mat_l=Q'*sigma_mat_g*Q;% 局部下应力张量 3x3
    
    y=zeros(12,4);
    for i=1:12
        sl=ts(i,:)';sl=sl/norm(sl);% 局部 3x1
        ml=tm(i,:)';ml=ml/norm(ml);% 局部 3x1
        sg=Q*sl;mg=Q*ml;           % 全局 3x1

        y(i,1)=(sg'*sigma_g)*(mg'*sigma_g); % 全局坐标系下由应力矢量计算
        y(i,2)=(sl'*sigma_l)*(ml'*sigma_l); % 局部坐标系下由应力矢量计算
        y(i,3)=sg'*(sigma_mat_g*mg);        % 全局坐标系下由应力张量计算
        y(i,4)=sl'*(sigma_mat_l*ml);        % 局部坐标系下由应力张量计算
    end

end