Obtaining wrong answer using solve(...)

1 view (last 30 days)
Kenneth Andersen
Kenneth Andersen on 6 Feb 2019
Answered: Junfei Li on 7 Feb 2022
I am trying to use the equation solver: solve(). The goal is to obtain the bulk and shear compressibility of a composite medium using a method derived by Devaney. I am confident (now, with some doubt) that the code worked before I installed the latest Matlab version 2018b; however, it is almost two years since I wrote it. The code is shown below. The obvious problem is that solve() returns unphysical, negative numbers.
%% 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
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 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);
Initially, I used solve this way
S = solve([K0 + C(ii)*x*dK == Ke, G0 + C(ii)*y*dG == Ge]);
and the shear (Ge) compressibility then seems then to be correct.
What am I doing wrong? Is there another way to "solve" this problem?
In advance, thanks for any help

Answers (1)

Junfei Li
Junfei Li on 7 Feb 2022
%% 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

Tags

Products


Release

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!