版主 (知名作家)
认真做事，踏实做人 ![]()
|
【答案】应助回帖
★ ★
obaica: 金币+2, ★★★★★最佳答案 2016-07-21 19:10:51
最后声明：材料转自http://jerkwin.github.io/2014/04/17/材料弹性模量各向异性的三维图示方法/ 如有不妥，望作者通知，我进行删除。我只是数据的搬运工，不是原创者。排版可能有乱码，请参照原文。
杨氏模量的极值
对于立方晶体, 我们可以用与[100], [010], [001]三个晶向的方向余弦来表示任意方向的杨氏模量, 结果如下
![杨氏模量投影图-12]()
其中 S11,S12,S44分别为立方晶体的三个独立的弹性柔顺系数, 柔顺矩阵 S 与弹性矩阵 C 的矩阵互为逆矩阵. l1,l2,l3为空间某一方向与晶体主轴的方向余弦. A 为各向异性值. 因此, 知道了三个柔顺弹性常数的值, 即可求得空间任一方向的杨氏模量 E, E 的大小只与方向有关.
由于 l21+l22+l23=1,l1,l2,l3∈[0,1] 可以知道杨氏模量的两个极值为
![杨氏模量投影图-13]()
前者对应于坐标轴方向, 后者对应于体对角线方法. 根据各向异性值 AA 与1的大小不同, 相应于极小或极大值.对于其他晶系杨氏模量的极值, 不易得到解析公式, 直接使用数值方法搜索即可.
杨氏弹性模量各向异性的图示
为了直观地表达弹性模量的各向异性, 人们常常将其用三维图来表示. 这种各向异性的直观图示方法具有一般性, 在科学数据可视化中经常遇到. 量子化学中常用的原子轨道的角度分布图就是一例. 具体原理是, 在球坐标系 (r,θ,ϕ 中, 对仅依赖于方向的函数 F(θ,ϕ 中, 做曲面 r=F(θ,ϕ. 显然, 当 F 为常数时, 此曲面为球面, 各个方向函数值相同, 不存在各向异性; 当 FF 随 (θ,ϕ 变化时, 曲面便可表示出函数值的变化.
mathematica中可使用球坐标绘图函数SphericalPlot3D来做出这种图, 很多文献中的图就是利用这个函数做的, 请参考这个函数的说明和相应的弹性模量示例.
考虑到Matlab使用更广泛些, 下面给出基于Matlab的绘图方法.
利用Matlab绘制各向异性图时, 有两种实现方法. 一种是利用球坐标绘图, 像mathematica那样. 虽然Matlab没有球坐标绘图函数, 但可以先将球坐标转换为直角坐标然后再绘图, 也不是很麻烦. 另一种方法是直接使用直角坐标, 利用等值面函数, 绘制函数 r−E=0r−E=0 的等值面.
下面的代码绘制几种金属的杨氏模量三维各向异性曲面, 弹性常数来源于这里.
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;
% for i=2:6; for j=1:i-1; C(i,j)=C(j,i); end; end
% 立方晶系
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=zeros(6);
C(1:3,1:3)=C12;
for i=1:3; C(i,i)=C11; end
for i=4:6; C(i,i)=C44; 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);
% 立方晶系极值公式
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);
%% 使用球坐标作图
[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;
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');
% 作模量的某一切面图
% [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, 而非俯视角
θθ, El+θ=π/2
![杨氏模量投影图-14]()
3. 不同晶系杨氏模量的表达式不同, 只要把代码里E的表达式修改成相应的方程即可.
4. 杨氏模量在某一平面内的截面图形可利用极坐标或直角坐标绘制, 原理类似.
为了让大家有一个更直观的了解, 我们把具有不同各向异性值的立方金属选取具有代表性几个, 列于下表
几种金属的弹性数据(单位: GPa)
![杨氏模量投影图-15]()
相应的三维杨氏模量图如下
![杨氏模量投影图-16]()
参考资料
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
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
Matlab绘图高级部分
科学计算可视化: 三维流场绘图
Applied Mechanics of Solids: Constitutive Models Relations between Stress and Strain
球坐标 |
|
送红花一朵 
回复此楼
