Differential equation

[MartonTimea]

Hi. I would like to find z from the differential equation below as a function of time. Can anyone help me with this? All the other terms are constant. I would really appreciate some help as I am new to mathcad. Thank you in advance.

[AlanStevens]

Here's an example. The attached is in MC15 format. In case you aren't using MC15 I've inserted an image below (if you're using Prime you don't need Given, you put the ode inside a solve block). Obviously, you will need to replace the data with your own!

Alan

[MartonTimea]

Hi thank you for the quick answer. This is very helpful. However i need to compute z for a time interval of 1680 seconds. And if i do this, the graphic does not compute anything. Could you tell em what is wrong?

I have attached the file as well.

[Fred_Kohlhepp]

Don't know why Odesolve is crashing; but the older versions still work.

[AlanStevens]

You forgot the power of 5/3 on the derivative term Fred. Put this on and rkfixed doesn't like it!

The problem is that z goes negative fairly quickly (less than 20 seconds) with the values and ODE specified. This would make z^(5/3) return a complex number and Mathcad's ODE solvers don't know what to do with that.

The data and equations need to be looked at more carefully (perhaps it's just that the initial value is far too small, perhaps the z term should be the absolute value of z, perhaps ...). Without knowing where the equation and data come from I can't be more specific.

Alan

[AlanStevens]

Ok, I think the main problem is that cp is 1 kJ/kg; but since you aren't actually using units you need to make sure everything is consistent with standard SI. This means here that you should have cp = 1000 (J/kg). (The intial starting height in the workfile is 1 metre.) With cp as 1000 the file works up to 1680 ok.

Alan

[ValeryOchkov]

Cp=J/kg/K

[Werner_E]

The reason for the odesolve block failing is the expression z(t)^(5/3). When z(t) gets neagtive, that expressions yields a complex result and odesolve fails.

Mathcad will give complex results for a^(1/3) if a is negative. While often undesired this is mathenatically correct, Mathcad yields the solution of x^3=a with the smallest argument(phase).

You could use sign(a)*|a|^(1/3) but Mathcad also provides as a compromise the root operator (Calculator palette) which defaults to real results if used on real arguments (mathematically not consequent if seen in a complex context, but often very convenient).

So change the expression in your equation and you get what you demanded. but I guess Alan is right to set cp:=1000 as you probably don't expect negative z values.

[MartonTimea]

Q represents the heat release and is given in kW, therfore cp should be in kJ. If i use cp=1000 i get values that are very high, which means that i have a smoke filling in a compartment of hundreds of meters. Which practically seems impossible. I have attached a small pdf if somebody is curious to see where the equation comes from. I assume the equation only works for a weak heat release rate. That is a low Q.

Thank you all for trying to help. I really appreciate it.

[AlanStevens]

Marton Timea wrote:

Q represents the heat release and is given in kW, therfore cp should be in kJ.

In that case you have over 4MW of heat going in to your compartment. Perhaps this leads to a rapid pressure rise such that the approximation of dropping the first term in equation 8.30 is no longer valid.

Alan

[Fred_Kohlhepp]

Solution per pdf: