Pressure-volume isotherms 
for a Bragg-Williams lattice gas
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :
T. L. Hill, An Introduction to Statistical Thermodynamics, Dover, New York, 1986.

p[u_] := -k T Log[1 - 1/u] + c w/(2 u^2)

plt1 = 
      Plot[-p[u]/(c w) /. 
          k T ->  (c w)/(-3), {u, 1, 10}, 
        PlotStyle ->  RGBColor[0, 0, 1], 
        DisplayFunction -> Identity]

plt2 = 
      Plot[-p[u]/(c w) /. 
          k T ->  (c w)/(-4), {u, 1, 10}, 
        PlotStyle ->  RGBColor[0, 1, 0], 
        DisplayFunction -> Identity]

plt3 = 
      Plot[-p[u]/(c w) /. 
          k T ->  (c w)/(-5), {u, 1, 10}, 
        PlotStyle ->  RGBColor[1, 0, 0], 
        DisplayFunction -> Identity]

Plotting Isotherms. We get result similar to the van der Waals equation. The Bragg-Williams approximation pedicts a first-order phase transition.

Show[plt1, plt2, plt3, 
      DisplayFunction -> $DisplayFunction]

We use the symmetry condition:  (Subscript[θ, l]) + (Subscript[θ, g]) = 1 to draw the horizontal stable equilibrium path. We can also use the thermodynamic equal-area theorem shown below.

sol = 
      FindRoot[{Log[(ug - 1)/(ul - 1)] ==  
                c w/(k T) (1/ug - 
                      1/ul) /. (c w)/(k T) ->  -5, 
1/ug + 1/ul ==  1}, {ug, 0.1}, {ul, 0.5}] // 
        Chop

-p[ug]/(c w) /. k T ->  (c w)/(-5) /. 
      sol

-p[ul]/(c w) /. k T ->  (c w)/(-5) /. 
      sol

sol2 = 
      FindRoot[(-p[u]/(c w) /. 
                  k T ->  (c w)/(-5)) ==  -p[
                    ug]/(c w) /. 
            k T ->  (c w)/(-5) /. sol, {u, 
          2}]

plt5 = 
      Show[plt3, 
        Epilog -> {Line[{{ug /. 
                    sol, -p[ug]/(c w) /. 
                      k T ->  (c w)/(-5) /. sol}, {ul /. 
                    sol, -p[ul]/(c w) /. 
                      k T ->  (c w)/(-5) /. 
                    sol}}], {PointSize[0.025], 
              Point[{ug /. 
                    sol, -p[ug]/(c w) /. 
                      k T ->  (c w)/(-5) /. sol}], 
              Point[{ul /. 
                    sol, -p[ul]/(c w) /. 
                      k T ->  (c w)/(-5) /. 
                    sol}], 
Point[{u /. 
                    sol2, -p[u]/(c w) /. 
                      k T ->  (c w)/(-5) /. sol2}]}}, 
        DisplayFunction -> $DisplayFunction, 
        PlotRange -> {{0, 10}, {-0.01, 0.1}}, 
        DisplayFunction -> $DisplayFunction, AspectRatio ->  1.1, 
        Frame ->  True, AxesOrigin -> {0, 0}]

<< Graphics`FilledPlot`

plt7 = 
      FilledPlot[{-p[ug]/(c w) /. 
              k T ->  (c w)/(-5) /. 
            sol, (-p[u]/(c w) /. 
              k T ->  (c w)/(-5))}, {u, ug /. sol, 
          u /. sol2}, DisplayFunction -> Identity]

plt8 = 
      FilledPlot[{-p[ug]/(c w) /. 
              k T ->  (c w)/(-5) /. 
            sol, (-p[u]/(c w) /. 
              k T ->  (c w)/(-5))}, {u, u /. sol2, 
          ul /. sol}, DisplayFunction -> Identity]

As expected, both shaded areas are equal.

Show[plt5, plt8, plt7, 
      PlotRange ->  {{0, 10}, {-0.01, 0.06}}, 
      AspectRatio ->  1.1, Frame ->  True, 
      DisplayFunction -> $DisplayFunction]

NIntegrate [-(-p[u]/(c w) /. 
              k T ->  (c w)/(-5)) + (-p[
                    ug]/(c w) /. 
              k T ->  (c w)/(-5) /. sol), {u, 
        ug /. sol, u /. sol2}]

NIntegrate [(-p[u]/(c w) /. 
                k T ->  (c w)/(-5)) + p[ug]/(c w) /. 
          k T ->  (c w)/(-5) /. sol, {u, u /. sol2, 
        ul /. sol}]