diff --git a/.gitignore b/.gitignore index b471c606..04c96837 100644 --- a/.gitignore +++ b/.gitignore @@ -28,7 +28,6 @@ simulator.log **/*.bbl **/*.blg **/*.log -**/*.pdf **/*.synctex.gz **/*.toc diff --git a/OpenHPL/Controllers/Governor.mo b/OpenHPL/Controllers/Governor.mo index 458d6503..2bc4dc07 100644 --- a/OpenHPL/Controllers/Governor.mo +++ b/OpenHPL/Controllers/Governor.mo @@ -126,10 +126,36 @@ equation //else // der(Y_gv) = u / T_g; //end if; - annotation ( - Documentation(info=" -
This is a simple model of the governor that controls the guide vane - opening in the turbine based on the reference power production.
+ annotation (preferredView="info", Documentation(info=" ++Here, a simple model of the governor that controls the guide vane opening in the turbine based on the reference power +production is described. The block diagram of this governor model is shown in the figure. +
+ +
+
+
Figure: Block Diagram of the governor.
+ ++Using the model in the figure and the standard Modelica blocks, the governor model is encoded in our library as the +Governor unit. This unit has inputs as the reference power production and generator frequency that are implemented +with the standard Modelica RealInput connector. This Governor unit also uses the standard Modelica +RealOutput connectors in order to provide output information about the turbine guide vane opening. +
+ ++In the Governor unit (note: in the text it mentions SynchGen but this appears to be a typo in the +original document - should be Governor), the user can specify the various time constants of this model (see +figure): pilot servomotor time constant Tp, primary servomotor integration time Tg, and transient +droop time constant Tr. The user should also provide the following parameters: droop value σ, transient droop δ, +and nominal values for the frequency and power generation. The information about the maximum, minimum, and initial guide +vane opening should also be specified. +
+The model is taken from [Sharefi2011].
")); end Governor; diff --git a/OpenHPL/Controllers/package.mo b/OpenHPL/Controllers/package.mo index a0848637..9d419049 100644 --- a/OpenHPL/Controllers/package.mo +++ b/OpenHPL/Controllers/package.mo @@ -1,7 +1,6 @@ within OpenHPL; package Controllers "Collection of different controllers" extends Modelica.Icons.Package; - extends Icons.Governor; end Controllers; diff --git a/OpenHPL/ElectroMech/BaseClasses/BaseValve.mo b/OpenHPL/ElectroMech/BaseClasses/BaseValve.mo index 389a0b21..1d713837 100644 --- a/OpenHPL/ElectroMech/BaseClasses/BaseValve.mo +++ b/OpenHPL/ElectroMech/BaseClasses/BaseValve.mo @@ -34,28 +34,24 @@ equation Vdot = mdot/data.rho; dp*(C_v_*max(epsilon, u^alpha))^2 = Vdot*abs(Vdot) "Valve equation for pressure drop"; dp = i.p - o.p "Link the pressure drop to the ports"; - annotation ( Documentation(info=" + annotation (preferredView="info", Documentation(info="This is a partial, simple model of hydraulic valve.
This model is based on the energy balance of a valve.
Specifically:
--dp*f(opening)=v|v| -+
$$ \\Delta p \\cdot f(\\mathrm{opening}) = \\nu \\cdot | \\nu | $$
The function f(opening) is expressed as:
--(C_v_*max(epsilon, u^alpha))^2 -+
$$ f(\\mathrm{opening}) = \\left( C_\\mathrm{v} \\cdot \\mathrm{max}(\\epsilon, u^\\alpha)\\right)^2 $$
-When alpha is 1, this implies a linear relation between closing and head loss. +When \\(\\alpha\\) is 1, this implies a linear relation between closing and head loss.
The valve capacity can either be specified
directly by the user by specifying C_v or it will be calculated from
diff --git a/OpenHPL/ElectroMech/Generators/SimpleGen.mo b/OpenHPL/ElectroMech/Generators/SimpleGen.mo
index a16d06d8..a54648b7 100644
--- a/OpenHPL/ElectroMech/Generators/SimpleGen.mo
+++ b/OpenHPL/ElectroMech/Generators/SimpleGen.mo
@@ -13,20 +13,51 @@ model SimpleGen "Model of a simple generator with mechanical connectors"
rotation=270,
origin={0,120})));
- annotation (
+ annotation (preferredView="info",
Documentation(info= "
-
Simple model of an ideal generator with friction based on angular momentum balance.
-This model based on the angular momentum balance, which depends on the turbine shaft power, - the friction loss in the unit rotation and the power taken up by the generator.
--The generator can be loaded either:
+The kinetic energy stored in the rotating generator is \\(K_a = \\frac{1}{2}J_a\\omega_a^2\\), +where ωa is angular velocity and Ja is moment of inertia.
+ +From energy balance:
+$$ \\frac{\\mathrm{d}K_a}{\\mathrm{d}t} = \\dot{W}_s - \\dot{W}_{f,a} - \\dot{W}_g $$
+where:
+Frictional power loss (mainly from bearings):
+$$ \\dot{W}_{f,a} = \\frac{1}{2}k_{f,b}\\omega_a^2 $$
+where kf,b is the bearing friction factor.
+ +Electric power available on grid:
+$$ \\dot{W}_e = \\eta_e \\dot{W}_g $$
+where ηe is electrical efficiency.
+ +The generator can be loaded either:
Pload should be set to 0 in this case.Pload input to 0 in this case.Pload specifying the connected electrical load.Pload) and shaft powerUser specifies: moment of inertia, electrical efficiency, bearing friction factor, number of poles, and initial angular velocity.
+
This is a model of the generator that is connected to the grid. - This model could give some transient results. However, it is better to use generator models from IPSL.
-More info about this model can be found in [Sharefi2011].
- ")); + annotation (preferredView="info", + Documentation(info= " +Detailed synchronous generator model connected to the grid, based on d-q decomposition.
+ +$$ \\left[\\begin{matrix}R_a+R_e & x_q'+x_e\\\\ -x_d'-x_e & R_a+R_e\\end{matrix}\\right]\\left[\\begin{matrix}I_d \\\\ I_q\\end{matrix}\\right]= \\left[\\begin{matrix}E_d'+V_s\\sin\\delta_e \\\\ E_q'-V_s\\cos\\delta_e\\end{matrix}\\right] $$
+where:
+$$ \\frac{\\mathrm{d}\\delta_e}{\\mathrm{d}t} = (\\omega - \\omega_s)\\frac{n_p}{2} $$
+where \\(n_p\\) is number of poles, \\(\\omega\\) and \\(\\omega_s\\) are generator and grid angular velocities.
+ +$$ \\frac{\\mathrm{d}\\omega}{\\mathrm{d}t}=\\frac{\\dot{W}_s-P_e}{J\\omega} $$
+ +$$ +\\begin{array}{c} +T_{qo}'\\frac{\\mathrm{d}E_d'}{\\mathrm{d}t} =-E_d' + (x_q' - x_q)I_q \\\\ +T_{do}'\\frac{\\mathrm{d}E_q'}{\\mathrm{d}t} = -E_q' + (x_d - x_d')I_d + E_f +\\end{array} +$$
+where \\(T_{do}'\\) and \\(T_{qo}'\\) are d-/q-axis transient open-circuit time constants.
+ +Field voltage dynamics:
+$$ \\frac{\\mathrm{d}E_f}{\\mathrm{d}t} = \\frac{-E_f + K_E\\left(V_{tr}-V_t-V_{stab}\\right)}{T_E} $$
+where \\(K_E\\) is excitation system gain, \\(T_E\\) is excitation time constant, \\(V_{tr}\\) is voltage reference set point, +and \\(V_t = \\sqrt{\\left(E_d'-R_aI_d-x_q'I_q\\right)^2+\\left(E_q'-R_aI_q+x_d'I_d\\right)^2}\\) is terminal voltage.
+ +$$ \\frac{\\mathrm{d}V_{stab}}{\\mathrm{d}t} = \\frac{-V_{stab} + K_F\\frac{\\mathrm{d}E_f}{\\mathrm{d}t}}{T_{FE}} $$
+where \\(K_F\\) is stabilizer gain and \\(T_{FE}\\) is stabilizer time constant.
+ +Active and reactive power:
+$$ +\\begin{array}{c} +P_e = 3\\left(E_d'I_d+E_q'I_q\\right)\\\\ +Q_e = \\sqrt{9V_t^2I_t^2-P_e^2} +\\end{array} +$$
+where \\(I_t=\\sqrt{I_d^2+I_q^2}\\) is terminate current.
+ +User specifies: nominal active/reactive powers, phase winding resistance, number of poles, network parameters +(equivalent resistance/reactance, RMS voltage, grid angular velocity), d-/q-axis reactances and time constants, +field voltage limits, excitation/stabilizer gains and time constants, moment of inertia, friction factor, and +initialization options.
+ +Note: For more advanced modeling, consider using generator models from OpenIPSL.
+More details in [Sharefi2011].
+")); end SynchGen; diff --git a/OpenHPL/ElectroMech/PowerSystem/Grid.mo b/OpenHPL/ElectroMech/PowerSystem/Grid.mo index 6e71a359..210c57df 100644 --- a/OpenHPL/ElectroMech/PowerSystem/Grid.mo +++ b/OpenHPL/ElectroMech/PowerSystem/Grid.mo @@ -49,7 +49,7 @@ equation connect(product.y, dP2.u1) annotation (Line(points={{-13.5,67},{-8,67},{-8,33.6},{-12.8,33.6}}, color={0,0,127})); connect(toHz.u, w_m2pu.y) annotation (Line(points={{5.2,-50},{84,-50},{84,-40},{78.6,-40}}, color={0,0,127})); connect(toHz.y, dF.u1) annotation (Line(points={{-8.6,-50},{-22,-50}}, color={0,0,127})); - annotation ( + annotation (preferredView="info", Documentation(info="-This is the Francis turbine model that gives possibilities for proper modelling of the Francis turbine. -
-The mechanistic model is based on Euler equations for the Francis turbine. -Besides hydraulic input and output, there are input as the control signal for the valve opening -and also output as the turbine shaft power and input as angular velocity. -
+This is the Francis turbine model that gives possibilities for proper modelling of the Francis turbine. +The mechanistic model is based on Euler equations for the Francis turbine.
+ +Besides hydraulic input and output, there are input as the control signal for the valve opening +and also output as the turbine shaft power and input as angular velocity.
+
Figure: Key quantities in the Francis turbine model showing inlet (1) and outlet (2) sections.
+ +The shaft power produced by the turbine is given by:
+$$ \\dot{W}_s = \\dot{m}\\omega \\left(R_1\\frac{\\dot{V}}{A_1}\\cot{\\alpha_1} - R_2\\left(\\omega R_2 + \\frac{\\dot{V}}{A_2}\\cot{\\beta_2}\\right)\\right) $$
+ +where ṁ and V̇ are mass and volumetric flow rates, ω is angular velocity, R1 and R2 +are inlet and outlet radii, A1 and A2 are cross-sectional areas, α1 is +inlet guide vane angle, and β2 is outlet blade angle.
+ +The total work rate is:
+$$ \\dot{W}_t = \\dot{W}_s + \\dot{W}_{ft} + \\Delta p_v \\dot{V} $$
+where Ẇft represents various friction losses (shock, whirl, wall friction), +and ΔpvV̇ accounts for guide vane pressure drop. Turbine efficiency \\(\\eta = \\dot{W}_s / \\dot{W}_t\\).
+ +There is also available the runner design algorithm that can define all geometrical -parameters based on the nominal parameters.
The turbine losses coefficients
-(k_ft1, k_ft2, k_ft3) can be also defined automatically.
-However, if some dynamic data from real turbine is available it is better to tune
-these parameters a bit more and use the defined values as a starting point.
-
A model for servo that that runs the guide vane opening is also available. -Furthermore it is possible to automatically generate all need parameters for the servo, - or simply specify them. -
--This mechanistic turbine model does not work really well for low loads (<10% guide vane opening). +parameters based on the nominal parameters (net head, flow rate, power, speed). +The algorithm determines: outlet blade angle β2, runner radii R1 and R2, +runner width w1, and inlet blade angle β1.
+ +The turbine losses coefficients (k_ft1, k_ft2, k_ft3) can be also
+defined automatically. However, if some dynamic data from real turbine is available it is better to tune
+these parameters a bit more and use the defined values as a starting point.
A model for servo that runs the guide vane opening is also available. A guide vane opening model relates +actuator position Y to guide vane angle α1 through geometric relationships. Furthermore it is +possible to automatically generate all needed parameters for the servo, or simply specify them.
+ +This mechanistic turbine model does not work really well for low loads (<10% guide vane opening).
However there is parameters that could be tuned for low load regimes.
These are u_min and k_ft4.
More info about the mechanistic turbine model can be found in: -[Vytvytskyi2018] and about -the servo (also turbine model) in: +[Vytvytskyi2018] and +[Vytvytskyi2019]. +Additional details about the servo (and turbine model) are in: Resources/Documents/Turbines_model.pdf.
")); end Francis; diff --git a/OpenHPL/ElectroMech/Turbines/Pelton.mo b/OpenHPL/ElectroMech/Turbines/Pelton.mo index 57870521..0b019562 100644 --- a/OpenHPL/ElectroMech/Turbines/Pelton.mo +++ b/OpenHPL/ElectroMech/Turbines/Pelton.mo @@ -56,18 +56,53 @@ equation mdot=i.mdot; connect(p_out, P_out) annotation (Line(points={{40,90},{40,110}}, color={0,0,127})); - annotation ( + annotation (preferredView="info", Documentation(info=" -This is a model of the Pelton turbine. -This model is based on the Euler turbine equation. -
--The model has not been tested.
+Mechanistic Pelton turbine model based on the Euler turbine equation and impulse turbine principles.
+
-
-
More info about the model can be found in:
-Resources/Documents/Turbines_model.pdf
+
Figure: Key concepts of the Pelton turbine.
+ +The shaft power \\(\\dot{W}_s\\) produced in the Pelton turbine is:
+$$ \\dot{W}_s=\\dot{m}v_R\\left[\\delta(u_\\delta)\\cdot v_1-v_R\\right]\\left(1-k\\cos\\beta\\right) $$
+where:
+Total work rate removed through the turbine:
+$$ {\\dot{W}_t} = {\\dot{W}_s+\\dot{W}_{ft}} $$
+Friction losses:
+$$ \\dot{W}_{ft}=K\\left(1-k\\cos\\beta\\right)\\dot{m}v_R^2 $$
+with friction coefficient \\(K=0.25\\).
+ +Pressure drop across the nozzle (positions \"0\" and \"1\"):
+$$ \\Delta p_n=\\frac{1}{2}\\rho\\dot{V}\\left[\\dot{V}\\left(\\frac{1}{A_1^2(Y)}-\\frac{1}{A_0^2}\\right)+k_f\\right] $$
+where \\(A_0\\) is cross-sectional area at nozzle beginning, \\(A_1(Y)\\) is area at nozzle end (function of needle position Y), +and \\(k_f\\) is the nozzle friction loss coefficient.
+ +User specifies: turbine runner radius, nozzle input diameter, runner bucket angle, friction factors and coefficients, +deflector mechanism coefficient.
+ +Note: This model has not been tested.
+More info in: Resources/Documents/Turbines_model.pdf
")); end Pelton; diff --git a/OpenHPL/ElectroMech/Turbines/Turbine.mo b/OpenHPL/ElectroMech/Turbines/Turbine.mo index 690df2c5..fca25597 100644 --- a/OpenHPL/ElectroMech/Turbines/Turbine.mo +++ b/OpenHPL/ElectroMech/Turbines/Turbine.mo @@ -43,35 +43,49 @@ equation connect(lossCorrection.u1, power.y) annotation (Line(points={{-48,80},{-88,80},{-88,30},{-81,30}},color={0,0,127})); connect(frictionLoss.power, lossCorrection.u2) annotation (Line(points={{-1,12},{-40,12},{-40,72}},color={0,0,127})); connect(u_t, u) annotation (Line(points={{-80,120},{-80,90},{0,90},{0,70}},color={0,0,127})); - annotation ( + annotation (preferredView="info", Documentation(info=" --This is a simple model of the turbine that give possibilities for simplified -modelling of the turbine unit. The model can use a constant efficiency or varying -efficiency from a lookup-table. -This model does not include any information about rotational speed of the runner. -
-
-This model is based on the energy balance and a simple valve-like expression.
-The guide vane 'valve capacity' should be used for this valve-like expression and can either be specified
-directly by the user by specifying C_v or it will be calculated from
-the Nominal turbinenet head H_n and nominal flow rate
-Vdot_n.
-
-The turbine efficiency is in per-unit values from 0 to 1, where 1 means that there are no losses in the turbine. -The output mechanical power is defined as multiplication of the turbine efficiency and the total possible power: -
---turbine_pressure_drop * turbine_flow_rate-
Besides hydraulic input and output, -there are inputs as the control signal for the valve opening and also output as the turbine shaft power. -
+This is a simple model of the turbine that gives possibilities for simplified +modelling of the turbine unit. The model is based on a look-up table for turbine efficiency +vs. guide vane opening.
-
More info about the model can be found in: Resources/Report/Report.docx
+ +Figure: Simple turbine model schematic.
+ +The mechanical turbine shaft power is:
+$$ \\dot{W}_\\mathrm{tr} = \\eta_\\mathrm{h} \\Delta p_\\mathrm{tr} \\dot{V}_\\mathrm{tr} $$
+where ηh is hydraulic efficiency (from lookup table or constant value), +Δptr is pressure drop, and V̇tr is volumetric flow rate.
+ +The turbine efficiency is in per-unit values from 0 to 1, where 1 means that there are +no losses. The model can use a constant efficiency or varying efficiency from a lookup-table.
+ +This model is based on the energy balance and a simple valve-like expression. +The turbine flow rate relates to pressure drop through:
+$$ \\dot{V}_\\mathrm{tr} = C_\\mathrm{v} u_\\mathrm{v} \\sqrt{\\frac{\\Delta p_\\mathrm{tr}}{p^\\mathrm{atm}}} $$
+where Cv is the guide vane \"valve capacity\", uv is guide vane opening +signal (0 to 1), and patm is atmospheric pressure.
+ +The guide vane valve capacity can either be specified directly by the user via C_v
+or it will be calculated from the nominal turbine net head H_n and nominal flow rate
+Vdot_n.
Besides hydraulic input and output, there are inputs as the control signal for the valve +opening and also output as the turbine shaft power.
+ +More info about the model can be found in: Resources/Report/Report.docx +and [Vytvytskyi2019].
"), Icon(graphics={Text( visible=enable_P_out, extent={{30,100},{50,80}}, diff --git a/OpenHPL/Functions/DarcyFriction/package.mo b/OpenHPL/Functions/DarcyFriction/package.mo index 0c874cf1..ca90306a 100644 --- a/OpenHPL/Functions/DarcyFriction/package.mo +++ b/OpenHPL/Functions/DarcyFriction/package.mo @@ -4,6 +4,100 @@ package DarcyFriction "Functions to define the Darcy friction factor and frictio annotation ( Documentation(info = " -Functions to define the Darcy friction factor and friction Force, which is based on this coefficient.
++First, the functions for defining the friction force in the waterway are described. The friction force +Ff is directed in the opposite direction of the velocity v (the linear velocity average across +the cross-section of the pipe) of the fluid. A common expression for friction force in the filled pipes +is the following: +
++$$ F_\\mathrm{f} = -\\frac{1}{8}\\pi\\rho LDf_\\mathrm{D}v|v| $$ +
++Here, L and D are related to the pipe length and diameter, respectively. fD is a Darcy friction +factor that is a function of Reynolds number NRe, with the roughness ratio ε/D as a parameter. +
+ +
+
+
Figure: Darcy friction factor as a function of the Reynolds number.
+ ++The turbulent region (NRe > 2.3×10³) is a flow regime where the velocity across the pipe has +a stochastic nature, and where the velocity v is relatively uniform across the pipe when we average the +velocity over some short period of time. The laminar region (NRe < 2.1×10³) is a flow regime +with a regular velocity v which varies as a parabola with the radius of the pipe, with zero velocity at the +pipe wall and maximal velocity at the centre of the pipe. +
+ ++Darcy friction factor varies with the roughness of the pipe surface, specified by roughness height ε. For +laminar flow in a cylindrical pipe (NRe < 2.1×10³), the Darcy friction factor fD +can be found using the following expression: +
++$$ f_\\mathrm{D} = \\frac{64}{N_\\mathrm{Re}} $$ +
++Here, the Reynolds number is found as follows: \\(N_\\mathrm{Re}=\\frac{\\rho|v|D}{\\mu}\\), where μ is the fluid viscosity. +
+ ++For turbulent flow (NRe > 2.3×10³), the Darcy friction factor is expressed as: +
++$$ f_\\mathrm{D} = \\frac{1}{\\left(2\\log_{10}\\left(\\frac{\\epsilon}{3.7D} + \\frac{5.74}{N_\\mathrm{Re}^{0.9}}\\right)\\right)^2} $$ +
+ ++In order to define the Darcy friction factor in a region between laminar and turbulent flow regimes, a +cubic polynomial interpolation is used between the laminar value at NRe=2100 and the turbulent +value at NRe=2300, with matching slopes at both endpoints to achieve global differentiability. +
+ +
+Based on the presented equations for calculation of the friction force in the waterway, two functions are
+encoded in class DarcyFriction:
+
fDarcy — calculates the Darcy friction factor. This function has the following inputs:
+Reynolds number NRe, pipe diameter D, and pipe roughness height ε. Returns the friction factor
+fD.Friction — calculates the actual friction force based on the response from the
+fDarcy function. This function has the following inputs: linear velocity v, pipe length and
+diameter L and D, liquid density and viscosity ρ and μ, and pipe roughness height ε. Returns friction force
+Ff.
+An example of a Modelica code for defining the Friction function:
+
+function Friction \"Friction force with Darcy friction factor\" + import Modelica.Constants.pi; + input Modelica.SIunits.Velocity v \"Flow velocity\"; + input Modelica.SIunits.Diameter D \"Pipe diameter\"; + input Modelica.SIunits.Length L \"Pipe length\"; + input Modelica.SIunits.Density rho \"Density\"; + input Modelica.SIunits.DynamicViscosity mu \"Dynamic viscosity of water\"; + input Modelica.SIunits.Height eps \"Pipe roughness height\"; + output Modelica.SIunits.Force F_f \"Friction force\"; +protected + Modelica.SIunits.ReynoldsNumber N_Re \"Reynolds number\"; + Real f \"friction factor\"; +algorithm + N_Re:= rho * abs(v) * D / mu; + f := fDarcy(N_Re, D, eps); + F_f := 0.5 * pi * f * rho * L * v * abs(v) * D / 4; +end Friction; +")); end DarcyFriction; diff --git a/OpenHPL/Functions/Fitting/FittingPhi.mo b/OpenHPL/Functions/Fitting/FittingPhi.mo index b75276e0..02402778 100644 --- a/OpenHPL/Functions/Fitting/FittingPhi.mo +++ b/OpenHPL/Functions/Fitting/FittingPhi.mo @@ -48,7 +48,7 @@ algorithm D_o, L); end if; - annotation ( + annotation (preferredView="info", Documentation(info = "
Define dimension factor φ for different types of fittings.
")); diff --git a/OpenHPL/Functions/Fitting/FittingVariants/Rounded.mo b/OpenHPL/Functions/Fitting/FittingVariants/Rounded.mo index 99baeac3..c7411e57 100644 --- a/OpenHPL/Functions/Fitting/FittingVariants/Rounded.mo +++ b/OpenHPL/Functions/Fitting/FittingVariants/Rounded.mo @@ -13,7 +13,7 @@ algorithm /* Rounded Expansion (same as Square Expansion) */ phi := Square(N_Re,p_eps,D_i,D_o); end if; - annotation (Documentation(info=" + annotation (preferredView="info", Documentation(info="Define dimension factor φ for Rounded Reduction. Rounded Expansion is the same as Squared Expansion.
Squared Reduction: |
diff --git a/OpenHPL/Functions/Fitting/FittingVariants/Tapered.mo b/OpenHPL/Functions/Fitting/FittingVariants/Tapered.mo
index 2f2dc02e..0f8d3d0b 100644
--- a/OpenHPL/Functions/Fitting/FittingVariants/Tapered.mo
+++ b/OpenHPL/Functions/Fitting/FittingVariants/Tapered.mo
@@ -36,7 +36,7 @@ algorithm
phi := Square(N_Re,p_eps,D_i,D_o);
end if;
end if;
- annotation (
+ annotation (preferredView="info",
Documentation(info="