ODE's error with if else function

[MC_14395780]

I am solving a system of 1D steady ODEs in Mathcad using RK4. The system works correctly until I add a heat transfer term representing two-way coupling between a gas phase and injected droplets. When this term is included, the solver no longer outputs results and reports that the gradient matrix must be scalar or a matrix.

The gas stagnation temperature gradient depends on Reynolds number, Nusselt number and a convective heat transfer coefficient that is itself a function of the ODE variables, which vary along z (length of nozzle)

All other gradients (Mach, pressure, velocity, etc.) are functions of the same state vector and work correctly when the heat term is removed. I suspect I am incorrectly defining a function within the gradient matrix or using If else wrong.

Any help would be appreciated, as I'm going round in circles.

below is screenshot of part of code:

Thanks very much

[Werner_E]

How about posting the worksheet itself instead of just pictures? And clarifying exactly which added term made the system fail.

What I noticed is that in your definitions of some grad... functions you divide by x7(y1) and in the vector x you provide x7 is zero.

Furthermore it looks like x8(T01) never is used, but I may have overlooked it in the pictures.

[MC_14395780]

Hi, 

I have now added the file above. I have removed x8 (T01) from the vector, and changed y1 (x7) to 0.1.

I think it may be something to do with my if else statement but I can't spot the mistake.

[Werner_E]

You can use Prime's error trace facility to trace back to the position the error shows up first.

In your case it traces back to the call of function F.drag.

You defined F.drag as a function with two arguments

but you call it with just one argument

The error message tells you exactly that 😉

Same problem with function C.d

Unfortunately fixing this

don't help for rkfixed.

It looks like the results do not depend much on z. The results are the same for z=0 and z=L.e=0.5.
Are the values expected to that high? (10^11)

Which was omitted when the system finally worked OK? What was added causing the problem?

[Werner_E]

I replaced gradT.dr and grady by

and I still get the same error in rkfixed ('this function cannot be used here').

So it looks that the problem are not the if..else..

Is it on purpose that these function are independent from z?