Bernard Thermal Instability: Marginal Mode Analysis

APPH 4200 Physics of Fluids
Columbia University
November, 2007

Introduction

This notebook solves for the marginal instability threshold for Bernard thermal instability. See Chapter 12, Section 3 of Kundu and Cohen's textbook.

Four steps are required. First, from the linearized equations for the fluid dynamics, a cubic equation is solved for the vertical variation of the vertical velocity. Second, these three solutions are used to satisfy three boundary conditions at the constant-temperature walls. This leads to a condition between Ra (the Rayleigh number) and (k^2)for marginal instability. Finally, we plot this condition and examine the modes.

Vertical Velocity Equation

The equation for vertical velocity takes the form

eqVert = ((((q^2) + (k^2)))^3) ==  (k^2)Ra

when the vertical velocity has the form, (Subscript[U, z]) ~ Cos[q z].

qSol = Solve[(eqVert /. (q^2) -> q2 /. (k^2) -> k2), q2]

Cos[z] is Even

First, we'll look for even modes...

Plot[Cos[z],{z,-(π/2),(π/2)}]

Plot[Evaluate[({Re[#1],Im[#1]}&)[Cos[z (1+I (Sqrt[3]))]]],{z,-(π/2),(π/2)},PlotStyle->{Blue,Red}]

An similiar analysis can be performed for the odd modes, but these modes have a higher critical Ra number.

Satisfy Boundary Conditions

uz[z_] = a (Cos[Sqrt[q2] z]/.qSol[[1]]) + b (Cos[Sqrt[q2] z]/.qSol[[2]])+ c (Cos[Sqrt[q2] z]/.qSol[[3]])

eq1 = uz[1/2] == 0

eq2 =(D[uz[z],z] == 0 /. z -> 1/2)

eq3 = (Collect[D[uz[z],{z,4}] - 2 k2 D[uz[z],{z,2}] + (k2^2) uz[z], {a,b,c}, Simplify]== 0 /. z -> 1/2)

These three equations can be satisfied simultaneous only if  the determinant of a characteristic matrix vanishes. This defines the marginal condition between (k^2)and Ra.

Marginal Instability Condition

The marginal instability criterion is found by simultaneously solving three boundary conditions. This possible when the determinant of a characteristic matrix vanishes.

First[eq1]

row1 = {Coefficient[First[eq1], a],Coefficient[First[eq1], b],Coefficient[First[eq1], c]}

row2 = {Coefficient[First[eq2], a],Coefficient[First[eq2], b],Coefficient[First[eq2], c]}

row3 = {Coefficient[First[eq3], a],Coefficient[First[eq3], b],Coefficient[First[eq3], c]}

Det[{row1,row2,row3}] // Simplify

marginal[k2_,Ra_] = (Det[{row1,row2,row3}]/(k2^(2/3))/(Ra^(2/3))) // Simplify

marginal[2.0,4000.0]

Plot[Evaluate[({Re[#1],Im[#1]}&)[marginal[k2,4000.`]]],{k2,0,50},PlotRange->All,AxesLabel->{"\!\(\*SuperscriptBox[\(k\), \(2\)]\)","Det"}]

Plot[Evaluate[({Re[#1],Im[#1]}&)[marginal[1.`,Ra]]],{Ra,0,5000},PlotRange->All,AxesLabel->{"Ra","Det"}]

FindRoot[marginal[1.0,Ra], {Ra, 2000.0}] // Chop

FindRoot[marginal[5.0,Ra], {Ra, 2000.0}] // Chop

FindRoot[marginal[50.0,Ra], {Ra, 2000.0}] // Chop

marginalRa[k2_] := Ra /. Chop[FindRoot[marginal[k2,Ra], {Ra, 2000.0}]]

marginalRa[2.0]

marginalBoundary = Table[{marginalRa[k2],Sqrt[k2]}, {k2, 1.0, 50.0,1.0}];

evenMarginalPlt=ListPlot[marginalBoundary,Joined->True,PlotRange->{{0,9000},{0,9}},PlotLabel->"Marginal Instability",AxesLabel->{"Ra","k"}]

Critical Mode

The critical mode and Rayleigh number were given in the textbook.

critRa = Min[First[Transpose[marginalBoundary]]]

critk = Select[marginalBoundary,(First[#] == critRa)&][[1,2]]

eq1Crit = eq1 /. k2 -> (critk^2) /. Ra -> critRa  /. a -> 1// Simplify

eq2Crit =eq2 /. k2 -> (critk^2) /. Ra -> critRa /. a -> 1 // Simplify

modeCrit = First[Solve[{eq1Crit,eq2Crit}, {b, c}]]

uzCrit[z_] = uz[z] /. a -> 1/. modeCrit/. k2 -> (critk^2) /. Ra -> critRa

uzCrit[0.0]

Plot[Evaluate[Chop[uzCrit[z]]],{z,-0.499`,0.499`},PlotLabel->"Critical Uz[z]",AxesLabel->{"z",""}]

Odd Mode (?)

Sin[z] is Even

First, we'll look for even modes...

Plot[Sin[z],{z,-π,π}]

Plot[Evaluate[({Re[#1],Im[#1]}&)[Sin[z (1+I (Sqrt[3]))]]],{z,-π,π},PlotStyle->{Blue,Red}]

An similiar analysis can be performed for the odd modes, but these modes have a higher critical Ra number.

Satisfy Boundary Conditions

uz[z_] = a (Sin[Sqrt[q2] z]/.qSol[[1]]) + b (Sin[Sqrt[q2] z]/.qSol[[2]])+ c (Sin[Sqrt[q2] z]/.qSol[[3]])

eq1 = uz[1/2] == 0

eq2 =(D[uz[z],z] == 0 /. z -> 1/2)

eq3 = (Collect[D[uz[z],{z,4}] - 2 k2 D[uz[z],{z,2}] + (k2^2) uz[z], {a,b,c}, Simplify]== 0 /. z -> 1/2)

These three equations can be satisfied simultaneous only if  the determinant of a characteristic matrix vanishes. This defines the marginal condition between (k^2)and Ra.

Marginal Instability Condition

The marginal instability criterion is found by simultaneously solving three boundary conditions. This possible when the determinant of a characteristic matrix vanishes.

First[eq1]

row1 = {Coefficient[First[eq1], a],Coefficient[First[eq1], b],Coefficient[First[eq1], c]}

row2 = {Coefficient[First[eq2], a],Coefficient[First[eq2], b],Coefficient[First[eq2], c]}

row3 = {Coefficient[First[eq3], a],Coefficient[First[eq3], b],Coefficient[First[eq3], c]}

Det[{row1,row2,row3}] // Simplify

marginal[k2_,Ra_] = (Det[{row1,row2,row3}]/(k2^(2/3))/(Ra^(2/3))) // Simplify

marginal[50.0,4000.0]
marginal[50.0,1000.0]

Plot[Evaluate[({Re[#1],Im[#1]}&)[marginal[k2,4000.`]]],{k2,0,70},PlotRange->All,AxesLabel->{"\!\(\*SuperscriptBox[\(k\), \(2\)]\)","Det"}]

Plot[Evaluate[({Re[#1],Im[#1]}&)[marginal[50.`,Ra]]],{Ra,0,5000},PlotRange->All,AxesLabel->{"Ra","Det"}]

Plot[Evaluate[({Re[#1],Im[#1]}&)[marginal[25.`,Ra]]],{Ra,0,5000},PlotRange->All,AxesLabel->{"Ra","Det"}]

FindRoot[marginal[25.0,Ra], {Ra, 3000.0}] // Chop

FindRoot[marginal[50.0,Ra], {Ra, 3000.0}] // Chop

FindRoot[marginal[100.0,Ra], {Ra, 3000.0}] // Chop

marginalRa[k2_,start_:3000.0] := Ra /. Chop[FindRoot[marginal[k2,Ra], {Ra, start}]]

globalStart = marginalRa[25.0,500.0]

marginalBoundary = Table[{globalStart=marginalRa[k2,globalStart],Sqrt[k2]}, {k2, 25.0, 100.0,2.0}];

oddMarginalPlt=ListPlot[marginalBoundary,Joined->True,PlotRange->{{0,9000},{0,9}},PlotLabel->"Marginal Instability",AxesLabel->{"Ra","k"}]

Show[evenMarginalPlt,oddMarginalPlt]

Critical Mode

The critical mode and Rayleigh number were given in the textbook.

critRa = Min[First[Transpose[marginalBoundary]]]

critk = Select[marginalBoundary,(First[#] == critRa)&][[1,2]]

eq1Crit = eq1 /. k2 -> (critk^2) /. Ra -> critRa  /. a -> 1// Simplify

eq2Crit =eq2 /. k2 -> (critk^2) /. Ra -> critRa /. a -> 1 // Simplify

modeCrit = First[Solve[{eq1Crit,eq2Crit}, {b, c}]]

uzCrit[z_] = uz[z] /. a -> 1/. modeCrit/. k2 -> (critk^2) /. Ra -> critRa

uzCrit[0.0]

Plot[Evaluate[{Re[uzCrit[z]],Im[uzCrit[z]]}],{z,-0.499`,0.499`},PlotLabel->"Critical Uz[z]",AxesLabel->{"z",""},PlotRange->All,PlotStyle->{Blue,Red}]

Summary

The linearized fluid dynamics equations were solved to find the marginal instability boundary for Bernard thermal instability.