J. Feicht 


This Mathematica notebook generates the beamtube pressure profile using the btwaterleak methodology of particle diffusion. The pressure  distribution is the solution to the set of n-coupled differential equations where “n” represents the number of  sections.

Options include the pump arraignments, leaks, temperature profile, etc.  A switch is included to force NDSolve to use the 4th order Runge-Kutta technique to solve the ODE’s, although doing this is not recommended. There are some asides in this version that discuss how to solve initial value problem, also setting up coupled differential equations.

Any data input or changes to defined values should use Mathematica floating point notation, for example use “2.0”, or “2.” instead of symbolic notation “2”, doing this makes the code run much faster because the code interprets this as a number.

Note: Deleting all code output before saving will greatly reduce file size.

Note: It is always a good idea to quit then restart the kernel after changing variables. Mathematica has a good (too good) memory.

Define tube geometry, pump speeds and pump locations. Warning mixed units!

ClearAll  (*See note on reseting the kernel, ClearAll etc. are not always throrough *)

ClearSystemCache

Tubelength=300*100.(*300 meter converted to centimeter*)

Tubediameter=10*2.54(*nominal 10" diameter converted to centimeter*)

Tubearea=((π (Tubediameter^2))/4.); (*(centimeter^2)*)

IonpumpS=1000*100. (*500 liter/sec ion pump speed converted to cc/sec*);
(*TempCryopumpS=1000.*6.5 (10^(3 ))(*speed cc/sec*)(*these are the temporary TempCryopumpS*)*)

Nopump=0.0(*speed cc/sec*) (*use this or code in 0.0 when temporary pumps are removed*);
Turbospeed=1000*250.0 (*speed of turbo in cc/sec*);

(*ChamberSpeed=500*1000(*speed cc/sec*)(*these are the end station speeds*)*)

Populate the Number of Ports (i.e. number of pumps)

Numberofports=7(*need to have at least 2 ports, one on each end, this value is the number of pumps too*)

Portpitch=Tubelength/(Numberofports-1)  (*centimeter*)

Populate the pump speeds and locations (comment this section out if you want to directly input pump speeds)

Do[Port[i]={(i-1)*Portpitch,
IonpumpS},{i,1,Numberofports}]; (*set pump vector and speeds*)

Do[Print[Port[i]],{i,1,Numberofports}];  (*output is a vector (ordered pair) reported as the distance from "one end", then what the pump speed is*)

Plot pump speeds (cc/sec) along beamtube

portspeedtable=Table[Port[j],{j,1,Numberofports}]; (*first need to make a table of the indexed ordered pairs*)

ListPlot[portspeedtable, {Joined->False}, AxesLabel->{"Distance", "Speed (cc/sec)"}, Filling->Bottom]

Define starting pressure, temperature and gas amu

startpress=1.0 (10^(-3)) (*Starting pressure in Torr*)

amu=18. (*water*)

Tkelvin=300. (*Kelvin*)

Define physical constants, conversion factors, variable declaration, etc.

alpha=0.5   (*accomodation coeff*)
rp=0.7   (*potential height fraction*)
sigma0=150.    (*number of monolayers at start*)
Tao=1.0 (10^(-12))   (*sec*)
tresmin=1.0 (10^(-6))   (*sec*)
tresmax=1.0 (10^10)   (*sec*)
emtmin=Log[tresmin]
emtmax=Log[tresmax]
vt=5.263 (10^4)Sqrt[(18/amu)]Sqrt[(Tkelvin/296)]  (*thermal velocity cm/sec*)
particlescm3=(296./Tkelvin)3.0 (10^16) Ptorr
particlescm2=1.0 (10^15)monolayers

pipeconductance=3.81*((Tubediameter^3)/Tubelength) ((((Tkelvin/amu)))^(1/2)) (*conductance for a long round pipe in liters/second*)

timeconstant=IntegerPart[(((Tubearea*Tubelength))/(1000pipeconductance))](*get an estimate of the system time constant (volume/(liter/sec))*)

Define how many sections to break the calculation up into, normally use 100 sections, higher values cost solve time but increase accuracy.

Sections=100(* this is the INTEGER number of calculation sections that the beamtube will be broken into, i.e. the number of coupled differential equations to solve, affects pump placement calc, etc.*)

Establish which tube sections have vacuum pumps,  then map pump location & speed to the “n” tube sections.

(*You may ask why we need this section? The reason is that because the tube can be broken up into many sections(user defined), the particular section where the placement of pumps, leaks, etc. can change. This code puts the pumps in their correct section.*)

Seclength=(Tubelength/Sections) (*length of each "calculation" section*)

Sections(*how many sections were requested*)

$MachineEpsilon

sectionrange[1]={1,Interval[{0,Seclength}]}(*defines the first section*)

sectionrange[Sections]={Sections,Interval[{(Tubelength-Seclength+.0001),Tubelength}]}(*defines the last section*)

Do[sectionrange[i]={i,Interval[{((i-1)*Seclength+.0001),(i*Seclength)}]},{i,2,Sections-1}] (*defines the intermediate sections*)(*the 0.0001 is a suck-ass way to stop MemberIntervalQ from duplicating pumps by not allowing the intervals to touch, there is probably a switch in the code to fix this but I havent found it yet*)

Table[sectionrange[i],{i,1,Sections}] ;(*take a look at the sections and their range*)

Do[findport[i,n]=Boole[IntervalMemberQ[sectionrange[i][[2]],Port[n][[1]]]],{n,1,Numberofports},{i,1,Sections}](* find the section (i) that has the pump port(n)*)

Do[tablel[j]=Table[findport[i,j],{i,1,Sections}],{j,1,Numberofports}] (*intentionally misspelled*)

portboole=Sum[tablel[i],{i,1,Numberofports}] ;(*make a boolean table of which section has the port*)

Check that the pump distribution is correct and no overlaps have occurred

Total[portboole]==Numberofports (*these numbers must be equal/True (which means the number of ports and the boolean sum are the same) or there is an overlap ant that port is being counted twice*)

Map the pump speeds to each port

Speedy[x_]:=IonpumpS*x   (*construct a function that multiplies by pump speed*)
Do[g=MapAt[Speedy,portboole,;;],{n,1,Sections}]; (*map the speed over the table*)

ListPlot[g,Filling->Bottom]

Do[speed[i]=g[[i]],{i,1,Sections}]  (*fyi, [[i]] is a compact notation for Part, it picks out that element in the list. We need these to be individual values (unlike Pumps) so they can go into the ODE by section. Speed should be in cc/sec*) 

(*>>>>>>>>>Note this will overwrite any earlier speed definitions so make a choice how you want to do this<<<<<<<<*)

speed[1001] (*speed in cc/sec*)

speed[101]  (*check if working correctly*)

speed[1]=500.*1000*0.0

speed[Sections]=500.*1000*0.0

speeds=Table[speed[i],{i,1
,Sections}] ;(*allows a quick look to be sure it's working*)

ListPlot[speeds,Filling->Bottom,Joined->True] (*speed plot versus section number in cc/sec*)

Setup beamtube leak(s). This section can be developed further as required. A similar section with programmed 
temperature also needs to be developed.

(*Leaklocation=2 (10^5)  (*centimeters from vertex*);*)

(*Sectionnumber=(Leaklocation/Tubelength)*Sections; (*finds tube section with leak*)*)

(*Leakrate=1 (10^(-5))(1000) (*leak rate in torr-liter/sec converted to torr-cc/sec*);*)

(*Leakstartday=500;*)

(*Leakstopday=501;*)

Do[qlk[i]=0.0,{i,1,Sections}] (*sets a leak rate in each section, qlk is an variable in each of the differential equations*)

(*qlk[1]=1.0 (10^(-5)) way to short-circuit the previous Do-loop and set the value of qlk[1] to do testing*)

Leaks=Table[qlk[i],{i,1
,Sections}];

ListPlot[Leaks,Filling->Bottom] (*speed plot versus section number in cc/sec*)

Set time step size

timestep=recordstep*24*3600/1440  (*timesttep in seconds >>>>>>>  1440 is not remarkable other than it gives 60 second intervals*)

RKstep=1 (*RungaKutta step size*);

Load solver

Needs["DifferentialEquations`NDSolveUtilities`"];

Set up the derivative coefficients

aa=((2 vt .5Tubediameter)/(3 (Seclength^2)))  (* units are (cm/sec)*(cm/(cm^2))=(1/sec)*)
bb=(4/Tubediameter) (* units are (1/cm)*)
dd=(1/(π (((.5Tubediameter))^2) Seclength)) (* units are (1/(cm^3))*)

Get the 4th order Runge-Kutta coefficients (it is better to let the code choose the solver, but we can force using RK if we want)

(*the following are the R-K 4th order coefficients*)

crkamat={{1/2},{0,1/2},{0,0,1}};
crkbvec={1/6,1/3,1/3,1/6};
crkcvec={1/2,1/2,1};ClassicalRungeKuttaCoefficients[4,p_]:=N[{crkamat,crkbvec,crkcvec},p];

Set  the outgassing coefficients of each section.

Do[(Subscript[ajj, n])=3.0 (10^(-7)),{n,1,Sections}]; (*outgassing per section in ((torr-(cm^3))/((cm^2)-sec)), typical values for unbaked stainless steel are on the order of 6(x10^(-9))((Torr-liter)/(sec-(cm^2))) or since 1000(cc/liter)= 6(x10^(-6))((Torr-cc)/(sec-(cm^2)))*)

(*(Subscript[ajj, 1])=1.0 (10^(-5))*)

Set up the n-coupled differential equations to get pressure profile.

diffeqns=Table[{(Subscript[p, i+1])'[t]==aa((Subscript[p, i])[t]+(Subscript[p, i+2])[t]-2(Subscript[p, i+1])[t])+bb (Subscript[ajj, i+1])+dd (qlk[i+1]-speed[i+1]*(Subscript[p, i+1])[t]),(Subscript[p, i+1])[0]==  startpress},{i,1,Sections-2}] ; (*this generates the ODE's for sections n+1 to n-1 of the n beam tube sections*)

firstsection={(Subscript[p, 1])'[t]==aa((Subscript[p, 2])[t]-(Subscript[p, 1])[t])+bb (Subscript[ajj, 1]) +dd (qlk[1]-speed[1]*(Subscript[p, 1])[t]),(Subscript[p, 1])[0]==startpress}  (*this generates the ODE for section 1 of the n beam tube sections*)

lastsection={(Subscript[p, Sections])'[t]==aa((Subscript[p, Sections-1])[t]-(Subscript[p, Sections])[t])+bb (Subscript[ajj, Sections]) +dd (qlk[Sections]-speed[Sections]*(Subscript[p, Sections])[t]),(Subscript[p, Sections])[0]==startpress} (*this generates the ODE for the last section of the n beam tube sections*) (* note that btwaterleak has a section 101, this needs more investigation, may be required here too*)

temp=Prepend[diffeqns,firstsection];  (* add all the sections to ODE list*)
diffeqns=Append[temp,lastsection];

vars2=Table[(Subscript[p, i])[t],{i,1,Sections}]; (*define the variables*)

Solve the coupled ODE’s. Can force solver to use 4th order RK  if desired (see notation). There is a subtlety here, the time that the solution runs over only needs to be the settling time of the system, not how long you want to simulate the pumpdown.  At sound speeds this is a relatively short time.  If you want to simulate a water desorption rate it needs to be wrapped around the ODE  solution with a Do statement or similar. For example the temperature, outgassing rate, or leak rate could all be parameterized around the ODE solution.

btw=NDSolve[diffeqns,vars2,{t,1,timeconstant}] ;(*let Mathematica choose the solver*)

(*btwRK=NDSolve[diffeqns,vars2,{t,1,timeconstant},Method->{"ExplicitRungeKutta","DifferenceOrder"->4,"Coefficients"->ClassicalRungeKuttaCoefficients},StartingStepSize->.01, MaxSteps->1000000]*) (*force using Runge Kutta method*)

Plot the beamtube pressure distribution. Example of pumpdown with constant outgassing rate/section

ListPlot3D[Table[vars2/.First[btw],{t,1,timeconstant,10}],AxesLabel->{"Section Number","Time","Pressure"}]

ListPlot3D[Log[Table[vars2/.First[btw],{t,1,timeconstant,10}]]]  (*try a log plot*)

(*ListPlot3D[Table[vars2/.First[btwRK],{t,1,timeconstant,10}],AxesLabel->{"Section Number","Time","Pressure"}]*)

(*ListPlot3D[Log[Table[vars2/.First[btwRK],{t,1,timeconstant,10}]]]  (*try a log plot*)*)

secpresvstime=Table[vars2/.First[btw],{t,1,timeconstant,1}]; (*these are the solutions, this can be a really big file, it contains a value for every time set from 1 to timeconstant, for each section (typically 100 sections)*)

Length[secpresvstime] (* determine the file size*)

First[secpresvstime]; (*take a look at first file elements. This is the pressure in each section at start. You should see a lot of initial conditions*)

Last[secpresvstime] ;(*take a look at last file elements*)

secpresvstime[[timeconstant]]; (*this is another way to get the last set of data in the file, or at any time interval*)

finalpressures=Table[secpresvstime[[i]][[Sections]];,{i,1,timeconstant}]; (*pull out the last sections pressure vs time at time=timeconstant*)

secpresvstime[[timeconstant]][[1Sections]]  (*pull out the final pressure in the last section at time=timeconstant*)

secpresvstime[[timeconstant]][[1]] (*pull out the final pressure in the first section at time=timeconstant*)

sec1pvst=Table[secpresvstime[[i]][[1]],{i,1,timeconstant}];  (*get the pressure vs time history of section 1*);
seclastpvst=Table[secpresvstime[[i]][[Sections]],{i,1,timeconstant}];(*get the pressure vs time history of the last section *);

ListLogPlot[{sec1pvst,seclastpvst},Joined->True, PlotLegends->{"First Section","Last Section"}]

secpresvstime[[timeconstant]][[1]] (* final pressure in the 1st section, change the 1 to whatever section you want to query*)

secpresvstime[[timeconstant]][[Sections]]   (* final pressure in the nth or last section*)

deltaP=secpresvstime[[timeconstant]][[1]]-secpresvstime[[timeconstant]][[Sections]]

conductance=((π Tubediameter*Seclength((Subscript[ajj, 1])))/deltaP)*0.001  (* (cm^2)x((torr (cm^3))/((cm^2)sec)) x (1/torr) x ((.001 liter)/(cm^3))=(liter/second) *)

pipeconductance (*this is what the handbook formula gives as the conductance in liter/sec, use for comparison*)

p8=LogPlot[(Subscript[p, 1])[t]/.First[btw],{t,1,timeconstant},PlotRange->All,GridLines->Automatic] (*this allows you to plot the pressure in a section (Subscript[p, n])*)

p50=LogPlot[(Subscript[p, 50])[t]/.First[btw],{t,1,timeconstant},PlotRange->All,GridLines->Automatic,PlotStyle->{Red,Dashed}] (*this allows you to plot the pressure in each section (Subscript[p, n])*)

Show[p8,p50]

ListLogPlot[Table[vars2/.First[btw],{t,1,timeconstant,10}],AxesLabel->{"Section Number","Pressure"},Joined->True]

ListLogPlot[Table[vars2/.First[btw],{t=timeconstant}],AxesLabel->{"Section Number","Pressure"},Joined->True]

Portpitch/100.  (*pump to pump distance in meter*)