rewrite some derivatives

modelica-3rdparty · Dec 17, 2013 · 3f6c4d4 · 3f6c4d4
1 parent f55bb0b
commit 3f6c4d4
Show file tree

Hide file tree

Showing 2 changed files with 91 additions and 8 deletions.
diff --git a/Extras/LaTeX/tables/RewriteDerivatives.m b/Extras/LaTeX/tables/RewriteDerivatives.m
@@ -0,0 +1,77 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Rewrite derivatives in terms of few basic derivatives
+% written by Matthis Thorade
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+clear;
+clc;
+
+% independent vars and base props
+syms D T P;
+% three basic first-order derivatives (wrt d and T)
+syms dP_D dP_T dU_T;
+% four basic second-order derivatives (wrt d and T)
+syms d2P_D2 d2P_T2 d2P_TD d2U_T2;
+
+% step 1: rewrite the first-order derivatives (wrt d and T)
+dU_D = -T/D^2*dP_T+P/D^2
+dS_T = 1/T*dU_T
+dS_D = -1/D^2*dP_T
+dH_T = dU_T+1/D*dP_T
+dH_D = 1/D*dP_D - T/D^2*dP_T
+
+% step 2: rewrite the second-order derivatives (wrt d and T)
+d2S_T2 = 1/T*d2U_T2 - 1/T^2*dU_T
+d2S_D2 = -1/D^2*d2P_TD + 2/D^3*dP_T
+d2S_TD = -1/D^2*d2P_T2
+
+
+%% engineering derivs
+dH_P_T = dH_D/dP_D
+dH_P_T = simplify(dH_P_T)
+dH_P_T = expand(dH_P_T)
+
+dH_T_P = dH_T -dH_D*dP_T/dP_D
+dH_T_P = simplify(dH_T_P)
+dH_T_P = expand(dH_T_P)
+
+dP_D_S = dP_D-dP_T*dS_D/dS_T
+dP_D_S = simplify(dP_D_S)
+dP_D_S = expand(dP_D_S)
+
+
+%% dynamic derivs (first)
+dT_P_H = 1/(dP_T - dP_D*dH_T/dH_D)
+dT_P_H = simplify(dT_P_H)
+dT_P_H = expand(dT_P_H)
+
+dT_H_P = 1/(dH_T - dH_D*dP_T/dP_D)
+dT_H_P = simplify(dT_H_P)
+dT_H_P = expand(dT_H_P)
+
+dD_P_H = 1/(dP_D - dP_T*dH_D/dH_T)
+dD_P_H = simplify(dD_P_H)
+dD_P_H = expand(dD_P_H)
+
+dD_H_P = 1/(dH_D - dH_T*dP_D/dP_T)
+dD_H_P = simplify(dD_H_P)
+dD_H_P = expand(dD_H_P)
+
+
+%% PVT second order derivatives
+d2T_D2_P = -(d2P_D2*dP_T - dP_D*d2P_TD)/dP_T^2 ...
+           +(d2P_TD*dP_T - dP_D*d2P_T2)*dP_D/dP_T^3
+d2T_D2_P = simplify(d2T_D2_P)
+d2T_D2_P = expand(d2T_D2_P)
+
+d2D_T2_P = -(d2P_T2*dP_D - dP_T*d2P_TD)/dP_D^2 ...
+           +(d2P_TD*dP_D - dP_T*d2P_D2)*dP_T/dP_D^3
+d2D_T2_P = simplify(d2D_T2_P)
+d2D_T2_P = expand(d2D_T2_P)
+
+
+%% derivative for fundamental derivative of gas dynamics
+d2P_D2_S = d2P_D2 -(dP_T*d2S_D2   + 2*d2P_TD*dS_D)/dS_T ...
+                  +(d2P_T2*dS_D^2 + 2*dP_T*dS_D*d2S_TD)/dS_T^2 ...
+                  -(dP_T  *dS_D^2*d2S_T2)/dS_T^3
+d2P_D2_S = simplify(d2P_D2_S)
+d2P_D2_S = expand(d2P_D2_S)
diff --git a/HelmholtzMedia/Examples/Validation/Derivatives_SinglePhase.mo b/HelmholtzMedia/Examples/Validation/Derivatives_SinglePhase.mo
@@ -50,7 +50,8 @@ model Derivatives_SinglePhase
   Medium.Types.Der2EntropyByTemperatureDensity d2sTd_analytical;
   Medium.Types.Der2EntropyByTemperatureDensity d2sTd_numerical;
 // Energy wrt. dT
-  Medium.Types.DerEnergyByDensity dudT_analytical;
+  Medium.Types.DerEnergyByDensity dudT_analytical1;
+  Medium.Types.DerEnergyByDensity dudT_analytical2;
   Medium.Types.DerEnergyByDensity dudT_numerical;
   Medium.Types.DerEnergyByTemperature duTd_analytical;
   Medium.Types.DerEnergyByTemperature duTd_numerical;
@@ -66,7 +67,8 @@ model Derivatives_SinglePhase
   Medium.Types.DerEnthalpyByDensity dhdT_analytical1;
   Medium.Types.DerEnthalpyByDensity dhdT_analytical2;
   Medium.Types.DerEnthalpyByDensity dhdT_numerical;
-  Medium.Types.DerEnthalpyByTemperature dhTd_analytical;
+  Medium.Types.DerEnthalpyByTemperature dhTd_analytical1;
+  Medium.Types.DerEnthalpyByTemperature dhTd_analytical2;
   Medium.Types.DerEnthalpyByTemperature dhTd_numerical;
 // Enthalpy wrt. dT 2nd order
   Medium.Types.Der2EnthalpyByDensity2 d2hd2T_analytical1;
@@ -297,10 +299,12 @@ equation
   Modelica.Utilities.Streams.print("====|====|====|====|====|====|====|====|====|====|====|====|====|====|====|====|"); // 80 characters
   Modelica.Utilities.Streams.print("internal Energy wrt. dT");
   // check (du/dd)@T=const
-  dudT_analytical = Medium.EoS.dudT(f);
+  dudT_analytical1 = Medium.EoS.dudT(f);
+  dudT_analytical2 = -state.T/state.d^2*Medium.EoS.dpTd(f) + state.p/state.d^2;
   dudT_numerical = (dplus_Tconst.u-dminus_Tconst.u)/(dplus_Tconst.d-dminus_Tconst.d);
-  Modelica.Utilities.Streams.print("  (du/dd)@T=const analytical= " + String(dudT_analytical));
-  Modelica.Utilities.Streams.print("  (du/dd)@T=const  numerical= " + String(dudT_numerical));
+  Modelica.Utilities.Streams.print("  (du/dd)@T=const analytical1= " + String(dudT_analytical1));
+  Modelica.Utilities.Streams.print("  (du/dd)@T=const analytical2= " + String(dudT_analytical2));
+  Modelica.Utilities.Streams.print("  (du/dd)@T=const   numerical= " + String(dudT_numerical));
   // check (du/dT)@d=const
   duTd_analytical = Medium.EoS.duTd(f);
   duTd_numerical = (Tplus_dconst.u-Tminus_dconst.u)/(Tplus_dconst.T-Tminus_dconst.T);
@@ -334,10 +338,12 @@ equation
   Modelica.Utilities.Streams.print("  (dh/dd)@T=const analytical2= " + String(dhdT_analytical2));
   Modelica.Utilities.Streams.print("  (dh/dd)@T=const   numerical= " + String(dhdT_numerical));
   // check (dh/dT)@d=const
-  dhTd_analytical = Medium.EoS.dhTd(f);
+  dhTd_analytical1 = Medium.EoS.dhTd(f);
+  dhTd_analytical2 = Medium.EoS.duTd(f) + 1/state.d*Medium.EoS.dpTd(f);
   dhTd_numerical = (Tplus_dconst.h-Tminus_dconst.h)/(Tplus_dconst.T-Tminus_dconst.T);
-  Modelica.Utilities.Streams.print("  (dh/dT)@d=const analytical= " + String(dhTd_analytical));
-  Modelica.Utilities.Streams.print("  (dh/dT)@d=const  numerical= " + String(dhTd_numerical));
+  Modelica.Utilities.Streams.print("  (dh/dT)@d=const analytical1= " + String(dhTd_analytical1));
+  Modelica.Utilities.Streams.print("  (dh/dT)@d=const analytical2= " + String(dhTd_analytical2));
+  Modelica.Utilities.Streams.print("  (dh/dT)@d=const   numerical= " + String(dhTd_numerical));
   // check (d2h/dd2)@T=const
   d2hd2T_analytical1 = Medium.EoS.d2hd2T(f);
   d2hd2T_analytical2 = -state.T/state.d^2*Medium.EoS.d2pTd(f) + 1/state.d*Medium.EoS.d2pd2T(f) - 1/state.d^2*Medium.EoS.dpdT(f) + 2*state.T/state.d^3*Medium.EoS.dpTd(f);