剪切模量各向异性

2023-08-11 09:34| 来源: 网络整理| 查看: 265

2014-10-03 10:10:58 硕士生伊国辉对此文有所贡献 2017-07-19 10:32:41 与弹性模量代码合并 2018-07-20 09:21:44 三斜晶系剪切模量代码, 感谢陈磊提供信息

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

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

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

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

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

\[E= S_{11}-2(S_{11}-S_{12}-{1 \over 2}S_{44})(l_1^2l_2^2+l_1^2l_3^2+l_2^2l_3^2)\]

相应的剪切模量为

\[\begin{align} G &=4S_{11}(l_1^2m_1^2+l_2^2m_2^2+l_3^2m_3^2) \\ &+ 8S_{12}(l_1l_2m_1m_2 +l_1l_3m_1m_3+l_2l_3m_2m_3) \\ &+ S_{44}[(l_1m_2+l_2m_1)^2+(l_1m_3+l_3m_1)^2+(l_2m_3+l_3m_2)^2] \\ &= A m_1^2 + B m_2^2 + C m_3^2 + 2D m_1m_2 + 2E m_1m_3 +2F m_2m_3 \end{align}\]

其中

\[l=\left( \begin{matrix} \sin \theta \cos \phi \\ \sin \theta \sin \phi \\ \cos \theta \\ \end{matrix} \right);\quad m=\left( \begin{matrix} \cos \theta \cos \phi \cos \chi -\sin \phi \sin \chi \\ \cos \theta \sin \phi \cos \chi +\cos \phi \sin \chi \\ -\sin \theta \cos \chi \\ \end{matrix} \right)\]

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

若令 $\a=4S_{11}-S_{44}, \b=4S_{12}+S_{44}, \g=S_{44}$, 则 $G$ 可写为

\[G= A m_1^2 + B m_2^2 + C m_3^2 + 2D m_1m_2 + 2E m_1m_3 +2F m_2m_3\]

其中

\[\begin{align} A &= \a l_1^2+\g=(4S_{11}-S_{44})l_1^2+S44 \\ B &= \a l_1^2+\g \\ C &= \a l_1^2+\g \\ D &= \b l_1l_2 = (4S_{12}+S_{44})l_1l_2 \\ E &= \b l_1l_3 \\ F &= \b l_2l_3 \end{align}\]

可见 $G$ 为 $m_1, m_2, m_3$ 的二次型, 可根据其正定性判断极值情况.

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

\[\begin{bmatrix} A & D & E \\ D & B & F \\ E & F & C \\ \end{bmatrix} \left[ \begin{array}{c} m_1 \\ m_2 \\ m_3 \end{array} \right] =0\]

其行列式

\[\D = (\a^3+2\b^3-3\a\b^2)l_1^2l_2^2l_3^2 + \g(\a^2-\b^2)(l_1^2l_2^2+l_1^2l_3^2+l_2^2l_3^2)+\a\g^2+\g^3\]

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

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

\[\alg m_1 &= \cos \theta \cos \phi \cos \chi -\sin \phi \sin \chi = P \cos\chi + Q \sin\chi \\ m_2 &= \cos \theta \sin \phi \cos \chi +\cos \phi \sin \chi = R \cos\chi + S \sin\chi \\ m_3 &= -\sin \theta \cos \chi = T \cos\chi \ealg\]

则

\[\alg G &= (APQ + BRS + DPS + DQR + EQT + FST) \sin(2x) \\ &+ (AP^2 + BR^2 + CT^2 + 2DPR + 2EPT + 2FRT ) \cos^2 x \\ &+ (AQ^2 + BS^2 + 2DQS ) \sin^2 x \\ &= C_1 \sin(2x) + C_2 \cos^2 x + C_3 \sin^2 x \\ &= C_1 \sin(2x) + (C_2-C3) \cos^2 x + C_3 \\ &= C_1 \sin(2x) + {C_2-C_3 \over 2} \cos(2x) + {C_2+C_3 \over 2} \\ &= \sqrt{C_1^2+{(C_2-C_3)^2 \over 4}} \sin(2x+\f) + {C_2+C_3 \over 2} , \f= \arctan {C_2-C_3 \over 2 C_1} \ealg\]

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

$G_{max}={C_2+C_3 \over 2}+\sqrt{C_1^2+{(C_2-C_3)^2 \over 4}}$

最小值

$G_{min}= {C_2+C_3 \over 2} - \sqrt{C_1^2+{(C_2-C_3)^2 \over 4}}$

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

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

推广一下, 对三斜晶系, 剪切模量计算公式为

\[\alg 1/G = &4 [S_{11} (l_1 m_1)^2 + S_{22} (l_2 m_2)^2 + S_{33} (l_3 m_3)^2] \\ +&S_{44} (l_2 m_3+l_3 m_2)^2 + S_{55} (l_1 m_3+l_3 m_1)^2 + S_{66} (l_1 m_2+l_2 m_1)^2 \\ +&8 [ S_{12} (l_1 l_2 m_1 m_2) + S_{13} (l_1 l_3 m_1 m_3) + S_{23} (l_2 l_3 m_2 m_3) ] \\ +&4 l_1 m_1 [S_{14} (l_2 m_3+l_3 m_2)+S_{15} (l_1 m_3+l_3 m_1)+S_{16} (l_1 m_2+l_2 m_1)] \\ +&4 l_2 m_2 [S_{24} (l_2 m_3+l_3 m_2)+S_{25} (l_1 m_3+l_3 m_1)+S_{26} (l_1 m_2+l_2 m_1)] \\ +&4 l_3 m_3 [S_{34} (l_2 m_3+l_3 m_2)+S_{35} (l_1 m_3+l_3 m_1)+S_{36} (l_1 m_2+l_2 m_1)] \\ +&2 \quad \quad \;S_{45} (l_2 m_3+l_3 m_2) (l_1 m_3+l_3 m_1) \\ +&2 \quad \quad \;S_{46} (l_2 m_3+l_3 m_2) (l_1 m_2+l_2 m_1) \\ +&2 \quad \quad \;S_{56} (l_1 m_3+l_3 m_1) (l_1 m_2+l_2 m_1) \\ = & S_{11} A_{11}^2 + S_{22} A_{22}^2 + S_{33} A_{33}^2 + S_{44} A_{23}^2 + S_{55} A_{13}^2 + S_{66} A_{12}^2 \\ +&2 [ S_{12} A_{11} A_{22} + S_{13} A_{11} A_{33}+ S_{23} A_{22} A_{33} ] \\ +&2 A_{11} [S_{14} A_{23} + S_{15} A_{13}+S_{16} A_{13}] \\ +&2 A_{22} [S_{24} A_{23} + S_{25} A_{13}+S_{26} A_{13}] \\ +&2 A_{33} [S_{34} A_{23} + S_{35} A_{13}+S_{36} A_{13}] \\ +&2 S_{45} A_{23} A_{13} + 2 S_{46} A_{23} A_{12} + 2 S_{56} A_{13} A_{12} \\ \\ A_{ij}=&l_i m_j+l_j m_i \ealg\]

因为 $\vec l$ 与 $\vec m$ 垂直, 因此其点积 $l_1m_1+l_2m_2+l_3m_3=0$, 利用此条件, 可以得到上式的另一种形式:

\[\alg 1/G= &4 [2 S_{12}-(S_{11}+S_{22}-S_{66})] l_1 m_1 l_2 m_2 \\ +&4 [2 S_{23}-(S_{22}+S_{33}-S_{44})] l_2 m_2 l_3 m_3 \\ +&4 [2 S_{13}-(S_{33}+S_{11}-S_{55})] l_1 m_1 l_3 m_3 \\ +&4 (l_1 m_2+l_2 m_1) [(S_{16}-S_{36}) l_1 m_1+(S_{26}-S_{36}) l_2 m_2] \\ +&4 (l_2 m_3+l_3 m_2) [(S_{24}-S_{14}) l_2 m_2+(S_{34}-S_{14}) l_3 m_3] \\ +&4 (l_1 m_3+l_3 m_1) [(S_{35}-S_{25}) l_3 m_3+(S_{15}-S_{25}) l_1 m_1] \\ +& \;\; S_{44} (l_2 m_3-l_3 m_2)^2 \\ +& \;\; S_{55} (l_3 m_1-l_1 m_3)^2 \\ +& \;\; S_{66} (l_1 m_2-l_2 m_1)^2 \\ +&2 S_{45} (l_2 m_3+l_3 m_2) (l_1 m_3+l_3 m_1) \\ +&2 S_{46} (l_2 m_3+l_3 m_2) (l_1 m_2+l_2 m_1) \\ +&2 S_{56} (l_1 m_3+l_3 m_1) (l_1 m_2+l_2 m_1) \ealg\]

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

作图代码 CrystalModulus.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 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211function CrystalModulus clc; clear; 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 %% 计算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); % 弹性模量 [V, A, Vmin, Vmax] = getG(S, theta, phi, -1); % 剪切模量 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(S, L1, L2, L3); % [V, A, Vmin, Vmax] = getG(S, theta, phi); % fprintf('%s: Aniso=%9.4f Min=%9.4f Max=%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默认的填色模式 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']) %% 剪切模量theta, phi方向随chi变化曲线 % theta=0.1*pi; phi=0.1*2*pi; % chi=linspace(0, 2*pi, 200); % G = Gchi(S, theta, phi, chi); % plot(chi,G); 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 * 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=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极值 %% minmax设定极值, -1: 最小值; 1: 最大值 function [G, A, Gmin, Gmax] = getG(S, theta, phi, minmax) A=2*(S(1,1)-S(1,2))/S(4,4); % 只适用于立方晶系 Gsize=size(theta); theta = theta(:); phi = phi(:); G = zeros(length(phi),1); for i=1:length(G) [x, ret] = fminbnd(@(x) -minmax*Gchi(S, theta(i), phi(i), x), 0,pi); G(i) = -minmax*ret; end G = reshape(G, Gsize); Gmin=min(G(:)); Gmax=max(G(:)); end %% 剪切模量计算公式 function G = Gchi(S, theta, phi, chi) 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); 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 + S22*(L2.*M2).^2 + S33*(L3.*M3).^2 ) ... + S44*(L2.*M3+L3.*M2).^2 + S55*(L1.*M3+L3.*M1).^2 + S66*(L1.*M2+L2.*M1).^2 ... +8*( S12*(L1.*L2.*M1.*M2) + S13*(L1.*L3.*M1.*M3) + S23*(L2.*L3.*M2.*M3) ) ... +4*(L1.*M1).*( S14*(L2.*M3+L3.*M2) + S15*(L1.*M3+L3.*M1) + S16*(L1.*M2+L2.*M1) ) ... +4*(L2.*M2).*( S24*(L2.*M3+L3.*M2) + S25*(L1.*M3+L3.*M1) + S26*(L1.*M2+L2.*M1) ) ... +4*(L3.*M3).*( S34*(L2.*M3+L3.*M2) + S35*(L1.*M3+L3.*M1) + S36*(L1.*M2+L2.*M1) ) ... +2* S45*(L2.*M3+L3.*M2).*(L1.*M3+L3.*M1) ... +2* S46*(L2.*M3+L3.*M2).*(L1.*M2+L2.*M1) ... +2* S56*(L1.*M3+L3.*M1).*(L1.*M2+L2.*M1); %% 立方晶系剪切模量公式, 用于和通用公式对比, 测试二者是否一致 % 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 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) %% 作模量的某一切面图, 使用前最好注释掉上面的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 辅助推导的代码

立方晶系剪切模量极值

Gcub.m 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19clc; syms A B C D E F G m1 m2 m3 P Q R S T x m1 = 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)); G pretty(G)

剪切模量公式一致性

dG.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 46clc; clear; syms G1 G2 dG L1 L2 L3 M1 M2 M3 ... S11 S12 S13 S14 S15 S16 S22 S23 S24 S25 S26 ... S33 S34 S35 S36 S44 S45 S46 S55 S56 S66 G1=4*(2*S12-(S11+S22-S66))*L1*M1*L2*M2 ... + 4*(2*S23-(S22+S33-S44))*L2*M2*L3*M3 ... + 4*(2*S13-(S33+S11-S55))*L3*M3*L1*M1 ... + 4*(L1*M2+L2*M1)*((S16-S36)*L1*M1+(S26-S36)*L2*M2) ... + 4*(L2*M3+L3*M2)*((S24-S14)*L2*M2+(S34-S14)*L3*M3) ... + 4*(L3*M1+L1*M3)*((S35-S25)*L3*M3+(S15-S25)*L1*M1) ... + S44*(L2*M3-L3*M2)^2 ... + S55*(L3*M1-L1*M3)^2 ... + S66*(L1*M2-L2*M1)^2 ... +2*S45*(L2*M3+L3*M2)*(L3*M1+L1*M3) ... +2*S56*(L3*M1+L1*M3)*(L1*M2+L2*M1) ... +2*S46*(L1*M2+L2*M1)*(L2*M3+L3*M2); G1=collect(expand(G1), ... [S11 S12 S13 S14 S15 S16 S22 S23 S24 S25 S26 S33 S34 S35 S36 S44 S45 S46 S55 S56 S66]); G1 G2=4*( S11*(L1*M1)^2 + S22*(L2*M2)^2 + S33*(L3*M3)^2 ) ... + S44*(L2*M3+L3*M2)^2 + S55*(L1*M3+L3*M1)^2 + S66*(L1*M2+L2*M1)^2 ... +8*( S12*(L1*L2*M1*M2) + S13*(L1*L3*M1*M3) + S23*(L2*L3*M2*M3) ) ... +4*(L1*M1)*( S14*(L2*M3+L3*M2) + S15*(L1*M3+L3*M1) + S16*(L1*M2+L2*M1) ) ... +4*(L2*M2)*( S24*(L2*M3+L3*M2) + S25*(L1*M3+L3*M1) + S26*(L1*M2+L2*M1) ) ... +4*(L3*M3)*( S34*(L2*M3+L3*M2) + S35*(L1*M3+L3*M1) + S36*(L1*M2+L2*M1) ) ... +2* S45*(L2*M3+L3*M2)*(L1*M3+L3*M1) ... +2* S46*(L2*M3+L3*M2)*(L1*M2+L2*M1) ... +2* S56*(L1*M3+L3*M1)*(L1*M2+L2*M1); G2=collect(expand(G2), ... [S11 S12 S13 S14 S15 S16 S22 S23 S24 S25 S26 S33 S34 S35 S36 S44 S45 S46 S55 S56 S66]); G2 dG=collect(expand(G2-G1), ... [S11 S12 S13 S14 S15 S16 S22 S23 S24 S25 S26 S33 S34 S35 S36 S44 S45 S46 S55 S56 S66]); % dG=subs(dG, ... % [S11 S12 S13 S14 S15 S16 S22 S23 S24 S25 S26 S33 S34 S35 S36 S44 S45 S46 S55 S56 S66], ... % [0*S11 0*S12 0*S13 0*S14 0*S15 0*S16 0*S22 0*S23 0*S24 0*S25 0*S26 0*S33 0*S34 0*S35 0*S36 0*S44 0*S45 0*S46 0*S55 0*S56 0*S66]); dG=subs(simplify(dG), [L1*M1+L2*M2+L3*M3], [0]) 常值讨论

评论中有人指出, 在某些方向上剪切模量为常值, 这是可能的, 我给了一个简单的说明. 详细的讨论请参考下面两篇文章:

T. C. T. Ting; J Elasticity 81(3):271-292, 2006; 10.1007/s10659-005-9016-2 Kevin M. Knowles, Philip R. Howie; J Elast 120(1):87-108, 2014; 10.1007/s10659-014-9506-1 在线工具

本来我还打算做一个在线工具, 但是发现早就有人做好了, 所以还是作罢吧.

如果你需要对比, 可以使用ELATE: Elastic tensor analysis

评论 2015-04-01 16:52:10 棒棒真是的太棒了

2015-04-01 20:06:30 Jerkwin 希望能提供更多信息.

2015-11-16 21:01:02 李润光程序很棒，但为何（111）剪切面上不同方向剪切模量是定值？可否给出推导过程或出处？非常感谢！ 2015-11-16 22:38:25 Jerkwin 不是很清楚你的问题, 请说得具体些. 2015-11-17 09:33:42 李润光比如Fe是面心立方结构，最易开动的滑移系应该是（111）,但根据你的公式计算（111）面上沿任意方向滑移其剪切模量均是同一定值，从最大最小剪切模量图也可以看出，（111）面和（100）面上沿不同方向滑移的剪切模量都是定值，我认为这是不合理的，但没有找到相关的推导，可否请您解释一下。 2015-11-17 11:27:41 Jerkwin 如果我没理解错的话, 你把图的含义搞混了. 那些图给出的是某一方向上剪切模量的极值, 而不是同chi对应的剪切模量. 如果那样的话, 对每一方向, 都可以做出剪切模量和chi的关系图, 可能那是你想要的. 2015-11-17 12:17:40 李润光我懂您的意思，图的含义我理解，是每个晶面上的模量极值。但仔细看一下，您的每一对图中，同种材料（111）面和（100）面上剪切模量的极大值和极小值均是相同的，我也用您的子程序画了这两个面上剪切模量和chi的关系图，是一条水平线，这是我困惑的地方，晶体学上，同一晶面沿不同剪切方向的模量不会是定值，我查找了相关文献，极少看到有关于此方面的推导，非常希望您给予进一步解释。 2015-11-17 13:02:20 Jerkwin 如果你说的(111)面对应的theta和phi都是PI/4的话, 执行主函数里面的[Gmin, Gmax] = pltGchi(.25pi, .25pi)语句会给出G和chi的关系图, 我得到的结果是正弦曲线, 不是直线, 和你说的不同. 2015-11-17 13:13:40 李润光（111）面theta并非PI/4。按以下语句执行， Plane = [1 1 1]; pltGchi ( acos( Plane(3)./norm(Plane)), atan( Plane(2)./Plane(1)) ); 输出图形即为（111）面对不同chi的shear modulus，是一条直线。 Plane = [1 0 0]时也是直线。请您验证一下。 2015-11-17 14:40:38 Jerkwin 你说的对. 对这个特殊的方向, 剪切模量确实是常数. 原因如下: 根据最后的公式, 当G为常数时, 只需要满足Gmax=Gmin即可, 这样就可以得到C1=0并且C2=C3. 这两个等式也可以直接从G的第一个表达式看出来. 根据这两个条件就可以算出那些方向上的G为常数. 我粗略地算了一下, C1=-(alpha-beta) sin^2(tht) cos(tht) sin(4phi)/4 C2-C3=(alpha-beta) sin^2(tht) { sin^2(tht) [ sin^2(2phi)-4 ] + cos(4phi)+3 } 所以, 当alpha=beta时, G恒为常量. 除此之外, 在特定方向上G也可以为常量.其中的一个解就是你所说的(111)方向, 此时phi=PI/4, C1=0, cos^2(tht)=1/3, C2-C3=0, 满足等式. 当然, 根据这两个等式, 你还可以算出其他的方向, 我就没有验证了. 你有兴趣的话可以验算一下, 告诉我结果, 以便更新下上面的博文. 2015-11-18 17:29:36 李润光恩，我想了下，（111）面模量为定值也是合理的。通常所说的最密排面最易滑移方向指实质上是指临界分切应力较小，势垒比较低。方向最易滑移并非因为模量小，而是原子间距小，平均伯氏矢量小，势垒低，这样看来你的计算在晶体学也是合理的。还是想看下原始推导过程，可能更有启发。 2015-11-18 23:45:34 Jerkwin 我不是研究晶体的, 所以晶体方面的知识不足. 希望你能把这个问题及其含义从晶体学方面论述一下, 我把数学部分补充上, 这样才好整理发布. 你写好后可以发到我的邮箱: [email protected]. 2016-03-15 15:02:34 敦煌的木桶下面这篇文章有推导，另外，对立方晶系而言，确实在(100)和(111)面上沿任意方向的剪切模量的值都是一样的。其他晶系的结果就不是这样的。 J Elast (2015) 120:87–108 Kevin M. Knowles · Philip R. Howie The Directional Dependence of Elastic Stiffness and Compliance Shear Coefficients and Shear Moduli in Cubic Materials

2016-03-17 06:02:36 Jerkwin 谢谢提供信息. 等我看过文章之后更新下博文.

2016-07-17 07:06:45 767676 必须点赞！良心好贴！

2016-07-25 16:04:52 Jerkwin 谢谢

2016-08-15 22:42:29 767676 本人在您的程序基础上，尝试做某切面轮廓，加代码如下： function SphPlt(X, Y, Z, F, File, Title) surf(X,Y,Z,F, ‘FaceColor’,’interp’, ‘EdgeColor’,’none’); hold on % 作模量的某一切面图 [X,Y,Z]=meshgrid(linspace(-Gmax,Gmax)); contourslice(X,Y,Z,X,x,y,z,[0 0]) axis tight; title(Title); view(30,30); daspect([1 1 1]); 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 但是没有成功，敬请版主点播！ 2016-08-16 23:51:27 Jerkwin 你要理解代码的含义. 你照抄的这几句放在这里不行, 因为大写的XYZ和小写的XYZ弄反了, 而且没有Gmax. 正确的差不多是这样 hold on % 作模量的某一切面图 [x,y,z]=meshgrid(linspace(-max(F(:)),max(F(:)))); contourslice(x,y,z,x,X,Y,Z,[0 0]) 2016-08-17 21:14:02 767676 多谢大神点拨！正如您所言，是我没有吃透代码含义，太心急了。敢问大神，能否给出您关于此段程序的引用文献。倘若有所小作，绝不敢忘引用您的原创。

2016-08-18 22:23:24 Jerkwin 如上所言, 我不是做材料的, 这些东西只不过是我业余帮人写的, 没有必要, 也不可能引用我的文献. 如果你真要引用相关文献的话, 引用上面ELAM的文献吧, 有些公式是我从它的源代码中抽出来的.

2016-08-17 17:40:47 斌楼主可否提供立方外其它晶系的剪切模量的计算公式

2016-08-18 22:20:13 Jerkwin 公式我有, 但太复杂了, 还没有测试正确性. 由于我本人不是做材料的, 所以也没有这种需求. 如果你确实需要计算公式的话, 请参考我的另一篇博文那个程序的代码中有你需要的公式.

2016-10-27 12:37:45 好友昵称么 “此外, 与杨氏模量不同的是，剪切模量在某一方向 l 上的值还与角度 χ 有关. 这就导致剪切模量的各向异性比杨氏模量更加复杂。”我想请问下，上面这句话的角度 χ是怎么定义的

2016-10-27 21:31:06 Jerkwin 上面图中标出了χ 角, 仔细看下.

2017-03-12 21:49:50 流星楼主其他晶系的剪切模量公式在哪有啊，能不能提供链接？

2017-03-13 11:07:58 Jerkwin 上面有人问了, 见我的回复.

2017-03-14 00:11:49 流星请问下楼主，求不同晶面的比如（111）的剪切模量，要改代码的哪个部分啊？ 2017-03-14 08:07:53 Jerkwin 不需要改什么, 使用Gchi函数就可以了 2017-03-14 19:10:52 流星楼主大大，刚学matlab,对于脚本看的不是很懂，需要画个（111）面的切变模量，你说的用Gchi函数具体怎么操作啊？要是有时间，请指导下，谢谢！ 2017-03-15 07:27:27 Jerkwin 给了theta和phi之后, 函数pltGchi(theta, phi)可以计算G和chi的关系, 并作图. 你试下就知道了. 2017-03-15 08:13:33 流星对于（111）面的剪切模量，是不是二维的平面图啊？那知道theta和phi，在哪里调用pltGchi(theta,phi)函数啊，是在最后的function SphPlt（）里面吗，那原来的 surf(X,Y,Z,F, ‘FaceColor’,’interp’, ‘EdgeColor’,’none’);是不是用#注释掉啊？哈哈，问题有点多，脑袋不灵光，希望楼主能回答。嗯，楼主大才，赞！！ 2017-03-15 09:44:26 Jerkwin 基本语法你都不明白的话, 我说了也是白白浪费时间. 你先去看看matlab的基础再说. 2017-03-15 12:09:57 流星好的，还是多谢大神啦

【本文地址】

公司简介

联系我们