Nozzle Design Error: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
10 views (last 30 days)
Show older comments
Bamelari Jovani
on 13 Aug 2024
Answered: Divyajyoti Nayak
on 21 Aug 2024
Greetings everyone, I am trying to look into a code found here corydodson/NozzleDesign: Design of supersonic nozzles (github.com). However, I run into an error whenever I run the code which reads "Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. " I belive the error comes from this internalnode.m function. x,y are the coordinates, t0 is the temperature and d can be 0 or 1 depending on the user. Anyway I can resolve this? Any help is appreciated!.
I have attached the relevant codes below which can also be found here: corydodson/NozzleDesign: Design of supersonic nozzles (github.com).
function [xOut,yOut,uOut,vOut] = internalnode(t0,species,moleFracs,x,y,u,v,d)
h0 = mixprop('h',species,moleFracs,t0);
r = mixprop('r',species,moleFracs)*1000;
L = length(x);
minus = logical([ones(L - 1,1);0]);
plus = logical([0;ones(L - 1,1)]);
xm = x(minus);
xmCalc = xm;
xp = x(plus);
ym = y(minus);
ymCalc = ym;
yp = y(plus);
ypCalc = yp;
um = u(minus);
umCalc = um;
up = u(plus);
upCalc = up;
vm = v(minus);
vmCalc = vm;
vp = v(plus);
vpCalc = vp;
sp = vm;
N = 20;
err = 1e-4;
notConv = true(1,4);
for i = 1:N
ymCalc = (ym + ymCalc)/2;
ypCalc = (yp + ypCalc)/2;
umCalc = (um + umCalc)/2;
upCalc = (up + upCalc)/2;
vmCalc = (vm + vmCalc)/2;
vpCalc = (vp + vpCalc)/2;
uCalc = [umCalc;upCalc(end)];
vCalc = [vmCalc;vpCalc(end)];
vMag = sqrt(uCalc.^2 + vCalc.^2);
h = h0 - vMag.^2/2000;
t = tempfromprop(species,moleFracs,'h',h);
g = mixprop('gamma',species,moleFracs,t);
a = sqrt(r*g.*t);
am = a(minus);
ap = a(plus);
mu = asind(a./vMag);
mum = mu(minus);
mup = mu(plus);
theta = atand(vCalc./uCalc);
thetam = theta(minus);
thetap = theta(plus);
lm = tand(thetam - mum);
lp = tand(thetap + mup);
q = uCalc.^2 - a.^2;
qm = q(minus);
qp = q(plus);
rm = 2*umCalc.*vmCalc - qm.*lm;
rp = 2*upCalc.*vpCalc - qp.*lp;
sm = d*am.^2.*vmCalc./ymCalc;
switch isempty(sp(ypCalc == 0))
case 1
sp = d*ap.^2.*vpCalc./ypCalc;
otherwise
sp(ypCalc ~= 0) = d*ap(ypCalc ~= 0).^2.*vpCalc(ypCalc ~= 0)./ypCalc(ypCalc ~= 0);
sp(ypCalc == 0) = sm(end);
end
A = zeros(L);
B = zeros(L,1);
for j = 1:L - 1
A(2*(j - 1) + 1:2*j,2*(j - 1) + 1:2*j) = [1,-lp(j);...
1,-lm(j)];
B(2*(j - 1) + 1:2*j) = [yp(j) - lp(j)*xp(j);...
ym(j) - lm(j)*xm(j)];
end
X = A\B;
ymCalc = X(1:2:end);
ypCalc = ymCalc;
xmCalc = X(2:2:end);
xpCalc = xmCalc;
for j = 1:L - 1
A(2*(j - 1) + 1:2*j,2*(j - 1) + 1:2*j) = [qp(j),rp(j);...
qm(j),rm(j)];
B(2*(j - 1) + 1:2*j) = [sp(j)*(xpCalc(j) - xp(j)) + qp(j)*up(j) + rp(j)*vp(j);...
sm(j)*(xmCalc(j) - xm(j)) + qm(j)*um(j) + rm(j)*vm(j)];
end
X = A\B;
umCalc = X(1:2:end);
upCalc = umCalc;
vmCalc = X(2:2:end);
vpCalc = vmCalc;
switch i ~= 1
case 1
notConv = abs([xmCalc,ymCalc,umCalc,vmCalc]./check0 - 1) > err;
end
check0 = [xmCalc,ymCalc,umCalc,vmCalc];
switch sum(sum(notConv)) == 0
case 1
break
end
switch isnan(xmCalc) | isnan(ymCalc) | isnan(umCalc) | isnan(vmCalc)
case 1
break
end
end
xOut = xmCalc;
yOut = ymCalc;
uOut = umCalc;
vOut = vmCalc;
end
4 Comments
Torsten
on 14 Aug 2024
Edited: Torsten
on 14 Aug 2024
If it is a test case supplied by the authors of the package, you should contact the authors for help.
It won't help if you find that the error message stems from line 92 where X = A\B is computed. You must be able to deduce how A and B are being built from your inputs and what finally makes the command X = A\B fail. For this, it will be necessary to understand the complete code and its algorithm.
Accepted Answer
Divyajyoti Nayak
on 21 Aug 2024
Hi,
The reason you are getting this warning is because some matrices used in the code contain NaN values. On debugging the code, I found that these lines in file ‘ivcurvekliegl.m’ are the causing issue.
% z(end) = interp1(r,z,0,'spline');
% u(end) = griddata(r,z,u,0,z(end));
% v(end) = griddata(r,z,v,0,z(end));
% t(end) = griddata(r,z,t,0,z(end));
The ‘griddata’ function tries to interpolate a 3d surface at the required query points but if it fails it returns NaN. On commenting out these lines the warnings disappear. Although more errors do pop up if you run the code for rectangular cross section (‘d’ = 0). The code runs well for axisymmetric cross section (‘d’ = 1) as given in the readme file of the repository.
In the case of rectangular cross section, the difference between exit mach number and design mach number doesn’t fall under tolerance, so the ‘maxTurnAngle’ does not get calculated. Instead ‘theta’ remains a vector which causes error. You can still remove the error by hard coding ‘theta’ to be a scalar without checking the tolerance. I chose the last value of the ‘theta’ vector, hardcoded it and was able to run the code with no problems.
switch abs(mExit(i) - mDesign) < err
case 1
%This line is not reached when d = 0
theta = theta(i);
break
end
end
if(d == 0)
theta = theta(end); %Hardcoded theta for d = 0
end
[x,y,u,v,nozzleCd] = kernel(t0,p0,species,moleFrac,rTd,yt,d,N,theta,1);
Results are in the attached image. Hope this helps!!
0 Comments
More Answers (0)
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!