晶体塑性有限元仿真入门(5)—欧拉角与晶体取向 

备注：网页排版有乱码，建议下载附件pdf查看

晶体取向是材料学科中的重要分支，当晶粒发生择优取向时，则导致材料性能（力学，物理和化学性能）的各向异性。各向异性会造成材料实际应用中的各种问题，如铝合金典型的制耳现象，再如取向硅钢中存在Goss织构时，有利于其磁学性能。在基础研究领域，织构的形成与演变是基本的科学问题。在工业应用领域，通过织构的设计和控制可以提高材料的性能。随着近年来EBSD和XRD等表征技术的发展，各种SCI期刊的发文都已离不开对材料晶体学取向的分析。这篇文章介绍晶体塑性有限元仿真过程中的欧拉角与晶体取向。  

   

图1 塑性变形过程导致的材料各向异性

全文包括以下几个部分：

1) 材料晶体结构

2) EBSD工作原理

3) 晶体取向分析

4) 晶体塑性材料模型

5) 织构演变结果

6) 参考资料

7) 附录

材料晶体结构 

在晶体学中，晶体结构是对晶体材料中原子、离子或分子有序排列的描述。有序结构由组成粒子的内在性质产生，形成沿物质三维空间的主要方向重复的对称模式，如图2所示。

   

图2 高分辨率透射电子显微镜图片的铁晶体，完美单晶的二维示意图  

     

构成这种重复图案的材料中最小的一组粒子是结构的晶胞。晶胞完全反映了整个晶体的对称性和结构，这是通过晶胞沿其主轴重复平移而建立的。平移向量定义布拉维点阵的节点，不同的晶体内部原子排列称为具有不同的晶格结构。各种晶格结构可以归纳为七大晶系，各种晶系分别与十四种空间格（称为Bravais晶格）相对应，如图3所示。

图3 三种具有立方体对称性的Bravais晶格：simple/primitive cubic (sc), the body centered cubic (bcc) and the face centered cubic (fcc)lattice

     

实际上的金属材料大多是多晶体，多晶体是由许多形状、大小、取向各不相同的单晶体晶粒所组成的，如图4所示。每个晶粒里面的晶体取向相同(并不完全相同，只是差别很小)，通过电子背散射衍射分析技术(EBSD)分析可以定量获得这些晶体取向信息。

   

图4 EBSD分析获得的反极图分布图(颜色代表取向)，多晶体组成示意图

       

EBSD工作原理 

如图5所示，EBSD利用从样品表面反弹回来的高能电子衍射，得到一系列的菊池花样。根据菊池花样的特点得出晶面间距d和晶面之间的夹角θ，从数据库中查出可能的晶体结构和晶胞参数。再利用化学成分等信息采用排除法确定该晶粒的晶体结构，并得出晶粒与膜面法向的取向关系。通过以上基本操作，可以得到样品表面测试点的Phase和Orientation实验数据。

  

图5 EBSD分析的工作原理(Phase discriminated, Orientation determined)

     

通过步进间距将样品表面划分为m*n个采样点，并依次获得这m*n个采样点菊池花样和匹配的Phase、Orientation实验数据，最后对m*n个采样点的这些数据进行整理、汇总、计算等，可以进行晶粒尺寸分析、织构分析、相分析、应变分析、再结晶分析、晶界分析、断裂分析等。如图6所示，EBSD获得多晶体试样表面的取向信息过程，类似于相机的成像过程，首先将画面分解成为m*n个像素点，然后依次获得各个像素点的信息，获得全部信息后进行整理、汇总、计算。

图6 EBSD分析的工作原理(多晶体分析)

       

晶体取向分析 

欧拉角 

通过EBSD实验，我们可以获得多晶体试样表面m*n个采样点的取向信息(Euler1、Euler2、Euler3)。与彩色图像每个像素点存在RGB三组信息类似，取向信息的每个采样点存在三组欧拉角角度信息，这三组角度信息代表了晶格在空间中的唯一取向信息。如图7所示，空间中任意取向的晶格，通过将全局坐标系依次采用三组欧拉角进行旋转后，都可以与晶体坐标系重合。第一次旋转是围绕Z轴旋转的，第二次旋转是围绕新的X轴，第三次旋转角度围绕新的Z轴。与三维空间中点的位置信息包含X、Y、Z三个自由度一样，三维空间中晶格的取向信息也包含三个自由度，记三个欧拉角为ϕ1、Φ、ϕ2。

图7 三次旋转与任意取向晶格重合的示意图

       

取向矩阵 

除了使用欧拉角表示空间中晶格的取向，还可以使用取向矩阵(Orientation matrix, Rotation matrix) g(ϕ1,Φ,ϕ2)进行描述。取向矩阵的表达式是是通过组合三个旋转得到的，满足：

g(ϕ1,Φ,ϕ2) = g(ϕ1)*g(Φ)*g(ϕ2)

那么：

由于取向矩阵的性质，满足：

假设在全局坐标系S下，局部坐标系s₁的取向矩阵为g₁，局部坐标系s₂的取向矩阵为g₂，那么，在局部坐标系s₁下的坐标系s₂的取向矩阵可以看作：在局部坐标系s₁下首先反方向旋转一组得到全局坐标系S，取向矩阵为g₂为，然后在得到全局坐标系S情况下继续旋转取向矩阵g₂，得到局部坐标系s₂，即：

例如，全局坐标系S下，欧拉角为[90,35,45]的局部坐标系s₁的取向矩阵g₁为：  

全局坐标系S下，欧拉角为[142.8,32.0,214.4]的局部坐标系s₂的取向矩阵g₂为：

那么，在局部坐标系s₁下的坐标系s₂的取向矩阵等于：

可以计算出（这些取向描述方法如何相互计算将在后面讨论）取向矩阵对应的欧拉角为[136.20,66.68,83.91]，其数值并不等于两个坐标系的欧拉角简单加减。

同理，假设在全局坐标系S下，局部坐标系s₁的取向矩阵为g₁，在局部坐标系s₁下，局部坐标系s₂的取向矩阵为g₂’。那么，在全局坐标系S下，依次旋转g₁和g₂’，得到局部坐标系s₂的取向矩阵g₂为：

同理，假设在全局坐标系S下，局部坐标系s₂的取向矩阵为g₂，在局部坐标系s₁下，局部坐标系s₂的取向矩阵为g₂’。那么，在全局坐标系S下，依次旋转g₂和g₂’^{-1，局部坐标系s}₁的取向矩阵为g₁：

     

     

弥勒指数 

除了使用欧拉角、取向矩阵表示空间中晶格的取向，还可以使用弥勒指数进行描述，表示格式为(hkl)[uvw]。如(-112)[1-11]，表示晶胞的(-112)面平行于轧板的轧面，[1-11]方向平行于轧向。已知一个晶面(hkl)和它所属的晶向[uvw]，就很容易得到二者之间的关系：hu+kv+lw=0，通常把这个关系式称为晶带定律。  

以上三种方法是描述晶格取向最常见的方法，由于三维空间中晶格的取向信息包含三个自由度，因此以上三种方法并不独立，而是存在如下联系。

图8 欧拉角、取向矩阵与弥勒指数之间的联系

     

晶体塑性材料模型 

晶体塑性材料模型在ABAQUS中作为用户材料子程序(Huang’s UMAT)实现，取向信息通过inp文件中材料参数部分里的local system和global system实现赋予。

图9 晶体取向信息的赋予  

     

变形过程中通过SDV将取向数据进行输出，SDV37~39和SDV73~75代表第一组滑移系(1,1,1)[0,-1,1]对应的滑移面法向和滑移方向，取向信息表示为g1。由于第一组滑移系的取向信息是固定的G1，因此，在全局坐标系下，晶体取向信息可表示为：G = G1*g1^-1。例如，晶体初始取向信息为[90,35,45]，全部滑移系输出的晶体取向信息可表示为如下：

表1 初始取向信息为[90,35,45]12组滑移系输出的晶体取向信息

h

k

l

u

v

w

H

K

L

U

V

W

Slip

0.406

0.406

0.819

-0.579

-0.579

0.574

0

0

1

1

0

0

0

-0.338

0.000

0.941

0.815

0.500

0.292

1

1

1

0

-1

1

1

-0.338

0.000

0.941

-0.815

0.500

-0.292

1

1

1

1

0

-1

2

-0.338

0.000

0.941

0.000

-1.000

0.000

1

1

1

-1

1

0

3

0.331

-0.816

0.473

-0.004

0.500

0.866

-1

1

1

1

0

1

4

0.331

-0.816

0.473

-0.819

0.000

0.574

-1

1

1

1

1

0

5

0.331

-0.816

0.473

0.815

0.500

0.292

-1

1

1

0

-1

1

6

0.331

0.816

0.473

-0.004

-0.500

0.866

1

-1

1

0

1

1

7

0.331

0.816

0.473

-0.819

0.000

0.574

1

-1

1

1

1

0

8

0.331

0.816

0.473

-0.815

0.500

-0.292

1

-1

1

1

0

-1

9

-1.000

0.000

-0.005

-0.004

-0.500

0.866

1

1

-1

0

1

1

10

-1.000

0.000

-0.005

-0.004

0.500

0.866

1

1

-1

1

0

1

11

-1.000

0.000

-0.005

0.000

-1.000

0.000

1

1

-1

-1

1

0

12

     

     

  

图8 274个晶粒的初始织构(0~180°随机织构)

     

建立完模型后对第一增量步的晶体取向(初始取向)进行验证，如图9所示，说明有限元模型被正确的赋予了这些随机取向，并验证了取向计算程序的正确。

     

     

  

图9 建立模型后对第一步晶体取向的验证

图11 织构演变模拟常见的边界条件

     

织构演变结果 

完成Abaqus构建有限元模型所有关键步骤后，输出inp文件并提交Job，查看织构演变结果如下(由于计算资源的限制，仅计算了simple compression和plane strain compression)：

  

 

simple compression

     

  

 

plane strain compression

     

以多晶体中一号节点为例，在塑性变形过程中它的织构演变如下：

  

1号节点织构取向演变

     

参考资料 

Polycrystalline Plasticity and the Evolution of Crystallographic Texture in FCC Metals

Texture evolution and mechanical behaviour of irradiated face-centred cubic metals

A User-Material Subroutine Incorporating Single Crystal Plasticity in the ABAQUS Finite Element Program

     

附件 

取向旋转矩阵计算 

% 取向旋转矩阵计算，本程序适用于将欧拉角(角度制)转换为取向旋转矩阵(G)

function Rotation = Euler_to_Rotation(Euler_Input)

euler_1 = Euler_Input(1,1)/180*pi;

euler_2 = Euler_Input(1,2)/180*pi;

euler_3 = Euler_Input(1,3)/180*pi;

 

u0=cos(euler_1)*cos(euler_3)-sin(euler_1)*sin(euler_3)*cos(euler_2);

v0=-cos(euler_1)*sin(euler_3)-sin(euler_1)*cos(euler_3)*cos(euler_2);

w0=sin(euler_1)*sin(euler_2);

uvw0 = [u0;v0;w0];

r0=sin(euler_1)*cos(euler_3)+cos(euler_1)*sin(euler_3)*cos(euler_2);

s0=-sin(euler_1)*sin(euler_3)+cos(euler_1)*cos(euler_3)*cos(euler_2);  

t0=-cos(euler_1)*sin(euler_2);

rst0 = [r0;s0;t0];

h0=sin(euler_3)*sin(euler_2);

k0=cos(euler_3)*sin(euler_2);

l0=cos(euler_2);

hkl0 = [h0;k0;l0];

Rotation = [uvw0,rst0,hkl0];

     

欧拉角计算 

% 欧拉角计算，本程序适用于将取向旋转矩阵(G)转换为欧拉角(角度制0~180)

function Euler = Rotation_to_Euler(Rotation_Input)

phi1 = atan2(-Rotation_Input(3,1), Rotation_Input(3,2))*180/pi;

phi2 = atan2(Rotation_Input(1,3), Rotation_Input(2,3))*180/pi;

phi = atan2(sqrt(Rotation_Input(1,3)^2+Rotation_Input(2,3)^2), ...

   Rotation_Input(3,3))*180/pi;

Euler_origin = [phi1,phi,phi2];

Euler_origin(Euler_origin<0)=180+Euler_origin(Euler_origin<0);

 

euler_1 = [Euler_origin(1),Euler_origin(1)+180];

euler_2 = [Euler_origin(2),Euler_origin(2)+180];

euler_3 = [Euler_origin(3),Euler_origin(3)+180];

Num = 500;

Euler_all_flag = randi(2,[Num,3]);

Rotation_Input_error = zeros(Num,1);

Euler_all = zeros(Num,3);

for i = 1 : Num

   Euler_all(i,:) = [euler_1(Euler_all_flag(i,1)),euler_2(Euler_all_flag(i,2)),euler_3(Euler_all_flag(i,3))];

   Rotation_Input_error(i,:) = sum(sum(Euler_to_Rotation(Euler_all(i,:)) - Rotation_Input)).^2;

end

Euler = Euler_all(find(Rotation_Input_error == min(Rotation_Input_error),1),:);

     

SDV提取保存 

# SDV提取，本程序适用于将odb输出文件的SDV提取导出保存

     

import odbAccess

import numpy as np

     

myOdb = odbAccess.openOdb(r"D:\CPFEM_Abaqus_V2\Course3_3\Job-1201-001.odb",readOnly=False)

myStep = myOdb.steps["Step-1"]

myFrame = myStep.frames

i = len(myFrame)-1

myField = myFrame[i].fieldOutputs

AA = 3;

myValue = myField["SDV37"].values

temp_h = []

for value in myValue:  

   myData = value.data

   temp_h.append(myData)

     

myValue = myField["SDV38"].values

temp_k = []

for value in myValue:

   myData = value.data

   temp_k.append(myData)

     

myValue = myField["SDV39"].values

temp_l = []

for value in myValue:

   myData = value.data

   temp_l.append(myData)

     

myValue = myField["SDV73"].values

temp_u = []

for value in myValue:

   myData = value.data

   temp_u.append(myData)

     

myValue = myField["SDV74"].values

temp_v = []

for value in myValue:

   myData = value.data

   temp_v.append(myData)

     

   myValue = myField["SDV75"].values

temp_w = []

for value in myValue:

   myData = value.data

   temp_w.append(myData)

     

Data = temp_h+temp_k+temp_l+temp_u+temp_v+temp_w

     

np.savetxt(r"D:\CPFEM_Abaqus_V2\Course3_3\Out_end.txt",Data,fmt="%.5f")

     

SDV数据转换 

% SDV数据转换，本程序适用于将输出SDV数据转换为欧拉角(角度制0~180)

function Euler = SDV_to_Euler(SDV_Input_hkl,SDV_Input_uvw)

hkl=[1,1,1]; uvw=[0,-1,1];

b_bf = transpose( uvw/sqrt(sum(uvw.^2)) );

n_bf = transpose( hkl/sqrt(sum(hkl.^2)) );

t_bf = cross(n_bf,b_bf) / sqrt(sum(cross(n_bf,b_bf).^2));

G_1 = [b_bf,t_bf,n_bf];

 

SDV_Input_rst = cross(SDV_Input_hkl,SDV_Input_uvw) /...

   sqrt(sum(cross(SDV_Input_hkl,SDV_Input_uvw).^2));

g_1 = [transpose(SDV_Input_uvw),transpose(SDV_Input_rst),transpose(SDV_Input_hkl)];

 

g_ij_new = G_1 * inv(g_1);

Euler = Rotation_to_Euler(g_ij_new);

       

     

SDV数据批处理 

% SDV数据批处理，本程序适用于将输出SDV数据转换为欧拉角(角度制0~180)

clc

clear

close all

% Data = load('D:\CPFEM_Abaqus_V2\Course3_3\Out_end_rolling.txt');

Data = load('D:\CPFEM_Abaqus_V2\Course3_3\Out_end_compression.txt');

HKLandUVW = reshape(Data,length(Data)/6,6);

data = HKLandUVW(1:1:end,:);

Unique_data = data(1:20:end,:);

Euler = zeros(length(Unique_data(:,1)),3);

for i = 1:length(Unique_data(:,1))

   Euler(i,:) = SDV_to_Euler(Unique_data(i,1:3),Unique_data(i,4:6));

end

Unique_Euler = Euler;

plot(Unique_Euler(:,1),'r.');hold on

plot(Unique_Euler(:,2),'k.');hold on

plot(Unique_Euler(:,3),'b.');hold on

     

     

多晶体取向可视化 

% 取向可视化，本程序适用于将多晶体取向可视化

cs = crystalSymmetry('cubic');        

ss = specimenSymmetry('triclinic');

ori = orientation.byEuler(Euler(:,1)*degree,Euler(:,2)*degree,Euler(:,3)*degree,cs,ss);

plotPDF(ori,Miller({1,0,0},{1,1,0},{1,1,1},cs),'all','MarkerSize',1)

pause(5)

     

Index = transpose(1:length(Euler(:,1)));

Phase = ones(length(Euler(:,1)),1);

X = transpose(1:length(Euler(:,1)));

Y = transpose(1:length(Euler(:,1)));

MAD = rand(length(Euler(:,1)),1);

BC = ceil(100*rand(length(Euler(:,1)),1));

BS = zeros(length(Euler(:,1)),1);

Bands = ceil(10*rand(length(Euler(:,1)),1));

Error = zeros(length(Euler(:,1)),1);

Reli = ones(length(Euler(:,1)),1);

Data = [Index, Phase, X, Y, Euler, MAD, BC, BS,...

   Bands, Error, Reli];

 

for i = 1 : length(Euler(:,1))

   fid = fopen('Euler_to_EBSD_File.txt','a+');

   if i == 1

       fprintf(fid,'Index Phase  X  Y  phi1  Phi  phi2  MAD  BC  BS  Bands Error Reli  \n');

   end

   fprintf(fid,'%4i %5i %6.1f %4.1f %10.3f %6.2f %6.2f %6.4f %2i %5i %8i %5i %10.4f \n',Data(i,:));  

   fclose(fid);

end

 

cs = crystalSymmetry('m-3m', 'mineral', 'Cu') ;

fname = fullfile('D:\Download_Files\Download_Doc_From_IE_Browser\Spider_Teacher','Euler_to_EBSD_File.txt');

ebsd = loadEBSD(fname, 'CS' ,cs,...

'interface', 'generic', 'Bunge', 'ignorePhase' ,[0 2],...

'ColumnNames', { 'Phase' 'x' 'y' 'phi1' 'Phi' 'phi2'},...

'Columns', [2 3 4 5 6 7]);

     

odf = calcODF(ebsd('Cu').orientations, 'halfwidth' ,15*degree);

h = [Miller(0,0,1,cs),Miller(0,1,1,cs),Miller(1,1,1,cs)];

plotPDF(odf,h,'points','all');

     

     

 

       

网上参考数据测试：

http://muchong.com/bbs/viewthread.php?tid=12562139

     

随机欧拉角数据测试：

  

[-180°~180°]

[-90°~90°]=[0°~180°] 1000seeds

       

[-90°~90°]=[0°~180°] 8000seeds

     

[-90°~90°]=[0°~180°] 80000seeds

     

     

  

[-45°~45°]

     

[0°~90°]

     

[0°~90°] 80000seeds

     

[-10°~10°]

     

单个欧拉角数据测试：