产品
之前帮研究第一性原理的同学在matlab中绘制过有关三维杨氏模量的图，代码是从论坛copy下来的，根据那位同学的实际要求把代码稍微改了一下。
关于原理的介绍：http://jerkwin.github.io/2014/04/17/%E6%9D%90%E6%96%99%E5%BC%B9%E6%80%A7%E6%A8%A1%E9%87%8F%E5%90%84%E5%90%91%E5%BC%82%E6%80%A7%E7%9A%84%E4%B8%89%E7%BB%B4%E5%9B%BE%E7%A4%BA%E6%96%B9%E6%B3%95/#opennewwindow
原代码：http://jerkwin.github.io/2014/10/03/%E5%89%AA%E5%88%87%E6%A8%A1%E9%87%8F%E5%90%84%E5%90%91%E5%BC%82%E6%80%A7/
根据需要修改过的代码：
function CrystalModulus
clc; clear; close all;
%% 设置C矩阵
C=zeros(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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 单斜晶系示例
Name='Fe1Mn2B4';
C(1,1)=440.4; C(1,2)= 149.6; C(1,3)= 149.9; C(1,4)= 0; C(1,5)=0; C(1,6)= 0;
C(2,2)=449.0; C(2,3)=150.1; C(2,4)= 0; C(2,5)=0; C(2,6)= 0;
C(3,3)=378.8; C(3,4)= 0; C(3,5)=0; C(3,6)= 0;
C(4,4)=176.7; C(4,5)= 0; C(4,6)=0;
C(5,5)=175.9; C(5,6)= 0;
C(6,6)=167.3;
C11=C(1,1); C12=C(1,2); C13=C(1,3); C14=C(1,4); C15=C(1,5); C16=C(1,6);
C22=C(2,2); C23=C(2,3); C24=C(2,4); C25=C(2,5); C26=C(2,6);
C33=C(3,3); C34=C(3,4); C35=C(3,5); C36=C(3,6);
C44=C(4,4); C45=C(4,5); C46=C(4,6);
C55=C(5,5); C56=C(5,6);
C66=C(6,6);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 赋值C矩阵对称元素
for i=2:6; for j=1:i-1; C(i,j)=C(j,i); end; end
%% 计算S矩阵
S=inv(C);
%% 球坐标(二维)表面作图方法
[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(S, L1, L2, L3); % 弹性模量
X=V.*L1; Y=V.*L2; Z=V.*L3;
% [X,Y,Z] = sph2cart(phi, pi/2-theta, V); % 或使用函数将球坐标转为直角坐标
fprintf('%s: Aniso=%9.4f\n Emin=%9.4f\n Emax=%9.4f\n', Name, A, Vmin, Vmax);
SphPlt(X, Y, Z, V);
%% 设置查看方式
axis tight; title(sprintf('%s Aniso=%9.4f',Name,A));
daspect([1 1 1]);
view(45,30); % 查看角度, 根据需要更改. (30,30)可能更合适
colormap jet; % 低版本Matlab默认的填色模式
caxis([Vmin Vmax]);
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(S, L1, L2, L3)
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; % 只适用于立方晶系
E=S11+(1-A)*S44*( (L1.*L2).^2+(L2.*L3).^2+(L3.*L1).^2 );
E=1./E;
%% 数值查找极值, 可用于任意晶系
Emin=min(E(:)); Emax=max(E(:));
end
%%% 绘制二维表面图
function SphPlt(X, Y,Z, V)
% % %% 作颜色映射曲面
% surf(X,Y,Z,V, 'FaceColor','interp','EdgeColor','none');
% 作曲面在做坐标面的投影
hold on
f=1.5; % 控制投影位置
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)
%% 作模量的某一切面图, 使用前最好注释掉上面的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]);
% %% 将投影面绘制于xy平面
% k = boundary(X(:),Y(:), .5);
% plot(X(k),Y(k), '-r', 'linewidth', 2);
% k = boundary(X(:),Z(:), 1);
% plot(X(k),Z(k), '-g', 'linewidth', 2);
% k = boundary(Y(:),Z(:), 1);
% plot(Y(k),Z(k), '-b', 'linewidth', 2);
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
六方三维表面图
六方二维剖面图
单斜三维图
欢迎评论区讨论
155-2731-8020
