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