MSU

LIBRARIES
.—:—.

RETURNING MATERIALS:

 

 

 

Place in book drop to
remove this checkout from
your record. flfl§§ wil]
be charged if book is
returned after the date
stamped below.

 

 

 

 

THE COMPUTER-AIDED DESIGN OF
MOLDBOARD PLOW SURFACES USING
THREE DIMENSIONAL GRAPHIC TECHNIQUES
By

Stephen Brian Richey

A THESIS

Submitted to
Michigan State University

in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Agricultural Engineering

1987

L./l ti,-

Approved by:

 

Major Professor

p1 / ‘7 ”‘5 5 1'.) I 4‘
3W3jyéir>w @131;

 

Department Chairman

 

ABSTRACT
THE COMPUTER-AIDED DESIGN OF

MOLDBOARD PLOW SURFACES USING
THREE DIMENSIONAL GRAPHIC TECHNIQUES

By
Stephen Brian Richey

A computer model has been developed to show that three
dimensional graphic techniques can be used to describe
moldboard plow surfaces. A three dimensional soil-tool
mechanics was developed to calculate the forces applied to
the plow by the soil. The forces exerted on the plow by the
soil are predicted without prior knowledge of the soil path
over the surface.

The plow shape is defined by two Bezier surfaces. Each
surface can be divided into a maximum of 100 elements.
Force calculations are made for each element and then the
forces and moments are summed over all the elements to
determine the totals for each of the x, y and 2 directions.
The soil path is determined by the surface Shape and the
forces exerted on the soil.

The force predictions were checked against measurements
made by Emmet and Gunkel (1983) on a 16 inch Ford moldboard
plow. Predicted draft forces were 15 percent larger than

those measured.

This thesis is dedicated to Mary.

11

ACKNOWLEDGMENTS

 

The author gratefully acknowledges Dr. AJit Srivastava,
major professor, for his guidance of the research effort and
his patience and understanding while this manuscript was
completed in the two years that followed.

Many thanks to Dr. Larry Segerlind and Dr. Robert
Wilkinson, committee members, for their support and advice
while working towards the completion of this graduate
degree.

There are many people that contributed to this research
in small but important ways. Their assistance made the
research easier and the daily activities enjoyable. The
author specifically wishes to recognize Cliff Zehr, Tim
Oliver, Wilbur Mahoney, Tom Esch and Luke Reese. I

Deepest thanks are expressed to the author's parents,
Fred and Arlene Richey, for their continual encouragement
and support in reaching for my goals.

The person who made the most personal sacrific for this
research was the author's wife, Mary. For her sacrific,
encouragement and assistance the author extends his deepest

heart-felt graditude.

iii

TABLE OF CONTENTS

 

List or Tables. 0 O O O C O O O O 0

List of F1

sures C 0 O O O O O O O 0

List of Symbols . . . . . . . . . .

Chapter
1

Appendices
A

DOD)

List of Re

Introduction. . . . . . .
1.1 Computer Graphics. .
1.2 Soil-tool Mechanics.
1.3 Research Objectives.

Literature Review . . . .
2.1 Plow Description . .
2.2 Soil-Tool Mechanics.

Soil-Tool Mechanics Model

3.1 Two Dimensional Model.
3.2 Three Dimensional Model.

Graphics Representation .

O O O O

A. 1 Three Dimensional Curves
A. 2 Three Dimensional Surfaces

Computer Model. . . . .
5.1 Moldboard Plow Model
5.2 Computer Program . .

Results and Verification.
6.1 Inclined Blade . . .
6.2 Moldboard. . . . . .

Conclusions . . . . . . .

Suggestions for Further Research.

PLOW Users Guide
PLOW Source Code
PLOW Data Files
Measured Data

ferences

iv

Page
. v

. vi
viii

sit-l:- NNHH

. 10
. 21

. 27
. 27
. 3A
. 37
. 37
. A1
. AS

. 50

6.2

6.3

LIST OF TABLES

Page
Inclined Blade Operating Conditions, Measured
Forces and Predicted Forces. . . . . . . . . 48

Soil and Operating Conditions Used for
MOdel verification O I O O O O O O O O O O O 51

Comparison Between the Predicted Plow Forces
and the Measured Forces. . . . . . . . . . . 58

3.2
3.3
30“

3.5

3.6
3.7

3.8

3.9

3010

h.3
A.A

5.1
5.2

5.3

LIST OF FIGURES

 

Thomas Jefferson's Framework for Designing a
Moldboard Plow. . . . . . . . . . . . .

Forces on a $011 Segment Reacting on an
Inclined Plane Tillage Tool . . . . . .

A Mohr Envelope of Soil Internal Stresses. .
A Mohr Envelope of Soil-Metal Stresses . . .

Formation of Soil Blocks by Repeated Shear
Failure 0 O O O O O O O O O O O O O O O

A Free Body Diagram of a Soil Slice for the
Three Dimensional Model . . . . . . . .

A Free Body Diagram of an Inclined Blade . .

Three Dimensional Representation of an
Inclined Plane Tillage Tool . . . . . .

Three Dimensional Representation of an
Inclined Plane Tillage Tool . . . . . .

Soil Path on the Surface of the Inclined
Plane Tillage Tool. . . . . . . . . . .

Soil and Tillage Tool Orientation in the
Three Dimensional Space . . . . . . . .

Hermite Curve Representation . . . . . . . .

Hermite Curve With Varying Tangent Vector
Length 0 O O O O O O O O O O O O O O O O

Bezier Curve Representation. . . . . . . . .

Sixteen Control Points for Bezier
Bicubic Patch . . . . . . . . . . . . .

Location of Four Elements on Moldboard . . .
PLOW menu 0 O O O O O O O O O O O O 0 O 0 O 0

Example of PLOW Force Output . . . . . . . .

vi

Page

11
11
1H

1H

15
20

20

22

22

25
28

28
30

36
A0
A2

A2

6.1a 63 mm Wide Inclined Blade. . . . . . . . . .
6.1b 12h mm Wide Inclined Blade . . . . . . . . .
6.1c 200 mm Wide Inclined Blade . . . . . . . . .
6.2 Sixteen Inch Ford Plow Descriptor Points . .
6.3 Representation of a Moldboard Plow Bottom by

Multiple Three Dimensional Inclined
Planes. . . . . . . . . . . . . . . . .

vii

(11>

O: Q D O 0’

LIST OF SYMBOLS

 

Adhesion of $011 on Metal
Acceleration Force of Soil
Tool Width

Soil Cohesion

Direction Cosine

Depth of Tool Cut

Distance the Soil Travels Along the Forward Failure
Shear Plane

Distance the Soil Travels on the Tool Over the
Described Path

Draft Force

Lateral Force

Vertical Force

Contact Area of Soil Block and Inclined Blade
Area of Forward Shear Failure Surface

Weight of Soil Segment

Acceleration Due to Gravity

Distance From One Shear Plane to Next

Mass of Accelerated 8011

Normal Load on Inclined Blade

Normal Load on Leading Forward Failure Surface
Normal Load on Trailing Forward Failure Surface
Time

Absolute Soil Velocity

viii

_ oo~<ufl

Q t 'R d

0-.

¢!

Velocity of Soil Along Inclined Blade
Velocity of Soil Along Shear Plane
Forward Velocity of the Inclined Blade
Initial Soil Velocity

Velocity of Accelerated Soil

Velocity of the Tool

Angle Described in Figure (3.7)

Angle of Forward Failure Surface

Bulk Density of Soil

Lift Angle of Tool

Angle Described in Figure (3.9)

Soil Internal Friction Coefficient
Soil-metal Friction Coefficient
Normal Stress

Shear Strength

Angle of Internal Friction

Angle Described in Figure (3.8)

ix

Chapter 1
INTRODUCTION

 

Ancient Egypt, Greece and Rome tilled their fields with
a "Y" shaped stick pulled by draft animals or slaves. This
plow design was practically unchanged until the 18th Century
(Carlson, 1961) when formed metals were used for the
moldboard (Gill and VandenBerg, 1967). Over many years
moldboard plow shapes have evolved through the 'trial and
error' design efforts of many individuals. The modern
moldboards are very complex surfaces that have required
elaborate and restrictive shape description techniques to

manufacture.

1.1 Computer Graphics

 

The aeronautics industry has developed graphic
techniques to describe complicated shapes of aerodynamically
efficient components (Foley and VanDam, 1982). These
techniques are directly applicable to the complex surfaces
of moldboard plows. They are very general and can represent
virtually any surface. The technique used in this study is
a Bezier bicubic parametric surface.

Graphic techniques by their nature also offer a very
effective commUnication media between the computer and the
Operator. The ancient Chinese proverb "one picture is worth
a thousand words" is also very applicable when communicating

with the computer. Graphic displays are more natural and

1

2
efficient than the 'traditional' alphanumeric displays used

on most computers (Foley and VanDam, 1982).

1.2 Soil-Tool Mechanics

 

Soil-tool mechanics offers a method of describing the
application of forces to the soil and how the soil reacts.
It was pioneered in the 1930's by M. L. Nichols and applied
to a simple two dimensional inclined blade in 1956 by W.
Soehne (Gill and VandenBerg, 1967). Other researchers have
developed mechanics that are dependent on emperically
derived soil paths (Emmet, 1983; O'Callaghan and McCoy,
1965).

This study modifies Soehne's two dimensional inclined
blade mechanics and expands it to three dimensions for
application to a moldboard plow. The soil path is
calculated based on the plow surface shape and forces

exerted on the soil.

1.3 Research Objectives

 

Research reported here sought to combine computer
graphics and three dimensional soil-tool mechanics to
predict the force required to move a moldboard plow through
the soil. The specific objectives were:

1. To develop a computer model that:

A. Uses three dimensional graphics techniques to
describe a share and moldboard surface.

B. Evaluates three dimensional force and moment
reactions of a single plow without prior
knowledge of the soil path over the plow
surface.

3

C. Allows the user to specify a plow shape, size
and set of operating conditions.

2. To verify the model results using existing force
data for a moldboard plow.

Chapter 2

LITERATURE REVIEW

 

2.1 Plow Description

 

One of the first moldboard description techniques for
construction purposes was developed by Thomas Jefferson
(Gill and VandenBerg, 1967). The rigid framework of Figure
(2.1) was its' basis. Using lines e-m and o-h as guides, a
straight edge is moved from e-o to m-h so that the straight
edge remains in a plane parallel to the y-z plane. The path
swept out by this line describes the surface. I

Several other moldboard shapes have been constructed
with a modified Jefferson framework. The guides and
straight edge were changed from straight lines to various
catenaries, arcs of circles, etc. Intuition was always the
basis for choosing one of these shapes over another (Gill
and VandenBerg, 1967).

White (1918) was the first to mathematically describe
the surface generated by this method. He achieved this by
moving a reference cartesian coordinate system around on the
surface, reducing the description to a series of hyperbolic
paraboloids (Gill and VandenBerg, 1967).

Several graphical techniques have been used by
researchers to describe the shape of the surface. The most

common method is to project the shape onto a gridded plane

4

 

 

 

 

Figure 2.1. Thomas Jefferson's Framework for Designing a
Moldboard Plow.

6

perpendicular to the direction of travel and then on to a
plane parallel to the direction of travel. For example, if
the relative position of a point in the perpendicular planes
is known, a contour line for constant x value can be
determined. Therefore, by plotting the contour lines in any
one of the coordinate directions and connecting the contour
line ends, a two dimensional representation of the surface
is drawn (Gill and VandenBerg, 1967).

Soehne (1959) photographed a slit of light projected on
the plow surface. By keeping the camera and light projector
stationary and moving the surface incrementally, the same
coordinate dimensions as described in the previous paragraph
could be found.

Although moldboard shapes had to be defined for
fabrication purposes, they were also defined in hopes of
relating plow shape to forces applied by the soil or the
resulting soil conditions. In this research, the actual
path the soil travels over the moldboard was determined.
Gill and VandenBerg (1967) reported Pfost, Corley and Soehne
in separate efforts, covered the moldboard surface with a
thin coat of lacquer. The scratches left in the lacquer
after plowing a short distance were recorded using
techniques previously described (Gill and VandenBerg, 1967).

Nichols and Kummer (1932) developed a method of
mathematically describing the surface using polar
coordinates and arcs independent of soil path. Their method

allowed them to look at three different functional areas and

7
evaluate the soil-tool mechanics for each area.

With the advent of the digital computer, several
researchers (Carlson, 1961; O'Callaghan and McCoy, 1965;
Emmett, 1983; Orlandea, et al., 1983) fit multi-order
polynomials to these scratch traces and the contour lines.
With the development of hardware that allows the computer to
read the x, y and z coordinate values directly (Emmett,
1983). a very accurate description of the surface and/or

soil path has been implemented.

2.2 Soil—Tool Mechanics

 

M. L. Nichols was the pioneer in soil tillage
mechanics. He worked with various researchers in studying
the dynamic properties of soil and their influence on
tillage operations. He attempted to explain soil forces
using colloid film theory (Nichols, 1931a). Nichols (1931b)
also researched the evaluation of soil-metal friction
coefficients.

Nichols and Kummer (1932) deduced there are three
different frictional areas of a moldboard plow. The lower
portion, forming a wedge for breaking the soil loose, a
central pulverization area and a turning or inversion area
on the upper portion. Using this theory they defined and
designed each of these areas to manipulate the soil as
desired.

Nichols and Reed (193A) found that two shear planes are
formed as a moldboard plow is pulled through the soil. The

primary shear plane is A5 degrees to the right of the

8
direction of travel and the secondary plane is 90 degrees
from the primary plane which is parallel to the leading edge
of the share.

Until the 1950's, the evaluation of moldboard plow
performance was largely based on visual observations.
Nichols, Reed and Reaves (1958) evaluated plow-share design
by observing the soil action in a glass sided soil bin and
then measuring the draft forces. The effects of shape,
wear, angle of approach and soil conditions were considered.

In 1956, Soehne analyzed the action of a simple
inclined plane (Gill and VandenBerg, 1967). Four behavior
equations described the tillage action. They are soil-metal
friction, shear failure, acceleration of a soil block after
forming shear planes and cutting resistance. Rowe and
Barnes (1961) added an adhesion term to soil-metal friction.

The inclined plane moving through the soil was assumed
to be in equilibrium. Using a force balance, the draft and
vertical loads could be computed (Rowe and Barnes, 1961;
G111 and VandenBerg, 1967). Several variations of this
analysis have been studied.

Rowe and Barnes (1961) concluded that the rate of
shearing of the soil was more important than acceleration
and should be included in the mathematical model.
Hettiaratchi and Reece (1967), Kawamura (Gill and
VandenBerg, 1967), McKyes (1977), Perumpral, et a1. (1983),
and Siemens, et a1. (1965) defined the shear planes

differently and completed analysis accordingly.

9
Carlson (1961) fit multi-order polynomial equations to
soil path scratches. He then differentiated the soil
displacements twice to get accelerations. Emmett and Gunkel
(1983) and O'Callaghan and McCoy (1965) also implemented
this technique to compute accelerations for draft prediction

models.

Chapter 3

SOIL-TOOL MECHANICS MODEL

 

3.1 Two Dimensional Model

 

W. Soehne developed the soil-tool mechanics model for
an inclined blade using four simple behavior equations.
They are soil-metal friction, shear failure, acceleration
force of each soil block and soil cutting resistance (Gill
and VandenBerg, 1967). Rowe and Barnes (1961) added an
adhesion term to the soil-metal interface stress equation.
The following assumptions have been made to realize the
mechanics (Rowe and Barnes, 1961).

l. The soil failure ahead of an implement is by a

series of shear failures.

2. The angle of the forward shear failure plane is
dependent only on the soil internal frictional
properties.

3. The soil internal shearing strength and the
soil-metal shearing strength can be approximated by
a linear function of the stress normal to the
surface of shear or of sliding.

A. The soil is considered to be a homogeneous and

isotropic material.

The following development of the two dimensional
soil-tool mechanics is from Gill and VandenBerg (1967).
Figure (3.1) shows the forces that act on a soil block as
the blade advances. Forces CF1 and‘uNi are soil cohesion
and internal frictional forces. They are due to soil shear.

10

11

 

 

 

Figure 3.1. Forces on a Soil Segment Reacting on an
Inclined Plane Tillage Tool.

slope=u 1],,

\
\

' 7‘
C } \: a

 

Figure 3.2. A Mohr Envelope of 8011 Internal Stresses.

12
CF1 is present only at the instant the shear failure plane
forms. Forces u'No and AFO are soil-metal friction and
adhesion forces. Force B is due to the soil acceleration
and G is the soil weight. The pure cutting resistance force
occurs at the edge of the tool. Soehne determined that the
pure cutting resistance of the soil is very small and only
important when the tool is dull or when there are stones or
organic matter present. Therefore, it is not included in

this analysis.

The shear strength of the soil can be expressed as

t = C +Au a (3-1)
where

t = soil internal shear strength

C = apparent soil cohesion

u = soil internal friction coefficient

a = normal stress on shear surface

Mohr, according to Gill and VandenBerg (1967). proposed
the following shear failure theory. Shear strength and
normal stress have a functional relationship on the failure
plane. If a series of different stress states are applied
to the same material and plotted as Mohr's circles, a line
drawn tangent to all the circles defines a shear failure
envelope (Figure (3.2)). The soil internal friction

coefficient can be expressed as

p.= tan¢ (3.2)

where ¢ is the angle of soil internal friction. The shear

13
stress at zero normal stress is the soil cohesion.
The parameter used to geometrically describe the soil
internal friction is beta. Figure (3.1) shows the angle

beta. Numerically it is

1

s = 1/2 (90 - tan- u) (3.3)

Like the shear strength of soil, the stress at the

soil-metal interface can be expressed as

t = A + u'0 (3.A)
where

t = shear stress at soil-metal interface

A = adhesion of soil on metal

u' = soil-metal friction coefficient

a = normal stress on soil-metal interface

Again plotting several stress states the failure envelope
can be drawn (Figure (3.3)). The adhesion and soil-metal
friction are the y-intercept and slope, respectively.

The reaction of soil to an inclined plane is assumed to
be repeating failure by shear forming small blocks of soil
(Figure (3.A)). As one soil block moves onto the incline
another soil block moves off. The sum of the forces the
soil blocks exert on the incline will be the force required
to propel the blade through the soil.

Using the notation in Figure (3.5), an equation can be

written for a soil block as follows

 

 

 

   
     

 

 

 

   
 

  

~\\\\\\\ 0
Figure 3.3. A Mohr Envelope of Soil-Metal Stresses.
forward
«— shear planes
h
i \
d v \
0 \
& —--
'/rr/lf//l‘/ fI/Ff’ WW" r//T' //f/:/’//:,7 rv'//7/
//// ////////// / /t’/ // /,{/: /////,///////l /////I// //:II /////I/
//./" ///// I // /// I ”/ I/I'l/ [VT/x. // A / I l// l ' //
inclined
blade

Figure 3.A. Formation of Soil Blocks by Repeated Shear
Failure.

15

mouse . . onzmfim
choamcofiwamonm < m m
.waouwswofio an
How a
ooaam H
on» Low

 

 

 

n-cm

 

N:
ocean cocadocs

 

soodm 30%
t

16

(B + CF1 +ILN1 -[LN2) cos(5 - p) +

(N2 - N1) sin(8 - 5) + G sin 8 - N0 = O (3.5)

where
B = acceleration force of soil
C = cohesion of soil
F1 = area of forward shear failure surface
# = coefficient of internal soil frictions
N1 = normal load on leading forward failure

surface

N2 = normal load on trailing forward failure
surface

G = weight of soil segment

N0 = normal load on inclined blade
5 = lift angle of tool
fl = angle of forward failure surface

Upon close inspection of the free-body diagram, it can be
shown by summing forces perpendicular to the forward failure
surface that the difference between N1 and N2 is a component

of the soil block weight G. As the soil blocks are reduced

in size N1 will equal N2, thus reducing Equation (3.5) to
N0 = (B + CFl) cos (5 - B) + G sinb (3.6)

The weight of soil may be calculated from the volume of

soil between shear planes. From Figure (3.”)

G __, thd (3.7)
sin.fi

17

where

= bulk density of soil

= tool width

distance from one shear plane to next
= depth of tool out

= angle of forward failure surface

In oatr o‘~<
u

The area of shear F1 can be computed from Figure (3.u)

F = bd (3.8)
sinli

 

The acceleration force B can be derived from Newton's

Second Law of Motion

B 3 m g! (3.9)
dt
where
m = mass of accelerated soil

velocity of accelerated soil
time

('1‘
II

If the periodic failure process is considered to be
continuous, the total work done for acceleration to occur
does not change significantly. Therefore, an average
constant force can be determined. The mass to be
accelerated in time t is assumed to be the volume of soil

disturbed in that time.

bdtV (3.10)

m= O

on|~¢

18

where
Y = wet bulk density
b = width of tool
d = depth of tool
t = time
g = acceleration due to gravity
Vo = velocity of the tool

The soil can be visualized to be at rest and in time t

it reaches a velocity Vs’ as shown in Figure (3.").

Therefore,

__.= _§____ =.J§ (3.11)

The magnitude of the velocity vectors are related so

that they form a closed triangle. From Figure (3.fl)

V0 = Vs cos B + Ve cos 8 (3-12)

and

VS sinfi = Ve sin 8 (3.13)

Eliminating Ve gives

sin!)
V - V (3.14)

8 ° sin(6+fl)

Substituting Equations (3.10), (3.11) and (3.14) into

Equation (3.9) and simplifying gives

19

2 sin 8
O (3.15)

8 sin(8 + a)

C and B are soil parameters determined by soil type. All
the variables of Equation (3.6) are now known, therefore, NO
can be computed.

With NO known, a free body diagram of the inclined
blade is considered. From Figure (3.6), forces in the x and

z direction can be summed.

’13
ll

(AF0 + uFNO) cos 6 (3.15)
F = (AF0 + u'NO) sin 8 + No cos 8 (3.17)

where

= draft force

F
Fz = vertical force
A = soil-metal adhesion
F

O = contact area of soil block and inclined blade
n' = coefficient of soil-metal friction

Forces computed in Equations (3.16) and (3.17) are due
to individual soil blocks. Total inclined blade forces are
summations over the number of soil blocks maintaining
contact. The thickness of the block, h in Figure (3.”), can
be selected for the desired resolution of the model.

The first soil block on the blade is different from the
others. The soil cohesion term is due to the force required

to form the forward shear plane. Once the shear plane is

20

X‘—

 

 

 

Figure 3.6. A Free Body Diagram of an Inclined Blade.

90-5
Xi

 

Figure 3.7. Three Dimensional Representation of an
Inclined Plane Tillage Tool.

21

formed, only the soil internal frictional force acts along

it.

3.2 Three Dimensional Model

 

The three dimensional model takes into account the
rotation of the inclined blade about z-axis, so that as the
blade moves through the soil, there is a side force parallel
to the y-axis. Figure (3.7) illustrates this orientation of
the blade as it moves in the negative x-direction.

Two angles are used to describe the blade orientation
in three space. Delta 8 is the angle between 0A and the
z-axis in Figure (3.7). CA is a line in the plane of the
inclined blade and in the x-z plane. Alpha a is the angle
between OB and the x-axis. OB is a line in the plane of the
inclined blade and in the x-y plane.

To aide in simplifying the force calculations, an angle
phi prime ¢' is defined in terms of alpha a and delta 8.

The relationship is

o' = tan-1 (sina tan5) (3.18)

Phi prime ¢' is the angle between CD and CE in Figure (3.8).
CD is a line in the plane of the inclined blade and
perpendicular to line OB. CE is parallel to the z-axis.

The forces considered in Equation (3.6) are also
considered here. The forward failure shear surface is
assumed to form perpendicular to line 0B of Figure (3.8) and
at an angle beta.p from the x-y plane. The acceleration

and cohesion force act along this plane.

22

 

 

 

Figure 3.8. Three Dimensional Representation of an
Inclined Plane Tillage Tool.

3 Path of 5011

f’ .. +

Inclined Blade

 

 

 

Y“!

 

0

Figure 3.9. Soil Path on the Surface of the Inclined
Plane Tillage Tool.

23
As in the two dimensional case, the soil weight may be

calculated from the volume between shear planes

G =Yh b d (3.19)
sina sinp

The area of the forward shear plane F1 is

F = ' bd (3.20)
sina sinB

 

Acceleration forces are computed from Newton's Second

Law of Motion

The mass of soil accelerated is

m = E (3.21)
S

For computation of the change in velocity with time,
the displacement of the soil must be considered. The
absolute displacement of the soil is along the forward shear
plane and in a direction perpendicular to OB of Figure
(3.7). The displacement of the soil relative to the incline
can be determined empirically with adequate facilities or
theoretically if all influencing forces can be quantified.
A direction of angle eta n from x-axis and in the plane of
the inclined blade was assumed (Figure (3.9)). Further, it

was assumed that forward velocity of tool is the velocity at

2“
which the soil traverses the blade over the above described

path. The time required to traverse the blade is

Vt
where
t = time to traverse blade

p distance the soil travels on the tool over
the described path

forward velocity of the inclined blade

Vt

The time required for the soil to traverse the path
over the blade (A to B of Figure (3.9)) is the same time
required for the soil to displace along the forward failure
surface (A to B Figure (3.10)). Therefore, the absolute

velocity is

v 3 .2 I (3.23)

where

<
II

A absolute soil velocity

distance the soil travels along the forward
failure shear plane

Q.
on
II

t = time to travel ds

The acceleration is

dv _ (Va - V

dt At

o) (3.21:)

 

25

    

z
\
\
\
\
\
‘\
‘\
‘\
‘\
\ 3
Soil Block \\
Inclined Blade
\
\
\
‘\
\
‘\
‘\
\

 

 

Figure 3.10. Soil and Tillage Tool Orientation in the
Three Dimensional Space.

26

where

V0 = initial soil velocity
At = time for velocity to go from V0 to VA

For this case the soil block starts at rest. Substituting
Equations (3.20) and (3.23) into Equation (3.8) gives

B = __e_ (3.25)

All the variables are now known in Equation (3.6),
therefore N0 can be solved for. As in the two dimensional
case, the free body diagram of the inclined blade can be
drawn. Note the direction of adhesion and soil-metal
frictional forces are in the opposite direction of soil

movement. Summation of forces in the x, y and z-directions

yields.
Fx = -(AFO + u'NO) cos 1| cos 5 + N0 cos ¢' sina(3.26)
Fy = -(AF0 + uJNo) sin'n + N0 cos ¢' cos a (3.27)
Fz = -(AFO + #‘No) cost] sin 8 + No sin ¢' (3.28)

CHAPTER u

GRAPHICS REPRESENTATION

 

The plow surface is modeled using Bezier parametric
bicubic surfaces. A single surface is derived by
representing a three dimensional curve and then generalizing
the mathematical development to three dimensional curved

surfaces. This development is from Foley and VanDam (1982).

h.1 Three Dimensional Curves

 

There“ are two ways to represent curves; as functions
of x, y and z or as functions of a parameter t, which is a
position along a curve. In the first case, it is difficult
to represent an infinite slope or a progression of points on
a segment if the segment includes a 100p. The latter case
is a parametric curve that allows closed and multi-valued
functions and replaces the slopes with tangent vectors.
This eliminates problems with infinite slopes.

Cubic curves are used because they are the lowest order
representation of curved segments that provide slope and
position continuity at segment endpoints. Higher order
representations also tend to have undesirable oscillations.
Lower order representations restrict the shape of the curve.

There are many ways to define a cubic parametric curve.
Only two are considered here. The Hermite curve

representation is shown in Figure (".1). It uses the

27

28

 

VJ SXY4
(X1.Y1) \\
SXYl \ny1
Figure “.1. Hermite Curve Representation.
increasing
constant

 

 

 

Figure h.2. Hermite Curve With varying Tangent Vector
Length.

29
positions of the curved segment endpoints (X1,Y1) and
(Xu,Yu); unit tangent vectors from each endpoint in the
direction of the other endpoint Sxyl and Sxyfl and tangent
vector lengths L and L

xyl
describe the 'fullness' of the curve segment. Figure (h.2)

xyu' The tangent vector lengths

shows how the curve changes as Lx u increases with Lx 1

y Y

constant.

The Bezier curve representation is shown in Figure
(B.3). It uses the endpoint positions (X1,Y1) and (Xu,Yu)
and the position of two control points (X2,Y2) and (X3,Y3)
which, when connected with the endpoints, define both the
unit tangent vectors and the tangent vector lengths. The
latter method makes the physical meaning of the tangent
vectors more apparent. This is important since this model
is intended for the design of plow shapes.

The previous description of the Hermite and Bezier
curve representations used only two dimensions for clarity.
This can be expanded to three dimensions by adding a z
coordinate to each point description in the Hermite and
Bezier forms and by adding two tangent vectors 3 and S

yz zx’

and two tangent vector lengths L and sz, to each endpoint

yz
in the Hermite form.

The derivation of the equations for calculating curve
point positions follows. The Hermite representation is used
in deriving the Bezier form, therefore it is considered

first. The derivation uses the subscript x to refer to the

x component. The y and z derivations are similar and are

30

 

(x2.Y2)

Figure “.3. Bezier Curve Representation.

31
not shown here.

A cubic curve parameterized with the variable t is

_ 3 2
x(t) — axt + bxt + cxt + dx (U.1)

In matrix form

[t3 t2 t 1]

a
x(t) = b (h.2)
c
d
or
x(t) = T Cx (“-3)
differentiating Equation (fl.1) with respect to t
Si—"—=3at2+2btz+c “H”
dt x x x
in matrix form
[:31:2 2t 1 O:l a
x'(t) = b (n.5)
c
d
Subject to Hermite form conditions
x(0) = P1x
X(1) 3 Pux (“06)
x'(O) = R1x
' =
x (1) th
P P

lx’ “x and Rlx’ Rux are endpoints and tangent vector

32
endpoints, respectively. Index points 1 and H are used
instead of 1 and 2 for consistency with the notation used
for Bezier curves. Substituting these conditions into

Equations (fl.2) and (h.5)

x(0) = P1x = [b o o {] cx
x(1) = Pux a [1 1 1 {] cx (u.7)
x'(O) = Rlx = [o o 1 Q] cx
x'(l) = Rux = [3 2 1 d] cx
OP
P1 0 o o 1
Pa = 1 1 1 1 cx (u.8)
R1 0 o 1 o
Ru 3 2 1 0

Solving for Cx by inverting the “x“ matrix in Equation (u.8)

yields
2 -2 1 1 P1
ox = -3 3 -2 -1 Pa (u.9)
o o 1 0 R1
1 o o o Ru
OP
Cx = Mn th

where Mh is the Hermite matrix and th is the Hermite
geometry vector. Substituting for Cx in Equation (fl.3)

gives

T M G (".10)

x(t) hx

Similarly

33

ll
0-3
3
Q

y(t) (“.11)

z(t) (“.12)

II
a
3
Q

As discussed previously, Bezier curve endpoint tangent
vectors are defined by two points. These points are defined

by the following relationships.

R

P'(0) (“.13)

1 3(P2 ' P1)

P'(1) (“.1“)

Ba

The relationship between the Hermite geometry matrix Gh

and the Bezier geometry matrix Gb is

P1 1 0 0 0 Pl ° “ )
G = P = 0 0 0 1 P = M G ( .15
h a? -3 3 o 0 P: “D b
Substituting for Gh in Equation (“.10) gives
x(t) = T Mh Mhb Gbx (“.16)
If we let
Mb = Mh Mhb (“.17)
then
x(t) = T Mb Gbx (“.18)

Mb is called the Bezier matrix. Numerically it is

-l 3 '3 1

Mb = 3 -6 3 0 (“.19)
-3 3 O O
1 O O O

“.2 Three Dimensional Surfaces

 

To move from the cubic curve to the bicubic surface an
additional variable 3 is used. Replacing t with s in
Equation (“.18) we obtain

x(3) = S M G

b (“.20)

bx

Instead of writing the Bezier geometry matrix as a

constant, it can be written as a function of t.

x(s,t) = 3 Mb Gbx(t) = S M

(t)

(t) (“.21)
(t)

(t)

Now let P1(t), P2(t), P3(t), Pu(t) each be cubes in Bezier

form, so that

P (t) p11 (a 22)
1x b p12 '

p13
pl“

ll
*3
3

p21
T Mb p22 (“.23)

P23
qu

P2x(t)

P (t) = T M p31 (u 2n)
3x b p32 °

p33
p31;

35

Pal
Pux(t) = T Mb pu2 (“.25)

91,3
Pun

As a row vector
P1(t) P2(t) P3(t) Pu(t) =

p11 p21 p31 p“1
b p12 p22 p32 Pu2
p13 p23 p33 ”“3
Plu Pan P3u Pun

(“.26)

Transposing both sides of the equation with the identity

(ABC)T=CT BT AT

P (t) p p p p
P1(t) 3 p11 p12 p13 p1u M T TT = P M T TT (u 27)
2 21 22 23 2a b x b
fi3‘t) p31 p32 p33 P3“
u‘t) Pu1 Pua Pu3 Pun
Substituting into Equation (“.21)
_ T T
x(s,t) - 3 Mb Px Mb T (n.28)
Similarily
- T T
y(s,t) - S Mb Py Mb T (“.29)
_ T T
z(s,t) - 3 Mb Pz Mb T (n.30)

The geometry matrix P consists of sixteen control

points as shown in Figure (“.“).

36

.nopwa cansoam
Loaumm com manaom Honucoo :oouxam

.=.= onsmfim

 

Chapter 5
COMPUTER MODEL

 

The moldboard plow model and graphics representation
were implemented on a Prime 750 computer equipped with a
Fortran V compiler, a Tektronix Plot-10 Terminal Control
System (T08) and a Tektronix “000 Series terminal. Also
used in this study was three dimensional plotting software
written by personnel of The Case Center for Computer Aided
Design, Michigan State University.

The moldboard plow was the primary application of this
program so plow terminology is used. However, the reader
should be aware that the program is also applicable to other
fixed surface tillage tools. These applications are

discussed in Chapter 6.

5.1 Moldboard Plow Model

 

Two Bezier surfaces are used to graphically depict a
moldboard plow. One for the share and one for the
moldboard. Each surface is divided into planar elements
with four corner points describing their bounds. Each
element has the forces computed in the x, y and 2 directions
as if each were an independent plane. The sum of forces in
each direction over the entire surface is the force required
to propel the plow through the soil. Moments created by the
forces on each element are also summed about a user

37

specified point.

location

The

38

of the plow to the implement frame.

following geometric assumptions were used

implement the moldboard plow model:

1.

The

important

The point of the share has (x,y,z) coordinates
of (0,0,0).

The plow moves in the negative x-direction.
Therefore, the control points used to describe
the surface shapes must also describe the
orientation of the plow as it moves through
the soil.

Each surface is divided into one to one
hundred elements; up to ten divisions in each
of the s and t directions. See Figure (5.1)
for the definition of the s and t directions.

The same number of elements in the s and t
directions are used. However, the share and
moldboard can have a different number of
elements.

All four edges of the share are straight.
Therefore, only the corner point coordinates
are used to describe the surface.

The edge of the moldboard adjacent to the
share is straight.

Elements in the curved portion of the surface
do not have all four corner points in a plane,
since only three points describe a plane.
Therefore, the point farthest from the plow
point on each element is not used to compute
the plane orientation.

The four center surface control points are not
user specified. They are computed from the
specified edge points such that the surface
has twist vectors of zero. This eliminates
local surface irregularities.

This point may correspond to the mounting

to

soil path over a share and moldboard is very

in determining the plow forces. This study

assumes the soil path is a function of force exerted on soil

39

by the element geometrically preceding it on the Bezier
surface. For example, in Figure (5.1) the soil path over
element two is dependent on the forces applied to the soil
by element one. The soil path over element three is
dependent on the forces applied to the soil by element two.
The soil path over element four is dependent on the forces
applied to the soil by element three, and so on.

The soil path dependency is determined by the y and z
force direction cosines and the geometric orientation of an

element (that is, the soil must be in contact with the

 

 

 

 

element).
F
Dy = y (5.1)
2 2 2
\/(Fx + 18y + F2 )
F
Dz = Z (5.2)
2 2 2
fix + Fy + F2 )
D
D = -—~Y— + D tan 8 (5’3)
x tan¢1 z
where

D = the direction cosine

F 8 force applied to the soil by the previous
element

a = as defined in Chapter 3
8 = as defined in Chapter 3

The x direction cosine is not used because this adds a

fourth constraint to a three dimensional system. Relaxing

“1
the x direction cosine constraint for a moldboard plow is
reasonable since the absolute soil movement is up and to the
right (2 and y directions) with very little absolute soil
movement in the direction of plow travel (x direction).

As stated previously, the soil path is assumed to be a
function of the force exerted on the soil by the previous
element. Before the soil moves on to the first row of share
elements, however, there has not been a force applied to the
soil. The path of the soil over the first row of elements
was assumed to be parallel to the direction of tool travel

(there is no soil movement in the y direction).

5.2 Computer Program

 

An interactive computer program named PLOW was
developed to prove that computer graphic techniques (i.e.,
Bezier surfaces) could be used to represent tillage tools
(primarily moldboard plow surfaces) and then be combined
with a soil-tool mechanics to calculate the forces exerted
on the tool by the soil. Several 'utilities' were written
to realize the implementation. To assure the user that the
plow is accurately represented, a picture of the plow is
plotted on the terminal using drawing functions. There are
also editing functions to change the data files. Lastly,
there is a function to calculate the forces on a tillage
tool.

Figure (5.2) is the menu displayed by PLOW. Procedures
1 thru 9 are drawing functions. Procedures 10 thru 13 are

editing functions. Procedure 1“ is the force calculation

“2

1) Plot

2) Change point of view

3) Change look at point

“) Change magnification factor

5) Change center position

6) Set plotting conditions to default
7) Label control points on plot

8) Plot measured points

9) Change number of surface elements
10) Change operating conditions
11) Change share, moldboard or scaling
12) Get new tool geometry data
13) Get new operating data
1“) Compute forces

15) End
Enter number of desired procedure (99 for menu).

Figure 5.2. PLOW Menu.

Tool totals X Y Z
Forces (N)
Share 2357“. -11936. -18053.
Moldboard 3988. -13“9. -951.
Total 27562. -13285. -1900“.

Moments (N-M)

Moldboard 253. 2613. -“158.
Total -8813. 15260. -2“358.
Location of zero moment 0.00 0.00 0.00

Figure 5.3. Example of PLOW Force Output.

“3
function and procedure 15 is the program termination
function.

PLOW is setup to run with or without a graphics
terminal. The graphics are Tektronix “000 series terminal
compatable. Procedure 1 is for plotting the tillage tool
any time the user is prompted for a procedure. Procedures 2
thru 6 are purely for drawing pictures. Procedure 7 is for
plotting the control points discussed in Chapter “.
Procedure 8 is for plotting any points the user wishes.
Procedure 9 is used to change the number of elements on the
tillage tool. This is also considered an editing function
for force calculations. The drawing functions are not
usable without a graphics terminal.

The editing functions are for changing the information
read from two data files. This information is the control
points of the tool geometry and operating conditions. As
stated above, procedure 9 is for changing the number of
elements on the tillage tool. Procedure 10 is for changing
the soil and tool operating conditions. Procedure 11 is for
changing the control points of the share or moldboard.
Procedure 11 is also used to change the overall geometric
scaling factors. Procedure 12 is for reading a different
tool geometry file and procedure 13 is for reading a
different operating conditions file.

Procedure 1“ computes the forces on the tillage tool.
Force calculations are made on each element and then summed

over all the elements of the tool. In addition, the user is

““

prompted for a point that has a zero moment. The distance
from each element centroid to the zero moment point is
calculated and multiplied by each of the component forces to
give a torque. The element torques are summed to give tool
torques about each of the component axes. The zero moment
point might correspond to the mounting location of the tool
to the implement frame. Figure (5.3) shows a typical
output.

Finally, procedure 15 ends the program gracefully. The
users guide for PLOW installed on the Prime in the Case

Center at Michigan State University is in Appendix A.

CHAPTER 6

RESULTS AND VERIFICATION

 

Verification of the computer program is accomplished by
comparing the predicted forces with actual force
measurements. The measured force data was extracted from
publications of authors with the equipment and resources for
the measurements. Two types of tools are considered for
comparisons. First, inclined blades of varying widths in
four soil conditions are discussed and second, a moldboard

plow in one soil type is discussed.

6.1 Inclined Blade

 

Plasse, Raghavan and McKyes (1983) reported that F. L.
Desir measured the draft force for narrow tines in varying
soil conditions. In the context of this discussion narrow
tines are inclined blades. Two soil types were used; clay
and sandy clay, each at two moisture contents. Three blade
widths were used; 63, 125 and 200 millimeters. The
operating depth and lift angle were held constant at 250
millimeters and 35 degrees, respectively.

Although the computer program was written for modeling
a conventional moldboard plow, the inclined blade can be
considered geometrically a very simple moldboard without
lateral or vertical curvature. For further simplification,

4S

“6

when a blade is considered rectangular, in a single plane,
and the leading edge is perpendicular to the direction of
travel; the soil will not apply a lateral force to the
blade. For this tool orientation the model is two
dimensional. The same results would be achieved with the
two dimensional model developed by Soehne discussed in
Chapter 3.

The description of the inclined blade for computer
representation was very simple. Like the moldboard plow,
the -blade was represented as two surfaces; the share and
moldboard. However, unlike the moldboard, the surfaces were
in the same plane with straight edges. The control points
on the moldboard edges were specified to be on the straight
lines connecting the corner points. By changing only the
y-direction scaling factor the three blades 'were
geometrically represented (Figures (6.1a, b and c)).

When the draft force of each blade was measured, the
operating conditions were also measured. These parameters
were soil cohesion, soil internal friction, soil bulk
density, soil-metal friction and soil-metal adhesion. Table
(6.1) shows the parameters and draft forces as they were
measured. The vertical and lateral force measurements were
not reported.

The right most column of Table (6.1) is the draft force
predicted by the model. Comparison of the predicted forces
to the measured forces (second from right of Table (6.1))

shows wide variation. This can be explained by considering

“7

 

Y

Figure 6.1a. 63 mm Wide Inclined Blade.

 

 

 

 

Figure 6.1c. 200 mm Wide Inclined Blade.

48

 

 

>m¢.m >>.m m.m o—m. mmm. oa m.wP oom NF
Fbm. m—.m o.m won. mmmP mm m.o— mm, FF
Foe. o>.m 0.? eon. mow— mm n.0— no 0—
emm. ho.fi m.m woe. mas, em o.m com o
.om._ mv.m v.m mom. «mm_ em m.>w mm, m
o_e. em._ m.m mes. men. mm m.m no N
mam.m mw.m m.o. eon. mvov em m.>m cow 0
mm_.m oa.o m.o owe. 0mm om m.mm mm? m
was. mm.m _.m_ mmm. . mmoi an m.m_ no a
nmm.¢ o>.o_ 0... son. mmo_ mm m.me 00m m
mF>._ oo.m w.m_ mmw. 500— me m.>_ mm. m
mmm.. m_.w m.. wee. emow mm m.mm no _
3.3
Azgv Azxv Amado is :0\mev mamg<
couch monom commosv< coflpoflpm hpwmcon :ofipOHum Asov
among pmmpm Hobos Hmpmz xasm HangoudH Ammxv new“: *
nopofivmpm cmuzmwmz uaflom aHHom HHom HHom :ofimmsoo madam ommo

.moouom Beneficopm was moouom umpsmwoe .mcofipflccoo mafipmpmmo woman omcflaodH —.w mgm<e

“9
the logic of the model.

Since the model was originally intended to predict the
forces exterted on a moldboard plow, those forces associated
with forming the shear planes, lifting and turning the soil,
and resistances to movement in the soil and at the soil-tool
interface were included. Three forces were excluded to
simplify the calculations.

The first is pure cutting resistance. Soehne
determined in his research (Gill and VandenBerg, 1967), that
the cutting resistance is only important when the tool is
dull or when there are stones or organic matter present.
Since these conditions were not present, excluding the pure
cutting resistance is appropriate for the inclined blade
also.

The other two forces excluded are due to the
interaction between the soil moving on the blade and the
stationary soil at each side. These forces are soil
cohesion and soil internal friction. A moldboard plow does
not have either of these forces on the dead furrow side.
The soil internal frictional force does not apply to the
land side since the soil's first movement is away from the
land. The soil cohesion force in the land side shear planes
was considered small compared to the cohesion in the forward
shear planes and the acceleration forces required to lift
and turn the furrow.

The latter two forces discussed may have a significant

effect on the total draft force of the inclined blade.

50
There is stationary soil on each side of the moving soil,
adding a soil internal frictional force. Also, as the width
of the blade is reduced, the area difference between the
cohesion forward and side shear surfaces becomes small,
giving the side cohesion a relative larger contribution. In
general, the relative contribution of smaller forces becomes
more significant when a planar tool is used because the
acceleration force, which is the largest force, becomes much

smaller.

6.2 Moldboard Plow

 

Emmet (1983) measured the forces on a moldboard plow at
the National Tillage Machinery Laboratory in Auburn,
Alabama. The test plow was a 16 inch Ford general purpose
bottom with a disposable-type share. The soil and operating
conditions are shown in Table (6.2). The soil-metal
adhesion is zero because an effective soil-metal friction
was measured for the plow material that included the
adhesion. Details of the test procedures and results were
reported by Emmet (1983).

The shape of the Bezier surfaces in the model had to be
manipulated to match the shape of the test plow. This was
accomplished by 1) measuring the x, y and z coordinates of
eight (8) descriptor points on the share and twenty nine
(29) descriptor points on the moldboard, 2) plotting them
over the Bezier surfaces on the computer terminal and 3)
changing the surfaces to match the descriptor points.

Figure (6.2) shows the descriptor points plotted on the

51

Table 6.2. Soil and Operating Conditions Used for Model

Verification.
Soil Type Norfolk Sandy Loam
Soil Bulk Density (kg/cu m) 1770
Soil Cohesion (kPa) 19.7
Soil Internal Friction Angle (deg) 3“.8
Soil-Metal Friction Coefficient 0.“
Soil-Metal Adhesion (kPa) 0.0
Forward Velocity (m/s) 1.3“

Tool Depth (m) .203

53

surfaces.

The plow on which the descriptor points were measured
had been used regularly for several years with little
maintenance, consequently the plow was worn considerably.
With this in mind, descriptor points on the share and
moldboard were selected to best describe the plow shape.
The plow was then mounted on the table of a Bridgeport
vertical spindle mill with the moldboard face up. Each
point was moved to a pointer chucked in spindle while the
table travel was recorded. Because of the limited travel on
the mill table, the plow had to be moved once to measure the
whole. plow. Datum points located before and after the move
were used to fit the two sets of descriptor points together.
The author estimated the accuracy of each measurement to be
plus or minus two (2) millimeters.

Once the descriptor points were modified to be in the
same coordinate system, they had to be manipulated so that
the moldboard sat upright as it would be during Operation.
This required writing a computer program to translate and
rotate the points by multiplying them by a transformation
matrix. Discussion of this technique can be found in Foley
and VanDam (1982) and Faux and Pratt (1979).

Matching the Bezier surfaces to the descriptor points
was a tedious trial and error process. Four of the
descriptor points were selected so that they fell on the
corners of the share surface and four were selected so they

fell on the corners of the moldboard surface. The control

5“
points were then manually manipulated to provide a visually
estimated minimum error between the descriptor points and
the Bezier surfaces. Additional interactive software could
be developed to greatly simplify this task. Software could
also be developed to mathematically quanitify and minimize
the error.

The number of elements on each Bezier surface can be
changed independent of the other surface. The greater the
number of elements, the smaller the size of each element and
the better the accuracy of geometric representation.
However, when the elements get too small and the soil 'sees'
a slightly curved surface section after a sharply curved
section, the soil deceleration force is much greater than
the soil gravitational force and the calculation of the
surface normal load becomes negative. Since this is
physically impossible, the normal load is set to zero.
However, the path of the soil over the next row of elements
is dependent on the forces on that row. Consequently, the
user is notified that the soil-tool model can not calculate
the forces. The one exception is when a zero normal load is
encountered on the last row of elements before the soil
leaves the plow.

To calculate the forces for a surface that has the
negative normal load problem, the number of surface elements
can be decreased until the element representation is course
enough to 'hide' the slightly curved surface section after

the sharply curved section. The maximum number of elements

55
for the Ford moldboard plow is 36 (that is, six along each

edge). Note that this model assumes an equal number of
elements along each edge for an individual surface.

As implied previously, when elements too large in size
are used, large geometric inaccuracies occur. Also, the
soil-tool mechanics eliminates forces with the assumption
that the elements are reduced in size. Therefore, the share
and moldboard surfaces with one or four elements were not
considered.

Figure (6.3) is a graphical representation of the share
with 16 elements and the moldboard with 25 elements. Force
results were generated for 9, 16, 25 and 36 element share
and moldboard combinations. The average of the forces in
each of the component directions is shown in Table (6.3).
For the x, y and z directions all the forces were within
3.5, 5.1 and 2.“ percent of the averages, respectively.

Table (6.3) also shows the comparisons between the
predicted forces and the measured forces. The predicted
draft force (Fx) is 15 percent larger than the measured
force; the predicted side force (Fy) is 53 percent larger
than that measured; and the predicted vertical force (Fz) is
115 percent larger than that measured.

The predictions could be improved by:

1. measuring the descriptor points more accurately.

2. fitting the surfaces to the descriptor points less
subjectively.

56

3. calculating the geometrical properties of each
element more accurately.

“. refining the soil-tool mechanics assumptions.

57

 

.mocsam vocaaocH Hocofimcosao Tosca sandpasz
an souuom zoam cpsonoHoz a no cofiuwucowopqom .m.o onsmam

 

 

 

 

58

Table 6.3. Comparison Between the Predicted Plow Forces
and the Measured Forces.

Mode Draft (n) Side (n) Vertical (n)
(Fx> (Ry) (F2)
Predicted 3895 -1712 ~3319

Measured 3389 -1118 -15“2

Chapter 7
CONCLUSIONS

 

Three dimensional graphics techniques (Bezier surfaces)
can be used to describe share and moldboard shapes.

The user should be aware of the assumptions of the model
when applying it to tool geometries other than a
moldboard plow. The computer model does not accurately
predict draft forces for inclined blades.

Three dimensional force reactions can be calculated
without prior knowledge of the soil path over the
moldboard plow surface. Moldboard plow draft force
predictions were within 15 percent of the draft force
measurements.

Force calculations for other moldboard plow shapes and

operating conditions should be checked against measured
forces to better determine the validity of the model.

R0

Chapter 8
SUGGESTIONS FOR FURTHER RESEARCH

 

This research illustrates only one application of
computer graphic techniques to modeling. Many others are
possible. These techniques are good for any application
where complicated and universal surface descriptions are
required. The soil tillage tool area has several more
applications. Modeling of shallow and deep narrow tines,
wide sweep or shovel cultivator tines and disk type tools
are all possible.

The computer program PLOW could be used for some of the
above applications. PLOW was written in modules to
facilitate replacing, adding or deleting functions as
required for a particular model. For example, a more
accurate inclined blade model could be developed by
replacing the force calculation subroutines with those more
closely representing that tool.

Another possible application of PLOW is using the
output as input to another program for designing an
implement hitching system. The three component dirction
forces and the three corresponding moments located along a
beam or tool bar can be summed to size the implement
structural members and hitch orientation to Optimize
available pulling power.

60

61

A very computer time intensive application would be
optimization of a tillage tool surface geometry to minimize
draft while maximizing the desired soil manipulation. There
would be a considerable amount of research required to
develop the optimization algorithms, but once in place only
computer time is needed to design a new implement.

There are many Computer Aided Design (CAD) software
packages available that use very elegant graphics
manipulations as part of their drawing functions. If these
were married with a force calculation function, the extent

of design using computers would be greatly expanded.

APPENDIX A

PLOW USERS GUIDE

PLOW USERS GUIDE

 

This document describes the operation of the computer
program PLOW running on a Prime Computer System. PLOW is
comprised of sixteen (16) procedures. This document
outlines the steps to execute the program and each procedure
within it.

The following conventions are used in this document:

1. Text or numbers enclosed in quotation marks " "
represents information displayed by PLOW.

2. Text or numbers enclosed in left and right carrots
< > represents information entered from the
keyboard by the user. Where more than one response
is possible, each response will be followed by an
explaination of that action.

3. Text or numbers not falling in areas 1 or 2 are
explainations.

To begin execution of the program enter:

<SEG #PLOW) SEG is the system command to execute a run
image. #PLOW is the run image. The rest
of PLOW is interactive with the user.

"Is a graphics terminal being used? (y or n)"
<y> is the normal response.

<n> will stop PLOW from sending graphics information
to the terminal. This severely limits the
usefulness of PLOW.

"Enter baud rate."
(1200) can be any of the standard baud rates. This

information is needed to initialize the three
dimensional graphics drivers.

A1

A2

"Enter file name for writing tool force information."

<D_FORCES> is any file name. The force calculations
will be stored here for later printing or
review.
"Enter name of geometry data file."
<D_FORD> is file name of geometry previously created
by PLOW.
"Enter name of operating conditions file."
<D_LMCLAY1> is file name of operating conditions file
previously created by PLOW.

The following main menu is displayed. The remainder of PLOW
will work from this menu.

"1) Plot

2) Change point of view

3) Change look at point

“) Change magnification factor

5) Change center position

6) Set plotting conditions to default
7) Label control points on plot

8) Plot measured points

9) Change number of surface elements
10) Change operating conditions

11) Change share, moldboard and scaling
12) Get new tool geometry data

13) Get new operating data

1“) Compute forces

15) End
Enter number of desired procedure. (99 for menu)"

<1> plots the active tool. Return to main menu
prompt.

<2) "Current view point is (-l.0,l.0,l.0). Enter new
X,Y,Z coordinates."

<1.0,.5,.2> The numbers identify the point

coordinates you are looking from.
Return to main menu prompt.

<3)

<“>

<5)

<6>

<7>

<8)

A3

"Current look at point is (0.0,0.0,0.0). Enter
new X,Y,Z coordinates."

<l.,0.,0.> The numbers identify the point
coordinates you are looking at.
Return to main menu prompt.

"Current magnification factor is 1.0. Enter new
factor."

<2.0> The screen is magnified about the screen
center by the value entered. Return to
main menu prompt.

"Current center point is (0.0,0.0). Enter new X,Y
coordinates."

<-1.,.5> The three dimensional image is projected
onto' a two dimensional plane, so the
value of the center point is relative to
the size and shape of the tool. The
user will have to 'get a feel' for this
procedure with each tool. Return to
main menu prompt.

"Plotting conditions set to default."

While the user is trying to get a better view of
a tool, he may modify several of the display
values. This procedure provides an easy way to
get back to the starting plot. Return to main
menu prompt.

This procedure labels the control points on the
tool. Return to main menu prompt.

"Enter measured points data file name."

<D_FORDPLOW> This file has data points that were
measured and entered into a file.
The file has four values per line.
The first value is an integer
identifying the point number and the
next three are the X,Y,Z
coordinates. There can be up to one
hundred points in the file. Return
to main menu prompt.

A“
<9> "Enter number of elements along the share edge."
<2) This number can range from 1 to 10.

"Enter number of elements along the moldboard
edge."

<5) This number can also range from 1 to 10. As
with the share, there is an equal number of
elements along each surface edge. Return to
main menu prompt.

<10) " Soil Parameters

1) Soil type: Norfolk Sandy Loam
2) Bulk Density (kg/cu m) 1770.

3) Cohesion (kPa) 19.7

“) Internal Friction Angle (deg) 3“.8

5) Soil-Metal Friction Coefficient 0.“0

6) Soil-Metal Adhesion (kPa) 0.0

Tool Parameters
7) Forward Velocity (m/s) 1.3“
8) Tool Depth (m) 0.203

Enter parameter number to be changed. (0 to stop)"

<2) "Enter new value for parameter number
2." '

(1690.)

"Enter parameter number to be changed. (0 to
stop)"

<0) "Store modified data? (y or n)"
<y> "Enter new file name."

<D_NEXT_SOIL> Return to main menu
prompt.

<n> Return to main menu prompt.

<11) "Do you want to change:
1) share
2) moldboard
3) scaling
“) stop"

A5
<1) Four share control points will be displayed.
"Enter point # to be edited. (0 to stop)"
(1) "Enter X,Y,Z."

<0.,0.,0.>

<0) Return to prompt at 11.

<2) Twelve moldboard control points will be
displayed.

"Enter point # to be edited. (0 to stop)"
<1) "Enter X,Y,Z."
<0.,0.,0.>
<0) Return to prompt at 11.
<3) "X-scale = l.
Y-scale 8 1.
Z-scale = 1.
Enter new X-scale, Y-scale, Z-scale."
<.5,.5,.5) This changes the size of the tool
in each of the component

directions by the factor entered.
Return to prompt at 11.

<“> "Do you want to store modified data in file?
(y or n)."

<y> "Enter name of tool geometry data file
to be written to."

<D_NEXT_GEOMETRY> Return to main menu
prompt.

<n> Return to main menu prompt.

<12) "Enter name of geometry data file to be read."

<D_NEXT_GEOMETRY) Return to main menu prompt.

<13) "Enter name of operating data file to be read."

<D_NEXT_OPERATING> Return to main menu prompt.

<1“)

(15>

A6

"Tool Totals X Y Z
Forces (N)
Moldboard 3988. -l3“9. -951.
Total 27562. -13285. -1900“.
Moments (N-m)
Share -9066. l26“7. -20200.
Moldboard 253. 2613. -“158.
Total -8813. 15260. -2“358.
Location of zero moment 0. O. ."
The numbers in the displayed table were just
calculated. The time taken to calculate depends
on the number of elements on the share and
moldboard and the load on the system. These
calculations can take several minutes. Return to
main menu prompt.
"End PLOW program."

APPENDIX B

PLOW SOURCE CODE

Bl

PLOW

This program interacts with the user to design a moldboard.
It was written to run on the Prime 750 computer at The
Case Center for Computer Aided Design, Michigan State
University. The computer was equipped with a Fortran-77
compiler and a Tektronix Plot-10 Terminal Control
System (TCS). In addition, three dimensional plotting
software developed by Case Center personnel was used.
This program has been run with a Tektronix 4100 series, a
Tektronix 4000 series, a DEC VTlOO with graphics expansion
card, and DEC VT240 terminals.

Every attempt was made to adhere to Fortran-77 conventions.
There are no guarantees expressed or implied.

C)C5C5C)C)C>C5C5C5C5CDC363C3C§C563

***********************************************************************
Written by:
Stephen B. Richey

January 1983

(2¢1(3(563(5C3C5

************************fi*******************************************fit

Array definitions:

PSS - coordinates of points blended in a direction describing
share surface.

PST - coordinates of points blended in t direction describing
share surface.

PMS - coordinates of points blended in 8 direction describing
moldboard surface.

PMT - coordinates of points blended in t direction describing

moldboard surface.

XS,YX,ZS - coordinates of share control points.
XM,YM,ZM - coordinates of moldboard control points.

DA - coordinates of measured points to be plotted.

********************

Variable definitions:

SOILTY - name of soil type (character string).
GAMMA - soil bulk density (kg/cu m).

COHES - soil cohesion (kpa).

PHI - soil internal friction (deg).

MUP - soil-metal friction.

ADHES - soil-metal adhesion (kpa).

VELTL - tool velocity (m/s).

DEPTHF - plowing depth (m).

C)CDC)(if)C5C5C3C3C5C3C563C3C3C)CDCDCDC5C3C3C3C3C5C)(if)C)

82

SIZEX,SIZEY,SIZEZ - scaling factors for controls points.
VX,VY,VZ - view point coordinates.

AX,AY,YZ - look at point coordinates.

MAG - magnification factor.

CX,CY - center position coordinates.

OLM,JJM - number of lines drawn on moldboard.
OLS,JJS - number of lines drawn on share.

FOFILE - name of data file for writing force info
(character string).

DATA - name of data file for reading measured data
(character string).

LUNTI - logical unit number for terminal input.
LUNTO - logical unit number for terminal output.
LUNDA,LUNDB,LUNDC - logical unit number for data input/output.

GRAPH - answer to question of graphics terminal
(character string).

0 - answer to question (character string).

01 - numerical answer to question.

BAUD - baud rate.

Subroutines called:

GETDAT - reads tool geometry data.
STODAT - stores tool geometry data.
OPERAT - reads soil and operating conditions.
DEFPLT - sets plotting parameters to default.

EDGE - calculates control points for a straight edge.
SURFACE - calculates point coordinates to be plotted.
SETPLT - sets up plot.

PLOT - plots surfaces.

LABPTS - labels control points on moldboard.

PARAM - edits soil and tool parameters.

NEWPS - edits share control points.

NEWPM - edits moldboard control points.

NEWSCA - edits scaling factors.

COMFOR - computes forces.

The following subroutines are on the Prime system at The
Case Center at Michigan State University.

INITT3 - initializes graphics file.
CLOSTK - closes graphics file.
MOVER3 - relative move (3-D).

MOVEA3 - absolute move (3-D).

DRAWR3 - relative draw (3-D).

RECOVR - recover to graphics mode.

HOME - move curser to screen upper left.
NEWPAG - page screen.

00000000000OOOOOOOOOCOOOGOOOOOOOOOOOOOOOOOOOOCOOOOOOOO

 

CDC)(3C)C)C3CJC)C)C5C)C5CDCJC3CJC3CDC)

¢~C3C3C3

B3
CARTN - generates carriage return.

CARTVP - sets cartisan view point.

LOOKAT - set lookat point.

PERSPT - turns on or off perspective viewing.

XUP,YUP,ZUP - defines which axis is up on plot.

PPCENT - sets distance to 2-D viewing plane from
look at point to view point as a
percentage.

DWINDO - sets bounds of area displayed on screen.

CENTER - sets center point of display window.

MAGNFY - sets display magnification factor.

MOVEA3 - move to specified coordinates (3-D).

VECTOR - draw line to specified coordinates and
draw an arrowhead (3-D).

CHARTK - write text on plot.

*********************

IMPLICIT REAL (A-H,O-Z)

IMPLICIT INTEGER (I-N)

COMMON /A/GAMMA,COMES,PHI,MUP,ADHES,VELTL,DEPTHF
COMMON /B/VX,VY,VZ,AX,AY,AZ,MAG,CX,CY

COMMON /LUN/LUNTI,LUNTO

COMMON /LUND/LUNDA,LUNDB

DIMENSION PMS(ll,3,ll),PMT(ll,3,ll),PST(ll,3,ll),PSS(11,3,ll),

+XM(4.“).YM(4.“).ZM(4.“).XS(4.“).YS(“.4).ZS(“.“)

DIMENSION DA(3,SO)

CMARACTER*1 Q,GRAPH

CHARACTER SOILTY*“O

CHARACTER*80 FOFILE,DATA

REAL VX,VY,VZ,AX,AY,AZ,MAG,CX,CY

REAL GAMMA,COHES,PHI,MUP,VELTL,DEPTHF
INTEGER OLM,OLS,BAUD,Q1

DATA LUNTI,LUNTO/l,l/
DATA LUNDA,LUNDB,LUNDC/12,11,13/

WRITE (LUNTI,*) 'Is a grpahics terminal being used?(y/n).'
READ (LUNTO,'(A3)') GRAPH

IF (GRAPH.EQ.'Y'.OR. GRAPH.EQ.'y') CALL NEWPAG

WRITE (LUNTI,*) 'Enter baud rate.‘

READ (LUNTO,*) BAUD

EAUD-EAUD/lo

Open file for writing tool info.

WRITE (LUNTI,*) ' Enter file name for writing tool force',

+ ' information.‘

READ (LUNTO,'(A)') FOFILE
OPEN (LUNDA,FILE-FOFILE,STATUSI'NEW',ERR'“)

COO

WOOD

000

@000

+

B“

Get tool geometry data, soil and tool operating conditions
and set plotting conditions to default.

CALL GETDAT (SIZEX,SIZEY,SIZEZ,XS,YS,ZS,XM,YM,ZM)

CALL OPERAT (SOILTY)

CALL DEFPLT(OLM,JJM,OLS,JJS)

Calculate points for share.

CALL
CALL
CALL
CALL
CALL
CALL
CALL
CALL
CALL
CALL
CALL
CALL
CALL

EDGE
EDGE
EDGE
EDGE
EDGE
EDGE
EDGE
EDGE
EDGE
EDGE
EDGE
EDGE

(xs,1,1)
(xs,1,4)
(XS,2,1)
(XS,2,“)
(YS,1,1)
(YS,1,“)
(YS,2,1)
(YS,2,“)
(23,1,1)
(23,1,4)
(28,2,1)
(23.2.4)

SURFAC (OLS,JJS,PSS,PST,XS,YS,ZS)

Calculate points for moldboard.

CALL EDGE (XM,1,1)
CALL EDGE (YM,1,1)
CALL EDGE (ZM,1,1)
CALL SURFAC (OLM,JJM,PMS,PMT,XM,YM,ZM)

Write menu to

WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE

(LUNTI,*)
(LUNTI,*)
(LUNTI,*)
(LUNTI,*)
(LUNTI,*)
(LUNTI,*)
(LUNTI,*)
(LUNTI.*)
(LUNTI,*)
(LUNTI,*)
(LUNTI,*)
(LUNTI,*)
(LUNTI,*)
(LUNTI,*)
(LUNTI.*)

lunti.

1)
2)
3)
4)
5)
6)
7)
8)
9)
10)
11)
12)
13)
14)
15)

Plot'

Change point of view.‘

Change look at point.’

Change magnification factor.‘

Change center position.‘

Set plotting conditions to default.’
Label control points on plot.'

Plot measured points.’

Change number of surface e1ements.’
Change operating conditions.’

Change share, moldboard, or scaling.‘
Get new tool geometry data.’
Get new operating data.‘
Compute forces.‘

Store plot in file and end.‘

IF (GRAPH.EQ.'Y'.0R. GRAPH.EQ.'y') CALL ANMODE

WRITE (LUNTI,*) 'Enter number of desired procedure',

'(99 for menu).'

READ (LUNTO,*,ERR-10) Q1

35
IF (Q1.EQ.99) COTO 99

GOTO (6,15,16,25,30,35,40,7S,“5,50,60,55,80,70,100),Ql
GOTO 10

Set up and plot tool.

@000

IF (GRAPH.EQ.'N'.OR.GRAPH.EQ.'n') GOTO 10
CALL INITT3 (BAUD)

7 CALL SETPLT (SIZEX,SIZEY,SIZEZ)

CALL PLOT (PSS,PST,OLS,JJS)

CALL PLOT (PMS,PMT,OLM,JJM)

CALL ANMODE

CALL HOME

GOTO 10

Change view point.

-c>c:c>

5 WRITE (LUNTI,'(A2“,F4.l,2(Al,F4.l),A32)')
+ ' Current view point is (',VX,',',VY,',',VZ,
+ '). Enter new X,Y,Z coordinates.‘
READ (LUNTO,*) vx,vy,vz
GOTO 10

Change lookat point.

0-‘000

6 WRITE (LUNTI,'(A26,F“.l,2(Al,F4.l),A32)')
+ ' Current lookat point is (',AX,’,',AY,',',AZ,
+ '). Enter new X,Y,Z coordinates.‘
READ (LUNTO,*) Ax,AY,Az
GOTO 10

Change magnification factor.

N000

5 WRITE (LUNTI,'(A33,F3.1,A20)')
+ ' Current magnification factor is ',MAG,
+ '. Enter new factor.‘
READ (LUNTO,*) MAG
GOTO 10

Change screen center point.

90000

0 WRITE (LUNTI,'(A25,F“.1,A1,F4.1,Al9)')
+ 'Current center point is (',CX,',',CY,
+ '). Enter new X,Y coordinates.‘
READ (LUNTO,*) CX,CY
GOTO 10

Reset plotting conditions

w000

5 WRITE (LUNTI,*) 'Plotting conditions set to default.‘
CALL DEFPLT (OLM,JJM,OLS,JJS)
GOTO 5

00

Label moldboard surface control points.

#000

5

“6

(13000 U'I000 W000

@000

61

IF (GRAPH.EQ.'N'
CALL RECOVR

CALL LABPTS (XM,

CALL CARTN
CALL ANMODE
GOTO 10

Change number

WRITE (LUNTI,*)
+ ' share edge.’

BO

.OR.GRAPH.EQ.'n') GOTO 10

YM,ZM)

of surface elements.

' Enter number of elements along the',

READ (LUNTO,'(I2)',ERR-“5) OLS

IF (OLS.LT.1.0R.
OLS‘OLS+1
JJSIOLS

WRITE (LUNTI , *)

OLS.GT.10) GOTO 45

' Enter number of elements along the',

+ ' moldboard edge.‘
READ (l,'(I2)',ERR=46) OLM

IF (OLM.LT.1.0R.
OLM-OLM+1
JJM-OLM

IF (GRAPH.EQ.'Y'
GOTO 10

OLM.GT.10) GOTO 46

.OR.CRAPH.EQ.'y') GOTO 5

Change operating conditions.

CALL PARAM (SOILTY)

IF (GRAPH.EQ.'Y'
GOTO 10

.OR.GRAPH.EQ.'y‘) CALL NEWPAG

Get new tool geometry data.

CALL GETDAT (SIZEX,SIZEY,SIZEZ,XS,YS,ZS,XM,YM,ZM)

IF (GRAPH.EQ.'Y'
GOTO 10

.OR.GRAPH.EQ.'y') GOTO 5

Get new operating data.

CALL OPERAT (SOILTY)

WRITE (LUNTI , *)
GOTO 10

Change share,

WRITE (LUNTI,*)
WRITE (LUNTI,*)
WRITE (LUNTI,*)
WRITE (LUNTI,*)
WRITE (LUNTI,*)

' New data read.’

moldboard, or scaling.

' Do you want to change:'
' 1) share'

' 2) moldboard'

' 3) scaling'

'4) stop'

READ (LUNTO,*) Ql
GOTO (61,62,63,64).Q1

GOTO 60

CALL NEWPS (XS,YS,ZS)

N000

N000

0U)

F000

B7

GOTO 60

CALL NEWPM (XM,YM,ZM)

GOTO 60

CALL NEWSCA (XS,YS,ZS,XM,YM,ZM,SIZEX,SIZEY,SIZEZ)
GOTO 60

WRITE (LUNTI,*) ' Do you want to store modified data',
+ ' in file?(y/n)'

READ(LUNTO, ' (A3) ') Q

IF (Q.EQ.'Y'.OR.Q.EQ.'y') THEN

CALL STODAT (SIZEX,SIZEY,SIZEZ,XS,YS,ZS,XM,YM,ZM)

ENDIF

IF (GRAPH.EQ.'Y'.OR.GRAPH.EQ.'y') GOTO 5

GOTO 10

Calculate forces.

IF (GRAPH.EQ.'Y'.OR.GRAPH.EQ.'y') CALL NEWPAG
CALL COMFOR (PMS,PSS,OLS,JJS,OLM,JJM,SOILTY)
GOTO 10

Plot measured points.

IF (GRAPH.EQ.'N'.OR.GRAPH.EQ.'n') GOTO 10
WRITE (LUNTI,*) ' Enter data file'name.’
READ (LUNTO,'(A10)',ERR-7S) DATA
OPEN (LUNDC,FILE=DATA,STATUS='OLD',ERR=75)
NUMP-l
READ (LUNDC,*,END-2) J,(DA(I,NUMP),I=1,3)
NUMP‘NUMP+1 '
GOTO l
NUMP-NUMP-l
CLOSE (LUNDC)
CALL RECOVR
CALL NEWPAG
CALL SETPLT (SIZEX,SIZEY,SIZEZ)
Plot points.
DO 3 I=l,NUMP
CALL MOVEA3 (DA(1,I),DA(2,I),DA(3,I))
CALL MOVER3 (-.OS*SIZEX,0.,0.)
CALL DRAWR3 (.l*SIZEX,0.,0.)
CALL MOVER3 (-.05*SIZEX,-.05*SIZEY,O.)
CALL DRAWR3 (0.,.l*SIZEY,O.)
CALL MOVER3 (0.,-.OS*SIZEY,-.OS*SIZEZ)
CALL DRAWR3 (O.,O.,.l*SIZEZ)
CONTINUE

GOTO 7

End program.
CALL CLOSTK (IERR)
CALL HOME

CALL ANMODE
WRITE (LUNTI,*) ' End PLOW program.’

88

CLOSE (LUNDA)
CLOSE (LUNTI)
CLOSE (LUNTO)

C)

END

INIT

********************

C)C)C§C3(5

SUBROUTINE DEFPLT (OLM,JJM,OLS,JJS)
This subroutine sets plotting parameters to default
Variable definitions:

VX,VY,VZ - view point coordinates.

AX,AY,AZ - look at point coordinates.

MAG - magnification factor.

CX,CY - center position coordinates.

OLM,OLS - number of horizontal lines drawn on moldboard
and share respectively.

JJM,JJS - number of vertical lines drawn on moldboard
and share respectively.

(if)C)C)C)C)C)C)C)C)C)CDC)C)

COMMON /B/VX,VY,VZ,AX,AY,AZ,MAG,CX,CY
REAL VX,VY,VZ,AX,AY,AZ,MAG,CX,CY
INTEGER OLM,OLS

VX'-l.

W. l e

VZ'l.

AX'O.

AY'O e

AZ‘O.

HAG’1.5

CX'-l.

CY'.5

0LM'4

JJM-4

0LS=2

JJS=2

RETURN

END
C
C********************
C

SUBROUTINE OPERAT (SOILTY)

This subroutine sets soil and operating conditions to
values read from user specified file.

C)C)C)C)C)

Variable definitions:

89

C
C DFILE - name of Operating conditions file
C (character string).
C
C SOILTY - name of soil type (character string).
C GAMMA - soil bulk density (kg/cu m).
C COHES - soil cohesion (kpa).
C PHI - soil internal friction (deg).
C MUP - soil-metal friction.
C ADHES - soil-metal adhesion (kpa).
C VELTL - tool velocity (m/s).
C DEPTHF - plowing depth (m).
C
C LUNTO - logical unit number for terminal output.
C LUNTI - logical unit number for terminal input.
C LUNDI - logical unit number for data input.
C
COMMON IA/GAMMA,COHES,PHI,MUP,ADHES,VELTL,DEPTHF
COMMON /LUN/LUNTI,LUNTO
CHARACTER SOILTY*“O
CHARACTER DFILE*80
REAL GAMMA,COHES,PHI,MUP,ADHES,VELTL,DEPTHF
DATA LUNDI/lO/
C

10 WRITE (LUNTI,*) 'ENTER NAME OF OPERATING CONDITIONS FILE.’
READ (LUNTO,'(A10)') DFILE

OPEN (UNIT-LUNDI,FILE=DFILE,STATUS='OLD',ERR=lO)
READ (LUNDI,'(A20)') SOILTY

READ (LUNDI,*) GAMMA

READ (LUNDI,*) COMES

READ (LUNDI,*) PHI

READ (LUNDI,*) MUP

READ (LUNDI,*) ADHES

READ (LUNDI ,*) VELTL

READ (LUNDI,*) DEPTHF

CLOSE(LUNDI)
RETURN
END
c
C********************
C
SUBROUTINE GETDAT (SIZEX,SIZEY,SIZEZ,XS,YS,ZS,XM,YM,ZM)

This subroutine reads tool geometry data from user
specified file.

Array definitions:

XS,YX,ZS - coordinates of share descritive points.
XM,YM,ZM - coordinates of moldboard descriptive points.

0000000000

Variable definitions:

C5C5C3C3CDC3C3C)

C

B10

DFILE - name of file to be read (character string).
SIZEX,SIZEY,SIZEZ - scaling factors for x,y,z data.

LUNTO - logical unit number for terminal output.
LUNTI - logical unit number for terminal input.
LUNDI - logical unit number for data input.

COMMON /LUN/LUNTI,LUNTO

DIMENSION XS(“,“),YS(4,4),ZS(4,“),XM(“,4),YM(4,4),ZM(“,“)
CHARACTER DFILE*80

DATA LUNDA/IO/

WRITE (LUNTI,*) 'ENTER NAME OF GEOMETRY DATA FILE TO BE READ.’
READ (LUNTO,'(A20)') DFILE

OPEN (LUNDI,FILE=DFILE,STATUSt'OLD',ERRBIO)
READ (LUNDI,'(3(F9.3))') SIZEX,SIZEY,SIZEZ

DO 20 I=l,“
Do 20 J=1,4
READ (LUNDI,'(3(F9.3))') XS(J,I),YS(J,I),ZS(J,I)
CONTINUE

DO 30 1'1,“
DO 30 J-l,“
READ (LUNDI,'(3(F9.3))') XM(J,I),YM(J,I),ZM(J,I)
CONTINUE

CLOSE (LUNDI)
RETURN
END

C********************

C

C3C3C5C3CIC3C3C363CDCDC)CDCDC)

SUBROUTINE STODAT (SIZEX,SIZEY,SIZEZ,XS,YS,ZS,XM,YM,ZM)

This subroutine writes current tool geometry data
to user specified file.

Array definitions:

XS,YX,ZS - coordinates of share descritive points.
XM,YM,ZM - coordinates of moldboard descriptive points.

Variable definitions:

DATA - name of file to be read (character string).
SIZEX,SIZEY,SIZEZ - scaling factors for x,y,z data.

LUNTO - logical unit number for terminal output.

Bll

LUNTI - logical unit number for terminal input.
LUNDI - logical unit number for data input.

000

COMMON /LUN/LUNTI,LUNTO

DIMENSION xs(4,4),Ys(4,4),zs(4,4),xn(4,4),YM(4,4),ZM(4,4)
CHARACTER DFILE*80

DATA LUNDI/lO/

lO WRITE (LUNTI,*) ' ENTER NAME OF TOOL GEOMETRY DATA FILE TO ',
+ 'BE WRITTEN TO.’
READ (LUNTO,'(AIO)') DFILE

C

OPEN (LUNDI,PILE-DPILE,STATUS='NEw',ERR-lo)
C

WRITE (LUNDI,'(3(F9.3))') SIZEX,SIZEY,SIZEZ
C
C

no 20 1:1,4

Do 20 J-l,“

WRITE (LUNDI,'(3(F9.3))') XS(J,I),YS(J,I),ZS(J,I)
20 CONTINUE

DO 30 I=l,“
DO 30 J=l,“
WRITE (LUNDI,'(3(F9.3))') XM(J,I),YM(J,I),ZM(J,I)
30 CONTINUE

C
C
CLOSE (LUNDI)
RETURN
END
C
C
C SURFACE
C
C********************
C

SUBROUTINE SURFAC (OL,JJ,PTS,PTT,PX,PY,PZ)

This subroutine calculates point coordinates to be plotted.

Array definitions:

PTS,PTT - coordinates of points for describing surface.
PX,PY,PZ - coordinates of control points.

Variable definitions:
OL,JJ - number of splines used to draw surface.

Subroutines called:

000000000000000

BlZ

INTNOD - calculates internal control points.
POINT - calculates point for describing surface.

CDCDC)

DIMENSION PTS(ll,3,ll),PTT(ll,3,ll),PX(“,4),PY(“,“),PZ(“,4)
REAL U,V
INTEGER OL,JJ

CALL INTNOD (PX)
CALL INTNOD (PY)
CALL INTNOD (PZ)

C BLEND IN T-DIRECTION

DO 10 I-1,OL

u-(I-1)/REAL(0L-1)

DO 10 J-1,JJ
v-(J-1)/REAL(JJ-1)
CALL POINT (PX,U,V,PTT(J,1,I))
CALL POINT (PY,U,V,PTT(J,2,I))
CALL POINT (PZ,U,V,PTT(J,3,I))

o CONTINUE

BLEND IN S-DIRECTION

C)(>C)h-

DO 20 I=1,0L
v=(I-1)/REAL(0L-1)
Do 20 J-l,JJ
U=(J-1)/REAL(JJ-1)
CALL POINT (Px,U,v,PTs(J,1,I))
CALL POINT (PY,U,V,PTS(J,2,I))
CALL POINT (PZ,U,V,PTS(J,3,I))
20 CONTINUE
RETURN
END
C
C*********************
C
SUBROUTINE INTNOD (P)

This subroutine calculates the four internal control points.
Array definition:

P - coordinates of control points

(363(363C5C3C)

DIMENSION P(4,4)
P(2,2)-P(2,l)+P(l,2)-P(l,l)
P(2,3)-P(l,3)-P(l,“)+P(2,“)
P(3,2)-P(3,l)-P(“,l)+P(“,2)
P(3,3)-P(“,3)+P(3,“)-P(“,4)
RETURN

END

C
C*********************

C)

CDCDC)C3C3C5C3C3C5C3C3CDCDC)C3C3C3C3C)

C

B13

SUBROUTINE POINT (P,S,T,Q)
This subroutine calculates a point for surface description.
Array definitions:

P - control points.

SMAT - 8 matrix.

TMAT - T matrix.
M - M matrix.

SM - S*M
SMP - S*M*P
MT - M*T

Q - calculated surface point.
Subroutine called:

MMLT - matrix multiplication (Case Center library
function).

DIMENSION P(“,4),SMAT(“),TMAT(“),SM(4),SMP(4),Q(l)
REAL M(“,4),MT(“)
DATA Ml-le ,3. ,‘30 ,1. ,3. ,“60 ,3. ’00 ,‘30 ,3. ,0. ,0. ,le ,0. ,0. ,0./

SMAT(1)-s*s*s
SMAT(2)=S*S
SMAI(3)-s
SMAT(4)-1

TMAT(l)-T*T*T
TMAT(2)=T*T
TMAT(3)=T
TMAT(“)-l

CALL MMLT (SM,SMAT,M,1,“,“)
CALL MMLT (MT,M,TMAT,4,“,1)
CALL MMLT (SMP,SM,P,1,“,“)
CALL MMLT (Q,SMP,MT,1,“,1)

RETURN
END

c********************

C

C)C)C)C)CIC)C)C)

SUBROUTINE EDGE (P,I,J)

This subroutine calculates control points along a
straight edge.

Array definitions:

P - control point coordinates.

Bl“

DIMENSION P(4,4)

IF (I.EQ.1) THEN
P(J.2)=P(J.l)+((P(J.4)-P(J.l))/3.)
P(J,3)-P(J,l)+((P(J,“)-P(J,1))*2./3.)

ELSE
P(2.J)-P(l.J)+((P(4,J)-P(1,J))/3.)
P(3.J)=P(1.J)+((P(4.J)-P(1.J))*2./3.)

ENDIF

RETURN

END

PLOT

*******************

09000

SUBROUTINE SETPLT (SIZEX,SIZEY,SIZEZ)

This subroutine sets up the plot and draws the axes for
subroutine PLOT.

Variable definitions:

SIZEX,SIZEY,SIZEZ - scaling factors for controls points.
VX,VY,VZ - view point coordinates.

AX,AY,YZ - look at point coordinates.

MAG - magnification factor.

CX,CY - center position coordinates.

SIZE - largest of SIZEX,SIZEY,SIZEZ.
SIZEA - average of SIZEX,SIZEY,SIZEZ.

The following subroutines are on the Prime system at The
Case Center at Michigan State University.

CARTVP - sets cartisan view point.

LOOKAT - set lookat point.

PERSPT - turns on or off perspective viewing.

XUP,YUP,ZUP - defines which axis is up on plot.

PPCENT - sets distance to 2-D viewing plane from
look at point to view point as a
percentage.

DWINDO - sets bounds of area displayed on screen.

CENTER - sets center point of display window.

MACNFY - sets display magnification factor.

MOVEA3 - move to specified coordinates (3-D).

VECTOR - draw line to specified coordinates and
draw an arrowhead (3-D).

CHARTK - write text on plot.

000000000000000000000000000000000

COMMON /B/VX,VY,VZ,AX,AY,AZ,MAG,CX,CY
REAL MAG

SIZE-~10000.

IF (SIZEX.CE.SIZE) SIZE-SIZEX

IF (SIZEY.GE.SIZE) SIZE-SIZEY

C5CDC)

C

B15

IF (SIZEZ.GE.SIZE) SIZE=SIZEZ
SIZEA=(SIZEX+SIZEY+SIZEZ)/3

SET UP PLOT

CALL CARTVP (VX,VY,VZ)
CALL LOOKAT (AX,AY,AZ)
CALL PERSPT (.FALSE.)

IF (VX.EQ.O..AND.VY.EQ.O.) THEN

CALL XUP

ELSE

CALL ZUP

ENDIF

CALL
CALL
CALL
CALL

CALL
CALL
CALL
CALL
CALL
CALL
CALL
CALL
CALL
CALL
CALL
CALL

PPCENT
DWINDO
CENTER
MAGNFY

(100.)

(-3.*SIZE,2.*SIZE,-2.*SIZE,3.*SIZE)

(CX*SIZE,CY*SIZE)
(MAG)

DRAW AXES

MOVEA3
VECTOR
MOVEA3
CHARTK
MOVEA3
VECTOR
MOVEA3
CHARTK
MOVEA3
VECTOR
MOVEA3
CHARTK

RETURN

END

(0.,0.,o.)
(2.*SIZEX,0.,O.)
(2.1*SIZEX,O.,O.)
('X',l.5)
(0.,0.,o.)
(O.,2.*SIZEY,O.)
(O.,2.l*SIZEY,O.)
('Y',l.5)
(0.,0.,o.)
(O.,0.,l.5*SIZEZ)
(O.,0.,l.6*SIZEZ)
('Z',1.5)

C*******************

C

CDC)C3C3C3C363CIC3CDC36363

SUBROUTINE PLOT (PS,PT,OL,JJ)

This subroutine plots the surfaces on the region set up

in subroutine SETPLT.

Array definitions:

PS,PT - contains point locations for 3-D plots
Variable definitions:

OL,JJ - number of splines used in PS and PT to
draw the surface.

DIMENSION PS(ll,3,ll),PT(ll,3,ll)
INTEGER 0L,JJ

CECDCEC)

non.—

20

C

816

DRAW SURFACE

[First draw splines in S direction

DO 10 I=1,0L
Do 10 J=1,JJ
IF (J.EQ.l) THEN
CALL MOVEA3 (PS(J,1,I),PS(J,2,I),PS(J,3,I))
ELSE
CALL DRAWA3 (PS(J,1,I),PS(J,2,I),PS(J,3,I))
ENDIF
CONTINUE

Then in T direction

00 20 I-1,0L
DO 20 J-1,JJ
IF (J.EQ.1) THEN
CALL MOVEA3 (PT(J,l,I),PT(J,2,I),PT(J,3,I))
ELSE
CALL DRAWA3 (PT(J,1,I),PT(J,2,I),PT(J,3,I))
ENDIF
CONTINUE
RETURN
END

C*******************

C

00000000000‘00

000

SUBROUTINE LABPTS (PX,PY,PZ)
This subroutine labels the moldboard control points.
Array definitions:

PX,PY,PZ - contains coordinates of control points.
Variable definitions:

Q - answer to question (character string).

Q1 - numbers to plotted (character string).

DIMENSION PX(“,“),PY(“,4),PZ(4,4)
CHARACTER Q*3,Ql*32
Ql-' 1 2 3 4 5 6 7 8 910111213141516'

LABEL POINTS

Do 2 1-1,4
DO 2 J-l,“
CALL MovEA3 (PX(J,I),PY(J,I),PZ(J,I))
Q-Ql(((((I-1)*4+J)*2)-l):(((I-l)*4+J)*2))
CALL CHARTK (Q,l.)
CONTINUE

Bl7

RETURN
END

TRANSFORM

*********************

00000

SUBROUTINE DISTAN (PS,I,J,X,Y,Z,ANGLE1,ANGLE2,DIST)

This subroutine calculates the distance the soil
travels over an element. These calculations
require translating and rotating elements.

Array definitions:
PS - coordinates of surface description points.
Variable definitions:

I,J - pointers indicating which element in P8 is
to be translated and rotated. Points
PS(J,*,I),PS(J+1,*,I),PS(J,*,I+1) and
PS(J+1,*,I+1) form the element.

X - x direction cosine of element orientation.

Y - y direction cosine of element orientation.

Z - z direction cosine of element orientation.

ANGLEl - the first angle of rotation.

ANGLEZ - the second angle of rotation.

DIST - the distance the soil travels over the element.

Subroutines called:

TRANSL - translates point.
ROTATE - rotates point.

000000000000000000000000000

DIMENSION PS(ll,3,ll)
REAL LEN

XS-X
Y5=Y
ZS-Z

CALL TRANSL (PS(J,1,I),PS(J,2,I),PS(J,3,I),X1,Y1,Zl,
+-PS(J,l,I),-PS(J,2,I),-PS(J,3,I))

CALL TRANSL (PS(J+l,l,I),PS(J+l,2,I),PS(J+1,3,I),X2,Y2,22,
+-PS(J,l,I),-PS(J,2,I),-PS(J,3,I))

CALL TRANSL (PS(J+l,1,I+l),PS(J+l,2,I+l),PS(J+l,3,I+l),X3,Y3,Z3,
+-PS(J,l,I),-PS(J,2,I),-PS(J,3,I))

CALL TRANSL (PS(J,1,I+l),PS(J,2,I+l),PS(J,3,I+1),X4,Y4,Z“,
+-PS(J,l,I),-PS(J,2,I),-PS(J,3,I))

CALL ROTATE (X1,Yl,Zl,Xl,Yl,Zl,'Z',ANGLEl)
CALL ROTATE (X2,Y2,ZZ,X2,Y2,22,'Z',ANGLEI)
CALL ROTATE (X3,Y3,Z3,X3,Y3,Z3,'Z',ANGLEl)

C

CALL
CALL

CALL
CALL
CALL
CALL
CALL

318

ROTATE (X4,Y4,Z“,X“,Y4,Z“,'Z',ANGLE1)
ROTATE (XS,Y5,ZS,XS,Y5,ZS,'Z',ANGLE1)

ROTATE (X1,Yl,Zl,Xl,Yl,Zl,'Y',ANGLEZ)
ROTATE (X2,Y2,22,X2,Y2,ZZ,'Y',ANGLEZ)
ROTATE (X3,Y3,Z3,X3,Y3,23,'Y',ANGLE2)
ROTATE (X4,Y4,Z“,X4,Y“,Z“,'Y',ANGLE2)
ROTATE (XS,Y5,ZS,X5,YS,ZS,'Y',ANGLEZ)

DIMY=(AES(Y1-T4)+AES(Y2-Y3))/2.
DIMZ-(ABS(Zl-ZZ)+ABS(23-Z“))/2.
LEN-SQRT(DIMY**2+DIMZ**2)
DIST-DIMz/zs

IF (DIST.GT.LEN) DIST-DIMY/YS
IF (DIST.GT.LEN) DIST-LEN

RETURN

END

C********************

C

0000000000000000000000

0

SUBROUTINE TRANEL (PS,I,J,ANGLE1,ANGLEZ,DIMY,DIMZ)

This subroutine calculates the y and z dimension of an
element.

Array definitions:

PS - surface desciption coordinates.

Variable definitions:

I,J - pointers indicating which element in PS is
to be translated and rotated.

ANGLEI - the first angle of rotation.

ANGLEZ - the second angle of rotation.

DIMY - calculated y dimension.

DIME - calculated 2 dimension.

Subroutines called:

TRANSL - translates element.
ROTATE - rotates element.

DIMENSION PS(ll,3,ll)

CALL

TRANSL (PS(J,1,I),PS(J,2,I),PS(J,3,I),X1,Yl,Zl,

+-PS(J,l,I),-PS(J,2,I),-PS(J,3,I))

CALL

TRANSL (PS(J+1,1,I),PS(J+1,2,I),PS(J+1,3,I),X2,Y2,Z2,

+-PS(J,l,I),-PS(J,2,I),-PS(J,3,I))

CALL

TRANSL (PS(J+1,1,I+l),PS(J+1,2,I+1),PS(J+1,3,I+l),X3,Y3,Z3,

+-PS(J,1,1),-PS(J,2,I),-PS(J,3,I))

CALL

TRANSL (PS(J,1,I+1),PS(J,2,I+1),PS(J,3,I+1),X4,Y“,Zé,

+-PS(J,1,1),-PS(J,2,I),-PS(J,3,I))

Bl9

(X1,Yl,Zl,Xl,Yl,Zl,'Z',ANGLEl)
(X2,Y2,Z2,X2,Y2,ZZ,'Z',ANGLEI)
(X3,Y3,23,X3,Y3,23,'Z',ANGLEl)
(X4,Y“,Z4,X4,Y“,Z“,'Z',ANGLE1)

(X1,Yl,Zl,Xl,Yl,Zl,'Y',ANGLEZ)
(X2,Y2,ZZ,X2,Y2,ZZ,'Y',ANGLEZ)
(X3,Y3,Z3,X3,Y3,Z3,'Y',ANGLE2)
(x4,Y4,z4,x4,Y4,z4,'Y',ANGLE2)

DIMY'(ABS(Yl—Y“)+ABS(Y2-Y3))/2.
DIMZ-(ABS(Zl-ZZ)+ABS(23-Z“))/2.

SUBROUTINE ROTATE (X,Y,Z,X1,Yl,Zl,AXIS,ANG)

This subroutine rotates a point.

Variable definitions:

point to be rotated.

Xl,Yl,Zl - point after rotation.

AXIS - axis point is rotated about (character

string).

ANG - angle point is to be rotated.

Subroutines called:

ROT - calculates T matrix.

MMLT - matrix multiplication (Case Center Math

library function).

C
CALL ROTATE
CALL ROTATE
CALL ROTATE
CALL ROTATE

C
CALL ROTATE
CALL ROTATE
CALL ROTATE
CALL ROTATE

C

C
RETURN
END

C

c********************

C

C

C

C

'C

C

C X,Y,Z -

C

C

C

C

C

C

C

C

C

C

C

DIMENSION T(4,4),P(4),RP(4)

CHARACTER AXIS*1

Do 10 1-1,4
DO 10 J-l,“

IF (I.EQ.J) THEN

T(I,J)-l.
ELSE
T(I.J)-0.

ENDIF
10 CONTINUE

P(1)-x
P(2)-Y
P(3)-z
P(4)-l.

IF (AXIS.EQ.'X') CALL ROT (T,2,3,ANG)

C

320

IF (AXIS.EQ.'Y') CALL ROT (T,l,3,ANG)
IF (AXIS.EQ.'Z') CALL ROT (T,l,2,ANG)

CALL MMLT(RP,P,T,1,4,“)
XI-RP(1)
Yl-RP(2)
Zl-RP(3)

RETURN
END

c*********************

C

(EC)C5C3CDCSC5C363C5C5CDCMO

C)

C

SUBROUTINE ROT (T,l,J,ANG)
This subroutine calculates T matrix.
Array definitions:

T - work array for assigning values to used in
subroutine ROTATE.

Variable definitions:

I,J - pointers in T matrix.
ANG - angle of rotation.

DIMENSION T(“,4)

T(I,I)-COS(ANG)
T(I,J)-SIN(ANC)
T(J,I)--SIN(ANG)
T(J,J)-COS(ANG)

IF (I.EQ.1.AND.J.EQ.3) THEN
T(I,J)--SIN(ANG)
T(J,I)-SIN(ANG)

ENDIF

RETURN
END

C********************

C

C)C)C)C)C)€)C)

SUBROUTINE TRANSL (X,Y,Z,Xl,Yl,Zl,x2,Y2,ZZ)
This subroutine translated a point.
Variable definitions:

X,Y,Z - point to be translated.
X1,Yl,Zl - translated point.

821
X2,Y2,22 - distance to translated.
Subroutine called:

MMLT - matrix multiplication (Case Center Math
library function).

(3C5C3C3C3C10

DIMENSION T(“,“),P(4),TP(4)

C5

D0 10 III,“
DO 10 J-l,“

IF (I.EQ.J) THEN
T(I,J)-1.

ELSE
T(I,J)=0.

ENDIF

10 CONTINUE

P(1)-x
P(2)=Y
P(3)-z
P(4)-l.

T(4,1)-x2

T(4,2)-Y2

T(4,3)-zz

CALL MMLT (TP,P,T,1,4,4)
x1-TP(1)

Yl-TP(2)

Zl-TP(3)

RETURN
END

NEWPTS

*************************

C1C3CDC3C)

SUBROUTINE NEWSCA (XS,YS,ZS,XM,YM,ZM,SIZEX,SIZEY,SIZEZ)

This subroutine changes the scaling factors of the
control points.

Array definitions:

XS,YX,ZS - coordinates of share control points.
XM,YM,ZM - coordinates of moldboard control points.

Variable definitions:

SIZEX,SIZEY,SIXEZ - scaling factors for control
points.

(3(5C3C5C3C3C3C3C3C3C)C3C3

322

SIZEXI,SIZEY1,SIZEZI - intermediate values for
SIZEX,SIZEY,SIZEZ.

LUNTO - logical unit number for terminal output.
LUNTI - logical unit number for terminal input.

000000

COMMON /LUN/LUNTI,LUNTO
DIMENSION xs(4,4),Ys(4,4),zs(4,4),XM(4,4),YM(4,4),ZM(4,4)
SIZEX1=SIZEX
SIZEYl-SIZEY
SIZEZl-SIZEZ
WRITE (LUNTI,'(A,F9.3)') 'x-SCALE- ',SIZEXI
WRITE (LUNTI,'(A,F9.3)') 'Y-SCALE- ',SIZEYl
WRITE (LUNTI,'(A,F9.3)') 'Z-SCLAE= ',SIZEZI
5 WRITE (LUNTI,*) 'ENTER NEW x-SCALE, Y-SCALE, z-SCALE.'
READ (LUNTO,*,ERR-5) SIZEX,SIZEY,SIZEZ

Do 10 1-1,4
Do 10 J-l,“
XS(I,J)-XS(I,J)*SIZEX/SIZEX1
YS(I,J)-YS(I,J)*SIZEY/SIZEY1
ZS(I,J)-ZS(I,J)*SIZEZ/SIZEZl
XM(I,J)-XM(I,J)*SIZEX/SIZEX1
YM(I,J)-YM(I,J)*SIZEY/SIZEY1
ZM(I,J)-2M(I,J)*SIZEZ/SIZE21
10 CONTINUE '

C
RETURN
END
C
C**********************
C
SUBROUTINE NEWPM (PX,PY,PZ)
C
C This subroutine changes the moldboard control points.
C
C Array definitions:
C
C PX,PY,PZ - coordinates of moldboard control
C points.
C
C LUNTO - logical unit number for terminal output.
C LUNTI - logical unit number for terminal input.
C
COMMON /LUN/LUNTI,LUNTO
DIMENSION PX(“,“),PY(“,“),PZ(4,4)
WRITE (LUNTI,*) ' MOLDBOARD DESCRIPTIVE POINTS'
WRITE (LUNTI,'(A2“)') 'PTI X Y 2'
C
DO 10 I-l,“
DO 10 J-l,“

IF (((I-l)*“+J).EQ.S) GOTO 10
IF (((I-l)*4+J).EQ.6) GOTO 10
IF (((I-1)*4+J).EQ.7) GOTO 10

10

15

20

C

823

IF (((I-l)*“+J).EQ.9) GOTO 10

IF (((I-1)*4+J).EQ.10) GOTO 10

IF (((I-l)*“+J).EQ.ll) GOTO 10

WRITE (LUNTI,'(IZ,3(3X,F6.3))') ((I-l)*“+J),PX(J,I),

+ PY(J,I),PZ(J,I)

CONTINUE

WRITE (LUNTI,'(/,A34)') 'ENTER PT# To BE EDITED.(0 To STOP)‘
READ (LUNTO,*,ERR=15) J

IF (J.GT.l6) GOTO 15

IF (J.EQ.0) GOTO 20

I-((J-l)/4)+l

J=MOD(J-l,“)+l

WRITE (LUNTI,'(A11)') 'ENTER X,Y,Z'

READ (LUNTO,*,ERR-15) PX(J,I),PY(J,I),PZ(J,I)
GOTO 15

CONTINUE

RETURN

END

C********************

C

00000000000

10

15

SUBROUTINE NEWPS (PX,PY,PZ)
This subroutine changes the share control points.
Array definitions:

PX,PY,PZ - coordinates of moldboard control
points.

LUNTO - logical unit number for terminal output.
LUNTI - logical unit number for terminal input.

COMMON /LUN/LUNTI,LUNTO

DIMENSION PX(“,4),PY(4,4),PZ(“,“)

WRITE (LUNTI,*) ' SHARE DESCRIPTIVE POINTS.‘
WRITE (LUNTI,'(A2“)') 'PT# x Y 2'

Do 10 I-l,“
Do 10 J-l,“
IF (((I-I)*4+J).EQ.1) GOTO 5
IF (((1—1)*4+J).EQ.4) GOTO 5
IF (((I-I)*4+J).EQ.13) GOTO 5
IF (((I-l)*“+J).EQ.l6) GOTO 5
GOTO 10
WRITE (LUNTI,'(12,3(3X,F6.3))') ((I—1)*4+J),Px(J,I),

+ PY(J,I),PZ(J,I)

CONTINUE

WRITE (LUNTI,'(/,A34)') 'ENTER PT# TO BE EDITED.(O TO STOP)‘
READ (LUNTO,*,ERR-15) J

IF (J.GT.16) GOTO 15

IR (J.EQ.0) GOTO 20

824

I=((J-l)/“)+l

J'MOD(J-l,4)+l

WRITE (LUNTI,'(A11)') 'ENTER X,Y,Z'

READ (LUNTO,*,ERR=15) PX(J,I),PY(J,I),PZ(J,I)

GOTO 15
20 CONTINUE
RETURN
END
C
c PARAM
C
C********************
C

SUBROUTINE PARAM (SOILTY)

This subroutine allow editing of the soil and tool
parameters.

Variable definitions:

DFILE - name of operating conditions file
(character string).

SOILTY - name of soil type (character string).
GAMMA - soil bulk density (kg/cu m).

COHES - soil cohesion (kpa).

PHI - soil internal friction (deg).

MUP - soil-metal friction.

ADHES - soil-metal adhesion (kpa).

VELTL - tool velocity (m/s).

DEPTHF - plowing depth (m).

LUNTI - logical unit number for terminal input.
LUNTO - logical unit number for terminal output.
LUNDI - logical unit number for data input.

Q - answer to question (character string).

00000000000000000000000

COMMON /A/GAMMA,COHES,PHI,MUP,ADHES,VELTL,DEPTHP
COMMON /LUN/LUNTI,LUNTO

CHARACTER SOILTY*“0

CHARACTER Q*3

CHARACTER DFILE*80

REAL GAMMA,COHES,PHI,MUP,ADHES,VELTL,DEPTHF

DATA LUNDI/lO/

CALL NEWPAG

WRITE (LUNTI,'(10X,A)') 'SOIL PARAMETERS'

WRITE (LUNTI,'(A,T60,A)')' 1) SOIL TYPE: ',SOILTY

WRITE (LUNTI,100)' 2) BULK DENSITY (RG/CU—M) ',GAMMA

WRITE (LUNTI,100)' 3) COHESION (kPA) ',COHES

WRITE (LUNTI,lOO)' 4) INTERNAL SOIL FRICTION ANGLE (DEG) ',PHI
WRITE (LUNTI,100)' 5) SOIL-METAL FRICTION COEFFICIENT ',MUP
WRITE (LUNTI,lOO)' 6) SOIL ADHESION (KPA)',ADHES

WRITE (LUNTI,'(10X,A)')'TOOL PARAMETERS'

10

100

C3C3C5C3C)

C5C3C3€3C3C36563C3

+

B25

WRITE (LUNTI,100)' 7) FORWARD VELOCITY (M/S) ',VELTL

WRITE (LUNTI,lOO)' 8) TOOL DEPTH (M) ',DEPTHF

WRITE (LUNTI,'(A)') 'ENTER PARAMETER NUMBER To BE CHANGED. ',
'(0 TO STOP)'

READ (LUNTO,*) I
IF (I.EQ.O) GOTO 20
WRITE (LUNTI,*) 'ENTER NEW VALUE FOR PARAMETER NUMBER ',I

IF (I.EQ.1) READ (LUNTO,'(A20)’) SOILTY
IF (I.EQ.2) READ (LUNTO,*) GAMMA

IF (I.EQ.3) READ (LUNTO,*) COHES

IF (I.EQ.4) READ (LUNTO,*) PHI

IF (I.EQ.5) READ (LUNTO,*) MUP

IF (I.EQ.6) READ (LUNTO,*) ADHES

IF (I.EQ.7) READ (LUNTO,*) VELTL

IF (I.EQ.8) READ (LUNTO,*) DEPTHF

GOTO 10

WRITE (LUNI,*) ' STORE MODIFIED DATA?(Y OR N)’
READ (LUNTO,'(A3)') Q

IF (Q.EQ.'N') GOTO 30

WRITE (LUNTI,*) 'ENTER NEW FILE NAME.’

READ (LUNTO,'(A20)') DPILE

OPEN (UNIT=LUNDI,FILE=DFILE,STATUS='NEW',ERR820)
WRITE (LUNDI,'(A20)') SOILTY

WRITE (LUNDI,*) GAMMA

WRITE (LUNDI,*) COMES

WRITE (LUNDI,*) PHI

WRITE (LUNDI ,*) MUP

WRITE (LUNDI,*) ADHES

WRITE (LUNDI,*) VELTL

WRITE (LUNDI,*) DEPTHF

CLOSE (LUNDI)

RETURN
FORMAT (A,T60,F9.3)
END

ANGLE

********************

SUBROUTINE ANGLE (PS,I,J,X,Y,ALPHA)

This subroutine calculates the angle between one of
axes and the element.

Array definitions:
PS - coordinates of surface description points.

Variable definitions:

00000000000

000

. 000

826

I,J - pointers indicating which element in P8 is
to have the orientation computed for.

X,Y - numbers indicating the plane from which the
angle is calculated.

ALPHA - calculated angle.

Subroutines called:

INTERP - interpolates between to points.

DIMENSION PS(11,3,11)
INTEGER X,Y,Z,O,P
DATA PI/3.14159265359/

Check which plane is to be considered.

IF (XOEQOIOANDOYOEQOZ) 233
IF (X.EQ.2.AND.Y.EQ.3) 2'1
IF (XOEQO3OANDOYOEQOI) 2'2

Find point closest to 0,0,0

IF (PS(J,Z,I).GE.PS(J,Z,I+l).AND.PS(J,Z,I).GE.PS(J+1,Z,I)) THEN
IF (PS(J,z,I+1).GE.PS(J+1,z,I)) THEN
K-J+l
L=J
M=J
N=I
o-I
P-I+1
ELSE
K-J
L-J
M-J+1
N-I+1
o-I
P-I
ENDIF
ELSE IF (PS(J,Z,I+1).GE.PS(J,Z,I).AND.PS(J,Z,I+1).GE.PS(J+1,Z,I))
+ THEN .
IF (PS(J,Z,I).GE.PS(J+1,Z,I)) THEN
K=J+1
L-J
M-J
NaI
o-I+1
P-I
ELSE
x-J
L-J
M-J+1
N-I
o=I+1

827

P81
ENDIF
ELSE IF (PS(J+l,Z,I).GE.PS(J,Z,I).AND.PS(J+1,Z,I).GE.PS(J,Z,I+1))
+ THEN
IF (PS(J,Z,I).GE.PS(J,Z,I+1)) THEN
K-J
L-J+l
M-J
N=I+l
0-I
P-I
ELSE
K-J
L-J+l
M-J
N-I
O-I
P-I+l
ENDIF
ENDIF
C
C Interpolate to find the point used to calculate the
C angle.
C
CALL INTERP (PS(K,Z,N),PS(L,Z,O),PS(M,Z,P),
+ PS(K,X,N),PS(L,X,O),XX)
CALL INTERP (PS(K,Z,N),PS(L,Z,O),PS(M,Z,P),
+ PS(K,Y,N),PS(L,Y,O),YY)
C
IF ((XX-PS(M,X,P)).LT..001.AND.(XX-PS(M,X,P)).GT.-.001) THEN
ALPHA-PI/Z.
ELSE
RATYx-((YY-PS(M,Y,P))/(xx-PS(M,x,P)))
IF ((YY-PS(M,Y,P))/(XX-PS(M,X,P)).LT.0.) THEN
ALPHAa-ATAN(-RATYX)
ELSE
ALPHA-ATAN (RATYX)
ENDIF
ENDIF
C
RETURN
END
C
C*******************
C
SUBROUTINE INTERP (MINI,MAXl,MIDI,MIN2,MAX2,MID2)
C
C This subroutine interpolates between to points.
C
C Variable definitions:
C
C MINI - minimum independent value.
C MAXI - maximum independent value.
C MIDI - middle independent value.

B28

C MINZ - minimum dependent value.
C MAXZ - maximum dependent value.
C MIDZ - middle dependent value (calculated).
C
REAL MINI,MAXl,MID1,MIN2,MAX2,MIDZ
C
IF ((MAXl-MINI).EQ.O.) THEN
MIDZ-MAXI
ELSE
MID2-MIN2+((MIDl-MIN1)/(MAXl-MIN1)*(MAX2-MIN2))
ENDIF
C
RETURN
END
C
C FORCES
C
c********************
C

SUBROUTINE COMFOR (PMS,PSS,OLS,JJS,OLM,JJM,SOILTY)

This subroutine calculates the forces the soil applies
to the tool.

Array definitions:

PMS - coordinates of moldboard surface description
points.

PSS - coordinates of share surface description
points.

SFORCE - Share forces.

MFORCE - Moldboard forces.

TOTFOR - Total forces.

TOTMOM - Total moments.

Variable definitions:

OLM,OLS - number of horizontal lines drawn on
moldboard and share respectively.
JJM,JJS - number of vertical lines drawn on
moldboard and share respectively.
SOILTY - soil type (character string).

LUNTI - logical unit number for terminal input.
LUNTO - logiéal unit number for terminal output.
LUNDA - logical unit number for file data storage.
LUNDB - logical unit number for file data storage.

Subroutines called:

ELEFOR - calculates element forces.
SUMFOR - sums the forces over all elements.

000000000000000000000000000000000

COMMON /A/GAMMA,COHES,PHI,MUP,ADHES,VELTL,DEPTHF

829

COMMON /LUN/LUNTI,LUNTO
COMMON /LUND/LUNDA,LUNDB

DIMENSION SFORCE(11,3,ll,lO),MFORCE(ll,3,ll,10)
+,PMS(11,3,11),PSS(11,3,ll),TOTFOR(3,2),TOTMOM(3,2)
+,XNS(11),YNS(11),ZNS(11),XNM(ll),YNM(ll),ZNM(11)

REAL MUP,MFORCE

000 0O‘Ui0

000

INTEGER OLS,JJS,OLM,JJM
CHARACTER SOILTY*40

IC=O

OPEN (LUNDB,FILE-'D_FORCES',STATUSs'UNKNOWN')
(LUNDB,'(//A)') 'SOIL AND OPERATING CONDITIONS'
(LUNDB,'(A,T20,A)') 'SOIL TYPE',SOILTY

WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE

WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE

(LUNDB,5)
(LUNDB,5)
(LUNDB,5)
(LUNDB,5)
(LUNDB,5)
(LUNDB,5)
(LUNDB,5)

(LUNDA.'(////A)') '

'SOIL BULK DENSITY (KG/CU M)',GAMMA
'SOIL COHESION (KPA)',COHES

'SOIL INTERNAL FRICTION ANGLE (DEG)',PHI
'SOIL-METAL FRICTION COEFFICIENT',MUP
'SOIL-METAL ADHESION (KPA)',ADHES
'FORWARD VELOCITY (M/S)',VELTL

'TOOL DEPTH (M)',DEPTHF

SOIL AND OPERATING CONDITIONS'

(LUNDA,'(A,T30,A)') 'SOIL TYPE',SOILTY

(LUNDA,S)
(LUNDA,5)
(LUNDA,5)
(LUNDA,5)
(LUNDA,S)
(LUNDA,5)
(LUNDA,5)
(LUNDA,6)
(LUNDA,6)

'SOIL BULK DENSITY (KG/CU M)',GAMMA

'SOIL COHESION (RPA)',COHES

'SOIL INTERNAL FRICTION ANGLE (DEG)',FHI
'SOIL-METAL FRICTION COEFFICIENT',MUP
'SOIL-METAL ADHESION (KPA)',ADHES

'FORWARD VELOCITY (M/S)',VELTL

'TOOL DEPTH (M)',DEPTHF

'NUMBER OF SHARE ELEMENTS',(JJS-l)*(OLS-l)
'NUMBER OF MOLDBOARD ELEMENTS',(JJM-1)*(0LM-1)

FORMAT (A,T“O,F10.“)
FORMAT (A,T“O,I3)

WRITE (LUNTI,*) 'ENTER LOCATION OF ZERO MOMENT.(X,Y,Z)'
READ (LUNTO,*) X,Y,Z

WRITE (LUNDB.'(////A)') '

COMPUTE FORCES ON SHARE

DO 10 J'1,JJS-l
DO 10 I-l,OLS-l

SHARE ELEMENT FORCES'

Assign force cosines.

IF (J.EQ.1) THEN

“-1 e

YN'O.

ZN'O e
ELSE

xN-XNs(I)

000

000

000

0001-0

000

15

B30

YNsYNs(I)
ZN-zNS(I)
ENDIF

IF (XNS(I).LE.0..AND.YNS(I).LE.O..AND.ZNS(I).LE.0.) THEN
II-I+1
JJ-J
GOTO 100

ENDIF

Calculate element forces.
CALL ELEFOR (PSS,I,J,SFORCE,'S',XN,YN,ZN,IC)

AVEX'(PSS(J,1,I)+PSS(J+1,1,I)+PSS(J+1,1,I+l)+PSS(J,l,I+l))/“.
AVEY-(PSS(J,2,I)+PSS(J+1,2,I)+PSS(J+1,2,I+1)+PSS(J,2,I+1))/“.
AVEZ-(PSS(J,3,I)+PSS(J+1,3,I)+PSS(J+1,3,I+l)+PSS(J,3,I+l))/“.

Sum forces.
CALL SUMFOR (SFORCE,I,J,X,Y,Z,AVEX,AVEY,AVEZ)
Calculate forces cosines for next row.

IF (SFORCE(J,1,I,lO).LE..OOOl.AND.SFORCE(J,2,I,10).LE..0001
+ .AND.SFORCE(J,3,I,lO).LE..OOOl) THEN
XNS(I)=O.
YNS(I)-0.
ZNS(I)-o.
ELSE
XNS(I)--(SFORCE(J,1,I,lO))/SQRT(SFORCE(J,1,I,lO)**2+
SFORCE(J,2,I,10)**2+SFORCE(J,3,I,lO)**2)
YNS(I)--(SFORCE(J,2,I,lO))/SQRT(SFORCE(J,1,I,10)**2+
+ SFORCE(J,2,I,10)**2+SFORCE(J,3,I,10)**2)
ZNS(I)--(SFORCE(J,3,I,10))ISQRT(SFORCE(J,1,I,lO)**2+
+ SFORCE(J,2,I,IO)**2+SFORCE(J,3,I,10)**2)
ENDIF

+

CONTINUE
Compute forces on moldboard.

WRITE (LUNDB,'(////A)') ' MOLDBOARD ELEMENT FORCES'
Do 20 J-1,JJM-1
DO 20 I-1,0LM-1

Assign force cosines.

IF (J.EQ.1) THEN
Ic-l
IF ((PMS(J,2,I).GE.PSS(JJS,2,IC).AND.PMS(J,2,I).LT.
+ PSS(JJS,2,IC+1)).OR.I.EQ.1) GOTO 16
IC=1C+1
GOTO 15

B31

16 XN=XNs(IC)
YN-YNS(IC)
ZN=ZNS(IC)
ELSE
XN-XNM(I)
YN-YNM(I)
ZN-ZNM(I)
ENDIF

Calculate element forces.

000

CALL ELEFOR (PMS,I,J,MFORCE,'M',XN,YN,ZN,IC)
AVEx-(PMS(J,1,I)+PMS(J+1,1,I)+PMS(J+1,1,I+1)+PMS(J,1,I+1))/4.
AVEY-(PMS(J,2,I)+PMS(J+1,2,I)+PMS(J+1,2,I+l)+PMS(J,2,I+1))/4.
AVEz-(PMS(J,3,I)+PMS(J+1,3,I)+PMS(J+1,3,I+1)+PMS(J,3,I+1))/4.
C Sum forces.

CALL SUMFOR (MFORCE,I,J,X,Y,Z,AVEX,AVEY,AVEZ)

Calculate force cosines for next row.

000

IF (MFORCE(J,1,I,lO).LE..OOOl.AND.MFORCE(J,2,I,10).LE.0.
+ .AND.MFORCE(J,3,I,10).LE..0001) THEN
XNM(I)-0.
YNM(I)-0.
ZNM(I)-o.
ELSE
XNM(I)--(MFORCE(J,1,I,10))/SQRT(MFORCE(J,1,I,10)**2+
+ MFORCE(J,2,I,lO)**2+MFORCE(J,3,I,10)**2)
YNM(I)--(MFORCE(J,2,I,10))/SQRT(MFORCE(J,1,I,10)**2+
+ MFORCE(J,2,I,10)**2+MFORCE(J,3,I,lO)**2)
ZNM(I)--(MFORCE(J,3,I,10))/SQRT(MFORCE(J,1,I,10)**2+
+ MFORCE(J,2,I,10)**2+MFORCE(J,3,I,lO)**2)
ENDIF

O CONTINUE

Sum forces over entire tool.

000N0

DO 30 K-l,3
TOTF0R(R,1)-o.
TOTMOM(K,l)-O.
TOTFOR(R,2)-o.
TOTMOM(K,2)-o.

DO 40 J-l,JJS-l
DO 40 I-l,OLS-l
TOTMOM(K,l)-TOTMOM(K,l)+SFORCE(J,K,I,9)
TOTFOR(K,1)-TOTFOR(K,1)+SFORCE(J,K, I, 10)
40 CONTINUE

DO 30 J'1,JJM-1

50
60

F‘C5C5C)

1000

+

B32

DO 30 I=l,OLM-l
TOTMOM(K,2)-TOTMOM(K,2)+MFORCE(J,K,I,9)
TOTFOR(K,2)=TOTFOR(K,2)+MFORCE(J,K,I,10)

CONTINUE

WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE

WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE

WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE

(LUNDB.'(////)')

(LUNDB,50) ' TOOL TOTALS','X','Y','Z'
(LUNDB,'(/,A)') ' FORCES (N)'

(LUNDB,60) 'SHARE',(TOTFOR(R,1),K=1,3)
(LUNDB,60) 'MOLDBOARD',(TOTFOR(K,2),R-1,3)

 

(LUNDB,'(T33,A)') '
(LUNDB,60) 'TOTAL',(TOTFOR(K,1)+TOTFOR(K,2),K-l,3)
(LUNDB,'(//A)') ' MOMENTS (N-M)‘

(LUNDB,60) 'SHARE',(TOTMOM(K,1),K-l,3)

(LUNDB,60) 'MOLDBOARD',(TOTMOM(R,2),R-1,3)

 

(LUNDB,'(T33,A)') '
(LUNDB,6O) 'TOTAL',(TOTMOM(R,1)+T01MOM(R,2),R-1,3)
(LUNDB,60) 'LOCATION OF ZERO MOMENT',x,Y,z

(LUNDA.'(//)')

(LUNDA,50) ' TOOL TOTALS','x','Y','z'
(LUNDA,'(/,A)') ' FORCES (N)'
(LUNDA,60) 'SHARE',(TOTFOR(R,1),K=1,3)

'(LUNDA,60) 'MOLDBOARD',(TOTFOR(K,2),K'1,3)

 

(LUNDA,'(T33,A)') '
(LUNDA,60) 'TOTAL',(TOTFOR(K,l)+TOTFOR(K,2),K=I,3)

(LUNTI,50) ' TOOL TOTALS','X',
(LUNTI.'(/,A)') ' FORCES (N)'
(LUNTI,60) 'SHARE',(TOTFOR(K,1),K-l,3)
(LUNTI,60) 'MOLDBOARD',(TOTFOR(K,2),K-1,3)

IYI’VZO

 

(LUNTI,'(T33,A)') '
(LUNTI,60) 'TOTAL‘,(TOTFOR(K,l)+TOTFOR(K,2),K-1,3)
(LUNTI,'(//A)') ' MOMENTS (N-M)'

(LUNTI,60) 'SHARE',(TOTMOM(K,1),K=1,3)

(LUNTI,60) 'MOLDBOARD',(TOTMOM(R,2),K-l,3)

 

(LUNTI,'(T33,A)') '
(LUNTI,60) 'TOTAL',(TOTMOM(K,l)+TOTMOM(K,2),K=1,3)
(LUNTI,60) 'LOCATION OF ZERO MOMENT',X,Y,Z

FORMAT (A,T28,3(TR11,A))
FORMAT (A,T33,3(F10.2,2X))

CLOSE

RETURN

(LUNDB)

Write message.

WRITE
WRITE
WRITE

FORHAT ( '******************** ATTENTION ********************v’/
' SOIL IS NOT IN CONTACT WITH SURFACE AT POINT ',12,12)

CLOSE

(LUNTI,1000) II,JJ
(LUNDA,1000) II,JJ
(LUNDB,1000) II,JJ

(LUNDB)

C

B33

RETURN

C*********************

C

000000000000000000000000000000000000000000000000

SUBROUTINE ELEFOR (PMS,I,J,FORCE,TOOLP,XN,YN,ZN,IC)
This subroutine calculates the soil forces on an element.

Array definitions:

PMS - coordinates of surface description points.
I,J - pointers indicating which element in PMS is
to have forces calculated for.
FORCE - calculated forces for elements where
lst and 3rd dimensions indicate element
2nd dimension indicates coordinate
l - x; 2 - y; 3 - 2
4th indicates which force
- normal load
- soil-metal friction
adhesion
- moments for each element
- total force for each element
VELM - velocity of soil over moldboard element.
VELS - velocity of soil over share element.

UIL‘UWNH
l

Variable definitions:

TOOLP - indicates which tool part is being calculated.
8 for share, M for moldboard (character string)

SOILTY - name of soil type (character string).

GAMMA - soil bulk density (kg/cu m).

COMES - soil cohesion (kpa).

PHI - soil internal friction (deg).

MUP - soil-metal friction.

ADHES - soil-metal adhesion (kpa).

VELTL - tool velocity (m/s).

DEPTHF - plowing depth (m).

XN,YN,ZN - force direction cosines.
IC - pointer to indicate which share element is before
moldboard element.

MU - internal soil friction coefficient.

F0 - contact area of soil block and tool.
Fl - area of forward shear failure surface.
WEIGHT - weight of soil block.

ACCEL - soil acceleration.

NO - normal load on tool.

LUNTI - logical unit number for terminal input.
LUNTO - logical unit number for terminal output.
LUNDA - logical unit number for data file.

00000000

000

000

834
LUNDB - logical unit number for data file.
Subroutines called:

ANGLE - computes angles of angles elements.
TRANEL - translates elements.
DISTAN - computes distance soil travels over element.

COMMON /A/GAMMA,COHES,PHI,MUP,ADHES,VELTL,DEPTHF
COMMON /LUN/LUNTI,LUNTO
COMMON /LUND/LUNDA,LUNDB

DIMENSION PMS(ll,3,ll),FORCE(ll,3,ll,5),VELM(ll),VELS(ll)
REAL MU,MUP,N0

CHARACTER TOOLP*1

SAVE VELM,VELS

DATA PI/3.l“lS9265“/

Calculate element orientation

CALL ANGLE (PMS,I,J,1,2,ALPHA)
CALL ANGLE (PMS,I,J,3,1,DELTA)
PHIP=ATAN(SIN(ALPHA)*TAN(DELTA))
BETA-.5*(90.-PHI)/(180./PI)

Compute other parameters

IF (J.EQ.1.AND.TOOLP.EQ.'S') THEN
COHESE-COHES*1000.

ELSE
COHESE-o.

ENDIF

Mu-TAN(PHI/(180./PI))

CALL TRANEL (PMS,I,J,(90./(180./PI)-ALPHA),-PHIP,DIMY1,DIMZl)
FO-DIMY1*DIMZI

CALL TRANEL (PMS,I,J,(90./(180./PI)-ALPHA),-BETA,DIMY,DIMZ)
Fl-DEPTHF/SIN(BETA)*DIMY

WEIGHTsGAMMA*9.81*F1*DIMz

Calculate force direction cosines.

IF (J.EQ.1.AND.TOOLP.EQ.'S') THEN
XNl-SIN(DELTA)

YNl'Oe
ZNl-Cos(DELTA)

ELSE
XN-YN/TAN(ALPHA)+ZN*TAN(DELTA)
XNl-XN/SQRT(XN**2+YN**2+ZN**2)
YNl-YN/SQRT(XN**2+YN**2+ZN**2)
ZN1a2N/SQRT(XN**2+YN**2+EN**2)

ENDIF

Compute accelerations.

B35

CALL DISTAN (PMS,I,J,XN1,YN1,ZN1,(90./(l80./PI)-ALPHA),
+ -PHIP,DIST)

000 000

000

SDIST-DIM21*COS(PHIP)/SIN(BETA)
TIME-DIST/VELTL
IF (J.EQ.1) THEN
IF (TOOLP.EQ.'S') THEN
VELI -0 e
VELs(I)-SDIST/TIME
VELz-VELS(I)
ELSE
VEL1-VELS(IC)
VELM(I)-SDIST/TIME
VEL2-VELM(I)
ENDIF
ELSE
IF (TOOLP.EQ.'S') THEN
VELl-VELS(I)
VELS(I)-SDIST/TIME
VEL2-VELS(I)
ELSE
VEL1-VELM(I)
VELM(I)-SDIST/TIME
VEL2-VELM(I)
ENDIF
ENDIF
ACCEL-GAMMA*Fl*DIMZ*(VEL2-VEL1)/TIME

NO'(COHESE*F1+ACCEL)*COS(PHIP-BETA)+WEIGHT*SIN(PHIP)
IF (NO.LE.O.) N030.

Calculate normal load.

FORCE(J,l,I,l)-N0*COS(PHIP)*SIN(ALPHA)
FORCE(J,2,I,1)I-NO*COS(PHIP)*COS(ALPHA)
FORCE(J,3,I,l)--NO*SIN(PHIP)

Calculate Soil-metal frictional force.

FORCE(J,1,I,2)-N0*MUP*(XN1)
FORCE(J,2,I,2)=N0*MUP*(YN1)
FORCE(J,3,I,Z)-NO*MUP*(ZN1)

Calculate adhesion force.
FORCE(J,1,I,3)-ADHES*F0*(XN1)
FORCE(J,2,I,3)-ADHES*F0*(YNI)
FORCE(J,3,I,3)=ADHES*F0*(ZN1)

Print out info to LUNDB.

WRITE (LUNDB,'(//A)') 'ELEMENT COORDINATES'
WRITE (LUNDB,'(11X,3(A,9X))') 'X','Y','Z'
WRITE (LUNDB,3S) I,J,PMS(J,I,I),PMS(J,2,I),PMS(J,3,I)

B36

WRITE (LUNDB,35) I,J+1,PMS(J+1,1,I),PMS(J+1,2,I),PMS(J+1,3,I)
WRITE (LUNDB,35) I+l,J+1,PMS(J+l,1,I+1),PMS(J+1,2,I+1),

+ PMS(J+1,3,I+1)

WRITE (LUNDB,35) I+1,J,PMS(J,1,I+l),PMS(J,2,I+1),PMS(J,3,I+1)

C
WRITE (LUNDB,'(A,T28,3(TR11,A))') ' FORCES (N)','X','Y','Z'
WRITE (LUNDB,50)'NORMAL LOAD',(FORCE(J,K,I,1),K=1,3)
WRITE (LUNDB,50)'SOIL-METAL FRICTION',(FORCE(J,K,I,2),K=1,3)
WRITE (LUNDB,50)'ADHESION',(FORCE(J,K,I,3),K=1,3)

C

35 FORMAT (12,3X,12,3(F9.3))
50 FORMAT (A,T33,3(F10.2,2X))

C

RETURN

END
C
C********************
C

SUBROUTINE SUMFOR (FORCE,I,J,X,Y,Z,AVEX,AVEY,AVEZ)
C
C This subroutine sums forces for an element.
C
C Array definitions:
C
C FORCE - forces values to be summed.
C
C Variable definitions:
c .
C I,J - pointers indication which element is being
C considered.
C X,Y,Z - coordinates of point force moments are being
C summed about.
C AVEX,AVEY,AVEZ - coordinates of center of point.
C
C LUNDA - logical unit number for data file.
C LUNDB - logical unit number for data file.
C

COMMON /LUND/LUNDA,LUNDB

DIMENSION FORCE(11,3,11,5)
C

FORCE(J,K,I,S)=0.

5 CONTINUE
C
C Sum forces.
C

no 10 Rs1,3
FORCE(J,1,I,5)-FORCE(J,1,I,5)+FORCE(J,1,I,K)
FORCE(J,2,I,S)-FORCE(J,2,I,5)+FORCE(J,2,I,K)
FORCE(J,3,I,5)-FORCE(J,3,I,5)+FORCE(J,3,I,K)
o CONTINUE

Calculate moments.

0001-

B37

FORCE(J,l,I,“)=(FORCE(J,2,I,5)*(Z-AVEZ))
+ -(FORCE(J,3,I,lO)*(Y-AVEY))

FORCE(J,2,I,“)--(FORCE(J,l,I,5)*(Z-AVEZ))
+ +(FORCE(J,3,I,S)*(X-AVEX))

FORCE(J,3,I,“)=(FORCE(J,l,I,5)*(Y-AVEY))
+ -(FORCE(J,2,I,5)*(X-AVEX))

 

WRITE (LUNDB,'(T33,36A)')' -------- '
WRITE (LUNDB,30) 'TOTALS',(FORCE(J,K,I,10),K-l,3)
WRITE (LUNDB,'(/)')
WRITE (LUNDB,30) 'DISTANCE TO ZERO MOMENT',x-AVEx,Y-AVEY,z-AVEZ
WRITE (LUNDB,30) 'MOMENT (N-M)',(FORCE(J,R,I,9),R=1,3)

3o FORMAT (A,T33,3(F10.2,2x))

RETURN
END

APPENDIX C

PLOW DATA FILES

PLOW DATA FILES

 

This section contains the data files used to generate

the information for this research. Each file will be

referred to by the name as it exists on the Prime 750

located at The Case Center for Computer Aided Design,

Michigan State University. A brief explanation of the

function of each file will also be included. Those

interested in a detailed description of each value are

encouraged to review the fortran source code for the

variable name the values are read into.
The first file is D_TINE. It is the geometry file for

the 63 mm wide inclined blade.

0.357 0.063 0.250
0.000 0.000 0.000
0.083 0.000 0.083
0.167 0.000 0.167
0.089 0.000 0.062
0.000 0.333 0.000
0.083 0.333 0.083
0.167 0.333 0.167
0.250 0.333 0.250
0.000 0.667 0.000
0.083 0.667 0.083
0.167 0.667 0.167
0.250 0.667 0.250
0.000 0.063 0.000
0.083 1.000 0.083
0.167 1.000 0.167
0.089 0.063 0.062
0.089 0.000 0.062
0.178 0.000 0.125
0.268 0.000 0.187
0.357 0.000 0.250
0.089 0.021 0.062
0.500 0.333 0.500

0.750
0.357
0.089
0.500
0.750
0-357
0.089
0.178
0.268
0-357

The

0.3“59
0.0000
0.0359
0.0719
0.1078
0.1153
0.1513
0.1658
0.2018
0.2306
0.2“52
0.2811
0.2957
0.3“59
0.3605
0.3751
0.3896
0.1078
0.18“2
0.2“76
0.2682
0.1958
0.2722
0.703“
0.7239
0.2838
0.“962
0.8185
0.8128
0.3718
0.58“2
0.8382
0.8325

0.333
0.021
0.0“2
0.667
0.667
0.0“2
0.063
0.063
0.063
0.063

0.“07“
0.0000
0.0000
0.0000
0.0000
0.1358
0.1358
0.106“
0.106“
0.2716
0.2“22
0.2“22
0.2128
0.“07“
0.3780
0.3“86
0.3192
0.0000
0.0000
0.0000

0.0991
0.0991
0.“597
0.“572
0.1981
0.3073
0.53“3
0.2972
0.“06“
0.“890
0.56“3

0.750
0.250
0.062
00500
0.750
0.250
0.062
0.125
0.187
0.250

0.259“
0.0000
0.0160
0.0320
0.0“80
0.0000
0.0160
0.0321
0.0“81
0.0000
0.0160
0.0321
0.0“81
0.0000
0.0160
0.0321
0.0“81
0.0“80
0.0826
0.1651
0.259“
0.0“7“
0.0819
0.“391
0.533“
0.0“67
0.0“67
0.1671
0.279“
0.0“61
0.0“61
0.1270

C2

following file is D_FORD.

It is the geometry file

for the 16 inch wide Ford moldboard plow.

C3

The following file is

D_FORDPLOW.

It contains the

descriptor points for the 16 inch wide Ford moldboard plow.

\OCDNGU'It'WNI—I

file.

moldboard plow model

0.198“
0.0000
0.1285
0.2301
0.3“59
0.3896
0.2886
0.1078
0.1723
0.212“
0.2680
0.360“
0.“627
0.6070
0.“955
0.3718
0.5638
0.639“
0.7817
0.2“92
0.2682
0.3033
0.3235
0.“3“0
0.5“53
0.6663
0.7501
0.7530
0.6“72
005081
0.“003
0.7990
0.8325
0.8181

.1035
.0000
.1502
.2673
.“07“
.3192
0.2053
0.0000
0.002“

000000

-0.0030

0.0635
0.1827
0.2861
0.“03“
0.2972
0.3059
0.5011

-0.0039
-0.0025

0.0“9“
0.0“75
0.1578
0.28“1
0.“223
0.5167
0.“896
0.3783
0.2“18
0.15“2
0.5605
0.56“3
005338

0.0“89
0.0000
0.001“
0.0006
0.0000
0.0“81
0.0“87
0.0“80
0.0918
0.1387
0.1266
0.1056
0.1029
0.0800
0.0551
0.0“61
0.2002
0.19“1
001581
0.2017
0.259“
0.2067
0.2830
0.3““0
0.3816
0.366“
0.3“02
0.2538
0.2608
0.2“18
0.1868
0.3096
0-2393
0.1893

The thirteen files that follow are operating conditions

The first file is the operating conditions for the

remainder

of the files,

verification

(Table

6.2). For the

the number at the end of each file

name corresponds to the case number of Figure 6.1.

C“

D_SANDYLOAM
Norfolk Sandy Loam
1770.

1907

3“.8

.“

0.
1.3“
.203

D_LMCLAYl

Low Moisture Clay-Case 1
103“.

35.9

38.

.““5

1.8

1.3“

.25

D_LMCLAY2

Low Moisture Clay-Case 2
1007.

17.2

“5.

.625

12.1

1.3“

.25

D_LMCLAY3
Low Moisture Clay-Case 3
1089. '
“8.3
35.
.36“
11.
1.3“
025

D_HMCLAY“

High Moisture Clay-Case “
1025.

15.5

39.

.625

12.1

1.3“

.25

C5

D_HMCLAYS

High Moisture Clay-Case 5
980.

55-2

26.

.“66

6.9

1.3“

.25

D_HMCLAY6

High Moisture Clay-Case 6
10“3.

37.9

2“.

.36“

10.5

1.3“

.25

D_LMSAND7

Low Moisture Sandy Clay-Case 7
1379.

6.9

36.

.““5

2.8

1.3“

.25

D_LMSAND8

Low Moisture Sandy Clay-Case 8
1361.

17.2

3“.

.306

5.”

1.3“

.25

D_LMSAND9
Low Moisture Sandy Clay-Case 9
1“15.

C6

D_HMSANDlO

High Moisture Sandy Clay-Case 10
1“06.

10.3

33.

.36“

“.6

103“

.25

D_HMSANDll
High Moisture Sandy Clay-Case 11
1288. -
10.3
32.
.36“
9.
1031'
.25

D_HMSAND12

High Moisture Sandy Clay-Case 12
1325.

17.2

“0.

.510

3.5

1.3“

.25

APPENDIX D

MEASURED DATA

MEASURED DATA

 

The following file contains the data points measured on

the 16 inch wide Ford moldboard plow.

the entire

times.

moves were used to piece together the data.

The

surface,

point

as D_FORDPLOW in

the

numbers

Appendix C is

manipulation of this file.

To adequately measure

plow had to be positioned three

measured before and after the
The file listed
of the

the result

Points 1-18 were taken at the first plow position.

1 0.000

2 8.605

3 0.999

“ -“.980

5 -11.951

6 -11.135

7 -5.253

8 5.305

9 3.296

10 2.237
11 -1.110
12 -6.875
13 -12.5“8
1“ -19089u
15 “150303
16 -9093?
17 -16.0“6
18 -20.123
Points

11 0.000
17 15.201
18 19.251
19 26.578
20 -2.035
21 “10305
22 0.936

0.000
2.327
2.111
2.016
1.901
-0.218
-0.109
0.1“9
-1.7“2
-3.699
-3.290
—2.566
-2.607
”10952
-0.732
-0.l23
-6.679
-60551

0.000
2.210
1.75“
-0.030
3.176
5-“38
3.216

5.200
6.568
8.203
9.“37
11.050
7.380
6.232
“.127
2.781
1.780
2.500
3.965
“.753
“.917
6.017
6.958
3.220
3.63“

11 and 17-31 were

8.000
7.291
6.882
6.500
9.“67
9.72“
9.077

01

taken at the second position.

23 1.658
2“ 7.852
25 1“.“07
26 21.“55
27 26.281
28 25.572
29 19.625
30 11.987
31 6.“31
Points
position.
17 -11.361
18 -70303
19 0.000
27 -o.281
28 -00980
32 2.227
33 3.253
3“ 1.99“

6.218
8.283
90380
8.372
7.061
3.757
“.360
“.022
2.120

17.19;

2.202
1.762
0.000
7.103
3.792
5.768
30005
1.13“

27,

D2

9.“l7
8.““1
7.050
5.658
“.771
5.918
6.855
7.855
8.169

28 and 32-3“ were taken at the third

7.288
6.880
6.500
“.780

.5.91“

“.653

5.520
6.260

LIST OF REFERENCES

LIST OF REFERENCES

 

Carlson, E.C., 1961. Plow and Computers. Agricultural
Engineering. “2(6):292-295, 307.

Emmett, W.T. Jr., 1983. Computer-Aided Performance
Predictions of Coated Moldboard Plows. M.S. Thesis
Cornell University, Ithaca, NY.

Emmett, W.T. Jr. and W.W. Gunkel, 1983. A Computer-Aided
Analysis of the Performance of Coated Moldboard Plow

Surfaces. American Society of Agricultural Engineers
Paper No. 83-10“6.

Faux, I.D. and M.J. Pratt, 1979. Computational Geometry for
Design and Manufacture. Ellis Horwood Limited, West
Sussex, England.

Foley, J.D. and A. VanDam, 1982. Fundamentals of
Interactive Computer Graphics. Addison Wesley
Publishing Company, Reading, Massachusetts.

Gill, W.R. and F.E. VandenBerg, 1967. Soil Dynamics in

Tillage and Traction. Agricultural Handbook No. 316.
U.S. Government Printing Office, Washington, D.C.

Hettiaratchi, D.R.P. and A.R. Reece, 1967. Symmetrical
Three-Dimensional Soil Failure. Journal of
Terramechanics. “(3):“5-67.

McKyes, E., 1977. The Calculation of Draft Forces and Soil
Failure Boundaries of Narrow Cutting Blades.

Transactions of American Society of Agricultural
Engineers. 15(2):111-116.

Nichols, M.L., 1931a. Dynamic Properties of Soil. I. An
Explanation of the Dynamic Properties of Soils by
Means of Colloidal Films. Agricultural Engineering.
12(7):259-26“.

Nichols, M.L., 1931b. The Dynamic Properties of Soil.
II. Soil and Metal Friction. Agricultural
Engineering. 12(8):32l-32“.

Nichols, M.L. and T.H. Kummer, 1932. The Dynamic Properties
of Soil. IV. A Method of Analysis of Plow Moldboard
Design Based upon Dynamic Properties of Soil.
Agricultural Engineering. 13(11):279-285.

Nichols, M.L. and I.F. Reed, 193“. Soil Dynamics: VI.
Physical Reactions of Soils to Moldboard Surfaces.
Agricultural Engineering. 15(6):187-190.

Nichols, M.L., I.F. Reed and C.A. Reaves, 1958. Soil
Reaction: To Plow Share Design. Agricultural
Engineering. 39(6):336-339-

O'Callaghan, J.R. and J.G McCoy, 1965. The Handling of Soil
by Mouldboard Ploughs. Journal of Agricultural
Engineering Research. 10(1):23-25.

Orlandea, N., S.H. Chen and T. Berenyi, 1983. A Study of
Soil-Tool Interactions. Transactions of American
Society of Agricultural Engineers. 21(5):1619—1625.

Perumpral, J.V., R.D. Grisso and 0.8. Desai, 1983. A
Soil-Tool Model Based on Limit Equilibrium Analysis.
Transactions of American Society of Agricultural
Engineers. 21(3):991-995.

Plasse, R., G.S.V. Raghavan and F. McKyes, 1983. Simulation
of Narrow Blade Performance in Different Soils.
American Society of Agricultural Engineers Paper No.

83-1539.

Rowe, R.J. and K.K. Barnes, 1961. Influence of Speed on
Elements of Draft of a Tillage Tool. Transactions of
the American Society of Agricultural Engineers.

6(1):55-57.

Siemens, J.C., J.A. Weber and T.H. Thorburn, 1965.
Mechanics of Soil as Influenced by Model Tillage
Tools. Transactions of American Society of
Agricultural Engineers. 8(1):1-7.

Soehne, W., 1959. [Investigations on the Shape of Plough
Bodies for High Speeds.) Grundlagen der Landtechnik.
11:22-39, (National Institute of Agricultural
Engineering., English Translation 87.)

White, E.A., 1918. A Study of the Plow Bottom and its
Action Upon the Furrow Slice. Journal of Agricultural
Research. 12:1“9-182, illus.