29 #ifndef BoundaryFunctionsH
30 #define BoundaryFunctionsH
50 stFESub(
const double iAA,
const double iAd,
const double ig,
const double iK,
const double iCC) :
51 AA(iAA), Ad(iAd), Gam(ig), K(iK), BC(iCC) {
58 double operator()(
const double A2) {
59 double xx = A2 * invAdAA;
60 double yy = pow(xx, 2 / Ga3);
61 yy =
pow2(K) * yy - 1.;
64 U2_2 = AdAA * sqrt((
pow2(xx) - 1.) / (yy * Ga3));
77 stFESup(
const double ig,
const double iK) :
81 Ga9 = (Gam - 1) / Ga8;
85 double operator()(
const double M) {
87 double yy = Ga8 * pow(K * M, Ga9);
88 double xx = Ga3 *
pow2(M);
106 stFSSub(
const double iAA,
const double iAd,
const double ig,
const double iK,
const double iCC,
const double iAc) :
107 AA(iAA), Ad(iAd), Ac(iAc), Gam(ig), K(iK), BC(iCC) {
109 invAdAA = 1 / (AA * Ad);
114 double operator()(
const double A2) {
116 if(pow2Ac >
pow2(A2)) {
117 U2 = sqrt((pow2Ac -
pow2(A2)) / Ga3);
119 double A1 = Ac * (BC + Ga3 * U2) * invAdAA;
120 double U1 = K * U2 *
pow2(A1 / A2);
121 return pow2(A1) + Ga3 *
pow2(U1) - pow2Ac;
138 stFSSup(
const double iAA,
const double iFcc,
const double ig,
const double iK,
const double iCC,
const double iAc) :
139 AA(iAA), Ac(iAc), Gam(ig), K(iK), BC(iCC), Fcc(iFcc) {
147 double operator()(
const double U2) {
149 double dif =
pow2(Ac) - Ga3 *
pow2(U2);
156 double A2_2 = sqrt(U2 * pow((BC + Ga3 * U2) * invAA, Gam / Ga3) * invFcc);
177 stRecover(
const double iAA,
const double iAd,
const double ig,
const double iCR,
const double iCC) :
178 AA(iAA), Ad(iAd), Gam(ig), CR(iCR), BC(iCC) {
180 invAdAA = 1 / (Ad * AA);
184 double operator()(
double A2) {
185 U2 = (BC - A2) / Ga3;
186 double Ga3U22 = Ga3 *
pow2(U2);
187 double Adep1 = sqrt(
pow2(A2) + Ga3U22);
188 A1 = A2 * Adep1 * invAdAA;
189 UThroat = U2 * CR *
pow2(A1 / A2);
190 double Adep2 = sqrt(
pow2(A1) + Ga3 *
pow2(UThroat));
191 A2_2 = sqrt(
pow2(Adep2) - Ga3U22);
228 stExpansion(
const double iCCS,
const double iCCE,
const double R_E,
const double R_A,
const double iGam) :
229 CCS(iCCS), CCE(iCCE), rel_entropia(R_E), rel_area(R_A), Gam(iGam) {
235 double operator()(
const double A1) {
237 U1 = (CCS - A1) / Ga3;
238 double xx1 =
pow2(A1) + Ga3 *
pow2(U1);
239 double xx2 = rel_area *
pow2(A1) + Gam *
pow2(U1);
240 if(fabs(U1 - A1) < 1e-15) {
241 if(Ga2 / (Gam + rel_area) >= 1) {
242 U2 = (xx2 / (Ga2 * U1));
244 U2 = (xx2 / (Ga2 * U1)) * (1 - sqrt(1 -
pow2(Ga2 / (Gam + rel_area))));
247 if((xx1 * 2. * Ga2 * (
pow2(U1) /
pow2(xx2))) >= 1) {
250 U2 = xx2 * (1 - sqrt(1. - (xx1 * 2. * Ga2 * (
pow2(U1) /
pow2(xx2))))) / Ga2;
259 A2 = sqrt(xx1 - Ga3 *
pow2(U2));
260 xx3 = CCE + Ga3 * U2;
261 double xx = pow(rel_entropia * xx3, Gam);
262 double xtx = (
pow2(A2) + Gam *
pow2(U2)) / (((Gam *
pow2(U1)) / rel_area +
pow2(A1)) *
pow2(A2));
264 double A1p = xx * xtx;
294 stContraction(
const double iCCS,
const double iCCE,
const double R_E,
const double R_A,
const double iGam) :
295 CCS(iCCS), CCE(iCCE), rel_entropia(R_E), rel_area(R_A), Gam(iGam) {
300 double operator()(
const double U1) {
303 double xx1 = rel_entropia * 4. * CCE / Ga1;
304 double xx2 =
pow2(rel_entropia) / Ga3 + 1.;
305 double xx3 =
pow2(CCE) / Ga3 - Ga3 *
pow2(U1) -
pow2(A1);
307 double b =
pow2(xx1) - xx2 * 4. * xx3;
309 A2 = xx1 / (2 * xx2);
311 A2 = (xx1 + sqrt(b)) / (2. * xx2);
318 U2 = (A2 * rel_entropia - CCE) / Ga3;
319 double U1p = A1 / A2;
320 double xx = 1. / Ga3;
321 U1p = rel_area * U2 / pow(U1p, xx);
356 stPerdPresAd(
const double iCC1,
const double iCC2,
const double iFK,
const double iGam,
const double iFRE) :
357 CC1(iCC1), CC2(iCC2), FK(iFK), Gam(iGam), FRE(iFRE) {
359 Ga5 = (Gam - 1) / 2 / Gam;
362 double operator()(
const double A1) {
364 U1 = (CC1 - A1) / Ga3;
366 double b1 = -FK *
pow2(U1) /
pow2(A1) + 1.;
367 double b =
pow2(A1) * b1;
369 double u2u1 = (sqrt(
pow2(b) - 4. * a * c) - b) / (2. * a);
378 xx3 = CC2 + Ga3 * U2;
379 A1p = xx3 * FRE * pow(b1, -Ga5);
404 stPerdPresAdL(
const double iCC1,
const double iCC2,
const double iFK,
const double iGam,
const double iFRE,
405 const double iARef) :
406 CC1(iCC1), CC2(iCC2), FK(iFK), Gam(iGam), FRE(iFRE), ARef(iARef) {
408 Ga5 = (Gam - 1) / 2 / Gam;
411 double operator()(
const double A1) {
413 U1 = (CC1 - A1) / Ga3;
414 double b1 = -FK * fabs(U1) /
pow2(A1) + 1.;
418 double b = b1 *
pow2(A1);
420 double u2u1 = QuadraticEqP(a, b, c);
429 xx3 = CC2 + U2 * Ga3;
431 A1p = xx3 * FRE / pow(b1, Ga5);
452 double Gasto_calculado;
470 stComprVol(
const double iAA,
const double iCC,
const double iGam,
const double iA,
const double iGS,
const double iF,
471 const double iPRef,
const double iARef) :
472 AA(iAA), BC(iCC), Gam(iGam), A(iA), Gasto_calculado(iGS), F(iF), PRef(iPRef), ARef(iARef) {
474 Ga4 = 2 * Gam / (Gam - 1);
479 double operator()(
const double Vel) {
480 double entropia = A * AA / (BC + Ga3 * Vel);
481 U = Gasto_calculado * pow(entropia, Ga4) / (CD * Gam * F * __units::BarToPa(PRef) / ARef * pow(A, 1 / Ga3));
499 const double iGam,
const double R,
const double Sec) :
500 Mass(iMass), Lambda(iCC), Gam(iGam), Section(Sec), Eff(iEff) {
501 A01 = sqrt(Gam * R * T01);
502 AA1 = A01 * pow(P01, (1 - Gam) / 2 / Gam);
506 double operator()(
const double Match) {
508 double A2 = Lambda / (1 - Ga3 * Match);
509 double k =
pow2(A2) * (1 + Ga3 * Match * Match) / (A01 * A01) - 1;
510 double AA2 = AA1 * sqrt(1 + k) / (k * Eff + 1);
511 double Match2 = Mass / (Gam * pow(A2 / pow(AA2, Gam), 1 / Ga3) * A2 * Section);
512 return Match - Match2;
529 stCharOrigin(
const double iW00,
const double iW10,
const double iW20,
const double iW01,
const double iW11,
530 const double iW21,
const double iG0,
const double iG1,
const double idtdx,
int isigno) :
531 W00(iW00), W10(iW10), W20(iW20), W01(iW01), W11(iW11), W21(iW21), G0(iG0), G1(iG1), dtdx(idtdx), signo(isigno) {
534 double operator()(
const double x) {
535 double W0p = Interpola(W00, W01, 1., x);
536 double W1p = Interpola(W10, W11, 1., x);
537 double W2p = Interpola(W20, W21, 1., x);
538 double Gp = Interpola(G0, G1, 1., x);
539 double Up = W1p / W0p;
540 double Ap = sqrt(Gp * (Gp - 1) * (W2p / W0p -
pow2(Up) / 2.));
541 return x + signo * (Up - signo * Ap) * dtdx;
553 stPathOrigin(
const double iW00,
const double iW10,
const double iW01,
const double iW11,
const double idtdx,
555 W00(iW00), W10(iW10), W01(iW01), W11(iW11), dtdx(idtdx), signo(isigno) {
558 double operator()(
const double x) {
559 double W0p = Interpola(W00, W01, 1., x);
560 double W1p = Interpola(W10, W11, 1., x);
561 double Up = W1p / W0p;
562 return x + signo * Up * dtdx;