ibpsa / modelica-ibpsa

Modelica library for building and district energy systems developed within IBPSA Project 1
https://ibpsa.github.io/project1
145 stars 84 forks source link

reformulation of basicFlowFunction_der to avoid square root of negative #723

Closed Mathadon closed 7 years ago

Mathadon commented 7 years ago

IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp_der is currently implemented as:

 m_flow_der := (if noEvent(dp>dp_turbulent) then 0.5*k/sqrt(dp)
                elseif noEvent(dp<-dp_turbulent) then 0.5*k/sqrt(-dp)
                else 1.25*k/m_k-0.75*k/m_k^5*dp^2)*dp_der;

this is fine for modelica tools, but this causes problems in casadi, which always evaluated all branches of the if-then-else. Since all branches are always evaluated, the square root of negative numbers can be computed, which leads to problems. A reformulation that solves this is:

 m_flow_der := (if noEvent(dp>dp_turbulent) then 0.5*k/sqrt(abs(dp))
                elseif noEvent(dp<-dp_turbulent) then 0.5*k/sqrt(abs(dp))
                else 1.25*k/m_k-0.75*k/m_k^5*dp^2)*dp_der;

I don't think this affects computation time so it make sense to put this formulation in the library.

Mathadon commented 7 years ago

Or even better:

 m_flow_der := (if noEvent(abs(dp)>dp_turbulent) 
                then 0.5*k/sqrt(abs(dp))
                else 1.25*k/m_k-0.75*k/m_k^5*dp^2)*dp_der;

and using even fewer operations:

 m_flow_der := (if noEvent(abs(dp)>dp_turbulent)
                then 0.5*k/sqrt(abs(dp))
                else k/m_k*(1.25-0.75*(dp/dp_turbulent)^2))*dp_der;