3.8 Numerical Solution of the Loop Closure Equations

The equations derived in the previous sections for position analysis can easily be used to obtain numerical results when the values of the link lengths are known. One can use a scientific calculator to compute the values of the dependant variables for a given value of the input parameter.  If the analysis for the whole range of the input variable is to be performed, some means of programming or a package programme may be useful. In this section examples for different ways of numerical solution will be shown using.

• Scientific calculator,
• Spread-sheet package program such as EXCEL®
• A programming language such as BASIC, PASCAL, FORTRAN or C,
• A package program for solving equations such as MathCAD® or  MATLAB®

It must be noted that the numerical solution of mechanisms is not restricted by the examples shown. Advances in the programming techniques and hardware has a great influence on the numerical methods and new software and hardware capabilities may change the approach while keeping the basic mathematics the same.

Teaching a programming language or a package program is our aim. Therefore it is strongly recommended to the reader to refer to any one of the well written texts for a programming language or package program for a complete understanding of commands and features of these programs.

Example 1:

For the mechanism shown: a₂ = 50 mm; a₃ = 250 mm; b₃ = 120 mm; a₁ = 20 mm and γ₃ = 30°. Determine the coordinates of point C when θ₁₂ = 60°.

Since:

sinθ₁₃ = (a₂sinθ₁₂ − a₁)/a₃

sinθ₁₃ = (50 sin60° − 20)/250 = 0.0932

 

which yields θ₁₃ = 5.3° or 174.7°

From the figure, we note that we must select θ₁₃ = 174.7°. Then:

s₁₄ = a₂ cosθ₁₂ − a₃ cosθ₁₃ = 50 cos60° − 250 cos 174.6° = 273.9 mm

x_C = x + b₃cos(θ₁₃ − γ₃) = 273.9 + 120 cos(174.7° − 30°) = 176.0 mm

y_C = a₁ + b₃sin(θ₁₃ − γ₃) = 20 + 120 sin(174.7° − 30°) = 89.4mm

If the calculator is programmable, the above equations can all be programmed as successive operations and the co-ordinates of point C can be obtained for different crank angles. Otherwise, the same calculation must be typed several times. The values of the fixed parameters and the input parameter must all be stored and recalled when necessary. The computation will very much depend on the calculator that is available.

Example 2:

Consider the four-bar mechanism shown which has the link lengths:

a₁ = 70 mm;  a₂ = 35 mm; a₃ = 62.3 mm; a₄ = 56 mm; |BC| = c₃ = 84.1mm; |AC| = b₃ = 126.6mm

You are to determine q₁₃, q₁₄ and the co-ordinates of point C for 0° < q₁₂< 360° with 5° increments. The necessary equations are:

x_A = s cosϕ = a₂cosθ₁₂ − a₁      (Horizontal component of AB₀)

y_A = s sinϕ = a₂sinθ₁₂            (Vertical component of AB₀)

s = \displaystyle \sqrt{{{{\text{x}}_{\text{A}}}^{2}+{{\text{y}}_{\text{A}}}^{2}}}            ( magnitude of AB₀)

ϕ = atan2(x_A, y_A)           (angular position of AB₀ with respect to positive x axis)

μ = ∠ABB₀ = cos^-1 \displaystyle \left[ {\frac{{{{\text{a}}_{3}}^{2}+{{\text{a}}_{4}}^{2}-{{\text{s}}^{2}}}}{{2{{\text{a}}_{3}}{{\text{a}}_{4}}}}} \right]

ψ = ∠AB₀B = cos^-1 \displaystyle \left[ {\frac{{{{\text{a}}_{4}}^{2}+{{\text{s}}^{2}}-{{\text{a}}_{3}}^{2}}}{{2{{\text{a}}_{4}}{{\text{s}}}}}} \right]

θ₁₄ = ϕ − ψ

θ₁₃ = θ₁₄ − μ

x_C = x_A + a₁ + b₃cos(γ + θ₁₄)            (Horizontal component of the distance A₀C)

y_C = y_A + b₃sin(γ + θ₁₃)            (Vertical component of the distance A₀C)

where:

γ = ∠ABC = cos^-1 \displaystyle \left[ {\frac{{{{\text{a}}_{3}}^{2}+{{\text{b}}_{3}}^{2}-{{\text{c}}_{3}}^{2}}}{{2{{\text{a}}_{3}}{{\text{b}}_{3}}}}} \right]

For this example we will use of EXCEL, a spread sheet program. In a spread sheet program you have a number of cells which can be considered as small boxes. Each cell is in a certain column (designated by letters) and a certain row (designated by a number). Thus a cell has a name such as A1 (first row first column) or C25 (25th row third column), etc. For each cell you can place some number or you can compute and place a new number generated from the numbers which are already placed  into other cells. During this computation besides the simple mathematical operations such as addition and multiplication, the spreadsheet has several built-in functions. For example you can square the content of the cell A1 and multiply it will the cosine of the content of the cell A2 and place the new computed number into cell A3. Now, this future can be used very effectively for the analysis of mechanisms.

Place the values of the fixed link lengths a₁, a₂, a₃, a₄ into cells A1 to A4 respectively. Into the cell A10, place the value of the angle θ₁₂ (0). Now multiply the value in cell A10 by π/180 to obtain θ₁₂ in radians and place it into the cell B10 by simply typing =A10 * PI ()/180 into the cell B10. (all the angles must be in radians). You can also convert the angle in degrees into radians by using the build-in function in Excel “RADIANS()” by typing “=RADIANS(A10)” into cell B10.  Eq (a) is typed into the cell C10 by typing “=$A$2*cos(B10)-$A$1”. The $ sign must be used for fixed lengths ($ sign in front of column or row designation fixes that row or column. If $ sign is placed on both column and row that cell is fixed. i.e. the cell address becomes an absolute address instead of relative cell address. This is important when the cells are to be copied into other cells Instead of typing the dollar sign, after the cell name is typed you can press the “F4” key). In a similar fashion, you can type in all the formulas from B10 to M10 as shown in row 8 of the spread-sheet , remembering that the formulas are to be solved from top to bottom and from left to right. Normally you will see the resulting numerical values in the cells. Once these equations are typed, if you change the value of θ₁₂ in cell A10, the remaining dependent variables will all change automatically. Thus, you can determine the value of θ₁₃, θ₁₄ , and the coordinates of point C  for any value of the input crank angle.

the value of all the dependent variables and the co-ordinates of point C for every crank angle you have selected!

In order to visualise the motion, one can use the chart option of the spread-sheet program and plot θ₁₃, θ₁₄ and μ as function of θ₁₂, or plot y_c versus x_c to determine the coupler path C as shown in the following figures

Note that one can use the closed form method instead of step-wise method as well.

path of coupler point C

To simulate the mechanism, let us also copy cells A10 to M10 into cells A6 to M6. We shall use the Control Toolbox feature of Excel. Usually the toolbox is not switched on. To switch this toolbox on, click “tools-customize” menu, place a check mark on “Control Toolbox” (You can also place a check mark on “Visual Basic” since we shall utilize this toolbox as well). First click onto the blue triangle-yellow ruler icon to switch to the design mode. There are several different buttons that can be used. We shall use the “Spin Button” which has two arrows up and down. If you are in the design mode, when you click onto the spin button the arrow cursor will change to plus sign. When you are near the cell A7 Right click on the mouse button and form a rectangular box . You will see that the spin button will be formed when you release the right button of the mouse. Right click onto this button shape you have formed and from the pop-up menu select “properties”. You can change the colour, size of the button. One important feature is that you can link this button to a cell. In the “Linked Cell” box type A6, the cell which contains the value of the independent variable (crank angle in degrees). Into the “Max” cell type “360” and into the “SmallChange” cell type an integer step size (say 2). Close the menu and on the Control ToolBox click onto the Blue Triangle icon to close the design mode. İf you click the left arrow on the spin button the value in cell A6 will decrease by 2 and if you click the right arrow, it will increase by 2 and the cells B6-M6 will also change accordingly.

Let us now go to the cell Q2 and R2 and type 0 for the coordinates of A₀. Next into cell Q3 and R3 type the coordinates of point A which is “=A2*Cos(B6)” and “=A2*Sin(B6)”. Into cells Q9, R9 type the coordinates of point B₀ (“=A1” and “0”). Type the coordinates of point B into cells Q4 and R4 as “=Q9+A4*Cos(J6)” and “R9+A4*Sin(J6)”. If you type “=L6” and “=M6” into the cells Q5, R5 you have copied the coordinates of point C into these cells. Into the cells Q7, R7 and Q8, R8 copy the coordinates of points A and B Thus obtaining a table as shown in the following figure (actually in the cells you will be seeing the result as numbers). These are the coordinates of the relevant points on the mechanism.

Next click on the chart icon. Select “XY Scatter” as the Chart type. As sub-type select “scatter with data points connected by lines” and click next. Select series and add series. For x values, select cells Q2 to Q8 which are the x values of points A₀ABB₀ and C. For y values select cells R2 to R8 which are the y values. You will obtain a chart similar to the one as shown below.

Although the four-bar mechanism is drawn, figure is not exactly satisfactory. First click on the chart and select Chart options from the pop-up menu. Select gridlines for both x and y axes. From the pop-up menu also select “format plot area”to change the colour of the plot area from grey to some other colour (say no colour). Next click on the x and y axes and on the scales menu change the minimum and maximum -40<x<100 and -40<y<140 and change major unit to 20 for both x and y axis. On the screen try to make the grid lines look as a square. In such a case you have made the x and y axis scales the same and the mechanism proportion is more real. Left Click on the mechanism link and select “Format Data Point”. In “Patterns” menu select the colour and the weight of the line that represents the links. For the marker shape select the circle and for the foreground select “No Colour”. These circles will represent the revolute joints between the links. You can also plot the coupler curve shown in Figure as another series on this chart. Now the mechanism will be seen as shown. When you click on the spin button, the input crank angle value will change and so will the positions of all the links of the mechanism drawn on the chart. You can also use different colours for each link by defining A₀A, ABCA and BB₀ as different series.

Click here to download the Excel file for this example.

The analysis and simulation of mechanisms using Excel is explained in more detail in Appendix II.

Example 3:

In this example we shall use MathCAD® which is a program for solving equations. There are other programs such as Mathematica®, Eureka® , MATLAB®, etc. which can essentially perform similar equation solving. In spread-sheet approach on the screen you see the cells and the numbers stored in the cells. In case of MathCAD the same basic idea is utilised. However the cells or boxes in which you store the numbers are not shown directly but  they exist as a label name that the user assigns. For example, when you type: a:=120 , b:=50, θ:=20, etc. The program generates boxes with labels a, b, attached to the cells in the computer which stores the values you have typed. (This assignment is done by typing the label and then typing “:”, from which you see “:=” on the screen and the computer expects an assignment to be made to the label typed).  There is also a cell with label p which already has the value 3.141592654 that you can use. Similar to the spread-sheet, besides the simple mathematical operations there are several built-in functions as well. Now, if you type: x_a:=a*cos(θ*π/180)+b on the screen you will see (to type x_a type “x.a” on the screen. To type the Greek characters you can use “Greek toolbar” or use [Ctrl] G to change any latin character to its Greek equivalent i.e typing “q [Ctrl] G “will result in “q“

A new cell with label x_a is made available and into this cell the numerical value is placed using the value of the cells a, b and q and the operations indicated. If you want to see the contents of the cell simply type x_a =. If you had a:=120, b:=50, θ:=20 assigned beforehand, on the screen you will see: x_a = 143.969.  Hence you can see the operations that you have performed as you would write it as a text, and you see the values in certain cells whenever necessary. The formulas that you type are shown in a more realistic form.

In figure (a) below, an adjustable pump is shown. Crank A₀A (which is constructed as an eccentric) is driven by an electric motor through the worm gear. The schematic drawing of the mechanism is shown in figure (b). The stroke adjustment is obtained by moving the location of the pivot point B₀ by means of an adjustment screw.

Note that A₀ABB₀ forms a four-bar mechanism. Unlike the previous examples the fixed link A₀B₀ is not horizontal and B₀ is to the left of A₀.

We can write the vector equation:

B₀A = B₀A₀ + A₀A

se^iϕ  = b₁ − is₁ + a₂e^iθ₁₂

or, equating the real and imaginary parts:

s cosϕ = b₁ + a₂cosθ₁₂
s sinϕ = s₁ + a₂sinθ₁₂

From which we can solve for s and ϕ. Then using the triangle B₀A B, the angle ψ can be determined using the cosine theorem:

ψ = cos^-1 \displaystyle \left[ {\frac{{{{\text{a}}_{4}}^{2}+{{\text{s}}^{2}}-{{\text{a}}_{3}}^{2}}}{{2{{\text{a}}_{4}}{{\text{s}}}}}} \right]

Now the angle θ₁₄ is:

θ₁₄ = ϕ − ψ − π/2

For the double slide (links 4 and 5), the loop equation is:

c₄e^iθ₁₄ + is₄e^iθ₁₄ = −is₁ − s₁₅ + b₁

 which yields:

s₄ = −(s₁ + c₄sinθ₁₄)/cosθ₁₄

s₁₅ = −c₄cosθ₁₄ + s₄sinθ₁₄ + b₁

In MathCAD you will write these equations exactly the same way. But before performing any computation, you must define the numerical values of the variables you are going to use either as input or as a result of some computation. The MathCAD sheet for the input crank angle θ₁₂ = 120° will look as follows:

---

MathCAD Output – 1
ADJUSTABLE STROKE PUMP

Determination of the output for a given crank angle.

conv := π/180 (to convert degrees into radians)

A₀A = a₂ := 55        AB = a₃ := 240        BB₀ = a₄ := 165        b₁ := 185        s₁ := 90        c₄ := 70

θ := 320        θ := θ^.conv       

x_a := b₁ + a₂^.cos(θ)        y_a := −s₁ + a₂^.sin(θ)       

s := \displaystyle\sqrt{{{{{{{{\text{x}}_{\text{a}}}}}}^{2}}+{{{{{{\text{y}}_{\text{a}}}}}}^{2}}}}

ϕ := angle(x_a, y_a)

ψ := acos \displaystyle \left[ {\frac{{{{\text{a}}_{4}}^{2}+{{\text{s}}^{2}}-{{\text{a}}_{3}}^{2}}}{{2{{\text{a}}_{4}}{{\text{s}}}}}} \right]

If we want to analyse the motion for the whole cycle, say for every 5º, then we must use indices for the variables so that the computer can allocate several cells to that variable. Secondly, besides obtaining numerical output we can also ask MathCAD to plot the result. In such a case, the solution will be as follows:

---

MathCAD Output – 2

ADJUSTABLE STROKE PUMP

Complete cycle analysis.

        

In this case i is not a subscript but an index that tells the computer to allocate an array of cells for that variable with index i. θ_i is written as q[i not q.i and refers to the ith element of the array q_i as in the case of a subscript. If we want to determine the effect of the distance s₁ which is adjusted by the screw, to the displacement of the pump, then we must perform the same analysis for different values of s. In this case we must use another index, k. The result is as follows:

---

MathCAD Output – 3

ADJUSTABLE STROKE PUMP

Analysis for complete cycle and for different strokes:

As it can be seen from the figure, as the distance s₁ is increased from 50 mm (k = 0) to 150 mm (k = 5), the stroke changes 63.8 mm (k = 0) to 155.4 mm.

You can download the Mathcad file by clicking here

Example 4:

The example solved in MathCad can also be solved in Excel. Since the number of steps in the calculation has increased, if we use the Excel sheet only, the number of columns and rows to be used will be large. Instead of using the standard build in functions, we can write a function that will evaluate s₁₆ for a given value of the input angle q₁₂ and adjustment position s₁Hence the output parameter value can be evaluated in a single step. From “Tools” menu select “Macro” and then “Visual Basic Editor” . Else you can use Alt+F11 keys. When you are in VBA editor, select “module” from the insert menu. On the module you can write ant VBA program. Let us write a program named AdjustablePump() which will utize the formulas that we have derived. This program may look like:

Global Const PI = 3.1415926
Function AdjustablePump(Crank;Coupler;Rocker;FixedLink_X;Eccentricity;B$7;$A8)
Dim Xa, Ya As Double
Dim S, fi, si, Theta2, Theta4 As Double
Dim S_4 As Double
Theta2 = ThetaDegrees * PI / 180
Xa = FixedLink_X + Crank * Cos(Theta2)
Ya = -AdjustmentLength + Crank * Sin(Theta2)
S = Mag(Xa, Ya)                                      ‘BoA distance
fi = Ang(Xa, Ya)                                       ‘Angle AB0 makes with the horizontal
Si = AngCos(Rocker, S, Coupler)        ‘Angle BB0A
Theta4 = fi – si – PI / 2                             ‘Rocker angle
S4 = -(AdjustmentLength + Eccentricity * Sin(Theta4)) / Cos(Theta4)
AdjustablePump = FixedLink_X – Eccentricity * Cos(Theta4) + S4 * Sin(Theta4)
End Function

 

You can easily follow this function routine. However this function uses three functions which are not available in Excel Library. These are Mag(X,Y) which determines the hypotenuse of a right angled triangle whose two other sides are X and Y. Ang(X,Y) is the same as ATAN2(X,Y) function. The reason we use it is that ATAN2() is not available in visual Basic. AngCos(U1,U2, Opposite) applies the cosine theorem to determine the angle opposite to the side called “Opposite” in a triangle whose other two sides are U1 and U2. Detailed explanation for these functions and more are given in Appendix-2. You must insert the Basic.Bas file into your module to use these functions .

Let us place the values of the link lengths into cells B1 to B5 and name these cells as Crank, Coupler, Rocker, FixedLink_X, and Eccentricity. Into cell B7let us place an adjustment length s₁ (s₁ = 20mm).

Now click on to the cell B8 and into this cell write. = AdjustablePump(Crank; Coupler; Rocker; FixedLink_X; Eccentricity; B$7; $A8 )(you can also use “Insert-Function-User Defined” commands to reach this function and fill in the required parameters) Once you enter this function the position of the plunger (s₁₅) will be determined for the given link lengths, the crank angle in A8 and the adjustmant length in cell B7. If you copy This cell (B8) up to B26, you will be able to determine the plunger position for every crank angle and for the given adjustmant length B7, since when during vertical copying the cell B/ does not change since we have fixed the row as B$7.

The distance s ₁ is changed by an adjustment screw shown. If we want to determine the effect of this adjustment, s ₁, on the displacement of the plunger, s₁₅ , into cells B7, C7, D7, E7, F7 and G7 we place different values of s ₁. Since the value of s ₁ in cell B7 is written as B$7 row 7 does not change (whereas in column C B$7 will be C$7). Now copy cell B8 horizontally from B8 to G8 and then verticall up to G26. Each column will refer to a particilar adjustmant s ₁ and each row will refer to a particular crank angle q₁₂ The result is given as a Excell sheet and sketch below.

You can download the Excel file.

If your company is producing such a pump with different size, everybody in the company can make use of such a file very easily!!!

	A	B	C	D	E	F	G
1	A2	55
2	A3	240
3	A4	165
4	B1	185
5	C4	70
6
7	θ12 / s1	20	50	100	150	200	250
8	0	252.0937	247.4078	249.6989	254.5738	251.8945	231.2036
9	20	252.1405	245.7816	248.1027	255.8408	258.3663	245.9161
10	40	252.2511	245.3691	249.6373	261.9564	270.8585	266.6617
11	60	252.2215	246.2394	254.6551	273.3395	289.4999	292.9677
12	80	252.0927	248.7137	263.8473	290.7546	314.6072	324.3855
13	100	252.2614	253.48	278.3194	315.0693	345.947	359.3292
14	120	253.5854	261.6704	299.246	346.0074	380.8984	393.3704
15	140	257.3465	274.4691	325.9924	379.1378	412.1455	418.5263
16	160	264.4316	291.0171	352.1259	403.739	429.014	426.5381
17	180	272.9428	305.118	365.8015	409.1448	425.2503	414.8999
18	200	277.6904	308.9373	361.3043	395.2032	404.2678	388.3058
19	220	276.148	302.3755	344.4052	370.4372	374.5207	354.1352
20	240	270.7803	291.0105	323.5058	343.1006	343.1258	318.3708
21	260	264.7198	279.2832	303.5609	317.7134	314.1433	284.8884
22	280	259.5959	269.1044	286.4366	296.0851	289.5916	256.4869
23	300	255.902	260.9577	272.5591	278.8124	270.5826	235.7459
24	320	253.6054	254.7785	261.8858	266.0913	257.8185	224.6681
25	340	252.4534	250.3422	254.2934	258.0019	251.617	223.5825
26	360	252.0937	247.4078	249.6989	254.5738	251.8945	231.2036

A note on the accuracy of computation: In the computers the numbers are stored using finite number of digits. For example in MathCAD, all the computations are performed in 15 significant digits. In certain other programs or compilers the number of digits may vary from 6 to 32. In mechanism analysis the number of significant digits used by the hardware or software is more than sufficient and very precise results can be obtained. The exceptions are when the mechanism is at or near a dead centre position, in which case the precision may decrease or an error command may be displayed (such as division by zero or negative value in the square root). This will also indicate a misbehaving mechanism in practice.

In engineering computations “No result can be more accurate than the accuracy of the input parameters”. Therefore, your result must not contain more significant digits than the number of significant digit for the  input parameters. For example, if you are to manufacture the links using classical tools in a machine shop and if the accuracy of the link lengths are within ±0.1 mm, an x displacement result cannot be x = 26.231456789 mm. This must be rounded to x = 26.2 mm.