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

类别:    标签: 数理 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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
function aniso
	clc; clear; clear all; close all;

	%% 设置C矩阵
	C=zeros(6);

	% 单斜晶系示例
	% Name='Tric';
	% 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;

	% 六方晶系示例
	Name='Hex';
	C(1,1)=581;    C(1,2)=55;     C(1,3)=121;    C(3,3)=445; C(4,4)=240;
	C(2,2)=C(1,1); C(2,3)=C(1,3); C(5,5)=C(4,4); C(6,6)=(C(1,1)-C(1,2))/2;

	% 正交晶系示例
	% Name='Orth';
	% 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; Name='Nb';
	% C11=522.40; C12= 160.80; C44=204.40; Name='W ';
	% C11=107.30; C12= 60.90; C44= 28.30; Name='Al';
	% C11=346.70; C12= 250.70; C44= 76.50; Name='Pt';
	% C11=231.40; C12= 134.70; C44=116.40; Name='Fe';
	% C11=124.00; C12=  93.40; C44= 46.10; Name='Ag';
	% C11= 49.50; C12=  42.30; C44= 14.90; Name='Pb';
	% C11= 13.50; C12=  11.44; C44=  8.78; Name='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

	% 赋值C矩阵对称元素
	for i=2:6; for j=1:i-1; C(i,j)=C(j,i); end; end

	%% 球坐标(二维)表面作图方法
	[theta, phi]=meshgrid( linspace(0,pi), linspace(0,2*pi) );
	L1=sin(theta).*cos(phi);
	L2=sin(theta).*sin(phi);
	L3=cos(theta);

	[V, A, Vmin, Vmax]=getE(C, L1, L2, L3); % 弹性模量
	% [V, A, Vmin, Vmax] = getG(C, theta, phi); % 剪切模量

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

	fprintf('%s: Aniso=%9.4f Min=%9.4f Max=%9.4f\n', Name, A, Vmin, Vmax);
	SphPlt(x, y, z, V);

	%% 直角坐标(三维)等值面作图方法, 慢, 仅供参考
	% [x,y,z]=meshgrid(linspace(-Vmax,Vmax));
	% r=sqrt(x.^2+y.^2+z.^2);
	% L1=x./r; L2=y./r; L3=z./r;
	% theta=acos(L3); phi=atan2(L2,L1);

	% [V, A, Vmin, Vmax] = getE(C, L1, L2, L3);
	% [V, A, Vmin, Vmax] = getG(C, theta, phi);
	% fprintf('%s: Aniso=%9.4f Emin=%9.4f Emax=%9.4f\n', Name, A, Vmin, Vmax);
	% CartPlt(x, y, z, V-r, V);

	%% 设置查看方式
	axis tight; title(sprintf('%s  Aniso=%9.4f',Name,A));
	daspect([1 1 1]);
	view(45,30); % 查看角度, 根据需要更改. (30,30)可能更合适

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

	camlight; lighting phong;

	%% 设置图片格式并输出
	set(gca,'position',[0.12,0.05, 0.6,0.85]);
	set(gcf,'position',[500,500, 380,350]); % 根据需要更改大小如[20,20,1000,900]
	set(gcf, 'PaperPositionMode', 'auto');
	print(gcf,'-dpng','-r600', [Name,'.png'])
end

%% 由C矩阵以及方向矢量L1, L2, L3计算弹性模量E
function [E, A, Emin, Emax] = getE(C, L1, L2, L3)
	S=inv(C); % 计算S矩阵
	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);

	% 三斜晶系杨氏模量通用公式, 可用于任意晶系
	A=0;
	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;
	% 立方晶系杨氏模量公式, 用于和通用公式对比, 测试二者是否一致
	% A=2*(S11-S12)/S44;
	% E=S11+(1-A)*S44*( (L1.*L2).^2+(L2.*L3).^2+(L3.*L1).^2 );
	% 六方晶系杨氏模量公式, 用于和通用公式对比, 测试二者是否一致
	% E=S11*(1-L3.^2).^2+S33*L3.^4+(S44+2*S13)*(1-L3.^2).*L3.^2;

	E=1./E;

	% 数值查找极值, 可用于任意晶系
	Emin=min(E(:)); Emax=max(E(:));
	% 立方晶系极值公式, 用于和数值极值对比, 测试二者是否一致
	% Emax=1/S11; Emin=1/(S11+(1-A)*S44/3);
	% if(A>1); Emin=1/S11; Emax=1/(S11+(1-A)*S44/3); end
end

%% 由C矩阵以及球坐标theta, phi计算剪切模量G极值, 只适用于立方晶系
function [G, A, Gmin, Gmax] = getG(C, theta, phi)
	S=inv(C); % 计算S矩阵
	A=2*(S(1,1)-S(1,2))/S(4,4);

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

	theta = theta(:); phi = phi(:);
	G = zeros(length(phi),1);
	for i=1:length(G)
		[x, Gmin] = fminbnd(@(x) Gchi(S, theta(i), phi(i), x), 0,pi); G(i) = Gmin;
		% [x, Gmax] = fminbnd(@(x) -Gchi(theta(i), phi(i), x), 0,pi); G(i) = -Gmax;
	end
	G = reshape(G, size(L1));
	Gmin=min(G(:)); Gmax=max(G(:));
end

%% 立方晶系剪切模量计算公式
function G = Gchi(S, theta, phi, chi)
	S11=S(1,1); S12=S(1,2); S44=S(4,4);

	L1=sin(theta).*cos(phi);
	L2=sin(theta).*sin(phi);
	L3=cos(theta);
	M1=cos(theta).*cos(phi).*cos(chi)-sin(phi).*sin(chi);
	M2=cos(theta).*sin(phi).*cos(chi)+cos(phi).*sin(chi);
	M3= -sin(theta).*cos(chi);

	G = 4*S11*( (L1.*M1).^2+(L2.*M2).^2+(L3.*M3).^2 ) ...
	  + 8*S12*( L1.*L2.*M1.*M2+L1.*L3.*M1.*M3+L2.*L3.*M2.*M3) ...
	  +   S44*( (L2.*M3+M2.*L3).^2+(L1.*M3+M1.*L3).^2+(L1.*M2+M1.*L2).^2 );
	G = 1./G;
end

%% 剪切模量theta, phi方向绘图, 极值
function [Gmin, Gmax] = pltGchi(theta, phi)
	chi=linspace(0, 2*pi, 200);
	G = Gchi(theta, phi, chi);
	Gmin=min(G(:));
	Gmax=max(G(:));
	plot(chi,G);
end

%% 绘制二维表面图
function SphPlt(x, y, z, v)
	% 作颜色映射曲面
	surf(x,y,z,v, 'FaceColor','interp', 'EdgeColor','none');

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

	% 作模量的某一切面图, 使用前最好注释掉上面的surf和mesh语句
	% 参考 https://www.zhihu.com/question/48734216
	% vmax=max(v(:));
	% [X,Y,Z]=meshgrid(linspace(-vmax,vmax));
	% contourslice(X,Y,Z,  X, x,y,z,[0 0]);
	% contourslice(X,Y,Z, -Y, x,y,z,[0 0]);
	% contourslice(X,Y,Z,  Z, x,y,z,[0 0]);
end

%% 绘制三维等值面图
function CartPlt(x, y, z, v, c)
	p=patch(isosurface(x,y,z,v,0));
	isocolors(x,y,z,c,p);
	isonormals(x,y,z,v,p);
	set(p,'FaceColor','interp', 'EdgeColor','none');
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-09-26 07:50:25 UTC | 版权所有 © 2008 - 2017 Jerkwin