"Modularization" of the program

[pw]

Hello all,

I work on modeling / simulation of some thermal energy system in Mathcad 15. The existing model works on the equation oriented approach, at the end all the mass and energy balance equations are gathered in the given-find block. Now I want to "modularize" the program - the components of the system are to be black-boxes, which have inputs and outputs (more precisely input and output streams, which are vectors of parameters). I was thinking about defining modules as functions, but the problem is I can not use boolean equal sign inside a function, that I need, since the whole system of equations / functions can only be solved numerically, in case of Mathcad in Given-Find block. So, every black-box would be expressed as

Function(X,U) = Y

where

X - vector of input streams (if particular black-box has more than one input stream this vector will often contain nested vectors)

U - vector of internal parameters

Y - vector of output streams (if particular black-box has more than one output stream this vector will often contain nested vectors)

As each module is connected to some other by stream(s), output stream of one module can be at the same time inlet stream of the subsequent module. For example, output stream of the superheater is simultaneously input stream of the turbine.

I would be very grateful to any ideas how to resolve the above described problem. Thanks in advance for help.

Piotr

[HarveyHensley · Accepted solution]

I think all you need to do is parameterize the solve block for a unit. The parameterized function that is created can then be used in a program function. You can build program functions for each unit and then a master function for a process loop.

I have built a process flowsheet similar to what you are attempting. You may get some ideas by going to my blog and selecting the "process simulation" category on the right side of the page.

[Fred_Kohlhepp]

It is possible to define a solve block as a function

[pw]

Thanks for reply. Harvey, do you think about the same as Fred? I suppose yes. I start dealing with it. One problem I see now are variables names, that in every given-find block must be unique, because they are globally visible. When dealing with functions, variables have local names, so there is no problem. Is there any way to solve this problem?

Harvey, I have visited your blog, it is very interesting, however I have not find any Mathcad files or pdfs. I would like to see the process flowsheet that you mentioned above - is it possible?

Also, I would like to ask you how would you model thermal energy system (combined cycle in my example). The purpose of the model is mass and energy balance. There are no chemical reactions and change of working media composition (except gas turbine combustion chamber). The working media are air, methane, combustion / flue gas, and water / steam. Air and combustion / flue gas are modeled as 5 components vectors (N2, O2, CO2, Ar and H2O(g) fractions). I have properties libraries for all the working media (dll files). Philosophy that I figured out is shortly described in my first post (black-boxes (modules) linked witch each other by streams). Every mass stream is characterized by 3 independent parameters of state: mass flow rate, pressure, and specific enthalpy, expressed as 3 element vector. In the system there are also heat rate streams and mechanical work streams. When using EO approach it takes too much time to re-model the program if configuration of the system is changed - that is why I want to change the "philosophy" of modeling, for example by closing the units in functions or solve blocks. Can you propose some elegant, coherent "mind model"? From your experience, what would you advise? I will be very grateful for your help.

If anybody has something to add, please feel free

Piotr

[HarveyHensley]

Piotr,

I don't think my example program is quite what you are looking for after all. I used a lot of Prode Properties functions for one thing. Also, I didn't use any solve blocks it appears. But I think the direction you plan to go looks sound, and in fact I might want to see if it would work for my example also. Here's how your case would work.

You have a solve block for a unit.

Given

equations

At the end of the block you have:

Unit-A(X,U) := Find(Y)

This just establishes the function Unit-A, it doesn't actually solve it.

So anytime you want the matrix Y from that type of unit, you have a statement like:

Y1 := Unit-A(X1,U1) which gives you the desired numerical matrix, Y1, for the particular values, X1 and U1.

X and Y are probably matrices because each stream has composition flows, and the other state variables. Each column in X and Y would represent a stream. You may be able to get by without nested arrays unless you want to contain multiple phases in a single stream.

I don't think there will be a problem with the global vs. local names provided you are careful with the new guess values prior to each solve block definition.

If you need additional programming for say a loop, you can then use the Unit-A() function inside the program function.

Hope this helps.

[pw]

Harvey,

In my program each stream has of course its composition, but for air / gas this composition is loaded as an argument in the property functions, and for water / steam properties functions it is not - the program knows if it is gas / air or water / steam by function name. So, when I have a stream that is defined by m, p and h, the composition is hidden in h. For example, h for air / gas can be expressed as h(T,C), for water / steam can be expressed as h(p,T), h(p,x), h(p,s) etc. Maybe your approach is better, but my works too. C is a composition vector.

Also, I am not sure if X and Y should be a matrix with the columns representing streams, instead of vectors. This approach would be good, if every unit/black-box/module (whatever name it has) would have only mass streams in and out. But there are also heat rate streams and mechanical work streams going in or out the unit. For example, pump has 2 input streams: mass stream and mechanical work stream, and one output stream (mass stream). So that is why I am going to deal with nested vectors. The structure of such a vector seems to be more consistent, if all the streams are "encapsulated". For pump example, X vector will contain one (nested) vector of mass stream (which has m, p and h properties), and one mechanical work stream, which is a scalar.

If my approach is not correct from your point of view, please let me know why. Thanks for your involvement in this topic.

Piotr

[HarveyHensley]

Piotr,

You are correct...the nested array is what you need since you have more than material streams. This is what I am also using for the integrated FCC model discussed in my blog. The solution to the giant solve block is a nested array with 10 sub-arrays (vectors and matrices) and scalars.

[pw]

Harvey,

I was wonder if it is possible to solve a unit expressing its function as

Unit-A(X,Y,U) = 0

And how would in this case Given-Find structure look like? Unit-A(X,Y,U) = Find(0)? Maybe there is some other way to find the roots of this function?

Piotr

[HarveyHensley]

I'm not sure what you are after. I will repeat what I said earlier and put it in a form similar to what you wrote.

Unit-A(X,Y,U) := define function

known: X:= xxxx U:= uuuuu

guesses Y:= something

Given

Unit-A(X,Y,U) = 0

Unit-A-Out(X,U) := Find(Y)

Y1:= Unit-A-Out(X1,U1) value of Y given X1 and U1

I have assumed that you are after the output, Y, of the unit. If you don't know any, then the Find would be:

Unit-A-Out(X,Y,U ):= Find(X,Y,U)

Soln := Unit-A-Out(X1,Y1,U1) Soln is a three element vector

X := Soln0 Y := Soln1 U := Soln2

[pw]

Thank you very much, Harvey,

I have started to implement the model, without parametrization - until it will work. The code for HRSG economiser (cross-counterflow) looks as follows:

When I trace error, it says "Initial guess must consist of numbers and matrices". h_pT and Tsat_p are steam tables function names.

Probably there is some problem with nested vectors elements.

Is this good way? And do you know how to fix / force it to work?

Piotr

[HarveyHensley]

I presume you have ORIGIN = 1 based on your indexes.

Try putting the guess values above the "Given".

If that doesn't work, please attach a worksheet. There is another way to write the solve block, but I think your method should work.

[pw]

ORIGIN is set to 1. Units are kg/s for mass flow rate, bar(a) for pressure, kJ/kg for specific enthalpy, and deg C for temperature. Please download steam tables from cdejonge.com - Download free steam tables for mathcad plugin and put the dlls to User EFI folder. Does not work for Prime, I suspect.

Piotr

[HarveyHensley]

Moving the guesses isn't the answer. Attached is the way I would do it.

I used a constant for the enthalpy equation since I don't have your two functions.

[pw]

I published post half an hour after you but it had been moderated for a few hours don't know why. Thanks for the attached worksheet. Nevertheless, I would like to ask you, if it is possible to use nested vectors instead of matrices, like in the original worksheet. Since, like I mentioned, I will also have heat and work streams.

[HarveyHensley]

Yes you can use nested vectors. I need to see your worksheet to see why your way wouldn't work. I suspect it may be due to incorrect subscripts, but I can't tell by looking at the picture.

[pw]

It is attached above your previous post. My post of Mar 28, 12:21 PM. Do you see it?

P.

[HarveyHensley]

Yes I see it now. Sorry I missed it earlier. I'll take a look at it now.

[HarveyHensley]

Mathcad apparently doesn't like nested arrays as initial guesses. I broke the Y array into two vectors and got the solution attached.

I still think it would be a good idea to keep your material streams in a matrix. This matrix can then be an element in a nested array along with the energy stream matrix and other results. The Find operation naturally creates the nested matrix as output, as shown by my attached solution. Note how Y1 and Y2 became elements in a nested array in S. By keeping the material streams in a matrix element of the solution, it will be easy to print a stream report by printing that element.

[pw]

Yes it works, you help me a lot. I am much closer form finding "final solution", I think. Now I want to go for expressing function by Unit-A(X,Y,U)=Find(X,Y,U) - this should be more universal, I suppose.

I would like to create a set of functions, that will be able to solve thermodynamic cycle or its part, for example HRSG. Suppose we have a single pressure HRSG, which is composed of economiser, evaporator, and superheater. Given data are gas turbine exhaust parameters (m, p, h), superheated steam p and h, feedwater temperature and internal parameters of the units (pressure drops, heat exchange efficiencies, evaporator pinch point, economiser approach point etc) - just enough to solve all the remaining parameters of all the streams. Some parameters in streams can be expressed implicitly, for example feedwater enthalpy before economiser is a function of p and T, where T is known, but p will be calculated according to pressure drops during computation. So we know only selected stream parameters, and they can be input or output params. I solved such a case in an equation oriented approach. Now I would like to solve it in "modular" approach, so I need to have 3 functions, one for every unit. Connections of the units by streams could be expressed for example by connectivity / incidence matrix. Can you tell me if it is possible to do this in such a manner in Mathcad?

I enclose two worksheets. First is a slightly changed version of your last, I changed names a bit that now looks more consistent from my point of view. Just for your information. The second is my first try in expressing function by Unit-A(X,Y,U)=Find(X,Y,U), but this does not work. Somewhere I went wrong.

I appreciate your involvement in this topic, Harvey.

Piotr

[HarveyHensley]

It appears that you want to build a simulator like Hysys that solves blocks backwards and forwards instead of just forwards like most sequential modular flowsheet programs. Can it be done? I always assume it can be done in Mathcad, until proven otherwise. Even then, I assume the objective can still be met, but maybe the implementation will be different than initially desired.

I think you will end up with a bunch of functions that you collect into one solve block and solve them simultaneously. Seems possible.

I have corrected your test 2 file. See the comments in the file.

[HarveyHensley]

There was still a problem with my previous file. There needs to be one variable in the Find() arguments for each equation. If you have more, it doesn't know which ones to keep constant. In the original version U was changed by the Find operation. The second, corrected method solves for only the Y's so the solution is different.

This has implications for what you eventually want to do...solve forward and backward. You will have to somehow indicate which variables are to be found and which are to be left as specs or input. That could get messy.

Note that sometimes you need to change the guesses for a variable that is to be solved. You are allowed to have that variable on both sides of the Find definition.

R(X,Y) := Find(Y)

An example where this can be used is for a reactor where Y is a vector of flow rates. In order to provide a good initial guess, Y needs to be changed according to the feed rate or else the problem might not converge.

[pw]

I have noticed that your previous file was broken, but since I was already working on new implementation, it did not bother me. I think about the final concept, and after all, it is not far from the presented one. I will post test case of a single pressure HRSG simulation, I already have written 3 needed functions, and they seems to work bidirectional, which is very important to me. If this concept will work (should be), it greatly shorten time of building simulation models, and they will be less error-prone as well. I will post the file (whether working or not) in a few days, when I find some time to complete it. I hope for an interesting discussion, as before, Harvey. Maybe you will also take advantage from this concept in your process flowsheeting.

Piotr

[HarveyHensley]

Piotr,

You will probably need to specify which of the variables are specs by adding constraint equations. For example, say you have three elements in each of X (the inlets) and Y (the outlets). Also assume you have three equations (e.g. material and energy balances). If you include both X and Y in the Find, eg (Find(X,Y)), you will have 6 unknowns. The program will not know which to adjust, but it appears that doesn't result in termination. It does result in an incorrect result.

To correct the problem, you need to create an equation for each spec variable. Lets assume that the 3rd element of X and the 1st and 2rd elements of Y are specs. Then you would add

X[3 = spec1 Y[1 = spec2 Y[2 = spec3

Now you have 6 equations (3 original + 3 spec constraints) and 6 variables. This will allow you to use Find(X,Y) and it should only change X[1, X[2 and Y[3. See what I mean when I said this can get messy?

You said in your last note that you are getting forward and backward results. Be sure to check to see that the results haven't changed some of your spec variables.

I will look forward to your progress. This is interesting.

[pw]

Harvey, thanks for trakcing this thread. My assumption that the functions seem to work bidirectional was too optimistic. Currently I try to force ECO function to work bidirectional, but my attempts are slipping down to nothing. The function Y_ECO(X,U) := Find(Y), or rather Y_ECO(X.1, X.2, U) := Find(Y.1, Y.2) works perfectly foreward, that is when I put known X and U into the function parameters list. Obtained Y values are correct. I have 6 unknowns, so I have written 6 equations in the given-find definition block: two for conservation of mass (one for water side and one for flue gas side), two for conservation of momentum (one for water side and one for flue gas side), one for conservation of energy, and one remaining equation - for specific enthalpy of outlet water. Can you please take a look on that? I am not sure if I understand your previous statement properly. I enclose two worksheets. The difference is highlited by red colour. By yellow colour I marked given (known) parameters of streams.

Piotr

[pw]

Correct files

[HarveyHensley]

Yes, you have 6 unknowns, but the Find argument list has 16 variables (6 X's, 6 Y's, 4 U's). You can get rid of 4 of them by removing U from the Find (). The others you have to include as constraints in the solve block as shown below.

I had to change the picture above because I had the stream names wrong. To me, it would be easier if the water stream, both in and out, were named X and the gas stream named Y. Then the index 1 and 2 would refer to in and out locations. To go further, you might consider establishing a naming convention for all types of utility streams. Examples: W = water/steam, TO = thermal oil, S = solids (e.g. catalyst), A = air. General process streams that undergo chemical change might be F for feeds, V for vapor outlet, L for liquid outlet. Just a suggestion.

[pw]

OK, I have 16 variables, so I have arranged such an equation system in G-F block:

This time I use Y_ECO function. The number of equations seems to be corresponding the number of variables, so I do not get why it still does not work.

Referring to your idea to assign X to water, and then index 1 and 2 would refer to in and out locations - I was thinking about it, but since I have in my cycle also mixing fluids heat exchanger (deaerator), it has 2 inlets and one outlet, so it complicates this idea.

You have written that there is no way to find the water pressure, but it is given for water outlet, and relative pressure drop is also given, so there should not be any problem to find feed water pressure.

Piotr

[HarveyHensley]

I have fixed your test 2 file, attached. You have to have all of the equations in one solve block. Also, the water enthalpy equation wasn't needed. Finally, there was no way to find water flow rate so I added a spec for the water inlet. Without it the program still would solve but it would use a different flow rate.

[pw]

Harvey, I think that your solution is not yet what we look for, because you have closed particular constraint values (from particular case) in ECO function definition block. The point is to achieve reusability, so the function definition block must be fixed. Solving particular cases, I mean for different data sets, everytime we should call ECO function, passing streams as an arguments. Therefore, one fixed G-F block should be a definition, and second to solve particular case. I posted this second block in my last post as an image, but something is wrong, please take a look. Also, from my point of view enthalpy equation of outlet water is needed, because it is the only one that uses approach point, which is a design parameter of every economiser. We have 16 variables, 10 of them are passed as given (known), so we need to define 6 equations in the function definition block: 2 for mass balance, 2 for conservation of momentum (pressures), 1 energy balance, and 1 is remaining. This remaining equation should use approach point, so the only one that seems to fit is enthalpy equation.

Please look at 20150506 ECO test3 file. If we write 10 (6 + 4) constraint equations + Y_ECO (which returns 6 variables), in fact we have 16 equations. 10 values are given and we look for 6, so there must be unique solution, regardless of the set of parameters known.

Looking foreward for your opinion, the discussion is very valuable.

Piotr

[HarveyHensley]

You misunderstood how U_ECO was being used. Neither it or U are needed to define the solve block. I have changed the test 2 file to make that more clear. Also, I have added a Spec vector to take care of the stream specs. This vector also isn't specified until after the solve block, so it can be changed for different units as long as the same X and Y variables are being specified.

The approach equation must not be written correctly because there are no unknowns, as I said earlier. Please check it.

I don't understand your test 3 file at all. Hopefully my updated test 2 file is what you need.