哲●科●文Jerkwin 类别 标签 链接 关于 订阅

剪切模量各向异性

类别: 科   标签: matlab 编程   阅读次数:   版权: (CC) BY-NC-SA ✕

以前的博文中讲过用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
211
function 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
19
clc;

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
46
clc; 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])

常值讨论

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

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

在线工具

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

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

评论

◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆
前一篇:　GROMACS表面张力单位,计算及其长程校正
后一篇:　Taggie：基于标签的桌面搜索工具
访问人次(2015年7月 9日起): | 最后更新: 2024-11-01 02:53:58 UTC | 版权所有 © 2008 - 2024 Jerkwin 个