Symbolic Math Toolbox 

This example uses Symbolic Math Toolbox to model the motion of a simple automotive piston.
Where:
TDC and BDC = top dead center and bottom dead center
B = bore (i.e., diameter of the cylinder)
L  length of the connecting rod
S  stroke length
a  crank radius
θ  crank angle
A schematic of the crank and connecting rod (including relevant dimensions) is below:
We would like to define piston height relative to the BDC position. The height relative to the crank origin can be calculated through trigonometry as . At the BDC position, height relative to the crank origin is .To calculate piston height relative to BDC position, we simply subtract the terms.
pistHeight := (L,a,`θ`) > a*cos(`θ`) +
sqrt(L^2  a^2*sin(`θ`)^2)  (La)
We plot piston height as a function of crank angle (θ) for one revolution, assuming a connecting rod length (L) of 150 mm and crank radius (a) of 50 mm.
plot(pistHeight(150,50,`θ`),`θ`=0..2*PI,
AxesTitles=["Crank angle (rad)", "Height (mm)"])
As expected, piston height is highest when crank angle equals 0 and 2 π.We can also create a surface plot to visualize how height changes with both crank angle and crank radius.
plot(#3D, pistHeight(150,a,`θ`),a=30..60, `θ`=0..2*PI,
AxesTitles=["Crank radius (mm)","Crank angle (rad)","Height (mm)"])
The maximum volume in the piston chamber occurs when the piston is at the BDC position. At that positon, volume can be expressed as . In general, volume can be expressed as , where H is piston height relative to BDC position (defined above). We define an expression for volume, substituting stroke length (S) with 2a:
pistVol:= (L,a,`θ`,B) > PI*(B/2)^2 * (2*a  pistHeight(L,a,`θ`))
We now plot volume as a function of crank angle (θ) for one revolution, assuming the following parameters:
L (length of connecting rod) = 150 mm
a (crank radius) = 50 mm
B (bore) = 85 mm
plot(pistVol(150,50,`θ`,85),`θ`=0..2*PI,
AxesTitles=["Crank angle (rad)","Volume (mm^3)"])
Surface area is defined similar to volume. The surface area of the cylinder head and the piston head are both approximately , while the area of the wall is . We define an expression for surface area, again substituting (S) with 2a:
pistSA:= (L,a,`θ`,B)> 2*PI*(B/2)^2 + PI*B*(2*a  pistHeight(L,a,`θ`))
We plot surface area for a piston with dimensions defined in the previous section.
plot(pistSA(150,50,`θ`,85),`θ`=0..2*PI,AxesTitles=["Crank angle (rad)",
"Surface area (mm^2)"])
Assume the crank rotates at 1000 rpm (105 rad/sec) for 1 second, then steadily increases from 1000 to 2000 rpm for 1 second, and then remains at 2000 rpm. We define a piecewise function to define angular velocity.
angVel := t > piecewise([t<=1, 105],[t>1 and t<=2, (210105)*t],
[t > 2, 210])
We calculate crank angle by integrating angular velocity.
angPos := int(angVel(t),t)
We substitute the resulting crank angle into our expression for piston height, and plot piston height from t=1 to t=2, where the angular velocity is steadily increasing. Note the increase in frequency.
pistHeight(150,50,angPos)
plot(pistHeight(150,50,angPos),t=1..2, AxesTitles=["time (sec)",
"Height (mm)"])
We animate piston motion for one full revolution (θ=0 to θ=2π) for a piston with connecting rod length (L) = 150 mm, crank radius (a) = 50 mm, and bore (B) = 85 mm. This animation was created in a separate notebook and copied into this notebook. The ability to copy graphics and animations between notebooks is useful when users want to document supporting analysis without having the code displayed in the notebook.
The code used to create this animation is available here.