TrioCFD 1.9.9_beta
TrioCFD documentation
Loading...
Searching...
No Matches
Lois_sodium.h
1/****************************************************************************
2* Copyright (c) 2023, CEA
3* All rights reserved.
4*
5* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
6* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
7* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
8* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
9*
10* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
11* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
12* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
13*
14*****************************************************************************/
15
16#ifndef lois_sodium_span_included
17#define lois_sodium_span_included
18
19#include <assert.h>
20#include <span.hpp>
21#include <vector>
22#include <cmath>
23
24/* Sodium physical laws derived from fits in DEN/DANS/DM2S/STMF/LMES/RT/12-018/A */
25
26// iterator index !
27#define i_it std::distance(res.begin(), &val)
28
29using VectorD = std::vector<double>;
30using SpanD = tcb::span<double>;
31
32const double Ksil = 1e-9;
33
34/* inverse of the liquid density (pressure-free) */
35inline void IRhoL(const SpanD T, SpanD res, int ncomp, int id)
36{
37 for (auto& val : res)
38 {
39 const double tk = T[i_it * ncomp + id] + 273.15;
40 val = (9.829728e-4 + tk * (2.641186e-7 + tk * (-3.340743e-11 + tk * 6.680973e-14)));
41 }
42}
43
44/* its derivative */
45inline void DTIRhoL(const SpanD T, SpanD res, int ncomp, int id)
46{
47 for (auto& val : res)
48 {
49 const double tk = T[i_it * ncomp + id] + 273.15;
50 val = (2.641186e-7 + tk * (-3.340743e-11 * 2 + tk * 6.680973e-14 * 3));
51 }
52}
53
54inline void RhoL(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
55{
56 assert ((int)T.size() == ncomp * (int)P.size() && (int)T.size() == ncomp * (int)res.size());
57 IRhoL(T,res,ncomp,id); // Note: res is being used here ...
58 for (auto& val : res) val = exp(Ksil * (P[i_it] - 1e5)) / val;
59}
60
61/* derivatives */
62inline void DPRhoL(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
63{
64 assert ((int)T.size() == ncomp * (int)P.size() && (int)T.size() == ncomp * (int)res.size());
65 RhoL(T,P,res,ncomp,id); // Note: res is being used here ...
66 for (auto& val : res) val = Ksil * val;
67}
68
69inline void DTRhoL(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
70{
71 assert ((int)T.size() == ncomp * (int)P.size() && (int)T.size() == ncomp * (int)res.size());
72 VectorD invRhoL((int)res.size()), DinvRhoL((int)res.size()); // res could be used for one of the two but leaving as is ...
73 IRhoL(T,SpanD(invRhoL),ncomp,id);
74 DTIRhoL(T,SpanD(DinvRhoL),ncomp,id);
75 for (auto& val : res) val = - exp(Ksil * (P[i_it] - 1e5)) * DinvRhoL[i_it] / (invRhoL[i_it] * invRhoL[i_it]);
76}
77
78/* liquid enthalpy */
79inline void HL(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
80{
81 assert ((int)T.size() == ncomp * (int)P.size() && (int)T.size() == ncomp * (int)res.size());
82 VectorD invRhoL((int)res.size()), DinvRhoL((int)res.size()); // res could be used for one of the two but leaving as is ...
83 IRhoL(T,SpanD(invRhoL),ncomp,id);
84 DTIRhoL(T,SpanD(DinvRhoL),ncomp,id);
85 for (auto& val : res)
86 {
87 const double tk = T[i_it * ncomp + id] + 273.15;
88 val = (2992600 / tk - 365487.2 + tk * (1658.2 + tk * (-.42395 + tk * 1.4847e-4)))+ (invRhoL[i_it] - tk * DinvRhoL[i_it]) * (1 - exp(Ksil*(1e5 - P[i_it]))) / Ksil;
89 }
90}
91
92/* derivatives */
93inline void DPHL(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
94{
95 assert ((int)T.size() == ncomp * (int)P.size() && (int)T.size() == ncomp * (int)res.size());
96 VectorD invRhoL((int)res.size()), DinvRhoL((int)res.size()); // res could be used for one of the two but leaving as is ...
97 IRhoL(T,SpanD(invRhoL),ncomp,id);
98 DTIRhoL(T,SpanD(DinvRhoL),ncomp,id);
99 for (auto& val : res)
100 {
101 const double tk = T[i_it * ncomp + id] + 273.15;
102 val = (invRhoL[i_it] - tk * DinvRhoL[i_it]) * exp(Ksil * (1e5 - P[i_it]));
103 }
104}
105
106inline void DTHL(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
107{
108 assert ((int)T.size() == ncomp * (int)res.size() && (int)T.size() == ncomp * (int)P.size());
109 for (auto& val : res)
110 {
111 const double tk = T[i_it * ncomp + id] + 273.15;
112 val = (-2992600 / (tk*tk) + 1658.2 + tk * (-.42395 * 2 + tk * 1.4847e-4 * 3))- tk * (-3.340743e-11 * 2 + tk * 6.680973e-14 * 6) * (1 - exp(Ksil*(1e5 - P[i_it]))) / Ksil;
113 }
114}
115
116inline void BETAL(const SpanD T, SpanD res, int ncomp, int id)
117{
118 assert ((int)T.size() == ncomp * (int)res.size());
119 VectorD invRhoL((int)res.size()), DinvRhoL((int)res.size()); // res could be used for one of the two but leaving as is ...
120 IRhoL(T,SpanD(invRhoL),ncomp,id);
121 DTIRhoL(T,SpanD(DinvRhoL),ncomp,id);
122 for (auto& val : res) val = DinvRhoL[i_it] / invRhoL[i_it];
123}
124
125/* liquid dynamic viscosity */
126inline void MuL(const SpanD T, SpanD res, int ncomp, int id)
127{
128 assert ((int)T.size() == ncomp * (int)res.size());
129 for (auto& val : res)
130 {
131 const double tk = T[i_it * ncomp + id] + 273.15;
132 val = exp( -6.4406 - 0.3958 * log(tk) + 556.835 / tk );
133 }
134}
135
136/* liquid thermal conductivity */
137inline void LambdaL(const SpanD T, SpanD res, int ncomp, int id)
138{
139 assert ((int)T.size() == ncomp * (int)res.size());
140 for (auto& val : res)
141 {
142 const double tk = T[i_it * ncomp + id] + 273.15;
143 val = std::max(124.67 + tk * (-.11381 + tk * (5.5226e-5 - tk * 1.1848e-8)), 0.045);
144 }
145}
146
147/* surface tension */
148inline void SigmaL(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
149{
150 assert ((int)P.size() == (int)res.size());
151 for (auto& val : res)
152 {
153 const double tk = T[i_it * ncomp + id] + 273.15;
154 val = 0.2405 * pow(1 - tk / 2503.7, 1.126);
155 }
156}
157
158/*
159 * ****************************
160 * ********** GAZ ***********
161 * ****************************
162 */
163
164/* saturation temperature */
165inline void Tsat_Na(const SpanD P, SpanD res, int ncomp, int ind)
166{
167 assert (ncomp * (int)P.size() == (int)res.size());
168 const double A = 7.8130, B = 11209, C = 5.249e5;
169 for (int i =0; i < (int)P.size(); i++) res[i * ncomp + ind] = 2 * C / ( -B + sqrt(B*B + 4 * A * C - 4 * C * log(P[i] / 1e6))) - 273.15;
170}
171
172/* its derivative */
173inline void DTsat_Na(const SpanD P, SpanD res, int ncomp, int ind)
174{
175 assert (ncomp * (int)P.size() == (int)res.size());
176 const double A = 7.8130, B = 11209, C = 5.249e5;
177 Tsat_Na(P,res,ncomp,ind); // XXX: if it gets complicated we will need to create a local array ...
178 for (int i =0; i < (int)P.size(); i++)
179 {
180 const double Tsk = res[i * ncomp + ind] + 273.15; // XXX : need in Kelvin
181 res[i * ncomp + ind] = Tsk * Tsk / ( P[i] * sqrt(B*B + 4 * A * C - 4 * C * log(P[i] / 1e6)));
182 }
183}
184
185/* saturation vapor pressure */
186inline void Psat_Na(const SpanD T, SpanD res, int ncomp, int ind)
187{
188 assert ((int)T.size() == (int)res.size());
189 const double A = 7.8130, B = 11209, C = 5.249e5;
190 for (int i =0; i < (int)T.size() / ncomp; i++)
191 {
192 const double tk = T[i * ncomp + ind] + 273.15;
193 res[i * ncomp + ind] = 1.e6 * exp(A - B / tk - C / (tk * tk));
194 }
195}
196
197/* its derivative */
198inline void DPsat_Na(const SpanD T, SpanD res, int ncomp, int ind)
199{
200 assert ((int)T.size() == (int)res.size());
201 const double B = 11209, C = 5.249e5;
202 // note: resu is being used here!
203 Psat_Na(T,res,ncomp,ind);
204 for (int i =0; i < (int)T.size() / ncomp; i++)
205 {
206 const double tk = T[i * ncomp + ind] + 273.15;
207 res[i * ncomp + ind] = (B / (tk * tk) + 2 * C / (tk * tk * tk)) * res[i * ncomp + ind];
208 }
209}
210
211/* saturation specific enthalpy */
212inline void Hsat(const SpanD P, SpanD res, int ncomp, int ind)
213{
214 assert (ncomp * (int)P.size() == (int)res.size());
215 VectorD Ts_((int)P.size()), HL_((int)P.size());
216 Tsat_Na(P,SpanD(Ts_),1,0); // 1D
217 HL(SpanD(Ts_),P,SpanD(HL_),1,0); // 1D
218 for (int i =0; i < (int)P.size(); i++) res[i * ncomp + ind] = HL_[i];
219}
220
221/* its derivative */
222inline void DHsat(const SpanD P, SpanD res, int ncomp, int ind)
223{
224 assert (ncomp * (int)P.size() == (int)res.size());
225 VectorD Ts_((int)P.size()), DTs_((int)P.size()), DTHL_((int)P.size()), DPHL_((int)P.size());
226 Tsat_Na(P,SpanD(Ts_),1,0); // 1D
227 DTsat_Na(P,SpanD(DTs_),1,0); // 1D
228 DTHL(SpanD(Ts_),P,SpanD(DTHL_),1,0);
229 DPHL(SpanD(Ts_),P,SpanD(DPHL_),1,0);
230 for (int i =0; i < (int)P.size(); i++) res[i * ncomp + ind] = DTs_[i] * DTHL_[i] + DPHL_[i];
231}
232
233/* vapor density: (ideal gas) * f1(Tsat_Na) * f2(DTsat_Na) */
234#define f1 (2.49121 + Tsk * (-5.53796e-3 + Tsk * (7.5465e-6 + Tsk * (-4.20217e-9 + Tsk * 8.59212e-13))))
235#define f2 (1 + dT * (-5e-4 + dT * (6.25e-7 - dT * 4.1359e-25)))
236#define Df1 (-5.53796e-3 + Tsk * (2 * 7.5465e-6 + Tsk * (-3 * 4.20217e-9 + Tsk * 4 * 8.59212e-13)))
237#define Df2 (-5e-4 + dT * (2 * 6.25e-7 - dT * 3 * 4.1359e-25))
238
239inline void RhoV(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
240{
241 assert ((int)T.size() == ncomp * (int)P.size() && (int)T.size() == ncomp * (int)res.size());
242 Tsat_Na(P,res,1,0); // XXX: note res is being used here
243 for (auto& val : res)
244 {
245 const double Tsk = val + 273.15; // XXX : need in Kelvin
246 const double tk = T[i_it * ncomp + id] + 273.15;
247 const double dT = tk - Tsk;
248 val = P[i_it] * 2.7650313e-3 * f1 * f2 / tk;
249 }
250}
251
252/* its derivatives */
253inline void DPRhoV(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
254{
255 assert ((int)T.size() == ncomp * (int)P.size() && (int)T.size() == ncomp * (int)res.size());
256 VectorD Tsk_((int)res.size()), DTsk_((int)res.size()); // res could be used for one of the two but leaving as is ...
257 Tsat_Na(P,SpanD(Tsk_),1,0);
258 DTsat_Na(P,SpanD(DTsk_),1,0);
259 for (auto& val : res)
260 {
261 const double Tsk = Tsk_[i_it] + 273.15; // XXX : need in Kelvin
262 const double tk = T[i_it * ncomp + id] + 273.15;
263 const double dT = tk - Tsk;
264 val = 2.7650313e-3 * (f1 * f2 + P[i_it] * DTsk_[i_it] * (Df1 * f2 - f1 * Df2)) / tk;
265 }
266}
267
268inline void DTRhoV(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
269{
270 assert ((int)T.size() == ncomp * (int)P.size() && (int)T.size() == ncomp * (int)res.size());
271 Tsat_Na(P,res,1,0); // XXX: note res is being used here
272 for (auto& val : res)
273 {
274 const double Tsk = val + 273.15; // XXX : need in Kelvin
275 const double tk = T[i_it * ncomp + id] + 273.15;
276 const double dT = tk - Tsk;
277 val = P[i_it] * 2.7650313e-3 * f1 * (Df2 - f2 / tk) / tk;
278 }
279}
280
281#undef f1
282#undef Df1
283#undef f2
284#undef Df2
285
286/* latent heat, to be evaluated at Tsat_Na */
287inline void Lvap_Na(const SpanD P, SpanD res, int ncomp, int ind)
288{
289 assert (ncomp * (int)P.size() == (int)res.size());
290 Tsat_Na(P,res,ncomp,ind); // XXX: note res is being used here. XXX : if this gets complicated we will need to create a local array ...
291 const double Tc = 2503.7;
292 for (int i =0; i < (int)P.size(); i++)
293 {
294 const double Tsk = res[i * ncomp + ind] + 273.15; // XXX : need in Kelvin
295 res[i * ncomp + ind] = 3.9337e5 * ( 1 - Tsk / Tc) + 4.3986e6 * pow( 1 - Tsk / Tc, .29302);
296 }
297}
298
299/* its derivative */
300inline void DLvap_Na(const SpanD P, SpanD res, int ncomp, int ind)
301{
302 assert (ncomp * (int)P.size() == (int)res.size());
303 const double Tc = 2503.7;
304 VectorD Tsk_((int)P.size()), DTsk_((int)P.size());
305 Tsat_Na(P,SpanD(Tsk_),1,0); // 1D
306 DTsat_Na(P,SpanD(DTsk_),1,0); // 1D
307 for (int i =0; i < (int)P.size(); i++)
308 {
309 const double Tsk = Tsk_[i] + 273.15; // XXX : need in Kelvin
310 res[i * ncomp + ind] = DTsk_[i] * (-3.9337e5 / Tc - 4.3986e6 * .29302 * pow( 1 - Tsk / Tc, .29302 - 1) / Tc);
311 }
312}
313
314/* saturation specific enthalpy (vapor = Hsat(P) + Lvap_Na(P)) */
315inline void HVsat(const SpanD P, SpanD res, int ncomp, int ind)
316{
317 assert (ncomp * (int)P.size() == (int)res.size());
318 VectorD Hs_((int)P.size()), Lv_((int)P.size());
319 Hsat(P,SpanD(Hs_),1,0);
320 Lvap_Na(P,SpanD(Lv_),1,0);
321 for (int i =0; i < (int)P.size(); i++) res[i * ncomp + ind] = Hs_[i] + Lv_[i];
322}
323
324/* its derivative (vapor = DHsat(P) + DLvap_Na(P)) */
325inline void DHVsat(const SpanD P, SpanD res, int ncomp, int ind)
326{
327 assert (ncomp * (int)P.size() == (int)res.size());
328 VectorD DHs_((int)P.size()), DLv_((int)P.size());
329 DHsat(P,SpanD(DHs_),1,0);
330 DLvap_Na(P,SpanD(DLv_),1,0);
331 for (int i =0; i < (int)P.size(); i++) res[i * ncomp + ind] = DHs_[i] + DLv_[i];
332}
333
334/* vapor enthalpy */
335inline void HV(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
336{
337 assert ((int)T.size() == ncomp * (int)P.size() && (int)T.size() == ncomp * (int)res.size());
338 const double Cp0 = 910, k = -4.6e-3;
339 VectorD Ts_((int)res.size()), Lvap_((int)res.size()); // res could be used for one of the two but leaving as is ...
340 Tsat_Na(P,SpanD(Ts_),1,0);
341 Lvap_Na(P,SpanD(Lvap_),1,0);
342 HL(SpanD(Ts_),P,res,1,0); // XXX: note res is being used here
343 for (auto& val : res)
344 {
345 const double dT = T[i_it * ncomp + id] - Ts_[i_it]; // XXX : not the same size, be careful
346 const double CpVs = (-1223.89 + Ts_[i_it] * (14.1191 + Ts_[i_it] * (-1.62025e-2 + Ts_[i_it] * 5.76923e-6)));
347 val = val + Lvap_[i_it] + Cp0 * dT + (Cp0 - CpVs) * (1 - exp(k * dT)) / k;
348 }
349}
350
351/* its derivatives */
352inline void DPHV(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
353{
354 assert ((int)T.size() == ncomp * (int)P.size() && (int)T.size() == ncomp * (int)res.size());
355 VectorD DTHL_((int)res.size()), DPHL_((int)res.size()), Ts_((int)res.size()), DTs_((int)res.size()), DLvap_((int)res.size()); // res could be used for one of the two but leaving as is ...
356 Tsat_Na(P,SpanD(Ts_),1,0);
357 DTsat_Na(P,SpanD(DTs_),1,0);
358 DLvap_Na(P,SpanD(DLvap_),1,0);
359 DTHL(SpanD(Ts_),P,SpanD(DTHL_),1,0);
360 DPHL(SpanD(Ts_),P,SpanD(DPHL_),1,0);
361
362 const double Cp0 = 910, k = -4.6e-3;
363 for (auto& val : res)
364 {
365 const double dT = T[i_it * ncomp + id] - Ts_[i_it]; // XXX : not the same size, be careful
366 const double CpVs = (-1223.89 + Ts_[i_it] * (14.1191 + Ts_[i_it] * (-1.62025e-2 + Ts_[i_it] * 5.76923e-6))); // TODO : FIXME : should be a method
367 const double DCpVs = (14.1191 + Ts_[i_it] * (- 2 * 1.62025e-2 + Ts_[i_it] * 3 * 5.76923e-6));
368 val = DPHL_[i_it] + DLvap_[i_it] + DTs_[i_it] * (DTHL_[i_it] - Cp0 - DCpVs * (1 - exp(k * dT)) / k + (Cp0 - CpVs) * exp(k * dT));
369 }
370}
371
372inline void DTHV(const SpanD T, const SpanD P, SpanD res, int ncomp, int id)
373{
374 assert ((int)T.size() == ncomp * (int)res.size() && (int)T.size() == ncomp * (int)P.size());
375 const double Cp0 = 910, k = -4.6e-3;
376 Tsat_Na(P,res,1,0); // XXX: note res is being used here
377 for (auto& val : res)
378 {
379 const double dT = T[i_it * ncomp + id] - val; // XXX : not the same size, be careful : val = res[i_it]
380 const double CpVs = (-1223.89 + val * (14.1191 + val * (-1.62025e-2 + val * 5.76923e-6)));
381 val = Cp0 - (Cp0 - CpVs) * exp(k * dT);
382 }
383}
384
385/* inverse of the vapor density */
386inline void IRhoV_vec(SpanD T, const SpanD P, SpanD res, int ncomp, int id)
387{
388 RhoV(T,P,res,ncomp,id); // XXX: note res is being used here
389 for (auto& val : res) val = 1. / val;
390}
391/* its derivatives */
392inline void DTIRhoV_vec(SpanD T, const SpanD P, SpanD res, int ncomp, int id)
393{
394 VectorD RhoV_((int)res.size()), DTRhoV_((int)res.size()); // res could be used for one of the two but leaving as is ...
395 RhoV(T,P,SpanD(RhoV_),ncomp,id);
396 DTRhoV(T,P,SpanD(DTRhoV_),ncomp,id);
397 for (auto& val : res) val = -DTRhoV_[i_it] / ( RhoV_[i_it] * RhoV_[i_it] );
398}
399
400inline void DPIRhoV_vec(SpanD T, const SpanD P, SpanD res, int ncomp, int id)
401{
402 VectorD RhoV_((int)res.size()), DPRhoV_((int)res.size()); // res could be used for one of the two but leaving as is ...
403 RhoV(T,P,SpanD(RhoV_),ncomp,id);
404 DPRhoV(T,P,SpanD(DPRhoV_),ncomp,id);
405 for (auto& val : res) val = -DPRhoV_[i_it] / ( RhoV_[i_it] * RhoV_[i_it] );
406}
407
408inline void BETAV(SpanD T, const SpanD P, SpanD res, int ncomp, int id)
409{
410 assert ((int)T.size() == ncomp * (int)res.size() && (int)T.size() == ncomp * (int)P.size());
411 VectorD IRhoV_((int)res.size()), DTIRhoV_((int)res.size()); // res could be used for one of the two but leaving as is ...
412 IRhoV_vec(T,P,SpanD(IRhoV_),ncomp,id);
413 DTIRhoV_vec(T,P,SpanD(DTIRhoV_),ncomp,id);
414 for (auto& val : res) val = DTIRhoV_[i_it] / IRhoV_[i_it];
415}
416
417/* vapor dynamic viscosity */
418inline void MuV(const SpanD T, SpanD res, int ncomp, int id)
419{
420 assert ((int)T.size() == ncomp * (int)res.size());
421 for (auto& val : res) val = 2.5e-6 + 1.5e-8 * (T[i_it * ncomp + id] + 273.15);
422}
423
424/* vapor thermal conductivity */
425inline void LambdaV(const SpanD T, SpanD res, int ncomp, int id)
426{
427 assert ((int)T.size() == ncomp * (int)res.size());
428 for (auto& val : res) val = 0.045;
429}
430
431#undef i_it
432#endif /* lois_sodium_span_included */
433
434//// TODO : FIXME : to rewrite using span ... not used for now
435///* Sodium physical laws derived from fits in DEN/DANS/DM2S/STMF/LMES/RT/12-018/A */
436//#define Tk ((T) + 273.15)
437//const double Ksil_d = 1e-9;
438///* inverses by Newton's method */
439//inline double TLh(double h, double T0, double P)
440//{
441// double T = T0, sec;
442// int i;
443// for (i = 0; i < 100 && std::abs( sec = HL(T, P) - h ) > 1e-8; i++)
444// T -= sec / DTHL(T, P);
445// return i < 100 ? T : -1e10;
446//}
447///* volumetric energy of the liquid */
448//inline double EL(double T, double P) { return HL(T, P) * RhoL(T, P); }
449///* derivatives */
450//inline double DTEL(double T, double P) { return DTHL(T, P) * RhoL(T, P) + HL(T, P) * DTRhoL(T, P); }
451//inline double DPEL(double T, double P) { return DPHL(T, P) * RhoL(T, P) + HL(T, P) * DPRhoL(T, P); }
452//inline double TLe(double e, double T0, double P)
453//{
454// double T = T0, sec;
455// int i;
456// for (i = 0; i < 100 && std::abs( sec = EL(T, P) - e ) > 1e-3; i++)
457// T -= sec / DTEL(T, P);
458// return i < 100 ? T : -1e10;
459//}
460///* Nusselt number of the liquid -> Skupinski */
461//inline double NuL(double T, double P, double w, double Dh)
462//{
463// double Pe = RhoL(T, P) * std::abs(w) * Dh * DTHL(T, P) / LambdaL(T);
464// return 4.82 + 0.0185 * pow(Pe, 0.827);
465//}
466//
467///* volumetric energy of the vapor */
468//inline double EV( double T, double P) { return RhoV(T, P) * HV(T, P); }
469///* its derivatives */
470//inline double DTEV(double T, double P) { return DTRhoV(T, P) * HV(T, P) + RhoV(T, P) * DTHV(T, P); }
471//inline double DPEV(double T, double P) { return DPRhoV(T, P) * HV(T, P) + RhoV(T, P) * DPHV(T, P); }
472//
473///* Nusselt number of the vapor -> Dittus/Boetler */
474//inline double NuV(double T, double P, double w, double Dh)
475//{
476// double Re = RhoV(T, P) * std::abs(w) * Dh / MuV(T), Pr = MuV(T) * DTHV(T, P) / LambdaV(T);
477// return 0.023 * pow(Re, 0.8) * pow(Pr, 0.4);
478//}
479//
480//inline double Rho( double T, double x, double P ) { return 1 / ( (1-x)/RhoL(T, P) + x / RhoV(T, P) ); }