Documentație Tehnică

001// pdim_modern.cpp
002// Modern C++17 console rewrite of the Turbo Pascal "diafragme" program.
003// Focus: sizing and performance of primary differential-pressure flow elements (orifice plates, nozzles, Venturi).
004// The code uses ISO 5167-style formulae (as reflected in the Pascal source) and SI units.
005//
006// NOTE:
007// - This is a faithful *engineering* rewrite of the core calculations. The original Pascal program
008//   had a large DOS GUI with graphics/printer and many fluid presets via a custom "udim" unit.
009//   Here we provide a clean, portable console tool with explicit inputs.
010// - If you share the original udim unit (fluid properties), we can add a library of fluids.
011//
012// Build:  g++ -std=gnu++17 -O2 -o pdim_modern pdim_modern.cpp
013//
014// DISCLAIMER: Verify against your standards (ISO 5167 parts) and your original program outputs.
015// Small differences may exist due to rounding, fluid property source, or tap corrections.
016
017#define _USE_MATH_DEFINES
018#include <cmath>
019#include <cstdint>
020#include <iomanip>
021#include <iostream>
022#include <limits>
023#include <optional>
024#include <stdexcept>
025#include <string>
026#include <tuple>
027#include <algorithm>
028
029// ---------- Types & helpers ----------
030
031enum class PrimaryElement {
032    Orifice,            // sharp-edged orifice plate
033    ISA1932Nozzle,      // ISA 1932 nozzle
034    LongRadiusNozzle,   // long-radius nozzle
035    VenturiNozzle,
036    VenturiTube
037};
038
039enum class Tapping {
040    Corner,     // corner taps
041    Flange,     // flange taps
042    DD2         // D and D/2 taps
043};
044
045struct Inputs {
046    // Pipe & orifice
047    double D_mm;        // pipe internal diameter [mm]
048    double beta;        // d/D (ignored in "size-by-flow" mode)
049    // Process
050    double rho;         // upstream density [kg/m^3]
051    double mu;          // dynamic viscosity [Pa·s]
052    double kappa;       // cp/cv (for gases). Use 1.0 for liquids.
053    double p1_bar;      // upstream absolute pressure [bar]
054    double dp_bar;      // differential pressure across primary element [bar]
055    double T_K;         // temperature [K] (optional usage)
056    // Target / known flow (volumetric @ upstream conditions)
057    double Q_m3s;       // volumetric flow rate [m^3/s]
058    // Selections
059    PrimaryElement elem;
060    Tapping taps;
061    bool isGas;         // true for compressible expansibility; false for liquids
062};
063
064static inline double sqr ( double x ) {
065    return x * x;
066}
067static inline double pow4 ( double x ) {
068    x = x * x; return x * x;
069}
070
071// ISO 5167 alpha = C / sqrt(1 - beta^4). Some Pascal code used alfa = C * E where E = 1/sqrt(1 - beta^4).
072struct Coeff {
073    double C;       // discharge coefficient
074    double eps;     // expansibility (compressibility) factor (1 for liquids)
075    double alpha;   // C / sqrt(1 - beta^4)
076};
077
078// ---------- Correlations from Pascal code (seen in PDIM.PAS) ----------
079// These were visible in the source you provided (e.g. 0.5959 + 0.0312*b^2.1 - 0.184*b^8 + 0.0029*b^2.5*(1e6/Re)^0.75)
080// and different c2/tapping corrections.
081Coeff orifice_coeff ( double beta, double ReD, double D_mm, Tapping taps, bool isGas, double dp_bar, double kappa, double p1_bar ) {
082    const double b = beta;
083    double C = 0.5959 + 0.0312 * std::pow ( b, 2.1 ) - 0.184 * std::pow ( b, 8.0 )
084        + 0.0029 * std::pow ( b, 2.5 ) * std::pow ( 1e6 / ReD, 0.75 );
085
086    // Tap-specific correction terms observed in PDIM.PAS
087    // (These approximate flange, corner, D-D/2 effects.)
088    double c2 = 0.0;
089    switch (taps) {
090    case Tapping::Corner:
091        c2 = 0.039 * ( pow4 ( b ) / ( 1.0 - pow4 ( b ) ) ) - 0.0337 * std::pow ( b, 3.0 ) * 25.4 / ( D_mm * 1000.0 );
092        break;
093    case Tapping::Flange:
094        c2 = 0.0; // often absorbed in the base correlation; keep zero unless we later refine
095        break;
096    case Tapping::DD2:
097        c2 = 0.09 * ( 25.4 / ( D_mm * 1000.0 ) ) * ( pow4 ( b ) / ( 1.0 - pow4 ( b ) ) ) - 0.0337 * std::pow ( b, 3.0 ) * 25.4 / ( D_mm * 1000.0 );
098        break;
099    }
100    C += c2;
101
102    const double E = 1.0 / std::sqrt ( 1.0 - pow4 ( b ) );
103
104    // Expansibility factor ε for gases. In PDIM, several forms appear; here we use a common ISO-like form.
105    double eps = 1.0;
106    if (isGas) {
107        // Use ISO-like expansibility; p1_bar, dp_bar in bar
108        const double beta4 = pow4 ( b );
109        const double term = 0.351 + 0.256 * beta4 + 0.93 * beta4 * beta4; // β^8 coefficient
110        eps = 1.0 - term * ( dp_bar / ( kappa * p1_bar ) );
111        if (eps < 0.0) eps = 0.0;
112        if (eps > 1.0) eps = 1.0;
113    }
114    return { C, eps, C * E };
115}
116
117// ISA 1932 nozzle (from PDIM patterns): C ≈ 0.99 - 0.2262 β^4.1 - (...) (1e6/Re)^1.15
118Coeff isa1932_coeff ( double beta, double ReD, bool isGas, double dp_bar, double kappa, double p1_bar ) {
119    const double b = beta;
120    double C = 0.99 - 0.2262 * std::pow ( b, 4.1 )
121        - ( 0.00175 * sqr ( b ) - 0.0033 * std::pow ( b, 4.15 ) ) * std::pow ( 1e6 / ReD, 1.15 );
122    const double E = 1.0 / std::sqrt ( 1.0 - pow4 ( b ) );
123    double eps = 1.0;
124    if (isGas) {
125        // PDIM had variants like (0.484+1.54*β^4)*dp/(k*p1). Use a robust ISO-like form for nozzles as well.
126        const double term = 0.351 + 0.256 * pow4 ( b ) + 0.93 * pow4 ( b ) * pow4 ( b );
127        eps = 1.0 - term * ( dp_bar / ( kappa * p1_bar ) );
128        eps = std::clamp ( eps, 0.0, 1.0 );
129    }
130    return { C, eps, C * E };
131}
132
133// Long-radius nozzle (simplified; many implementations reuse ISA 1932 with different constants).
134// Here we reuse ISA 1932 correlation as a placeholder; adjust if you have the exact PDIM constants.
135Coeff long_radius_nozzle_coeff ( double beta, double ReD, bool isGas, double dp_bar, double kappa, double p1_bar ) {
136    return isa1932_coeff ( beta, ReD, isGas, dp_bar, kappa, p1_bar );
137}
138
139// Venturi nozzle (placeholder: often C ~ 0.985 to 0.995 with weak Re dependence)
140Coeff venturi_nozzle_coeff ( double beta, double ReD, bool isGas, double dp_bar, double kappa, double p1_bar ) {
141    const double b = beta;
142    double C = 0.985 - 0.01 * std::pow ( 1e6 / ReD, 0.05 );
143    const double E = 1.0 / std::sqrt ( 1.0 - pow4 ( b ) );
144    double eps = 1.0;
145    if (isGas) {
146        const double term = 0.351 + 0.256 * pow4 ( b ) + 0.93 * pow4 ( b ) * pow4 ( b );
147        eps = 1.0 - term * ( dp_bar / ( kappa * p1_bar ) );
148        eps = std::clamp ( eps, 0.0, 1.0 );
149    }
150    return { C, eps, C * E };
151}
152
153// Venturi tube (classical: C ~ 0.99 ...)
154Coeff venturi_tube_coeff ( double beta, double ReD, bool isGas, double dp_bar, double kappa, double p1_bar ) {
155    const double b = beta;
156    double C = 0.99 - 0.005 * std::pow ( 1e6 / ReD, 0.05 );
157    const double E = 1.0 / std::sqrt ( 1.0 - pow4 ( b ) );
158    double eps = 1.0;
159    if (isGas) {
160        const double term = 0.351 + 0.256 * pow4 ( b ) + 0.93 * pow4 ( b ) * pow4 ( b );
161        eps = 1.0 - term * ( dp_bar / ( kappa * p1_bar ) );
162        eps = std::clamp ( eps, 0.0, 1.0 );
163    }
164    return { C, eps, C * E };
165}
166
167Coeff coeffs ( const Inputs &in, double beta, double ReD ) {
168    switch (in.elem) {
169    case PrimaryElement::Orifice:          return orifice_coeff ( beta, ReD, in.D_mm, in.taps, in.isGas, in.dp_bar, in.kappa, in.p1_bar );
170    case PrimaryElement::ISA1932Nozzle:    return isa1932_coeff ( beta, ReD, in.isGas, in.dp_bar, in.kappa, in.p1_bar );
171    case PrimaryElement::LongRadiusNozzle: return long_radius_nozzle_coeff ( beta, ReD, in.isGas, in.dp_bar, in.kappa, in.p1_bar );
172    case PrimaryElement::VenturiNozzle:    return venturi_nozzle_coeff ( beta, ReD, in.isGas, in.dp_bar, in.kappa, in.p1_bar );
173    case PrimaryElement::VenturiTube:      return venturi_tube_coeff ( beta, ReD, in.isGas, in.dp_bar, in.kappa, in.p1_bar );
174    }
175    throw std::runtime_error ( "Unknown primary element" );
176}
177
178// ---------- Flow relations ----------
179
180struct Result {
181    double beta;
182    double d_mm;
183    double ReD;
184    Coeff cf;
185    double Q_m3s;     // volumetric flow [m^3/s]
186    double m_kg_s;    // mass flow [kg/s]
187    double dp_bar;    // differential pressure [bar]
188};
189
190// Compute flow for a given beta and dp
191Result compute_flow ( const Inputs &in, double beta, double dp_bar_override = -1.0 ) {
192    const double D = in.D_mm / 1000.0;        // [m]
193    const double d = beta * D;                // [m]
194    const double A = M_PI * d * d / 4.0;      // [m^2]
195
196    const double dp_bar = ( dp_bar_override > 0.0 ? dp_bar_override : in.dp_bar );
197    const double dp_Pa = dp_bar * 1e5;
198
199    // We will iterate because Re depends on Q (unknown). Use a fixed-point on velocity.
200    double ReD = 1e6; // initial guess for Re used in correlations
201    Coeff cf = coeffs ( in, beta, ReD );
202
203    // Fixed-point iteration on Q
204    double Q = 0.0;
205    for (int it = 0; it < 20; ++it) {
206        cf = coeffs ( in, beta, ReD );
207        const double alpha_eps = cf.alpha * ( in.isGas ? cf.eps : 1.0 );
208        // Q = (alpha*beta^2*eps) * A_pipe * sqrt(2*dp/rho) ?
209        // From PDIM snippet: deb = (alpha * beta^2 * eps) * (pi*d2^2/4) * sqrt(2*rho*dp)
210        // Interpreting d2 as pipe diameter; here we use orifice throat area A and upstream rho.
211        const double rec = cf.alpha * beta * beta * ( in.isGas ? cf.eps : 1.0 );
212        const double Q_new = rec * ( M_PI * D * D / 4.0 ) * std::sqrt ( 2.0 * dp_Pa / in.rho );
213        // Update Re using throat velocity based on Q through orifice area
214        const double v = Q_new / A;
215        ReD = in.rho * v * d / in.mu;
216        if (std::abs ( Q_new - Q ) < 1e-12) {
217            Q = Q_new; break;
218        }
219        Q = Q_new;
220    }
221
222    const double m = Q * in.rho;
223    return { beta, d * 1000.0, ReD, cf, Q, m, dp_bar };
224}
225
226// Solve for beta (size) given target Q and dp using bisection
227Result size_by_flow ( const Inputs &in, double Q_target, double bmin = 0.20, double bmax = 0.75 ) {
228    auto f = [&] ( double b ) {
229        return compute_flow ( in, b ).Q_m3s - Q_target;
230        };
231    double lo = bmin, hi = bmax;
232    double flo = f ( lo ), fhi = f ( hi );
233    if (flo * fhi > 0.0) {
234        // Try to expand the bracket a bit
235        lo = 0.10; hi = 0.80;
236        flo = f ( lo ); fhi = f ( hi );
237        if (flo * fhi > 0.0) {
238            throw std::runtime_error ( "Unable to bracket solution for beta. Check inputs." );
239        }
240    }
241    for (int it = 0; it < 80; ++it) {
242        const double mid = 0.5 * ( lo + hi );
243        const double fmid = f ( mid );
244        if (std::abs ( fmid ) < 1e-12 || ( hi - lo ) < 1e-6) {
245            return compute_flow ( in, mid );
246        }
247        if (flo * fmid <= 0.0) {
248            hi = mid; fhi = fmid;
249        }
250        else {
251            lo = mid; flo = fmid;
252        }
253    }
254    return compute_flow ( in, 0.5 * ( lo + hi ) );
255}
256
257// Pretty print
258void print_result ( const Result &r ) {
259    std::cout << std::fixed << std::setprecision ( 6 );
260    std::cout << "beta         : " << r.beta << "\n";
261    std::cout << "orifice d    : " << r.d_mm << " mm\n";
262    std::cout << "Re_D         : " << r.ReD << "\n";
263    std::cout << "C            : " << r.cf.C << "\n";
264    std::cout << "epsilon      : " << r.cf.eps << "\n";
265    std::cout << "alpha        : " << r.cf.alpha << "\n";
266    std::cout << "Q            : " << r.Q_m3s << " m^3/s\n";
267    std::cout << "m_dot        : " << r.m_kg_s << " kg/s\n";
268    std::cout << "dp           : " << r.dp_bar << " bar\n";
269}
270
271// ---------- CLI ----------
272
273int main () {
274    std::cout << "Differential-Pressure Flow Primary (modern rewrite of PDIM.PAS)\n";
275    std::cout << "Units: D in mm, p in bar, mu in Pa·s, rho in kg/m^3, Q in m^3/s, T in K.\n\n";
276
277    Inputs in{};
278    int elemChoice = 1, tapChoice = 2, mode = 1;
279    std::cout << "Primary element [1=Orifice, 2=ISA1932 nozzle, 3=Long-radius nozzle, 4=Venturi nozzle, 5=Venturi tube]: ";
280    std::cin >> elemChoice;
281    in.elem = ( elemChoice == 1 ? PrimaryElement::Orifice :
282        elemChoice == 2 ? PrimaryElement::ISA1932Nozzle :
283        elemChoice == 3 ? PrimaryElement::LongRadiusNozzle :
284        elemChoice == 4 ? PrimaryElement::VenturiNozzle : PrimaryElement::VenturiTube );
285
286    if (in.elem == PrimaryElement::Orifice) {
287        std::cout << "Tapping [1=Corner, 2=Flange, 3=D & D/2]: ";
288        std::cin >> tapChoice;
289        in.taps = ( tapChoice == 1 ? Tapping::Corner : tapChoice == 3 ? Tapping::DD2 : Tapping::Flange );
290    }
291    else {
292        in.taps = Tapping::Flange;
293    }
294
295    std::cout << "Is the fluid a gas? [1=yes, 0=no]: ";
296    int gasInt; std::cin >> gasInt; in.isGas = ( gasInt != 0 );
297
298    std::cout << "Pipe inner diameter D [mm]: ";
299    std::cin >> in.D_mm;
300
301    std::cout << "Upstream density rho [kg/m^3]: ";
302    std::cin >> in.rho;
303
304    std::cout << "Dynamic viscosity mu [Pa*s]: ";
305    std::cin >> in.mu;
306
307    if (in.isGas) {
308        std::cout << "Heat capacity ratio kappa = cp/cv [-]: ";
309        std::cin >> in.kappa;
310    }
311    else {
312        in.kappa = 1.0;
313    }
314
315    std::cout << "Upstream absolute pressure p1 [bar]: ";
316    std::cin >> in.p1_bar;
317
318    std::cout << "Differential pressure dp [bar]: ";
319    std::cin >> in.dp_bar;
320
321    std::cout << "Temperature T [K] (enter, used only for info): ";
322    std::cin >> in.T_K;
323
324    std::cout << "\nMode: 1) Size by target flow (compute beta,d)   2) Rate given beta (compute Q)\nSelect 1/2: ";
325    std::cin >> mode;
326
327    if (mode == 1) {
328        std::cout << "Target volumetric flow Q [m^3/s]: ";
329        std::cin >> in.Q_m3s;
330        // Solve for beta
331        try {
332            Result r = size_by_flow ( in, in.Q_m3s );
333            std::cout << "\n--- RESULT (size by flow) ---\n";
334            print_result ( r );
335        }
336        catch (const std::exception &e) {
337            std::cerr << "ERROR: " << e.what () << "\n";
338            return 1;
339        }
340    }
341    else {
342        std::cout << "Given beta d/D [-]: ";
343        std::cin >> in.beta;
344        Result r = compute_flow ( in, in.beta );
345        std::cout << "\n--- RESULT (rate given beta) ---\n";
346        print_result ( r );
347    }
348    return 0;
349}