以前的博文中讲过用matlab作弹性模量(杨氏模量)的各向异性图，这里更进一步, 讲讲剪切模量各向异性的作图。

我们知道，剪切模量定义为: 在某个晶面上沿着某一特定方向施加剪切应力后剪切应力和应变的比值。它的值是和两个方向有关的.

在图中，与向量 l 垂直的平面即是晶体要发生剪切的平面. 由于弹性模量的各向异性，在这一特定平面的不同方向， 施加相同的剪切应力会得到不同的剪切应变，故会得到不同的剪切模量的值。此外, 与杨氏模量不同的是，剪切模量在某一方向 l 上的值还与角度 χ 有关. 这就导致剪切模量的各向异性比杨氏模量更加复杂。$;$

虽然剪切模量的各项异性比杨氏模量复杂, 但是作图的方法没有什么本质区别. 作图的关键是它的计算.

以立方晶系为例, 其杨氏模量

相应的剪切模量为

其中

通常对剪切模量进行表征时使用其在每一方向上最大值或最小值, 这就需要我们计算剪切模量某一方向的极值.

若令 α=4S11−S44,β=4S12+S44,γ=S44, 则 G 可写为

其中

可见 G 为 m1,m2,m3 的二次型, 可根据其正定性判断极值情况.

根据极值的必要条件, 也可对 G 进行求导, 根据 m1,m2,m3 方程组解的情况进行判断极值情况. 求导后

其行列式

这两种情形下, 得到的极值点都为 (0,0,0), 但显然无法满足 m21+m22+m23=1 的约束条件, 所以这两种方法都不可行.

换用最笨的方法, 直接将 mi 的表达式代入, 化简后得到 G 满足的方程. 设

则

最终表示为简单的三角函数, 容易知道, G 的最大值

Gmax=C2+C32+C21+(C2−C3)24−−−−−−−−−−−√

最小值

Gmin=C2+C32−C21+(C2−C3)24−−−−−−−−−−−√

上面的推导过程可借助matlab的符号计算功能完成.

虽然有了解析的表达式, 但实际编码时很繁琐, 且大多数情况下我们并不能得到解析的公式. 所以这只是作为求极值的一种方法. 常用的数值方法中最简单的就是分割点直接计算所有值, 然后求其最大值或最小值, 复杂一点的就是利用优化的各种方法.

照例, 下面给出几种典型金属的剪切模量最大值与最小值的图像, 前者为最大值, 后者为最小值.

作图代码# Language: matlabfunction Aniso    clc; clear; close all;    global S11 S12 S44%     C11=240.20; C12=125.60; C44=28.20; Name='Nb'; Aniso='0.49';%     C11=522.40; C12=160.80; C44=204.40; Name='W';  Aniso='1.13';%     C11=107.30; C12=60.90; C44=28.30; Name='Al'; Aniso='1.22';%     C11=346.70; C12=250.70; C44=76.50; Name='Pt'; Aniso='1.59';%     C11=231.40; C12=134.70; C44=116.40; Name='Fe'; Aniso='2.41';%     C11=124.00; C12=93.40; C44=46.10; Name='Ag'; Aniso='3.01';%     C11=49.50; C12=42.30; C44=14.90; Name='Pb'; Aniso='4.14';    C11=13.50; C12=11.44; C44=8.78; Name='Li'; Aniso='8.52';    C=zeros(6);C(1:3,1:3)=C12;for i=1:3;C(i,i)=C11; endfor i=4:6;C(i,i)=C44; end    S=inv(C);    S11=S(1,1); S12=S(1,2); S44=S(4,4);%[X, Y, Z, E]= Eval;%SphPlt(X, Y, Z, E,[Name,'_E.png'],[Name,' Aniso=', Aniso]);%[Gmin, Gmax]=pltGchi(.25*pi,.25*pi);[X, Y, Z, G]= Gval;SphPlt(X, Y, Z, G,[Name,'_Gmin.png'],[Name,' Aniso=', Aniso]);end%%function[X, Y, Z, E]= Eval    global S11 S12 S44    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); endfprintf('Aniso=%9.4f Emin=%9.4f Emax=%9.4fn', 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+(1-A)*S44*((L1.*L2).^2+(L2.*L3).^2+(L3.*L1).^2);    E=1./E;    X=E.*L1; Y=E.*L2; Z=E.*L3;end%%function[X, Y, Z, G]= Gval[theta, phi]=meshgrid(linspace(0,pi,50),linspace(0,2*pi));    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(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));    X=G.*L1; Y=G.*L2; Z=G.*L3;end%%function G =Gchi(theta, phi, chi)    global S11 S12 S44    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%%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%%functionSphPlt(X, Y, Z, F, File, Title)surf(X,Y,Z,F,'FaceColor','interp','EdgeColor','none');    axis tight;title(Title);view(30,30);daspect([111]);    camlight; lighting phong;    cbar=colorbar;title(cbar,'GPa');set(gca,'position',[0.12,0.05,0.6,0.85]);set(gcf,'position',[500,500,380,350]);set(gcf,'PaperPositionMode','auto');print(gcf,'-dpng','-r300', File)end

辅助推导的代码

# Language: matlabclc;syms A B C D E F G m1 m2 m3 P Q R S T xm1 = P*cos(x)+Q*sin(x);m2 = R*cos(x)+S*sin(x);m3 = T*cos(x);G = A*m1^2+ B*m2^2+ C*m3^2+2*D*m1*m2 +2*E*m1*m3 +2*F*m2*m3;G =expand(G);G =simple(G);G=collect(G,cos(x));G=collect(G,sin(x));G=collect(G,sin(2*x));Gpretty(G)

◆图片/表格/公式/代码完整版请参看：剪切模量各向异性◆