Skip to content

Commit

Permalink
Merge pull request #4358 from HansOlsson/UseSqrt
Browse files Browse the repository at this point in the history 
Use sqrt(x) instead of x^0.5 in several places.
  • Loading branch information
@hubertus65
hubertus65 authored Mar 27, 2024
2 parents 1b88744 + be7c326 commit 691f3f0
Show file tree
Hide file tree
Showing 3 changed files with 32 additions and 32 deletions.
32 changes: 16 additions & 16 deletions Modelica/Fluid/Dissipation.mo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3683,7 +3683,7 @@ Calculation of pressure loss in edged bends with sharp corners at overall flow r
PI*IN_con.d_cir else if IN_con.geometry == TYP.Elliptical then PI*(
IN_con.a_ell + IN_con.b_ell) else if IN_con.geometry == TYP.Rectangular then
2*(IN_con.a_rec + IN_con.b_rec) else if IN_con.geometry == TYP.Isosceles then
IN_con.a_tri + 2*((IN_con.h_tri)^2 + (IN_con.a_tri/2)^2)^0.5 else 0)
IN_con.a_tri + 2*sqrt((IN_con.h_tri)^2 + (IN_con.a_tri/2)^2) else 0)
"Perimeter";
SI.Diameter d_hyd=4*A_cross/perimeter "Hydraulic diameter";
Real beta=IN_con.beta*180/PI "Top angle";
Expand Down Expand Up @@ -3810,7 +3810,7 @@ Generally this function is numerically best used for the <strong>incompressible
PI*IN_con.d_cir else if IN_con.geometry == TYP1.Elliptical then PI*
(IN_con.a_ell + IN_con.b_ell) else if IN_con.geometry == TYP1.Rectangular then
2*(IN_con.a_rec + IN_con.b_rec) else if IN_con.geometry == TYP1.Isosceles then
IN_con.a_tri + 2*((IN_con.h_tri)^2 + (IN_con.a_tri/2)^2)^0.5 else 0)
IN_con.a_tri + 2*sqrt((IN_con.h_tri)^2 + (IN_con.a_tri/2)^2) else 0)
"Perimeter";
SI.Diameter d_hyd=4*A_cross/perimeter "Hydraulic diameter";
Real beta=IN_con.beta*180/PI "Top angle";
Expand Down Expand Up @@ -4548,7 +4548,7 @@ This record is used as <strong>input record</strong> for the pressure loss funct
"Volume flow rate";
SI.Pressure dp_min=max(Modelica.Constants.eps, abs(IN_con.dp_min))
"Start of approximation for decreasing pressure loss";
SI.VolumeFlowRate V_flow_smooth=if a > 0 and b <= 0 then (dp_min/a)^0.5 else 0
SI.VolumeFlowRate V_flow_smooth=if a > 0 and b <= 0 then sqrt(dp_min/a) else 0
"Start of approximation for decreasing volume flow rate";

//Documentation
Expand Down Expand Up @@ -5909,8 +5909,8 @@ Calculation of pressure loss for <strong>two phase flow</strong> in a horizontal
SI.Area Av=if IN_con.valveCoefficient == TYP1.AV then IN_con.Av else if
IN_con.valveCoefficient == TYP1.KV then IN_con.Kv*27.7e-6 else if IN_con.valveCoefficient
== TYP1.CV then IN_con.Cv*24e-6 else if IN_con.valveCoefficient == TYP1.OP then
IN_con.m_flow_nominal/max(MIN, IN_con.opening_nominal*(IN_con.rho_nominal
*IN_con.dp_nominal)^0.5) else MIN
IN_con.m_flow_nominal/max(MIN, IN_con.opening_nominal*sqrt(IN_con.rho_nominal
*IN_con.dp_nominal)) else MIN
"Av (metric) flow coefficient [Av]=m^2";

TYP.PressureLossCoefficient zeta_bal=SMOOTH(
Expand Down Expand Up @@ -5955,12 +5955,12 @@ Calculation of pressure loss for <strong>two phase flow</strong> in a horizontal
== TYP2.Gate then zeta_gat else if IN_con.geometry == TYP2.Sluice then
zeta_slu else 0 "Pressure loss coefficient of chosen valve";

Real valveCharacteristic=(2/min(IN_con.zeta_TOT_max, max(MIN, max(IN_con.zeta_TOT_min,
abs(zeta_TOT)))))^0.5
Real valveCharacteristic=sqrt(2/min(IN_con.zeta_TOT_max, max(MIN, max(IN_con.zeta_TOT_min,
abs(zeta_TOT)))))
"Valve characteristic considering opening of chosen valve";

SI.MassFlowRate m_flow_small=valveCharacteristic*Av*(IN_var.rho)^0.5*(IN_con.dp_small)
^0.5 "Mass flow rate at linearisation";
SI.MassFlowRate m_flow_small=valveCharacteristic*Av*sqrt(IN_var.rho)*sqrt(IN_con.dp_small)
"Mass flow rate at linearisation";

//Documentation

Expand Down Expand Up @@ -6018,8 +6018,8 @@ Generally this function is numerically best used for the <strong>incompressible
SI.Area Av=if IN_con.valveCoefficient == TYP1.AV then IN_con.Av else if
IN_con.valveCoefficient == TYP1.KV then IN_con.Kv*27.7e-6 else if IN_con.valveCoefficient
== TYP1.CV then IN_con.Cv*24e-6 else if IN_con.valveCoefficient == TYP1.OP then
IN_con.m_flow_nominal/max(MIN, IN_con.opening_nominal*(IN_con.rho_nominal
*IN_con.dp_nominal)^0.5) else MIN
IN_con.m_flow_nominal/max(MIN, IN_con.opening_nominal*sqrt(IN_con.rho_nominal
*IN_con.dp_nominal)) else MIN
"Av (metric) flow coefficient [Av]=m^2";

TYP.PressureLossCoefficient zeta_bal=SMOOTH(
Expand Down Expand Up @@ -6064,14 +6064,14 @@ Generally this function is numerically best used for the <strong>incompressible
== TYP2.Gate then zeta_gat else if IN_con.geometry == TYP2.Sluice then
zeta_slu else 0 "Pressure loss coefficient of chosen valve";

Real valveCharacteristic=(2/min(IN_con.zeta_TOT_max, max(MIN, max(IN_con.zeta_TOT_min,
abs(zeta_TOT)))))^0.5
Real valveCharacteristic=sqrt(2/min(IN_con.zeta_TOT_max, max(MIN, max(IN_con.zeta_TOT_min,
abs(zeta_TOT)))))
"Valve characteristic considering opening of chosen valve";

//Documentation

algorithm
M_FLOW := valveCharacteristic*Av*(IN_var.rho)^0.5*
M_FLOW := valveCharacteristic*Av*sqrt(IN_var.rho)*
Modelica.Fluid.Dissipation.Utilities.Functions.General.SmoothPower(
dp,
IN_con.dp_small,
Expand Down Expand Up @@ -6112,8 +6112,8 @@ Generally this function is numerically best used for the <strong>compressible ca
group="Valve", enable= valveCoefficient == 3));
SI.Pressure dp_nominal=1e3 "Nominal pressure loss" annotation (Dialog(group=
"Valve", enable= valveCoefficient == 4));
SI.MassFlowRate m_flow_nominal=opening_nominal*Av*(rho_nominal*dp_nominal)^
0.5 "Nominal mass flow rate" annotation (Dialog(group="Valve", enable=
SI.MassFlowRate m_flow_nominal=opening_nominal*Av*sqrt(rho_nominal*dp_nominal)
"Nominal mass flow rate" annotation (Dialog(group="Valve", enable=
valveCoefficient == 4));
SI.Density rho_nominal=1000 "Nominal inlet density" annotation (Dialog(group=
"Valve", enable= valveCoefficient == 4));
Expand Down
8 changes: 4 additions & 4 deletions Modelica/Media/R134a.mo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,9 +176,9 @@ package R134a "R134a: Medium model for R134a"
sat.cv := f.R_s*(-f.tau*f.tau*f.ftautau);
sat.pt := f.R_s*f.d*(f.delta*(f.fdelta - f.tau*f.fdeltatau));
sat.pd := f.R_s*f.T*(f.delta*(2.0*f.fdelta + f.delta*f.fdeltadelta));
sat.a := abs(f.R_s*f.T*(2*f.delta*f.fdelta + f.delta*f.delta*f.fdeltadelta
sat.a := sqrt(abs(f.R_s*f.T*(2*f.delta*f.fdelta + f.delta*f.delta*f.fdeltadelta
- ((f.delta*f.fdelta - f.delta*f.tau*f.fdeltatau)*(f.delta*f.fdelta -
f.delta*f.tau*f.fdeltatau))/(f.tau*f.tau*f.ftautau)))^0.5;
f.delta*f.tau*f.fdeltatau))/(f.tau*f.tau*f.ftautau))));
sat.kappa := 1/(f.d*f.R_s*f.T*f.delta*(2.0*f.fdelta + f.delta*f.fdeltadelta));
sat.beta := f.R_s*f.d*f.delta*(f.fdelta - f.tau*f.fdeltatau)*sat.kappa;
sat.gamma := sat.a^2/f.R_s/f.T;
Expand Down Expand Up @@ -1591,9 +1591,9 @@ Proceedings of the Joint Meeting of IIR Commissions B1, B2, E1, and E2, Padua, I
// assert(getPhase_ph(state.p, state.h)==1, "Function for velocity of sound is only valid for one-phase regime!");
else
f := f_R134a(state.d, state.T);
a := abs(R134aData.data.R_s*state.T*(2*f.delta*f.fdelta + f.delta*f.delta
a := sqrt(abs(R134aData.data.R_s*state.T*(2*f.delta*f.fdelta + f.delta*f.delta
*f.fdeltadelta - ((f.delta*f.fdelta - f.delta*f.tau*f.fdeltatau)*(f.delta
*f.fdelta - f.delta*f.tau*f.fdeltatau))/(f.tau*f.tau*f.ftautau)))^0.5;
*f.fdelta - f.delta*f.tau*f.fdeltatau))/(f.tau*f.tau*f.ftautau))));
end if;
annotation (Documentation(info="<html>
<p>This function calculates the velocity of sound of R134a from the state record (e.g., use setState_phX function for input). The velocity of sound is modelled by the fundamental equation of state of Tillner-Roth and Baehr (1994).</p>
Expand Down
24 changes: 12 additions & 12 deletions Modelica/Media/package.mo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7398,8 +7398,8 @@ critical pressure.
pro.cp := -pro.R_s*g.tau*g.tau*g.gtautau;
pro.cv := pro.R_s*(-g.tau*g.tau*g.gtautau + (g.gpi - g.tau*g.gtaupi)*(g.gpi
- g.tau*g.gtaupi)/(g.gpipi));
pro.a := abs(g.R_s*g.T*(g.gpi*g.gpi/((g.gpi - g.tau*g.gtaupi)*(g.gpi - g.tau
*g.gtaupi)/(g.tau*g.tau*g.gtautau) - g.gpipi)))^0.5;
pro.a := sqrt(abs(g.R_s*g.T*(g.gpi*g.gpi/((g.gpi - g.tau*g.gtaupi)*(g.gpi - g.tau
*g.gtaupi)/(g.tau*g.tau*g.gtautau) - g.gpipi))));
vt := g.R_s/g.p*(g.pi*g.gpi - g.tau*g.pi*g.gtaupi);
vp := g.R_s*g.T/(g.p*g.p)*g.pi*g.pi*g.gpipi;
pro.kappa := -1/(pro.d*g.p)*pro.cp/(vp*pro.cp + vt*vt*g.T);
Expand Down Expand Up @@ -7461,8 +7461,8 @@ critical pressure.
vp := g.R_s*g.T/(g.p*g.p)*g.pi*g.pi*g.gpipi;
pro.kappa := -1/((g.p/(pro.R_s*g.T*g.pi*g.gpi))*g.p)*pro.cp/(vp*pro.cp + vt
*vt*g.T);
pro.a := abs(g.R_s*g.T*(g.gpi*g.gpi/((g.gpi - g.tau*g.gtaupi)*(g.gpi - g.tau
*g.gtaupi)/(g.tau*g.tau*g.gtautau) - g.gpipi)))^0.5;
pro.a := sqrt(abs(g.R_s*g.T*(g.gpi*g.gpi/((g.gpi - g.tau*g.gtaupi)*(g.gpi - g.tau
*g.gtaupi)/(g.tau*g.tau*g.gtautau) - g.gpipi))));

d := g.p/(pro.R_s*g.T*g.pi*g.gpi);
pro.dudT := (pro.p - g.T*vt/vp)/(d*d);
Expand Down Expand Up @@ -7492,8 +7492,8 @@ critical pressure.
vt := g.R_s/g.p*(g.pi*g.gpi - g.tau*g.pi*g.gtaupi);
vp := g.R_s*g.T/(g.p*g.p)*g.pi*g.pi*g.gpipi;
pro.kappa := -1/(pro.d*g.p)*pro.cp/(vp*pro.cp + vt*vt*g.T);
pro.a := abs(g.R_s*g.T*(g.gpi*g.gpi/((g.gpi - g.tau*g.gtaupi)*(g.gpi - g.tau
*g.gtaupi)/(g.tau*g.tau*g.gtautau) - g.gpipi)))^0.5;
pro.a := sqrt(abs(g.R_s*g.T*(g.gpi*g.gpi/((g.gpi - g.tau*g.gtaupi)*(g.gpi - g.tau
*g.gtaupi)/(g.tau*g.tau*g.gtautau) - g.gpipi))));
pro.ddpT := -(pro.d*pro.d)*vp;
pro.ddTp := -(pro.d*pro.d)*vt;
pro.duTp := pro.cp - g.p*vt;
Expand Down Expand Up @@ -7530,9 +7530,9 @@ critical pressure.
pro.cv := f.R_s*(-f.tau*f.tau*f.ftautau);
pro.kappa := 1/(f.d*f.R_s*f.d*f.T*f.delta*f.fdelta)*((-pv*pro.cv + pt*pt*f.T)
/(pro.cv));
pro.a := abs(f.R_s*f.T*(2*f.delta*f.fdelta + f.delta*f.delta*f.fdeltadelta
pro.a := sqrt(abs(f.R_s*f.T*(2*f.delta*f.fdelta + f.delta*f.delta*f.fdeltadelta
- ((f.delta*f.fdelta - f.delta*f.tau*f.fdeltatau)*(f.delta*f.fdelta -
f.delta*f.tau*f.fdeltatau))/(f.tau*f.tau*f.ftautau)))^0.5;
f.delta*f.tau*f.fdeltatau))/(f.tau*f.tau*f.ftautau))));
pro.ddph := (f.d*(pro.cv*f.d + pt))/(f.d*f.d*pd*pro.cv + f.T*pt*pt);
pro.ddhp := -f.d*f.d*pt/(f.d*f.d*pd*pro.cv + f.T*pt*pt);
pro.duph := -1/pro.d + p/(pro.d*pro.d)*pro.ddph;
Expand Down Expand Up @@ -7575,9 +7575,9 @@ critical pressure.
pro.cv := f.R_s*(-f.tau*f.tau*f.ftautau);
pro.kappa := 1/(f.d*f.R_s*f.d*f.T*f.delta*f.fdelta)*((-pv*pro.cv + pt*pt*f.T)
/(pro.cv));
pro.a := abs(f.R_s*f.T*(2*f.delta*f.fdelta + f.delta*f.delta*f.fdeltadelta
pro.a := sqrt(abs(f.R_s*f.T*(2*f.delta*f.fdelta + f.delta*f.delta*f.fdeltadelta
- ((f.delta*f.fdelta - f.delta*f.tau*f.fdeltatau)*(f.delta*f.fdelta -
f.delta*f.tau*f.fdeltatau))/(f.tau*f.tau*f.ftautau)))^0.5;
f.delta*f.tau*f.fdeltatau))/(f.tau*f.tau*f.ftautau))));
pro.ddTp := -pt/pd;
pro.ddpT := 1/pd;
//problem with units in last two lines
Expand Down Expand Up @@ -7609,9 +7609,9 @@ critical pressure.
*f.fdeltatau)^2/(2*f.delta*f.fdelta + f.delta*f.delta*f.fdeltadelta));
pro.cv := f.R_s*(-f.tau*f.tau*f.ftautau);
pro.kappa := 1/(f.d*pro.p)*((-pv*pro.cv + pt*pt*f.T)/(pro.cv));
pro.a := abs(f.R_s*f.T*(2*f.delta*f.fdelta + f.delta*f.delta*f.fdeltadelta
pro.a := sqrt(abs(f.R_s*f.T*(2*f.delta*f.fdelta + f.delta*f.delta*f.fdeltadelta
- ((f.delta*f.fdelta - f.delta*f.tau*f.fdeltatau)*(f.delta*f.fdelta -
f.delta*f.tau*f.fdeltatau))/(f.tau*f.tau*f.ftautau)))^0.5;
f.delta*f.tau*f.fdeltatau))/(f.tau*f.tau*f.ftautau))));
pro.dudT := (pro.p - f.T*pt)/(f.d*f.d);
end helmholtzToProps_dT;

Expand Down

0 comments on commit 691f3f0

Please sign in to comment.