晶體塑性有限元仿真入門(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°]

     

單個歐拉角數據測試：