Simulation ODE with time functions

[AD_10780197]

Good morning, everyone,

I'm currently working on a Mathcad sheet. I have two time-dependent functions that I want to combine with a differential equation. However, I keep getting an error message, and I'm not sure what I'm doing wrong.

I managed to get it working in MATLAB, but I would like to finish this properly in Mathcad. The issue occurs on the last page, where the solve block isn't isolated from the rest of the calculations 

Does anyone know how I can solve this? Thanks in advance!

[LucMeekes]

1st. You should state that you are using Prime 5, to prevent that you get an answer in prime 6...11 which you cannot open.

2nd. See what happens if you remove the " x''(t)= " from the "else" part in your constraint.

Phrased as it is now, your constraint sets x''(t) to be either 0 when the IF condition is true, or it is the result of the boolean expression in the else part.

Success!
Luc

[AD_10780197]

Hi Luc,

Thanks for the quick reply. I’m indeed working in Prime 5 (although I recently gained access to Prime 10 as well). I removed the x''(t)= part from the else branch exactly as you suggested, but unfortunately, I’m still getting the same “Argument not specified” error.

[Werner_E]

1) You should use x1 and x2 in the odesolve command, not their derivatives.

2) p.gas is a variable (value 8*10^6) but is used as if it was a function here on the LHS.

This can't work

3) k.gas is here used like a function but its nowhere defined, neither as a function nor as a variable

The same undefined variable (k.gas) also occurs in the else branch in the definition of x.2''

Simply deleting the undefined expressions makes the solve block work, but sure that's not what you ar looking for. Of course x2 is constant zero.

I am not sure if the solve block will work using the branching (if ... else) on the RHS but as a first step you have to fix the problems with p.gas and k.gas.

[AD_10780197]

Hi Werner,

Thanks again for the detailed feedback!

Since my last message, I made this model in the past in Python  (slightly different parameter set), and that version solves without errors, so the core equations do work. I linearized the orifice equation using the chain rule, but the overall structure remains the same.I’ll walk through the sheet once more and, line by line, compare it with the Python code to spot the mistake and your remarks below

I’ve already:

1. switched the odesolve vector to use x1 and x2 (not their derivatives sorry );

2. renamed the scalar p_gas to p_gas0 and added a proper p_gas(x) function;

3. introduced the missing k_gas(x) definition. ( Gas stiffness as a function of the position of x2)

With those changes the solve block runs, but the results still look off. I’ll keep checking and let you know what I find to solve this, maybe you can see my mistake in one eye blink. The function shared by you was my initial starting point see sheet ref 12 working fine. The add on was the extra motion equation for the accumulator. Thanks a lot for your help so far it’s really appreciated! 

[Werner_E]

---

@AD_10780197 wrote:

I’ve already:

1. switched the odesolve vector to use x1 and x2 (not their derivatives sorry );

2. renamed the scalar p_gas to p_gas0 and added a proper p_gas(x) function;

3. introduced the missing k_gas(x) definition. ( Gas stiffness as a function of the position of x2)

Hmm, I loaded your "A14" sheet and you did NOT changed the variables/function you solve for in the odesolve block
You renamed p_gas to p_gas0 but you did not defined functions p_gas and k_kas. You used boolean expressions inside the solve block and not function definitions in front of it as shown in my previous answer below.

See the attached Prime 10 sheet.

PS: I noticed that in the A12 sheet you provided you tried to use units in the plot. To do so you should use units right from the start!

Or you redefine like x1(t):=x1(t)*1/bar, but using units throughout sure is preferable.

[Werner_E]

Could it be that you had something like the following in mind?

I defined a function k.gas in front of the solve block so it could be used in the definition of x2''. I had no idea what you had in mind with your "function" p.gas, so I deleted it and just used the already defined variable p.gas as you had set it up.

vgas and v1 are defined after the solve block using the results for x1 and x2. x2 still seems to be constant zero.

[AD_10780197]

Dear Werner,

Thanks again — I’ve made good progress due to you, and the system is almost working as expected. When I add an extra differential equation based on the conservation law (mass balance), everything still works correctly and behaves as expected. See the named working as the latest extension name in the attachment

However, when I try to make the gas pressure and gas stiffness explicitly dependent on the dynamic motion of x2(t), the model starts to fail. It simply stops working as expected. I’ve double-checked the free-body diagram. I know that in steady-state, the pressures should be balanced, so this behavior is unexpected

 I will try to solve this, maybe you have some tips. Perhaps the issue is due to the highly nonlinear nature of the gas-related functions (as shown in the appendix). I’m considering linearising the expressions for gas stiffness and gas volume to improve stability. It’s also unclear whether this inconsistency only affects the transient phase or if it also has an impact in steady-state conditions.

Best regards,

Albert 

[Werner_E]

---

However, when I try to make the gas pressure and gas stiffness explicitly dependent on the dynamic motion of x2(t), the model starts to fail.

---

I looked at your "almost" file but found no failing calculation ?

[AD_10780197]

Dear Werner,

I apologize for attaching the wrong file earlier. Please find the correct file attached. This one is with gas and pressure as a function of the movement of the accumulator  x2(t)

Thank you in advance

With kind regards

Albert

[Werner_E]

You have to DEFINE the functions (using := ) - not using the boolean comparison fat equal. 
And the name of the formal argument does not matter at all when you are defining a function

And you have to do this IN FRONT of the solve block. Your expressions were positioned somewhere underneath(!) the solve block.

As already written it doesn't make much sense to put unit "bar" in the unit placeholder of the plot. Your whole calculation is done dimensionless because for some reason you stripped the units (you may consider not doing so).

You see that Prime realizes that the results of sol.p1 are dimensionless and so corrects the wrong pressure unit "bar" which you typed in by dividing by the default pressure unit (which is Pascal in the unit system SI you ar using).

You could achieve the very same scaling effect by typing 10^5 in the unit placeholder

or by dividing the values by 10^5 (or multiplying them by 10^-5)

While I still would prefer not stripping the units and using them throughout you could also add the units later. If you know that the dimensionless results are in Pascal, then add that unit. Later in the plot you can demand the display in bar.

BTW, somewhere in your sheet you write "Program has problems with units" and you start stripping all units.
Better approach would be to find the reason for the problem (guess you ran into some kind of unit mismatch) and fix it rather than to get rid of the units altogether.

In Prime you may use units for most calculations - one exception are the 3D plots where no units are allowed at the x- and y-axis. Here a workaround is necessary, but that should be no reason to get rid of all units for the whole (or for most) of the calculations.

Prime 10 sheet attached

[Werner_E]

See the attached file and the remarks therein - it uses units wherever possible.

Its a crude fix and I simply replaced assignments likebyor added the unit m³ in the definition of A.accu.

I often wasn't sure about the correct/desired units (m³ or L, J or kW*hr, etc.).

There also may be leftovers from my various attempts.

So the sheet sure still needs quite some rework.

[AD_10780197]

Hi Werner,

First of all, thanks again 🙏 for your earlier help — your tips have helped me understand what was going wrong and given me a clear direction to continue. I still plan to clean up the sheet properly myself, but right now I'm stuck on one specific issue related to the solve block and a position-dependent function. I have one small question, then I will stop stalking. If I want to link with the second mass spring damper, how can I do that? You already shared with a x, but this x is a vector linear and not linked with x2(t). Also, tried to use x2 and seems not to work. I removed the stiffness block and changed the < to > to test the ODE x2(t)

The gas pressure in the accumulator depends on the position of the piston x2(t) , and I’ve defined the gas pressure as a separate function outside the solve block:

​

But when I do this, Mathcad returns the following error:

"The target of an invocation has caused an exception."

Thanks in advance, maybe you’ll spot it right away.

[Werner_E · Accepted solution]

You do it exactly as you did 😉

And the name of the formal argument (currently x2) when you define function p.gas does not matter at all.

The error message is not helpful at all. The reason is a unit mismatch.

You defined p.gas0 to be 8*10^6, but it should be 80 bar. Actually I forgot to add the unit in the sheet where I tried to use units allover.

Not sure why this unit mismatch did not show in the sheet I posted - maybe because the else branch in the definition for x2'' never was called!?

I also noticed in that older sheet that the usage of function k.gas also would result in a unit mismatch even with the correct unit for p.gas0:

Anyway, in your new sheet k.gas is not used and once you add the correct unit for the initial gas pressure the solve block works OK.

Prime 10 sheet attached

[AD_10780197]

Dear Werner, I am super happy with your support. I learned a lot, and now I am tidying up the mess I made, but without you, I wouldn't have succeeded. I have learned a lot.

Thank you so much.

[vnamboodheri]

Hello @AD_10780197,

It looks like you have some responses from our community champion. If any of these replies helped you solve your question please mark the appropriate reply as the Accepted Solution. 

Of course, if you have more to share on your issue, please let the Community know so other community members can continue to help you.

Thanks,
Vivek N.
Community Moderation Team.