Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Pump kinematic viscosity/dynamic viscosity #168

Draft
wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

Pump kinematic viscosity/dynamic viscosity #168

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
099ed41
Corrected error dynamic viscosity vs. kinematic viscosity and put an …
RaphaelGebhart Mar 1, 2024
f0bd48e
Resolved UTF8 issue
RaphaelGebhart Mar 1, 2024
eed14b2
Updated package.order
RaphaelGebhart Mar 5, 2024
b1b95a3
Changed mu to nu
RaphaelGebhart Apr 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions ThermofluidStream/Processes/Internal/TurboComponent/dp_tau_centrifugal.mo
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,8 @@ protected
Real gamma(unit="1") = V_flow_D/V_flow_D_ref "Flow scaling factor";

SI.SpecificVolume v_in = 1/max(rho_min, Medium.density(state_in)) "Specific volume at inlet";
SI.SpecificVolume mu_in = Medium.dynamicViscosity(state_in) "Specific volume at inlet";
SI.DynamicViscosity eta_in = Medium.dynamicViscosity(state_in) "Dynamic viscosity at inlet";
SI.KinematicViscosity nu_in = eta_in*v_in "Kinematic viscosity at inlet";
SI.SpecificVolume v_ref = 1/rho_ref;

SI.Power W_t "Technical work going into pump";
Expand Down Expand Up @@ -92,9 +93,9 @@ protected

algorithm
//limit abs(omega) to effectiveley limit the Re_mod to Re_mod_min
omega_hat := max(Re_mod_min/(omega_s^1.5*f_q^0.75)/r^2*mu_in, abs(omega));
omega_hat := max(Re_mod_min/(omega_s^1.5*f_q^0.75)/r^2*nu_in, abs(omega));
V_flow_BEP := K_D*omega_hat;
Re_mod := (omega_hat*r^2)/(mu_in)*(omega_s^1.5*f_q^0.75);
Re_mod := (omega_hat*r^2)/(nu_in)*(omega_s^1.5*f_q^0.75);

// compute corresponding volume flow of water and factors
f_Q := Re_mod^(-6.7/(Re_mod^0.735));
Expand Down
302 changes: 302 additions & 0 deletions ThermofluidStream/PumpTFSExplicit.mo
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,302 @@
within ThermofluidStream;
model PumpTFSExplicit "Inserted PartialTurboComponent and function dp_tau_centrifugal into Pump to have access to all variables"
extends ThermofluidStream.Interfaces.SISOFlow(m_flowStateSelect=StateSelect.prefer, final clip_p_out=true);

import ThermofluidStream.Utilities.Types.InitializationMethods;
import Quantity = ThermofluidStream.Sensors.Internal.Types.MassFlowQuantities;

parameter Real max_rel_volume(min=0, max=1, unit="1") = 0.05 "Maximum relative volume change"
annotation(Dialog(tab="Advanced"));

Real eta(unit="1") = if noEvent(abs(W_t)>1e-4) then dp*v_in*m_flow/W_t else 0.0;

protected
SI.SpecificVolume v_out = 1/max(rho_min, Medium.density(inlet.state));

public
parameter Boolean omega_from_input = false "Direct omega input"
annotation(Dialog(group="Input/Output"));
parameter Boolean enableOutput = false "Include output for selectable quantity"
annotation(Dialog(group="Input/Output"));
parameter Quantity outputQuantity=Quantity.m_flow_kgps "Quantity to output"
annotation(Dialog(group="Input/Output", enable=enableOutput));
parameter Boolean enableAccessHeatPort = false "Include access heatport"
annotation(Dialog(group="Input/Output"));
parameter SI.MomentOfInertia J_p = 5e-4 "Moment of inertia"
annotation(Dialog(group="Parameters", enable=not omega_from_input));
parameter SI.AngularVelocity omega_reg = dropOfCommons.omega_reg "Angular velocity used for normalization"
annotation(Dialog(tab="Advanced", group="Regularization"));
parameter SI.MassFlowRate m_flow_reg = dropOfCommons.m_flow_reg "Mass flow used for normalization"
annotation(Dialog(tab="Advanced", group="Regularization"));
parameter SI.Density rho_min = dropOfCommons.rho_min "Minimum for rho (to make model stable for rho=0 @ p=0)"
annotation(Dialog(tab="Advanced", group="Regularization"));
parameter StateSelect omegaStateSelect = if omega_from_input then StateSelect.default else StateSelect.prefer "State select for m_flow"
annotation(Dialog(tab="Advanced"));
parameter InitializationMethods initOmega = ThermofluidStream.Utilities.Types.InitializationMethods.none "Initialization method for omega"
annotation(Dialog(tab= "Initialization", group="Angular", enable=not omega_from_input));
parameter SI.AngularVelocity omega_0 = 0 "Initial value for omega"
annotation(Dialog(tab= "Initialization", group="Angular", enable=(initOmega == InitializationMethods.state)));
parameter SI.AngularAcceleration omega_dot_0 = 0 "Initial value for der(omega)"
annotation(Dialog(tab= "Initialization", group="Angular", enable=(initOmega == InitializationMethods.derivative)));
parameter Boolean initPhi = false "If true phi is initialized"
annotation(Dialog(tab= "Initialization", group="Angular", enable=not omega_from_input));
parameter SI.Angle phi_0 = 0 "Initial value for phi"
annotation(Dialog(tab= "Initialization", group="Angular", enable=initPhi));

Modelica.Mechanics.Rotational.Interfaces.Flange_a flange(phi=phi, tau=tau) if not omega_from_input "Flange to receive mechanical power"
annotation (Placement(transformation(extent={{-20,-20},{20,20}}, origin={0,-100}, rotation=-90),
iconTransformation(extent={{-20,-20},{20,20}}, origin={0,-100}, rotation=-90)));
Modelica.Blocks.Interfaces.RealInput omega_input(unit = "rad/s") = omega if omega_from_input "Input to directly set pump speed [rad/s]"
annotation (Placement(transformation(extent={{-20,-20},{20,20}}, origin={0,-100}, rotation=-90),
iconTransformation(extent={{-20,-20},{20,20}}, origin={0,-100}, rotation=90)));
Modelica.Thermal.HeatTransfer.Interfaces.HeatPort_a heatport(Q_flow = Q_t) if enableAccessHeatPort "Access-heat dumping port"
annotation (Placement(transformation(extent={{-20,-20},{20,20}}, origin={0,100}, rotation=90),
iconTransformation(extent={{-20,-20},{20,20}}, origin={0,100}, rotation=90)));
Modelica.Blocks.Interfaces.RealOutput output_val(unit=ThermofluidStream.Sensors.Internal.getFlowUnit(
outputQuantity)) = getQuantity(inlet.state, m_flow, outputQuantity, rho_min) if enableOutput "Measured value [variable]"
annotation (Placement(transformation(extent={{-20,-20},{20,20}}, origin={-80,-100}, rotation=270)));

function getQuantity = ThermofluidStream.Sensors.Internal.getFlowQuantity (redeclare package Medium = Medium)
"Function to compute a selectable quantity"
annotation (
Documentation(info="<html>
<p>Function to compute a selectable quantity to output. The quantity is associated to the mass flow. </p>
</html>"));

SI.Power W_t "technichal work performed on fluid";
SI.Power Q_t "work that could not be performed on the fluid, and is dumped to heat port";
SI.Torque tau_st "moment that would result in static operation";

protected
SI.SpecificVolume v_in = 1/max(rho_min, Medium.density(inlet.state));

// intermediate quantities related to thermodynamical state
SI.SpecificEnergy dh "specific technichal work performed under the assumption of positive massflow";

// quantities related to inertia of component and mechanical power
SI.Angle phi "component angular position";
SI.Torque tau "moment from flange";
SI.Torque tau_normalized "moment after zero massflow normalization";
SI.AngularVelocity omega(stateSelect=omegaStateSelect) "component angular velocity";

// extends ThermofluidStream.Processes.Internal.TurboComponent.partial_dp_tau;

public
parameter Boolean parametrizeByScaling = true "If true scale by three parameters"
annotation(Dialog(enable=true));
parameter SI.Height TDH_D = 3.6610 "Design pressure head"
annotation(Dialog(group="Scaling", enable=parametrizeByScaling));
parameter SI.VolumeFlowRate V_flow_D = 3.06e-3 "Design volume flow"
annotation(Dialog(group="Scaling", enable=parametrizeByScaling));
parameter SI.AngularVelocity omega_D = 314.2 "Design angular velocity"
annotation(Dialog(group="Scaling", enable=parametrizeByScaling));
parameter Real K_D_input(unit="m3/rad") = 9.73e-06 "Vflow_D / omega_D"
annotation(Dialog(group="General", enable=not parametrizeByScaling));
parameter Integer f_q_input = 1 "Number of floods"
annotation(Dialog(group="General", enable=not parametrizeByScaling));
parameter SI.Radius r_input = 1.60e-2 "Radius of pump"
annotation(Dialog(group="General", enable=not parametrizeByScaling));
parameter SI.Density rho_ref_input = 1.00e3 "Reference Density"
annotation(Dialog(group="General", enable=not parametrizeByScaling));
parameter Real a_h_input(unit= "m.s2/rad2", displayUnit="m.min2") = 4.864e-5 "HQ factor 1"
annotation(Dialog(group="HQ characteristic", enable=not parametrizeByScaling));
parameter Real b_h_input(unit="s2/(m2.rad)", displayUnit="m.h.min/l") = -2.677 "HQ factor 2"
annotation(Dialog(group="HQ characteristic", enable=not parametrizeByScaling));
parameter Real c_h_input(unit="s2/m5", displayUnit="m.h2/l2") = 3.967e+5 "HQ factor 3"
annotation(Dialog(group="HQ characteristic", enable=not parametrizeByScaling));
parameter Real a_t_input(unit="N.m.s2/(rad.m3)", displayUnit="N.m.min.h/l") = 5.427e-1 "TQ factor 1"
annotation(Dialog(group="TQ characteristic", enable=not parametrizeByScaling));
parameter Real b_t_input(unit="N.m.s2/m6", displayUnit="N.m.h2/l2") = 2.777e+4 "TQ factor 2"
annotation(Dialog(group="TQ characteristic", enable=not parametrizeByScaling));
parameter Real v_i_input(unit="N.m.s2/rad2", displayUnit="N.m.min2") = 1.218e-6 "TQ factor 4"
annotation(Dialog(group="TQ characteristic", enable=not parametrizeByScaling));
parameter Real v_s_input(unit="N.m.s/rad", displayUnit="N.m.min") = 1.832e-4 "TQ factor 3"
annotation(Dialog(group="TQ characteristic", enable=not parametrizeByScaling));
parameter Real Re_mod_min(unit="1") = 1e2 "Minimum modified Reynolds number"
annotation(Dialog(tab="Advanced", enable=true));

protected
Medium.ThermodynamicState state_in = inlet.state "Thermodynamic state at inlet";
Medium.MassFlowRate m_flow_norm = m_flow_reg "Mass flow used for normalization";
SI.AngularVelocity omega_norm = omega_reg "Angular velocity used for normalization";

Real alpha(unit="1") = TDH_D/TDH_D_ref "Pressure scaling factor";
Real beta(unit="1") = omega_D/omega_D_ref "Speed scaling factor";
Real gamma(unit="1") = V_flow_D/V_flow_D_ref "Flow scaling factor";

SI.Density rho_in = max(rho_min, Medium.density(state_in)) "Density at inlet";
SI.DynamicViscosity eta_in = Medium.dynamicViscosity(state_in) "Dynamic viscosity at inlet";
SI.KinematicViscosity mu_in = eta_in/rho_in "Kinematic viscosity at inlet";

// SI.Power W_t_f "Technical work going into pump";
SI.Length TDH "Total dynamic head";
SI.VolumeFlowRate V_flow "Volume flow through pump";

constant SI.AngularVelocity omega_D_ref = 314.2;
constant SI.VolumeFlowRate V_flow_D_ref = 3.06e-3;
constant SI.Height TDH_D_ref = 3.6610;

constant Real a_h_ref(unit="m.s2/rad2") = 4.864e-5;
constant Real b_h_ref(unit="s2/(m2.rad)") = -2.677;
constant Real c_h_ref(unit="s2/m5") = 3.967e+5;

constant Real a_t_ref(unit="N.m.s2/(rad.m3)") = 5.427e-1;
constant Real b_t_ref(unit="N.m.s2/m6") = 2.777e+4;
constant Real v_i_ref(unit="N.m.s2/rad2") = 1.218e-6;
constant Real v_s_ref(unit="N.m.s/rad") = 1.832e-4;

constant Real f_q_ref = 1;
constant Real K_D_ref(unit="m3/rad") = 9.73e-06;
constant SI.Density rho_ref_ref = 1.00e3;
constant SI.Length r_ref = 1.60e-2;

Real a_h(unit="m.s2/rad2") = if parametrizeByScaling then alpha/(beta^2)* a_h_ref else a_h_input;
Real b_h(unit="s2/(m2.rad)") = if parametrizeByScaling then alpha/(beta*gamma)* b_h_ref else b_h_input;
Real c_h(unit="s2/m5") = if parametrizeByScaling then alpha/(gamma^2)* c_h_ref else c_h_input;

Real a_t(unit="N.m.s2/(rad.m3)") = if parametrizeByScaling then alpha/(beta^2)* a_t_ref else a_t_input;
Real b_t(unit="N.m.s2/m6") = if parametrizeByScaling then alpha/(beta*gamma)* b_t_ref else b_t_input;
Real v_i(unit="N.m.s2/rad2") = if parametrizeByScaling then sqrt(alpha)*gamma/beta* v_i_ref else v_i_input;
Real v_s(unit="N.m.s/rad") = if parametrizeByScaling then alpha^(-1/4)*beta^(1/2)*gamma^(3/2)* v_s_ref else v_s_input;

Real f_q = if parametrizeByScaling then f_q_ref else f_q_input;
Real K_D(unit="m3/rad") = if parametrizeByScaling then gamma/beta* K_D_ref else f_q_input;
SI.Density rho_ref = if parametrizeByScaling then rho_ref_ref else rho_ref_input;
SI.Length r = if parametrizeByScaling then sqrt(alpha)/beta* r_ref else r_input;

Real omega_s(unit="1") = ((K_D/f_q)^0.5)/((Modelica.Constants.g_n*(a_h+b_h*K_D+c_h*K_D^2)) ^0.75);
SI.VolumeFlowRate V_flow_BEP "optimal volume flow at omega";
SI.AngularVelocity omega_hat "abs omega clipped Re_mod_min";
Real Re_mod(unit="1");
Real f_Q(unit="1");
Real f_H(unit="1");
Real f_eta(unit="1");

initial equation
if initOmega == InitializationMethods.state then
omega = omega_0;
elseif initM_flow == InitializationMethods.derivative then
der(omega) = omega_dot_0;
elseif initM_flow == InitializationMethods.steadyState then
der(omega) = 0;
end if;

if initPhi then
phi = phi_0;
end if;

equation

algorithm
omega_hat := max(Re_mod_min/(omega_s^1.5*f_q^0.75)/r^2*mu_in, abs(omega));
//omega_hat := abs(omega);
V_flow_BEP := K_D*omega_hat;
Re_mod := (omega_hat*r^2)/(mu_in)*(omega_s^1.5*f_q^0.75);

// compute corresponding volume flow of water and factors

f_Q := Re_mod^(-6.7/(Re_mod^0.735));
f_eta := Re_mod^(-19/(Re_mod^0.705));
V_flow := m_flow/f_Q/rho_ref;
f_H := 1 - (1 - f_Q)*abs(V_flow/V_flow_BEP)^0.75;

/*
f_Q := 1;
f_eta := 1;
V_flow := m_flow/f_Q/rho_ref;
f_H := 1;
*/
// corrected characteristic curves for TDH and tau_st
TDH := f_H *(a_h*omega*abs(omega) - b_h*omega*abs(V_flow) - c_h*V_flow*abs(V_flow));

// original version
// tau_st := (f_Q*f_H/f_eta) * (rho_in/rho_ref*(a_t*V_flow*abs(omega) - b_t*V_flow*abs(V_flow) + v_i*omega*abs(omega))) + v_s*abs(omega); // v_s is mechanical and does not scale with medium
// docu version
// tau_st := (f_Q*f_H/f_eta)*(rho_in/rho_ref*a_t*abs(omega)*V_flow - rho_in/rho_ref*b_t*abs(V_flow)*V_flow + v_i*abs(omega)*omega + v_s*abs(omega));

// there could be several versions....
tau_st := (f_Q*f_H/f_eta)*rho_in/rho_ref*((a_t*V_flow*abs(omega) - b_t*V_flow*abs(V_flow) + v_i*omega*abs(omega)) + v_s*abs(omega)); // v_s is mechanical and does not scale with medium

//compute dp
dp := rho_in*Modelica.Constants.g_n*TDH;
equation
// compute dp, tau_st from characteristic curve
//(dp, tau_st) = dp_tau_pump(m_flow, omega, inlet.state, m_flow_reg, omega_reg, rho_min);

h_out = h_in + dh;
Xi_out = Xi_in;

// regularize tau_st/W_t such that the equations dont break for m_flow->0
W_t = tau_st*omega;
dh = (m_flow*W_t)/(m_flow^2 + (m_flow_reg)^2);
if noEvent(W_t >= 0) then
// if work is given to fluid, dump access heat to port
Q_t = W_t - m_flow*dh;
tau_normalized = tau_st;
else
// if work is taken from fluid, reduce tau, so that the work can be taken from the fluid
Q_t = 0;
tau_normalized = dh*m_flow/(noEvent(if abs(omega)>omega_reg then omega else (if omega < 0 then -omega_reg else omega_reg)));
end if;

// omega is either a state or directly given by omega_input (depending on omega_from_input)
if noEvent(omega_from_input) then
tau = 0;
phi = 0;
else
J_p * der(omega) = tau - tau_normalized;
omega = der(phi);
end if;

// Pump test for incompressibility
assert(abs(v_in- v_out)/v_in < max_rel_volume, "Medium in pump is assumed to be incompressible, but check failed", dropOfCommons.assertionLevel);

annotation (Icon(coordinateSystem(preserveAspectRatio=false), graphics={
Text(visible=displayComponentName,
extent={{-150,140},{150,100}},
textString="%name",
textColor={0,0,255}),
Ellipse(
extent={{-56,54},{64,-66}},
lineColor={28,108,200},
lineThickness=0.5,
fillColor={215,215,215},
fillPattern=FillPattern.Solid,
pattern=LinePattern.None),
Line(
points={{-100,0},{100,0}},
color={28,108,200},
thickness=0.5),
Ellipse(
extent={{-60,60},{60,-60}},
lineColor={28,108,200},
lineThickness=0.5,
fillColor={255,255,255},
fillPattern=FillPattern.Solid),
Line(
points={{0,60},{60,0}},
color={28,108,200},
thickness=0.5),
Line(
points={{0,-60},{60,0}},
color={28,108,200},
thickness=0.5),
Rectangle(
extent={{-60,20},{66,-20}},
lineColor={28,108,200},
fillColor={255,255,255},
fillPattern=FillPattern.Solid),
Text(
extent={{-52,16},{60,-14}},
textColor={28,108,200},
textString="test version")}),
Diagram(coordinateSystem(preserveAspectRatio=false)),
Documentation(info="<html>
<p>This model is&nbsp;working&nbsp;under&nbsp;the&nbsp;assumption&nbsp;of&nbsp;incompressible&nbsp;fluids and performs an assert for this. </p>
<p>It can be chosen between</p>
<ul>
<li>Nominal-flow pump, where a nominal flow is computed and the difference between it and the actual flow is linearly producing a pressure</li>
<li>Centrifugal pump, which implements the equations of a scalable centrifugal pump.</li>
</ul>
</html>"));
end PumpTFSExplicit;
Loading
Loading