3.016 Problem Set #3, 2010

I3-1

Before starting, let's just make some matrices to use in this and following problems

Matrix 1 (matOne) below is from the first set of equations:

matOne = {
{1, -3/2,-1},
{1/3,-1/3,1},
{-2 , -1, 2}};
MatrixForm[matOne]

Now Matrix 1 (matOne) is listed above in matrix form.

Matrix 2 (matTwo) below is from the second set of equations:

matTwo = {
{1, -3/2,-1},
{1/3,-1/3,1},
{7 , -9, 5}};
MatrixForm[matTwo]

Now Matrix 2 (matTwo) is listed above in matrix form.

Solve the two sets of equations in the traditional way. In each case, we have three equations and three unknowns.

Now create the unknown vector, and the two right hand side vectors:

unknown = {p,q,r};
rhs1 = 33{3,1,3};
rhs1a=33{2,1,2};
rhs2 = {0,0,0};

Below are the three equations for the first problem:

equations1 = Table[rhs1[[i]] == matOne[[i]].unknown,{i,1,3}]

Below are the three equations for the problem 1a (matching the case to be studied in I3-2):

equations1a = Table[rhs1a[[i]] == matOne[[i]].unknown,{i,1,3}]

Below are the three equations for the second problem:

equations2 = Table[rhs2[[i]] == matTwo[[i]].unknown,{i,1,3}]

Solving for the first problem:

sol1 = Solve[equations1,unknown]
unknown/.sol1

Solving for the problem 1a (matching the case to be studied in I3-2):

sol1a = Solve[equations1a,unknown]
unknown/.sol1a

Solving for the second problem:

sol2 = Solve[equations2,unknown]
unknown/.sol2

Note that there are infinitely many solutions for the second problem.  Geometrically all the solutions form a line and not just one point.

I3-2

Generate the inverse matrix:

invMat1 = Inverse[matOne];
invMat1//MatrixForm

Demonstrate that we get the same solution above (problem 1a):

invMat1.rhs1a

Graphical representation of the solution:

Obtain parametric expressions of the three equations (planes) in terms of p and q:

threeplanes = Table[ r/.Solve[equations1a[[i]],r],{i,1,3}]

Plot the three planes with plotting ranges include and near the solution point (experimented with different colors, p/q ranges and turned mesh on/off to get a better graphical output):

plotthreeplanes = Plot3D[threeplanes,{p,0,10},{q,-55,-45},PlotStyle-> {Red,Blue,Yellow},Mesh-> False]

Showing the solution with a green sphere, together with the three planes (experimented with different colors and sizes of the sphere, and rotated the plot to see it from different angles):

Show[plotthreeplanes,Graphics3D[{Green,Sphere[invMat1.rhs1a,.25]}]]

I3-3

The inverse doesn't exist because the determinant is zero

First trying to invert the matrix:

invMat2 = Inverse[matTwo]

... and found that the matrix is singular. Now checking the determinant:

Det[matTwo]

Indeed the determinant is zero.
Therefore we cannot multiply it by the right hand side vector.

However, the problem still has a graphical representation...

Obtain parametric expressions of the three equations (planes) in terms of p and q:

threeplanes2 = Table[ r/.Solve[equations2[[i]],r],{i,1,3}]

Plot the three planes with plotting ranges include and near the solution line (experimented with different colors, p/q ranges and turned mesh on/off to get a better graphical output):

plotthreeplanes3 = Plot3D[threeplanes2,{p,-80,0},{q,-35,-45},PlotStyle-> {Red,Blue,Yellow},Mesh-> False]

Obtain all the solutions (line):

line=unknown/.sol2[[1]]

Find two points on the line for plotting:

points = {line/.r->20,line/.r-> -30}

Showing the solution line with a green cylinder, together with the three planes (experimented with different colors and sizes of the sphere, and rotated the plot to see it from different angles):

Show[plotthreeplanes3,Graphics3D[{Green,Cylinder[points,2]}]]

I3-4

Exploring to understand better the "equation of state"

Read from http://en.wikipedia.org/wiki/Equation_of_state regarding "equation of state", "van der Waals equation of state", and the physical meaning of "b". Found out that "b" represents "the repulsion parameter or the effective molecular volume.

4-1: conditions on a and b

Write down expressions for the equation of state in terms of P and T:

Pvdw = R T/(V-b) - a/V^2

Tvdw = (P+ a/V^2)(V-b)/R

Now we will find the conditions in a couple steps...

First specify that the equation of state for T must be positive (note that there are implicit assumptions that R>0 and V>0):

cond1 =FullSimplify[Reduce[Tvdw >0,{a,b},Reals],Assumptions -> P>0 && R> 0 && V> 0]

We can do the same thing specifying that P must be postive:

cond2 =FullSimplify[Reduce[Pvdw >0,{a,b},Reals]/.T-> Tvdw,Assumptions -> T>0 && R> 0 && V> 0]

Note that, in the above we replaced T using the the expression of Tvdw, so the conditions would be expressed in terms of P, V, R, a, b, so that cond1 and cond2 can be directly compared.

Trying to simply jointly cond1 and cond2:

FullSimplify[Reduce[cond1 && cond2 ,{a,b},Reals],Assumptions -> R> 0 && V> 0 && P > 0]

This is the same as cond1. This confirms that cond1 and cond2 are equivalent!

Thus, we should have a < -P (V^2) and b > V (this is unlikely because b is supposed to represent the effective volume of one molecule...) or b < V  and a > -P (V^2)

Then the physically relevant condition should be b < V and a > -P (V^2).

4-2: water isotherms

First, find out typical P-V behavior online to understand ideal gases and real gases. As an example, http://www.brighthub.com/engineering/mechanical/articles/35392.aspx shows good example plots. Here http://www.daviddarling.info/encyclopedia/V/van_der_Waals_equation_of_state.html gives anthor schematic plot.

Generating the specific equation of state, at three given temperatures:

PvdwWater = Pvdw/.{a-> 558 , b-> 3.05 10^(-5), R-> 8.03}

threeIsotherms = {PvdwWater/.T-> 700, PvdwWater/.T-> 647,PvdwWater/.T-> 600}

Plotting the results, and exploring different ranges of V:

Plot[threeIsotherms,{V,0,0.8}, AxesLabel->{V,P}]

The behavior does not look like the correct behavior for water vapor!!!

Trying to find out if the parameters given are really the typical/reasonable values. In fact, looking at http://hyperphysics.phy-astr.gsu.edu/hbase/kinetic/waal.html, it looks like  a = 0.558 pascal . (m^6)  is a more reasonable parameter with the correct order of maganitude and correct units. The other parameters look alright.

Generating the specific equation of state again (with a = 0.558), at three given temperatures:

PvdwWater = Pvdw/.{a-> .558 , b-> 3.05 10^(-5), R-> 8.03}

threeIsotherms = {PvdwWater/.T-> 700, PvdwWater/.T-> 647,PvdwWater/.T-> 600}

Plotting the new results (with a = 0.558), and exploring different ranges of V:

One of the difficulties of working with numbers like this is that one has to "hunt" for the right region to plot (this is why it is useful to non-dimensionalize. Here is a result of hunting around...

Plot[threeIsotherms,{V,0.00004,.0005}, AxesLabel->{V,P}]

Now the behavior does look like physically meaningful behavior for water!

Here are more isotherms just to visualize...

manyIsotherms =Table[ PvdwWater/.T-> isoT,{isoT,100,1000,100}];

Plot[manyIsotherms,{V,0.00004,.0005}, AxesLabel->{V,P}]

It would be much easier to calculate the critical point. The conditions are that the isotherm's slope and inflection point appear at the same conditions.

Pvdw = R T/(V-b) - a/V^2

Tvdw = (P+ a/V^2)(V-b)/R

cond1 = 0 == D[Pvdw,V]==0
cond2  = 0 == D[Pvdw,{V,2}]==0

This will give the temperature and volume at the critical point in terms of a and b

critPointConds = Solve[{cond1,cond2},{T,V}][[1]]

The pressure at the critical point can be determined by back-substituting these volumes and pressures.

Pvdw/.critPointConds

Thus, for water,

critPointConds/.{a-> .558 , b-> 3.05 10^(-5), R-> 8.03}

The predicted molar volume at the critical point is 0.0000915 cubic meters.

And the critical point pressure is (in pascals)

Pvdw/.critPointConds/.{a-> .558 , b-> 3.05 10^(-5), R-> 8.03}

4-3: Non-Dimensionalizing

The van der Waals equation of state:

vdwEq = R T== (P + a/V^2)(V-b)

Critical conditions:

{Tcrit,Vcrit} = {T,V}/.critPointConds

Pcrit  = Pvdw/.{T-> Tcrit, V-> Vcrit}

Now using non-dimensionalized parameters Θ, Ω and Π to express the van der Waals equation of state:

vdwEqScaled = vdwEq/.{T-> Θ Tcrit , V-> Ω Vcrit,P -> Π Pcrit}

Solve the non-dimensionalized equation vdwEqScaled for Π:

ΠvdwRule= Solve[vdwEqScaled,Π][[1]]

Πvdw =  Π/.ΠvdwRule

Solve the non-dimensionalized equation vdwEqScaled for Θ:

ΘvdwRule = Solve[vdwEqScaled,Θ][[1]]

Θvdw=Θ/.ΘvdwRule

4-4: Visualizing for the case of a generic gas

This is now easy because if we plot in these normalized coordinates, the critical transition point appears at Π=Ω=Θ=1...

Generate non-dimensionalized isothermal functions near critical temperature:

isotherms = Table[Πvdw/.Θ->temp,{temp,.9,1.1, .025}];

Plot the generated functions:

P1 = Plot[isotherms,{Ω,.5,3},PlotRange-> {0.1,2},ImageSize-> Large, BaseStyle-> Large, AxesLabel-> {Ω,Π}]

Plot the functions, overlapping with the critical curve at Θ = 1:

Show[P1,Plot[Πvdw/.Θ->1,{Ω,0.1,3},PlotStyle-> {Thickness[0.01],Cyan},PlotRange-> All]]

4-5: bulk modulus

Bulk modulus β = -V(((∂P)/(∂V))) with T= contant, so we know β has the units of pressure (pascal).

We therefore can define non-dimensional bulk modulus as B = β/Pcrit, 
and immediately we can obtain B = -Ω ∂Π/∂Ω.

BulkModDimensionless =Simplify[ -Ω D[Πvdw,Ω]]

To see how this behaves near the critical point, the most straightforward thing to do is plot it.

isothermsForBulkModulus = Table[BulkModDimensionless/.Θ-> temp1,{temp1,0.9,1.1,.05}]

Plot[isothermsForBulkModulus,{Ω,0.9,1.1},AxesLabel-> {"Ω", "Dimensionless β"}, PlotLabel-> "Isotherms for Bulk Modulus near Critical Point", ImageSize-> Large, BaseStyle-> Large]

The critical point is where the bulk modulus becomes negative.  In fact, a negative bulk modulus creats an unstable material.

4-6: coefficient of linear expansion

α = (1/V)(((∂V)/(∂T))) with P= contant, so the coefficient of thermal expansion α has the units of inverse temperature.

We therefore can define non-dimensional coefficient of thermal expansion as A = α Tcrit, 
and immediately we can obtain A = (1/Ω) ∂Ω/∂Θ.

It will be simpler to compute the inverse of A. To do this, first find Θ(Π,Ω)

Θvdw = Θ/.Solve[Π == Πvdw,Θ][[1]]

Now compute 1/A:

inverseCoTEDimensionless =Simplify[ Ω D[Θvdw,Ω]]

Finally obtain A:

CoTEDimensionless =1/inverseCoTEDimensionless

To see how this behaves near the critical point, the most straightforward thing to do is plot it.

isoPressuresForCoTE = Table[CoTEDimensionless/.Π-> p,{p,0.9,1.1,.05}]

Plot[isoPressuresForCoTE,{Ω,0.9,1.1},AxesLabel-> {"Ω", "Dimensionless α"}, PlotLabel-> Style["Iso-Pressures for Coefficient of Thermal Expansion near Critical Point",FontFamily->"ArialNarrow",FontWeight->Bold,20], ImageSize-> Large, BaseStyle-> Large,PlotStyle-> {Red,Orange,Green,Cyan,Blue}]

The critical point is where the coefficient of thermal expansion goes from being negative (unstable) and then has a singularity as it goes to infinity at the critical point.

4-7: Cp - Cv

The dimensionless form will be ΘΩA^2B:

CPminusCV =Simplify[Θ Ω CoTEDimensionless^2 BulkModDimensionless]

Use the van der Waals equation to make CpminusCv a function only of Π and Ω...

CPminusCV =Simplify[ CPminusCV/.Θ-> Θvdw]

To see how this behaves near the critical point, the most straightforward thing to do is plot it.

isoPressuresForCpmCv = Table[CPminusCV/.Π-> p,{p,0.9,1.1,.05}]

Plot[isoPressuresForCpmCv,{Ω,0.9,1.1},AxesLabel-> {"Ω", "Dimensionless Cp-Cv"}, PlotLabel-> "Iso-Pressures for Cp-Cv near Critical Point", ImageSize-> Large, BaseStyle-> Large,PlotStyle-> {Red,Orange,Green,Cyan,Blue}]

The critical point for this case has the same type of behavior as that for the coefficient of thermal expansion. The system goes from unstable to stable as the pressure is increased past the critical point.