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

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

[Torsten]

If you don't give us code to reproduce the error, how should we find the problem ?

[Torsten]

Thank you for the files.

I guess there are test case for the code. Does it work for these test cases ? Or is it a test case you try to make work ?

[Bamelari Jovani]

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]

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.

Answer 1

[Divyajyoti Nayak]

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