Modeling a Mechanical System by Numerical Solutions to
Differential Equations

Description

This Notebook shows how Mathematica can be used in the modeling of mechanical systems; it involves a simple example representing automobile shock absorbers.

By: John M. Novak
Wolfram Research, Inc.

Model

The system we investigate here is a simple fourth-order mechanical network.  Naturally, there are many approaches to working with such a system.  We chose to use the numeric differential equation solver in Mathematica to investigate this situation.

Show[Graphics[{
	GrayLevel[0],Disk[{0,0},1],
	GrayLevel[1],Disk[{0,0},.4],
	GrayLevel[0],
	Line[{{1.1,0},{2,0},{1.9,1.3},{-1.4,1.7},
		{-1.7,3.0},{-2.2,3.0}}],
	Line[{{-1.1,0},{-2.2,0}}],
	Line[{{0,-1.1},{0,-2.0}}],
	Polygon[{{0,-1.1},{.2,-1.6},{-.2,-1.6},
		{0,-1.1}}],
	Text["F",{0,-2.3}]}],
	AspectRatio->Automatic]

Consider a simplistic view of an automobile suspension.  We have a pair of masses, the car and the wheel, coupled by a spring and a damper (the shock absorbers).  In addition, the wheel is coupled to the road by a spring and a damper (the tire).  We choose here to model the road as an oscillatory force on the tire.  We then observe the oscillations of the wheel and car, with a view to minimizing the effect on the car, while not straining the systems coupling the wheel to the car and road.  (See diagrams.)  Again,  it must be emphasized that this is a simplistic model-- much is  omitted or just plain not right in this model, but it does give the general idea.

Mathematically, we can now express the equations of motion of the wheel and car by summing the forces on each mass.  Here, the resting length of each of the springs is .5 m.  The first equation represents the forces on the upper mass, the second on the lower.  f[t] is our driving force, the force the road exerts on the wheel.

eqn1 =  m1 x1''[t] - k1(.5 - (x1[t] - x2[t])) +
	d1(x1'[t] - x2'[t]) == 0;

eqn2 =  m2 x2''[t] - k1((x1[t] - x2[t])- .5) -
	d1(x1'[t] - x2'[t]) - k2 (.5 - x2[t]) +
	d2 x2'[t] == f[t];

We also need to specify boundary values for the equations.  Here, assume there is  no initial acceleration or velocity, with the wheel  at 1/2 m and the car at 1 m.

bvals = {x1[0] ==1, x1'[0] == 0, x1''[0] == 0,
	x2[0] == .5,x2'[0] == 0, x2''[0] == 0};

These equations we just assembled are entirely symbolic.  To use them, we need a simple mechanism for substituting the values for our constants, that will at the same time allow us to easily modify the equations.  This mechanism can also be used to pass in the equations representing the driving force.  We can do this in the following fashion:

eqns[em1_,em2_,ek1_,ek2_,ed1_,ed2_,force_] :=
	Flatten[{eqn1,eqn2,bvals}]/.
	{m1->em1,m2->em2,k1->ek1,k2->ek2,
		d1->ed1,d2->ed2,f[t]->force}

This is all we need to start investigating solutions.  As a note, I assume that consistent units are used everywhere here - these were generated with MKS units in mind, so our force, for instance, ranges approximately between positive and negative 1/3 of the force of gravity.

For convenience, we can define a simple function to represent an impulse--a square pulse starting at some particular time, and lasting for a defined time.

Pulse[t_?NumberQ,tmin_,length_] :=
	If[t >= tmin && t <= tmin + length,
		1, 0]

Solutions

Now that we have our equations, we want to see solutions.  Here is a simple example:

soln1 = NDSolve[eqns[1000,15,1000,8000,1000,800,
			1000 Pulse[t,0,.5]],
	{x1[t],x2[t]}, {t,0,20}];

This results in a pair of substitution rules, representing the solutions (this format is best if multiple solutions are generated).  I suppressed the actual print form for space saving purposes.  To display these, then, we need to extract the functions from the rule form, then plot them.  This can be done as follows.

Plot[Evaluate[First[{x1[t],x2[t]}/.soln1]],{t,0,20},
	Frame->True,PlotRange->{0,1.3}]

Here, we can generate a number of solutions, and then display them together with a scaled representation of the impulse:

gr1 = Plot[Evaluate[Flatten[{(1000 Pulse[t,0,.5])/3000,
		First[{x1[t],x2[t]}/.soln1]}]],{t,0,20},
		PlotStyle->{Hue[0],GrayLevel[0],
			GrayLevel[0]},
		Frame->True,PlotRange->{0,1.3},
		DisplayFunction->Identity];

soln2 = NDSolve[eqns[1000,15,1000,8000,1000,800,
			2000 Pulse[t,0,.5] + 3000 Pulse[t,3,.6] +
			1500 Pulse[t,5,.4]],
	{x1[t],x2[t]}, {t,0,20},
		MaxSteps->3000];

gr2 = Plot[Evaluate[Flatten[{
	(2000 Pulse[t,0,.5] + 3000 Pulse[t,3,.6] +
	1500 Pulse[t,5,.4])/3000,
	First[{x1[t],x2[t]}/.soln2]}]],
	{t,0,20},
	PlotStyle->{Hue[0],GrayLevel[0],
			GrayLevel[0]},
	Frame->True,PlotRange->{0,1.3},
	DisplayFunction->Identity];

Show[GraphicsArray[{{gr1},{gr2}}],
	DisplayFunction->$DisplayFunction]