Nonlinear Dynamics of a Barometer: Effect of Compressibility of Air Above Mercury
​Authors: Housam Binous and Ahmed Bellagi

In an ideal mercury barometer, the space above mercury is nearly vacuum . But if a small air pocket remains, its compressibility can modify the transient dynamics significantly.
​
We propose to study the transient behavior of a mercury barometer subjected to a sudden atmospheric pressure increase (equal to 5000 Pa). To do so, we account for the variable inertia of the mercury column and the compressibility of a trapped air pocket above the fluid. As opposed to routine hydrostatic analyses, we capture the full nonlinear dynamics, including height-dependent viscous damping, adiabatic air compression, and the variable-mass inertial term. For narrow tubes, the evolving Reynolds number indicates that the flow regime is laminar throughout the process.
​
​Numerical Application: The initial height is 749.46 mm. The final height if the air compressibility is neglected is 789.44 mm. The final height is equal to 785.61 mm when air compressibility is accounted for.
In[]:=
rho=13600;(*mercurydensitykg/m^3*)​​mu=1.5*10^-3;(*dynamicviscosityPa.s*)
In[]:=
R=0.001;(*tuberadiusm*)​​A=PiR^2;(*cross-sectionalarea*)
In[]:=
g=9.81;(*gravitym/s^2*)
In[]:=
gamma=1.4;(*adiabaticindexforair*)
In[]:=
Patm0=101325;(*initialatmosphericpressure*)​​Patm1=105325;(*steppressure*)
In[]:=
H0=Patm0/(rhog)(*initialequilibriumheight*)
Out[]=
0.759467
H1=Patm1/(rhog)(*finalequilibriumheightwhenoneneglectsaircompressibility*)
Out[]=
0.789448
In[]:=
V0=5*10^-6;(*initialtrappedairvolume*)​​Pair0=500;(*initialairpressure*)
In[]:=
tmax=20;
In[]:=
(*AtmosphericForcing---weapplyastepchangeofmagnitude5000Pa*)​​Patm[t_]:=Patm0+(Patm1-Patm0)UnitStep[t]
In[]:=
(*AirCompressibility*)​​V[x_]:=V0-Ax
In[]:=
Pair[x_]:=Pair0(V0/Max[V[x],10^-12])^gamma(*adiabaticcompressionfollowsP*V^gamma=constant*)​​​​(*ViscousForce---assumeLaminarflow---Reynoldsnumberlessthan2000*)​​Fvisc[x_,xdot_]:=(8mu(H0+x)/R^2)xdotA​​​​(*Newton'ssecondlawofmotion*)​​eq=rhoA((H0+x[t])x''[t]+x'[t]^2)==A(Patm[t]-Pair[x[t]])-rhogA(H0+x[t])-Fvisc[x[t],x'[t]];​​​​(*InitialConditions*)​​ics={x[0]==0,x'[0]==0};​​​​(*NumericalSolutionusingNDSolve*)​​sol=NDSolve[{eq,ics,WhenEvent[V[x[t]]<=10^-12,"StopIntegration"]},x,{t,0,tmax},MaxSteps->Infinity];​​​​(*Totalheightofmercurycolumnequalsintialheight(H0)+rise(x)*)​​H[t_]:=H0+x[t]/.sol[[1]]​​
First[1000H[t]/.sol/.t->tmax](*exactfinalequilibriumheight*)
Out[]=
785.61
In[]:=
Reynolds[t_]:=2RrhoAbs[x'[t]/.sol[[1]]]/mu(*Reynoldsnumbers*)
In[]:=
PairvsTime[t_]:=Pair[x[t]/.sol[[1]]](*compressedairpressure*)​​
In[]:=
heightPlot=Plot[Evaluate[1000H[t]],{t,0,tmax},FrameLabel->{"t (s)","Total Height (m)"},PlotLabel->"Mercury Column Height",PlotRange->All,Frame->True,AspectRatio->1,GridLines->Automatic,PlotStyle->Darker@Red]
Out[]=
In[]:=
velocityPlot=Plot[Evaluate[x'[t]/.sol],{t,0,tmax},FrameLabel->{"t (s)","Velocity (m/s)"},PlotLabel->"Rise Velocity",PlotRange->All,Frame->True,AspectRatio->1,GridLines->Automatic,PlotStyle->Darker@Green]
Out[]=
rePlot=Plot[Evaluate[Reynolds[t]],{t,0,tmax},FrameLabel->{"t (s)","Reynolds Number"},PlotLabel->"Instantaneous Reynolds Number",PlotRange->All,Frame->True,AspectRatio->1,GridLines->Automatic,PlotStyle->Darker@Blue](*thisplotconfirmsthelaminarflowassumption​​usedtoderivethevisocusdampingtermduetotubewall*)
Out[]=
In[]:=
airPressurePlot=Plot[Evaluate[PairvsTime[t]],{t,0,tmax},FrameLabel->{"t (s)","Air Pressure (Pa)"},PlotLabel->"Trapped Air Pressure",PlotRange->All,Frame->True,AspectRatio->1,GridLines->Automatic,PlotStyle->Darker@Magenta]​​
Out[]=