root/branches/maintenance/3.0/Modelica/Media/Incompressible.mo

Revision 1083, 32.4 kB (checked in by dietmarw, 9 months ago)

updated the copyright notices to 2008 and corrected some Umlauts to HTML entities when present in the HTML environment (closes ticket:47)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1within Modelica.Media;
2package Incompressible
3  "Medium model for T-dependent properties, defined by tables or polynomials"
4
5  import SI = Modelica.SIunits;
6  import Cv = Modelica.SIunits.Conversions;
7  import Modelica.Constants;
8  import Modelica.Math;
9
10  package Common "Common data structures"
11
12    // Extended record for input to functions based on polynomials
13    record BaseProps_Tpoly "Fluid state record"
14      extends Modelica.Icons.Record;
15      annotation(Documentation(info="<html></html>"));
16      SI.Temperature T "temperature";
17      SI.Pressure p "pressure";
18      //    SI.Density d "density";
19    end BaseProps_Tpoly;
20
21    //     record BaseProps_Tpoly_old "fluid state record"
22    //       extends Modelica.Media.Interfaces.PartialMedium.ThermodynamicState;
23    //       //      SI.SpecificHeatCapacity cp "specific heat capacity";
24    //       SI.Temperature T "temperature";
25    //       SI.Pressure p "pressure";
26    //       //    SI.Density d "density";
27    //       parameter Real[:] poly_rho "polynomial coefficients";
28    //       parameter Real[:] poly_Cp "polynomial coefficients";
29    //       parameter Real[:] poly_eta "polynomial coefficients";
30    //       parameter Real[:] poly_pVap "polynomial coefficients";
31    //       parameter Real[:] poly_lam "polynomial coefficients";
32    //       parameter Real[:] invTK "inverse T [1/K]";
33    //     end BaseProps_Tpoly_old;
34  end Common;
35
36  package TableBased "Incompressible medium properties based on tables"
37    import Poly = Modelica.Media.Incompressible.TableBased.Polynomials_Temp;
38
39    extends Modelica.Media.Interfaces.PartialMedium(
40       final reducedX=true,
41       final fixedX = true,
42       mediumName="tableMedium",
43       redeclare record ThermodynamicState=Common.BaseProps_Tpoly,
44       singleState=true);
45    // Constants to be set in actual Medium
46    constant Boolean enthalpyOfT=true
47      "true if enthalpy is approximated as a function of T only, (p-dependence neglected)";
48    constant Boolean densityOfT = size(tableDensity,1) > 1
49      "true if density is a function of temperature";
50    constant Temperature T_min "Minimum temperature valid for medium model";
51    constant Temperature T_max "Maximum temperature valid for medium model";
52    constant Temperature T0=273.15 "reference Temperature";
53    constant SpecificEnthalpy h0=0 "reference enthalpy at T0, reference_p";
54    constant SpecificEntropy s0=0 "reference entropy at T0, reference_p";
55    constant MolarMass MM_const=0.1 "Molar mass";
56    constant Integer npol=2 "degree of polynomial used for fitting";
57    constant Integer neta=size(tableViscosity,1)
58      "number of data points for viscosity";
59    constant Real[:,:] tableDensity "Table for rho(T)";
60    constant Real[:,:] tableHeatCapacity "Table for Cp(T)";
61    constant Real[:,:] tableViscosity "Table for eta(T)";
62    constant Real[:,:] tableVaporPressure "Table for pVap(T)";
63    constant Real[:,:] tableConductivity "Table for lambda(T)";
64    //    constant Real[:] TK=tableViscosity[:,1]+T0*ones(neta) "Temperature for Viscosity";
65    constant Boolean TinK "true if T[K],Kelvin used for table temperatures";
66    constant Boolean hasDensity = not (size(tableDensity,1)==0)
67      "true if table tableDensity is present";
68    constant Boolean hasHeatCapacity = not (size(tableHeatCapacity,1)==0)
69      "true if table tableHeatCapacity is present";
70    constant Boolean hasViscosity = not (size(tableViscosity,1)==0)
71      "true if table tableViscosity is present";
72    constant Boolean hasVaporPressure = not (size(tableVaporPressure,1)==0)
73      "true if table tableVaporPressure is present";
74    final constant Real invTK[neta] = if size(tableViscosity,1) > 0 then 
75        invertTemp(tableViscosity[:,1],TinK) else fill(0,0);
76  annotation(__Dymola_keepConstant = true, Documentation(info="<HTML>
77<p>
78This is the base package for medium models of incompressible fluids based on
79tables. The minimal data to provide for a useful medium description is tables
80of density and heat capacity as functions of temperature.
81</p>
82
83<h4>Using the package TableBased</h4>
84<p>
85To implement a new medium model, create a package that <b>extends</b> TableBased
86and provides one or more of the constant tables:
87</p>
88
89<pre>
90tableDensity        = [T, d];
91tableHeatCapacity   = [T, Cp];
92tableConductivity   = [T, lam];
93tableViscosity      = [T, eta];
94tableVaporPressure  = [T, pVap];
95</pre>
96
97<p>
98The table data is used to fit constant polynomials of order <b>npol</b>, the
99temperature data points do not need to be same for different properties. Properties
100like enthalpy, inner energy and entropy are calculated consistently from integrals
101and derivatives of d(T) and Cp(T). The minimal
102data for a useful medium model is thus density and heat capacity. Transport
103properties and vapor pressure are optional, if the data tables are empty the corresponding
104function calls can not be used.
105</p>
106</HTML>"));
107    final constant Real poly_rho[:] = if hasDensity then 
108                                         Poly.fitting(tableDensity[:,1],tableDensity[:,2],npol) else 
109                                           zeros(npol+1) annotation(keepConstant = true);
110    final constant Real poly_Cp[:] = if hasHeatCapacity then 
111                                         Poly.fitting(tableHeatCapacity[:,1],tableHeatCapacity[:,2],npol) else 
112                                           zeros(npol+1) annotation(keepConstant = true);
113    final constant Real poly_eta[:] = if hasViscosity then 
114                                         Poly.fitting(invTK, Math.log(tableViscosity[:,2]),npol) else 
115                                           zeros(npol+1) annotation(keepConstant = true);
116    final constant Real poly_pVap[:] = if hasVaporPressure then 
117                                         Poly.fitting(tableVaporPressure[:,1],tableVaporPressure[:,2],npol) else 
118                                            zeros(npol+1) annotation(keepConstant= true);
119    final constant Real poly_lam[:] = if size(tableConductivity,1)>0 then 
120                                         Poly.fitting(tableConductivity[:,1],tableConductivity[:,2],npol) else 
121                                           zeros(npol+1) annotation(keepConstant = true);
122    function invertTemp "function to invert temperatures"
123      input Real[:] table "table temperature data";
124      input Boolean Tink "flag for Celsius or Kelvin";
125      output Real invTable[size(table,1)] "inverted temperatures";
126    algorithm
127      for i in 1:size(table,1) loop
128        invTable[i] := if TinK then 1/table[i] else 1/Cv.from_degC(table[i]);
129      end for;
130    end invertTemp;
131
132    redeclare model extends BaseProperties(
133      final standardOrderComponents=true,
134      p_bar=Cv.to_bar(p),
135      T_degC(start = T_start-273.15)=Cv.to_degC(T),
136      T(start = T_start,
137        stateSelect=if preferredMediumStates then StateSelect.prefer else StateSelect.default))
138      "Base properties of T dependent medium"
139    //  redeclare parameter SpecificHeatCapacity R=Modelica.Constants.R,
140
141      annotation(Documentation(info="<html>
142<p>
143Note that the inner energy neglects the pressure dependence, which is only
144true for an incompressible medium with d = constant. The neglected term is
145p-reference_p)/rho*(T/rho)*(partial rho /partial T). This is very small for
146liquids due to proportionality to 1/d^2, but can be problematic for gases that are
147modeled incompressible.
148</p>
149 
150<p>                               
151Enthalpy is never a function of T only (h = h(T) + (p-reference_p)/d), but the
152error is also small and non-linear systems can be avoided. In particular,
153non-linear systems are small and local as opposed to large and over all volumes.
154</p>                               
155   
156<p>
157Entropy is calculated as
158</p>
159<pre>
160  s = s0 + integral(Cp(T)/T,dt)
161</pre>
162<p>
163which is only exactly true for a fluid with constant density d=d0.
164</p>
165</html>
166      "));
167      SI.SpecificHeatCapacity cp "specific heat capacity";
168      parameter SI.Temperature T_start = 298.15 "initial temperature";
169    equation
170      assert(hasDensity,"Medium " + mediumName +
171                        " can not be used without assigning tableDensity.");
172      assert(T >= T_min and T <= T_max, "Temperature T (= " + String(T) +
173             " K) is not in the allowed range (" + String(T_min) +
174             " K <= T <= " + String(T_max) + " K) required from medium model \""
175             + mediumName + "\".");
176      R = Modelica.Constants.R;
177      cp = Poly.evaluate(poly_Cp,if TinK then T else T_degC);
178      h = if enthalpyOfT then h_T(T) else  h_pT(p,T,densityOfT);
179      if singleState then
180        u = h_T(T) - reference_p/d;
181      else
182        u = h - p/d;
183      end if;
184      d = Poly.evaluate(poly_rho,if TinK then T else T_degC);
185      state.T = T;
186      state.p = p;
187      MM = MM_const;
188    end BaseProperties;
189
190    redeclare function extends setState_pTX
191      "Returns state record, given pressure and temperature"
192    algorithm
193      state := ThermodynamicState(p=p,T=T);
194    end setState_pTX;
195
196    redeclare function extends setState_dTX
197      "Returns state record, given pressure and temperature"
198    algorithm
199      assert(false, "for incompressible media with d(T) only, state can not be set from density and temperature");
200    end setState_dTX;
201
202    function setState_pT "returns state record as function of p and T"
203      input AbsolutePressure p "pressure";
204      input Temperature T "temperature";
205      output ThermodynamicState state "thermodynamic state";
206    algorithm
207      state.T := T;
208      state.p := p;
209    end setState_pT;
210
211    redeclare function extends setState_phX
212      "Returns state record, given pressure and specific enthalpy"
213    algorithm
214      state :=ThermodynamicState(p=p,T=T_ph(p,h));
215    end setState_phX;
216
217    function setState_ph "returns state record as function of p and h"
218      input AbsolutePressure p "pressure";
219      input SpecificEnthalpy h "specific enthalpy";
220      output ThermodynamicState state "thermodynamic state";
221    algorithm
222      state :=ThermodynamicState(p=p,T=T_ph(p,h));
223    end setState_ph;
224
225    redeclare function extends setState_psX
226      "Returns state record, given pressure and specific entropy"
227    algorithm
228      state :=ThermodynamicState(p=p,T=T_ps(p,s));
229    end setState_psX;
230
231    function setState_ps "returns state record as function of p and s"
232      input AbsolutePressure p "pressure";
233      input SpecificEntropy s "specific entropy";
234      output ThermodynamicState state "thermodynamic state";
235    algorithm
236      state :=ThermodynamicState(p=p,T=T_ps(p,s));
237    end setState_ps;
238
239    redeclare function extends specificHeatCapacityCv
240      "Specific heat capacity at constant volume (or pressure) of medium"
241
242     annotation(smoothOrder=2);
243    algorithm
244      assert(hasHeatCapacity,"Specific Heat Capacity, Cv, is not defined for medium "
245                                             + mediumName + ".");
246      cv := Poly.evaluate(poly_Cp,if TinK then state.T else state.T - 273.15);
247    end specificHeatCapacityCv;
248
249    redeclare function extends specificHeatCapacityCp
250      "Specific heat capacity at constant volume (or pressure) of medium"
251
252     annotation(smoothOrder=2);
253    algorithm
254      assert(hasHeatCapacity,"Specific Heat Capacity, Cv, is not defined for medium "
255                                             + mediumName + ".");
256      cp := Poly.evaluate(poly_Cp,if TinK then state.T else state.T - 273.15);
257    end specificHeatCapacityCp;
258
259    redeclare function extends dynamicViscosity
260      "Return dynamic viscosity as a function of the thermodynamic state record"
261
262     annotation(smoothOrder=2);
263    algorithm
264      assert(size(tableViscosity,1)>0,"DynamicViscosity, eta, is not defined for medium "
265                                             + mediumName + ".");
266      eta := Math.exp(Poly.evaluate(poly_eta, 1/state.T));
267    end dynamicViscosity;
268
269    redeclare function extends thermalConductivity
270      "Return thermal conductivity as a function of the thermodynamic state record"
271
272     annotation(smoothOrder=2);
273    algorithm
274      assert(size(tableConductivity,1)>0,"ThermalConductivity, lambda, is not defined for medium "
275                                             + mediumName + ".");
276      lambda := Poly.evaluate(poly_lam,if TinK then state.T else Cv.to_degC(state.T));
277    end thermalConductivity;
278
279    function s_T "compute specific entropy"
280      input Temperature T "temperature";
281      output SpecificEntropy s "specific entropy";
282    algorithm
283      s := s0 + (if TinK then 
284        Poly.integralValue(poly_Cp[1:npol],T, T0) else 
285        Poly.integralValue(poly_Cp[1:npol],Cv.to_degC(T),Cv.to_degC(T0)))
286        + Modelica.Math.log(T/T0)*
287        Poly.evaluate(poly_Cp,if TinK then 0 else Modelica.Constants.T_zero);
288    end s_T;
289
290    redeclare function extends specificEntropy "Return specific entropy
291 as a function of the thermodynamic state record"
292
293    protected
294      Integer npol=size(poly_Cp,1)-1;
295     annotation(smoothOrder=2);
296    algorithm
297      assert(hasHeatCapacity,"Specific Entropy, s(T), is not defined for medium "
298                                             + mediumName + ".");
299      s := s_T(state.T);
300    end specificEntropy;
301
302    function h_T "Compute specific enthalpy from temperature"
303      import Modelica.SIunits.Conversions.to_degC;
304      extends Modelica.Icons.Function;
305      input SI.Temperature T "Temperature";
306      output SI.SpecificEnthalpy h "Specific enthalpy at p, T";
307     annotation(derivative=h_T_der);
308    algorithm
309      h :=h0 + Poly.integralValue(poly_Cp, if TinK then T else Cv.to_degC(T), if TinK then 
310      T0 else Cv.to_degC(T0));
311    end h_T;
312
313    function h_T_der "Compute specific enthalpy from temperature"
314      import Modelica.SIunits.Conversions.to_degC;
315      extends Modelica.Icons.Function;
316      input SI.Temperature T "Temperature";
317      input Real dT "temperature derivative";
318      output Real dh "derivative of Specific enthalpy at T";
319    algorithm
320      dh :=Poly.evaluate(poly_Cp, if TinK then T else Cv.to_degC(T))*dT;
321    end h_T_der;
322
323    function h_pT "Compute specific enthalpy from pressure and temperature"
324      import Modelica.SIunits.Conversions.to_degC;
325      extends Modelica.Icons.Function;
326      input SI.Pressure p "Pressure";
327      input SI.Temperature T "Temperature";
328      input Boolean densityOfT = false
329        "include or neglect density derivative dependence of enthalpy" annotation(Evaluate);
330      output SI.SpecificEnthalpy h "Specific enthalpy at p, T";
331     annotation(smoothOrder=2);
332    algorithm
333      h :=h0 + Poly.integralValue(poly_Cp, if TinK then T else Cv.to_degC(T), if TinK then 
334      T0 else Cv.to_degC(T0)) + (p - reference_p)/Poly.evaluate(poly_rho, if TinK then 
335              T else Cv.to_degC(T))
336        *(if densityOfT then (1 + T/Poly.evaluate(poly_rho, if TinK then T else Cv.to_degC(T))
337      *Poly.derivativeValue(poly_rho,if TinK then T else Cv.to_degC(T))) else 1.0);
338    end h_pT;
339
340    function density_T "Return density as function of temperature"
341
342      input Temperature T "temperature";
343      output Density d "density";
344    algorithm
345      d := Poly.evaluate(poly_rho,if TinK then T else Cv.to_degC(T));
346    end density_T;
347
348    redeclare function extends temperature
349      "Return temperature as a function of the thermodynamic state record"
350
351    algorithm
352     T := state.T;
353    end temperature;
354
355    redeclare function extends pressure
356      "Return pressure as a function of the thermodynamic state record"
357    algorithm
358     p := state.p;
359    end pressure;
360
361    redeclare function extends density
362      "Return density as a function of the thermodynamic state record"
363    algorithm
364      d := Poly.evaluate(poly_rho,if TinK then state.T else Cv.to_degC(state.T));
365    end density;
366
367    redeclare function extends specificEnthalpy
368      "Return specific enthalpy as a function of the thermodynamic state record"
369    algorithm
370      h := if enthalpyOfT then h_T(state.T) else h_pT(state.p,state.T);
371    end specificEnthalpy;
372
373    redeclare function extends specificInternalEnergy
374      "Return specific internal energy as a function of the thermodynamic state record"
375    algorithm
376      u := if enthalpyOfT then h_T(state.T) else h_pT(state.p,state.T)
377        - (if singleState then  reference_p/density(state) else state.p/density(state));
378    end specificInternalEnergy;
379
380    function T_ph "Compute temperature from pressure and specific enthalpy"
381      input AbsolutePressure p "pressure";
382      input SpecificEnthalpy h "specific enthalpy";
383      output Temperature T "temperature";
384    protected
385      package Internal
386        "Solve h(T) for T with given h (use only indirectly via temperature_phX)"
387        extends Modelica.Media.Common.OneNonLinearEquation;
388
389        redeclare record extends f_nonlinear_Data
390          "superfluous record, fix later when better structure of inverse functions exists"
391            constant Real[5] dummy = {1,2,3,4,5};
392        end f_nonlinear_Data;
393
394        redeclare function extends f_nonlinear "p is smuggled in via vector"
395        algorithm
396          y := if singleState then h_T(x) else h_pT(p,x);
397        end f_nonlinear;
398
399        // Dummy definition has to be added for current Dymola
400        redeclare function extends solve
401        end solve;
402      end Internal;
403    algorithm
404     T := Internal.solve(h, T_min, T_max, p, {1}, Internal.f_nonlinear_Data());
405    end T_ph;
406
407    function T_ps "Compute temperature from pressure and specific enthalpy"
408      input AbsolutePressure p "pressure";
409      input SpecificEntropy s "specific entropy";
410      output Temperature T "temperature";
411    protected
412      package Internal
413        "Solve h(T) for T with given h (use only indirectly via temperature_phX)"
414        extends Modelica.Media.Common.OneNonLinearEquation;
415
416        redeclare record extends f_nonlinear_Data
417          "superfluous record, fix later when better structure of inverse functions exists"
418            constant Real[5] dummy = {1,2,3,4,5};
419        end f_nonlinear_Data;
420
421        redeclare function extends f_nonlinear "p is smuggled in via vector"
422        algorithm
423          y := s_T(x);
424        end f_nonlinear;
425
426        // Dummy definition has to be added for current Dymola
427        redeclare function extends solve
428        end solve;
429      end Internal;
430    algorithm
431     T := Internal.solve(s, T_min, T_max, p, {1}, Internal.f_nonlinear_Data());
432    end T_ps;
433
434    package Polynomials_Temp
435      "Temporary Functions operating on polynomials (including polynomial fitting); only to be used in Modelica.Media.Incompressible.TableBased"
436      extends Modelica.Icons.Library;
437
438      annotation (Documentation(info="<HTML>
439<p>
440This package contains functions to operate on polynomials,
441in particular to determine the derivative and the integral
442of a polynomial and to use a polynomial to fit a given set
443of data points.
444</p>
445<p>
446
447<p><b>Copyright &copy; 2004-2008, Modelica Association and DLR.</b></p>
448
449<p><i>
450This package is <b>free</b> software. It can be redistributed and/or modified
451under the terms of the <b>Modelica license</b>, see the license conditions
452and the accompanying <b>disclaimer</b> in the documentation of package
453Modelica in file \"Modelica/package.mo\".
454</i>
455</p>
456
457</HTML>
458",     revisions="<html>
459<ul>
460<li><i>Oct. 22, 2004</i> by Martin Otter (DLR):<br>
461       Renamed functions to not have abbrevations.<br>
462       Based fitting on LAPACK<br>
463       New function to return the polynomial of an indefinite integral<li>
464<li><i>Sept. 3, 2004</i> by Jonas Eborn (Scynamics):<br>
465       polyderval, polyintval added<li>
466<li><i>March 1, 2004</i> by Martin Otter (DLR):<br>
467       first version implemented
468</li>
469</ul>
470</html>"));
471      function evaluate "Evaluate polynomial at a given abszissa value"
472        extends Modelica.Icons.Function;
473        input Real p[:]
474          "Polynomial coefficients (p[1] is coefficient of highest power)";
475        input Real u "Abszissa value";
476        output Real y "Value of polynomial at u";
477        annotation(derivative(zeroDerivative=p)=evaluate_der);
478      algorithm
479        y := p[1];
480        for j in 2:size(p, 1) loop
481          y := p[j] + u*y;
482        end for;
483      end evaluate;
484
485      function derivative "Derivative of polynomial"
486        extends Modelica.Icons.Function;
487        input Real p1[:]
488          "Polynomial coefficients (p1[1] is coefficient of highest power)";
489        output Real p2[size(p1, 1) - 1] "Derivative of polynomial p1";
490      protected
491        Integer n=size(p1, 1);
492      algorithm
493        for j in 1:n-1 loop
494          p2[j] := p1[j]*(n - j);
495        end for;
496      end derivative;
497
498      function derivativeValue
499        "Value of derivative of polynomial at abszissa value u"
500        extends Modelica.Icons.Function;
501        input Real p[:]
502          "Polynomial coefficients (p[1] is coefficient of highest power)";
503        input Real u "Abszissa value";
504        output Real y "Value of derivative of polynomial at u";
505        annotation(derivative(zeroDerivative=p)=derivativeValue_der);
506      protected
507        Integer n=size(p, 1);
508      algorithm
509        y := p[1]*(n - 1);
510        for j in 2:size(p, 1)-1 loop
511          y := p[j]*(n - j) + u*y;
512        end for;
513      end derivativeValue;
514
515      function secondDerivativeValue
516        "Value of 2nd derivative of polynomial at abszissa value u"
517        extends Modelica.Icons.Function;
518        input Real p[:]
519          "Polynomial coefficients (p[1] is coefficient of highest power)";
520        input Real u "Abszissa value";
521        output Real y "Value of 2nd derivative of polynomial at u";
522      protected
523        Integer n=size(p, 1);
524      algorithm
525        y := p[1]*(n - 1)*(n - 2);
526        for j in 2:size(p, 1)-2 loop
527          y := p[j]*(n - j)*(n - j - 1) + u*y;
528        end for;
529      end secondDerivativeValue;
530
531      function integral "Indefinite integral of polynomial p(u)"
532        extends Modelica.Icons.Function;
533        input Real p1[:]
534          "Polynomial coefficients (p1[1] is coefficient of highest power)";
535        output Real p2[size(p1, 1) + 1]
536          "Polynomial coefficients of indefinite integral of polynomial p1 (polynomial p2 + C is the indefinite integral of p1, where C is an arbitrary constant)";
537      protected
538        Integer n=size(p1, 1) + 1 "degree of output polynomial";
539      algorithm
540        for j in 1:n-1 loop
541          p2[j] := p1[j]/(n-j);
542        end for;
543      end integral;
544
545      function integralValue "Integral of polynomial p(u) from u_low to u_high"
546        extends Modelica.Icons.Function;
547        input Real p[:] "Polynomial coefficients";
548        input Real u_high "High integrand value";
549        input Real u_low=0 "Low integrand value, default 0";
550        output Real integral=0.0
551          "Integral of polynomial p from u_low to u_high";
552        annotation(derivative(zeroDerivative=p)=integralValue_der);
553      protected
554        Integer n=size(p, 1) "degree of integrated polynomial";
555        Real y_low=0 "value at lower integrand";
556      algorithm
557        for j in 1:n loop
558          integral := u_high*(p[j]/(n - j + 1) + integral);
559          y_low := u_low*(p[j]/(n - j + 1) + y_low);
560        end for;
561        integral := integral - y_low;
562      end integralValue;
563
564      function fitting
565        "Computes the coefficients of a polynomial that fits a set of data points in a least-squares sense"
566        extends Modelica.Icons.Function;
567        input Real u[:] "Abscissa data values";
568        input Real y[size(u, 1)] "Ordinate data values";
569        input Integer n(min=1)
570          "Order of desired polynomial that fits the data points (u,y)";
571        output Real p[n + 1]
572          "Polynomial coefficients of polynomial that fits the date points";
573        annotation (Documentation(info="<HTML>
574<p>
575Polynomials.fitting(u,y,n) computes the coefficients of a polynomial
576p(u) of degree \"n\" that fits the data \"p(u[i]) - y[i]\"
577in a least squares sense. The polynomial is
578returned as a vector p[n+1] that has the following definition:
579</p>
580<pre>
581  p(u) = p[1]*u^n + p[2]*u^(n-1) + ... + p[n]*u + p[n+1];
582</pre>
583</HTML>"));
584      protected
585        Real V[size(u, 1), n + 1] "Vandermonde matrix";
586      algorithm
587        // Construct Vandermonde matrix
588        V[:, n + 1] := ones(size(u, 1));
589        for j in n:-1:1 loop
590          V[:, j] := {u[i] * V[i, j + 1] for i in 1:size(u,1)};
591        end for;
592
593        // Solve least squares problem
594        p :=Modelica.Math.Matrices.leastSquares(V, y);
595      end fitting;
596
597      function evaluate_der "Evaluate polynomial at a given abszissa value"
598        extends Modelica.Icons.Function;
599        input Real p[:]
600          "Polynomial coefficients (p[1] is coefficient of highest power)";
601        input Real u "Abszissa value";
602        input Real du "Abszissa value";
603        output Real dy "Value of polynomial at u";
604      protected
605        Integer n=size(p, 1);
606      algorithm
607        dy := p[1]*(n - 1);
608        for j in 2:size(p, 1)-1 loop
609          dy := p[j]*(n - j) + u*dy;
610        end for;
611        dy := dy*du;
612      end evaluate_der;
613
614      function integralValue_der
615        "Time derivative of integral of polynomial p(u) from u_low to u_high, assuming only u_high as time-dependent (Leibnitz rule)"
616        extends Modelica.Icons.Function;
617        input Real p[:] "Polynomial coefficients";
618        input Real u_high "High integrand value";
619        input Real u_low=0 "Low integrand value, default 0";
620        input Real du_high "High integrand value";
621        input Real du_low=0 "Low integrand value, default 0";
622        output Real dintegral=0.0
623          "Integral of polynomial p from u_low to u_high";
624      algorithm
625        dintegral := evaluate(p,u_high)*du_high;
626      end integralValue_der;
627
628      function derivativeValue_der
629        "time derivative of derivative of polynomial"
630        extends Modelica.Icons.Function;
631        input Real p[:]
632          "Polynomial coefficients (p[1] is coefficient of highest power)";
633        input Real u "Abszissa value";
634        input Real du "delta of abszissa value";
635        output Real dy
636          "time-derivative of derivative of polynomial w.r.t. input variable at u";
637      protected
638        Integer n=size(p, 1);
639      algorithm
640        dy := secondDerivativeValue(p,u)*du;
641      end derivativeValue_der;
642    end Polynomials_Temp;
643  annotation (
644    Documentation(info="<HTML>
645<h4>Table based media</h4>
646<p>
647This is the base package for medium models of incompressible fluids based on
648tables. The minimal data to provide for a useful medium description is tables
649of density and heat capacity as functions of temperature.
650</p>
651<h4>Using the package TableBased</h4>
652<p>
653To implement a new medium model, create a package that <b>extends</b> TableBased
654and provides one or more of the constant tables:
655<pre>
656tableDensity        = [T, d];
657tableHeatCapacity   = [T, Cp];
658tableConductivity   = [T, lam];
659tableViscosity      = [T, eta];
660tableVaporPressure  = [T, pVap];
661</pre>
662The table data is used to fit constant polynomials of order <b>npol</b>, the
663temperature data points do not need to be same for different properties. Properties
664like enthalpy, inner energy and entropy are calculated consistently from integrals
665and derivatives of d(T) and Cp(T). The minimal
666data for a useful medium model is thus density and heat capacity. Transport
667properties and vapor pressure are optional, if the data tables are empty the corresponding
668function calls can not be used.
669</HTML>"));
670
671  end TableBased;
672