材料弹性模量各向异性的三维图示方法

类别:    标签: 数理 matlab   阅读次数:   版权: (CC) BY-NC-SA

背景原理简介

各向异性广泛存在各种材料之中. 什么是各向异性呢? 简单来说, 就是晶体在不同方向有不同的性能. 我们很容易把各向异性和不均匀性混淆, 其实这是两个完全不同的概念. 各向异性就是说材料的性能与方向有关, 而不均匀性是指材料的性能与部位有关. 拿单晶来说, 内部任一点, 结构与性能都是相同的, 但不同方向却有不同的性能.

晶体是各向异性的, 不同方向性能不一样, 而且有着严格的对称性. 这里要说一下, 对称性是一个很了不起的性质, 在数学, 物理学中都有着广泛的应用, 而且在材料性质的分析中也常用到. 多晶材料存在择优取向, 也有一定的各向异性.

各种材料都有弹性, 大多数材料的弹性性质也具有各向异性. 例如, 在立方晶体中[111]方向通常比[100]方向更难压缩(stiff). 当我们对材料施加载荷, 材料会发生相应的形变, 在弹性范围内, 形变遵循胡克(Hook)定律, 即应力与应变是线性关系, 可以表示为 $\s = C \ve$, 其中, $\s$ 是应力, $\ve$ 是应变, $C$ 为杨氏模量(或称弹性模量), 也常用 $E$ 或 $Y$ 来表示.

材料不同方向上弹性模量不同, 我们怎么描述这种不同呢? 最好用数学方法, 建立数学框架, 准确直观地将弹性各项异性描述出来. 下面我们就进行这种数学的描述. 本人数学水平有限, 不能一步一步推导, 但我们可以简要理解一下推导过程. 弹性各项异性的推导就是利用张量和群论推广胡克定律. 我们先只考虑低阶弹性常数. 考虑二阶的 $C_{ijkl}$, 原来的胡克定律 $\s = C \ve$ 可推广为为矩阵形式

可以证明, 刚度矩阵 $\mathbf C$ 为对称阵, $C_{ij}=C_{ji}$. 因此, 独立张量元数目至多只有21个. 晶系的对称性越高, 独立的张量元数目就越少. 需要指出的是, $C_{ij}$ 的数目只与晶系有关, 而与晶系中具体的对称类型无关.

$\mathbf C$ 的逆矩阵 $\mathbf S$ 称为柔顺矩阵. 利用 $\mathbf S$ 可得到杨氏弹性模量的一般表达式. 我们用与[100], [010], [001]三个晶向的方向余弦来表示任意方向的杨氏模量. 设 $l_1, l_2, l_3$ 为空间某一方向与晶体主轴的方向余弦, 空间任一方向的杨氏模量 $E$ 的大小只与方向有关, 具体表达式如下

写为矩阵元素加和的形式为

一种较对称, 方便推导的形式为

这三种不同的表达形式, 可根据需要选择使用.

不同晶系的杨氏弹性模量

上面杨氏弹性模量的公式有些复杂, 好在除三斜晶系外, 大多数晶体都具有对称性. 考虑到晶体的对称性, 某些弹性常数必定为零, 而某些则相等, 所以对具有对称性的晶体, 相应的的杨氏模量公式简单些. 下面两张图总结了不同晶系刚度矩阵和柔顺矩阵的特点, 以及不同晶系杨氏弹性模量的公式, 后者为许多文献所引用.

1 三斜晶系(Triclinic system)

三斜晶系是所有七大晶系中对称性最低的晶系, 因此拥有最多的独立矩阵元, 其形式为:

共有21个独立的矩阵元, 杨氏模量

2 单斜晶系(Monoclinic system)

考虑对称性后, 单斜晶系有13个独立的矩阵单元:

3 正交晶系(Orthorhombic system)

正交晶系拥有相当高的对称性, 其独立矩阵元的数目为9个.

4 四方晶系(Tetragonal system)

4.1 四方晶系 $4, \bar 4, 4/m$

对于具有 $4, \bar4, 4/m$ 对称操作的四方晶系, 其独立矩阵元的数目为7个:

4.2 四方晶系 $422, 4mm, \bar 42m, 4/mmm$

对于具有 $422, 4mm, \bar 42m, 4/mmm$ 对称操作的四方晶系, 独立矩阵元的数目仅为6个:

5 三方晶系(Trigonal system)

5.1 三方晶系 $3, \bar 3$

三方晶系 $3, \bar 3$ 的独立矩阵元的数目为7个.

5.2 三方晶系 $32, 3m, \bar 3m$

三方晶系 $32, 3m, \bar 3m$ 独立矩阵元的数目为6个.

6 六方晶系(Hexagonal system)

六方晶系共有5个独立的矩阵元.

7 立方晶系(Cubic system)

立方晶系是所有晶系中对称度最高的晶系, 其独立矩阵元数目仅为3个, $C_{11}, C_{12}, C_{44}$

杨氏模量的极值

对于立方晶体, 我们可以用与[100], [010], [001]三个晶向的方向余弦来表示任意方向的杨氏模量, 结果如下

其中 $S_{11}, S_{12}, S_{44}$ 分别为立方晶体的三个独立的弹性柔顺系数, 柔顺矩阵 $S$ 与弹性矩阵 $C$ 的矩阵互为逆矩阵. $l_1, l_2, l_3$ 为空间某一方向与晶体主轴的方向余弦. $A$ 为各向异性值. 因此, 知道了三个柔顺弹性常数的值, 即可求得空间任一方向的杨氏模量 $E$, $E$ 的大小只与方向有关.

由于 $l_1^2+l_2^2+l_3^2=1, l_1, l_2, l_3 \in [0,1]$, 可以知道杨氏模量的两个极值为

前者对应于坐标轴方向, 后者对应于体对角线方法. 根据各向异性值 $A$ 与1的大小不同, 相应于极小或极大值.

对于其他晶系杨氏模量的极值, 不易得到解析公式, 直接使用数值方法搜索即可.

杨氏弹性模量各向异性的图示

为了直观地表达弹性模量的各向异性, 人们常常将其用三维图来表示. 这种各向异性的直观图示方法具有一般性, 在科学数据可视化中经常遇到. 量子化学中常用的原子轨道的角度分布图就是一例. 具体原理是, 在球坐标系 $(r,\q, \f)$ 中, 对仅依赖于方向的函数 $F(\q, \f)$ 中, 做曲面 $r=F(\q, \f)$. 显然, 当 $F$ 为常数时, 此曲面为球面, 各个方向函数值相同, 不存在各向异性; 当 $F$ 随 $(\q,\f)$ 变化时, 曲面便可表示出函数值的变化.

mathematica中可使用球坐标绘图函数SphericalPlot3D来做出这种图, 很多文献中的图就是利用这个函数做的, 请参考这个函数的说明和相应的弹性模量示例.

考虑到Matlab使用更广泛些, 下面给出基于Matlab的绘图方法.

利用Matlab绘制各向异性图时, 有两种实现方法. 一种是利用球坐标绘图, 像mathematica那样. 虽然Matlab没有球坐标绘图函数, 但可以先将球坐标转换为直角坐标然后再绘图, 也不是很麻烦. 另一种方法是直接使用直角坐标, 利用等值面函数, 绘制函数 $r-E=0$ 的等值面.

下面的代码绘制几种金属的杨氏模量三维各向异性曲面, 弹性常数来源于这里.

Aniso.m
  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
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
function Aniso
clc;clear;close all;

%% 处理数据, 计算矩阵以及弹性模量的极值

C=zeros(6);

% 单斜晶系测试
% C(1,1)=125;C(1,2)=87; C(1,3)=90; C(1,4)=0; C(1,5)=-9;C(1,6)=0;
%            C(2,2)=169;C(2,3)=105;C(2,4)=0; C(2,5)=-7;C(2,6)=0;
%                       C(3,3)=128;C(3,4)=0; C(3,5)=11;C(3,6)=0;
%                                  C(4,4)=53;C(4,5)=0; C(4,6)=-0.6;
%                                            C(5,5)=36;C(5,6)=0;
%                                                      C(6,6)=48;

% 正交晶系测试
C(1,1)=115.9; C(1,2)=35.3;  C(1,3)=46.8;
              C(2,2)=174.1; C(2,3)=38.7;
                            C(3,3)=153.1;
C(4,4)=50.9;  C(5,5)=70.2;  C(6,6)=26.6;

% 立方晶系
% C11=240.20; C12= 125.60; C44= 28.20; % Nb
% C11=522.40; C12= 160.80; C44=204.40; % W
% C11=107.30; C12=  60.90; C44= 28.30; % Al
% C11=346.70; C12= 250.70; C44= 76.50; % Pt
% C11=231.40; C12= 134.70; C44=116.40; % Fe
% C11=124.00; C12=  93.40; C44= 46.10; % Ag
% C11= 49.50; C12=  42.30; C44= 14.90; % Pb
% C11= 13.50; C12=  11.44; C44=  8.78; % Li
% C(1:3,1:3)=C12;
% for i=1:3; C(i,i)=C11; end
% for i=4:6; C(i,i)=C44; end

for i=2:6; for j=1:i-1; C(i,j)=C(j,i); end; end

S=inv(C);
S11=S(1,1); S12=S(1,2); S13=S(1,3); S14=S(1,4); S15=S(1,5); S16=S(1,6);
            S22=S(2,2); S23=S(2,3); S24=S(2,4); S25=S(2,5); S26=S(2,6);
                        S33=S(3,3); S34=S(3,4); S35=S(3,5); S36=S(3,6);
                                    S44=S(4,4); S45=S(4,5); S46=S(4,6);
                                                S55=S(5,5); S56=S(5,6);
                                                            S66=S(6,6);

%% 使用球坐标作图

[theta, phi]=meshgrid( linspace(0,pi), linspace(0,2*pi) );

L1=sin(theta).*cos(phi);
L2=sin(theta).*sin(phi);
L3=cos(theta);

% 三斜晶系杨氏模量公式, 可用于任意晶系
E=S11 * L1.^4 + S22 * L2.^4 + S33 * L3.^4 ...
 + (S44+2*S23) * (L2.*L3).^2 + (S55+2*S13) * (L1.*L3).^2 + (S66+2*S12) * (L1.*L2).^2 ...
 + 2*((S14+S56) * L1.^2 +       S24 * L2.^2 +       S34 * L3.^2) .* L2.*L3 ...
 + 2*(      S15 * L1.^2 + (S25+S46) * L2.^2 +       S35 * L3.^2) .* L1.*L3 ...
 + 2*(      S16 * L1.^2 +       S26 * L2.^2 + (S36+S45) * L3.^2) .* L1.*L2;

% 立方晶系
% E=S11+(1-A)*S44*( (L1.*L2).^2+(L2.*L3).^2+(L3.*L1).^2 );

E=1./E;

% 查找极值, 可用于任意晶系
A=0; Emin=min(E(:)); Emax=max(E(:));

% 立方晶系极值公式
% A=2*(S11-S12)/S44;
% Emax=1/S11; Emin=1/(S11+(1-A)*S44/3);
% if(A>1); Emin=1/S11; Emax=1/(S11+(1-A)*S44/3); end

fprintf('A=%9.4f Emin=%9.4f Emax=%9.4f\n', A, Emin, Emax);

x=E.*L1; y=E.*L2; z=E.*L3;
% 或使用函数转为直角坐标
% [x,y,z] = sph2cart(v, pi/2-u,E);

% 作颜色映射曲面
surf(x,y,z, E, 'FaceColor','interp', 'EdgeColor','none');

% 作曲面在做坐标面的投影, f 控制投影位置
hold on
xr=xlim; yr=ylim; zr=zlim; f=1.2;
mesh(zeros(size(x))+f*xr(1), y, z, E)
mesh(x, zeros(size(y))-f*yr(1), z, E)
mesh(x, y, zeros(size(z))+f*zr(1), E)
hold off

% 作模量的某一切面图, 使用前先注释掉上面的 surf 语句
% [X,Y,Z]=meshgrid(linspace(-Emax,Emax));
% contourslice(X,Y,Z,X,x,y,z,[0 0])

%% 或使用直角坐标等值面方法作图

% [x,y,z]=meshgrid(linspace(-Emax,Emax));
% r=sqrt(x.^2+y.^2+z.^2);
% L1=x./r; L2=y./r; L3=z./r;
%
% % 立方晶系
% E=S11+(1-A)*S44*( (L1.*L2).^2+(L2.*L3).^2+(L3.*L1).^2 );
% E=1./E;
% v=r-E;
%
% p=patch(isosurface(x,y,z,v,0));
% isocolors(x,y,z,E,p);
% isonormals(x,y,z,v,p);
% set(p,'FaceColor','interp','EdgeColor','none');

%% 设置图片格式, 输出图片

axis tight; title 'Nb A=0.5';
view(45,30); daspect([1 1 1]);
camlight; lighting phong;

colormap jet; % 低版本Matlab默认的填色模式
cbar=colorbar; title(cbar, 'GPa');

set(gca,'position',[0.12,0.05, 0.6,0.85]);
set(gcf,'position',[20,20, 1000,900]);
set(gcf, 'PaperPositionMode', 'auto');
print(gcf,'-dpng','-r300','Nb.png')

end

注意

  1. matlab默认的渲染颜色取决于Z轴大小, 这不符合我们的要求, 因为我们需要用颜色表示E的大小, 这样图形更直观.
  2. matlab球坐标转换函数使用的球坐标采用数学约定, 与物理上常用的不同, 使用仰角El, 而非俯视角 $\q$, $El+\q=\p/2$
  3. 不同晶系杨氏模量的表达式不同, 只要把代码里E的表达式修改成相应的方程即可.
  4. 杨氏模量在某一平面内的截面图形可利用极坐标或直角坐标绘制, 原理类似.

为了让大家有一个更直观的了解, 我们把具有不同各向异性值的立方金属选取具有代表性几个, 列于下表

几种金属的弹性数据(单位: GPa)
金属 C11 C12 C44 S11 S12 S44 A Emin Emax
铌Nb 240.20 125.60 28.20 0.0065 -0.0022 0.0355 0.49 80.01 153.95
钨W 522.40 160.80 204.40 0.0025 -0.0007 0.0062 1.13 446.71 493.65
铝Al 107.30 60.90 28.30 0.0158 -0.0057 0.0353 1.22 63.20 75.57
铂Pt 346.70 250.70 76.50 0.0073 -0.0031 0.0131 1.59 136.29 210.51
铁Fe 231.40 134.70 116.40 0.0076 -0.0028 0.0086 2.41 132.28 283.34
银Ag 124.00 93.40 46.10 0.0229 -0.0098 0.0217 3.01 43.75 120.44
铅Pb 49.50 42.30 14.90 0.0951 -0.0438 0.0671 4.14 10.52 40.23
锂Li 13.50 11.44 8.78 0.3328 -0.1526 0.1139 8.52 3.00 21.23

相应的三维杨氏模量图如下

参考资料

  1. T. C. T. Ting; On Anisotropic Elastic Materials for which Young’s Modulus E(n) is Independent of n or the Shear Modulus G(n,m) is Independent of n and m; J Elasticity 81(3):271-292, 2006; 10.1007/s10659-005-9016-2
  2. Kevin M. Knowles, Philip R. Howie; The Directional Dependence of Elastic Stiffness and Compliance Shear Coefficients and Shear Moduli in Cubic Materials; J Elast 120(1):87-108, 2014; 10.1007/s10659-014-9506-1
  3. Matlab绘图高级部分
  4. 科学计算可视化: 三维流场绘图
  5. Applied Mechanics of Solids: Constitutive Models Relations between Stress and Strain
  6. 球坐标

评论

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: gnuplot颜色设置
后一篇: 希腊字母与拉丁字母的对应

访问人次(2015年7月 9日起): | 最后更新: 2017-06-16 03:43:39 UTC | 版权所有 © 2008 - 2017 Jerkwin