3
$\begingroup$

I'm working on a psychrometric chart on Mathematica. The code is as follows:

(*psat is saturation pressure calculated by Hyland Wexler correlation*)
psat[T_?NumericQ] := 
 Exp[(c8/(T + 273.15)) + 
   c9 + (c10*(T + 
       273.15)) + (c11*((T + 273.15)^2)) + (c12*((T + 
         273.15)^3)) + (c13*Log[T + 273.15])]
c8 = -5800.2206;
c9 = 1.3914993;
c10 = -0.04860239;
c11 = 0.41764768*10^-4;
c12 = -0.14452093*10^-7;
c13 = 6.5459673;
patm = 101325;(*Pa*)(*atm pressure*)
rn = 0.622;(*kg/kg*);(*molar ratio*)
cp = 1.004;(*KJ/kg K*);(*air specific heat*)
hvref = 2501;(*KJ/kg*);(*vapour enthalpy*)
cpda = 1.006;(*KJ/kg K*);(*specific heat dry air*)
cpwv = 1.84;(*KJ/kg K*);(*specific heat water vapour*)

(*Y is humidity ratio as function of relative humidity and temperature*)
Y[phi_?NumericQ, T_?NumericQ] := (rn*phi*psat[T])/(patm - psat[T])

phiPlot = Table[Plot[Y[phi, T], {T, 0, 50}, PlotStyle -> {Thin, Blue}], {phi,0.1, 1, 0.1}];

phiSatPlot = Plot[Y[1, T], {T, 0, 50}, PlotStyle -> Green, PlotRange -> {{0, 50}, {0, 0.050}}];

Tref = 0;

(*YB is humidity ratio as function of enthalpy and temperature*)
YB[h_?NumericQ,T_?NumericQ] := (h - (cpda*(T + 273.15 - (Tref + 273.15))))/((cpwv*(T + 273.15 - (Tref + 273.15))) + hvref)

h1Plot = Table[Plot[YB[h, T], {T, 0, 55}, PlotStyle -> {Thin, Red},PlotRange -> {{0, 50}, {0, 0.050}}], {h, 0, 200, 5}];

Show[h1Plot, phiSatPlot, phiPlot, Frame -> True, 
 FrameTicks -> {{None, Range[0, 0.050, 0.002]}, {Range[0, 50, 10], 
    None}}]

So the following graphic shows humidity ratio (Y axis), dry bulb temperature (x axis), relative humidity (blue lines) and enthaply

psychrometric chart

The issue is that I can't erase the enthalpy (red) lines above the saturation line (green curve). Does anybody knows how I could do that?

$\endgroup$
3
  • 1
    $\begingroup$ Welcome to the Mathematica Stack Exchange. Could you please include definition for psat? $\endgroup$
    – Syed
    Commented Nov 13, 2022 at 5:12
  • $\begingroup$ phiSatPlot = Plot[Y[1, T], ... What is Y? Please add all this info to your post and delete your comments. Make sure you copy from this web page to a new notebook on a fresh kernel. If you can't see the plot, we can't see it either. $\endgroup$
    – Syed
    Commented Nov 13, 2022 at 10:50
  • $\begingroup$ Thank you @Syed! I did the alterations and tried the code in a new notebook to check the plot. $\endgroup$ Commented Nov 13, 2022 at 11:29

3 Answers 3

4
$\begingroup$
Clear["Global`*"]

c8 = -5800.2206;
c9 = 1.3914993;
c10 = -0.04860239;
c11 = 0.41764768*10^-4;
c12 = -0.14452093*10^-7;
c13 = 6.5459673;
patm = 101325;
Tref = 0;

(*Pa*)(*atm pressure*)rn = 0.622;(*kg/kg*);(*molar ratio*)cp \
    = 1.004;(*KJ/kg K*);(*air specific heat*)hvref = \
    2501;(*KJ/kg*);(*vapour enthalpy*)cpda = 1.006;(*KJ/kg K*);(*specific \
    heat dry air*)cpwv = 1.84;(*KJ/kg K*);(*specific heat water \
    vapour*)(*Y is humidity ratio as function of relative humidity and \
    temperature*)

    Y[phi_?NumericQ, T_?NumericQ] := (rn*phi*psat[T])/(patm - psat[T])
    
    (*psat is saturation pressure calculated by Hyland Wexler correlation*)


psat[T_?NumericQ] := 
 Exp[(c8/(T + 273.15)) + 
   c9 + (c10*(T + 
       273.15)) + (c11*((T + 273.15)^2)) + (c12*((T + 
         273.15)^3)) + (c13*Log[T + 273.15])]

(*YB is humidity ratio as function of enthalpy and temperature*)
YB[h_?NumericQ, 
  T_?NumericQ] := (h - (cpda*(T + 
        273.15 - (Tref + 273.15))))/((cpwv*(T + 
        273.15 - (Tref + 273.15))) + hvref)

h1Plot = ListLinePlot[
   Table[YB[h, T], {h, 0, 200, 5}, {T, 0, 55}]
   , PlotStyle -> {{Thin, Red}}
   ];

phiPlot = 
  Plot[Evaluate@Table[Y[phi, T], {phi, 0.1, 0.9, 0.1}], {T, 0, 50}
   , PlotStyle -> Blue
   , PlotRange -> {0, 0.05}
   , PlotRangePadding -> None
   , Frame -> True
   , FrameTicks -> {
     {None, Range[0, 0.050, 0.002]}
     , {Range[0, 50, 10], None}
     }
   , ImageSize -> 500
   ];

phiSatPlot = Plot[Y[1, T]
   , {T, 0, 50}
   , PlotStyle -> {Thick, Darker@Green}
   , PlotRange -> {{0, 50}, {0, 0.050}}
   , FillingStyle -> Directive[Red, HatchFilling[-10 Degree, 0.2, 7]]
   , Filling -> {1 -> {Top, Directive[Opacity[1], White]}}
   
   ];

Show[phiPlot, h1Plot, phiSatPlot]

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ Thank you @Syed! It helped me a lot! $\endgroup$ Commented Nov 13, 2022 at 19:49
4
$\begingroup$

Add a white domain erasePlot to the plot.

erasePlot = 
  Plot[Y[1, T], {T, 0, 50}, PlotStyle -> Green, 
   PlotRange -> {{0, 50}, {0, 0.050}}, Filling -> Top, 
   FillingStyle -> White];
Show[h1Plot, phiSatPlot, phiPlot, erasePlot, Frame -> True, 
 FrameTicks -> {{None, Range[0, 0.050, 0.002]}, {Range[0, 50, 10], 
    None}}]

enter image description here

(*psat is saturation pressure calculated by Hyland Wexler correlation*)
psat[T_?NumericQ] := 
 Exp[(c8/(T + 273.15)) + 
   c9 + (c10*(T + 
       273.15)) + (c11*((T + 273.15)^2)) + (c12*((T + 
         273.15)^3)) + (c13*Log[T + 273.15])]
c8 = -5800.2206;
c9 = 1.3914993;
c10 = -0.04860239;
c11 = 0.41764768*10^-4;
c12 = -0.14452093*10^-7;
c13 = 6.5459673;
patm = 101325;(*Pa*)(*atm pressure*)
rn = 0.622;(*kg/kg*);(*molar ratio*)
cp = 1.004;(*KJ/kg K*);(*air specific heat*)
hvref = 2501;(*KJ/kg*);(*vapour enthalpy*)
cpda = 1.006;(*KJ/kg K*);(*specific heat dry air*)
cpwv = 1.84;(*KJ/kg K*);(*specific heat water vapour*)

(*Y is humidity ratio as function of relative humidity and \
temperature*)
Y[phi_?NumericQ, T_?NumericQ] := (rn*phi*psat[T])/(patm - psat[T])

phiPlot = 
  Table[Plot[Y[phi, T], {T, 0, 50}, PlotStyle -> {Thin, Blue}], {phi, 
    0.1, 1, 0.1}];

phiSatPlot = 
  Plot[Y[1, T], {T, 0, 50}, PlotStyle -> Green, 
   PlotRange -> {{0, 50}, {0, 0.050}}];

Tref = 0;

(*YB is humidity ratio as function of enthalpy and temperature*)
YB[h_?NumericQ, 
  T_?NumericQ] := (h - (cpda*(T + 
        273.15 - (Tref + 273.15))))/((cpwv*(T + 
        273.15 - (Tref + 273.15))) + hvref)

h1Plot = 
  Table[Plot[YB[h, T], {T, 0, 55}, PlotStyle -> {Thin, Red}, 
    PlotRange -> {{0, 50}, {0, 0.050}}], {h, 0, 200, 5}];

erasePlot = 
  Plot[Y[1, T], {T, 0, 50}, PlotStyle -> Green, 
   PlotRange -> {{0, 50}, {0, 0.050}}, Filling -> Top, 
   FillingStyle -> White];
Show[h1Plot, phiSatPlot, phiPlot, erasePlot, Frame -> True, 
 FrameTicks -> {{None, Range[0, 0.050, 0.002]}, {Range[0, 50, 10], 
    None}}]
$\endgroup$
3
  • $\begingroup$ Thanks @cvgmt! I don´t think I can remove h1Plot entirely since it shows enthalpy lines that are important to the chart. $\endgroup$ Commented Nov 13, 2022 at 11:57
  • $\begingroup$ @JulioAraujoDosSantos updated. $\endgroup$
    – cvgmt
    Commented Nov 13, 2022 at 12:07
  • $\begingroup$ thank you very much my friend! It solves the problem. $\endgroup$ Commented Nov 13, 2022 at 12:56
4
$\begingroup$

You can use RegionFunction option with your boundary function Y[1,T]

Note that you don't need to plot each contour separately. The method used below for newh1Plot will save memory.

With definitions as in OP then

newh1Plot =
 Plot[
  Table[YB[h, T], {h, 0, 200, 5}]
  , {T, 0, 55}
  , PlotStyle -> {Thin, Red}
  , PlotRange -> {{0, 50}, {0, 0.050}}
  , RegionFunction -> Function[{x, y}, y <= Y[1, x]]
  ]

Mathematica graphics

and

Show[newh1Plot, phiSatPlot, phiPlot
 , Frame -> True
 , FrameTicks -> {{None, Range[0, 0.050, 0.002]}, {Automatic, None}}
 ]

Mathematica graphics

Hope this helps.

$\endgroup$
1
  • 1
    $\begingroup$ thank you very much, it helps a lot! And also thank you for the memory optmization tip. $\endgroup$ Commented Nov 13, 2022 at 12:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.