%% Solve Eq.(8a) and (8b) symbolically of Devaney's:
% Effective elastic parameters of random composites
% K: bulk compressibility
% G: shear compressibility
% C: volume fraction
% Subscript 0 refers to matric medium, i.e. epoxy
% Subscript e refers to effective medium
% no subscript refers to the particles being added, i.e. aluminium oxide
clear; clc;
syms Ke Ge
C = 0:0.1:1;
% Epoxy material parameters
K0 = 5.1655e09; % Bulk compressibility of epoxy matrix
G0 = 1.5464e09; % Shear compressibility of epoxy matrix
rho0 = 1094; % Density of epoxy matrix
% Al2O3 material parameters
K = 228e9; % Bulk compressibility of Al2O3
G = 152e9; % Shear compressibility of Al2O3
rho = 3950; % Density of Al2O3
%% Calculate
dK = K - K0;
dG = G - G0;
drho = rho - rho0;
x = (3*Ke + 4*Ge)./(3*Ke + 4*Ge + 3*dK);
y = ( 5*(3*Ke + 4*Ge)*Ge )./( (15*Ke + 20*Ge)*Ge + 6*(Ke + 2*Ge)*dG );
Kv = zeros(1, length(C));
Gv = zeros(1, length(C));
rhov = zeros(1, length(C));
for ii = 1:length(C)
S = solve([K0 + C(ii)*x*dK == Ke, G0 + C(ii)*y*dG == Ge, Ke>0, Ge>0], [Ke Ge]);
Kv(ii) = vpa(S.Ke);
Gv(ii) = vpa(S.Ge);
clear S
rhov(ii) = rho0 + C(ii)*drho;
end
% Phase and shear velocities Eqs.(9a,b)
Vc = sqrt( (Kv + 4/3*Gv)./rhov );
Vs = sqrt(Gv./rhov);
Z=rhov.*Vc;
figure;
subplot(3,1,1)
plot(C,Vc)
ylabel('Velocity')
subplot(3,1,2)
plot(C,rhov)
ylabel('Density')
subplot(3,1,3)
plot(C,Z)
ylabel('Impedance')
Hi Kenneth,
You may not be interested in this question anymore, but I am posting my modification here. (Your code helped me hope it can help others in the future!) I added two ineualities into the equations and got a seemingly reasonable answer.
S = solve([K0 + C(ii)*x*dK == Ke, G0 + C(ii)*y*dG == Ge, Ke>0, Ge>0], [Ke Ge]);
I have checked the plot with previous literature, and the curve look indeed similar. Hope this is still useful to you...
Best,
Junfei
