cancel
Showing results for 
Search instead for 
Did you mean: 
cancel
Showing results for 
Search instead for 
Did you mean: 

Help to solve a non-linear ODE in mathCAD 15

Pyros1111
5-Regular Member

Help to solve a non-linear ODE in mathCAD 15

Hi,

 

I am attempting to dynamically model the behaviour of granular solids in a rotary kiln using known equations, however i am having trouble setting up the following ODEs in mathCAD for the model....can anyone help. I have included the key points i am trying to model in Mathcad.

 

 

 

bed height equation.PNG

boundary condition: boundary condition.PNG

where h0 is a fixed value (say 10% of the value of R).

 

From the above i'm looking to be able to determine the bed height (h) at any axial location (z) along the kiln tube.

 

Subsequently equations for Cross section area of the granular bed and residence time of the granular transport through the kiln (as a function of bed depth (h)) follow on from this...

 

Cross section area:

 

Cross section Area equation.PNG

Residence time

residence time equation.PNG

So essentially i am looking to use MatCAD to find the residence time (t), given fixed values for: n (rotation speed in rpm), R (radius of kiln in m), D (diameter of kiln in m), angle of repose of the material (beta), angle of inclination of the kiln (theta), mass flow rate (ms) , bulk density of material (rho), axial length of the kiln tube (z).

 

My initial attempt using ODEsolver block did not get that far and with limited success (see attached) and seems very temperamental depending on what values are used and cannot utilise units for the variables (hence why they are all unitless.

 

 

 

29 REPLIES 29
-MFra-
21-Topaz I
(To:Pyros1111)

Hi,

1111.jpg

Pyros1111
5-Regular Member
(To:-MFra-)

Thanks very much !!

Pyros1111
5-Regular Member
(To:Pyros1111)

is there anyway to do it using the specific units for the different variables?

Units in ODE's in Mathcad 15 pose a significant challenge: take them out incorrectly and the problem is not properly scaled.  One solution I found was to use the function "UnitsOf". 

If you define a variable with a unit (L := 2 ft), to use it in an ODE it must have the units stripped.  Clearly, dividing by inches, feet, or meters will create different values.  Mathcad has a consistent unit system (when I type L = my version 15 would have defaulted to meters,) to get all the scaling correct, divide each expression by the function.  So L becomes "L/UnitsOf(ft)".  Note that "L/UnitsOf(L) would work as well EXCEPT if you have a unit mismatch the first expression would catch it, the second one would not.  (If L had dimensions of velocity the first expression (UnitsOf(ft) would crash, the second would happily just strip the wrong units.

 

But ODE's in Prime can use units. . .

Pyros1111
5-Regular Member
(To:Fred_Kohlhepp)

Hmmm interesting... trying your example out for myself i get this:

 

Capture x.PNG

is this correct? it still looks like it has defaulted to metres i.e 2ft = 0.61m?

is this correct? it still looks like it has defaulted to metres i.e 2ft = 0.61m?

 

You are correct-- 2 ft has become 0.61, but there are no attached units.  So your equation (ODE) will be rendered unitless but with proper scaling.

 

The unitless values are dependent on the unit system chosen.

With unit system "U.S."you get:

LucMeekes_0-1578431431077.png

But with unit system set to "S.I." you get:

LucMeekes_1-1578431514978.png

Success!
Luc

Thanks Luc!

 

The point I was trying to make is that to solve an ODE correctly it must be properly scaled; using the UnitsOf function ensures that Mathcad will scale the equation correctly into whichever default system is selected.  The solution can then be "reunitized" by multiplying by UnitsOf(the correct type of dimension).

Pyros1111
5-Regular Member
(To:Fred_Kohlhepp)

applying what your units workaround to my current Mathcad 15 sheet i cant get it to run - see attached? what have i missed?

There are a couple things I've noticed, tried hard to address; see comments pdf attached

LucMeekes
23-Emerald I
(To:Pyros1111)

Your rotational speed may be a litle low. 1/3 of a rotation per minute apparently isn't enough, 333 rotations per minute does give a solution.

Mind that you don't have to use the UnitsOf function to cancel units, and 'rad' is a unitless 'unit' anyway. You can also divide the units out directly.

See below.

LM_20200107_Odesolve1.png

Then, with Dw reduced to 0.5 m,and the minus sign in front of the derivative removed, I get:

LM_20200107_Odesolve2.png

Success!
Luc

Pyros1111
5-Regular Member
(To:-MFra-)

thanks. - will this work using specific units with the variables?

LucMeekes
23-Emerald I
(To:Pyros1111)

In Mathcad 15 and before you cannot have units with an Odesolve (unless they cancel out).

One of the advantages of Prime is that you can.

 

Success!
Luc

Pyros1111
5-Regular Member
(To:LucMeekes)

tried putting the same thing into Mathcad prime (so i could use units as suggested) but i just get various errors... clearly its done a bit differently in prime and i am missing something - any help?

Pyros1111
5-Regular Member
(To:LucMeekes)

tried putting the same thing into Mathcad prime (so i could use units as suggested) but i just get various errors... clearly its done a bit differently in prime and i am missing something - any help?

Prime is different than 15.

 

Did you translate the 15 file, or start a new file?

 

Can you post your Prime sheet?  (Usually gets prompter and more focused help.

Pyros1111
5-Regular Member
(To:Fred_Kohlhepp)

tried to start a new one... see attached.

LucMeekes
23-Emerald I
(To:Pyros1111)

You don't need guess values for an odesolve. Put all variable definitions above, out of the guess part of the solve block.

 

Success!
Luc

Concerning the odesolve block it should make no difference if you define your constants inside the solve block or above. But if you define them inside, the values are not available below and outside of the solve block. So it may be a good idea anyway to define those constants above the block.

The reason for your solve block failing is the value of your variable R.w.

Compared to your Mathcad sheet you used different values for m.s and R.w in your Prime attempt.

Especially the low value 0.25 for R.w compared to the value 2 you had in Mathcad seems to make the solve block fail. It looks like (with the other values unchanged) R.w must be at least 0.502 or larger.

Werner_E_0-1577134407264.png

 

Pyros1111
5-Regular Member
(To:Werner_E)

Hi,

 

when i tried to put the relevant units in (ant altered the values accordingly) it didn't want to work - various warnings (see attached). Any pointers?

  1. The initial conditions must be given with the appropriate units. Its not necessary to write 0m instead of just 0, but its good habit to do so, too
  2. The solve block fails if R is lower than 1.97 m. You'll have to check your equation or the other values assigned if you think that R=0,5m should yield a valid result. I notice that you changed several values compared to the file I posted (where R must had a value of at least 0.502).
  3. To plot with units on the abscissa you need to define a range variable. I have chosen z.plot as the name for this variable which is only to be used for plotting, but you can chose any name you like (even simply z :-).

Werner_E_0-1578322799185.png

 

Pyros1111
5-Regular Member
(To:Werner_E)

Hi,

 

thanks for this - i see how that works now.

 

I was also trying to alter the initial condition for Hz as its not correct and should be Hz(L)=0.01 rather than Hz(0)=0.01 however this doesn't seem to work? does something else need to be changed?


@Pyros1111 wrote:

Hi,

 

thanks for this - i see how that works now.

 

I was also trying to alter the initial condition for Hz as its not correct and should be Hz(L)=0.01 rather than Hz(0)=0.01 however this doesn't seem to work? does something else need to be changed?


Alan already solved that. Its as simple as writing the desired inital condition in the solve block and set the "end" value to zero. No R can be a higher value, too.

Werner_E_0-1578342963850.png

 

LucMeekes
23-Emerald I
(To:Pyros1111)

One other thing:

Theta and beta are angles I guess. You provided a large value for Theta that you did not intend: 40, that is 40 radians, not degrees!

Sad to say that when I change to 40 degrees and 3 degrees respectively, no solution is found...

Success!

Luc

Pyros1111
5-Regular Member
(To:-MFra-)

sorry to be a bother again but when i try and put this in it still doesn't seem to work properly for the last bit - says variable undefined??? any pointers on what im doing wrong?

 

Capture.PNG

Pyros1111
5-Regular Member
(To:Pyros1111)

my fault - was using a Boolean equals sign.

-MFra-
21-Topaz I
(To:Pyros1111)

If constants and variables are dimensioned, they must be rendered dimensionless before the Odesolve block. After it, you can redimension them again.

Pyros1111
5-Regular Member
(To:-MFra-)

hi again,

 

thanks again for your previous reply - i have been looking over it and trying to use the more accurate values to see what comes out and noticed i had initially described the inmitial condition for the ODEsolver as Hx(0)=0.01, whereas what i wanted to do was define it as Hx(L) = 0.01. However, when i do this in the Mathcad (See attached) it doesn't work and says the endpoints cannot be the same in the solve block?

Capture 4.PNG

can you explain why and how to make it work for this end constraint please?

Try integrating in reverse from Lretort to zero:

 

reverse.JPG 

See attached.

Alan

 

Announcements
Check out the latest
Mathcad Tip
"PTC Mathcad 15 / Prime 1-6 Update."