Clear Filters
Clear Filters

Nozzle Design Error: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

10 views (last 30 days)
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
Bamelari Jovani
Bamelari Jovani on 14 Aug 2024
It is a test case and I am trying to make it work. However, I encounter the error as described above. I still believe the error lies in the internalnode.m function. Maybe in line 92 where X = A\B
Torsten
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.

Sign in to comment.

Accepted Answer

Divyajyoti Nayak
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!!

More Answers (0)

Products


Release

R2024a

Community Treasure Hunt

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

Start Hunting!