Energy not conserved in non-linear systems containing Annex60.Fluid.Interfaces.StaticTwoPortConservationEquation #282

Mathadon · 2015-06-30T13:55:02Z

When a MixingVolume is configured to steady state (i.e. Annex60.Fluid.Interfaces.StaticTwoPortConservationEquation is used) energy is in some cases not conserved at zero flow. The problem arises when the following equation is used in a non-linear system:

port_b.h_outflow = inStream(port_a.h_outflow) + Q_flow * m_flowInv;

A resulting non-linear algebraic loop may look like:

algorithm // Torn part
    vol.heatPort.T := 273.15+0.0002390057361376673*vol.ports[2].h_outflow;
    prescribedTemperature.Q_flow := 0.01*((290-vol.heatPort.T)*sin(time));
    vol.heatPort.Q_flow := prescribedTemperature.Q_flow*(1+prescribedTemperature.alpha
      *(vol.heatPort.T-prescribedTemperature.T_ref));

  equation // Residual equations
    0 = vol.ports[2].h_outflow-(vol.steBal.port_a.h_outflow_inStream+
      vol.heatPort.Q_flow*vol.steBal.m_flowInv);

In the last equation vol.steBal.m_flowInv is zero around zero flow. This means Q_flow can take any value in the last equation and the equation will still be satisfied.

An example:

results in:

ie Q_flow is not 0 when time > 1s. Note that this problem exists regardless of if prescribedHeatFlow is used.

Note that this is the same problem as #280 , and this problem can again be solved by setting use_safeDivision=false. The only difference is that the singularity of this equation does not lead to a division by zero since the equation is used in a non-linear loop, where less symbolic manipulations are performed. The implication is however that the singularity is not detected and that this leads to non-physical results.

I think we'll need to think of an approach for handling this properly, or provide good parameters to the user, like suggested in #280 .

The text was updated successfully, but these errors were encountered:

…ound zero flow this is for issue ibpsa#282

Mathadon · 2015-06-30T14:08:08Z

See https://github.com/Mathadon/modelica-annex60/tree/issue282_singularMixVol for an example
Annex60.Fluid.MixingVolumes.Examples.MixingVolumeNonLinSysZeroFlow

Mathadon · 2015-06-30T14:18:09Z

Setting use_safeDivision=false can lead to other problems for the same model so it's getting tricky. I'll try to look into it further but any input is welcome.

mwetter · 2015-06-30T22:45:41Z

m_flowInv approaches zero only near zero flow. This is required for m_flowInv to be smooth and to be equal to 1/m_flow for reasonably large flow rates. I do not know how to otherwise robustly handle the problem described by equation (1) in http://www.tandfonline.com/doi/pdf/10.1080/19401493.2013.765506

If m_flow is that small to be inside the region where this regularization happens, then no heat should be exchanged anyway as any small perturbation would yield to very large outlet temperatures in such steady-state models. All we attempt to do is to get robustly through the zero flow region.

If you want to exchange heat at or near zero mass flow rate, you will need a dynamic model which then won't use this equation.

The parameter use_safeDivision is only for testing and should not be changed for normal use.

Mathadon · 2015-07-01T06:27:01Z

I agree that no heat should be exchanged when m_flow is small, but the problem is that it does exchange heat.

I think I found a solution. I'll try to explain.

In any model vol.heatPort.Q_flow and vol.heatPort.T will need to be computed somehow. Let's consider 3 cases how this can happen, depending on what is connected at the HeatPort. Depending on the case a different model formulation may be required. All these cases are for m_flow=0.

Case 1: `vol.heatPort.Q_flow` is determined by a `prescribedHeatFlow` and the input value is computed using only 'external' variables

Note that the example above does not classify under this category since the input of the prescribedHeatFlow requires vol.heatPort.T.

In this case Q_flow is known and h_outflow needs to be computed from this value. We need equation:

port_b.h_outflow = inStream(port_a.h_outflow) + Q_flow * m_flowInv;

since this can be solved for port_b.h_outflow = inStream(port_a.h_outflow). I.e. the heat flow entering the mixing volume is lost, but at least the model does not crash. The user should set vol.prescribedHeatFlow=true.

Case 2: `vol.heatPort.T` is determined by a `prescribedTemperature` (or a state value?) and the input value is computed using only external variables.

T (and thus also port_b.h_outflow) is known and Q_flow needs to be computed from this value. In this case

port_b.h_outflow = inStream(port_a.h_outflow) + Q_flow * m_flowInv;

will not work since it leads to a division by zero when solving for Q_flow. Therefore we should use

port_a.m_flow * (inStream(port_b.h_outflow) - port_a.h_outflow) = +Q_flow;

Now a second problem arises since the 'second' outlet temperature (port_a instead of port_b) still needs to be computed from:

port_a.m_flow * (inStream(port_a.h_outflow) - port_b.h_outflow) = -Q_flow;  (1)

Since Q_flow is known by now we again get a division by zero. It is important to note at this point that vol.prescribedHeatFlow = false should be set by the user. Therefore the model is only used when allowFlowReversal = false because dynBal would be used otherwise:

useSteadyStateTwoPort=(nPorts == 2) and 
      (prescribedHeatFlowRate or (not allowFlowReversal)) and (
      energyDynamics == Modelica.Fluid.Types.Dynamics.SteadyState) and (
      massDynamics == Modelica.Fluid.Types.Dynamics.SteadyState) and (
      substanceDynamics == Modelica.Fluid.Types.Dynamics.SteadyState) and (
      traceDynamics == Modelica.Fluid.Types.Dynamics.SteadyState)

Since allowFlowReversal=false we can use this to our advantage. I propose two options for changing equation (1):

Option 1:

Use:

port_a.h_outflow = inStream(port_b.h_outflow) - Q_flow * m_flowInv;

since Q_flow is already known at this time we can resort to the solution of case 1.

Option 2:

As suggested in #281 we can remove the calculation completely since its result should not be used anyway:

port_a.h_outflow = Medium.h_default;

Case 3: `vol.heatPort.T` and `vol.heatPort.Q_flow` need to be computed from a (non-)linear algebraic loop

This is the case in the example on the top.

non-linear case

All equations are solved sequentially and iteratively. Assuming that no symbolic manipulation occurs we should use

port_a.m_flow * (inStream(port_a.h_outflow) - port_b.h_outflow) = -Q_flow;

since this always leads to Q_flow = 0 when m_flow=0. Other equations can then be used to identify the value of port_b.h_outflow and vol.heatPort.T.

linear case

Symbolic manipulation does occur. I'm not sure what the impact of this is but I think we should also use

port_a.m_flow * (inStream(port_a.h_outflow) - port_b.h_outflow) = -Q_flow;

In both the linear and non-linear case Q_flow and port_b.h_outflow are now known but port_a.h_outflow still needs to be computed. This results in the same problems as in case 2 and the same equations can be used to overcome this.

Note that the solution of case 2 is the same as the solution of case 3.

Summary

I propose to keep the current equations only when prescribedHeatFlow=true and to use the alternative implementation otherwise. Ie this can be implemented by setting use_safeDivision=prescribedHeatFlowRate. Then secondly we should change the equation that computes port_a.h_outflow using one of the two options I proposed.

I'll make examples and unit tests for all these cases and demonstrate the solution.

Mathadon · 2015-07-01T06:45:09Z

An important note for debugging this:

Dymola does not always generate errors when divisions by zero occur, this seems to be the case when the division is 'consistent', i.e. equation (1) can be evaluated as x=0/0. Therefore in some cases it's difficult to 'prove' that an implementation generates errors, especially with simple examples since they can be solved quite accurately. A way to overcome this is to set Advanced.Define.AimForHighAccuracy = false. This will cause non-linear systems to be solved less accurately, which results in an inconsistent division, for example x = 1E-5/0.

Mathadon · 2015-07-01T07:06:00Z

An additional advantage of this approach is that more systems will be linear and that they can therefore be solved more easily.

around zero flow updated documentation this is for ibpsa#282

Mathadon · 2015-07-01T07:59:13Z

The pull request contains an example model (Annex60.Fluid.MixingVolumes.Examples.MixingVolumeZeroFlow) that implements the changes outlined above. For each mixing volume prescribedHeatFlowRate can be changed and this will result in an error. This demonstrates how the previous implementation produced errors.

I did not implement port_a.h_outflow = Medium.h_default since we should treat this separately in #281 .

mwetter · 2015-07-01T16:25:52Z

Case 1 is OK.

For Case 2, I suggest to go with option 2, e.g., use

port_a.m_flow * (inStream(port_b.h_outflow) - port_a.h_outflow) = +Q_flow;
port_a.h_outflow = Medium.h_default; // Provided https://trac.modelica.org/Modelica/ticket/1742 shows that this is OK.

Please also write in the code and the info section that this relies on

useSteadyStateTwoPort=(nPorts == 2) and 
      (prescribedHeatFlowRate or (not allowFlowReversal)) and (
      energyDynamics == Modelica.Fluid.Types.Dynamics.SteadyState) and (
      massDynamics == Modelica.Fluid.Types.Dynamics.SteadyState) and (
      substanceDynamics == Modelica.Fluid.Types.Dynamics.SteadyState) and (
      traceDynamics == Modelica.Fluid.Types.Dynamics.SteadyState)

Also, at the above statement, a remark should be made to avoid that someone removes the (not allowFlowReversal).

For Case 3, I think we should stay with the current formulation that uses m_flowInv.
We cannot know whether computing m_flow and Q_flow requires iteration, in which
only approximate values will be computed. This is the case described in (1) in the paper
http://www.tandfonline.com/doi/pdf/10.1080/19401493.2013.765506
Your current suggestion for Case 3 also is based on "Assuming that no symbolic manipulation occurs".
Such an assumption cannot be made when developing a component model.

Using use_safeDivision will confuse others (and us in a few years). I suggest
to use an enumeration instead, with the following values:

configuration 
 - prescribedHeatFlow "heatPort.Q_flow is prescribed, i.e., it is a control input and not computed using heatPort.T";
 - prescribedTemperature "heatPort.T is prescribed, for example using a control signal or the solution of a differential equation";
 - computedHeatFlow "heatPort.Q_flow is computed based on the difference between an ambient temperature and heatPort.T";

I also don't understand the title for Case 3: You wrote "vol.heatPort.T and vol.heatPort.Q_flow need to be computed...". But in fact, only Q_flow is computed. If T were computed, it would be assigned, in which case we have Case 2.

Do you agree with Case 3?
Could you implement the above using the enumeration (and document accordingly)? I can then also run the large tests on the Buildings library before moving it to the master as this is a significant change.

Mathadon · 2015-07-01T17:14:49Z

Case 1: Ok
Case 2: Will fix this
Case 3: I disagree:

Note that keeping mFlow_inv for case 3 is simply not an option since energy is not conserved. I've had a simulation where a storage tank was charged by itself when m_flow=0 because of this.

I'm not sure what you mean but T and Q_flow both do need to be computed from the algebraic loop as is illustrated by the following piece of code (1):

  // Nonlinear system of equations number = 1.
    // Start values for iteration variables of non-linear system of 1 equations: 
    //   volNonLinSys.ports[2].h_outflow(start = 83680.0)
  algorithm // Torn part
    volNonLinSys.heatPort.T := 273.15+0.0002390057361376673*volNonLinSys.ports[2].h_outflow;
    preHea2.Q_flow := 0.01*((290-volNonLinSys.heatPort.T)*volNonLinSys.heatPort.T
      ^0.6*time);
    volNonLinSys.heatPort.Q_flow := preHea2.Q_flow*(1+preHea2.alpha*(
      volNonLinSys.heatPort.T-preHea2.T_ref));

  equation // Residual equations
    0 = volNonLinSys.ports[1].m_flow*(sou2.ports[1].h_outflow-volNonLinSys.ports[2].h_outflow)
      +volNonLinSys.heatPort.Q_flow;

It's important that the current (updated) implementation of the residual equation is used since this guarantees that Q_flow=0 when m_flow=0.

Using the old implementation with m_flowInv you get the following piece of code (2):

algorithm // Torn part
    volNonLinSys.heatPort.T := 273.15+0.0002390057361376673*volNonLinSys.ports[2].h_outflow;
    preHea2.Q_flow := 0.01*((290-volNonLinSys.heatPort.T)*volNonLinSys.heatPort.T
      ^0.6*time);
    volNonLinSys.heatPort.Q_flow := preHea2.Q_flow*(1+preHea2.alpha*(
      volNonLinSys.heatPort.T-preHea2.T_ref));

  equation // Residual equations
    0 = volNonLinSys.ports[2].h_outflow-(sou2.ports[1].h_outflow+
      volNonLinSys.heatPort.Q_flow*volNonLinSys.steBal.m_flowInv);

where Q_flow can take on every value since it is multiplied with 0 anyway.
My statement that no symbolic manipulation should occur was a bit misleading and I'd like to disregard that for now.

Regarding the paper: you there assume that delta h is computed as delta h= Q/m. This is however not the case as illustrated above by piece of code 1. The actual equation is Q= m * delta h. I also see no reason why Dymola would invert this equation except for in 'case 1', exactly because it could introduce a division by zero.

So I'd like to stick with the proposed implementation. Note that the examples show that it works.

Regarding the enumeration. I would stick with a parameter since case 2 and case 3 have identical code and therefore only 2 different cases need to exist. Moreover case 1 is pretty easy to identify and quite rare whereas the other cases are not and this could confuse the user. I do propose to rename use_safeDivision to prescribedHeatFlowRate as this is a lot clearer.

Mathadon · 2015-07-01T17:28:01Z

Note that the above problem was only introduced recently when we started using StaticTwoPortConservationEquation also when allowFlowReversal = false. Before dynBal was used when allowFlowReversal = false so the issue did not exist since this model was only used for case 1, where use_safeDivision = true was the right approach.

mwetter · 2015-07-01T21:37:12Z

I see why you get the wrong results in Case 3. So let's use your proposed formulation for Case 3, e.g.,

port_a.m_flow * (inStream(port_a.h_outflow) - port_b.h_outflow) = -Q_flow;

Renaming use_safeDivision to prescribedHeatFlowRate is good (although this is a non-backward compatible change, but it is clearer).

Please also update the info section in StaticTwoPortConservationEquation so that users don't need to read github tickets. The info section should be self-contained. It is however fine to refer to the ticket in the revision notes, but users will be scared away if they need to read pages of tickets in order to understand the info section.

The proper explanation is "prescribedHeatFlow should be true if heatPort.Q_flow is prescribed, i.e., it is computed from a control input and not from a heat transfer in which heatPort.T is part of the driving temperature difference"
("Q_flow is an input" is not clear as input is a Modelica keyword.)

Mathadon · 2015-07-01T22:21:49Z

Fixed.

I also updated and documented Annex60.Fluid.Interfaces.StaticTwoPortHeatMassExchanger since it was affected by the renaming of use_safeDivision. Unit tests give correct results.

mwetter · 2015-07-02T22:16:26Z

@Mathadon
I did some more revisions on issue282_singularMixVol

I revised the comments and refactored StaticTwoPortConservationEquation to set default outlet values if allowFlowReversal=false.

I also set a start value for prescribedHeatFlow but not a default value. This will give a warning if no value is provided.

All unit tests are fine on the Annex60 library. However, in the Buildings library, I get different translations statistics:

*** Warning: SpaceCoolingSystem3.mat: Translation statistics for initialization changed for linear, but results are unchanged.
 Old = 0, 0
 New = 0, 2

*** Warning: DataCenterContinuousTimeControl.mat: Translation statistics for initialization changed for nonlinear, but results are unchanged.
 Old = 1, 3, 1, 1, 17
 New = 1, 3, 1, 1, 19

*** Warning: DataCenterDiscreteTimeControl.mat: Translation statistics for initialization and results changed for nonlinear.
 Old = 1, 3, 1, 1, 17
 New = 1, 3, 1, 1, 19

*** Warning: DataCenterRenewables.mat: Translation statistics for initialization and results changed for nonlinear.
 Old = 1, 3, 1, 1, 17
 New = 1, 3, 1, 1, 19

I also updated BuildingsPy on the lbl repository so that it issues a warning if a parameter has a start value but no value has been assigned.

These are on the branch Annex60_issue282_singularMixVol of the Buildings library. I have not yet looked at why the nonlinear systems become bigger.
I will run larger tests over the week-end.

Mathadon · 2015-07-03T04:37:25Z

That's interesting. I'll have a look at it!

Mathadon · 2015-07-03T09:53:22Z

I cannot reproduce these results using branch buildings/Annex60_issue282_singularMixVol using Dymola 2015FD01 or Dymola2016 using Evaluate=true or false.

Can you verify their result when doing a manual translation?

This is for ibpsa/modelica-ibpsa#282

mwetter · 2015-07-17T16:26:23Z

@Mathadon : The large tests in the Buildings library work fine. I also verified the above statistics for SpaceCoolingSystem3, DataCenterDiscreteTimeControl and DataCenterRenewables and I get the smaller system (as described under 'Old').

@Mathadon This is ready to merge from my point of view after these issues are fixed:

Annex60.Fluid.Examples.Performance.Example4 fails to translate.
There are many warnings, see below.

*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Delays/Examples/Delay.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Examples/ReverseBuoyancy.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Validation/ThreeRoomsContam.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Examples/Performance/Example1v2.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Experimental/Benchmarks/AirFlow/Examples/MultipleFloorsVectors.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Sensors/Conversions/Examples/To_VolumeFraction.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Experimental/Benchmarks/AirFlow/Examples/ZoneStepResponse.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Storage/Examples/ExpansionVessel.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Experimental/Benchmarks/AirFlow/Examples/TwoFloors.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/MassExchangers/Examples/ConstantEffectiveness.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/MixingVolumes/Examples/MixingVolumeMoistAir.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/MixingVolumes/Examples/MixingVolumeMassFlow.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/MixingVolumes/Examples/MixingVolumeMFactor.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/MixingVolumes/Examples/MixingVolume.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Examples/OneOpenDoor.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/MixingVolumes/Examples/MixingVolumeInitialization.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Examples/ReverseBuoyancy3Zones.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Examples/ZonalFlow.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Sensors/Examples/TraceSubstances.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Sources/Examples/TraceSubstancesFlowSource.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Examples/ChimneyShaftWithVolume.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Sensors/Examples/Density.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Examples/NaturalVentilation.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Interfaces/Examples/ReverseFlowMassExchanger.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Examples/ClosedDoors.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Interfaces/Examples/Humidifier_u.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Examples/OneRoom.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Sensors/Examples/MassFraction.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Interfaces/Examples/HeaterCooler_u.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Examples/OneEffectiveAirLeakageArea.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Examples/ResistanceVolumeFlowReversal.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Fluid/Interfaces/Examples/ConservationEquation.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Examples/ChimneyShaftNoVolume.mos");'.
*** Warning: Parameter with start value only in 'RunScript("Resources/Scripts/Dymola/Airflow/Multizone/Examples/CO2TransportStep.mos");'.

mwetter · 2015-07-17T17:25:34Z

@Mathadon : Can you please check why Annex60.Fluid.Examples.Performance.Example4 no longer translates after I merged the master to this branch.

A related problem is that after merging this branch into https://github.com/lbl-srg/modelica-buildings/tree/issue428_fmu_stat, the model
Buildings.Fluid.MassExchangers.Humidifier_u fails the model check if allowFlowReversal=false.

In version 2.1.0 of the Buildings library, this model worked.

Revision:
The commit a909a52 corrects the bug that caused Example4 and Humidifier_u to fail during translation.
The only remaining item is to address the missing start values.

The model Fluid/Interfaces/Examples/StaticTwoPortConservationEquation.mo does not translate

mwetter · 2015-07-21T21:08:46Z

@Mathadon : I think we should look at setting prescribedHeatFlowRate=false by default by assigning the parameter in the base class. Models such as a chiller or furnace could set it to true.

Using true as a default value can give strange results, see vol[1].T when running

simulateModel("Annex60.Fluid.Examples.ResistanceVolumeFlowReversal(vol(each prescribedHeatFlowRate=true))", stopTime=10000, method="dassl", tolerance=1e-05, resultFile="ResistanceVolumeFlowReversal");

While the results are correct (because the volumes are steady-state and hence the temperature is not defined if m_flow=0), the formulation with prescribedHeatFlowRate=false does not give this behavior and leads to fewer operations.
This change will also require updating the info section and the documentation of the parameter.

Changed parameter to constant as we can assume that users won't need to change this, and hence should not be shown the complexity of having this parameter when they instanciate a volume. Update documentation. Removed non-required assignments

Mathadon · 2015-08-02T17:39:28Z

I agree that we should set prescribedHeatFlowRate=false by default since this will be the correct setting in most of the cases.

I'm not sure if it should be a constant or a parameter. Not being able to use the gui to change/see the variable seems a bit impractical. I do agree that its value should never need to be changed after compilation.

mwetter · 2015-08-11T22:09:22Z

@Mathadon If you agree with the current implementation, I (or you) can merge
#289

I merged this already to the master of the Buildings library and see no issue. In fact, lbl-srg/modelica-buildings#448 leads to about 5 to 10% faster simulation across a few large test models.

Mathadon · 2015-08-12T08:21:18Z

That's good news :) I'll do the merge!

Mathadon · 2015-08-12T08:30:27Z

This is closed by #289

Mathadon · 2017-09-22T16:14:57Z

@firasomran I cannot see the file that you attached.
Can you elaborate why you want to read the equations and why state that you cannot read ' its cmoplete system of equations'? Are some equations missing?

Mathadon · 2017-09-22T18:37:46Z

See here (http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.571.9335&rep=rep1&type=pdf) for more information about tearing / torn equations.
Essentially, all these equations are solved. But iteration variables are added only for the residual equations. This reduces the complexity that the Newton solver has to deal with.

Instead of looking at the dsmodel.mof, you can also just look at the model equations of course :)

But most equations should be in the dsmodel.mof

Mathadon · 2017-09-23T11:20:54Z

In the Modelica text view of the individual models. E.g. IBPSA.Fluid.Interfaces.ConservationEquation contains the equations for conservation of energy and conservation of mass of MixingVolumes.

Mathadon · 2017-09-25T07:19:13Z

An example based on your file:

  // Nonlinear system of equations
  // Tag: simulation.nonlinear[2]
    // It depends on the following parameters: 
    //   Kondensator.shell.geo.A_hor
    //   Kondensator.shell.geo.height_fill
    //   Kondensator.shell.phaseBorder.table.columns[1]
    //   Kondensator.shell.phaseBorder.table.tableID
    // It depends on the following timevarying variables: 
    //   Kondensator.shell.volume_liq
    // Start values for iteration variables of non-linear system of 1 equations: 
    //   Wasserstand_des_Kondensators
  algorithm // Torn part
    Kondensator.shell.phaseBorder.table.y[1] := ClaRa.Components.Utilities.Blocks.ParameterizableTable1D.tableIpo
      (Kondensator.shell.phaseBorder.table.tableID, Kondensator.shell.phaseBorder.table.columns[1],
       Wasserstand_des_Kondensators);
    Kondensator.shell.phaseBorder.A_hor_act := Kondensator.shell.geo.A_hor*
      Kondensator.shell.phaseBorder.table.y[1];
    Kondensator.shell.phaseBorder.level_abs := min(Kondensator.shell.geo.height_fill,
       max(1E-006, Kondensator.shell.volume_liq/Kondensator.shell.phaseBorder.A_hor_act));

  equation // Residual equations
    0 = Wasserstand_des_Kondensators-Kondensator.shell.phaseBorder.level_abs/
      Kondensator.shell.geo.height_fill;
  // Analytic Jacobian was produced, but it is not listed here.
  // To have it listed, set
  //   Advanced.OutputModelicaCodeWithJacobians = true
  // before translation. May give much output,
  // because common subexpression elimination is not activated.
  // End of nonlinear system of equations

This is nonlinear system number 2 with 1 iteration variable and 3+1 equations. So in the Dymola statistics this would be listed as a non-linear system of equations with size 4 before reduction and size 1 after reduction.

The torn equations are the equations that could be solved for the variable in the left hand side of the equation. These variables can be evaluated using the iteration variable(s) and 'inputs' to the algebraic loop only. I.e. no other variable values of the algebraic loop are required.
Since not all variables could be solved like this, one or more equations remain that could not be solved for a variable. These are the residual equations. The residual equations are solved by a Newton solver such that their residuals (i.e. the left hand side minus the right hand side of the equation) become zero. The algebraic loop is solved when this is the case.

So the model equations are all four equations (3 torn and 1 residual equation) that are listed in the above excerpt.
I hope this helps!

thorade · 2017-09-25T08:53:07Z

@firasomran
Dymola 2018 introduced some new debugging and diagnostics tools that might be helpful for you as well.
https://www.3ds.com/products-services/catia/products/dymola/latest-release/

Mathadon · 2017-09-25T13:18:44Z

I think they should be included in the translation statistics, but they simply don't have a number in dsmodel.mof?

Mathadon added a commit to Mathadon/modelica-annex60 that referenced this issue Jun 30, 2015

added example model that demonstrates that energy is not conserved ar…

363add6

…ound zero flow this is for issue ibpsa#282

Mathadon mentioned this issue Jul 1, 2015

Division by zero in Annex60.Fluid.Interfaces.StaticTwoPortConservationEquation #280

Closed

Mathadon added a commit to Mathadon/modelica-annex60 that referenced this issue Jul 1, 2015

fixed implementation of StaticTwoPortConservationEquation.mo

e067435

around zero flow updated documentation this is for ibpsa#282

Mathadon mentioned this issue Jul 1, 2015

Issue282 singular mix vol #283

Closed

Mathadon mentioned this issue Jul 2, 2015

Issue282 singular mix vol #289

Merged

Mathadon mentioned this issue Jul 3, 2015

Separate energy dynamics for heat exchanger solid in Buildings.Fluid.Storage.StratifiedEnhancedInternalHex lbl-srg/modelica-buildings#434

Closed

mwetter added a commit to lbl-srg/modelica-buildings that referenced this issue Jul 17, 2015

Added prescribedHeatFlowRate

3941ec9

This is for ibpsa/modelica-ibpsa#282

mwetter added a commit that referenced this issue Jul 18, 2015

Added example that shows that #282 has a bug

33bd6f9

The model Fluid/Interfaces/Examples/StaticTwoPortConservationEquation.mo does not translate

mwetter added a commit that referenced this issue Jul 18, 2015

Bug fix for #282

a909a52

Mathadon closed this as completed Aug 12, 2015

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Energy not conserved in non-linear systems containing Annex60.Fluid.Interfaces.StaticTwoPortConservationEquation #282

Energy not conserved in non-linear systems containing Annex60.Fluid.Interfaces.StaticTwoPortConservationEquation #282

Mathadon commented Jun 30, 2015

Mathadon commented Jun 30, 2015

Mathadon commented Jun 30, 2015

mwetter commented Jun 30, 2015

Mathadon commented Jul 1, 2015

Mathadon commented Jul 1, 2015

Mathadon commented Jul 1, 2015

Mathadon commented Jul 1, 2015

mwetter commented Jul 1, 2015

Mathadon commented Jul 1, 2015

Mathadon commented Jul 1, 2015

mwetter commented Jul 1, 2015

Mathadon commented Jul 1, 2015

mwetter commented Jul 2, 2015

Mathadon commented Jul 3, 2015

Mathadon commented Jul 3, 2015

mwetter commented Jul 17, 2015

mwetter commented Jul 17, 2015

mwetter commented Jul 21, 2015

Mathadon commented Aug 2, 2015

mwetter commented Aug 11, 2015

Mathadon commented Aug 12, 2015

Mathadon commented Aug 12, 2015

Mathadon commented Sep 22, 2017

Mathadon commented Sep 22, 2017

Mathadon commented Sep 23, 2017

Mathadon commented Sep 25, 2017

thorade commented Sep 25, 2017

Mathadon commented Sep 25, 2017

Energy not conserved in non-linear systems containing Annex60.Fluid.Interfaces.StaticTwoPortConservationEquation #282

Energy not conserved in non-linear systems containing Annex60.Fluid.Interfaces.StaticTwoPortConservationEquation #282

Comments

Mathadon commented Jun 30, 2015

Mathadon commented Jun 30, 2015

Mathadon commented Jun 30, 2015

mwetter commented Jun 30, 2015

Mathadon commented Jul 1, 2015

Case 1: vol.heatPort.Q_flow is determined by a prescribedHeatFlow and the input value is computed using only 'external' variables

Case 2: vol.heatPort.T is determined by a prescribedTemperature (or a state value?) and the input value is computed using only external variables.

Option 1:

Option 2:

Case 3: vol.heatPort.T and vol.heatPort.Q_flow need to be computed from a (non-)linear algebraic loop

non-linear case

linear case

Summary

Mathadon commented Jul 1, 2015

Mathadon commented Jul 1, 2015

Mathadon commented Jul 1, 2015

mwetter commented Jul 1, 2015

Mathadon commented Jul 1, 2015

Mathadon commented Jul 1, 2015

mwetter commented Jul 1, 2015

Mathadon commented Jul 1, 2015

mwetter commented Jul 2, 2015

Mathadon commented Jul 3, 2015

Mathadon commented Jul 3, 2015

mwetter commented Jul 17, 2015

mwetter commented Jul 17, 2015

mwetter commented Jul 21, 2015

Mathadon commented Aug 2, 2015

mwetter commented Aug 11, 2015

Mathadon commented Aug 12, 2015

Mathadon commented Aug 12, 2015

Mathadon commented Sep 22, 2017

Mathadon commented Sep 22, 2017

Mathadon commented Sep 23, 2017

Mathadon commented Sep 25, 2017

thorade commented Sep 25, 2017

Mathadon commented Sep 25, 2017

Case 1: `vol.heatPort.Q_flow` is determined by a `prescribedHeatFlow` and the input value is computed using only 'external' variables

Case 2: `vol.heatPort.T` is determined by a `prescribedTemperature` (or a state value?) and the input value is computed using only external variables.

Case 3: `vol.heatPort.T` and `vol.heatPort.Q_flow` need to be computed from a (non-)linear algebraic loop