TwoWayValve model division by 0 error

[EttoreZ]

I am working on a small hydronic floor heating system for two thermal zones. It is a closed circuit that branches into two floor heating slabs (in the example approximated as resistances). The control is a simple on-off that in the example is given by a boolean pulse function. While working on the model I encountered the problem shown below:

Error: The following error was detected at time: XXX
Model error - division by zero: (1.0) / (ValDayZ.kVal^2) = (1) / (0)

The equation that raises the issue seems to be:
k = sqrt(1/(1/kFixed^2 + 1/kVal^2)); in the PartialTwoWayValveKv model where the variable k is computed as a function of kVal^2.

I tried to narrow down the possible causes but I am a stuck. A few things to mention:

• The error happens only when the circuit has two valves and the flow splitters (if you try to remove one of the two valves it disappears)

• The error sometimes happens on one valve or the other. However plotting the variables involved the values seem exactly the same.

I attach the total model to showcase the encountered error.
IssueTwoWayValveKValTotal.zip

[Mathadon]

@EttoreZ kVal=0 should not be possible since the valve leakage l should always be larger than 0. The only way I can think of that causes k=0 is when the control signal becomes negative. This can happen when you filter the valve inlet filter with use_inputFilter=true. If the time integrator makes an integration error such that the control signal becomes slightly negative instead of 0, this could happen. But even then it is weird that k would be exactly zero since k should be come negative when the control signal is too small. If you set a smaller tolerance then the problem could be solved.

General debugging approach:

• check your model equations in dsmodel.mof to see that they are as you expect
• check how k could become exactly zero

[mwetter]

@EttoreZ : Please let me know if you are still stuck after following the advice of @Mathadon .
If you are still stuck, please verify if your zip file is valid. I can't unzip it:

$ unzip IssueTwoWayValveKValTotal.zip 
Archive:  IssueTwoWayValveKValTotal.zip
  End-of-central-directory signature not found.  Either this file is not
  a zipfile, or it constitutes one disk of a multi-part archive.  In the
  latter case the central directory and zipfile comment will be found on
  the last disk(s) of this archive.
unzip:  cannot find zipfile directory in one of IssueTwoWayValveKValTotal.zip or
        IssueTwoWayValveKValTotal.zip.zip, and cannot find IssueTwoWayValveKValTotal.zip.ZIP, period.

[EttoreZ]

@mwetter I downloaded it and it unpacked fine. I made it with winRar though. This is the same file made with WinZip IssueTwoWayValveKValTotal_fin.zip .
I am trying following the advice of @Mathadon. The equations on the .mof file look fine to me. I agree with him that the issue could be caused by an error in the filter integration that leads to a negative value of y_actual. This leads to phi = 0 and kVal to 0 as consquence. It can be derived from the equations below:

ValDayZ.phi := max(0, homotopy(IBPSA.Fluid.Actuators.BaseClasses.equalPercentage
(ValDayZ.y_actual, ValDayZ.R, ValDayZ.l, ValDayZ.delta0), ValDayZ.l+
ValDayZ.y_actual*(1-ValDayZ.l)));
ValDayZ.kVal := 0.008695819912499183*ValDayZ.phi;
ValDayZ.k := sqrt(1/(13224.489795918365+1/ValDayZ.kVal^2)); .

However from the simulations results it does not seem that y_actual approaches 0. If you want we can talk about it before/after today's coordination meeting

[PMehrfeld]

I also receive very often a similar error.
In my case, the issue lays in m_flow_turbulent/k (see below).
I assume this error corresponds with this issue.

However, it seems as the solver is able to overcome the error through further iterations and the simulation continues.

Error: The following error was detected at time: 4385.110851261505

Model error - division by zero: (m_flow_turbulent) / (k) = (0.000160471) / (0)

The stack of functions is:
AixLib.Fluid.BaseClasses.FlowModels.basicFlowFunction_m_flow
AixLib.Fluid.BaseClasses.FlowModels.basicFlowFunction_m_flow(hydronicsRadiatorsValves.volRet.ports[9].m_flow, hydronicsRadiatorsValves.val[9].kVal, hydronicsRadiatorsValves.val[9].m_flow_turbulent)
First evaluation failed for non-linear solver.

[EttoreZ]

@PMehrfeld in my case the solver manages to get past it too. However I wonder in complex models that have to run for a while (in my case the two zone apartment model has to run for a yearly simulation), how much does this recurring error slows down the simulation?

@Mathadon Just a quick update after briefly talking with @mwetter. I was not able to figurate out how to extract the failed step data (in the results I have only the results where y_actual never becomes negative). However by modifiying the equation for phi from :
ValDayZ.phi := max(0, homotopy(IBPSA.Fluid.Actuators.BaseClasses.equalPercentage
(ValDayZ.y_actual, ValDayZ.R, ValDayZ.l, ValDayZ.delta0), ValDayZ.l+
ValDayZ.y_actual*(1-ValDayZ.l)));

to
ValDayZ.phi := max(ValDayZ.l, homotopy(IBPSA.Fluid.Actuators.BaseClasses.equalPercentage
(ValDayZ.y_actual, ValDayZ.R, ValDayZ.l, ValDayZ.delta0), ValDayZ.l+
ValDayZ.y_actual*(1-ValDayZ.l)));

The error disappears. @mwetter raised the issue that putting the value at exactly the leakage value would be a problem for the derivative estimation around the leakage point due to max oscillating near the leakage value. He suggested it could be still changed but to a smaller value (0.1*l) for example. It should be tested on larger models to see if there is an actual reduction of the simulation time, else the documentation could be updated to just mention that is a normal error that the user should not be too bothered by.

[Mathadon]

I'd be in favor of adding a max() that enforces the k to be positive and indeed, it's better to use max(0.1 * l, .) than max(l, .).

[mwetter]

@EttoreZ : The development branch is issue1376_valveLeakage. Please make a pull request against that branch for further testing.

[EttoreZ]

I report below a .zip file containing preliminary test on the buildings library with the suggested modification
ResultsBuildings.zip. The change does not seem to have a huge impact on the existing models of the Buildings. To further test it on the other libraries I will make a pull request on the branch issue1376_valveLeakage.