@@ -34,7 +34,7 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
3434 int maxiter , iter , maxiterx , iterx , CDNoz ;
3535 int interpErr = 0 ;
3636
37-
37+
3838 /* Determine Nozzle Type */
3939 /*prm->SwitchType: */
4040 /* 1: Convergent Nozzle */
@@ -54,9 +54,9 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
5454 /* Where gas constant is R = f(FAR), but NOT P & T */
5555 Rt = interp1Ac (prm -> Y_N_FARVec ,prm -> T_N_RtArray ,FARcIn ,prm -> A ,& interpErr );
5656 if (interpErr == 1 && * (prm -> IWork + Er1 )== 0 ){
57- #ifdef MATLAB_MEX_FILE
57+ #ifdef MATLAB_MEX_FILE
5858 printf ("Warning in %s, Error calculating Rt1. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
59- #endif
59+ #endif
6060 * (prm -> IWork + Er1 ) = 1 ;
6161 }
6262 Rs = Rt ;
@@ -67,19 +67,19 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
6767 if (Ptin <= PambIn ) {
6868 Ptin = PambIn + 0.1 ;
6969 if (* (prm -> IWork + Er2 )== 0 ){
70- #ifdef MATLAB_MEX_FILE
70+ #ifdef MATLAB_MEX_FILE
7171 printf ("Warning in %s, Backflow warning PtIn <= Pamb\n" , prm -> BlkNm );
72- #endif
72+ #endif
7373 * (prm -> IWork + Er2 ) = 1 ;
7474 }
7575 }
7676 /* Determine ideal velocity defined by perfect expansion to Pambient */
7777 PcalcStat (Ptin , PambIn , TtIn , htin , FARcIn , Rt , & Sin , & Ts , & hs , & rhos , & V );
7878 gammas_s = interp2Ac (prm -> Y_N_FARVec ,prm -> X_N_TtVec ,prm -> T_N_MAP_gammaArray ,FARcIn ,Ts ,prm -> A ,prm -> B ,& interpErr );
7979 if (interpErr == 1 && * (prm -> IWork + Er3 )== 0 ){
80- #ifdef MATLAB_MEX_FILE
80+ #ifdef MATLAB_MEX_FILE
8181 printf ("Warning in %s, Error calculating gammas. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
82- #endif
82+ #endif
8383 * (prm -> IWork + Er3 ) = 1 ;
8484 }
8585 MN_s = V * divby (sqrtT (gammas_s * Rs * Ts * C_GRAVITY * JOULES_CONST ));
@@ -93,9 +93,9 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
9393 MNg = 1 ;
9494 gammatg = interp2Ac (prm -> Y_N_FARVec ,prm -> X_N_TtVec ,prm -> T_N_MAP_gammaArray ,FARcIn ,TtIn ,prm -> A ,prm -> B ,& interpErr );
9595 if (interpErr == 1 && * (prm -> IWork + Er4 )== 0 ){
96- #ifdef MATLAB_MEX_FILE
96+ #ifdef MATLAB_MEX_FILE
9797 printf ("Warning in %s, Error calculating gammatg. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
98- #endif
98+ #endif
9999 * (prm -> IWork + Er4 ) = 1 ;
100100 }
101101 /* use isentropic equations for a first cut guess */
@@ -106,9 +106,9 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
106106 PcalcStat (Ptin , PsMNg , TtIn , htin , FARcIn , Rt , & Sin , & TsMNg , & hsg , & rhosg , & Vg );
107107 gammasg = interp2Ac (prm -> Y_N_FARVec ,prm -> X_N_TtVec ,prm -> T_N_MAP_gammaArray ,FARcIn ,TsMNg ,prm -> A ,prm -> B ,& interpErr );
108108 if (interpErr == 1 && * (prm -> IWork + Er4 )== 0 ){
109- #ifdef MATLAB_MEX_FILE
109+ #ifdef MATLAB_MEX_FILE
110110 printf ("Warning in %s, Error calculating gammasg. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
111- #endif
111+ #endif
112112 * (prm -> IWork + Er4 ) = 1 ;
113113 }
114114 MNg = Vg * divby (sqrtT (gammasg * Rs * TsMNg * C_GRAVITY * JOULES_CONST ));
@@ -132,9 +132,9 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
132132 PcalcStat (Ptin , PsMNg , TtIn , htin , FARcIn , Rt , & Sin , & TsMNg , & hsg , & rhosg , & Vg );
133133 gammasg = interp2Ac (prm -> Y_N_FARVec ,prm -> X_N_TtVec ,prm -> T_N_MAP_gammaArray ,FARcIn ,TsMNg ,prm -> A ,prm -> B ,& interpErr );
134134 if (interpErr == 1 && * (prm -> IWork + Er5 )== 0 ){
135- #ifdef MATLAB_MEX_FILE
135+ #ifdef MATLAB_MEX_FILE
136136 printf ("Warning in %s, Error calculating iteration gammasg. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
137- #endif
137+ #endif
138138 * (prm -> IWork + Er5 ) = 1 ;
139139 }
140140 MNg = Vg * divby (sqrtT (gammasg * Rs * TsMNg * C_GRAVITY * JOULES_CONST ));
@@ -146,9 +146,9 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
146146 iter = iter + 1 ;
147147 }
148148 if (iter == maxiter && * (prm -> IWork + Er6 )== 0 ){
149- #ifdef MATLAB_MEX_FILE
149+ #ifdef MATLAB_MEX_FILE
150150 printf ("Warning in %s, Error calculating Ps at MN = 1.\n" , prm -> BlkNm );
151- #endif
151+ #endif
152152 * (prm -> IWork + Er6 ) = 1 ;
153153 }
154154 /* MN = 1 parameters */
@@ -180,9 +180,9 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
180180 MNth = 1 ;
181181 gammasth = interp2Ac (prm -> Y_N_FARVec ,prm -> X_N_TtVec ,prm -> T_N_MAP_gammaArray ,FARcIn ,Tsth ,prm -> A ,prm -> B ,& interpErr );
182182 if (interpErr == 1 && * (prm -> IWork + Er7 )== 0 ){
183- #ifdef MATLAB_MEX_FILE
183+ #ifdef MATLAB_MEX_FILE
184184 printf ("Warning in %s, Error calculating iteration gammasg. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
185- #endif
185+ #endif
186186 * (prm -> IWork + Er7 ) = 1 ;
187187 }
188188 Vth = MNth * sqrtT (gammasth * Rs * Tsth * C_GRAVITY * JOULES_CONST );
@@ -191,9 +191,9 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
191191
192192 /* error('Nozzle Error: Negative Mach number!!') */
193193 if (MNth < 0 && * (prm -> IWork + Er8 )== 0 ){
194- #ifdef MATLAB_MEX_FILE
194+ #ifdef MATLAB_MEX_FILE
195195 printf ("Error in %s: negative throat mach number, MN = %f.\n" , prm -> BlkNm , MNth );
196- #endif
196+ #endif
197197 * (prm -> IWork + Er8 ) = 1 ;
198198 }
199199
@@ -207,16 +207,16 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
207207 /* look up Flow Coefficient */
208208 CdTh = interp1Ac (prm -> X_N_PEQPaVec ,prm -> T_N_CdThArray ,PQPaMap ,prm -> B1 ,& interpErr );
209209 if (interpErr == 1 && * (prm -> IWork + Er9 )== 0 ){
210- #ifdef MATLAB_MEX_FILE
210+ #ifdef MATLAB_MEX_FILE
211211 printf ("Warning in %s, Error calculating CdTh. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
212- #endif
212+ #endif
213213 * (prm -> IWork + Er9 ) = 1 ;
214214 }
215215 Therm_growth = interp1Ac (prm -> X_N_TtVecTG ,prm -> T_N_TGArray ,TtIn ,prm -> C ,& interpErr );
216216 if (interpErr == 1 && * (prm -> IWork + Er10 )== 0 ){
217- #ifdef MATLAB_MEX_FILE
217+ #ifdef MATLAB_MEX_FILE
218218 printf ("Warning in %s, Error calculating Therm_growth. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
219- #endif
219+ #endif
220220 * (prm -> IWork + Er10 ) = 1 ;
221221 }
222222
@@ -225,9 +225,9 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
225225 if (prm -> IDes < 0.5 ) {
226226 Ath = WIn * C_PSItoPSF * divby (Therm_growth * (1 - prm -> flowLoss /100 )* CdTh * rhosth * Vth );
227227 if (choked == 0 && * (prm -> IWork + Er11 )== 0 ){
228- #ifdef MATLAB_MEX_FILE
228+ #ifdef MATLAB_MEX_FILE
229229 printf ("Warning in %s, Calculating prm->IDes Area with un-choked nozzle.\n" , prm -> BlkNm );
230- #endif
230+ #endif
231231 * (prm -> IWork + Er11 ) = 1 ;
232232 }
233233 }
@@ -303,31 +303,43 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
303303 iterx = iterx + 1 ;
304304 }
305305 if (iterx == maxiterx && * (prm -> IWork + Er12 )== 0 ){
306- #ifdef MATLAB_MEX_FILE
306+ #ifdef MATLAB_MEX_FILE
307307 printf ("Warning in %s, Error calculating Ps at exit.\n" , prm -> BlkNm );
308- #endif
308+ #endif
309309 * (prm -> IWork + Er12 ) = 1 ;
310310 }
311311 /* Collect data from the expansion to exit area: */
312- Tsx = Ts ;
313- Vx = V ;
314- Psx = Psxg ;
315- rhosx = rhos ;
316- gammasx = interp2Ac (prm -> Y_N_FARVec ,prm -> X_N_TtVec ,prm -> T_N_MAP_gammaArray ,FARcIn ,Ts ,prm -> A ,prm -> B ,& interpErr );
317- if (interpErr == 1 && * (prm -> IWork + Er13 )== 0 ){
318- #ifdef MATLAB_MEX_FILE
319- printf ("Warning in %s, Error calculating gammas. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
320- #endif
321- * (prm -> IWork + Er13 ) = 1 ;
312+ if (Psxg > PambIn ) {
313+ /* if Psxg is greater than Pamb, use calculated exit values */
314+ Tsx = Ts ;
315+ Vx = V ;
316+ Psx = Psxg ;
317+ rhosx = rhos ;
318+ gammasx = interp2Ac (prm -> Y_N_FARVec ,prm -> X_N_TtVec ,prm -> T_N_MAP_gammaArray ,FARcIn ,Ts ,prm -> A ,prm -> B ,& interpErr );
319+ if (interpErr == 1 && * (prm -> IWork + Er13 )== 0 ){
320+ #ifdef MATLAB_MEX_FILE
321+ printf ("Warning in %s, Error calculating gammas. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
322+ #endif
323+ * (prm -> IWork + Er13 ) = 1 ;
324+ }
325+ MNx = Vx * divby (sqrtT (gammasx * Rs * Tsx * C_GRAVITY * JOULES_CONST ));
326+ }
327+ else {
328+ /* if Psxg is less than Pamb, use exit values consistant with ambient pressure and issue warning */
329+ Psx = PambIn ;
330+ Tsx = Ts_s ;
331+ Vx = V_s ;
332+ rhosx = rhos_s ;
333+ MNx = MN_s ;
334+ gammasx = gammas_s ;
322335 }
323- MNx = Vx * divby (sqrtT (gammasx * Rs * Tsx * C_GRAVITY * JOULES_CONST ));
324336 }
325337
326338 /* error('Nozzle Error: Negative Mach number!!') */
327339 if (MNx < 0 && * (prm -> IWork + Er14 )== 0 ){
328- #ifdef MATLAB_MEX_FILE
340+ #ifdef MATLAB_MEX_FILE
329341 printf ("Error in %s: negative exit mach number, MN = %f.\n" , prm -> BlkNm , MNx );
330- #endif
342+ #endif
331343 * (prm -> IWork + Er14 ) = 1 ;
332344 }
333345
@@ -338,19 +350,19 @@ void Nozzle_TMATS_body(double* y, const double* u, const NozzleStruct* prm)
338350 if (prm -> CfgEn < 0.5 ){
339351 Cv = interp1Ac (prm -> X_N_PEQPaVec ,prm -> T_N_CvArray ,PQPaMap ,prm -> B1 ,& interpErr );
340352 if (interpErr == 1 && * (prm -> IWork + Er15 )== 0 ){
341- #ifdef MATLAB_MEX_FILE
353+ #ifdef MATLAB_MEX_FILE
342354 printf ("Warning in %s, Error calculating Cv. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
343- #endif
355+ #endif
344356 * (prm -> IWork + Er15 ) = 1 ;
345357 }
346358 Cfg = 1 ;
347359 }
348360 else {
349361 Cfg = interp1Ac (prm -> X_N_PEQPaVec ,prm -> T_N_CfgArray ,PQPaMap ,prm -> B1 ,& interpErr );
350362 if (interpErr == 1 && * (prm -> IWork + Er16 )== 0 ){
351- #ifdef MATLAB_MEX_FILE
363+ #ifdef MATLAB_MEX_FILE
352364 printf ("Warning in %s, Error calculating Cfg. Vector definitions may need to be expanded.\n" , prm -> BlkNm );
353- #endif
365+ #endif
354366 * (prm -> IWork + Er16 ) = 1 ;
355367 }
356368 Cv = 1 ;
0 commit comments