Pulsatile Flow

1) Replace the right hand side (RHS) of the "p" equation with your list of pressure values enclosed in {}.

2) Replace the RHS of the "twave" equation with the twave time (the integer which corresponds to the pressure data point where the twave ends).

3) Replace the RHS of the "period" equation with the actual period of the heartbeat in seconds.

4) Replace the RHS of the "volume" equation with your volume function, which describes the total volume of blood expelled into the aorta as a function of time.

5) If desired, replace the RHS of the "k" equation with the value of the spring constant at minimum radius.

6) If desired, replace the RHS of the "rmin" equation with the minimum aortal radius.

Note: The above may change for arythmias with extrordinarily large systolic or extrordinarily small diastolic pressures.

7) For truncated or multiple ventricular contraction arythmias, replace the RHS of the "stroke" equation with the stroke volume.

8) Push the "enter" key after clicking on each cell you want to evaluate.

If you are not happy with the "smoothed pressure" graph, increase the RHS of the "terms" equation.

Smooth[x_List,v_,t_]:=Block[

	{y,l,loy,hiy,ser,loser,hiser},

	y=Fourier[x];

	l=Length[x];
	
	loy=Min[x];
	
	hiy=Max[x];

	ser=N[Pi Re[Sum[y[[i]]*
					Exp[-2*Pi*I*(v-1)*(i-1)/l],
				{i,1,t}]]/Sqrt[l]];

	y=Table[ser,{v,num}];
	
	loser=Min[y];
	
	hiser=Max[y];
	
	N[loy+(hiy-loy)(ser-loser)/(hiser-loser)]];

p=Table[N[80+30((t/6)^3)(E^-(t/6))],{t,64}];

twave=32;

period=60/72;

k=1.05*(10^6);

rmin=1.15;

stroke=80;

terms=6;

Use one of the following cells to define and plot your volume function (to be sure it's what you want).  Delete the other cells before evaluating the notebook.  Adjustments to the numerical constants may be necessary to make your volume fit the expected shape.  It is important that the peak(s) of your volume function occur at the end of the twave(s).

volume=If[t<twave/3,
			stroke(1-(Abs[t-twave/3]/(twave/3))^3)/2,
		If[t<2twave/3,
			3stroke(1-(Abs[t-2twave/3]/(twave/3))^3)/10,
			2stroke(1-(Abs[t-twave]/(twave/3))^3)/10]];

Plot[volume,{t,0,twave},PlotRange->All,
	PlotLabel->"Your Volume Function"];

volume=If[t<twave/2,
			2stroke(1-(Abs[t-twave/2]/(twave/2))^3)/3,
			stroke(1-(Abs[t-twave]/(twave/2))^3)/3];

Plot[volume,{t,0,twave},PlotRange->All,
	PlotLabel->"Your Volume Function"];

volume=stroke(1-(Abs[t-twave]/twave)^3);

Plot[volume,{t,0,twave},PlotRange->All,
	PlotLabel->"Your Volume Function"];

The following does all of the work - Don't change anything here!

dvol=N[volume-(volume /. t->t-1)];

vol:=If[dvol>0,dvol,0];

Now we create a smooth function of time from your pressure data.

ListPlot[p,PlotRange->All,
		PlotLabel->"Your Pressure Data"];
		
num=Length[p];

diastolic=Min[p];

pr=Smooth[p,t,terms];

Plot[pr,{t,0,num},PlotRange->All,
	PlotLabel->"Smoothed Pressure Data"];

We now fit a non-Hookean Young's Modulus to empirical curves for your value of the initial "spring constant" k (at minimum radius).

ymt={{rmin,0},{1.25rmin,1},{1.3rmin,3}};

ym[rd_]=k*(E^Fit[ymt,{1,rd},rd]);

Plot[ym[rd],{rd,rmin,rmin+.1},
	PlotRange->{0,4*10^6},
	PlotLabel->"Young's Modulus of Aorta Wall"];

Laplace's Law for the circumferential tension in a cylindrical wall of radius r and infinitesimal thickness is T = Pr.  If the wall thickness is smaller than the radius by a factor of at least 10, a good approximation is T = Pr / thickness.  The thickness of the aorta is about 2 mm.  The resulting tension divided by the Young's Modulus gives the strain (fractional extension), from which we compute the radius.

rad[r_]:=N[rmin(1+(pr-diastolic)1333r/(0.2ym[r]))];

radii=Table[Nest[rad,rmin,Round[t]],{t,1,num}];

PrependTo[radii,N[rmin]];

We now find a smooth function of t for the radius.

r=Smooth[radii,t,terms];

Plot[r,{t,0,num},PlotRange->All,
	PlotLabel->"Aortal Radius"];

Finally, we find the cross sectional area of the aortal opening, the pressure difference across and the effective force applied to that area,  and plot the velocity (power divided by force) of our parcel of blood.  vcrit is the velocity at which the Reynold's Number exceeds 2000.

a=Pi*r^2;

dp=pr-(pr /. t->t-1);

f=dp*a;

vcrit=2000*.04/(2rmin*1.05);

Plot[{pr,dp*vol/(f*period/num),vcrit},
	{t,1,twave},PlotRange->All,
	PlotStyle->{RGBColor[0,1,0],RGBColor[0,0,1],
					RGBColor[1,0,0]},
	PlotLabel->"Pressure (green) and Velocity (blue)"];

To complete the project, find the length of time (from the graph above) which the velocity exceeded the critical value.  If the velocity curve is "bumpy", treat each bump separately.  Replace the constant below with your value.

tn=15;

Now we compute the actual time of nonzero velocity, the blood volume ejected during that time, the kinetic energy it had, the eddy advection time, the eddy length scale, and the expected eddy length:

Print["delta time = ",dt=N[period*tn/num],
		" seconds."];

Print["volume = ",eddyvol=N[Sum[vol,{t,1,tn}]],
		" cubic centimeters."];

Print["kinetic energy = ",
		kin=0.5*1.05N[Sum[vol*(dp*vol/(f*period/num))^2,
						{t,1,tn}]],
		" ergs."];

Print["eddy advection time = ",
		ta=Sqrt[.04dt*eddyvol/kin],
		" seconds."];

Print["eddy length scale = ",
		le=Sqrt[Sqrt[(.04^3)dt*
		eddyvol/((1.05^2)kin)]],
		" centimeters."];

Print["expected eddy length = ",
		lexp=le*Sqrt[dt/ta],
		" centimeters."];