求顺度系数,需注意应变与剪切应变的2倍关系
clear;
syms c11 c12 c44;
c11 = 10;c12 = 60;c44 = 20;
epsilon_11 = 0.1;epsilon_22=0.2;shear_strain_12 = 0.3;epsilon_12 = 1/2 * shear_strain_12;
e = [epsilon_11 epsilon_22 0 shear_strain_12 0 0]';%此处e(6)=0.3表示剪切应变xy,应变xy=1/2*剪切应变=0.15
e11 = epsilon_11; e22 =epsilon_22; e12 = epsilon_12;e33=0;%plant strain problem
Cijkl = [c11 c12 c12 0 0 0;
c12 c11 c12 0 0 0;
c12 c12 c11 0 0 0;
0 0 0 c44 0 0;
0 0 0 0 c44 0;
0 0 0 0 0 c44];
Cijkl_inv = inv(Cijkl);
s_origin = Cijkl * e; %此处s12表示剪切应力
e_origin = Cijkl \ s_origin; % Cijkl^-1 * sigma
s11 = c11 * e11 + c12 * e22;
s22 = c12 * e11 + c11 * e22;
s33 = c12 * e11 + c12 * e22;
s12 = 2 * c44 * e12;
e11_byhand = (c11 + c12)/(c11^2 + c11*c12 - 2*c12^2) * s11 -c12/(c11^2 + c11*c12 - 2*c12^2) * (s22+s33);
e22_byhand = (c11 + c12)/(c11^2 + c11*c12 - 2*c12^2) * s22 -c12/(c11^2 + c11*c12 - 2*c12^2) * (s11+s33);
e33_byhand = -c12/(c11^2 + c11*c12 - 2*c12^2) * (s11+s22) + (c11 + c12)/(c11^2 + c11*c12 - 2*c12^2) * s33;
e12_byhand = 1.0 / c44 * s12;