Matlab中有關(guān)三維楊氏模量的繪圖
? ? ? ?之前幫研究第一性原理的同學(xué)在matlab中繪制過(guò)有關(guān)三維楊氏模量的圖，代碼是從論壇copy下來(lái)的，根據(jù)那位同學(xué)的實(shí)際要求把代碼稍微改了一下。
? ? ? ? 關(guān)于原理的介紹：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/
根據(jù)需要修改過(guò)的代碼：
function CrystalModulus
clc; clear; close all;
%% 設(shè)置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矩陣對(duì)稱元素
for i=2:6; for j=1:i-1; C(i,j)=C(j,i); end; end
%% 計(jì)算S矩陣
S=inv(C);
%% 球坐標(biāo)(二維)表面作圖方法
[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); % 或使用函數(shù)將球坐標(biāo)轉(zhuǎn)為直角坐標(biāo)
fprintf('%s: Aniso=%9.4f\n Emin=%9.4f\n Emax=%9.4f\n', Name, A, Vmin, Vmax);
SphPlt(X, Y, Z, V);
%% 設(shè)置查看方式
axis tight; title(sprintf('%s? Aniso=%9.4f',Name,A));
daspect([1 1 1]);
view(45,30); % 查看角度, 根據(jù)需要更改. (30,30)可能更合適
colormap jet; % 低版本Matlab默認(rèn)的填色模式
? ? caxis([Vmin Vmax]);
cbar=colorbar; title(cbar, 'GPa');
camlight; lighting phong;
%% 設(shè)置圖片格式并輸出
set(gca,'position',[0.12,0.05, 0.6,0.85]);
set(gcf,'position',[500,500, 380,350]); % 根據(jù)需要更改大小如[20,20,1000,900]
set(gcf, 'PaperPositionMode', 'auto');
print(gcf,'-dpng','-r600', [Name,'.png'])
end
%% 由C矩陣以及方向矢量L1, L2, L3計(jì)算彈性模量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;
? ? %% 數(shù)值查找極值, 可用于任意晶系
Emin=min(E(:)); Emax=max(E(:));
end
%%% 繪制二維表面圖
function SphPlt(X, Y,Z, V)
% %? %% 作顏色映射曲面
% surf(X,Y,Z,V, 'FaceColor','interp','EdgeColor','none');
%? 作曲面在做坐標(biāo)面的投影
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語(yǔ)句
% %% 參考 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
? ? ? ? ? 歡迎評(píng)論區(qū)討論
