Symbolic Math Toolbox

Modeling the Motion of an Automotive Piston

This example uses Symbolic Math Toolbox to model the motion of a simple automotive piston. 

image

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

Calculate Piston Height

A schematic of the crank and connecting rod (including relevant dimensions) is below:    

image

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 math.  At the BDC position, height relative to the crank origin is math .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) - (L-a)

math

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)"])

MuPAD graphics

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)"])

MuPAD graphics

Calculate Volume of Piston Cylinder

The maximum volume in the piston chamber occurs when the piston is at the BDC position.  At that positon, volume can be expressed as math.  In general, volume can be expressed as math, 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,`θ`))

math

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)"])

MuPAD graphics

Calculate Surface Area of Piston Cylinder

Surface area is defined similar to volume. The surface area of the cylinder head and the piston head are both approximately math, while the area of the wall is math. 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,`θ`))

math

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)"])

MuPAD graphics

Evaluate Piston Motion for Changing Angular Velocities

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, (210-105)*t],
[t > 2, 210])

math

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)

math

plot(pistHeight(150,50,angPos),t=1..2, AxesTitles=["time (sec)",
"Height (mm)"])

MuPAD graphics

Animate Piston Motion

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.

MuPAD graphics

The code used to create this animation is available here.