Obtaining wrong answer using solve(...)
조회 수: 1 (최근 30일)
이전 댓글 표시
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
댓글 수: 0
답변 (1개)
Junfei Li
2022년 2월 7일
%% 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
댓글 수: 0
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!