OpenWAM
NewCompSolver.h
1 /*--------------------------------------------------------------------------------*\
2 ==========================|
3  \\ /\ /\ // O pen | OpenWAM: The Open Source 1D Gas-Dynamic Code
4  \\ | X | // W ave |
5  \\ \/_\/ // A ction | CMT-Motores Termicos / Universidad Politecnica Valencia
6  \\/ \// M odel |
7  ----------------------------------------------------------------------------------
8  License
9 
10  This file is part of OpenWAM.
11 
12  OpenWAM is free software: you can redistribute it and/or modify
13  it under the terms of the GNU General Public License as published by
14  the Free Software Foundation, either version 3 of the License, or
15  (at your option) any later version.
16 
17  OpenWAM is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  GNU General Public License for more details.
21 
22  You should have received a copy of the GNU General Public License
23  along with OpenWAM. If not, see <http://www.gnu.org/licenses/>.
24 
25 
26  \*--------------------------------------------------------------------------------*/
27 
28 //---------------------------------------------------------------------------
29 #ifndef NewCompSolverH
30 #define NewCompSolverH
31 //---------------------------------------------------------------------------
32 
33 //#include "roots.h"
34 #include "Constantes.h"
35 
37 
38  double CC2;
39  double Gam;
40  double AA1ini;
41  double AA2ini;
42  double RelComp;
43  double A1ini;
44  double Sec1;
45  double Sec2;
46  double V1ini;
47 
48  double RelSec;
49  double V2;
50  double Gaa;
51  double Gae;
52  double n;
53  double k;
54  double AA1new;
55  double AA2new;
56  double CP;
57 
58  stCompSolverA2(double c2, double g, double RC, double cp0, double aa1, double aa2, double a1, double s1, double s2,
59  double v1) :
60  Gam(g), CC2(c2), CP(cp0), AA1ini(aa1), AA2ini(aa2), RelComp(RC), A1ini(a1), Sec1(s1), Sec2(s2), V1ini(v1) {
61  Gaa = 2 / (Gam - 1);
62  Gae = (Gam - 1) / Gam;
63  RelSec = s1 / s2;
64  }
65  double operator()(const double A2) {
66  V2 = (CC2 * __cons::ARef - A2) * (-Gaa);
67  if(V2 < 0) {
68  AA1new = AA2ini * pow(RelComp, Gae / 2.) / sqrt(1 + CP);
69  AA2new = AA2ini;
70  } else {
71  AA2new = sqrt(1 + CP) * AA1ini / pow(RelComp, Gae / 2.);
72  AA1new = AA1ini;
73  }
74  return pow(A1ini / pow(AA1new, Gam), Gaa) * V1ini * RelSec - pow(A2 / pow(AA2new, Gam), Gaa) * V2;
75  }
76 };
77 
79 
80  double CC1;
81  double CC2;
82  double AA1ini;
83  double AA2ini;
84  double Sec1;
85  double Sec2;
86  double RelComp;
87  double Rend;
88  double Gam;
89  double Dt;
90  double CPAnt;
91  double LCar;
92 
93  double RelSec;
94  double V1;
95  double V2;
96  double A2;
97  double CP;
98  double n;
99  double AA1new;
100  double AA2new;
101  double Gaa;
102  double Gae;
103  double CPfin;
104  double A2Max;
105  double A2Min;
106  double A2MaxLim;
107  double A2MinLim;
108 
109  stCompSolverA1(double c1, double c2, double s1, double s2, double g, double aa1, double aa2, double RC, double Rd,
110  double Lc, double cp0, double inct) :
111  CC1(c1), CC2(c2), AA1ini(aa1), AA2ini(aa2), Sec1(s1), Sec2(s2), RelComp(RC), Gam(g), Rend(Rd), Dt(inct), CPAnt(cp0),
112  LCar(Lc) {
113 
114  Gaa = 2 / (Gam - 1);
115  Gae = (Gam - 1) / Gam;
116  CP = (pow(RC, Gae) - 1) / Rd;
117  RelSec = s1 / s2;
118  }
119 
120  double operator()(const double A1) {
121  V1 = (CC1 * __cons::ARef - A1) * Gaa;
122  if(V1 < 0) {
123  CPfin = CPAnt + (1 / LCar * 10) * (-CPAnt) * 0.5 * fabs(V1) * Dt;
124  A2Min = CC2 * __cons::ARef + V1 * RelSec / Gaa;
125  A2Max = CC2 * __cons::ARef;
126  A2MinLim = CC2 * __cons::ARef * 2 / (Gam + 1);
127  A2MaxLim = CC2 * __cons::ARef;
128  } else if(V1 > 0) {
129  CPfin = CPAnt + (1 / LCar * 10) * (CP - CPAnt) * 0.5 * fabs(V1) * Dt;
130  A2Min = CC2 * __cons::ARef;
131  A2Max = CC2 * __cons::ARef + V1 * RelSec / Gaa;
132  A2MinLim = CC2 * __cons::ARef;
133  A2MaxLim = CC2 * __cons::ARef * 2 / (3 - Gam);
134  } else {
135  A2 = CC2 * __cons::ARef;
136  V2 = 0.;
137  return (A2 * A2) - (A1 * A1) * (1 + CPfin);
138  }
139  stCompSolverA2 FunA2(CC2, Gam, RelComp, CPAnt, AA1ini, AA2ini, A1, Sec1, Sec2, V1);
140 
141  bool Verdad = zbrac2(FunA2, A2Min, A2Max, A2MinLim, A2MaxLim);
142  if(Verdad) {
143  A2 = FindRoot(FunA2, A2Min, A2Max);
144  } else {
145  if(fabs(FunA2(A2Min)) < fabs(FunA2(A2Max))) {
146  double Error = FunA2(A2Min);
147  A2 = A2Min;
148  } else {
149  double Error = FunA2(A2Max);
150  A2 = A2Max;
151  }
152  }
153  V2 = FunA2.V2;
154 
155  AA1new = FunA2.AA1new;
156  AA2new = FunA2.AA2new;
157 
158  return (A2 * A2 + V2 * V2 / Gaa) - (A1 * A1 + V1 * V1 / Gaa) * (1 + CPfin);
159  }
160 };
161 #endif
Constantes.h
FindRoot
double FindRoot(T &func, const double x1, const double x2)
Finds the root of a function.
Definition: Math_wam.h:459
stCompSolverA1
Definition: NewCompSolver.h:78
stCompSolverA2
Definition: NewCompSolver.h:36