6
$\begingroup$

A Maxwell construction is just an integral of area and then solving for a root to find vapor pressure.

This must be exactly the same value as where the Gibbs energy plot crosses at a point. I found that to be p = 0.81.

The problem is that my areadifferential function doesn't work properly; it doesn't return p = 0.81;

What mistake I have made in my coding?

t = 0.95; 

p[v_] := (8*t)/(3*v - 1) - 3/v^2;

Plot[p[v], {v, 0.5, 3}, 
  PlotRange -> {{0, 3}, {0.6, 1}}, 
  AxesLabel -> {V/Vc, P/Pc}]

g[v_] := (-t)*Log[3*v - 1] + 0.95/(3*v - 1) - 9/(4*v); 

ParametricPlot[{p[v], g[v]}, {v, .65, 2.25},
  AxesLabel -> {P/Pc, G/NKT},
  PlotLabel -> "Gibbs Free Energy Vs. P/Pc"]

That returns pressure p = P/Pc = 0.81 as a plot which is correct

enter image description here

This is the error Maxwell construction "area differential" part

pint[v_] := (8/3)*t*Log[3*v - 1] + 3/v;

areadifferential[p0_, v1guess_, v2guess_] := 
  (v1 = FindRoot[p[v] == p0, {v, v1guess}][[1,2]]; 
   v2 = FindRoot[p[v] == p0, {v, v2guess}][[1,2]];
pint[v2] - pint[v1] - p0*(v2 - v1))

FindRoot[areadifferential[p0, 0.7, 2] == 0, {p0, 0.8, 0.82}]
$\endgroup$
9
$\begingroup$

Good job on describing your problem.

Copying your definitions from your question

p[v_] := (8 t)/(3 v - 1) - 3/v^2
g[v_] := t/(3 v - 1) - t Log[-1 + 3 v] - 9/(4 v)
pint[v_] := 8/3 t Log[-1 + 3 v] + 3/v
t = 0.95

I found the value 0.81188 to be the value that minimized your areadifferential function (shown at the bottom).

Look at the plot of P/Pc vs V/Vc and note the three points corresponding to three different V/Vc values resulting in p[v] = 0.8118.

Module[
 {
  v1 = v /. FindRoot[p[v] == 0.8118, {v, 0.5}],
  p1,
  v2 = v /. FindRoot[p[v] == 0.8118, {v, 1}],
  p2,
  v3 = v /. FindRoot[p[v] == 0.8118, {v, 2}],
  p3
  },
 p1 = p[v1];
 p2 = p[v2];
 p3 = p[v3];
 Show[
  Plot[p[v], {v, 0.65, 2.25}, PlotRange -> {{0, 3}, {0.6, 1}}, 
   AxesLabel -> {"V/Vc", "P/Pc"}],
  Graphics[
   {
    PointSize[0.03],
    Red,
    Point[{v1, p1}],
    Green,
    Point[{v2, p2}],
    Black,
    Point[{v3, p3}]
    }
   ]
  ]
 ]

Mathematica graphics

The problem that you are experiencing is because the value of 0.7 for v1guess is just to the right of the valley so it converges on the wrong point (i.e., the green point).

Change v1guess to 0.5 and you will be fine, it will converge on the red point. Using v2guess of 2.0 works fine causing convergence on the black point.

I slightly modified the areadifferential function to use Module rather than parenthesis. One needs to constrain the input arguments to be numeric so FindRoot doesn't complain (FindFoot first tries to work with symbolic arguments which doesn't work with your problem).

areadifferential[
p0_?NumericQ,
v1guess_?NumericQ,
v2guess_?NumericQ] := Module[
  {
   v1 = FindRoot[p[v] == p0, {v, v1guess}][[1, 2]],
   v2 = FindRoot[p[v] == p0, {v, v2guess}][[1, 2]]
   },
  pint[v2] - pint[v1] - p0*(v2 - v1)
  ]

and then

FindRoot[areadifferential[p0, 0.5, 2] == 0, {p0, 0.8, 0.82}]
(* {p0 -> 0.811879} *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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