390 lines
14 KiB
C#
390 lines
14 KiB
C#
///
|
|
/// Copyright (c) 2013-2018 Sensus Slovensko a.s.
|
|
///
|
|
using System;
|
|
using System.Collections.Generic;
|
|
using log4net;
|
|
|
|
namespace Config
|
|
{
|
|
public static class Formulas
|
|
{
|
|
private static readonly ILog log = LogManager.GetLogger(typeof(Formulas));
|
|
|
|
///
|
|
/// Tables to calculate specific enthalpy
|
|
///
|
|
private static readonly int[] Ii;
|
|
private static readonly int[] Ji;
|
|
private static readonly double[] ni;
|
|
|
|
private static readonly double[] Di;
|
|
|
|
/// <summary>
|
|
/// Constructor
|
|
/// </summary>
|
|
static Formulas()
|
|
{
|
|
///
|
|
/// Initialize tables to calculate specific enthalpies
|
|
///
|
|
Ii = new int[34] { 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32 };
|
|
Ji = new int[34] { -2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41 };
|
|
ni = new double[34] {
|
|
0.14632971213167, /// 1
|
|
-0.84548187169114, /// 2
|
|
-0.37563603672040E1, /// 3
|
|
0.33855169168385E1, /// 4
|
|
-0.95791963387872, /// 5
|
|
0.15772038513228, /// 6
|
|
-0.16616417199501E-1, /// 7
|
|
0.81214629983568E-3, /// 8
|
|
0.28319080123804E-3, /// 9
|
|
-0.60706301565874E-3, /// 10
|
|
-0.18990068218419E-1, /// 11
|
|
-0.32529748770505E-1, /// 12
|
|
-0.21841717175414E-1, /// 13
|
|
-0.52838357969930E-4, /// 14
|
|
-0.47184321073267E-3, /// 15
|
|
-0.30001780793026E-3, /// 16
|
|
0.47661393906987E-4, /// 17
|
|
-0.44141845330846E-5, /// 18
|
|
-0.72694996297594E-15, /// 19
|
|
-0.31679644845054E-4, /// 20
|
|
-0.28270797985312E-5, /// 21
|
|
-0.85205128120103E-9, /// 22
|
|
-0.22425281908000E-5, /// 23
|
|
-0.65171222895601E-6, /// 24
|
|
-0.14341729937924E-12, /// 25
|
|
-0.40516996860117E-6, /// 26
|
|
-0.12734301741641E-8, /// 27
|
|
-0.17424871230634E-9, /// 28
|
|
-0.68762131295531E-18, /// 29
|
|
0.14478307828521E-19, /// 30
|
|
0.26335781662795E-22, /// 31
|
|
-0.11947622640071E-22, /// 32
|
|
0.18228094581404E-23, /// 33
|
|
-0.93537087292458E-25, /// 34
|
|
};
|
|
|
|
///
|
|
/// Initialize a table to calculate temperature of a platinum thermometer from resistance
|
|
///
|
|
Di = new double[]
|
|
{
|
|
439.932854,
|
|
472.418020,
|
|
37.684494,
|
|
7.472018,
|
|
2.920828,
|
|
0.005184,
|
|
-0.963864,
|
|
-0.188732,
|
|
0.191203,
|
|
0.049025,
|
|
};
|
|
}
|
|
|
|
/// <summary>
|
|
/// Calculate density of distilled water from temperature
|
|
/// </summary>
|
|
/// <param name="t">ITS-90 temperature in [°C]</param>
|
|
/// <returns>Density in [kg/m3]</returns>
|
|
public static double DistilledWaterDensityFromTemp(double t)
|
|
{
|
|
if (t <= 40)
|
|
{
|
|
const double c0 = 999.839564;
|
|
const double c1 = 0.067998613;
|
|
const double c2 = -0.0091101468;
|
|
const double c3 = 0.00010058299;
|
|
const double c4 = -0.0000011275659;
|
|
const double c5 = 6.5985371e-09;
|
|
|
|
return ((((c5 * t + c4) * t + c3) * t + c2) * t + c1) * t + c0;
|
|
}
|
|
else
|
|
{
|
|
const double a0 = 9.9983952E2;
|
|
const double a1 = 1.6952577E1;
|
|
const double a2 = -7.9905127E-3;
|
|
const double a3 = -4.6241757E-5;
|
|
const double a4 = 1.0584601E-7;
|
|
const double a5 = -2.8103006E-10;
|
|
const double b = 1.6887236E-2;
|
|
|
|
return (((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t + a0) / (1.0 + b * t);
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Calculate density of distilled water from temperature (obsolete)
|
|
/// </summary>
|
|
/// <param name="t">IPTS-68 temperature in [°C]</param>
|
|
/// <returns>Density in [kg/m3]</returns>
|
|
public static double DistilledWaterDensityFromTempIPTS68(double t)
|
|
{
|
|
const double a0 = 999.842594;
|
|
const double a1 = 0.06793952;
|
|
const double a2 = -0.009095290;
|
|
const double a3 = 0.0001001685;
|
|
const double a4 = -0.000001120083;
|
|
const double a5 = 6.536332e-09;
|
|
|
|
return ((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t + a0;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Calculate density by comparing calculated data and data from a certificate
|
|
/// </summary>
|
|
/// <param name="realDensity">Density from a certificate in [kg/m3]</param>
|
|
/// <param name="atTemperature">Temperature from a certificate in [°C]</param>
|
|
/// <returns>Density correction in [kg/m3]</returns>
|
|
public static double DensityCorrection(double realDensity, double atTemperature)
|
|
{
|
|
/// Calculated data
|
|
double calculatedDensity = DistilledWaterDensityFromTemp(atTemperature);
|
|
|
|
return realDensity - calculatedDensity;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Calculate corrected (real) water density from temperature
|
|
/// </summary>
|
|
/// <param name="t">Temperature in [°C]</param>
|
|
/// <returns>Density in [kg/m3]</returns>
|
|
public static double WaterDensityFromTemp(double t, double realDensity, double atTemperature)
|
|
{
|
|
return DistilledWaterDensityFromTemp(t) + DensityCorrection(realDensity, atTemperature);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Calculate corrected (real) water density from temperature
|
|
/// </summary>
|
|
/// <param name="t">Temperature in [°C]</param>
|
|
/// <returns>Density in [kg/m3]</returns>
|
|
public static double WaterDensityFromTempPress(double temp, double pressure, double realDensity, double atTemperature)
|
|
{
|
|
double x0 = 5.08821E-10;
|
|
double x1 = 1.2639418;
|
|
double x2 = 0.2660269;
|
|
double x3 = 0.3734838;
|
|
double x4 = 2.0205242;
|
|
double theta = temp / 100.0;
|
|
double B = x0 * (((x3 * theta + x2) * theta + x1) * theta + 1) / (1 + x4 * theta);
|
|
|
|
return WaterDensityFromTemp(temp, realDensity, atTemperature) * (1 + B * Config.Units.ConvertTo(Config.Unit.Pa, pressure));
|
|
}
|
|
|
|
|
|
public static float AirDensityFromAmbientVales(float tempC, float pressureBar, float humiPct)
|
|
{
|
|
double pressurePa = 100000.0 * (double)pressureBar; /// [Pa]
|
|
double tempKelvin = 273.15 + (double)tempC;
|
|
double coef1 = 1.2811805 / 10000.0 * tempKelvin * tempKelvin
|
|
- 1.950987 / 100.0 * tempKelvin
|
|
+ 34.04926034
|
|
- 6.353631 * 1000.0 / tempKelvin;
|
|
double coef3 = humiPct / 100.0 * System.Math.Exp(coef1) / pressurePa;
|
|
double airDensityKgm3 = 0.00348353 * pressurePa * (1.0 - 0.378 * coef3) / tempKelvin; /// kg/m3
|
|
return (float)airDensityKgm3;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Convert 'pulses' to 'volume', prevent division by zero
|
|
/// </summary>
|
|
public static double VolumeFromPulses(int pulses, double pulsesPerLiter)
|
|
{
|
|
if (pulsesPerLiter <= double.Epsilon) return 0;
|
|
return Convert.ToDouble(pulses) / pulsesPerLiter;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Calculate the error in % from 'measured' and 'true' volume, prevent division by zero
|
|
/// </summary>
|
|
public static double ErrorFromVolumes(double measuredVolume, double trueVolume)
|
|
{
|
|
if (trueVolume <= float.Epsilon)
|
|
{
|
|
if (measuredVolume <= float.Epsilon)
|
|
{
|
|
log.WarnFormat("ErrorFromVolumes({0},{1}) returns {2}", measuredVolume, trueVolume, -100.0);
|
|
return -100.0;
|
|
}
|
|
|
|
log.WarnFormat("ErrorFromVolumes({0},{1}) returns {2}", measuredVolume, trueVolume, 99.0);
|
|
return 99.0;
|
|
}
|
|
|
|
double error = 100.0 * (measuredVolume - trueVolume) / trueVolume;
|
|
log.InfoFormat("ErrorFromVolumes({0},{1}) returns {2}", measuredVolume, trueVolume, error);
|
|
return error;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Calculates corrected value form a list of corrections by interpolation.
|
|
/// It is assumed that values in the list 'corrections' are sorted.
|
|
/// </summary>
|
|
/// <param name="rawMeasurement">Raw uncorrected value</param>
|
|
/// <param name="corrections">Sorted (value, correction) pairs</param>
|
|
/// <returns>Corrected value</returns>
|
|
public static double CorrectedValue(double rawValue, IList<Config.Entities.MeasurementCorrection> corrections)
|
|
{
|
|
if (corrections == null || corrections.Count == 0) return rawValue;
|
|
|
|
if (rawValue < corrections[0].Measurement)
|
|
{
|
|
return rawValue + corrections[0].Correction;
|
|
}
|
|
|
|
int count = corrections.Count;
|
|
for (int i = 1; i < count; i++)
|
|
{
|
|
if (rawValue < corrections[i].Measurement)
|
|
{
|
|
double d1 = rawValue - corrections[i-1].Measurement;
|
|
double d2 = corrections[i].Measurement - rawValue;
|
|
|
|
double corr;
|
|
if (d1 + d2 <= float.Epsilon)
|
|
{
|
|
corr = (corrections[i - 1].Correction + corrections[i].Correction) / 2.0;
|
|
}
|
|
else
|
|
{
|
|
corr = (corrections[i - 1].Correction * d2 + corrections[i].Correction * d1) / (d1 + d2);
|
|
}
|
|
return rawValue + corr;
|
|
}
|
|
}
|
|
|
|
return rawValue + corrections[count - 1].Correction;
|
|
}
|
|
|
|
|
|
/// <summary>
|
|
/// Converts a measurement error to a correction (used when preparing correction tables).
|
|
/// </summary>
|
|
/// <param name="measuredValue">Measured value (in arbitrary units)</param>
|
|
/// <param name="error">Measurement error in %</param>
|
|
/// <returns>Correction in the same units as the measured value</returns>
|
|
public static double CorrectionFromError(double measuredValue, double error)
|
|
{
|
|
double trueValue = measuredValue / (1 + error/100);
|
|
double correction = trueValue - measuredValue;
|
|
return correction;
|
|
}
|
|
|
|
|
|
/// <summary>
|
|
/// Calculates the heat coefficient for water
|
|
/// </summary>
|
|
/// <param name="pressure">Pressure [bar]</param>
|
|
/// <param name="T_in">Inlet temperature [°C]</param>
|
|
/// <param name="T_out">Outlet temperature [°C]</param>
|
|
/// <param name="flowMeasuredAtInlet">true = flow measured @inlet, false = flow measured @outlet</param>
|
|
/// <returns> Heat coefficient for water [J/(m3 K)]</returns>
|
|
public static double HeatCoefficientWater(double pressure, double T_in, double T_out, bool flowMeasuredAtInlet)
|
|
{
|
|
if (T_in == T_out) return 0;
|
|
|
|
const double R = 461.526; /// [J kg^-1 K^-1]
|
|
const double p_star_Pa = 16.53E6; /// [Pa] (=16.53 MPa)
|
|
const double T_star = 1386.0; /// [K]
|
|
|
|
double T_in_K = Config.Units.ConvertTo(Config.Unit.K, T_in);
|
|
double T_out_K = Config.Units.ConvertTo(Config.Unit.K, T_out);
|
|
double tau_in = T_star / T_in_K;
|
|
double tau_out = T_star / T_out_K;
|
|
double pi = Config.Units.ConvertTo(Config.Unit.Pa, pressure) / p_star_Pa;
|
|
|
|
double h_in = tau_in * GammaTau(pi, tau_in) * R * T_in_K;
|
|
double h_out = tau_out * GammaTau(pi, tau_out) * R * T_out_K;
|
|
|
|
double ni = flowMeasuredAtInlet ? GammaPi(pi, tau_in) * R * T_in_K / p_star_Pa
|
|
: GammaPi(pi, tau_out) * R * T_out_K / p_star_Pa;
|
|
|
|
return (h_in - h_out) / (ni * (T_in - T_out));
|
|
}
|
|
|
|
/// <summary>
|
|
/// gamma(pi) see also STN EN 1434-1 Annex A (A.4)
|
|
/// </summary>
|
|
/// <param name="pi">pi = p / p* where p* = 16.53 MPa</param>
|
|
/// <param name="tau">tau = T* / T where T* = 1386 K</param>
|
|
/// <returns>gamma(pi)</returns>
|
|
static double GammaPi(double pi, double tau)
|
|
{
|
|
double result = 0;
|
|
for (int i = 0; i < 34; i++)
|
|
{
|
|
result -= ni[i] * Ii[i] * Math.Pow(7.1 - pi, Ii[i] - 1) * Math.Pow(tau - 1.222, Ji[i]);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// <summary>
|
|
/// gamma(tau) see also STN EN 1434-1 Annex A (A.7)
|
|
/// </summary>
|
|
/// <param name="pi">pi = p / p* where p* = 16.53 MPa</param>
|
|
/// <param name="tau">tau = T* / T where T* = 1386 K</param>
|
|
/// <returns>gamma(tau)</returns>
|
|
static double GammaTau(double pi, double tau)
|
|
{
|
|
double result = 0;
|
|
for (int i = 0; i < 34; i++)
|
|
{
|
|
result += ni[i] * Math.Pow(7.1 - pi, Ii[i]) * Ji[i] * Math.Pow(tau - 1.222, Ji[i] - 1);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Conversion of measured resistance of a platinum thermometer to temperature according to ITS-90
|
|
/// </summary>
|
|
/// <param name="R">Measured resistance in [°C]</param>
|
|
/// <param name="R001C">Calibrated resistance in Ohm at 0.01°C</param>
|
|
/// <param name="a7">Calibrated ITS-90 coefficient a7</param>
|
|
/// <param name="b7">Calibrated ITS-90 coefficient b7</param>
|
|
/// <param name="c7">Calibrated ITS-90 coefficient c7</param>
|
|
/// <returns>Temperature in [°C]</returns>
|
|
public static double PlatinumResistanceTM_ITS90_R2T(double R, double R001C, double a7, double b7, double c7)
|
|
{
|
|
double w = R / R001C; /// ratio
|
|
double r1 = w - 1.0;
|
|
double dw = r1 * (a7 + r1 * (b7 + r1 * c7)); /// = a7*r1 + b7*r1^2 + c7*r1^3
|
|
double wr = w - dw;
|
|
double x = (wr - 2.64) / 1.64;
|
|
|
|
double sum = 0;
|
|
for (int i = Di.Length - 1; i >= 0; i--)
|
|
{
|
|
sum = sum * x + Di[i];
|
|
}
|
|
|
|
return sum;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Conversion of measured resistance of a platinum thermometer to temperature using Callendar-Van Dusen equations
|
|
/// </summary>
|
|
/// <param name="R">Measured resistance in [°C]</param>
|
|
/// <param name="R0">Calibrated resistance in Ohm at 0°C</param>
|
|
/// <param name="A">Calibration coefficient a</param>
|
|
/// <param name="B">Calibration coefficient b</param>
|
|
/// <returns>Temperature in [°C]</returns>
|
|
public static double PlatinumResistanceTM_ITS27_R2T(double R, double R0, double A, double B)
|
|
{
|
|
if (R0 * R0 * A * A - 4 * R0 * B * (R0 - R) <= 0) return 0; /// Out of range
|
|
|
|
return (-(R0 * A) + Math.Sqrt(R0 * R0 * A * A - 4 * R0 * B * (R0 - R))) / (2 * R0 * B);
|
|
}
|
|
|
|
|
|
public static double Buoyancy()
|
|
{
|
|
return 1.00103;
|
|
}
|
|
}
|
|
}
|